text:
module MachineLearning.NeuralNetworkTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus

import MachineLearning.DataSets (dataset2)

import qualified Control.Monad.Random as RndM
import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning as ML
import qualified MachineLearning.Optimization as Opt
import MachineLearning.Model
import MachineLearning.NeuralNetwork
import qualified MachineLearning.NeuralNetwork.TopologyMaker as TM

(x, y) = ML.splitToXY dataset2

gradientCheckingEps = 0.1

checkGradientTest eps activation loss lambda = do
  let nnt = TM.makeTopology activation loss (LA.cols x) 2 [10]
      model = NeuralNetwork nnt
      thetas = initializeTheta 1511197 nnt
      diffs = take 5 $ map (\e -> Opt.checkGradient model lambda x y thetas e) [0.005, 0.0051 ..]
      diff = minimum $ filter (not . isNaN) diffs
  assertApproxEqual (show thetas) eps 0 diff

xPredict = LA.matrix 2 [ -0.5, 0.5
                       , 0.2, -0.2
                       , 1, 1
                       , 1, 0
                       , 0, 0]
yExpected = LA.vector [1, 1, 0, 0, 1]

learnTest activation loss minMethod nIters =
  let lambda = L2 $ 0.5 / (fromIntegral $ LA.rows x)
      x1 = ML.mapFeatures 2 x
      nnt = TM.makeTopology activation loss (LA.cols x1) 2 [10]
      model = NeuralNetwork nnt
      xPredict1 = ML.mapFeatures 2 xPredict
      initTheta = initializeTheta 5191711 nnt
      (theta, optPath) = Opt.minimize minMethod model 1e-7 nIters lambda x1 y initTheta
      yPredicted = hypothesis model xPredict1 theta
      js = (LA.toColumns optPath) !! 1
  in do
    assertVector (show js) 0.01 yExpected yPredicted

tests = [ testGroup "gradient checking" [
            testCase "Sigmoid - Logistic: non-zero lambda" $ checkGradientTest 0.1 TM.ASigmoid TM.LLogistic (L2 0.01)
            , testCase "Sigmoid - Logistic: zero lambda" $ checkGradientTest 0.1 TM.ASigmoid TM.LLogistic (L2 0)
            , testCase "ReLU - Softmax: non-zero lambda" $ checkGradientTest 0.1 TM.ARelu TM.LSoftmax (L2 0.01)
            , testCase "ReLU - Softmax: zero lambda" $ checkGradientTest 0.1 TM.ARelu TM.LSoftmax (L2 0)
            , testCase "Tanh - MultiSvm: non-zero lambda" $ checkGradientTest 0.1 TM.ATanh TM.LMultiSvm (L2 0.01)
            , testCase "Tanh - MultiSvm: zero lambda" $ checkGradientTest 0.1 TM.ATanh TM.LMultiSvm (L2 0)
            , testCase "Tanh - MultiSvm: no reg" $ checkGradientTest 0.1 TM.ATanh TM.LMultiSvm RegNone
            ]
        , testGroup "learn" [
            testCase "Sigmoid - Logistic: BFGS" $ learnTest TM.ASigmoid TM.LLogistic (Opt.BFGS2 0.01 0.7) 50
            , testCase "ReLU - Softmax: BFGS" $ learnTest TM.ARelu TM.LSoftmax (Opt.BFGS2 0.1 0.1) 50
            , testCase "Tanh - MultiSvm: BFGS" $ learnTest TM.ATanh TM.LMultiSvm (Opt.BFGS2 0.1 0.1) 50
            ]
        ]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Isumi.Math.FFT.Image
  ( fftWithImage
  ) where

import           Codec.Picture         (Image, Pixel8)
import           Data.Complex
import           Isumi.Math.FFT.Bruun  (fft2dBruun)
import           Isumi.Math.MatrixUtil (imageToMatrix)
import           Numeric.LinearAlgebra (Matrix)

fftWithImage :: Image Pixel8 -> Maybe (Matrix (Complex Double))
fftWithImage = fft2dBruun . imageToMatrix

ext: hs

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun


text:
module Newton.Render where

-- this module comes from JuicyPixels package
import qualified Codec.Picture            as J
import           Codec.Picture.ColorQuant (defaultPaletteOptions)
import           Data.Complex
import           Newton.Utils             (root)

-- | Add brightness to an RGB8 Image
brightnessRGB8 :: Int -> J.Palette -> J.Palette
brightnessRGB8 add = J.pixelMap brightFunction
  where
    up v = fromIntegral (fromIntegral v + add)
    brightFunction (J.PixelRGB8 r g b) = J.PixelRGB8 (up r) (up g) (up b)

-- | Save an PNG image as a file
saveBasinsPic ::
     FilePath -- ^ Where to write the image file.
  -> J.Palette
  -> IO ()
saveBasinsPic = J.writePng

saveBasinsAnim ::
     FilePath -- ^ Filename of an image
  -> J.GifDelay -- ^ Delay between frames
  -> [J.Palette] -- ^ List of pictures to store within the GIF
  -> IO ()
saveBasinsAnim path delay images =
  case J.writeGifImages path J.LoopingForever finalImgs of
    Left msg  -> error msg
    Right ans -> ans
  where
    palettized = map (J.palettize defaultPaletteOptions) images
    finalImgs = map (\(img, plt) -> (plt, delay, img)) palettized

-- | Render attraction basins for zeros (roots) of a complex function
-- as a PNG image.
renderBasins ::
     (Fractional a, RealFloat a)
  => Int -- ^ Size of the square image in pixels
  -> (Complex a -> Complex a) -- ^ A function f to find zero for.
  -> (Complex a -> Complex a) -- ^ A derivative of function f.
  -> (Complex a -> J.PixelRGB8) -- ^ How to color roots.
  -> J.Palette
renderBasins commonSize f f' color = image
  where
    image = J.generateImage render width height
    width = commonSize -- width in pixels
    height = commonSize -- height in pixels
    (xfrom, xto) = (-2, 2) -- range for the real axis
    (yfrom, yto) = (-2, 2) -- range for the imaginary axis
    -- | Convert pixel coordinates to a complex number.
    pixelToComplex (i, j) = re :+ im
      where
        re = xfrom + (xto - xfrom) * fromIntegral i / fromIntegral width
        im = yfrom + (yto - yfrom) * fromIntegral j / fromIntegral height
    n = 100 -- maximum number of iterations
    p e = magnitude e < 1e-7 -- precision tolerance
    -- | Convert actual number of interations into
    -- into a grey pixel using logarithmic scale
    -- (this helps bring out more details in the image).
    logscale :: Integral a => Int -> a -> a
    logscale k c = floor (fromIntegral c * (1 - logBase n' k'))
      where
        n' = fromIntegral n + 1
        k' = fromIntegral k + 1
    -- | Render one pixel.
    render :: Int -> Int -> J.PixelRGB8
    render i j =
      case root n p f f' (pixelToComplex (i, j)) of
        Nothing -> J.PixelRGB8 0 0 0 -- black
        Just (k, xr) ->
          let J.PixelRGB8 r g b = color xr
           in J.PixelRGB8 (logscale k r) (logscale k g) (logscale k b)

-- | Simple (unstructured) complex number coloring.
simpleColoring :: Complex Double -> J.PixelRGB8
simpleColoring z = J.PixelRGB8 r g b
  where
    r = floor ((2 ^ 8 - 1) * r')
    g = floor ((2 ^ 8 - 1) * g')
    b = floor ((2 ^ 8 - 1) * b')
    (r', g', b') = hsl (phase z) 1 (1 - 0.5 ** magnitude z)

-- | Convert HSL color to RGB.
-- This function is copied with slight changes
-- from http://hackage.haskell.org/package/colour-2.3.4
hsl :: (RealFrac a, Floating a, Ord a) => a -> a -> a -> (a, a, a)
hsl h s l = (component tr, component tg, component tb)
  where
    hk = h / (2 * pi)
    (tr, tg, tb) = (mod1 (hk + 1 / 3), mod1 hk, mod1 (hk - 1 / 3))
    q
      | l < 0.5 = l * (1 + s)
      | otherwise = l + s - l * s
    p = 2 * l - q
    component t
      | t < 1 / 6 = p + ((q - p) * 6 * t)
      | t < 1 / 2 = q
      | t < 2 / 3 = p + ((q - p) * 6 * (2 / 3 - t))
      | otherwise = p
    mod1 x =
      if pf < 0
        then pf + 1
        else pf
      where
        (_, pf) = properFraction x

max_stars_repo_name: unaimillan/NewtonFractals

max_issues_repo_name: unaimillan/NewtonFractals

max_forks_repo_name: unaimillan/NewtonFractals


text:
{-# language MagicHash #-}

module Prelude
  ( module P
  ) where

import Data.Either as P (Either(..))
import Streaming.Prelude as P (Stream, Of)
import Numeric as P (pi)
import Control.Monad as P (Monad(..))
import Data.Primitive.Instances ()
import Data.Complex as P (Complex(..))
import Control.Monad.ST as P (ST, runST)
import GHC.Err as P (error)
import GHC.Exts as P (Double(..),Int(..),Int#)
import Data.Function as P (($), id)
import Data.Bool as P (otherwise, Bool(..))
import Data.Bits as P (Bits(..))
import Data.Semigroup as P (Semigroup(..))
import Data.Monoid as P (Monoid(..))
import Control.Applicative as P (Applicative(..))
import Data.Semiring as P (Semiring(..), Ring(..), (+),(*),(-))
import Data.Int as P (Int)
import Data.Word as P (Word)
import GHC.Real as P (fromIntegral, (/), floor)
import Data.Eq as P (Eq(..))
import Data.Ord as P (Ord(..))
import Control.Monad.Primitive as P (PrimMonad(..))
import Data.Primitive.Types as P (Prim(..))
import Data.Primitive.PrimArray as P (PrimArray,MutablePrimArray)
import Data.Primitive.Contiguous as P (Contiguous,Element,Mutable)

max_stars_repo_name: chessai/streaming-fft

max_issues_repo_name: chessai/streaming-fft

max_forks_repo_name: chessai/streaming-fft


text:
{-# LANGUAGE
  DataKinds,
  DeriveDataTypeable,
  FlexibleInstances,
  GADTs,
  MultiParamTypeClasses,
  TemplateHaskell,
  TypeFamilies #-}

{- |
Module      :  Data.Real.Constructible
Description :  Constructible real numbers
Copyright   :  © Anders Kaseorg, 2013
License     :  BSD-style

Maintainer  :  Anders Kaseorg <andersk@mit.edu>
Stability   :  experimental
Portability :  Non-portable (GHC extensions)

The constructible reals, 'Construct', are the subset of the real
numbers that can be represented exactly using field operations
(addition, subtraction, multiplication, division) and positive square
roots.  They support exact computations, equality comparisons, and
ordering.

>>> [((1 + sqrt 5)/2)^n - ((1 - sqrt 5)/2)^n :: Construct | n <- [1..10]]
[sqrt 5,sqrt 5,2*sqrt 5,3*sqrt 5,5*sqrt 5,8*sqrt 5,13*sqrt 5,21*sqrt 5,34*sqrt 5,55*sqrt 5]

>>> let f (a, b, t, p) = ((a + b)/2, sqrt (a*b), t - p*((a - b)/2)^2, 2*p)
>>> let (a, b, t, p) = f . f . f . f $ (1, 1/sqrt 2, 1/4, 1 :: Construct)
>>> floor $ ((a + b)^2/(4*t))*10**40
31415926535897932384626433832795028841971

>>> let qf (p, q) = ((p + sqrt (p^2 - 4*q))/2, (p - sqrt (p^2 - 4*q))/2 :: Construct)
>>> let [(v, w), (x, _), (y, _), (z, _)] = map qf [(-1, -4), (v, -1), (w, -1), (x, y)]
>>> z/2
-1/16 + 1/16*sqrt 17 + 1/8*sqrt (17/2 - 1/2*sqrt 17) + 1/4*sqrt (17/4 + 3/4*sqrt 17 - (3/4 + 1/4*sqrt 17)*sqrt (17/2 - 1/2*sqrt 17))

Constructible complex numbers may be built from constructible reals
using 'Complex' from the complex-generic library.

>>> (z/2 :+ sqrt (1 - (z/2)^2))^17
1 :+ 0
-}

module Data.Real.Constructible (
  Construct,
  deconstruct,
  fromConstruct,
  ConstructException (..)) where

import Control.Applicative ((<$>), (<*>), (<|>), empty)
import Control.Exception (Exception, ArithException (..), throw)
import Data.Complex.Generic (Complex (..))
import Data.Complex.Generic.TH (deriveComplexF)
import Data.Ratio ((%), numerator, denominator)
import Data.Typeable (Typeable)
import Math.NumberTheory.Roots (exactSquareRoot)
import Numeric.Search.Integer (search)
import Text.Read (Lexeme (..), Read (..), lexP, parens, prec, readListPrecDefault, step)
import Text.Read.Lex (numberToInteger)

data FieldShape = QShape | SqrtShape !FieldShape deriving Show

data Field k where
  Q :: Field 'QShape
  Sqrt :: !(Field k) -> !(Elt k) -> Field ('SqrtShape k)

instance Show (Field k) where
  showsPrec _ Q = showString "Q"
  showsPrec d (Sqrt k r) =
    showParen (d > 9) $ showsPrec 9 k . showString "[sqrt " . showsPrecK k 10 r . showString "]"

type family Elt (k :: FieldShape)
type instance Elt 'QShape = Rational
data SqrtElt k = SqrtZero | SqrtElt !(Elt k) !(Elt k)
type instance Elt ('SqrtShape k) = SqrtElt k

sqrtElt :: Field k -> Elt k -> Elt k -> SqrtElt k
sqrtElt k a b | isZeroK k a && isZeroK k b = SqrtZero
              | otherwise = SqrtElt a b

sqrtLift :: Field k -> Elt k -> SqrtElt k
sqrtLift k a = sqrtElt k a (zeroK k)

addK :: Field k -> Elt k -> Elt k -> Elt k
addK Q a b = a + b
addK Sqrt{} SqrtZero a = a
addK Sqrt{} a SqrtZero = a
addK (Sqrt k _) (SqrtElt a b) (SqrtElt c d) = sqrtElt k (addK k a c) (addK k b d)

mulK :: Field k -> Elt k -> Elt k -> Elt k
mulK Q a b = a * b
mulK Sqrt{} SqrtZero _ = SqrtZero
mulK Sqrt{} _ SqrtZero = SqrtZero
mulK (Sqrt k r) (SqrtElt a b) (SqrtElt c d) =
  SqrtElt (addK k (mulK k a c) (mulK k r (mulK k b d))) (addK k (mulK k a d) (mulK k b c))

sqK :: Field k -> Elt k -> Elt k
sqK Q a = a*a
sqK Sqrt{} SqrtZero = SqrtZero
sqK (Sqrt k r) (SqrtElt a b) =
  let c = mulK k a b
  in SqrtElt (addK k (sqK k a) (mulK k r (sqK k b))) (addK k c c)

subK :: Field k -> Elt k -> Elt k -> Elt k
subK Q a b = a - b
subK k@Sqrt{} SqrtZero a = negateK k a
subK Sqrt{} a SqrtZero = a
subK (Sqrt k _) (SqrtElt a b) (SqrtElt c d) = sqrtElt k (subK k a c) (subK k b d)

negateK :: Field k -> Elt k -> Elt k
negateK Q a = negate a
negateK Sqrt{} SqrtZero = SqrtZero
negateK (Sqrt k _) (SqrtElt a b) = SqrtElt (negateK k a) (negateK k b)

absK :: Field k -> Elt k -> Elt k
absK Q a = abs a
absK k a = if sgnK k a == LT then negateK k a else a

signumK :: Field k -> Elt k -> Rational
signumK Q a = signum a
signumK k a = case sgnK k a of LT -> -1; EQ -> 0; GT -> 1

divK :: Field k -> Elt k -> Elt k -> Elt k
divK Q a b = a / b
divK k a b = mulK k a (recipK k b)

recipK :: Field k -> Elt k -> Elt k
recipK Q a = recip a
recipK Sqrt{} SqrtZero = throw DivideByZero
recipK (Sqrt k r) (SqrtElt a b) =
  let c = recipK k (subK k (sqK k a) (mulK k r (sqK k b)))
  in SqrtElt (mulK k a c) (mulK k (negateK k b) c)

eqK :: Field k -> Elt k -> Elt k -> Bool
eqK Q a b = a == b
eqK Sqrt{} SqrtZero SqrtZero = True
eqK (Sqrt k _) (SqrtElt a b) (SqrtElt c d) = eqK k a c && eqK k b d
eqK Sqrt{} SqrtZero SqrtElt{} = False
eqK Sqrt{} SqrtElt{} SqrtZero = False

isZeroK :: Field k -> Elt k -> Bool
isZeroK Q a = a == 0
isZeroK Sqrt{} SqrtZero = True
isZeroK Sqrt{} SqrtElt{} = False

compareK :: Field k -> Elt k -> Elt k -> Ordering
compareK Q a b = compare a b
compareK k a b = sgnK k (subK k a b)

sgnK :: Field k -> Elt k -> Ordering
sgnK Q a = compare a 0
sgnK Sqrt{} SqrtZero = EQ
sgnK (Sqrt k r) (SqrtElt a b) = case (sgnK k a, sgnK k b) of
  (o, EQ) -> o
  (EQ, o) -> o
  (GT, GT) -> GT
  (LT, LT) -> LT
  (GT, LT) -> sgnK k (subK k (sqK k a) (mulK k r (sqK k b)))
  (LT, GT) -> sgnK k (subK k (mulK k r (sqK k b)) (sqK k a))

zeroK :: Field k -> Elt k
zeroK Q = 0
zeroK Sqrt{} = SqrtZero

fromRationalK :: Field k -> Rational -> Elt k
fromRationalK Q a = a
fromRationalK (Sqrt k _) a = sqrtLift k (fromRationalK k a)

sqrtK :: Field k -> Elt k -> Maybe (Elt k)
sqrtK Q a = (%) <$> exactSquareRoot (numerator a) <*> exactSquareRoot (denominator a)
sqrtK Sqrt{} SqrtZero = return SqrtZero
sqrtK (Sqrt k r) (SqrtElt a b)
  | isZeroK k b = sqrtLift k <$> sqrtK k a <|> SqrtElt (zeroK k) <$> sqrtK k (divK k a r)
  | otherwise = do
    n <- sqrtK k $ subK k (sqK k a) (mulK k r (sqK k b))
    let half = fromRationalK k (1 % 2)
        p = mulK k half (addK k a n)
        q = mulK k half b
        y1 = do
          c <- sqrtK k p
          return (SqrtElt c (divK k q c))
        y2 = do
          d <- sqrtK k $ divK k p r
          return (SqrtElt (divK k q d) d)
    y1 <|> y2

negateS, sqrtS :: (Int -> ShowS) -> Int -> ShowS
(-!), (+!), (*!), (/!) :: (Int -> ShowS) -> (Int -> ShowS) -> Int -> ShowS
infixl 6 +!, -!
infixl 7 *!, /!
negateS s d = showParen (d > 6) $ showChar '-' . s 7
(+!) s1 s2 d = showParen (d > 6) $ s1 6 . showString " + " . s2 7
(-!) s1 s2 d = showParen (d > 6) $ s1 6 . showString " - " . s2 7
(*!) s1 s2 d = showParen (d > 7) $ s1 7 . showChar '*' . s2 8
(/!) s1 s2 d = showParen (d > 7) $ s1 7 . showChar '/' . s2 8
sqrtS s d = showParen (d > 9) $ showString "sqrt " . s 10

mulSqrtS :: Field k -> Elt k -> Elt k -> Int -> ShowS
mulSqrtS k b r
  | eqK k b (fromRationalK k 1) = sqrtS (flip (showsPrecK k) r)
  | otherwise = flip (showsPrecK k) b *! sqrtS (flip (showsPrecK k) r)

showsPrecK :: Field k -> Int -> Elt k -> ShowS
showsPrecK Q d x
  | q == 1 = showsPrec d p
  | p < 0 = negateS (flip showsPrec (-p) /! flip showsPrec q) d
  | otherwise = (flip showsPrec p /! flip showsPrec q) d
  where
    p = numerator x
    q = denominator x
showsPrecK Sqrt{} _ SqrtZero = showChar '0'
showsPrecK (Sqrt k r) d (SqrtElt a b) = case sgnK k b of
  EQ -> showsPrecK k d a
  GT | isZeroK k a -> mulSqrtS k b r d
     | otherwise -> (flip (showsPrecK k) a +! mulSqrtS k b r) d
  LT | isZeroK k a -> negateS (mulSqrtS k (negateK k b) r) d
     | otherwise -> (flip (showsPrecK k) a -! mulSqrtS k (negateK k b) r) d

fromRatioK :: Floating a => Field k -> Elt k -> Elt k -> a
fromRatioK Q = \a b -> fromRational (a/b)
fromRatioK (Sqrt k r) = er where
  e = fromRatioK k
  s = sqrt (e r (fromRationalK k 1))
  er SqrtZero _ = 0
  er _ SqrtZero = throw DivideByZero
  er (SqrtElt a0 b0) (SqrtElt a1 b1) = case (sgnK k a, sgnK k b) of
    (_, EQ) -> e a n1
    (EQ, _) -> e b n1*s
    (GT, GT) -> x1
    (LT, LT) -> x1
    (GT, LT) -> x2
    (LT, GT) -> x2
    where
      a = subK k (mulK k a0 a1) (mulK k r (mulK k b0 b1))
      b = subK k (mulK k b0 a1) (mulK k a0 b1)
      n0 = subK k (sqK k a0) (mulK k r (sqK k b0))
      n1 = subK k (sqK k a1) (mulK k r (sqK k b1))
      x1 = e a n1 + e b n1*s
      x2 = recip $ e a n0 - e b n0*s

fromConstructK :: Floating a => Field k -> Elt k -> a
fromConstructK k a = fromRatioK k a (fromRationalK k 1)

-- |The type of constructible real numbers.
data Construct where
  C :: !(Field k) -> !(Elt k) -> Construct

deconstructK :: Field k -> Elt k -> Either Rational (Construct, Construct, Construct)
deconstructK Q a = Left a
deconstructK Sqrt{} SqrtZero = Left 0
deconstructK (Sqrt k r) (SqrtElt a b)
  | isZeroK k b = deconstructK k a
  | otherwise = Right (C k a, C k b, C k r)

{- |
Deconstruct a rational constructible number as a 'Rational', or an
irrational constructible number as a triple @(a, b, r)@ of simpler
constructible numbers representing @a + b*sqrt r@ (with @b /= 0@ and
@r > 0@).  Recursively calling 'deconstruct' on all triples will yield
a finite tree that terminates in 'Rational' leaves.

Note that two irrational constructible numbers that compare as equal
may deconstruct in different ways.
-}
deconstruct :: Construct -> Either Rational (Construct, Construct, Construct)
deconstruct (C k a) = deconstructK k a

data JoinK k1 k2 where
  JoinK :: !(Field k) -> (Elt k1 -> Elt k) -> (Elt k2 -> Elt k) -> JoinK k1 k2

joinK :: Field k1 -> Field k2 -> JoinK k1 k2
joinK Q k = JoinK k (fromRationalK k) id
joinK k Q = JoinK k id (fromRationalK k)
joinK k1 (Sqrt k2 r) = case joinK k1 k2 of
  JoinK k f1 f2 -> let r' = f2 r in case sqrtK k r' of
    Nothing ->
      let f2' SqrtZero = SqrtZero
          f2' (SqrtElt a b) = SqrtElt (f2 a) (f2 b)
      in JoinK (Sqrt k r') (sqrtLift k . f1) f2'
    Just s ->
      let f2' SqrtZero = zeroK k
          f2' (SqrtElt a b) = addK k (f2 a) (mulK k (f2 b) s)
      in JoinK k f1 f2'

instance Show Construct where
  showsPrec d (C k a) = showsPrecK k d a

instance Read Construct where
  readPrec =
    parens $
    pNum <|>
    prec 6 (pNegate <|> (step readPrec >>= pAddSub)) <|>
    prec 7 (step readPrec >>= pMulDiv) <|>
    prec 10 pSqrt
    where
      pNum = do {Number n <- lexP; maybe empty (return . fromInteger) (numberToInteger n)}
      pNegate = do {Symbol "-" <- lexP; a <- negate <$> step readPrec; return a <|> pAddSub a}
      pAddSub a = do {Symbol "+" <- lexP; b <- (a +) <$> step readPrec; return b <|> pAddSub b} <|>
                  do {Symbol "-" <- lexP; b <- (a -) <$> step readPrec; return b <|> pAddSub b}
      pMulDiv a = do {Symbol "*" <- lexP; b <- (a *) <$> step readPrec; return b <|> pMulDiv b} <|>
                  do {Symbol "/" <- lexP; b <- (a /) <$> step readPrec; return b <|> pMulDiv b}
      pSqrt = do {Ident "sqrt" <- lexP; a <- step readPrec; return (sqrt a)}
  readListPrec = readListPrecDefault

instance Eq Construct where
  C k1 a1 == C k2 a2 = case joinK k1 k2 of JoinK k f1 f2 -> eqK k (f1 a1) (f2 a2)

instance Ord Construct where
  compare (C k1 a1) (C k2 a2) = case joinK k1 k2 of JoinK k f1 f2 -> compareK k (f1 a1) (f2 a2)

instance Num Construct where
  C k1 a1 + C k2 a2 = case joinK k1 k2 of JoinK k f1 f2 -> C k (addK k (f1 a1) (f2 a2))
  C k1 a1 * C k2 a2 = case joinK k1 k2 of JoinK k f1 f2 -> C k (mulK k (f1 a1) (f2 a2))
  C k1 a1 - C k2 a2 = case joinK k1 k2 of JoinK k f1 f2 -> C k (subK k (f1 a1) (f2 a2))
  negate (C k x) = C k (negateK k x)
  abs (C k x) = C k (absK k x)
  signum (C k x) = C Q (signumK k x)
  fromInteger = C Q . fromInteger

instance Fractional Construct where
  C k1 a1 / C k2 a2 = case joinK k1 k2 of JoinK k f1 f2 -> C k (divK k (f1 a1) (f2 a2))
  recip (C k a) = C k (recipK k a)
  fromRational = C Q

-- |The type of exceptions thrown by impossible 'Construct' operations.
data ConstructException =
  -- |'toRational' was given an irrational constructible number.
  ConstructIrrational |
  -- |'sqrt' was given a negative constructible number.
  ConstructSqrtNegative |
  -- |'**' was given an exponent that is not a dyadic rational, or a transcendental function was called.
  Unconstructible String
  deriving (Eq, Ord, Typeable)

instance Show ConstructException where
  showsPrec _ ConstructIrrational = showString "cannot convert irrational Construct to rational"
  showsPrec _ ConstructSqrtNegative = showString "Construct sqrt: negative argument"
  showsPrec _ (Unconstructible s) = showString s . showString " is not constructible"

instance Exception ConstructException

{- |
This partial 'Floating' instance only supports 'sqrt' and '**' where
the exponent is a dyadic rational.  Passing a negative number to
'sqrt' will throw the 'ConstructSqrtNegative' exception.  All other
operations will throw the 'Unconstructible' exception.
-}
instance Floating Construct where
  sqrt (C k a)
    | sgnK k a == LT = throw ConstructSqrtNegative
    | otherwise = case sqrtK k a of
        Nothing -> C (Sqrt k a) (SqrtElt (zeroK k) (fromRationalK k 1))
        Just b -> C k b
  pi = throw (Unconstructible "pi")
  exp = throw (Unconstructible "exp")
  log = throw (Unconstructible "log")
  a ** b = go (numerator b') (denominator b') where
    b' = toRational b
    go p q = let (n, p') = divMod p q in a^^n*go' p' q
    go' 0 _ = 1
    go' p q = case divMod q 2 of
      (q', 0) -> sqrt (go p q')
      _ -> throw (Unconstructible "(** non-dyadic-rational)")
  logBase = throw (Unconstructible "logBase")
  sin = throw (Unconstructible "sin")
  tan = throw (Unconstructible "tan")
  cos = throw (Unconstructible "cos")
  asin = throw (Unconstructible "asin")
  atan = throw (Unconstructible "atan")
  acos = throw (Unconstructible "acos")
  sinh = throw (Unconstructible "sinh")
  tanh = throw (Unconstructible "tanh")
  cosh = throw (Unconstructible "cosh")
  asinh = throw (Unconstructible "asinh")
  atanh = throw (Unconstructible "atanh")
  acosh = throw (Unconstructible "acosh")

{- |
This 'Real' instance only supports 'toRational' on numbers that are in
fact rational.  'toRational' on an irrational number will throw the
'ConstructIrrational' exception.
-}
instance Real Construct where
  toRational = either id (\_ -> throw ConstructIrrational) . deconstruct

instance RealFrac Construct where
  properFraction (C Q x) = (m, C Q y) where (m, y) = properFraction x
  properFraction x = (fromInteger m, x - fromInteger m)
    where m = search ((> x) . fromInteger) - 1

instance Enum Construct where
  succ = (+ 1)
  pred = subtract 1
  toEnum = fromIntegral
  fromEnum = fromInteger . truncate
  enumFrom n = n `seq` (n : enumFrom (n + 1))
  enumFromThen n m = n `seq` m `seq` (n : enumFromThen m (m + m - n))
  enumFromTo n m = takeWhile (<= m) (enumFrom n)
  enumFromThenTo e1 e2 e3 = takeWhile predicate (enumFromThen e1 e2) where
    predicate | e2 >= e1 = (<= e3)
              | otherwise = (>= e3)

mk :: a -> a -> Complex a
mk = (:+)

toPair :: Complex a -> (a, a)
toPair (x :+ y) = (x, y)

deriveComplexF ''Complex ''Construct 'mk 'toPair

{- |
Evaluate a floating-point approximation for a constructible number.

To improve numerical stability, addition of numbers with different
signs is avoided using quadratic conjugation.

>>> fromConstruct $ sum (map sqrt [7, 14, 39, 70, 72, 76, 85]) - sum (map sqrt [13, 16, 46, 55, 67, 73, 79])
1.8837969820815017e-19
-}
fromConstruct :: Floating a => Construct -> a
fromConstruct (C k a) = fromConstructK k a

max_stars_repo_name: andersk/haskell-constructible

max_issues_repo_name: andersk/haskell-constructible

max_forks_repo_name: andersk/haskell-constructible


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE FlexibleContexts #-}

module Taiji.Pipeline.SC.ATACSeq.Functions.Feature.Peak
    ( mkFeatMat
    , mkPeakMat
    , findPeaks
    , mergePeaks
    , getPeakEnrichment
    , mkCellClusterBed
    , computePeakRAS
    ) where

import Control.Arrow (first)
import           Bio.Pipeline
import qualified Data.HashSet as S
import Bio.Data.Bed hiding (NarrowPeak)
import qualified Bio.Data.Bed as Bed
import Data.Conduit.Internal (zipSinks, zipSources)
import qualified Data.Text as T
import qualified Data.ByteString.Char8 as B
import qualified Data.Matrix as Mat
import Data.Singletons.Prelude (SingI)
import Shelly hiding (FilePath)
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U
import Control.DeepSeq (force)
import Data.Conduit.Zlib (gzip)

import Statistics.Distribution (quantile)
import Statistics.Distribution.ChiSquared (chiSquared)

import Taiji.Prelude
import Taiji.Pipeline.SC.ATACSeq.Types
import Taiji.Pipeline.SC.ATACSeq.Functions.Utils
import Taiji.Pipeline.SC.ATACSeq.Functions.QC (streamTN5Insertion)
import qualified Taiji.Utils.DataFrame as DF
import Taiji.Utils

-- | Make the read count matrix.
mkPeakMat :: SCATACSeqConfig config
          => FilePath
          -> SCATACSeq S (TagsAligned, File '[Gzip] 'NarrowPeak)
          -> ReaderT config IO (SCATACSeq S (File '[Gzip] 'Matrix))
mkPeakMat dir input = do
    passedQC <- getQCFunction
    let output = printf "%s/%s_rep%d_peak.mat.gz" dir (T.unpack $ input^.eid)
            (input^.replicates._1)
    input & replicates.traverse.files %%~ liftIO . (\(tagFl, regionFl) -> do
        regions <- runResourceT $ runConduit $
            streamBedGzip (regionFl^.location) .| sinkList :: IO [BED3]
        runResourceT $ runConduit $ streamTN5Insertion passedQC tagFl .|
            mapC (\xs -> (fromJust $ head xs^.name, xs)) .| mkCountMat regions .|
            sinkRows' (length regions) (fromJust . packDecimal) output
        return $ emptyFile & location .~ output )

-- | Make the read count matrix.
mkFeatMat :: SCATACSeqConfig config
          => FilePath
          -> SCATACSeq S (TagsAligned, File '[Gzip] 'NarrowPeak)
          -> ReaderT config IO ( SCATACSeq S
              ( File '[RowName, Gzip] 'Tsv
              , File '[ColumnName, Gzip] 'Tsv
              , File '[Gzip] 'Matrix ))
mkFeatMat dir input = do
    passedQC <- getQCFunction
    let output = printf "%s/%s_rep%d_peak.mat.gz" dir (T.unpack $ input^.eid)
            (input^.replicates._1)
        rownames = printf "%s/%s_rep%d_rownames.txt.gz" dir (T.unpack $ input^.eid)
            (input^.replicates._1)
        features = printf "%s/%s_rep%d_features.txt.gz" dir (T.unpack $ input^.eid)
            (input^.replicates._1)
        bcPrefix = B.pack $ T.unpack (input^.eid) <> "_" <> show (input^.replicates._1) <> "+"
    input & replicates.traverse.files %%~ liftIO . (\(tagFl, regionFl) -> do
        regions <- runResourceT $ runConduit $
            streamBedGzip (regionFl^.location) .| sinkList :: IO [BED3]
        let rowSink = mapC f .| unlinesAsciiC .| gzip .| sinkFile rownames
              where
                f (nm, xs) = force $ nm <> "\t" <> fromJust (packDecimal $ foldl1' (+) $ map snd xs)
            colSink = do
                vec <- colSumC $ length regions
                let bs = B.unlines $ zipWith showBed regions $ U.toList vec
                yield bs .| gzip .| sinkFile features
            showBed (BED3 chr s e) x = B.concat
                [ chr, ":", fromJust $ packDecimal s, "-"
                , fromJust $ packDecimal e, "\t", fromJust $ packDecimal x]
            sink = (,,) <$> ZipSink (sinkRows' (length regions) (fromJust . packDecimal) output)
                <*> ZipSink rowSink <*> ZipSink colSink
        _ <- runResourceT $ runConduit $ streamTN5Insertion passedQC tagFl .|
            mapC (\xs -> (fromJust $ head xs^.name, xs)) .|
            mkCountMat regions .| mapC (first (bcPrefix <>)) .|
            getZipSink sink
        return ( location .~ rownames $ emptyFile
               , location .~ features $ emptyFile
               , location .~ output $ emptyFile )
        )

-- | Call Peaks for aggregated files
findPeaks :: (SingI tags, SCATACSeqConfig config)
          => FilePath
          -> CallPeakOpts
          -> (B.ByteString, File tags 'Bed)
          -> ReaderT config IO (B.ByteString, File '[Gzip] 'NarrowPeak)
findPeaks prefix opts (cName, bedFl) = do
    tmpdir <- asks _scatacseq_tmp_dir
    dir <- asks _scatacseq_output_dir >>= getPath . (<> asDir prefix)
    let output = dir ++ "/" ++ B.unpack cName ++ ".narrowPeak.gz"
    asks _scatacseq_blacklist >>= \case
        Nothing -> liftIO $ do
            r <- callPeaks output bedFl Nothing opts
            return (cName, r)
        Just blacklist -> liftIO $ withTemp tmpdir $ \tmp -> do
            _ <- callPeaks tmp bedFl Nothing opts
            blackRegions <- readBed blacklist :: IO [BED3]
            let bedTree = bedToTree const $ map (\x -> (x, ())) blackRegions
            runResourceT $ runConduit $
                (streamBedGzip tmp :: ConduitT () Bed.NarrowPeak (ResourceT IO) ()) .|
                filterC (not . isIntersected bedTree) .| sinkFileBedGzip output
            return (cName, location .~ output $ emptyFile)

-- | Merge peaks
mergePeaks :: SCATACSeqConfig config
           => FilePath
           -> [(B.ByteString, File '[Gzip] 'NarrowPeak)]
           -> ReaderT config IO (Maybe (File '[Gzip] 'NarrowPeak))
mergePeaks _ [] = return Nothing
mergePeaks dir input = do
    tmpdir <- asks _scatacseq_tmp_dir
    let output = dir <> "/merged.narrowPeak.gz"
    liftIO $ withTemp tmpdir $ \tmp1 -> withTemp tmpdir $ \tmp2 -> do
        runResourceT $ runConduit $ mapM_ (streamBedGzip . (^._2.location)) input .|
            mapC resize .| sinkFileBed tmp1
        shelly $ escaping False $ bashPipeFail bash_ "cat" $
            [T.pack tmp1, "|", "sort", "-k1,1", "-k2,2n", "-k3,3n", ">", T.pack tmp2]
        let source = streamBed tmp2 .|
                mergeSortedBedWith iterativeMerge .| concatC
        runResourceT $ runConduit $ zipSources (iterateC succ (0 :: Int)) source .|
            mapC (\(i, p) -> name .~ Just ("p" <> B.pack (show i)) $ p) .|
            sinkFileBedGzip output
    return $ Just $ location .~ output $ emptyFile
  where
    iterativeMerge [] = []
    iterativeMerge peaks = bestPeak : iterativeMerge rest
      where
        rest = filter (\x -> sizeOverlapped x bestPeak == 0) peaks
        bestPeak = maximumBy (comparing (fromJust . (^.npPvalue))) peaks
    resize pk = chromStart .~ max 0 (summit - halfWindowSize) $
        chromEnd .~ summit + halfWindowSize$
        npPeak .~ Just halfWindowSize $ pk
      where
        summit = pk^.chromStart + fromJust (pk^.npPeak)
    halfWindowSize = 200

-- | Extract BEDs for each cluster.
mkCellClusterBed :: SCATACSeqConfig config
                 => SCATACSeq S ( File '[NameSorted, Gzip] 'Bed
                                , [CellCluster] )  -- ^ clusters
                 -> ReaderT config IO
                    (SCATACSeq S [(B.ByteString, File '[Gzip] 'Bed, Int)])
mkCellClusterBed input = do
    let idRep = asDir $ "/Bed/" <> T.unpack (input^.eid) <>
            "_rep" <> show (input^.replicates._1)
    dir <- asks _scatacseq_output_dir >>= getPath . (<> idRep)
    input & replicates.traverse.files %%~ liftIO . ( \(bed, cs) -> do
        let sinks = sequenceConduits $ flip map cs $ \CellCluster{..} -> do
                let output = dir ++ "/" ++ B.unpack _cluster_name ++ ".bed.gz"
                    cells = S.fromList $ map _cell_barcode _cluster_member
                    fl = location .~ output $ emptyFile
                (_, depth) <- filterC (f cells) .|
                    zipSinks (sinkFileBedGzip output) lengthC
                return (_cluster_name, fl, depth)
        runResourceT $ runConduit $
            streamBedGzip (bed^.location) .| sinks )
  where
    f :: S.HashSet B.ByteString -> BED -> Bool
    f cells x = fromJust (x^.name) `S.member` cells

-- | Get signal enrichment for each peak
getPeakEnrichment :: SCATACSeqConfig config
                  => ( [(B.ByteString, File '[Gzip] 'NarrowPeak)]
                     , Maybe (File '[Gzip] 'NarrowPeak) )
                  -> ReaderT config IO (File '[] 'Tsv, File '[] 'Tsv)
getPeakEnrichment (peaks, Just refPeak) = do
    dir <- asks _scatacseq_output_dir >>= getPath . (<> "/Feature/Peak/")
    let output1 = dir <> "peak_signal.tsv"
        output2 = dir <> "peak_pvalue.tsv"
    liftIO $ do
        list <- runResourceT $ runConduit $
            streamBedGzip (refPeak^.location) .| sinkList
        let col = map toString list
        (row, val) <- fmap unzip $ forM peaks $ \(nm, p) -> do
            peak <- runResourceT $ runConduit $ streamBedGzip (p^.location) .| sinkList
            return (T.pack $ B.unpack nm, getValue peak list)
        let (signal, pval) = DF.unzip $ DF.transpose $ DF.mkDataFrame row col val
        DF.writeTable output1 (T.pack . show) signal
        DF.writeTable output2 (T.pack . show) pval
        return (location .~ output1 $ emptyFile, location .~ output2 $ emptyFile)
  where
    toString x = T.pack (B.unpack $ x^.chrom) <> ":" <>
        T.pack (show $ x^.chromStart) <> "-" <> T.pack (show $ x^.chromEnd)
    getValue :: [Bed.NarrowPeak]  -- ^ query
            -> [BED3]  -- ^ Peak list
            -> [(Double, Double)]
    getValue query peakList = runIdentity $ runConduit $
        yieldMany peakList .| intersectBedWith f query .| sinkList
      where
        f _ [] = (0, 0)
        f _ [x] = (x^.npSignal, quantile (chiSquared 1) $ 10**(negate $ fromJust $ x^.npPvalue))
        f _ _ = undefined
getPeakEnrichment _ = undefined

-- | Compute the relative accessibility score.
computePeakRAS :: SCATACSeqConfig config
               => FilePath
               -> ( Maybe (File '[Gzip] 'NarrowPeak)
                  , [SCATACSeq S (File '[Gzip] 'Matrix)] )
               -> ReaderT config IO (Maybe (FilePath, FilePath, FilePath))
computePeakRAS _ (Nothing, _) = return Nothing
computePeakRAS prefix (peakFl, inputs) = do
    dir <- asks ((<> asDir prefix) . _scatacseq_output_dir) >>= getPath
    let output1 = dir <> "relative_accessibility_scores.tsv"
        output2 = dir <> "cell_specificity_score.tsv"
        output3 = dir <> "/cell_specificity_pvalue.tsv"
    liftIO $ do
        peaks <- fmap (map mkName) $ runResourceT $ runConduit $
            streamBedGzip (fromJust peakFl^.location) .| sinkList

        (names, cols) <- fmap unzip $ forM inputs $ \input -> do
            vec <- computeRAS (input^.replicates._2.files.location)
            return (input^.eid, V.convert vec)

        let ras = DF.map (*2.5) $ DF.fromMatrix peaks names $ Mat.fromColumns cols
            ss = computeSS $ DF.map (logBase 2 . (+1)) ras
        DF.writeTable output1 (T.pack . show) ras
        DF.writeTable output2 (T.pack . show) ss
        cdf <- computeCDF $ DF.map (logBase 2 . (+1)) ras
        DF.writeTable output3 (T.pack . show) $ DF.map (lookupP cdf) ss
        return $ Just (output1, output2, output3)
  where
    mkName :: BED3 -> T.Text
    mkName p = T.pack $ B.unpack (p^.chrom) <> ":" <> show (p^.chromStart) <>
        "-" <> show (p^.chromEnd)
    lookupP (vec, res, n) x | p == 0 = 1 / n
                            | otherwise = p
      where
        p = vec V.! i
        i = min (V.length vec - 1) $ truncate $ x / res

max_stars_repo_name: Taiji-pipeline/Taiji-scATAC-seq

max_issues_repo_name: Taiji-pipeline/Taiji-scATAC-seq

max_forks_repo_name: Taiji-pipeline/Taiji-scATAC-seq


text:
module Plot (
    plot,
    plotDensity
) where

import Graphics.Rendering.Chart
import Graphics.Rendering.Chart.Backend.Cairo
import Control.Lens
import Data.Colour
import Data.Colour.Names
import Data.Default.Class

import Plot.Utils
import Statistics.Types
import Statistics.Distribution
import qualified Data.Vector.Unboxed as V

width, height :: Int
width = 800
height = 400
file_opts :: FileOptions
file_opts = fo_size .~ (width, height)
          $ fo_format .~ SVG
          $ def

plot :: Sample -> String -> FilePath -> IO (PickFn ())
plot vs title fn = renderableToFile file_opts fn (toRenderable layout)
  where
    n = V.length vs

    tv, ev, dv :: [(Double, Double)]
    tv = [(fromIntegral i, i /. n) | i <- [1 .. n]]
    ev = zip (map fromIntegral [1 .. n]) (V.toList vs)
    dv = zipWith (\ (x, t) (_, e) -> (x, t - e)) tv ev

    layout :: LayoutLR Double Double Double
    layout = layoutlr_title      .~ title
           $ layoutlr_background .~ solidFillStyle (opaque white)
           $ layoutlr_plots      .~ [ Left (toPlot theoretical),
                                      Left (toPlot empirical),
                                      Right (toPlot difference)]
           $ layoutlr_right_axis .~ yaxis
           $ layoutlr_foreground .~ opaque black
           $ def

    theoretical :: PlotLines Double Double
    theoretical = theoreticalWith tv

    empirical :: PlotLines Double Double
    empirical = empiricalWith ev

    difference :: PlotLines Double Double
    difference = differenceWith dv

yaxis :: LayoutAxis Double
yaxis  = laxis_title    .~ "theoretical - empirical"
       $ laxis_override .~ (axisGridHide . rStyleTicks)
       $ def

rStyleTicks     :: AxisData x -> AxisData x
rStyleTicks ad  = ad{ _axis_ticks = map invert_tick (_axis_ticks ad) }
    where
        invert_tick :: (x,Double) -> (x, Double)
        invert_tick (x, t) = (x, -t)

plotDensity :: (ContDistr d) => d -> Sample -> String -> FilePath -> IO (PickFn ())
plotDensity gd vs title fn = renderableToFile file_opts fn (toRenderable layout)
  where

    n  = V.length vs
    mx = V.maximum vs

    -- divide domain into bins
    bins :: Double
    bins = fromIntegral $ max (n `div` 100) 1000

    xs, ys :: [Double]
    xs = map (* (mx / bins)) [0 .. bins]
    -- count how many at each level
    ys =  let counts = count (V.toList vs) xs
              scale'  = fromIntegral n * (mx / bins)
          in  map ((/ scale') . fromIntegral) counts

    layout :: LayoutLR Double Double Double
    layout = stdLayout title
               [ Left (toPlot (  empiricalWith (zip xs ys)))
               , Left (toPlot (theoreticalWith [(x, density gd x) | x <- xs]) )]

(/.) :: (Real a, Real b, Fractional c) => a -> b -> c
(/.) x y = fromRational $ toRational x / toRational y

count :: [Double] -> [Double] -> [Int]
count samples breaks = map length (foldr go (:[]) breaks samples)
  where
    go x r s = let (ys, zs) = break (>= x) s in ys : r zs

max_stars_repo_name: finlay/random-dist

max_issues_repo_name: finlay/random-dist

max_forks_repo_name: finlay/random-dist


text:
module Statistics where

import qualified Data.Vector       as V
import qualified Statistics.Sample as S

data Description = Description
  { min'   :: Double
  , mean   :: Double
  , stdDev :: Double
  , max'   :: Double
  } deriving (Show)

describe :: V.Vector Double -> Description
describe xs =
  Description
  {min' = minimum xs, mean = S.mean xs, stdDev = S.stdDev xs, max' = maximum xs}

max_stars_repo_name: anthony-khong/data-analysis

max_issues_repo_name: anthony-khong/data-analysis

max_forks_repo_name: anthony-khong/data-analysis


text:
-- This program is intentionally poorly documented; fixing the
-- documentation and actually creating the program are left as exercises
-- for the reader.

import qualified Data.List as DL;
import qualified Data.Function as DF;
import qualified Numeric.FFT as FT;
import Data.Complex;

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

windowSize :: Int;
windowSize = 65536;

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

type Giraffe = [(Int, Pitch)];  -- LIST OF GRAPH COORDINATES
type Pitch = Double;
type Sample = Double;
type Frequency = Double;
type Multiplier = Double;
type FFTOutput = [Complex Double];
data X11Winda = ExerciseForTheReader;
data WaveType = Sin | Cos;

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

-- | createNewWinda creates a new X11 window and returns the tag of this
-- X11 window.
createNewWinda :: IO X11Winda;
createNewWinda = return ExerciseForTheReader;

fft :: [Sample] -> FFTOutput;
fft = FT.fft . map cis;

-- | For all FFT outputs k, fftToPitch equals the loudest pitch of k.
-- Implementing fftToPitch is left as an exercise for the reader.
-- Shut up and hack!
fftToPitch :: FFTOutput -> Pitch;
fftToPitch k = 0;

-- | readMicData is a stream of microphone samples.
-- Implementing readMicData is left as an exercise for the reader.
-- The author recommends using the SDL library or an operating system-
-- specific approach.
readMicData :: IO [Sample];
readMicData = return $ map (\a -> a * (sin a)) [1.0,1.1..];

-- | For all natural numbers n, readSamp n equals n microphone samples.
readSamp :: Int -> IO [Sample];
readSamp k = readMicData >>= return . take k;

-- | For all [Sample] n, seqToPitch n equals the output of the sequence-
-- to-pitch function on the input n.
seqToPitch :: [Sample] -> Pitch;
seqToPitch = fftToPitch . fft;

-- | genCoords yields the coordinates of the graph.
genCoords :: [Pitch] -> Giraffe;
genCoords x = zip [0..length x - 1] x;

-- | For all X11Winda k, for all Giraffe g, genGraph k g draws
-- g to k, returning the pitch/y values of g.
-- Implementing genGraph is left as an exercise for the reader.
-- The author recommends using Chart or creating a custom library.
genGraph :: X11Winda -> Giraffe -> IO [Pitch];
genGraph w k = putStrLn "TODO: IMPLEMENT GRAPHING." >> return (map snd k);

-- | mane is the main program loop.
mane :: X11Winda -> [Double] -> IO ();
mane w k = readSamp windowSize >>= doTheGraph w >>= loop w
  where
  doTheGraph a b = genGraph a (genCoords $ append k $ seqToPitch b)
  append a b = a ++ [b]
  loop a b = mane a (take 50 b);

-- | main is the program's entry point.  main initialises mane.
main :: IO ();
main = createNewWinda >>= \winda -> mane winda [];

max_stars_repo_name: varikvalefor/tthheeppaarrttyy

max_issues_repo_name: varikvalefor/tthheeppaarrttyy

max_forks_repo_name: varikvalefor/tthheeppaarrttyy


text:
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# LANGUAGE FlexibleInstances, OverlappingInstances, ScopedTypeVariables,
    ViewPatterns #-}
module Tests.Distribution (tests) where

import Control.Applicative ((<$), (<$>), (<*>))
import Data.Binary (Binary, decode, encode)
import Data.List (find)
import Data.Typeable (Typeable)
import Statistics.Distribution
import Statistics.Distribution.Beta (BetaDistribution, betaDistr)
import Statistics.Distribution.Binomial (BinomialDistribution, binomial)
import Statistics.Distribution.CauchyLorentz
import Statistics.Distribution.ChiSquared (ChiSquared, chiSquared)
import Statistics.Distribution.Exponential (ExponentialDistribution, exponential)
import Statistics.Distribution.FDistribution (FDistribution, fDistribution)
import Statistics.Distribution.Gamma (GammaDistribution, gammaDistr)
import Statistics.Distribution.Geometric
import Statistics.Distribution.Hypergeometric
import Statistics.Distribution.Normal (NormalDistribution, normalDistr)
import Statistics.Distribution.Poisson (PoissonDistribution, poisson)
import Statistics.Distribution.StudentT
import Statistics.Distribution.Transform (LinearTransform, linTransDistr)
import Statistics.Distribution.Uniform (UniformDistribution, uniformDistr)
import Test.Framework (Test, testGroup)
import Test.Framework.Providers.QuickCheck2 (testProperty)
import Test.QuickCheck as QC
import Test.QuickCheck.Monadic as QC
import Tests.ApproxEq (ApproxEq(..))
import Tests.Helpers (T(..), testAssertion, typeName)
import Tests.Helpers (monotonicallyIncreasesIEEE)
import Text.Printf (printf)
import qualified Control.Exception as E
import qualified Numeric.IEEE as IEEE

-- | Tests for all distributions
tests :: Test
tests = testGroup "Tests for all distributions"
  [ contDistrTests (T :: T BetaDistribution        )
  , contDistrTests (T :: T CauchyDistribution      )
  , contDistrTests (T :: T ChiSquared              )
  , contDistrTests (T :: T ExponentialDistribution )
  , contDistrTests (T :: T GammaDistribution       )
  , contDistrTests (T :: T NormalDistribution      )
  , contDistrTests (T :: T UniformDistribution     )
  , contDistrTests (T :: T StudentT                )
  , contDistrTests (T :: T (LinearTransform StudentT) )
  , contDistrTests (T :: T FDistribution           )

  , discreteDistrTests (T :: T BinomialDistribution       )
  , discreteDistrTests (T :: T GeometricDistribution      )
  , discreteDistrTests (T :: T GeometricDistribution0     )
  , discreteDistrTests (T :: T HypergeometricDistribution )
  , discreteDistrTests (T :: T PoissonDistribution        )

  , unitTests
  ]

----------------------------------------------------------------
-- Tests
----------------------------------------------------------------

-- Tests for continous distribution
contDistrTests :: (Param d, ContDistr d, QC.Arbitrary d, Typeable d, Show d, Binary d, Eq d) => T d -> Test
contDistrTests t = testGroup ("Tests for: " ++ typeName t) $
  cdfTests t ++
  [ testProperty "PDF sanity"              $ pdfSanityCheck     t
  , testProperty "Quantile is CDF inverse" $ quantileIsInvCDF   t
  , testProperty "quantile fails p<0||p>1" $ quantileShouldFail t
  , testProperty "log density check"       $ logDensityCheck    t
  ]

-- Tests for discrete distribution
discreteDistrTests :: (Param d, DiscreteDistr d, QC.Arbitrary d, Typeable d, Show d, Binary d, Eq d) => T d -> Test
discreteDistrTests t = testGroup ("Tests for: " ++ typeName t) $
  cdfTests t ++
  [ testProperty "Prob. sanity"         $ probSanityCheck       t
  , testProperty "CDF is sum of prob."  $ discreteCDFcorrect    t
  , testProperty "Discrete CDF is OK"   $ cdfDiscreteIsCorrect  t
  , testProperty "log probabilty check" $ logProbabilityCheck   t
  ]

-- Tests for distributions which have CDF
cdfTests :: (Param d, Distribution d, QC.Arbitrary d, Show d, Binary d, Eq d) => T d -> [Test]
cdfTests t =
  [ testProperty "C.D.F. sanity"        $ cdfSanityCheck         t
  , testProperty "CDF limit at +inf"    $ cdfLimitAtPosInfinity  t
  , testProperty "CDF limit at -inf"    $ cdfLimitAtNegInfinity  t
  , testProperty "CDF at +inf = 1"      $ cdfAtPosInfinity       t
  , testProperty "CDF at -inf = 1"      $ cdfAtNegInfinity       t
  , testProperty "CDF is nondecreasing" $ cdfIsNondecreasing     t
  , testProperty "1-CDF is correct"     $ cdfComplementIsCorrect t
  , testProperty "Binary OK"            $ p_binary t
  ]

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

-- CDF is in [0,1] range
cdfSanityCheck :: (Distribution d) => T d -> d -> Double -> Bool
cdfSanityCheck _ d x = c >= 0 && c <= 1
  where c = cumulative d x

-- CDF never decreases
cdfIsNondecreasing :: (Distribution d) => T d -> d -> Double -> Double -> Bool
cdfIsNondecreasing _ d = monotonicallyIncreasesIEEE $ cumulative d

-- cumulative d +∞ = 1
cdfAtPosInfinity :: (Param d, Distribution d) => T d -> d -> Bool
cdfAtPosInfinity _ d
  = cumulative d (1/0) == 1

-- cumulative d - ∞ = 0
cdfAtNegInfinity :: (Param d, Distribution d) => T d -> d -> Bool
cdfAtNegInfinity _ d
  = cumulative d (-1/0) == 0

-- CDF limit at +∞ is 1
cdfLimitAtPosInfinity :: (Param d, Distribution d) => T d -> d -> Property
cdfLimitAtPosInfinity _ d =
  okForInfLimit d ==> counterexample ("Last elements: " ++ show (drop 990 probs))
                    $ Just 1.0 == (find (>=1) probs)
  where
    probs = take 1000 $ map (cumulative d) $ iterate (*1.4) 1000

-- CDF limit at -∞ is 0
cdfLimitAtNegInfinity :: (Param d, Distribution d) => T d -> d -> Property
cdfLimitAtNegInfinity _ d =
  okForInfLimit d ==> counterexample ("Last elements: " ++ show (drop 990 probs))
                    $ case find (< IEEE.epsilon) probs of
                        Nothing -> False
                        Just p  -> p >= 0
  where
    probs = take 1000 $ map (cumulative d) $ iterate (*1.4) (-1)

-- CDF's complement is implemented correctly
cdfComplementIsCorrect :: (Distribution d) => T d -> d -> Double -> Bool
cdfComplementIsCorrect _ d x = (eq 1e-14) 1 (cumulative d x + complCumulative d x)

-- CDF for discrete distribution uses <= for comparison
cdfDiscreteIsCorrect :: (DiscreteDistr d) => T d -> d -> Property
cdfDiscreteIsCorrect _ d
  = counterexample (unlines badN)
  $ null badN
  where
    -- We are checking that:
    --
    -- > CDF(i) - CDF(i-e) = P(i)
    --
    -- Apporixmate equality is tricky here. Scale is set by maximum
    -- value of CDF and probability. Case when all proabilities are
    -- zero should be trated specially.
    badN = [ printf "N=%3i    p[i]=%g\tp[i+1]=%g\tdP=%g\trelerr=%g" i p p1 dp ((p1-p-dp) / max p1 dp)
           | i <- [0 .. 100]
           , let p      = cumulative d $ fromIntegral i - 1e-6
                 p1     = cumulative d $ fromIntegral i
                 dp     = probability d i
                 relerr = ((p1 - p) - dp) / max p1 dp
           ,  not (p == 0 && p1 == 0 && dp == 0)
           && relerr > 1e-14
           ]

logDensityCheck :: (ContDistr d) => T d -> d -> Double -> Property
logDensityCheck _ d x
  = counterexample (printf "density    = %g" p)
  $ counterexample (printf "logDensity = %g" logP)
  $ counterexample (printf "log p      = %g" (log p))
  $ counterexample (printf "eps        = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))
  $ or [ p == 0     && logP == (-1/0)
       , p < 1e-308 && logP < 609
       , eq 1e-14 (log p) logP
       ]
  where
    p    = density d x
    logP = logDensity d x

-- PDF is positive
pdfSanityCheck :: (ContDistr d) => T d -> d -> Double -> Bool
pdfSanityCheck _ d x = p >= 0
  where p = density d x

-- Quantile is inverse of CDF
quantileIsInvCDF :: (Param d, ContDistr d) => T d -> d -> Double -> Property
quantileIsInvCDF _ d (snd . properFraction -> p) =
  p > 0 && p < 1  ==> ( counterexample (printf "Quantile     = %g" q )
                      $ counterexample (printf "Probability  = %g" p )
                      $ counterexample (printf "Probability' = %g" p')
                      $ counterexample (printf "Error        = %e" (abs $ p - p'))
                      $ abs (p - p') < invQuantilePrec d
                      )
  where
    q  = quantile   d p
    p' = cumulative d q

-- Test that quantile fails if p<0 or p>1
quantileShouldFail :: (ContDistr d) => T d -> d -> Double -> Property
quantileShouldFail _ d p =
  p < 0 || p > 1 ==> QC.monadicIO $ do r <- QC.run $ E.catch
                                              (False <$ (return $! quantile d p))
                                              (\(_ :: E.SomeException) -> return True)
                                       QC.assert r

-- Probability is in [0,1] range
probSanityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Bool
probSanityCheck _ d x = p >= 0 && p <= 1
  where p = probability d x

-- Check that discrete CDF is correct
discreteCDFcorrect :: (DiscreteDistr d) => T d -> d -> Int -> Int -> Property
discreteCDFcorrect _ d a b
  = counterexample (printf "CDF   = %g" p1)
  $ counterexample (printf "Sum   = %g" p2)
  $ counterexample (printf "Delta = %g" (abs (p1 - p2)))
  $ abs (p1 - p2) < 3e-10
  -- Avoid too large differeneces. Otherwise there is to much to sum
  --
  -- Absolute difference is used guard againist precision loss when
  -- close values of CDF are subtracted
  where
    n  = min a b
    m  = n + (abs (a - b) `mod` 100)
    p1 = cumulative d (fromIntegral m + 0.5) - cumulative d (fromIntegral n - 0.5)
    p2 = sum $ map (probability d) [n .. m]

logProbabilityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Property
logProbabilityCheck _ d x
  = counterexample (printf "probability    = %g" p)
  $ counterexample (printf "logProbability = %g" logP)
  $ counterexample (printf "log p          = %g" (log p))
  $ counterexample (printf "eps            = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))
  $ or [ p == 0     && logP == (-1/0)
       , p < 1e-308 && logP < 609
       , eq 1e-14 (log p) logP
       ]
  where
    p    = probability d x
    logP = logProbability d x

p_binary :: (Eq a, Show a, Binary a) => T a -> a -> Bool
p_binary _ a = a == (decode . encode) a

----------------------------------------------------------------
-- Arbitrary instances for ditributions
----------------------------------------------------------------

instance QC.Arbitrary BinomialDistribution where
  arbitrary = binomial <$> QC.choose (1,100) <*> QC.choose (0,1)
instance QC.Arbitrary ExponentialDistribution where
  arbitrary = exponential <$> QC.choose (0,100)
instance QC.Arbitrary GammaDistribution where
  arbitrary = gammaDistr <$> QC.choose (0.1,10) <*> QC.choose (0.1,10)
instance QC.Arbitrary BetaDistribution where
  arbitrary = betaDistr <$> QC.choose (1e-3,10) <*> QC.choose (1e-3,10)
instance QC.Arbitrary GeometricDistribution where
  arbitrary = geometric <$> QC.choose (0,1)
instance QC.Arbitrary GeometricDistribution0 where
  arbitrary = geometric0 <$> QC.choose (0,1)
instance QC.Arbitrary HypergeometricDistribution where
  arbitrary = do l <- QC.choose (1,20)
                 m <- QC.choose (0,l)
                 k <- QC.choose (1,l)
                 return $ hypergeometric m l k
instance QC.Arbitrary NormalDistribution where
  arbitrary = normalDistr <$> QC.choose (-100,100) <*> QC.choose (1e-3, 1e3)
instance QC.Arbitrary PoissonDistribution where
  arbitrary = poisson <$> QC.choose (0,1)
instance QC.Arbitrary ChiSquared where
  arbitrary = chiSquared <$> QC.choose (1,100)
instance QC.Arbitrary UniformDistribution where
  arbitrary = do a <- QC.arbitrary
                 b <- QC.arbitrary `suchThat` (/= a)
                 return $ uniformDistr a b
instance QC.Arbitrary CauchyDistribution where
  arbitrary = cauchyDistribution
                <$> arbitrary
                <*> ((abs <$> arbitrary) `suchThat` (> 0))
instance QC.Arbitrary StudentT where
  arbitrary = studentT <$> ((abs <$> arbitrary) `suchThat` (>0))
instance QC.Arbitrary (LinearTransform StudentT) where
  arbitrary = studentTUnstandardized
           <$> ((abs <$> arbitrary) `suchThat` (>0))
           <*> ((abs <$> arbitrary))
           <*> ((abs <$> arbitrary) `suchThat` (>0))
instance QC.Arbitrary FDistribution where
  arbitrary =  fDistribution
           <$> ((abs <$> arbitrary) `suchThat` (>0))
           <*> ((abs <$> arbitrary) `suchThat` (>0))

-- Parameters for distribution testing. Some distribution require
-- relaxing parameters a bit
class Param a where
  -- Precision for quantileIsInvCDF
  invQuantilePrec :: a -> Double
  invQuantilePrec _ = 1e-14
  -- Distribution is OK for testing limits
  okForInfLimit :: a -> Bool
  okForInfLimit _ = True

instance Param a

instance Param StudentT where
  invQuantilePrec _ = 1e-13
  okForInfLimit   d = studentTndf d > 0.75

instance Param (LinearTransform StudentT) where
  invQuantilePrec _ = 1e-13
  okForInfLimit   d = (studentTndf . linTransDistr) d > 0.75

instance Param FDistribution where
  invQuantilePrec _ = 1e-12

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

unitTests :: Test
unitTests = testGroup "Unit tests"
  [ testAssertion "density (gammaDistr 150 1/150) 1 == 4.883311" $
      4.883311418525483 =~ (density (gammaDistr 150 (1/150)) 1)
    -- Student-T
  , testStudentPDF 0.3  1.34  0.0648215  -- PDF
  , testStudentPDF 1    0.42  0.27058
  , testStudentPDF 4.4  0.33  0.352994
  , testStudentCDF 0.3  3.34  0.757146   -- CDF
  , testStudentCDF 1    0.42  0.626569
  , testStudentCDF 4.4  0.33  0.621739
    -- Student-T General
  , testStudentUnstandardizedPDF 0.3    1.2  4      0.45 0.0533456  -- PDF
  , testStudentUnstandardizedPDF 4.3  (-2.4) 3.22 (-0.6) 0.0971141
  , testStudentUnstandardizedPDF 3.8    0.22 7.62   0.14 0.0490523
  , testStudentUnstandardizedCDF 0.3    1.2  4      0.45 0.458035   -- CDF
  , testStudentUnstandardizedCDF 4.3  (-2.4) 3.22 (-0.6) 0.698001
  , testStudentUnstandardizedCDF 3.8    0.22 7.62   0.14 0.496076
    -- F-distribution
  , testFdistrPDF  1  3   3     (1/(6 * pi)) -- PDF
  , testFdistrPDF  2  2   1.2   0.206612
  , testFdistrPDF  10 12  8     0.000385613179281892790166
  , testFdistrCDF  1  3   3     0.81830988618379067153 -- CDF
  , testFdistrCDF  2  2   1.2   0.545455
  , testFdistrCDF  10 12  8     0.99935509863451408041
  ]
  where
    -- Student-T
    testStudentPDF ndf x exact
      = testAssertion (printf "density (studentT %f) %f ~ %f" ndf x exact)
      $ eq 1e-5  exact  (density (studentT ndf) x)
    testStudentCDF ndf x exact
      = testAssertion (printf "cumulative (studentT %f) %f ~ %f" ndf x exact)
      $ eq 1e-5  exact  (cumulative (studentT ndf) x)
    -- Student-T General
    testStudentUnstandardizedPDF ndf mu sigma x exact
      = testAssertion (printf "density (studentTUnstandardized %f %f %f) %f ~ %f" ndf mu sigma x exact)
      $ eq 1e-5  exact  (density (studentTUnstandardized ndf mu sigma) x)
    testStudentUnstandardizedCDF ndf mu sigma x exact
      = testAssertion (printf "cumulative (studentTUnstandardized %f %f %f) %f ~ %f" ndf mu sigma x exact)
      $ eq 1e-5  exact  (cumulative (studentTUnstandardized ndf mu sigma) x)
    -- F-distribution
    testFdistrPDF n m x exact
      = testAssertion (printf "density (fDistribution %i %i) %f ~ %f [got %f]" n m x exact d)
      $ eq 1e-5  exact d
      where d = density (fDistribution n m) x
    testFdistrCDF n m x exact
      = testAssertion (printf "cumulative (fDistribution %i %i) %f ~ %f [got %f]" n m x exact d)
      $ eq 1e-5  exact d
      where d = cumulative (fDistribution n m) x

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
module FourierPinwheel.Hypergeo1F1 where

import           Data.Complex
import           Data.List              as L
import           Math.Sequence.Converge

{-# INLINE hypergeom #-}
hypergeom ::
     (Eq a, RealFloat a, Enum a)
  => Complex a
  -> Complex a
  -> Complex a
  -> Complex a
hypergeom a b z =
  converge .
  scanl (+) 0 .
  scanl (*) 1 .
  L.map (\i -> (a + ((i - 1) :+ 0)) / (b + ((i - 1) :+ 0)) * z / (i :+ 0)) $
  [1 ..]

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE DeriveAnyClass, DeriveGeneric, TemplateHaskell, OverloadedStrings #-}
module Numeric.Sundials.Types
  ( OdeRhsCType
  , OdeRhs(..)
  , UserData
  , Jacobian
  , StepControl(..)
  , ODEOpts(..)
  , SundialsDiagnostics(..)
  , emptyDiagnostics
  , SundialsSolution(..)
  , EventInfo(..)
  , CrossingDirection(..)
  , EventSpec(..)
  , SunVector(..)
  , SunIndexType
  , SunRealType
  , sunContentLengthOffset
  , sunContentDataOffset
  , sunCtx
  )
  where

import           Data.Int (Int32)
import qualified Data.Vector.Storable as VS
import qualified Data.Map.Strict as Map
import qualified Language.Haskell.TH as TH

import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix)
import           Control.DeepSeq (NFData)
import           GHC.Generics (Generic)
import           Foreign.C.Types
import           Foreign.Ptr
import           Language.C.Types as CT
import           Language.C.Inline.Context
import           Numeric.Sundials.Arkode (SunVector(..), SunMatrix(..),
                                          SunIndexType, SunRealType,
                                          sunContentLengthOffset,
                                          sunContentDataOffset)

-- | The type of the C ODE RHS function.
type OdeRhsCType = CDouble -> Ptr SunVector -> Ptr SunVector -> Ptr UserData -> IO CInt

data UserData

-- | The right-hand side of the ODE system.
--
-- Can be either a Haskell function or a pointer to a C function.
data OdeRhs
  = OdeRhsHaskell (CDouble -> VS.Vector CDouble -> VS.Vector CDouble)
  | OdeRhsC (FunPtr OdeRhsCType) (Ptr UserData)

type Jacobian = Double -> Vector Double -> Matrix Double

-- | Adaptive step-size control
-- functions.
--
-- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control)
-- allows the user to control the step size adjustment using
-- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where
-- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\)
-- is the required relative error, \(s_i\) is a vector of scaling
-- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and
-- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\).
--
-- [CVode](https://computation.llnl.gov/projects/sundials/cvode)
-- and [ARKode](https://computation.llnl.gov/projects/sundials/arkode)
-- allow the user to control the step size adjustment using
-- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with
-- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl),
-- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no
-- effect.
data StepControl = X     Double Double -- ^ absolute and relative tolerance for \(y\); in GSL terms, \(a_{y} = 1\) and \(a_{dy/dt} = 0\); in ARKode terms, the \(\eta^{abs}_i\) are identical
                 | X'    Double Double -- ^ absolute and relative tolerance for \(\dot{y}\); in GSL terms, \(a_{y} = 0\) and \(a_{dy/dt} = 1\); in ARKode terms, the latter is treated as the relative tolerance for \(y\) so this is the same as specifying 'X' which may be entirely incorrect for the given problem
                 | XX'   Double Double Double Double -- ^ include both via relative tolerance
                                                     -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
                 | ScXX' Double Double Double Double (Vector Double) -- ^ scale absolute tolerance of \(y_i\); in ARKode terms, \(a_{{dy}/{dt}}\) is ignored, \(\eta^{abs}_i = s_i \epsilon^{abs}\) and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
  deriving (Eq, Ord, Show, Read)

data ODEOpts method = ODEOpts {
    maxNumSteps :: Int32
  , minStep     :: Double
  , maxFail     :: Int32
  , odeMethod   :: method
  , stepControl :: StepControl
  , initStep    :: Maybe Double
    -- ^ initial step size - by default, CVode
    -- estimates the initial step size to be the
    -- solution \(h\) of the equation
    -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
    -- \(\ddot{y}\) is an estimated value of the second
    -- derivative of the solution at \(t_0\)
  } deriving (Read, Show, Eq, Ord)

data SundialsDiagnostics = SundialsDiagnostics {
    odeGetNumSteps               :: Int
  , odeGetNumStepAttempts        :: Int
  , odeGetNumRhsEvals_fe         :: Int
  , odeGetNumRhsEvals_fi         :: Int
  , odeGetNumLinSolvSetups       :: Int
  , odeGetNumErrTestFails        :: Int
  , odeGetNumNonlinSolvIters     :: Int
  , odeGetNumNonlinSolvConvFails :: Int
  , dlsGetNumJacEvals            :: Int
  , dlsGetNumRhsEvals            :: Int
  } deriving Show

emptyDiagnostics :: SundialsDiagnostics
emptyDiagnostics = SundialsDiagnostics 0 0 0 0 0 0 0 0 0 0

data SundialsSolution =
  SundialsSolution
  { actualTimeGrid :: VS.Vector Double    -- ^ actual time grid returned by the solver (with duplicated event times)
  , solutionMatrix :: Matrix Double       -- ^ matrix of solutions: each column is an unknwown
  , eventInfo      :: [EventInfo]         -- ^ event infos, as many items as triggered events during the simulation
  , diagnostics    :: SundialsDiagnostics -- ^ usual Sundials diagnostics
  }

data EventInfo =
  EventInfo
  { eventTime     :: !Double            -- ^ time at which event was triggered
  , eventIndex    :: !Int               -- ^ which index was triggered
  , rootDirection :: !CrossingDirection -- ^ in which direction ((+)->(-) or (-)->(+)) the root is crossed
  }
  deriving (Generic, Show, NFData)

-- | The direction in which a function should cross the x axis
data CrossingDirection = Upwards | Downwards | AnyDirection
  deriving (Generic, Eq, Show, NFData)

data EventSpec = EventSpec
  { eventCondition :: Double -> VS.Vector Double -> Double
  , eventDirection :: CrossingDirection
  , eventUpdate :: Double -> VS.Vector Double -> VS.Vector Double
  }

sunTypesTable :: Map.Map TypeSpecifier TH.TypeQ
sunTypesTable = Map.fromList
  [
    (TypeName "sunindextype", [t| SunIndexType |] )
  , (TypeName "SunVector",    [t| SunVector |] )
  , (TypeName "SunMatrix",    [t| SunMatrix |] )
  , (TypeName "UserData",     [t| UserData |] )
  ]

-- | Allows to map between Haskell and C types
sunCtx :: Context
sunCtx = mempty {ctxTypesTable = sunTypesTable}

max_stars_repo_name: idontgetoutmuch/hmatrix-sundials

max_issues_repo_name: idontgetoutmuch/hmatrix-sundials

max_forks_repo_name: idontgetoutmuch/hmatrix-sundials