text:
{-# LANGUAGE FlexibleContexts #-}

module Blash where

import qualified Data.Vector.Generic as V
import Data.Vector.Generic ((!))
import Data.Complex
import BlashImpl (fZERO, stride, isample, map_reduce, Size, Inc, Index, copyM, axpyM)

copy :: (V.Vector v a)
         => Size
         -> v a
         -> Inc
         -> v a
         -> Inc
         -> v a
copy n dx incx dy incy = V.modify modifier dy
  where
    modifier ys = copyM n dx incx ys incy

axpy :: (Ord a, Num a, V.Vector v a)
         => Size
         -> a
         -> v a
         -> Inc
         -> v a
         -> Inc
         -> v a
axpy n da dx incx dy incy = V.modify modifier dy
  where
    modifier ys = axpyM n da dx incx ys incy

-- /* dot product dx dot dy. */
dot :: (Num a, V.Vector v a)
        => Size
        -> v a
        -> Inc
        -> v a
        -> Inc
        -> a
dot n _ _ _ _ | n <= 0 = fZERO
dot n dx incx dy incy = sum prods
  where
    prods = flip map [0..(n-1)] f
    f i = (dx ! ix i) * (dy ! iy i)
    ix = stride n incx
    iy = stride n incy

-- /* dot product dx dot dy where dx & dy are complex. */
dotc :: (RealFloat a, V.Vector v (Complex a))
        => Size
        -> v (Complex a)
        -> Inc
        -> v (Complex a)
        -> Inc
        -> Complex a
dotc n _ _ _ _ | n <= 0 = fZERO :+ fZERO
dotc n dx incx dy incy = sum prods
  where
    prods = flip map [0..(n-1)] f
    f i = (conjugate (dx ! ix i)) * (dy ! iy i)
    ix = stride n incx
    iy = stride n incy

-- /* compute the L2 norm of array DX of length N, stride INCX */
-- TODO should really be (Num a) - or split between Complex and Real
nrm2 :: (Ord a, Floating a, V.Vector v a)
        => Size
        -> v a
        -> Inc
        -> a
nrm2 n _  incx | n <= 0 || incx <= 0 = fZERO
nrm2 n dx _ | n == 1 = abs $ dx ! 0
nrm2 n dx incx =
  let dx' = isample n dx incx
      xmax = V.maximum $ V.map abs dx'
  in
    case xmax == fZERO of
    True -> fZERO
    False ->
      let scale = 1.0 / xmax
          scaled_dx' = V.map (* scale) dx'
          n' = V.length scaled_dx'
          scaled_drnm2 = sqrt $ dot n' scaled_dx' 1 scaled_dx' 1
      in
        xmax * scaled_drnm2

-- DASUM takes the sum of the absolute values.
asum :: (Num a, V.Vector v a)
        => Size
        -> v a
        -> Inc
        -> a
asum = map_reduce abs V.sum

-- DAMAX takes the max of the absolute values.
iamax :: (Ord a, Num a, V.Vector v a)
        => Size
        -> v a
        -> Inc
        -> Index
iamax = map_reduce abs V.maxIndex

max_stars_repo_name: wyn/blash

max_issues_repo_name: wyn/blash

max_forks_repo_name: wyn/blash


text:
{-# LANGUAGE BangPatterns, FlexibleContexts #-}
-- |
-- Module    : Statistics.Transform
-- Copyright : (c) 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Fourier-related transformations of mathematical functions.
--
-- These functions are written for simplicity and correctness, not
-- speed.  If you need a fast FFT implementation for your application,
-- you should strongly consider using a library of FFTW bindings
-- instead.

module Transform ( kPointFFT ) where

import Control.Monad         (when)
import Control.Monad.ST      (ST)
import Data.Bits             (shiftL, shiftR, (.&.), (.|.))
import Data.Complex          (Complex(..), conjugate, realPart, magnitude)
import qualified Data.Vector.Generic         as G
import qualified Data.Vector.Generic.Mutable as M
import qualified Data.Vector.Unboxed         as U

type CD = Complex Float

-- | /O(log n)/ Compute the logarithm in base 2 of the given value.
log2 :: Int -> Int
log2 v0
    | v0 <= 0   = error "Statistics.Math.log2: invalid input"
    | otherwise = go 5 0 v0
  where
    go !i !r !v | i == -1        = r
                | v .&. b i /= 0 = let si = U.unsafeIndex sv i
                                   in go (i-1) (r .|. si) (v `shiftR` si)
                | otherwise      = go (i-1) r v
    b = U.unsafeIndex bv
    !bv = U.fromList [0x2, 0xc, 0xf0, 0xff00, 0xffff0000, 0xffffffff00000000]
    !sv = U.fromList [1,2,4,8,16,32]

hanning :: Float -> Float -> Float
hanning m n = 0.5 - 0.5 * cos(2 * pi * n / m)

kPointFFT :: U.Vector Float -> U.Vector Float
kPointFFT !vec = mag where
    winlen = realToFrac $ G.length vec - 1
    window = G.imap (\i x -> ((hanning winlen (realToFrac i)) * x) :+ 0.0) vec
    fftres = halfFFT window
    mag    = G.map (\x -> let val = 20 * (logBase 10 $ magnitude x) in if(val < 0) then 0 else val) fftres

halfFFT :: U.Vector CD -> U.Vector CD
halfFFT !vec = U.take (U.length vec `div` 2) (fft vec)

-- | Radix-2 decimation-in-time fast Fourier transform.
fft :: U.Vector CD -> U.Vector CD
fft !v = G.create $ do
          mv <- G.thaw v
          mfft mv
          return mv

mfft :: (M.MVector v CD) => v s CD -> ST s ()
mfft vec
    | 1 `shiftL` m /= len = error "Statistics.Transform.fft: bad vector size"
    | otherwise           = bitReverse 0 0
 where
  bitReverse i j | i == len-1 = stage 0 1
                 | otherwise  = do
    when (i < j) $ M.swap vec i j
    let inner k l | k <= l    = inner (k `shiftR` 1) (l-k)
                  | otherwise = bitReverse (i+1) (l+k)
    inner (len `shiftR` 1) j
  stage l !l1 | l == m    = return ()
              | otherwise = do
    let !l2 = l1 `shiftL` 1
        !e  = -6.283185307179586/fromIntegral l2
        flight j !a | j == l1   = stage (l+1) l2
                    | otherwise = do
          let butterfly i | i >= len  = flight (j+1) (a+e)
                          | otherwise = do
                let i1 = i + l1
                xi1 :+ yi1 <- M.read vec i1
                let !c = cos a
                    !s = sin a
                    d  = (c*xi1 - s*yi1) :+ (s*xi1 + c*yi1)
                ci <- M.read vec i
                M.write vec i1 (ci - d)
                M.write vec i (ci + d)
                butterfly (i+l2)
          butterfly j
    flight 0 0
  len = M.length vec
  m   = log2 len

fi :: Int -> CD
fi = fromIntegral

halve :: Int -> Int
halve = (`shiftR` 1)

ext: hs

max_stars_repo_name: matthewleon/haskell-portaudio

max_issues_repo_name: matthewleon/haskell-portaudio

max_forks_repo_name: matthewleon/haskell-portaudio


text:
{-# LANGUAGE TemplateHaskell #-}

module Scheme.Types where
import Data.DeriveTH
import Data.Maybe
import Data.IORef

import Control.Monad.Error

import qualified Data.Complex as C (Complex(..), realPart)
import Data.Ratio (Ratio, (%), numerator, denominator)
import Data.Vector as V (Vector, foldl)

import Text.Parsec.Error
import Text.Parsec.Pos

import System.IO
import DeriveGetMaybe

type ThrowsErrorIO = ErrorT LispError IO

data Env = Env
    { readOnly :: [(String, LispVal)]
    , definitions :: IORef [(String, IORef LispVal)] }
    deriving (Eq)

data LispVal
    = Atom String
    | Bool Bool
    | Num LispNum
    | Character Char
    | String String
    | List [LispVal]
    | DottedList [LispVal] LispVal
    | Vector (V.Vector LispVal)
    | Func [String] (Maybe String) [LispVal] Env
    | PrimitiveFunc String
    | Port Handle
    deriving (Eq)

data LispNum
    = Int Integer
    | Rational (Ratio Integer)
    | Double Double
    | Complex (C.Complex Double)
    deriving (Show, Eq)

data LispError
    = ParserError ParseError
    | NumArgs Int [LispVal]
    | BadExpr String LispVal
    | TypeMismatch String String
    | UnboundVar String String
    | Default String

instance Show LispVal where
    show (Atom s)            = s
    show (Bool s)            = if s then "#t" else "#f"
    show (Num n)             = show n
    show (Character c)       = "\\#" ++ [c]
    show (String s)          = "\"" ++ s ++ "\""
    show (PrimitiveFunc f)   = "<primitive : " ++ f ++ ">"
    show (List vs)           = "(" ++ unwords (map show vs) ++ ")"
    show (DottedList vs dot) = "(" ++ unwords (map show vs)
                                   ++ " . " ++ show dot ++ ")"
    show (Vector vs)         = "#(" ++ V.foldl (\acc e -> show e ++ " " ++ acc) "" vs ++ ")"
    show (Func params vararg body _) =
        "(lambda (" ++ unwords params
        ++ maybe "" (" . "++) vararg ++ ") " ++ unwords (map show body)++ ")"
    show (Port h)            = "<port: " ++ show h ++ ">"

$( derive makeIs ''LispVal )
$( derive makeIs ''LispNum )
$( derive makeGet ''LispVal )

getInt :: LispNum -> Maybe Integer
getInt (Int n) = Just n
getInt _ = Nothing

getExact :: LispNum -> Maybe (Ratio Integer)
getExact (Int n) = Just (n % 1)
getExact (Rational n) = Just n
getExact _ = Nothing

getAsDouble :: LispNum -> Double
getAsDouble (Int n) = fromInteger n
getAsDouble (Rational n) = fromInteger (numerator n) / fromInteger (denominator n)
getAsDouble (Complex n) = C.realPart n
getAsDouble (Double n) = n

getAsComplex :: LispNum -> C.Complex Double
getAsComplex (Complex n) = n
getAsComplex n = getAsDouble n C.:+  0

---------AHHHHHHHHH!!!!!!!!
--the horror

--engulfs

instance Fractional LispNum where
    a / b = if isComplex a && isComplex b
            then Complex $ getAsComplex a / getAsComplex b
            else Double $ getAsDouble a / getAsDouble b
    recip a = (Complex $ 1 C.:+ 0) / a
    fromRational = Rational

instance Num LispNum where
    (Int a) + (Int b) = Int $ a + b
    (Rational a) + (Rational b) = Rational $ a + b
    (Double a) + (Double b) = Double $ a + b
    (Complex a) + (Complex  b) = Complex $ a + b
    a + b = uncurry (+) $ coerceNums a b

    (Int a) * (Int b) = Int $ a * b
    (Rational a) * (Rational b) = Rational $ a * b
    (Double a) * (Double b) = Double $ a * b
    (Complex a) * (Complex  b) = Complex $ a * b
    a * b = uncurry (*) $ coerceNums a b

    a - b = (\(a', b') -> a' + negate b') (coerceNums a b)
    abs (Int n) = Int           $ abs n
    abs (Rational n) = Rational $ abs n
    abs (Double n) = Double     $ abs n
    abs (Complex n) = Complex   $ abs n
    negate (Int n) = Int           $ negate n
    negate (Rational n) = Rational $ negate n
    negate (Double n) = Double     $ negate n
    negate (Complex n) = Complex   $ negate n
    signum (Int n) = Int           $ signum n
    signum (Rational n) = Rational $ signum n
    signum (Double n) = Double     $ signum n
    signum (Complex n) = Complex   $ signum n
    fromInteger = Int

coerceNums :: LispNum -> LispNum -> (LispNum, LispNum)
coerceNums al bl = fromJust $ msum
    [ do
         a <- getInt al
         b <- getInt bl
         return (Int a, Int b)
    , do
         a <- getExact al
         b <- getExact bl
         return (Rational a, Rational b)
    , if isComplex al && isComplex bl
        then Just (Complex $ getAsComplex al, Complex $ getAsComplex bl)
        else Nothing
    , Just (Double $ getAsDouble al, Double $ getAsDouble bl) ]

instance Show LispError where
         show (ParserError e)       = show e
         show (NumArgs n vs)        = "expected " ++ show n ++ " args, found " ++ show (length vs) ++ "\n\t" ++ unwords (map show vs)
         show (TypeMismatch want v) = "type mismatch. expected " ++ want ++ ", found " ++ show v
         show (BadExpr msg expr)    = msg ++ " : " ++ show expr
         show (UnboundVar msg var)  = msg ++ " : " ++ var
         show (Default msg)         = "error : " ++ msg

instance Error LispError where
    noMsg = Default "mysterious failure experienced."
    strMsg = Default

nullEnv readOnlys = newIORef [] >>= \e -> return (Env readOnlys e)

setVar :: Env -> String -> LispVal -> ThrowsErrorIO ()
setVar e@(Env readOnlys env) var lisp = liftIO (readIORef env) >>= \envList ->
       maybe (defineVar e var lisp)
             (const $ throwError $ Default "setting an existing variable") $ lookup var envList

getVar (Env readOnlys env) var = maybe findInDefinitions return (lookup var readOnlys)
    where findInDefinitions = liftIO (readIORef env) >>= maybe
                    (throwError $ Default $ "undefined atom \"" ++ var ++ "\"" )
                    (liftIO . readIORef) . lookup var

defineVar (Env readOnlys env) var lisp =
    do envList <- liftIO $ readIORef env
       maybe (return ()) (throwError . Default . ("redefining read-only "++) . show) (lookup var readOnlys)
       maybe (liftIO $ newIORef lisp >>= \wrapped -> writeIORef env ((var, wrapped):envList))
             (\oldLisp -> liftIO $ writeIORef oldLisp lisp) $ lookup var envList

showIOThrows :: (MonadIO m, Show e) => ThrowsErrorIO e -> m String
showIOThrows = (liftIO . runErrorT) >=> return . either show show

-- introspection!
errToLisp :: LispError -> LispVal
errToLisp (ParserError p) =
    List [Atom "parse-error"
        , List [Atom "source-name", String . sourceName $ errorPos p]
        , List [Atom "source-line", Num . Int . fromIntegral . sourceLine $ errorPos p]
        , List [Atom "source-column", Num . Int . fromIntegral . sourceColumn $ errorPos p]
        , List [Atom "message", String $ showErrorMessages "or" "unknown parse error"
                                    "expecting" "unexpected" "end of input"
                                    (errorMessages p)]]
errToLisp err = String . show $ err

max_stars_repo_name: oroguh/SCMinHS

max_issues_repo_name: oroguh/SCMinHS

max_forks_repo_name: oroguh/SCMinHS


text:
{-# LANGUAGE FlexibleContexts #-}

module Main (main) where

import Control.Monad (mzero)
import Control.Monad.IO.Class (liftIO)
import Data.Complex (Complex)
import Data.Foldable (toList)
import Data.Typeable (Typeable)
import System.Console.GetOpt
import Text.PrettyPrint.Mainland
import Text.PrettyPrint.Mainland.Class

import Spiral
import Spiral.Backend.C
import Spiral.Config
import Spiral.Driver
import Spiral.Exp
import Spiral.FFT.CooleyTukey
import Spiral.Monad
import Spiral.OpCount
import Spiral.Program
import Spiral.SPL
import Spiral.SPL.Run
import Spiral.Search
import Spiral.Search.FFTBreakdowns
import Spiral.Util.Uniq

main :: IO ()
main = defaultMainWith' options mempty $ \fs args -> do
    useComplexType <- asksConfig $ testDynFlag UseComplex
    n <- case args of
           [s] -> return (read s)
           _   -> return 4
    f <- formula fs n
    pprint f
    if useComplexType
      then toProgram ("hspiral_dft_" ++ show n) f >>= go
      else toProgram ("hspiral_dft_" ++ show n) (Re f) >>= go
  where
    go :: (Typed a, Num (Exp a)) => Program a -> Spiral ()
    go prog = do
      pprint prog
      ops <- countProgramOps prog
      resetUnique
      defs <- evalCg $ cgProgram prog
      outp <- asksConfig output
      case outp of
        Nothing -> return ()
        Just{}  -> writeOutput (toList defs)
      liftIO $ putDocLn $
          text "Multiplications:" <+> ppr (mulOps ops) </>
          text "      Additions:" <+> ppr (addOps ops) </>
          text "          Total:" <+> ppr (allOps ops)

-- The SPL formula for which we generate code and count operations.
formula :: MonadSpiral m => [Flag] -> Int -> m (SPL (Exp (Complex Double)))
formula fs n =
  case fs of
    [Dif]                -> return $ dif n
    [Dit]                -> return $ dit n
    [SplitRadix]         -> runSearch () splitRadixSearch (DFT n)
    [ConjPairSplitRadix] -> runSearch () conjSplitRadixSearch (DFT n)
    [ImpSplitRadix]      -> runSearch () impSplitRadixSearch (DFT n)
    _                    -> fail "Must specify exactly on of --dif, --dit, --split-radix, conj-split-radix, or --improved-split-radix"
  where
    splitRadixSearch :: (Typeable a, Typed a, MonadSpiral m)
                     => SPL (Exp a)
                     -> S s m (SPL (Exp a))
    splitRadixSearch (F n w) = splitRadixBreakdown n w
    splitRadixSearch _       = mzero

    conjSplitRadixSearch :: (Typeable a, Typed a, MonadSpiral m)
                         => SPL (Exp a)
                         -> S s m (SPL (Exp a))
    conjSplitRadixSearch (F n w) = conjPairSplitRadixBreakdown n w
    conjSplitRadixSearch _       = mzero

    impSplitRadixSearch :: (Typeable a, Typed a, MonadSpiral m)
                        => SPL (Exp a)
                        -> S s m (SPL (Exp a))
    impSplitRadixSearch (F n w) = improvedSplitRadixBreakdown n w
    impSplitRadixSearch _       = mzero

data Flag = Dif
          | Dit
          | SplitRadix
          | ConjPairSplitRadix
          | ImpSplitRadix
  deriving (Eq, Ord, Show)

options :: [OptDescr Flag]
options =
    [ Option [] ["dif"] (NoArg Dif)                             "Use DIF"
    , Option [] ["dit"] (NoArg Dit)                             "Use DIT"
    , Option [] ["split-radix"] (NoArg SplitRadix)              "Use split radix"
    , Option [] ["conj-split-radix"] (NoArg ConjPairSplitRadix) "Use conjugate pair split radix"
    , Option [] ["imp-split-radix"] (NoArg ImpSplitRadix)       "Use improved split radix"
    ]

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}

{-# LANGUAGE CPP                 #-}
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeOperators       #-}

module Test.Grenade.Layers.LeakyRelu where

import           Data.Proxy

import           GHC.TypeLits

import           Grenade.Core
import           Grenade.Layers.LeakyRelu

import qualified Numeric.LinearAlgebra        as LA
import           Numeric.LinearAlgebra.Static (L)
import qualified Numeric.LinearAlgebra.Static as H

import           Hedgehog
import qualified Hedgehog.Range               as Range

import           Test.Hedgehog.Compat
import           Test.Hedgehog.Hmatrix

prop_leaky_relu_correct_on_one_element = property $ do
  input :: S ('D1 1) <- forAll genOfShape
  alpha              <- forAll $ genRealNum (Range.constant 0 0.9)

  let layer   = LeakyRelu alpha :: LeakyRelu
      S1D out = snd $ runForwards layer input :: S ('D1 1)
      x       = (\(S1D a) -> H.extract a LA.! 0) input

  H.extract out LA.! 0 === if x >= 0 then x else x * alpha

prop_leaky_relu_correct_on_3D_array = property $ do
  height   :: Int <- forAll $ choose 2 100
  width    :: Int <- forAll $ choose 2 100
  channels :: Int <- forAll $ choose 2 100
  alpha           <- forAll $ genRealNum (Range.constant 0 0.9)

  case (someNatVal (fromIntegral height), someNatVal (fromIntegral width), someNatVal (fromIntegral channels)) of
    (Just (SomeNat (Proxy :: Proxy h)), Just (SomeNat (Proxy :: Proxy w)), Just (SomeNat (Proxy :: Proxy c))) -> do
      input :: S ('D3 h w c) <- forAll genOfShape

      let layer    = LeakyRelu alpha :: LeakyRelu
          S3D out  = snd $ runForwards layer input :: S ('D3 h w c)
          inp' = (\(S3D x) -> H.dmmap (\a -> if a >= 0 then a else a * alpha) x) input :: L (h * c) w

      H.extract inp' === H.extract out

tests :: IO Bool
tests = checkParallel $$(discover)

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-|
Module      : Qubism.StateVec
Description : Types and functions for quantum state vectors
Copyright   : (c) Keith Pearson, 2018
License     : MIT
Maintainer  : keith@qubitrot.org
-}

{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}

module Qubism.StateVec
  ( StateVec (UnsafeMkStateVec)
  , mkStateVec
  , mkStateVec'
  , mkQubit
  , normalize
  , tensor
  , collapse
  , measureQubit
  , measure
  , dimension
  ) where

-- For dependent typing
import GHC.TypeLits
import Data.Singletons
import Data.Singletons.TypeLits
import Data.Finite

import           Control.Monad.Random
import           Control.Monad.Trans.State.Strict
import           Data.Complex
import           Numeric.LinearAlgebra ((|>))
import qualified Numeric.LinearAlgebra as LA
import           Text.Printf

import Qubism.Algebra
import Qubism.CReg

newtype StateVec (n :: Nat) =
  UnsafeMkStateVec (LA.Vector C)

-- | To provide "close enough" equality testing.
instance Eq (StateVec n) where
  UnsafeMkStateVec zs == UnsafeMkStateVec ws =
     LA.norm_2 (zs - ws) < 0.000001

instance KnownNat n => VectorSpace (StateVec n) where
  zero = UnsafeMkStateVec $ internalLen (sing :: Sing n) |> repeat 0
  z .: (UnsafeMkStateVec a) = UnsafeMkStateVec $ LA.scalar z * a
  (UnsafeMkStateVec a) +: (UnsafeMkStateVec b) = UnsafeMkStateVec $ a + b
  neg (UnsafeMkStateVec a) = UnsafeMkStateVec $ -a

instance KnownNat n => HilbertSpace (StateVec n) where
  (UnsafeMkStateVec zs) <.> (UnsafeMkStateVec ws) = zs LA.<.> ws

instance KnownNat n => Show (StateVec n) where
  show (UnsafeMkStateVec zs) = foldl show' "" $ zip [0..] (LA.toList zs)
    where show' str (i,z) = str ++ pf z ++ "  " ++ basis i ++ "\n"
          pf z    =  printf "% 6.4f" (realPart z) ++ "  + "
                  ++ printf "% 6.4f" (imagPart z) ++ "i"
          basis i = let bit j = if i `quot` 2^(n-j-1) `mod` 2 == 0
                                then '0' else '1'
                    in  "|" ++ fmap bit (take n [0..]) ++ ">"
          n = fromIntegral $ fromSing (sing :: Sing n)

internalLen :: (KnownNat n, Num a) => Sing n -> a
internalLen = (2 ^) . fromIntegral . fromSing

-- | Number of qubits in the StateVec
dimension :: forall (n :: Nat) a. (KnownNat n, Num a) => StateVec n -> a
dimension _ = fromIntegral $ fromSing (sing :: SNat n)

-- | StateVec's are intialized to |0>
mkStateVec :: forall n . KnownNat n => StateVec n
mkStateVec = UnsafeMkStateVec $ internalLen (sing :: Sing n) |> (1 : repeat 0)

-- | Make a statevector with an explicity specified size.
-- StateVec's are intialized to |0>
mkStateVec' :: Sing n -> StateVec n
mkStateVec' sn = UnsafeMkStateVec $ l |> (1 : repeat 0)
  where l = (2^) . fromIntegral $ fromSing sn

-- | A qubit is just a StateVec 1, initalized to |0>
mkQubit :: StateVec 1
mkQubit = UnsafeMkStateVec $ LA.fromList [1, 0]

normalize :: StateVec n -> StateVec n
normalize (UnsafeMkStateVec zs) = UnsafeMkStateVec $ LA.normalize zs

adjoint :: KnownNat n => StateVec n -> (StateVec n -> C)
adjoint qr = (qr <.>)

-- | The tensor product on elements of our Hilbert space
tensor :: StateVec n -> StateVec m -> StateVec (n + m)
tensor (UnsafeMkStateVec zs) (UnsafeMkStateVec ws) =
  UnsafeMkStateVec $ LA.flatten (zs `LA.outer` ws)

-- | Collapse a StateVec to a state compatable with qubit <n> being measured <bit>.
-- Only occurs physically with measurement, so internal use only.
collapse
  :: forall n . KnownNat n
  => Finite n -> Bit -> StateVec n -> StateVec n
collapse i b (UnsafeMkStateVec zs) = normalize . UnsafeMkStateVec $ zs * mask
  where
    mask   = l |> altseq
    altseq = replicate m ifZero ++ replicate m ifOne ++ altseq
    ifZero = if b == Zero then 1 else 0
    ifOne  = if b == One  then 1 else 0
    l      = internalLen (sing :: Sing n)
    m      = l `quot` (2 ^ (getFinite i + 1))

-- | Preforms a measurement of an induvidual qubit (indexed Finite n)
-- in a quantum register.
measureQubit
  :: (MonadRandom m, KnownNat n)
  => Finite n -> StateT (StateVec n) m Bit
measureQubit i = do
  qr <- get
  r  <- getRandomR (0, 1)
  let qrZero = collapse i Zero qr
      qrOne  = collapse i One qr
      pOne   = realPart $ qrOne <.> qr -- guaranteed to be real, so this
  if r < pOne                          -- is just a type-cast.
    then put qrOne  >> pure One
    else put qrZero >> pure Zero

-- | Preform a measurement on all qubits in a quantum register, returning a
-- computational basis state.
measure
  :: forall m n . (MonadRandom m, KnownNat n)
  => StateT (StateVec n) m CReg
measure = mkCReg <$> traverse measureQubit (take n [0 ..])
  where n = fromIntegral $ fromSing (sing :: Sing n)

max_stars_repo_name: qubitrot/qubism

max_issues_repo_name: qubitrot/qubism

max_forks_repo_name: qubitrot/qubism


text:
{-# OPTIONS_GHC -Wall #-}

-- | Ring classes. A distinguishment is made between Rings and Commutative Rings.
module Plankton.Ring
  ( Semiring
  , Ring
  , CRing
  ) where

import Data.Complex (Complex(..))
import Plankton.Additive
import Plankton.Distribution
import Plankton.Multiplicative
import Protolude (Bool(..), Double, Float, Int, Integer)

-- | Semiring
class (MultiplicativeAssociative a, MultiplicativeUnital a, Distribution a) =>
      Semiring a

instance Semiring Double

instance Semiring Float

instance Semiring Int

instance Semiring Integer

instance Semiring Bool

instance (AdditiveGroup a, Semiring a) => Semiring (Complex a)

-- | Ring
-- a summary of the laws inherited from the ring super-classes
--
-- > zero + a == a
-- > a + zero == a
-- > (a + b) + c == a + (b + c)
-- > a + b == b + a
-- > a - a = zero
-- > negate a = zero - a
-- > negate a + a = zero
-- > a + negate a = zero
-- > one `times` a == a
-- > a `times` one == a
-- > (a `times` b) `times` c == a `times` (b `times` c)
-- > a `times` (b + c) == a `times` b + a `times` c
-- > (a + b) `times` c == a `times` c + b `times` c
-- > a `times` zero == zero
-- > zero `times` a == zero
class ( AdditiveGroup a
      , MultiplicativeAssociative a
      , MultiplicativeUnital a
      , Distribution a
      ) =>
      Ring a

instance Ring Double

instance Ring Float

instance Ring Int

instance Ring Integer

instance (Ring a) => Ring (Complex a)

-- | CRing is a Ring with Multiplicative Commutation.  It arises often due to '*' being defined as a multiplicative commutative operation.
class (Multiplicative a, Ring a) =>
      CRing a

instance CRing Double

instance CRing Float

instance CRing Int

instance CRing Integer

instance (CRing a) => CRing (Complex a)

max_stars_repo_name: chessai/plankton

max_issues_repo_name: chessai/plankton

max_forks_repo_name: chessai/plankton


text:
module ImageProcess where

import Codec.Picture
import Data.Vector.Storable (Vector, length, toList)
import Numeric.LinearAlgebra
--import Data.ByteString

imagePath :: FilePath
imagePath = "./test.jpg"

getDynamicImage :: FilePath -> IO(Maybe DynamicImage)
getDynamicImage filepath = do
  res <- readImage filepath
  case res of
    Left err -> print err >> return Nothing
    Right image -> return (Just image)

readImageInfo :: DynamicImage -> Image PixelRGB8
readImageInfo image =
  let imagePixel8@(Image w h rgbs) = convertRGB8 image in
  imagePixel8

getImageInfo :: FilePath -> IO(Maybe (Image PixelRGB8))
getImageInfo filepath = do
  res <- readImage filepath
  case res of
    Left err -> print err >> return Nothing
    Right image -> return $ Just (readImageInfo image)

getImageWidth :: Image PixelRGB8 -> Int
getImageWidth image =
  let (Image w _ _) = image in
  w

getImageHeight :: Image PixelRGB8 -> Int
getImageHeight image =
  let (Image _ h _) = image in
  h

getImageRGBs :: Image PixelRGB8 -> Data.Vector.Storable.Vector(PixelBaseComponent Pixel8)
getImageRGBs image =
  let (Image _ _ rgbs) = image in
  rgbs

printPixels :: Show a => [a] -> IO()
printPixels [] = return()
printPixels (xr:xg:xb:xs) =
  print [xr, xg, xb] >>
  printPixels xs

genRGBList :: (Show a, Integral a, Fractional b) => [a] -> ([b], [b], [b])
genRGBList [] = ([], [], [])
genRGBList l =
  let func [] rlist glist blist = ((reverse rlist), (reverse glist), (reverse blist))
      func (xr:xg:xb:xs) rlist glist blist = func xs (((fromIntegral xr)/256.0):rlist) (((fromIntegral xg)/256.0):glist) (((fromIntegral xb)/256.0):blist) in
  func l [] [] []

genMatrixFromList :: [R] -> Int -> Int -> Maybe (Matrix R)
genMatrixFromList [] _ _ = Nothing
genMatrixFromList l w h = Just ((w >< h) l)

testImageProcess :: FilePath -> IO()
testImageProcess filepath = do
    imageM <- getImageInfo filepath
    case imageM of
      Nothing -> putStrLn "nothing.."
      Just imagePixel8 ->
        putStrLn "the width of image: " >>
        print (getImageWidth imagePixel8) >>
        putStrLn "the height ofimage: " >>
        print (getImageHeight imagePixel8) >>
        print (Data.Vector.Storable.length (getImageRGBs imagePixel8)) >>
        let rgbs = getImageRGBs imagePixel8 in
        let w = getImageWidth imagePixel8 in
        let h = getImageHeight imagePixel8 in
        let rgbsList = Data.Vector.Storable.toList rgbs in
        let (rlist, glist, blist) = genRGBList rgbsList in
        print (Prelude.length rlist) >>
        print (Prelude.length glist) >>
        print (Prelude.length blist) >>
        let Just rMat = genMatrixFromList rlist w h in
        let Just gMat = genMatrixFromList glist w h in
        let Just bMat = genMatrixFromList blist w h in
        print rMat

{-

decodeImageWidth :: IO(Maybe (Image PixelRGB8)) -> IO(Maybe Int))
decodeImageWidth imageM = do
  maybeImage <- imageM
  case maybeImage of
    Nothing -> return Nothing
    Just (Image w h rgbs) -> return (Just w)

decodeImageHeight :: IO(Maybe(Image PixelRGB8)) -> IO(Maybe Int)
decodeImageHeight imageM = do
  maybeImage <- imageM
  case maybeImage of
    Nothing -> return Nothing
    Just (Image w h rgbs) -> return (Just h)

decodeImageRGBs :: IO(Maybe (Image PixelRGB8)) -> IO(Maybe (Vector(PixelBaseComponent Pixel8)))
decodeImageRGBs imageM = do
  maybeImage <- imageM
  case maybeImage of
    Nothing -> return Nothing
    Just (Image w h rgbs) -> return (Just rgbs)

readRGBsPixel8 :: Maybe DynamicImage -> IO(Maybe (Vector(PixelBaseComponent Pixel8)))
readRGBsPixel8 Nothing = return Nothing
readRGBsPixel8 (Just image) =
  let imagePixel8@(Image _ _ rgbs) = convertRGB8 image in return (Just rgbs)

getRGBs :: FilePath -> IO(Maybe (Vector(PixelBaseComponent Pixel8)))
getRGBs filepath = do
  dynamicImage <- getDynamicImage filepath
  rgbs <- readRGBsPixel8 dynamicImage
  return rgbs

-}
{-}
readRGBsPixel8 :: DynamicImage -> Maybe (Vector (PixelBaseComponent Pixel8))
readRGBsPixel8 (ImageY8 image@(Image _ _ rgbs)) = Just rgbs
readRGBsPixel8 _ = Nothing

readRGBsPixel16 :: DynamicImage -> Maybe (Vector (PixelBaseComponent Pixel16))
readRGBsPixel16 :: (ImageY16 image@(Image _ _ rgbs)) = Just rgbs
readRGBsPixel16 _ = Nothing

readRGBsPixelYF :: DynamicImage -> Maybe (Vector (PixelBaseComponent PixelF))
readRGBsPixelYF :: (ImageF image@(Image _ _ rgbs)) = Just rgbs
readRGBsPixelYF _ = Nothing

readRGBsPixelYA8 :: DynamicImage -> Maybe (Vector (PixelBaseComponent PixelYA8))
readRGBsPixelYA8 :: (ImageYA8 image@(Image _ _ rgbs)) = Just rgbs
readRGBsPixelYA8 _ = Nothing

readRGBsPixelYA16 :: DynamicImage -> Maybe (Vector (PixelBaseComponent PixelYA16))
readRGBsPixelYA16 :: (ImageYA16 image@(Image _ _ rgbs)) = Just rgbs
readRGBsPixelYA16 _ = Nothing

getRGBsImage :: forall a. Pixel a => DynamicImage -> Vector (PixelBaseComponent a)
getRGBsImage (ImageY8 image@(Image _ _ rgbs@(MVector (Word8 _)))) = rgbs
getRGBsImage (ImageY16 image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageYF image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageYA8 image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageRGB8 image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageRGB16 image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageRGBF image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageRGBA8 image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageRGBA16 image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageYCbCr8 image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageCMYK8 image@(Image _ _ rgbs)) = rgbs
getRGBsImage (ImageCMYK16 image@(Image _ _ rgbs)) = rgbs

readRGBs :: forall a. Pixel a => FilePath  -> IO(Maybe (Vector (PixelBaseComponent a)))
readRGBs filepath = do
  res <- getDynamicImage filepath
  case res of
    Nothing -> print "we got nothing!" >> return Nothing
    Just image -> return (Just (getRGBsImage image))

-}
{-
getRGBs :: Image -> Vector(PixelBaseComponent Word8)
getRGBs (Image _ _ rgbs) = rgbs
-}

justH :: Maybe Char
justH = do
  (x:xs) <- Just "Hello"
  return x

max_stars_repo_name: neutronest/vortex

max_issues_repo_name: neutronest/vortex

max_forks_repo_name: neutronest/vortex


text:
-- |
-- Module      : Test.LeanCheck.Function.ShowFunction
-- Copyright   : (c) 2015-2018 Rudy Matela
-- License     : 3-Clause BSD  (see the file LICENSE)
-- Maintainer  : Rudy Matela <rudy@matela.com.br>
--
-- This module is part of LeanCheck,
-- a simple enumerative property-based testing library.
--
-- This module exports the 'ShowFunction' typeclass,
-- its instances and related functions.
--
-- Using this module, it is possible to implement
-- a 'Show' instance for functions:
--
-- > import Test.LeanCheck.ShowFunction
-- > instance (Show a, Listable a, ShowFunction b) => Show (a->b) where
-- >   show = showFunction 8
--
-- This shows functions as a case pattern with up to 8 cases.
--
-- It will only work for functions whose ultimate return value is an instance
-- of 'ShowFunction'.  This module provides instances for most standard data
-- types ('Int', 'Bool', 'Maybe', ...).  Please see the 'ShowFunction'
-- typeclass documentation for how to declare istances for user-defined data
-- types.
--
-- The modules
-- "Test.LeanCheck.Function"
-- and
-- "Test.LeanCheck.Function.Show"
-- exports an instance like the one above.
{-# LANGUAGE CPP #-}
module Test.LeanCheck.Function.ShowFunction
  (
  -- * Showing functions
    showFunction
  , showFunctionLine

  -- * Support for user-defined algebraic datatypes on return values
  , ShowFunction (..)
  , bindtiersShow

  -- * Listing functional bindings
  , Binding
  , bindings

  -- * Pipeline for explaining, describing and clarifying bindings
  , explainedBindings
  , describedBindings
  , clarifiedBindings

  -- * Re-exports
  , Listable
  )
where

import Test.LeanCheck.Core
import Test.LeanCheck.Error (errorToNothing)
import Test.LeanCheck.Utils.Types
import Test.LeanCheck.Stats (classifyOn)
import Data.List (intercalate, sortBy)
import Data.Maybe
import Data.Function (on)
import Data.Word
import Data.Int
import Data.Ratio
import Data.Complex
import Data.Char (GeneralCategory)
import System.Exit (ExitCode)
import System.IO (IOMode, BufferMode, SeekMode)
import Foreign.C

-- | A functional binding in a showable format.
--   Argument values are represented as a list of strings.
--   The result value is represented by 'Just' a 'String' when defined
--   or by 'Nothing' when 'undefined'.
type Binding = ([String], Maybe String)

-- | 'ShowFunction' values are those for which
--   we can return a list of functional bindings.
--
-- Instances for 'show'able algebraic datatypes are defined using
-- 'bindtiersShow':
--
-- > instance ShowFunction Ty where bindtiers = bindtiersShow
class ShowFunction a where
  bindtiers :: a -> [[Binding]]

-- | Given a 'ShowFunction' value, return a list of 'Binding's.
--   If the domain of the given argument function is infinite,
--   the resulting list is infinite.
--
-- Some examples follow.  These are used as running examples in the definition
-- of 'explainedBindings', 'describedBindings' and 'clarifiedBindings'.
--
-- * Defined return values are represented as 'Just' 'String's:
--
--     > > bindings True
--     > [([],Just "True")]
--
-- * Undefined return values are represented as @Nothing@:
--
--     > > bindings undefined
--     > [([],Nothing)]
--
-- * Infinite domains result in an infinite bindings list:
--
--     > > bindings (id::Int->Int)
--     > [ (["0"], Just "0")
--     > , (["1"], Just "1")
--     > , (["-1"], Just "-1")
--     > , ...
--     > ]
--
-- * Finite domains result in a finite bindings list:
--
--     > > bindings (&&)
--     > [ (["False","False"], Just "False")
--     > , (["False","True"], Just "False")
--     > , (["True","False"], Just "False")
--     > , (["True","True"], Just "True")
--     > ]
--
--     > > bindings (||)
--     > [ (["False","False"], Just "False")
--     > , (["False","True"], Just "True")
--     > , (["True","False"], Just "True")
--     > , (["True","True"], Just "True")
--     > ]
--
-- * Even very simple functions are represented by an infinite list of bindings:
--
--     > > bindings (== 0)
--     > [ (["0"], Just "True")
--     > , (["1"], Just "False")
--     > , (["-1"], Just "False")
--     > , ...
--     > ]
--
--     > > bindings (== 1)
--     > [ (["0"], Just "False")
--     > , (["1"], Just "True")
--     > , (["-1"], Just "False")
--     > , ...
--     > ]
--
-- * Ignored arguments are still listed:
--
--     > > bindings ((\_ y -> y == 1) :: Int -> Int -> Bool)
--     > [ (["0","0"], Just "False")
--     > , (["0","1"], Just "True")
--     > , (["1","0"], Just "False")
--     > , ...
--     > ]
--
-- * Again, undefined values are represented as 'Nothing'.
--   Here, the 'head' of an empty list is undefined:
--
--     > > bindings (head :: [Int] -> Int)
--     > [ (["[]"], Nothing)
--     > , (["[0]"], Just "0")
--     > , (["[0,0]"], Just "0")
--     > , (["[1]"], Just "1")
--     > , ...
--     > ]
bindings :: ShowFunction a => a -> [Binding]
bindings = concat . bindtiers

-- | A drop-in implementation of 'bindtiers' for 'show'able types.
--
-- Define instances for 'show'able algebraic datatypes as:
--
-- > instance ShowFunction Ty where bindtiers = bindtiersShow
bindtiersShow :: Show a => a -> [[Binding]]
bindtiersShow x = [[([],errorToNothing $ show x)]]

instance ShowFunction ()       where bindtiers = bindtiersShow
instance ShowFunction Bool     where bindtiers = bindtiersShow
instance ShowFunction Int      where bindtiers = bindtiersShow
instance ShowFunction Word     where bindtiers = bindtiersShow
instance ShowFunction Integer  where bindtiers = bindtiersShow
instance ShowFunction Char     where bindtiers = bindtiersShow
instance ShowFunction Float    where bindtiers = bindtiersShow
instance ShowFunction Double   where bindtiers = bindtiersShow
instance ShowFunction Ordering where bindtiers = bindtiersShow
instance Show a => ShowFunction [a]       where bindtiers = bindtiersShow
instance Show a => ShowFunction (Maybe a) where bindtiers = bindtiersShow
instance (Show a, Show b) => ShowFunction (Either a b) where bindtiers = bindtiersShow
instance (Show a, Show b) => ShowFunction (a,b) where bindtiers = bindtiersShow

-- instance for functional value type --
instance (Show a, Listable a, ShowFunction b) => ShowFunction (a->b) where
  bindtiers f = concatMapT bindtiersFor tiers
    where bindtiersFor x = mapFst (show x:) `mapT` bindtiers (f x)
          mapFst f (x,y) = (f x, y)

paren :: String -> String
paren s = "(" ++ s ++ ")"

showTuple :: [String] -> String
showTuple [x]  =  x
showTuple xs | all (== "_") xs  =  "_"
             | otherwise        =  paren $ intercalate "," xs

showBindings :: [Binding] -> [String]
showBindings bs = [ showTuple as ++ " -> " ++ r | (as, Just r) <- bs ]

showNBindings :: Bool -> Int -> [Binding] -> [String]
showNBindings infinite n bs' = take n bs ++ ["..." | infinite || length bs > n]
  where
  bs = showBindings bs'

isValue :: ShowFunction a => a -> Bool
isValue f = case bindings f of
              [([],_)] -> True
              _        -> False

showValueOf :: ShowFunction a => a -> String
showValueOf x = case snd . head . bindings $ x of
                  Nothing -> "undefined"
                  Just x' -> x'

-- | Given the number of patterns to show, shows a 'ShowFunction' value.
--
-- > > putStrLn $ showFunction undefined True
-- > True
-- >
-- > > putStrLn $ showFunction 3 (id::Int->Int)
-- > \x -> case x of
-- >       0 -> 0
-- >       1 -> 1
-- >       -1 -> -1
-- >       ...
-- >
-- > > putStrLn $ showFunction 4 (&&)
-- > \x y -> case (x,y) of
-- >         (True,True) -> True
-- >         _ -> False
-- >
--
-- In the examples above, "@...@" should be interpreted literally.
--
-- This can be used as an implementation of 'show' for functions:
--
-- > instance (Show a, Listable a, ShowFunction b) => Show (a->b) where
-- >   show = showFunction 8
--
-- See 'showFunctionLine' for an alternative without line breaks.
showFunction :: ShowFunction a => Int -> a -> String
showFunction n = showFunctionL False (n*n+1) n

-- | Same as 'showFunction', but has no line breaks.
--
-- > > putStrLn $ showFunctionLine 3 (id::Int->Int)
-- > \x -> case x of 0 -> 0; 1 -> 1; -1 -> -1; ...
-- > > putStrLn $ showFunctionLine 3 (&&)
-- > \x y -> case (x,y) of (True,True) -> True; _ -> False
--
-- This can be used as an implementation of 'show' for functions:
--
-- > instance (Show a, Listable a, ShowFunction b) => Show (a->b) where
-- >   show = showFunction 8
showFunctionLine :: ShowFunction a => Int -> a -> String
showFunctionLine n = showFunctionL True (n*n+1) n

-- | isUndefined checks if a function is totally undefined
--   for the given maximum number of values
isUndefined :: ShowFunction a => Int -> a -> Bool
isUndefined m = all (isNothing . snd) . take m . bindings

-- | checks if a function is constant
--   for the given maximum number of values
isConstant :: ShowFunction a => Int -> a -> Bool
isConstant m f = case take m $ bindings f of
                 []          -> False -- uninhabited type?
                 ((_,r'):bs) -> all (\(_,r) -> r == r') bs

-- | shows a constant function
showConstant :: ShowFunction a => Int -> a -> String
showConstant m f = "\\" ++ unwords vs ++ " -> " ++ fromMaybe "undefined" r
  where
  (as,r) = head $ bindings f
  vs = replicate (length as) "_"

-- The first boolean parameter tells if we are showing
-- the function on a single line
showFunctionL :: ShowFunction a => Bool -> Int -> Int -> a -> String
showFunctionL singleLine m n f | isValue f = showValueOf f
showFunctionL singleLine m n f | isConstant m f = showConstant m f
--showFunctionL singleLine m n f | canName m f = showName m f
showFunctionL singleLine m n f | otherwise = lambdaPat ++ caseExp
  where
    lambdaPat = "\\" ++ unwords vs ++ " -> "
    casePat = "case " ++ showTuple (filter (/= "_") vs) ++ " of"
    (vs, bindings) = clarifiedBindings m n f
    bs = showNBindings (length bindings >= m) n bindings
    sep | singleLine = " "
        | otherwise = "\n"
    cases | singleLine = intercalate "; " bs
          | otherwise  = unlines
                       $ (replicate (length lambdaPat) ' ' ++) `map` bs
    caseExp = if isUndefined m f
                then "undefined"
                else casePat ++ sep ++ cases

-- instances for further tuple arities --
instance (Show a, Show b, Show c)
      => ShowFunction (a,b,c) where bindtiers = bindtiersShow
instance (Show a, Show b, Show c, Show d)
      => ShowFunction (a,b,c,d) where bindtiers = bindtiersShow
instance (Show a, Show b, Show c, Show d, Show e)
      => ShowFunction (a,b,c,d,e) where bindtiers = bindtiersShow
instance (Show a, Show b, Show c, Show d, Show e, Show f)
      => ShowFunction (a,b,c,d,e,f) where bindtiers = bindtiersShow
instance (Show a, Show b, Show c, Show d, Show e, Show f, Show g)
      => ShowFunction (a,b,c,d,e,f,g) where bindtiers = bindtiersShow
instance (Show a, Show b, Show c, Show d, Show e, Show f, Show g, Show h)
      => ShowFunction (a,b,c,d,e,f,g,h) where bindtiers = bindtiersShow
instance ( Show a, Show b, Show c, Show d
         , Show e, Show f, Show g, Show h
         , Show i )
      => ShowFunction (a,b,c,d,e,f,g,h,i) where bindtiers = bindtiersShow
instance ( Show a, Show b, Show c, Show d
         , Show e, Show f, Show g, Show h
         , Show i, Show j )
      => ShowFunction (a,b,c,d,e,f,g,h,i,j) where bindtiers = bindtiersShow
instance ( Show a, Show b, Show c, Show d
         , Show e, Show f, Show g, Show h
         , Show i, Show j, Show k )
      => ShowFunction (a,b,c,d,e,f,g,h,i,j,k) where bindtiers = bindtiersShow
instance ( Show a, Show b, Show c, Show d
         , Show e, Show f, Show g, Show h
         , Show i, Show j, Show k, Show l )
      => ShowFunction (a,b,c,d,e,f,g,h,i,j,k,l) where bindtiers = bindtiersShow

-- Data.Ratio
instance (Integral a, Show a) => ShowFunction (Ratio a) where bindtiers = bindtiersShow

-- Data.Complex
instance (RealFloat a, Show a) => ShowFunction (Complex a) where bindtiers = bindtiersShow

-- instance for types from Data.Int and Data.Word
instance ShowFunction Int8  where bindtiers = bindtiersShow
instance ShowFunction Int16 where bindtiers = bindtiersShow
instance ShowFunction Int32 where bindtiers = bindtiersShow
instance ShowFunction Int64 where bindtiers = bindtiersShow
instance ShowFunction Word8  where bindtiers = bindtiersShow
instance ShowFunction Word16 where bindtiers = bindtiersShow
instance ShowFunction Word32 where bindtiers = bindtiersShow
instance ShowFunction Word64 where bindtiers = bindtiersShow

-- instance for types from Test.LeanCheck.Utils.Types
instance ShowFunction Nat   where bindtiers = bindtiersShow
instance ShowFunction Nat1  where bindtiers = bindtiersShow
instance ShowFunction Nat2  where bindtiers = bindtiersShow
instance ShowFunction Nat3  where bindtiers = bindtiersShow
instance ShowFunction Nat4  where bindtiers = bindtiersShow
instance ShowFunction Nat5  where bindtiers = bindtiersShow
instance ShowFunction Nat6  where bindtiers = bindtiersShow
instance ShowFunction Nat7  where bindtiers = bindtiersShow
instance ShowFunction Int1  where bindtiers = bindtiersShow
instance ShowFunction Int2  where bindtiers = bindtiersShow
instance ShowFunction Int3  where bindtiers = bindtiersShow
instance ShowFunction Int4  where bindtiers = bindtiersShow
instance ShowFunction Word1 where bindtiers = bindtiersShow
instance ShowFunction Word2 where bindtiers = bindtiersShow
instance ShowFunction Word3 where bindtiers = bindtiersShow
instance ShowFunction Word4 where bindtiers = bindtiersShow

instance ShowFunction Natural where bindtiers = bindtiersShow

instance ShowFunction Letter    where bindtiers = bindtiersShow
instance ShowFunction AlphaNum  where bindtiers = bindtiersShow
instance ShowFunction Digit     where bindtiers = bindtiersShow
instance ShowFunction Alpha     where bindtiers = bindtiersShow
instance ShowFunction Upper     where bindtiers = bindtiersShow
instance ShowFunction Lower     where bindtiers = bindtiersShow
instance ShowFunction Space     where bindtiers = bindtiersShow

instance ShowFunction Spaces    where bindtiers = bindtiersShow
instance ShowFunction Lowers    where bindtiers = bindtiersShow
instance ShowFunction Uppers    where bindtiers = bindtiersShow
instance ShowFunction Alphas    where bindtiers = bindtiersShow
instance ShowFunction Digits    where bindtiers = bindtiersShow
instance ShowFunction AlphaNums where bindtiers = bindtiersShow
instance ShowFunction Letters   where bindtiers = bindtiersShow

instance Show a => ShowFunction (X a) where bindtiers = bindtiersShow
instance Show a => ShowFunction (Xs a) where bindtiers = bindtiersShow
instance Show a => ShowFunction (Set a) where bindtiers = bindtiersShow
instance Show a => ShowFunction (Bag a) where bindtiers = bindtiersShow
instance Show a => ShowFunction (NoDup a) where bindtiers = bindtiersShow
instance (Show a, Show b) => ShowFunction (Map a b) where bindtiers = bindtiersShow

-- misc instances
instance ShowFunction ExitCode   where bindtiers = bindtiersShow
instance ShowFunction SeekMode   where bindtiers = bindtiersShow
instance ShowFunction IOMode     where bindtiers = bindtiersShow
instance ShowFunction BufferMode where bindtiers = bindtiersShow
instance ShowFunction GeneralCategory where bindtiers = bindtiersShow

-- instances for Foreign.C types
instance ShowFunction CChar      where bindtiers = bindtiersShow
instance ShowFunction CSChar     where bindtiers = bindtiersShow
instance ShowFunction CUChar     where bindtiers = bindtiersShow
instance ShowFunction CShort     where bindtiers = bindtiersShow
instance ShowFunction CUShort    where bindtiers = bindtiersShow
instance ShowFunction CInt       where bindtiers = bindtiersShow
instance ShowFunction CUInt      where bindtiers = bindtiersShow
instance ShowFunction CLong      where bindtiers = bindtiersShow
instance ShowFunction CULong     where bindtiers = bindtiersShow
instance ShowFunction CPtrdiff   where bindtiers = bindtiersShow
instance ShowFunction CSize      where bindtiers = bindtiersShow
instance ShowFunction CWchar     where bindtiers = bindtiersShow
instance ShowFunction CSigAtomic where bindtiers = bindtiersShow
instance ShowFunction CLLong     where bindtiers = bindtiersShow
instance ShowFunction CULLong    where bindtiers = bindtiersShow
instance ShowFunction CIntPtr    where bindtiers = bindtiersShow
instance ShowFunction CUIntPtr   where bindtiers = bindtiersShow
instance ShowFunction CIntMax    where bindtiers = bindtiersShow
instance ShowFunction CUIntMax   where bindtiers = bindtiersShow
instance ShowFunction CClock     where bindtiers = bindtiersShow
instance ShowFunction CTime      where bindtiers = bindtiersShow
instance ShowFunction CFloat     where bindtiers = bindtiersShow
instance ShowFunction CDouble    where bindtiers = bindtiersShow
#if __GLASGOW_HASKELL__ >= 802
instance ShowFunction CBool      where bindtiers = bindtiersShow
#endif
#if __GLASGOW_HASKELL__
instance ShowFunction CUSeconds  where bindtiers = bindtiersShow
instance ShowFunction CSUSeconds where bindtiers = bindtiersShow
#endif

-- | Returns a set of variables and a set of bindings
--   describing how a function works.
--
-- Some argument values are generalized to "@_@" when possible.
-- If one of the function arguments is not used altogether, it is ommited in
-- the set of bindings and appears as "_" in the variables list.
--
-- This is the /last/ function in the clarification pipeline.
--
-- It takes two integer arguments:
--
-- 1. @m@: the maximum number of cases considered for computing the description;
-- 2. @n@: the maximum number of cases in the actual description.
--
-- As a general rule of thumb, set @m=n*n+1@.
--
-- Some examples follow:
--
-- * When all arguments are used, the result is the same as 'describedBindings':
--
--     > > clarifiedBindings 100 10 (==1)
--     > ( ["x"], [ (["1"],Just "True"),
--     >          , (["_"],Just "False") ] )
--
-- * When some arguments are unused, they are omitted in the list of bindings and
--   appear as @"_"@ in the list of variables.
--
--     > > clarifiedBindings 100 10 (\_ y -> y == 1)
--     > ( ["_", "y"], [ (["1"],Just "True")
--     >               , (["_"],Just "False") ] )
clarifiedBindings :: ShowFunction a => Int -> Int -> a -> ([String],[Binding])
clarifiedBindings m n = clarifyBindings . describedBindings m n

clarifyBindings :: [Binding] -> ([String],[Binding])
clarifyBindings bs  =  (varnamesByUsage used, map (mapFst $ select used) bs)
  where
  mapFst f (x,y) = (f x, y)
  used = usedArgs bs

varnamesByUsage :: [Bool] -> [String]
varnamesByUsage = zipWith used varnames
  where
  used s False = "_"
  used s True  = s
  varnames = ["x","y","z","w"] ++ map (++"'") varnames

usedArgs :: [Binding] -> [Bool]
usedArgs = foldr1 (zipWith (||))
         . map (map (/= "_") . fst)

-- | Returns a set of bindings describing how a function works.
-- Some argument values are generalized to "@_@" when possible.
-- It takes two integer arguments:
--
-- 1. @m@: the maximum number of cases considered for computing description;
-- 2. @n@: the maximum number of cases in the actual description.
--
-- As a general rule of thumb, set @m=n*n+1@.
--
-- This is the /second/ function in the clarification pipeline.
--
-- This function processes the result of 'explainedBindings'
-- to sometimes return shorter descriptions.
-- It chooses the shortest of the following (in order):
--
-- * regular unexplained-undescribed 'bindings';
-- * regular 'explainedBindings';
-- * 'explainedBindings' with least occurring cases generalized first;
--
-- Here are some examples:
--
-- * Sometimes the result is the same as 'explainedBindings':
--
--     > > describedBindings 100 10 (||)
--     > [ (["False","False"],Just "False")
--     > , (["_","_"],Just "True") ]
--
--     > > describedBindings 100 10 (==0)
--     > [ (["0"],Just "True")
--     > , (["_"],Just "False") ]
--
-- * but sometimes it is shorter because we consider generalizing least
--   occurring cases first:
--
--     > > describedBindings 100 10 (&&)
--     > [ ( ["True","True"],Just "True")
--     > , ( ["_","_"],Just "False") ]
--
--     > > describedBindings 100 10 (==1)
--     > [ (["1"],Just "True"),
--     > , (["_"],Just "False") ]
--
--     > > describedBindings 100 10 (\_ y -> y == 1)
--     > [ (["_","1"],Just "True")
--     > , (["_","_"],Just "False") ]
describedBindings :: ShowFunction a => Int -> Int -> a -> [Binding]
describedBindings m n f
  | length bs1 <= n  =  bs1
  | otherwise        =  bs0
  where
  bs0  =  take m $ bindings f
  bs1  =  describeBindings bs0

describeBindings :: [Binding] -> [Binding]
describeBindings bs = head $ sortOn length $
  [ bs
  , explainBindings bs
  , explainBindings . concat . sortOn length $ classifyOn snd bs
  ]

-- | Returns a set of bindings explaining how a function works.
--   Some argument values are generalized to "@_@" when possible.
--   It takes as argument the maximum number of cases
--   considered for computing the explanation.
--
-- A measure of success in this generalization process is if this function
-- returns less values than the asked maximum number of cases.
--
-- This is the /first/ function in the clarification pipeline.
--
-- * In some cases, 'bindings' cannot be "explained"
--   an almost unchanged result of 'bindings' is returned
--   with the last binding having variables replaced by "@_@":
--
--     > > explainedBindings 4 (id::Int->Int)
--     > [ (["0"],Just "0")
--     > , (["1"],Just "1")
--     > , (["-1"],Just "-1")
--     > , (["_"],Just "2") ]
--
-- * When possible, some cases are generalized using @_@:
--
--     > > explainedBindings 10 (||)
--     > [ (["False","False"],Just "False")
--     > , (["_","_"],Just "True") ]
--
--     but the resulting "explanation" might not be the shortest possible
--     (cf. 'describedBindings'):
--
--     > > explainedBindings 10 (&&)
--     > [ ( ["False","_"],Just "False")
--     > , (["_","False"],Just "False")
--     > , (["_","_"],Just "True") ]
--
-- * Generalization works for infinite domains (heuristically):
--
--     > > explainedBindings 10 (==0)
--     > [ (["0"],Just "True")
--     > , (["_"],Just "False") ]
--
-- * Generalization for each item is processed in the order they are generated by 'bindings'
--   hence explanations are not always the shortest possible (cf. 'describedBindings').
--   In the following examples, the first case is redundant.
--
--     > > explainedBindings 10 (==1)
--     > [ (["0"],Just "False")
--     > , (["1"],Just "True"),
--     > , (["_"],Just "False") ]
--
--     > > explainedBindings 10 (\_ y -> y == 1)
--     > [ (["_","0"],Just "False")
--     > , (["_","1"],Just "True")
--     > , (["_","_"],Just "False") ]
explainedBindings :: ShowFunction a => Int -> a -> [Binding]
explainedBindings m = explainBindings . take m . bindings

explainBindings :: [Binding] -> [Binding]
explainBindings = explain []
  where
  explain :: [Binding] -> [Binding] -> [Binding]
  explain bs' []           =  reverse bs'
  explain bs' ((as,r):bs)  =  explain (bs''++bs') [b | b <- bs, none (b <~~) bs'']
    where
    bs'' = discardLater (<~~)
         [ (gas,r) | gas <- generalizations as
                   , and [r' == r | (as',r') <- bs, as' <~ gas] ]

generalizations :: [String] -> [[String]]
generalizations []     = [[]]
generalizations (v:vs) = map ("_":) gvs ++ map (v:) gvs
  where
  gvs = generalizations vs

-- | Should be read as "is generalized by":
--
-- > > ["1","2","3"] <~ ["_","_","_"]
-- > True
-- > > ["_","_","_"] <~ ["1","2","3"]
-- > False
-- > > ["1","3"] <~ ["_","3"]
-- > True
-- > > ["_","3"] <~ ["_","4"]
-- > False
(<~) :: [String] -> [String] -> Bool
[]     <~ []       =  True
(v:vs) <~ ("_":ws) =  vs <~ ws
(v:vs) <~ (w:ws)   =  v == w && vs <~ ws
_      <~ _        =  False

-- | Should be read as "is generalized by".
(<~~) :: Binding -> Binding -> Bool
(as,r) <~~ (as',r') = as <~ as' && r == r'

-- general auxiliary functions

discard :: (a -> Bool) -> [a] -> [a]
discard p = filter (not . p)

discardLater :: (a -> a -> Bool) -> [a] -> [a]
discardLater (?>) = dl
  where
  dl []     = []
  dl (x:xs) = x : discard (?> x) (dl xs)

none :: (a -> Bool) -> [a] -> Bool
none p = not . any p

-- sortOn is only available on GHC > 7.8
sortOn :: Ord b => (a -> b) -> [a] -> [a]
sortOn f = sortBy (compare `on` f)

select :: [Bool] -> [a] -> [a]
select [] _ = []
select _ [] = []
select (p:ps) (x:xs) = if p then x : xs' else xs' where xs' = select ps xs

max_stars_repo_name: jwaldmann/leancheck

max_issues_repo_name: jwaldmann/leancheck

max_forks_repo_name: jwaldmann/leancheck


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP          #-}
{-# LANGUAGE ViewPatterns #-}
-- | Tests for Statistics.Math
module Tests.SpecFunctions (
  tests
  ) where

import Control.Monad
import Data.List
import Data.Maybe
import qualified Data.Vector as V
import           Data.Vector   ((!))
import qualified Data.Vector.Unboxed as U

import Test.QuickCheck  hiding (choose,within)
import Test.Tasty
import Test.Tasty.QuickCheck   (testProperty)
import Test.Tasty.HUnit

import Tests.Helpers
import Tests.SpecFunctions.Tables
import Numeric.SpecFunctions
import Numeric.SpecFunctions.Internal   (factorialTable)
import Numeric.MathFunctions.Comparison (within,relativeError,ulpDistance)
import Numeric.MathFunctions.Constants  (m_epsilon,m_tiny)

erfTol,erfcTol,erfcLargeTol :: Int
#if USE_SYSTEM_ERF && !defined(__GHCJS__)
erfTol       = 1
erfcTol      = 2
erfcLargeTol = 2
#else
erfTol       = 4
erfcTol      = 4
erfcLargeTol = 64
#endif

isGHCJS :: Bool
#if defined(__GHCJS__)
isGHCJS = True
#else
isGHCJS = False
#endif

tests :: TestTree
tests = testGroup "Special functions"
  [ testGroup "erf"
    [ -- implementation from numerical recipes loses presision for
      -- large arguments
      testCase "erfc table" $
        forTable "tests/tables/erfc.dat" $ \[x, exact] ->
          checkTabularPure erfcTol (show x) exact (erfc x)
    , testCase "erfc table [large]" $
        forTable "tests/tables/erfc-large.dat" $ \[x, exact] ->
          checkTabularPure erfcLargeTol (show x) exact (erfc x)
      --
    , testCase "erf table" $
        forTable "tests/tables/erf.dat" $ \[x, exact] -> do
          checkTabularPure erfTol (show x) exact (erf x)
    , testProperty "id = erfc . invErfc" invErfcIsInverse
    , testProperty "id = invErfc . erfc" invErfcIsInverse2
    , testProperty "invErf  = erf^-1"    invErfIsInverse
    ]
  --
  , testGroup "log1p & Co"
    [ testCase "expm1 table" $
        forTable "tests/tables/expm1.dat" $ \[x, exact] ->
          checkTabularPure 2 (show x) exact (expm1 x)
    , testCase "log1p table" $
        forTable "tests/tables/log1p.dat" $ \[x, exact] ->
          checkTabularPure 1 (show x) exact (log1p x)
    ]
  ----------------
  , testGroup "gamma function"
    [ testCase "logGamma table [fractional points" $
        forTable "tests/tables/loggamma.dat" $ \[x, exact] -> do
          checkTabularPure 2 (show x) exact (logGamma x)
    , testProperty "Gamma(x+1) = x*Gamma(x)" $ gammaReccurence
    , testCase     "logGamma is expected to be precise at 1e-15 level" $
        forM_ [3..10000::Int] $ \n -> do
          let exact = logFactorial (n-1)
              val   = logGamma (fromIntegral n)
          checkTabular 8 (show n) exact val
    ]
  ----------------
  , testGroup "incomplete gamma"
    [ testCase "incompleteGamma table" $
        forTable "tests/tables/igamma.dat" $ \[a,x,exact] -> do
          let err | a < 10    = 16
                  | a <= 101  = if isGHCJS then 64 else 32
                  | a == 201  = 200
                  | otherwise = 32
          checkTabularPure err (show (a,x)) exact (incompleteGamma a x)
    , testProperty "incomplete gamma - increases" $
        \(abs -> s) (abs -> x) (abs -> y) -> s > 0 ==> monotonicallyIncreases (incompleteGamma s) x y
    , testProperty "0 <= gamma <= 1"               incompleteGammaInRange
    , testProperty "gamma(1,x) = 1 - exp(-x)"      incompleteGammaAt1Check
    , testProperty "invIncompleteGamma = gamma^-1" invIGammaIsInverse
    ]
  ----------------
  , testGroup "beta function"
    [ testCase "logBeta table" $
        forTable "tests/tables/logbeta.dat" $ \[p,q,exact] ->
          let errEst
                -- For Stirling approx. errors are very good
                | b > 10          = 2
                -- Partial Stirling approx
                | a > 10 = case () of
                    _| b >= 1    -> 4
                     | otherwise -> 2 * est
                -- sum of logGamma
                | otherwise = case () of
                    _| a <= 1 && b <= 1 -> 8
                     | a >= 1 && b >= 1 -> 8
                     | otherwise        -> 2 * est
                where
                  a = max p q
                  b = min p q
                  --
                  est = ceiling
                      $ abs (logGamma a) + abs (logGamma b) + abs (logGamma (a + b))
                      / abs (logBeta a b)
          in checkTabularPure errEst (show (p,q)) exact (logBeta p q)
    , testCase "logBeta factorial" betaFactorial
    , testProperty "beta(1,p) = 1/p"   beta1p
    -- , testProperty "beta recurrence"   betaRecurrence
    ]
  ----------------
  , testGroup "incomplete beta"
    [ testCase "incompleteBeta table" $
        forM_ tableIncompleteBeta $ \(p,q,x,exact) ->
          checkTabular 64 (show (x,p,q)) (incompleteBeta p q x) exact
    , testCase "incompleteBeta table with p > 3000 and q > 3000" $
        forM_ tableIncompleteBetaP3000 $ \(x,p,q,exact) ->
          checkTabular 7000 (show (x,p,q)) (incompleteBeta p q x) exact
    --
    , testProperty "0 <= I[B] <= 1" incompleteBetaInRange
    , testProperty "ibeta symmetry" incompleteBetaSymmetry
    -- XXX FIXME DISABLED due to failures
    -- , testProperty "invIncompleteBeta  = B^-1" $ invIBetaIsInverse
    ]
  ----------------
  , testGroup "digamma"
    [ testAssertion "digamma is expected to be precise at 1e-14 [integers]"
        $ digammaTestIntegers 1e-14
      -- Relative precision is lost when digamma(x) ≈ 0
    , testCase "digamma is expected to be precise at 1e-12" $
      forTable "tests/tables/digamma.dat" $ \[x, exact] ->
        checkTabularPure 2048
          (show x) (digamma x) exact
    ]
  ----------------
  , testGroup "factorial"
    [ testCase "Factorial table" $
      forM_ [0 .. 170] $ \n -> do
        checkTabular 1
          (show n)
          (fromIntegral (factorial' n))
          (factorial (fromIntegral n :: Int))
      --
    , testCase "Log factorial from integer" $
      forM_ [2 .. 170] $ \n -> do
        checkTabular 1
          (show n)
          (log $ fromIntegral $ factorial' n)
          (logFactorial (fromIntegral n :: Int))
    , testAssertion "Factorial table is OK"
    $ U.length factorialTable == 171
    , testCase "Log factorial table" $
      forTable "tests/tables/factorial.dat" $ \[i,exact] ->
        checkTabularPure 3
          (show i) (logFactorial (round i :: Int)) exact
    ]
  ----------------
  , testGroup "combinatorics"
    [ testCase "choose table" $
      forM_ [0 .. 1000] $ \n ->
        forM_ [0 .. n]  $ \k -> do
          checkTabular 2048
            (show (n,k))
            (fromIntegral $ choose' n k)
            (choose (fromInteger n) (fromInteger k))
    --
    , testCase "logChoose == log . choose" $
      forM_ [0 .. 1000] $ \n ->
        forM_ [0 .. n]  $ \k -> do
          checkTabular 2
            (show (n,k))
            (log $ choose n k)
            (logChoose n k)
    ]
    ----------------------------------------------------------------
    -- Self tests
  , testGroup "self-test"
    [ testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0
    ]
  ]

----------------------------------------------------------------
-- efr tests
----------------------------------------------------------------

roundtrip_erfc_invErfc,
  roundtrip_invErfc_erfc,
  roundtrip_erf_invErf
  :: (Double,Double)
#if USE_SYSTEM_ERF && !defined(__GHCJS__)
roundtrip_erfc_invErfc = (2,2)
roundtrip_invErfc_erfc = (2,2)
roundtrip_erf_invErf   = (1,1)
#else
roundtrip_erfc_invErfc = (2,8)
roundtrip_invErfc_erfc = (8,4)
roundtrip_erf_invErf   = (128,128)
#endif

-- id ≈ erfc . invErfc
invErfcIsInverse :: Double -> Property
invErfcIsInverse ((*2) . range01 -> x)
  = (not $ isInfinite x) ==>
  ( counterexample ("x        = " ++ show x )
  $ counterexample ("y        = " ++ show y )
  $ counterexample ("x'       = " ++ show x')
  $ counterexample ("calc.err = " ++ show (delta, delta-e'))
  $ counterexample ("ulps     = " ++ show (ulpDistance x x'))
  $ ulpDistance x x' <= round delta
  )
  where
    (e,e') = roundtrip_erfc_invErfc
    delta  = e' + e * abs ( y / x  *  2 / sqrt pi * exp( -y*y ))
    y  = invErfc x
    x' = erfc y

-- id ≈ invErfc . erfc
invErfcIsInverse2 :: Double -> Property
invErfcIsInverse2 x
  = (not $ isInfinite x') ==>
    (y > m_tiny)          ==>
    (x /= 0) ==>
    counterexample ("x        = " ++ show x )
  $ counterexample ("y        = " ++ show y )
  $ counterexample ("x'       = " ++ show x')
  $ counterexample ("calc.err = " ++ show delta)
  $ counterexample ("ulps     = " ++ show (ulpDistance x x'))
  $ ulpDistance x x' <= delta
  where
    (e,e') = roundtrip_invErfc_erfc
    delta  = round
           $ e' + e * abs (y / x  /  (2 / sqrt pi * exp( -x*x )))
    y  = erfc x
    x' = invErfc y

-- id ≈ erf . invErf
invErfIsInverse :: Double -> Property
invErfIsInverse a
  = (x /= 0) ==>
    counterexample ("x        = " ++ show x )
  $ counterexample ("y        = " ++ show y )
  $ counterexample ("x'       = " ++ show x')
  $ counterexample ("calc.err = " ++ show delta)
  $ counterexample ("ulps     = " ++ show (ulpDistance x x'))
  $ ulpDistance x x' <= delta
  where
    (e,e') = roundtrip_erf_invErf
    delta  = round
           $ e + e' * abs (y / x  *  2 / sqrt pi * exp ( -y * y ))
    x  | a < 0     = - range01 a
       | otherwise =   range01 a
    y  = invErf x
    x' = erf y

----------------------------------------------------------------
-- QC tests
----------------------------------------------------------------

-- B(p,q) = (x - 1)!(y-1)! / (x + y - 1)!
betaFactorial :: IO ()
betaFactorial = do
  forM_ prod $ \(p,q,facP,facQ,facProd) -> do
    let exact = fromIntegral (facQ * facP)
              / fromIntegral facProd
    checkTabular 16 (show (p,q))
      (logBeta (fromIntegral p) (fromIntegral q))
      (log exact)
  where
    prod    = [ (p,q,facP,facQ, factorial' (p + q - 1))
              | (p,facP) <- facList
              , (q,facQ) <- facList
              , p + q < 170
              , not (p == 1 && q== 1)
              ]
    facList = [(p,factorial' (p-1)) | p <- [1 .. 170]]

-- B(1,p) = 1/p
beta1p :: Double -> Property
beta1p (abs -> p)
  = p > 2 ==>
    counterexample ("p    = " ++ show p)
  $ counterexample ("logB = " ++ show lb)
  $ counterexample ("err  = " ++ show d)
  $ d <= 24
  where
    lb = logBeta 1 p
    d  = ulpDistance lb (- log p)

{-
-- B(p+1,q) = B(p,q) · p/(p+q)
betaRecurrence :: Double -> Double -> Property
betaRecurrence (abs -> p) (abs -> q)
  = p > 0  &&  q > 0  ==>
    counterexample ("p          = " ++ show p)
  $ counterexample ("q          = " ++ show q)
  $ counterexample ("log B(p,q) = " ++ show (logBeta p q))
  $ counterexample ("log B(p+1,q) = " ++ show (logBeta (p+1) q))
  $ counterexample ("err        = " ++ show d)
  $ d <= 128
  where
    logB  = logBeta p q + log (p / (p + q))
    logB' = logBeta (p + 1) q
    d     = ulpDistance logB logB'
-}

-- Γ(x+1) = x·Γ(x)
gammaReccurence :: Double -> Property
gammaReccurence x
  = x > 0  ==>  err < errEst
    where
      g1     = logGamma x
      g2     = logGamma (x+1)
      err    = abs (g2 - g1 - log x)
      -- logGamma apparently is not as precise for small x. See #59 for details
      errEst = max 1e-14
             $ 2 * m_epsilon * sum (map abs [ g1 , g2 , log x ])

-- γ(s,x) is in [0,1] range
incompleteGammaInRange :: Double -> Double -> Property
incompleteGammaInRange (abs -> s) (abs -> x) =
  x >= 0 && s > 0  ==> let i = incompleteGamma s x in i >= 0 && i <= 1

-- γ(1,x) = 1 - exp(-x)
-- Since Γ(1) = 1 normalization doesn't make any difference
incompleteGammaAt1Check :: Double -> Bool
incompleteGammaAt1Check (abs -> x) =
  ulpDistance (incompleteGamma 1 x) (-expm1(-x)) < 16

-- invIncompleteGamma is inverse of incompleteGamma
invIGammaIsInverse :: Double -> Double -> Property
invIGammaIsInverse (abs -> a) (range01 -> p) =
  a > m_tiny && p > m_tiny && p < 1 && x > m_tiny  ==>
    ( counterexample ("a    = " ++ show a )
    $ counterexample ("p    = " ++ show p )
    $ counterexample ("x    = " ++ show x )
    $ counterexample ("p'   = " ++ show p')
    $ counterexample ("err  = " ++ show (ulpDistance p p'))
    $ counterexample ("est  = " ++ show est)
    $ ulpDistance p p' <= est
    )
  where
    x  = invIncompleteGamma a p
    f' = exp ( log x * (a-1) - x - logGamma a)
    p' = incompleteGamma    a x
    -- FIXME: Test should be rechecked when #42 is fixed
    (e,e') = (32,32)
    est    = round
           $ e' + e * abs (x / p * f')

-- I(x;p,q) is in [0,1] range
incompleteBetaInRange :: Double -> Double -> Double -> Property
incompleteBetaInRange (abs -> p) (abs -> q) (range01 -> x) =
  p > 0 && q > 0  ==> let i = incompleteBeta p q x in i >= 0 && i <= 1

-- I(0.5; p,p) = 0.5
incompleteBetaSymmetry :: Double -> Property
incompleteBetaSymmetry (abs -> p) =
    p > 0 ==>
    counterexample ("p   = " ++ show p)
  $ counterexample ("ib  = " ++ show ib)
  $ counterexample ("err = " ++ show d)
  $ counterexample ("est = " ++ show est)
  $ d <= est
  where
    est | p < 1     = 80
        | p < 10    = 200
        | otherwise = round $ 6 * p
    d  = ulpDistance ib 0.5
    ib = incompleteBeta p p 0.5

-- invIncompleteBeta is inverse of incompleteBeta
invIBetaIsInverse :: Double -> Double -> Double -> Property
invIBetaIsInverse (abs -> p) (abs -> q) (range01 -> x) =
  p > 0 && q > 0  ==> ( counterexample ("p   = " ++ show p )
                      $ counterexample ("q   = " ++ show q )
                      $ counterexample ("x   = " ++ show x )
                      $ counterexample ("x'  = " ++ show x')
                      $ counterexample ("a   = " ++ show a)
                      $ counterexample ("err = " ++ (show $ abs $ (x - x') / x))
                      $ abs (x - x') <= 1e-12
                      )
  where
    x' = incompleteBeta    p q a
    a  = invIncompleteBeta p q x

-- Table for digamma function:
--
-- Uses equality ψ(n) = H_{n-1} - γ where
--   H_{n} = Σ 1/k, k = [1 .. n]     - harmonic number
--   γ     = 0.57721566490153286060  - Euler-Mascheroni number
digammaTestIntegers :: Double -> Bool
digammaTestIntegers eps
  = all (uncurry $ eq eps) $ take 3000 digammaInt
  where
    ok approx exact = approx
    -- Harmonic numbers starting from 0
    harmN = scanl (\a n -> a + 1/n) 0 [1::Rational .. ]
    gam   = 0.57721566490153286060
    -- Digamma values
    digammaInt = zipWith (\i h -> (digamma i, realToFrac h - gam)) [1..] harmN

----------------------------------------------------------------
-- Unit tests
----------------------------------------------------------------

-- Lookup table for fact factorial calculation. It has fixed size
-- which is bad but it's OK for this particular case
factorial_table :: V.Vector Integer
factorial_table = V.generate 2000 (\n -> product [1..fromIntegral n])

-- Exact implementation of factorial
factorial' :: Integer -> Integer
factorial' n = factorial_table ! fromIntegral n

-- Exact albeit slow implementation of choose
choose' :: Integer -> Integer -> Integer
choose' n k = factorial' n `div` (factorial' k * factorial' (n-k))

-- Truncate double to [0,1]
range01 :: Double -> Double
range01 = abs . (snd :: (Integer, Double) -> Double) . properFraction

readTable :: FilePath -> IO [[Double]]
readTable
  = fmap (fmap (fmap read . words) . lines)
  . readFile

forTable :: FilePath -> ([Double] -> Maybe String) -> IO ()
forTable path fun = do
  rows <- readTable path
  case mapMaybe fun rows of
    []   -> return ()
    errs -> assertFailure $ intercalate "---\n" errs

checkTabular :: Int -> String -> Double -> Double -> IO ()
checkTabular prec x exact val =
  case checkTabularPure prec x exact val of
    Nothing -> return ()
    Just s  -> assertFailure s

checkTabularPure :: Int -> String -> Double -> Double -> Maybe String
checkTabularPure prec x exact val
  | within prec exact val = Nothing
  | otherwise             = Just $ unlines
      [ " x         = " ++ x
      , " expected  = " ++ show exact
      , " got       = " ++ show val
      , " ulps diff = " ++ show (ulpDistance exact val)
      , " err.est.  = " ++ show prec
      ]

max_stars_repo_name: Shimuuar/math-functions

max_issues_repo_name: Shimuuar/math-functions

max_forks_repo_name: Shimuuar/math-functions


text:
module Hw4
    (
    Vec(..),
    range,
    absolute,
    rd,
    dft,
    fft) where

import Data.Complex
import GHC.Base

range :: Double -> Double -> Double -> [Double]
absolute :: Vec (Complex Double) -> Vec Double
rd :: Int -> Vec Double -> Vec Double
dft :: [Double] -> Vec (Complex Double)
fft :: [Double] -> Vec (Complex Double)

newtype Vec a = Vec {runVec :: [a]}

instance Show a => Show (Vec a) where
    show (Vec a) = "[" ++ show' a ++ "]"

show' :: Show a => [a] -> [Char]
show' [] = ""
show' [x] = show x
show' (x:xs) = show x ++ " " ++ show' xs

instance Functor Vec where
    fmap g x = Vec $ g <$> runVec x

instance Applicative Vec where
    pure x = Vec $ repeat x
    liftA2 f a b =  Vec $ zipWith f (runVec a) (runVec b)

instance Semigroup (Vec a) where
    (<>) (Vec a) (Vec b) = Vec $ a ++ b

instance Monoid (Vec a) where
    mempty = Vec []

instance Num a => Num (Vec a) where
    fromInteger x = Vec $ pure $ fromInteger x
    (+) = liftA2 (+)
    (-) = liftA2 (-)
    (*) = liftA2 (*)
    negate (Vec a) = Vec $ fmap negate a
    abs (Vec a) = Vec $ fmap abs a
    signum (Vec a) = Vec $ fmap signum a

instance Fractional a => Fractional (Vec a) where
    (/) = liftA2 (/)
    fromRational x = Vec $ pure $ fromRational x

instance Floating a => Floating (Vec a) where
    exp (Vec a) = Vec (fmap exp a)
    log (Vec a) = Vec (fmap log a)
    sin (Vec a) = Vec (fmap sin a)
    cos (Vec a) = Vec (fmap cos a)
    asin (Vec a) = Vec (fmap asin a)
    acos (Vec a) = Vec (fmap acos a)
    atan (Vec a) = Vec (fmap atan a)
    sinh (Vec a) = Vec (fmap sinh a)
    cosh (Vec a) = Vec (fmap cosh a)
    asinh (Vec a) = Vec (fmap asinh a)
    acosh (Vec a) = Vec (fmap acosh a)
    atanh (Vec a) = Vec (fmap atanh a)
    pi = Vec $ pure pi

range from to count = fmap (\x -> from + x * step) [0..count-1]
    where step = (to - from)/count

absolute = fmap (\x -> sqrt(realPart x * realPart x + imagPart x * imagPart x))

rd n = fmap (\x -> fromIntegral (round $ c * x) / c)
  where c = 10^n

dft x =
    let n = fromIntegral $ length x
        index = range 0 n n
        xn = x `zip` index

        f k = sum y
                where
                    y = map factor xn
                    factor (xi, j) = let y = 2 * pi * j * k / n
                                        in xi * cos y :+ (-xi * sin y)
    in
        Vec $ fmap f index

fft x
    | n <= 16 = dft x
    | otherwise =
        let (even, odd) = split x
            (e, o) = (fft even, fft odd)
            t = -2 * pi / fromIntegral n

            g (e', o', k) =
                        let
                            y = t * fromIntegral k
                            (fr, fi) = (cos y, sin y)
                            (pr, pi) = (fr * realPart o' - fi * imagPart o', fr * imagPart o' + fi * realPart o')
                        in
                            ((realPart e' + pr) :+ (imagPart e' + pi), (realPart e' - pr) :+ (imagPart e' - pi))

            (lower, upper) = unzip $ fmap g $ zip3 (runVec e) (runVec o) [0..]
        in
            Vec (lower ++ upper)

    where
            n = length x
            split [] = ([], [])
            split [a] = ([a], [])
            split (a:b:c) = (a:x, b:y)
                where (x,y) = split c

max_stars_repo_name: NikolaPeevski/Haskell-stuff

max_issues_repo_name: NikolaPeevski/Haskell-stuff

max_forks_repo_name: NikolaPeevski/Haskell-stuff


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE BangPatterns #-}
module AI.Singularity.Training where
import GHC.TypeLits
import qualified Numeric.LinearAlgebra.Data as LD
import qualified Numeric.LinearAlgebra.Devel as LDev
import qualified Numeric.LinearAlgebra as LA
import qualified Control.Monad.State as S
import Control.Lens
import Data.Proxy
import Control.Monad.ST
import System.Random
import Control.Monad
import Control.Arrow
import qualified Data.Vector.Storable as VS
import Data.Conduit
import System.IO (hPrint, stderr)
import qualified Data.Conduit.Binary as CB
import qualified Data.Conduit.List   as CL
import Data.Conduit.Internal (zipSources)
import Numeric.AD
import Control.Monad.Trans.Class (lift)

import AI.Singularity.Utils
import AI.Singularity.Data.Matrix
import AI.Singularity.Data.Vector
import AI.Singularity.Data.Network

data TrainConfig = TrainConfig { _learningRate :: !Float,
                                 _lrUpdate     :: Float -> Float,
                                 _batchsize    :: !Int,
                                 _momentumRate :: !Float }
makeLenses ''TrainConfig

newtype Loss = Loss { getLoss :: forall a. Floating a => a -> a -> a}

lossGrad :: Loss -> (forall a. Floating a => a -> a -> a)
lossGrad (Loss f) = f'
  where fList :: forall a. Floating a => [a] -> a
        fList [y,y'] = f y y'
        f' y y'      = head . tail . grad fList $ [y,y']

data TrainState i o = TrainState { _config  :: !TrainConfig,
                                   _network :: !(TrainedNetwork i o),
                                   _iter    :: !Int,
                                   _loss    :: Loss,
                                   _epochs  :: !Int,
                                   _evalSet :: [(i,o)]}

makeLenses ''TrainState

addGrad :: (KnownNat m, KnownNat n) => TrainConfig -> InTrainMat m n -> Identity (InTrainMat m n)
addGrad tr mp = Identity $ GeneralInTrainingMat newWeights newGradient hGr
  where momentumRatio = view momentumRate tr
        lr            = view learningRate tr
        newWeights    = checkNan $ view trMatr mp + (mmap (*lr) $ view trGrad mp) / mmap sqrt (hGr + 0.00000001)
        newGradient   = checkNan . over matr (LD.cmap (*momentumRatio)) $ view trGrad mp
        hGr           = 0.9 * view hGrad mp + 0.1 * (view trGrad mp ^ 2)

applyGrad :: TrainConfig -> TrainedNetwork inp out -> TrainedNetwork inp out
applyGrad tr = runIdentity . changeMatType (addGrad tr)

calcGrads :: forall n m. (KnownNat n, KnownNat m) => Vector n -> Vector m -> TrainedNetwork (Vector n) (Vector m) -> (Matrix (n+1) m, Vector n)
calcGrads !inp !outGrad (FFLayer !mpair !f) = (gradMat, inpGrad)
  where
    m       = view trMatr mpair
    gradMat = dy `outer` consV 1 inp
    dy      = outGrad * diffM
    diffM   = cmap (diff f) . (#>) m . consV 1 $ inp
    inpGrad = tailV $ transposeM m #> dy

trainWithErr :: forall a b. TrainConfig -> a -> b -> TrainedNetwork a b -> (TrainedNetwork a b, a)
trainWithErr tc !inp !outGrad n@(FFLayer !mpair !f) = (FFLayer  m' f, inpGrad)
  where lr = view learningRate tc
        m' = over trGrad (subtract gradMat) mpair
        (gradMat,inpGrad) = calcGrads inp outGrad n
trainWithErr tc !inp !outGrad (FFSeq !n1 !n2) = (FFSeq n1' n2', inpGrad)
  where (n1', inpGrad)  = trainWithErr tc inp inpGrad' n1
        inp'            = conductSignalInTrain inp n1
        (n2', inpGrad') = trainWithErr tc inp' outGrad n2
trainWithErr tc !inp !outGrad (SplitNet n) = (SplitNet *** splitV) unSpl
  where unSpl = trainWithErr tc (uncurry appendV inp) outGrad n
trainWithErr tc !inp !outGrad (FromBinFunc f) = (FromBinFunc f, grads)
  where grads = (\(x1,x2) -> ((!!0) &&& (!!1)). map (*outGrad) . fgrad $ [x1,x2]) inp
        fwrap :: forall a. Fractional a => [a] -> a
        fwrap [x1,x2] = f x1 x2
        fgrad = grad fwrap
trainWithErr tc !inp !outGrad (Split n1 n2) = (Split n1' n2', inGr)
  where (n1', inGr1) = trainWithErr tc inp (fst outGrad) n1
        (n2', inGr2) = trainWithErr tc inp (snd outGrad) n2
        inGr         = inGr1 + inGr2
trainWithErr tc !inp !outGrad (FromFunc f) = (FromFunc f, inGr)
  where inGr = diff f inp * outGrad
trainWithErr tc !inp !outGrad (FFDiv !n1 !n2) = (FFDiv n1' n2', (inGr1, inGr2))
  where (n1', inGr1) = trainWithErr tc (fst inp) (fst outGrad) n1
        (n2', inGr2) = trainWithErr tc (snd inp) (snd outGrad) n2
trainWithErr tc !inp !outGrad (Recurse n)   = if null inp then (Recurse n, []) else (Recurse n', [fst inGrad])
   where inputs      = scanl (\acc sym -> (sym, conductSignalInTrain acc n)) (head inp,0) (tail inp)
         (n',inGrad) = foldr (\input (net, outG) -> trainWithErr tc input (snd outG) net) (n,(undefined, outGrad)) inputs
trainWithErr tc !inp !outGrad (First n)     = (First n', inGrad)
  where (n', in1Grad) = trainWithErr tc (fst inp) outGrad n
        inGrad        = (in1Grad, 0)
trainWithErr tc !inp !outGrad (Second n)     = (Second n', inGrad)
  where (n', in1Grad) = trainWithErr tc (snd inp) outGrad n
        inGrad        = (0, in1Grad)

trainSample :: forall a b. TrainConfig -> [(a,b)] -> (b -> b -> b) -> Int -> Network a b -> Network a b
trainSample tc !corpus gLoss !epochs = fromTrain . (trainHelp tc 1 . concat $ replicate epochs corpus) . toTrain
  where
    trainHelp :: TrainConfig -> Int -> [(a, b)] -> TrainedNetwork a b -> TrainedNetwork a b
    trainHelp tc i ((!x,!y):corp) !net = next . fst $ trainWithErr tc x l net
      where l    = gLoss y (conductSignalInTrain x net)
            next = if i < b
                   then trainHelp tc (i+1) corp
                   else trainHelp nTc 1 corp . applyGrad tc
            nTc  = tc -- over learningRate (view lrUpdate tc) tc
            b    = view batchsize tc
    trainHelp _ _ [] !net = net

eps :: forall n. Floating n => n
eps = 0.00000001

defaultConfig = TrainConfig 0.05 id 1 0

endOfBatch :: TrainState i o -> Bool
endOfBatch = view iter &&& view (config.batchsize) >>> uncurry (>=)

startState conf net = defaultState conf net logloss
  where logloss y y' = negate $ y * log (y' + eps) + (1 - y) * log (1 - y' + eps)

defaultState :: TrainConfig -> Network i o -> (forall a. Floating a => a -> a -> a) -> Int -> [(i,o)] -> TrainState i o
defaultState conf net f ep = TrainState conf (toTrain net) 1 (Loss f) ep

-- applyTrainingExample :: (Floating o, Show o) => TrainState i o -> (i,o) -> IO (TrainState i o)
applyTrainingExample tr (i,o) = endBatch . over iter (+1) . over network (fst . trainWithErr trConf i outLoss) $ tr
  where
    outLoss  = lossGrad (view loss tr) o (conductSignalInTrain i (view network tr))
    trConf   = view config tr
    endBatch = if endOfBatch tr
               then \x -> do
                           let up = view (config.lrUpdate) x
                               nX = over (config.learningRate) up . set iter 1 . over network (applyGrad trConf) $ x
                               s  = view evalSet tr
                               ac = sum . map (uncurry. flip . getLoss . view loss $ tr) . map (first $ flip conductSignalInTrain (view network nX)) $ s
                           unless (null s) . lift $ do
                             hPrint stderr "LOSS:"
                             hPrint stderr (ac/fromIntegral (length s))
                           return nX
               else return

trainC cond conf net ep = trainStateC cond . startState conf net ep

trainStateC cond trState =
  if view epochs trState >= 1
  then
    runConduitRes (cond .| CL.foldM applyTrainingExample trState) >>= (over epochs (subtract 1) >>> trainStateC cond)
  else
    fromTrain . applyGrad (view config trState) . view network <$> runConduitRes (cond .| CL.foldM applyTrainingExample trState)

trainExample :: forall a b. Floating b => TrainConfig -> [(a, b)] -> Int -> Network a b -> Network a b
trainExample conf x = trainSample conf x logloss'
  where
    logloss [y, y'] = y * log (y' + eps) + (1 - y) * log (1 - y' + eps)
    logloss' y y'   = head . tail . grad logloss $ [y, y']

max_stars_repo_name: Antystenes/Memetic-Predictor

max_issues_repo_name: Antystenes/Memetic-Predictor

max_forks_repo_name: Antystenes/Memetic-Predictor


text:
-- Mandelbrot.hs
-- Caluclate Mandelbrot set in Haskell
--
-- vim: ft=haskell sw=2 ts=2 et
--

module Mandelbrot
  ( Pixel
  , mandelbrot ) where

import Data.Complex

-- Basic type defs
type Pixel = (Int, Int, Int)

-- Calculate pixel value
pixelValue :: (Integral a) => Double -> Double -> a -> a
pixelValue x y nMax =
  let z0 = x :+ y
      z = 0.0 :+ 0.0 in
    loop z z0 nMax
  where
    loop :: (Integral a) => Complex Double -> Complex Double -> a -> a
    loop z z0 nMax =
      nMax

valueToColor :: (Integral a) => a -> Pixel
valueToColor val =
  (0, 0, 0)

mandelbrot :: (Integral a, Floating b) => a -> a -> b -> b -> b -> [Pixel]
mandelbrot xMax yMax xCenter yCenter pixelSize =
  [(0, 0, 0)]

max_stars_repo_name: kkirstein/proglang-playground

max_issues_repo_name: kkirstein/proglang-playground

max_forks_repo_name: kkirstein/proglang-playground


text:
module DistanceMatrixRank where

import Data.List as L
import Data.Vector as V
import Data.Matrix as M

import Statistics.Sample

data DistanceMatrix = DistanceMatrix
  {
    distanceMatrix :: Matrix Double,
    origins :: [String],
    destinations :: [String]
  }
  deriving (Show)

rankDestinations :: DistanceMatrix -> [(String, Double)]
rankDestinations m = L.zip (destinations m) (L.map stdDev $ columns (distanceMatrix m))

columns :: Matrix a -> [Vector a]
columns m = L.map (flip M.getCol m) [1..cols]
  where
    cols = M.ncols m

max_stars_repo_name: WilliamDo/optimal-prime

max_issues_repo_name: WilliamDo/optimal-prime

max_forks_repo_name: WilliamDo/optimal-prime


text:
-- | Functionality for graphing 2-dimensional matrices.
module System.Console.Ansigraph.Internal.Matrix (
    matShow
  , displayMat
  , displayCMat
) where

import System.Console.Ansigraph.Internal.Core

import Data.Complex
import Data.List (intersperse)
import Control.Monad.IO.Class (MonadIO, liftIO)
-- for GHC <= 7.8
import Control.Applicative

---- Matrices ----

mmap :: (a -> b) -> [[a]] -> [[b]]
mmap = map . map

mmax :: (Num a, Ord a) => [[a]] -> a
mmax = maximum . map maximum . mmap abs

densityChars = "█▓▒░"

densityVals :: [Double]
densityVals = (+ 0.125) . (/4) <$> [3,2,1,0]
         -- = [7/8, 5/8, 3/8, 1/8]

blocks :: [(Double,Char)]
blocks = zip densityVals densityChars

data MatElement = MatElement !Bool {-# UNPACK #-} !Char

elemChar :: MatElement -> Char
elemChar (MatElement _ c) = c

putRealElement :: MonadIO m => GraphSettings -> MatElement -> m ()
putRealElement s (MatElement b c) = colorStr clring (c : " ")
  where clr    = if b then realNegColor s else realColor s
        clring = mkColoring clr (realBG s)

putImagElement :: MonadIO m => GraphSettings -> MatElement -> m ()
putImagElement s (MatElement b c) = colorStr clring $ c : " "
  where clr    = if b then imagNegColor s else imagColor s
        clring = mkColoring clr (imagBG s)

selectMatElement :: Double -> MatElement
selectMatElement x = let l = filter (\p -> fst p < abs x) blocks in case l of
  []    -> MatElement False   ' '
  (p:_) -> MatElement (x < 0) (snd p)

matElements :: [[Double]] -> [[MatElement]]
matElements m = let mx = mmax m
                in  mmap (selectMatElement . (/ mx)) m

-- | Given a matrix of Doubles, return the list of strings illustrating the absolute value
--   of each entry relative to the largest, via unicode chars that denote a particular density.
--   Used for testing purposes.
matShow :: [[Double]] -> [String]
matShow = mmap elemChar . matElements

newline = putStrLn' ""

intersperse' x l = intersperse x l ++ [x]

-- | Use ANSI coloring (specified by an 'GraphSettings') to visually display a Real matrix.
displayMat :: MonadIO m => GraphSettings -> [[Double]] -> m ()
displayMat s = liftIO . sequence_ . concat . intersperse' [newline] . displayRealMat s

-- | Use ANSI coloring (specified by a 'GraphSettings') to visually display a Real matrix.
displayRealMat :: MonadIO m => GraphSettings -> [[Double]] -> [[m ()]]
displayRealMat s = mmap (putRealElement s) . matElements

-- | Use ANSI coloring (specified by a 'GraphSettings') to visually display a Real matrix.
displayImagMat :: MonadIO m => GraphSettings -> [[Double]] -> [[m ()]]
displayImagMat s = mmap (putImagElement s) . matElements

-- | Use ANSI coloring (specified by an 'GraphSettings') to visually display a Complex matrix.
displayCMat :: MonadIO m => GraphSettings -> [[Complex Double]] -> m ()
displayCMat s m = liftIO . sequence_ . concat . intersperse' [newline] $
  zipWith (\xs ys -> xs ++ putStr " " : ys)
          (displayRealMat s $ mmap realPart m)
          (displayImagMat s $ mmap imagPart m)

max_stars_repo_name: fieldstrength/ansigraph

max_issues_repo_name: fieldstrength/ansigraph

max_forks_repo_name: fieldstrength/ansigraph


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE Arrows, FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE DeriveGeneric #-}

module Examples.MultiCamera where

import Control.Arrow hiding ((|||))
import Control.Monad (replicateM, forM)
import Control.Monad.Trans (MonadIO, liftIO, lift)
import Control.Monad.Bayes.Class (MonadSample, MonadCond, logCategorical)
import Control.Monad.Bayes.Sampler (SamplerIO, sampleIO)

import Data.Aeson hiding (Result)
import Data.Maybe (catMaybes)

import GHC.Generics (Generic)

import qualified Data.ByteString.Lazy.Char8 as BS (putStrLn)

import Numeric.LinearAlgebra.Static

import Inference (zdsparticles, zunheap)
import DelayedSampling (MDistr (..), children, state, State, MarginalT (..), SMarginalT(..), distr, typeOfDSDistr, DelayedInfer)
import qualified SymbolicDistr as DS
import DSProg (DeepForce (..), Result (..), Expr' (..), Expr, marginal, zdeepForce)
import Distributions (normal, poisson, sample, observe, bernoulli, replicateIID, passert, factor, uniformIntRange)
import MVDistributions (shuffleList, mvNormal, uniformNSphere)
import Util.ZStream (ZStream)
import qualified Util.ZStream as ZS
import Util.Ref (MonadState, Heap, emptyHeap, readRef)

import Numeric.SpecFunctions (logGamma)
import Numeric.Log (Log (Exp), ln)

import qualified Data.Vector.Storable as V (fromList)

import qualified Metaprob as MP

numCameras :: Int
numCameras = 2

numPoses :: Int
numPoses = 10

type CameraID = Int
type Pose = Int

data TrackG pv = Track
  { posWH :: pv
  , startTime :: Double
  , trackID :: Int
  , identity :: Int
  , camera :: CameraID
  }
  deriving (Generic, Show)

instance ToJSON pv => ToJSON (TrackG pv)

instance DeepForce pv => DeepForce (TrackG pv) where
  type Forced (TrackG pv) = TrackG (Forced pv)
  deepForce (Track pv h t id c) = Track <$> deepForce pv <*> deepForce h <*> pure t <*> pure id <*> pure c
  deepConst (Track pv h t id c) = Track (deepConst pv) (deepConst h) t id c

type STrack = TrackG (Expr (R 4))
type Track = TrackG (R 4)
type MarginalTrack = TrackG (Result (R 4))

type Appearance = Expr (R 10)

tdiff :: Double
tdiff = 1

posWHCovBlocks :: Sym 2 -> Sym 2 -> Sym 4
posWHCovBlocks pcov whcov = sym $
  ((unSym pcov ||| (konst 0 :: Sq 2))
        ===
  ((konst 0 :: Sq 2) ||| unSym whcov))

newTrack'' :: MonadState Heap m => MonadSample m => Double -> Int -> Int -> Int -> m STrack
newTrack'' t tid id c = do
  pv <- DS.sample (DS.mvNormal (Const mu) cov)
  pure (Track pv t tid id c)
  where
  mu = (konst 0 :: R 2) # (konst 1 :: R 2)
  cov = posWHCovBlocks (sym eye :: Sym 2) (0.2 * sym eye :: Sym 2)

newTrack :: MonadState Heap m => MonadSample m => [Appearance] -> Double -> Int -> m (STrack, [Appearance])
newTrack appearances t tid = do
  camera <- sample (uniformIntRange numCameras)
  newAppearance <- sample (bernoulli (1 / (fromIntegral (length appearances) + 1)))
  (appearanceID, appearances') <- if newAppearance
    then do
         app <- DS.sample (DS.mvNormal (Const 0) (sym eye))
         pure (length appearances, appearances ++ [app])
    else do
         appID <- sample (uniformIntRange (length appearances))
         pure (appID, appearances)
  newTrack <- newTrack'' t tid appearanceID camera
  pure (newTrack, appearances')

trackMotion :: MonadState Heap m => MonadSample m => Double -> STrack -> m STrack
trackMotion tdiff track = do
  pv' <- DS.sample (DS.mvNormal (MVMul (Const motionMatrix) (posWH track)) motionCov)
  pure $ track { posWH = pv' }
  where
  motionMatrix = eye :: Sq 4
  motionCov :: Sym 4
  motionCov = sym (konst tdiff) * posWHCovBlocks (sym eye) (sym (konst 0.1 * eye))

trackSurvivalMotion :: MonadState Heap m => MonadSample m => Double -> STrack -> m (Maybe STrack)
trackSurvivalMotion tdiff track = do
  survived <- sample (bernoulli (exp (- tdiff * deathRate)))
  if survived
    then Just <$> trackMotion tdiff track
    else pure Nothing
  where
  deathRate = 0.02

trackMeasurement :: [Appearance] -> STrack -> DS.Distr (Expr (R 4), (R 10, Expr (R 1)))
trackMeasurement appearances track =
  DS.indep (DS.mvNormal (posWH track) (sym eye)) $
  DS.bind uniformNSphere $ \pose ->
    DS.mvNormal (MVMul (Const (row pose)) (appearances !! identity track)) (sym eye)

-- Naive association
associationWithClutter :: forall a b. DeepForce b => DS.Distr [b] -> (a -> DS.Distr b) -> [a] -> DS.Distr [b]
associationWithClutter cluttersD obsD tracks = DS.Distr is iobs undefined where
  is :: MonadState Heap m => MonadSample m => m [b]
  is = do
    clutter <- DS.sample cluttersD
    observations <- mapM (DS.sample . obsD) tracks
    shuffleList (observations ++ clutter)
  defaultObserve :: DelayedInfer m => [Forced b] -> m ()
  defaultObserve observations = do
    passert (length tracks <= length observations)
    observations' <- shuffleList observations
    let (observations'', clutter) = splitAt (length tracks) observations'
    mapM_ (\(t, o) -> DS.observe (obsD t) o) (zip tracks observations'')
    DS.observe cluttersD clutter
  iobs :: DelayedInfer m => [Forced b] -> m ()
  iobs = if True
    then assocWithClutterCustomProposal cluttersD obsD tracks
    else defaultObserve

assocWithClutterCustomProposal :: DelayedInfer m => DeepForce b => DS.Distr [b] -> (a -> DS.Distr b) -> [a] -> [Forced b] -> m ()
assocWithClutterCustomProposal cluttersD obsD allTracks allObservations = go allTracks allObservations where
  numTracks = fromIntegral (length (allTracks))
  numObservations = fromIntegral (length (allObservations))
  go [] clutter = do
    let numClutter = fromIntegral (length clutter)
    factor (- (logGamma (numTracks + numClutter + 1) - logGamma (numClutter + 1)))
    DS.observe cluttersD clutter
  go (track : tracks) [] =
    passert False
  go (track : tracks) observations = do
    let obsDistr = obsD track
    lls <- map Exp <$> mapM (DS.score obsDistr) observations
    let totalLL = sum lls
    let adjLLs = map (/ totalLL) lls
    i <- logCategorical (V.fromList adjLLs)
    let (os1, obs : os2) = splitAt i observations
    DS.observe obsDistr obs
    factor (ln (recip (adjLLs !! i))) -- proposal correction
    go tracks (os1 ++ os2)

tracksMeasurement :: [Appearance] -> [STrack] -> DS.Distr [(Expr (R 4), (R 10, Expr (R 1)))]
tracksMeasurement appearances = associationWithClutter cluttersDistr (trackMeasurement appearances)
  where
  cluttersDistr = DS.replicateIID (poisson amtClutter) clutterDistr
  amtClutter :: Double
  amtClutter = 1
  clutterDistr :: DS.Distr (Expr (R 4), (R 10, Expr (R 1)))
  clutterDistr = DS.indep (DS.mvNormal (Const ((konst 0 :: R 2) # (konst 1 :: R 2))) (sym eye))
    $ DS.bind (mvNormal 0 (sym eye)) (\_ -> DS.mvNormal (Const 0) (sym eye))

-- NOTE: Currently, it is possible that we will generate a new track
-- with the same identity and camera as an existing track.
-- Should consider how to change the model so that this doesn't happen.
tracksMotion :: MonadState Heap m => MonadSample m => Double -> Double -> [STrack] -> Int -> [Appearance] -> m ([STrack], Int, [Appearance])
tracksMotion t tdiff tracks numTracks appearances = do
  liveTracks <- catMaybes <$> mapM (trackSurvivalMotion tdiff) tracks
  numNewTracks <- sample (poisson (birthRate * tdiff))
  (newTracks, appearances') <- generateNewTracks appearances (numTracks + 1) numNewTracks
  pure (liveTracks ++ newTracks, numTracks + numNewTracks, appearances')
  where
  generateNewTracks apps nt 0 = pure ([], apps)
  generateNewTracks apps nt n = do
    (newT, apps') <- newTrack apps t nt
    (\(x, y) -> (newT : x, y)) <$> generateNewTracks apps' (nt + 1) (n - 1)
  birthRate = 0.1

-- Note how extraneous state about the number of tracks and the time is kept hidden
zstepGen :: MonadState Heap m => MonadSample m => ZStream (MP.Gen m) () ([STrack], [(CameraID, R 10, Expr (R 4), Expr (R 1))])
zstepGen = ZS.fromStep stepf initState
  where
  initState :: (Double, [TrackG pv], Int, [Appearance])
  initState = (0, [], 0, [])
  stepf (t, tracks, numTracks, appearances) () = do
    (tracks', newNumTracks, appearances') <- lift (tracksMotion t tdiff tracks numTracks appearances)
    observations <- forM [0 .. numCameras - 1] $ \camID -> do
      obs <- ("obs." ++ show camID) MP.~~ MP.dsPrim (tracksMeasurement appearances' (filter (\t -> camera t == camID) tracks'))
      return [ (camID, pose, pvwh, obsAppear) | (pvwh, (pose, obsAppear)) <- obs ]
    return ((t + tdiff, tracks', newNumTracks, appearances'), (tracks', concat observations))

processObservationsStream :: DelayedInfer m => ZStream m [(CameraID, R 10, R 4, R 1)] [MarginalTrack]
processObservationsStream = proc observations -> do
  (tracks', _) <- MP.zobserving zstepGen -< ((), "obs" MP.|-> MP.obs observations)
  ZS.run -< mapM trackf tracks'
  where
  trackf :: DelayedInfer m => STrack -> m MarginalTrack
  trackf track = do
    Just pv <- marginal (posWH track)
    pure $ track { posWH = pv }

generateGroundTruth :: MonadState Heap m => MonadSample m => ZStream m () ([STrack], [(CameraID, R 10, Expr (R 4), Expr (R 1))])
generateGroundTruth = ZS.liftM MP.sim zstepGen

runMTTPF :: Int -> ZStream SamplerIO () ([Track], [[MarginalTrack]], [(CameraID, R 10, R 4, R 1)])
runMTTPF numParticles = proc () -> do
  (groundTruth, obs) <- zdeepForce generateGroundTruth -< ()
  particles <- zdsparticles numParticles processObservationsStream -< obs
  returnA -< (groundTruth, particles, obs)

runExample :: IO ()
runExample = sampleIO $ ZS.runStream (liftIO . BS.putStrLn . encode) (runMTTPF 100)

max_stars_repo_name: psg-mit/probzelus-haskell

max_issues_repo_name: psg-mit/probzelus-haskell

max_forks_repo_name: psg-mit/probzelus-haskell