text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE RecordWildCards #-}
module Bench where
import           Graphics.Rendering.Chart.Easy
import           Graphics.Rendering.Chart.Backend.Diagrams
import           System.FilePath
import           System.Directory
import qualified Data.Vector                   as Vec
import qualified Statistics.Sample             as Stat
import qualified Data.Csv                      as Csv
import qualified Data.ByteString.Lazy          as B
import qualified Data.Map                      as Map
import           Data.Type.Natural
import           Data.List                      ( intercalate )
import qualified Data.Store                    as Store
                                         hiding ( size )
import qualified Data.ByteString               as Bs
import           Options.Applicative
import           Data.Semigroup                 ( (<>) )

import           CodeGen.Type
import           ParPattern
import           CodeGen
import           Example.Mergesort
import           Example.DotProd
import           Example.Wordcount
import           CodeGen.Data

class RandomGenData a where
    genData :: Int -> Int -> IO a

class GenBenchData a where
    mkBenchData :: Int -> BenchRawData a

instance (GenBenchData a, GenBenchData b) => GenBenchData (a, b) where
    mkBenchData len = RPair (mkBenchData len) (mkBenchData len)

instance GenBenchData [a] where
    mkBenchData len = RList len

instance RandomGenData a => RandomGenData (a, a) where
    genData x y = do
        first  <- genData x y
        second <- genData x y
        return (first, second)

instance RandomGenData [Int] where
    genData x _y = return [x]

type Benchable a = (Serialise a, GenBenchData a)

benchmarkCompile
    :: Benchable a
    => FilePath
    -> [Int]
    -> [Int]
    -> [Int]
    -> (Int -> ArrowPipe a b)
    -> IO ()
benchmarkCompile path seeds unrolls sizes aexpr = mapM_ benchmarkUnRolls sizes
  where
    benchmarkUnRolls size = mapM_ (benchmarkSeeds size) unrolls
    benchmarkSeeds size unroll = mapM_ (benchmarkList size unroll) seeds
    benchmarkList size unroll seed = do
        let x = mkBenchData (2 ^ size)
        codeGenCompileWithHeader
            x
            (runPipe1 zero (aexpr unroll))
            (path </> (intercalate "_" [show size, show unroll, show seed]))

benchmarkRun :: FilePath -> Int -> [Int] -> [Int] -> [Int] -> IO ()
benchmarkRun path time seeds unrolls sizes =
    rawData >>= (Bs.writeFile (path </> rawDataName) . Store.encode)
  where
    rawData = fmap (zip sizes) $ mapM benchmarkUnRolls sizes
    benchmarkUnRolls size =
        fmap (zip unrolls) $ mapM (benchmarkSeeds size) unrolls
    benchmarkSeeds size unroll =
        toBenchData <$> (concat <$> mapM (benchmarkList size unroll) seeds)
    benchmarkList size unroll seed = mapM
        (const $ codeGenBenchRun
            (path </> (intercalate "_" [show size, show unroll, show seed]))
        )
        [1 .. time]

benchmarkRun2 :: FilePath -> Int -> [Int] -> [Int] -> [Int] -> IO ()
benchmarkRun2 path time seeds unrolls sizes =
    (toRawData <$> rawData)
        >>= (B.writeFile (path </> rawDataName2) . Csv.encode)
  where
    toRawData xs =
        [ RawData size' 0 unroll' (fst val)
        | (size'  , rest) <- xs
        , (unroll', val ) <- rest
        ]

    rawData = fmap (zip sizes) $ mapM benchmarkUnRolls sizes
    benchmarkUnRolls size =
        fmap (zip unrolls) $ mapM (benchmarkSeeds size) unrolls
    benchmarkSeeds size unroll =
        toBenchData <$> (concat <$> mapM (benchmarkList size unroll) seeds)
    benchmarkList size unroll seed = mapM
        (const $ codeGenBenchRun
            (path </> (intercalate "_" [show size, show unroll, show seed]))
        )
        [1 .. time]

rawDataName :: String
rawDataName = "rawData"

rawDataName2 :: String
rawDataName2 = "rawData.csv"

benchmarkCollectData :: FilePath -> IO ()
benchmarkCollectData path = do
    rawDataRead <- Bs.readFile (path </> rawDataName)
    let rawData = case Store.decode rawDataRead of
            Left  x -> error $ show x
            Right y -> y
    removeFile (path </> rawDataName)
    writeBenchmarkCSV path rawData
    writeBenchmarkPlot path rawData
    writeBenchmarkSpeedUpPlot path rawData

benchmarkCollectData2 :: FilePath -> IO ()
benchmarkCollectData2 path = do
    rawData <- rawDataDecode path
    -- removeFile (path </> rawDataName2)
    -- writeBenchmarkCSV path rawData
    writeBenchmarkPlot path rawData
    writeBenchmarkSpeedUpPlot path rawData

benchmarking
    :: Benchable a
    => FilePath
    -> [Int]
    -> [Int]
    -> [Int]
    -> (Int -> ArrowPipe a b)
    -> IO ()
benchmarking path seeds unrolls sizes expr = do
    benchmarkCompile path seeds unrolls sizes expr
    benchmarkRun path 1 seeds unrolls sizes
    benchmarkCollectData path

mean :: [Double] -> Double
mean = (Stat.mean . Vec.fromList)

stdDev :: [Double] -> Double
stdDev = (Stat.stdDev . Vec.fromList)

toBenchData :: [Double] -> (Double, Double)
toBenchData xs = (mean xs, stdDev xs)

writeBenchmarkPlot :: FilePath -> [(Int, [(Int, (Double, Double))])] -> IO ()
writeBenchmarkPlot path sourceData = toFile def imagePath $ do
    layout_title .= "benchmark"
    (plotLines . toPlottableData) sourceData
    where imagePath = path </> "bench.svg"

-- the first int is size
writeBenchmarkSpeedUpPlot
    :: FilePath -> [(Int, [(Int, (Double, Double))])] -> IO ()
writeBenchmarkSpeedUpPlot path sourceData = toFile def imagePath $ do
    layout_title .= "Speedup"
    (plotLines . toPlottableData . toSpeedUpData) sourceData
    where imagePath = path </> "benchSpeedup.svg"

-- the first int is k
plotLines :: [(Int, [(Int, Double)])] -> EC (Layout Int Double) ()
plotLines = mapM_ (uncurry plotLine)
  where
    plotLine :: Int -> [(Int, Double)] -> EC (Layout Int Double) ()
    plotLine title xs = plot (line ("k" ++ show title) [xs])

toSpeedUpData
    :: [(Int, [(Int, (Double, Double))])] -> [(Int, [(Int, (Double, Double))])]
toSpeedUpData = fmap (\(x, y) -> (x, helper helper2 y))
  where
    helper2
        :: (Int, (Double, Double))
        -> (Int, (Double, Double))
        -> (Int, (Double, Double))
    helper2 (size, parTime) (_, serTime) =
        (size, (fst serTime / fst parTime, snd parTime))

    helper f xs = fmap (\x -> f x (head xs)) (tail xs)

toPlottableData
    :: [(Int, [(Int, (Double, Double))])] -> [(Int, [(Int, Double)])]
toPlottableData = helper . flattenMap
  where
    helper :: (Ord a, Ord b) => [(a, b, Double)] -> [(b, [(a, Double)])]
    helper source = Map.toList $ foldl
        (\map (key1, key2, val) -> Map.insertWith (++) key2 [(key1, val)] map)
        (Map.empty)
        source

    flattenMap :: [(Int, [(Int, (Double, Double))])] -> [(Int, Int, Double)]
    flattenMap source =
        [ (key1, key2, d1) | (key1, val) <- source, (key2, (d1, _d2)) <- val ]

writeBenchmarkCSV :: FilePath -> [(Int, [(Int, (Double, Double))])] -> IO ()
writeBenchmarkCSV path sourceData = B.writeFile csvPath csvRaw
  where
    records :: [BenchData]
    records = concatMap (\(x, y) -> toRecordHelper x y) sourceData

    toRecordHelper :: Int -> [(Int, (Double, Double))] -> [BenchData]
    toRecordHelper x xs = fmap
        (\(unroll, (meanVal, stdVal)) -> BenchData (x, unroll, meanVal, stdVal))
        xs

    csvRaw  = Csv.encodeByName csvTitle records

    csvPath = path </> "bench.csv"

newtype BenchData = BenchData (Int, Int, Double, Double)

instance Csv.ToNamedRecord BenchData where
    toNamedRecord (BenchData (size, unroll, meanVal, stdVal)) = Csv.namedRecord
        [ "size" Csv..= size
        , "unroll" Csv..= unroll
        , "mean" Csv..= meanVal
        , "std" Csv..= stdVal
        ]

data RawData = RawData { size :: Int, seed :: Int, unroll :: Int, time :: Double }

instance Csv.FromRecord RawData where
    parseRecord v =
        RawData <$> v Csv..! 0 <*> v Csv..! 1 <*> v Csv..! 2 <*> v Csv..! 3

instance Csv.ToRecord RawData where
    toRecord RawData {..} = Csv.record
        [ Csv.toField size
        , Csv.toField seed
        , Csv.toField unroll
        , Csv.toField time
        ]

rawDataDecode :: FilePath -> IO [(Int, [(Int, (Double, Double))])]
rawDataDecode path = do
    rawData <- B.readFile filePath
    let decoder =
            (Csv.decode Csv.NoHeader rawData) :: Either
                    String
                    (Vec.Vector RawData)
    return $ rawDatadecode' decoder
  where
    filePath = path </> rawDataName2
    rawDatadecode'
        :: Either String (Vec.Vector RawData)
        -> [(Int, [(Int, (Double, Double))])]
    rawDatadecode' x = case x of
        Left  msg        -> error $ "failed " ++ msg
        Right sourceData -> (toList . toMap'' . toMap' . toMap) sourceData

    toMap = Vec.foldl
        (\map record ->
            Map.insertWith (++) (size record, unroll record) [time record] map
        )
        Map.empty

    toMap'' = Map.foldrWithKey
        (\(mySize, myUnroll) val acc ->
            Map.insertWith (Map.union) mySize (Map.singleton myUnroll val) acc
        )
        Map.empty

    toMap' = Map.map toBenchData

    toList = (Map.toList . Map.map Map.toList)

csvTitle :: Vec.Vector Csv.Name
csvTitle = Vec.fromList ["size", "unroll", "mean", "std"]

runParser :: Parser (Bool, Bool)
runParser =
    (,)
        <$> switch (long "collect" <> short 'c' <> help "whether to collect")
        <*> switch (long "run" <> short 'r' <> help "whether to run or codegen")

genBenchmarkFile
    :: Benchable a
    => FilePath
    -> [Int]
    -> [Int]
    -> (Int -> ArrowPipe a b)
    -> IO ()
genBenchmarkFile path unrolls sizes expr =
    benchmarkEntry path 1 [1] unrolls sizes expr

benchmarkEntry
    :: Benchable a
    => FilePath
    -> Int
    -> [Int]
    -> [Int]
    -> [Int]
    -> (Int -> ArrowPipe a b)
    -> IO ()
benchmarkEntry path time seeds unrolls sizes expr = execParser opts >>= \x ->
    case x of
        (True , _    ) -> benchmarkCollectData2 path
        (False, True ) -> benchmarkRun2 path time seeds unrolls sizes
        (False, False) -> benchmarkCompile path seeds unrolls sizes expr
    where opts = info (runParser <**> helper) mempty

benchmarkNewEntry :: FilePath -> Int -> [Int] -> [Int] -> [Int] -> IO ()
benchmarkNewEntry mainPath time seeds unrolls sizes =
    execParser opts >>= \x -> case x of
        (True, _) ->
            mapM_ (\path -> benchmarkCollectData (mainPath </> path)) names
        (False, True) -> mapM_
            (\path -> benchmarkRun (mainPath </> path) time seeds unrolls sizes)
            names
        (False, False) ->
            mapM_ (\name -> myHelper name (mainPath </> name)) names
  where
    names = ["intcount", "mergesort", "dotprod"]
    opts  = info (runParser <**> helper) mempty

    myHelper "intcount" path =
        benchmarkCompile path seeds unrolls sizes wordCount
    myHelper "mergesort" path =
        benchmarkCompile path seeds unrolls sizes mergeSort
    myHelper "dotprod" path = benchmarkCompile path seeds unrolls sizes dotProd
    myHelper _         _    = error "benchmark example not implemented"

max_stars_repo_name: Richard-zhang/SPar

max_issues_repo_name: Richard-zhang/SPar

max_forks_repo_name: Richard-zhang/SPar


text:
module HLearn.Optimization.Conic
    where

import Control.DeepSeq
import Control.Monad
import Control.Monad.Random
import Control.Monad.ST
import Data.List
import Data.List.Extras
import Debug.Trace
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as VM
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import qualified Data.Vector.Generic as VG
import qualified Data.Vector.Generic.Mutable as VGM
import qualified Data.Vector.Algorithms.Intro as Intro
import Numeric.LinearAlgebra hiding ((<>))
import qualified Numeric.LinearAlgebra as LA
import Data.Random.Normal

import HLearn.Algebra
import qualified HLearn.Optimization.Common as Recipe
import qualified HLearn.Optimization.LineMinimization as LineMin

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

conicprojection :: Matrix Double -> Matrix Double
conicprojection m = cmap realPart $ u LA.<> lambda' LA.<> trans u
    where
        lambda' = cmap (\x -> if realPart x < 0 then 0 else x) lambda
        lambda = diagRect 0 l (VG.length l) (VG.length l)
        (l,u) = eig m

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

data RandomConicPersuit a = RandomConicPersuit
    { _stdgen :: !StdGen
    , _soln :: !a
    , _fx :: !(Scalar a)
    , _solnlast :: !a
    }

-- data OptInfo a = OptInfo
--     { _stdgen :: !StdGen
--     , _x :: !a
--     , _fx ::
--     , _xold :: !a
--     }

itr :: Int -> (tmp -> a) -> (tmp -> Bool) -> (tmp -> tmp) -> tmp -> a
itr i result stop step init = --trace ("i="++show i++"; fx="++show (_fx init)) $ --trace ("i="++show i++"; init="++show init) $
  if i==0 || stop init
    then result init
    else itr (i-1) result stop step (step init)

randomConicPersuit f x0 = _soln $ itr 10 id (\x -> False) (step_RandomConicPersuit f) init
    where
        init = RandomConicPersuit
            { _stdgen = mkStdGen $ round $ sumElements x0
            , _soln = x0
            , _fx = f x0
            , _solnlast = x0
            }

conicProjection f x0 = _soln $ argmin _fx [itr 1 id (\x -> False) (step_RandomConicPersuit f) (init i) | i <- [0..10]]
    where
        init i = RandomConicPersuit
            { _stdgen = mkStdGen $ i+(round $ sumElements x0)
            , _soln = x0
            , _fx = f x0
            , _solnlast = x0
            }

stop_RandomConicPersuit (RandomConicPersuit stdgen soln fx solnlast) = undefined

step_RandomConicPersuit f (RandomConicPersuit stdgen soln fx solnlast) = --trace ("lambda = "++show lambda++ "; phi="++show phi) $
    ret
    where
--         (x', stdgen') = runRand ((LA.fromList . take (rows soln)) `liftM` getRandomRs (-1,1)) stdgen
        normL 0 (xs,g) = (xs,g)
        normL i (xs,g) = normL (i-1) (x:xs,g')
            where
                (x,g') = normal g

        (x'std, stdgen') = normL (rows soln) ([],stdgen)
        x' = VS.fromList x'std
--         (lambda,phi::Matrix Double) = eigSH soln
--         q = (diag $ cmap (sqrt.abs) lambda) LA.<> phi
--         x' = (scale (1-kappa) q - scale kappa (ident (rows q))) LA.<> LA.fromList x'std
--         kappa = 1e-4

        y' = asColumn x' LA.<> asRow x'
--         y' = asColumn (LA.fromList x') LA.<> asRow (LA.fromList x')

        g_alpha alpha = f $ scale alpha y' + soln
        alpha_hat = error "step_randomConic" -- LineMin._x $ runOptimization $ LineMin.brent g_alpha (LineMin.lineBracket g_alpha 0 1)

        alpha_hat_y' = scale alpha_hat y'
        g_beta beta = f $ alpha_hat_y' + scale beta soln
        beta_hat = error "step_randomConic" -- LineMin._x $ LineMin.brent g_beta (LineMin.lineBracket g_beta 0 1)

        soln' = alpha_hat_y' + scale beta_hat soln

        ret = RandomConicPersuit
            { _stdgen = stdgen'
            , _soln = soln'
            , _fx = f soln'
            , _solnlast = soln
            }

ext: hs

max_stars_repo_name: Heather/HLearn

max_issues_repo_name: Heather/HLearn

max_forks_repo_name: Heather/HLearn


text:
module Statistics.GLM
    ( glm
    , logisticReg
    ) where

import Numeric.LinearAlgebra hiding ((<>))

data RegFamily = RegFamily { g :: Matrix Double -> Vector Double -> Vector Double
                           , gprime :: Matrix Double -> Vector Double -> Vector Double
                           , variance :: Matrix Double -> Vector Double -> Vector Double
                           }

logisticReg :: RegFamily
logisticReg = RegFamily g gprime variance
  where
    g a x = 1 / (exp(-(a #> x)) + 1)
    gprime a x = exp(-(a #> x)) / (exp(-(a #> x)) + 1)**2
    variance a x = g a x * (1 - g a x)

doIteration :: Matrix Double -> Vector Double -> Vector Double -> RegFamily -> Vector R
doIteration a x b (RegFamily g gprime variance) =
  let gp = gprime a x
      w = gp ** 2 / variance a x
      z = a #> x + (b - g a x) / gp
      eqleft = tr a <> (matrix 1 (toList w) * a)
      eqright = tr a #> (w * z)
  in eqleft <\> eqright

glm :: RegFamily -> Matrix Double -> Vector Double -> Vector Double -> [Vector Double]
glm rf a x b = x : glm rf a (doIteration a x b rf) b


text:
{-# LANGUAGE TemplateHaskell #-}

module FractalsTest
( runTests
) where

import Test.QuickCheck
import Data.Complex

import Fractals
import Palette

prop_toGrayscaleLength :: [Double] -> Bool
prop_toGrayscaleLength ws = length (toBitmap grayscale ws) == 4 * (length ws)

inMandelbrotSet :: Int -> Complex Double -> Bool
inMandelbrotSet nMax z0 = nMax == (fst $ order (mandelbrot z0) 2.0 nMax z0)

prop_orderMandelbrotOutside :: Positive Int -> Complex Double -> Property
prop_orderMandelbrotOutside (Positive nMax) z0 =
    (mag2 z0 > 4.0) ==> not $ inMandelbrotSet nMax z0

toUnitCircle :: RealFloat a => Complex a -> Complex a -> Complex a
toUnitCircle w z = if magW>0.0 && magZ>0.0 then mu else (0.0 :+ 0.0)
    where magW = magnitude w
          magZ = magnitude z
          mu = if magW < magZ then w/(magZ:+0.0) else z/(magW:+0.0)

prop_orderMandelbrotCardioid :: Positive Int    -- ^ the max number of iterations
                             -> Complex Double  -- ^ required to generate a value in the unit disk
                             -> Complex Double  -- ^ required to generate a value in the unit disk
                             -> Bool
prop_orderMandelbrotCardioid (Positive nMax) w z = inMandelbrotSet nMax z0
    where z0 = mandelbrotCardioid $ toUnitCircle w z

prop_orderMandelbrotCircle :: Positive Int    -- ^ the max number of iterations
                           -> Complex Double  -- ^ required to generate a value in the unit disk
                           -> Complex Double  -- ^ required to generate a value in the unit disk
                           -> Bool
prop_orderMandelbrotCircle (Positive nMax) w z = inMandelbrotSet nMax z0
    where z0 = (-1.0) + mu * 0.25
          mu = toUnitCircle w z

return []
runTests :: IO Bool
runTests = $quickCheckAll

max_stars_repo_name: arekfu/fractals-haskell

max_issues_repo_name: arekfu/fractals-haskell

max_forks_repo_name: arekfu/fractals-haskell


text:
module MeasurementTests
  ( testMeasurements
  ) where

import Data.Maybe (fromJust)
import qualified Data.Vector as V
import HABQTlib.Data
import HABQTlib.MeasurementProcessing
import Numeric.LinearAlgebra (Complex(..), (><))
import qualified Numeric.LinearAlgebra as LA
import StateGenTests (arbPropPure)
import qualified Test.QuickCheck as QC
import Test.QuickCheck (Property, (.&&.))
import TestHelpers

povmValidProp :: PurePOVM -> Property
povmValidProp wsts = normedContents .&&. sumsToUnity
  where
    dim :: Int
    dim =
      (\(WeighedPureStateVector (_, PureStateVector sv)) -> LA.rows sv) $
      V.head wsts
    normedContents =
      QC.conjoin . V.toList $
      V.map (\(WeighedPureStateVector (_, sv)) -> arbPropPure sv) wsts
    msum :: LA.Matrix (LA.Complex Double)
    msum =
      V.foldl'
        (\acc (WeighedPureStateVector (w, PureStateVector sv)) ->
           acc + LA.scalar (w :+ 0) * (sv LA.<> LA.tr sv))
        ((dim >< dim) (repeat (0 :+ 0)))
        wsts
    sumsToUnity = LA.norm_Frob (msum - LA.ident dim) <= 1e-6

povmProbabilityNorm :: PurePOVM -> DensityMatrix -> Property
povmProbabilityNorm povm dm =
  QC.property $ abs (1 - V.sum (measurementProbs povm dm)) <= 1e-8

blochHelper :: Double -> Double -> SingleQbParam
blochHelper = (,)

coordProp :: Property
coordProp =
  let prop sv = fidelity sv sv' > 0.99999
        where
          sv' = blochAnglesToSV . fromJust . svToAngles $ sv
   in QC.forAll (QC.resize 2 QC.arbitrary) prop

testBlochCoords :: IO ()
testBlochCoords = do
  putStrLn "Testing Bloch coordinate transofrmations:"
  QC.quickCheckWith QC.stdArgs {QC.maxSuccess = 1000} coordProp

testAntipodal1QbPOVM :: IO ()
testAntipodal1QbPOVM = do
  putStrLn "Testing antipodal single qubit POVM:"
  QC.quickCheckWith QC.stdArgs {QC.maxSuccess = 1000} $
    ((povmValidProp . mkAntipodalPOVM) .) . blochHelper

genProdPOVM :: Int -> QC.Gen PurePOVM
genProdPOVM qbnum = do
  paramVs' <- QC.vectorOf qbnum QC.arbitrary
  let paramVs = fmap (uncurry blochHelper) paramVs'
      svs = map mkAntipodalPOVM paramVs
  return $ productPOVM svs

testProductPOVM :: IO ()
testProductPOVM = do
  putStrLn "Testing productPOVMs"
  let gen = do
        qbnum <- QC.suchThat QC.arbitrary (\n -> n >= 1 && n <= 6)
        let dim :: Dim
            dim = 2 ^ qbnum
        rank <- QC.suchThat QC.arbitrary (\n -> n <= dim && n > 0)
        povm <- genProdPOVM qbnum
        dm <- arbDM dim rank
        return (povm, dm)
  QC.quickCheckWith QC.stdArgs {QC.maxSuccess = 10000} $
    QC.forAll
      gen
      (\(povm, dm) -> povmValidProp povm .&&. povmProbabilityNorm povm dm)

testMeasurements :: IO ()
testMeasurements = do
  testBlochCoords
  testAntipodal1QbPOVM
  testProductPOVM

max_stars_repo_name: Belinsky-L-V/HABQT

max_issues_repo_name: Belinsky-L-V/HABQT

max_forks_repo_name: Belinsky-L-V/HABQT


text:
module Numeric.Fourier
    ( fft
    , ifft

    -- * Testing
    , prop_fourier
    )
  where

import Data.Complex ( cis, realPart, Complex )
import Data.Vector.Unboxed ( Unbox, Vector, (!) )

import qualified Data.Vector.Unboxed as V

{-----------------------------------------------------------------------------
    Fast Fourier transformation
------------------------------------------------------------------------------}
-- | Compute the discrete Fourier transform (DFT) of a vector
-- using a fast Fourier transform (FFT) algorithm.
--
-- The discrete Fourier transform is defined as
--
-- > (fft xs) ! k = sum [ exp(- j*k * 2*pi*i/n) * (xs ! j) | j <- [0..(n-1)] ]
-- >   where n = length xs
--
-- Let $n$ be the length of the input vector.
-- The algorithm recursively divides this vector into parts of equal lengths.
-- Consequently, the running time is $O(n \log n)$ if $n$ is a power of two.
-- However, the running time will degrade to $O(n^2)$
-- if $n$ contains prime factors other than $2$.
fft :: Vector (Complex Double) -> Vector (Complex Double)
fft xs
    | n <= 1    = xs
    | even n    = fftStep xs
    | otherwise = dftSlow xs
  where
    n = V.length xs

-- | Perform an fft step on a vector of even, nonzero length.
--
-- TODO: Perform more operations in-place.
fftStep :: Vector (Complex Double) -> Vector (Complex Double)
fftStep xs =
        V.izipWith (\k ev od -> ev + xi k * od) dftL dftR
    <>  V.izipWith (\k ev od -> ev - xi k * od) dftL dftR
  where
    n  = V.length xs
    n2 = n `div` 2

    -- exp(- 2πi·k (2j+1)/n) = exp(- 2πi k/n) * exp(- 2πi·kj/(2n))
    xi k = cis (- 2*pi* fromIntegral k/fromIntegral n)

    xsEven = V.generate n2 $ \j -> xs ! (2*j)
    xsOdd  = V.generate n2 $ \j -> xs ! (2*j+1)

    dftL = fft xsEven
    dftR = fft xsOdd

-- | Inverse of the discrete fourier transform.
--
-- Implemented using 'fft'.
ifft :: Vector (Complex Double) -> Vector (Complex Double)
ifft xs = V.generate n $ \k ->
    (if k == 0 then ys ! 0 else ys ! (n-k)) / fromIntegral n
  where
    ys = fft xs
    n = V.length xs

{-----------------------------------------------------------------------------
    Testing
------------------------------------------------------------------------------}
-- | Definition of the discrete Fourier transform (DFT).
-- Slow, for testing only.
dftSlow :: Vector (Complex Double) -> Vector (Complex Double)
dftSlow xs = V.generate n $ \k ->
    V.sum $ V.imap (\j x -> cis (- 2*pi * double (k*j) / double n) * x) xs
  where
    double = fromIntegral
    n = V.length xs

-- | Test whether the 'fft' function really is a Fourier transform,
-- i.e. gives the same result as 'dftSlow'.
prop_fourier :: Vector (Complex Double) -> Bool
prop_fourier xs = fft xs `approx` dftSlow xs
  where
    dist xs ys = realPart $ V.sum $ V.zipWith (\x y -> abs (x-y)) xs ys
    approx xs ys = dist xs ys / dist xs (V.map (const 0) ys) < eps
    eps = 1e-15

max_stars_repo_name: HeinrichApfelmus/chebyshev

max_issues_repo_name: HeinrichApfelmus/chebyshev

max_forks_repo_name: HeinrichApfelmus/chebyshev


text:
module Main where

import Numeric.LinearAlgebra

data Model = Model { w::Matrix R, beta::Vector R} deriving Show

-- sigmoid
sigmoid :: R -> R
sigmoid x = 1 / (1 + exp (-x))

elm :: Matrix R -> Vector R -> Int -> IO Model
elm x y n = do _w <- randn n $ cols x
               let hidden = cmap sigmoid $ _w <> (tr x)
                   _beta = tr (pinv hidden) #> y
               return Model { w = _w, beta = _beta}

classify :: Model -> Matrix R -> Vector R
classify m x =
  let _w = w m
      _beta = beta m
  in cmap signum $ tr (cmap sigmoid $ _w <> (tr x)) #> _beta

main :: IO ()
main = do
  model <- elm x y n
  print "class label"
  print y
  print $ "result: hidden_neuron = " ++ show n
  print $ classify model x
    where x = matrix 3 [0, 0, 1,
                        0, 1, 1,
                        1, 0, 1,
                        1, 1, 1]
          y = vector [-1, 1, 1, -1]
          n = 10

max_stars_repo_name: masaponto/haskell-elm

max_issues_repo_name: masaponto/haskell-elm

max_forks_repo_name: masaponto/haskell-elm


text:
{-# LANGUAGE ApplicativeDo                            #-}
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE EmptyCase                                #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE KindSignatures                           #-}
{-# LANGUAGE LambdaCase                               #-}
{-# LANGUAGE PartialTypeSignatures                    #-}
{-# LANGUAGE QuantifiedConstraints                    #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TupleSections                            #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE ViewPatterns                             #-}
{-# OPTIONS_GHC -fno-warn-partial-type-signatures     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Extra.Solver    #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

import           Backprop.Learn
import           Control.Applicative
import           Control.DeepSeq
import           Control.Exception
import           Control.Monad
import           Control.Monad.IO.Class
import           Control.Monad.Primitive
import           Control.Monad.Trans.Class
import           Control.Monad.Trans.State
import           Data.Char
import           Data.Conduit
import           Data.Default
import           Data.Kind
import           Data.Proxy
import           Data.Singletons
import           Data.Singletons.TypeLits
import           Data.Time
import           Data.Type.Equality
import           Data.Type.Tuple
import           GHC.TypeLits.Compare
import           GHC.TypeNats
import           Numeric.LinearAlgebra.Static.Backprop hiding ((<>))
import           Numeric.Natural
import           Numeric.Opto
import           Options.Applicative
import           Statistics.Distribution
import           Statistics.Distribution.Normal
import           Text.Printf
import qualified Data.Conduit.Combinators                     as C
import qualified Data.Vector.Sized                            as SV
import qualified System.Random.MWC                            as MWC

data Mode = CSimulate (Maybe Int)
          | forall s p. ( Floating s
                        , Floating p
                        , NFData s
                        , NFData p
                        , Metric Double s
                        , Metric Double p
                        , Show s
                        , Show p
                        , Initialize p
                        , Initialize s
                        , Backprop p
                        , Backprop s
                        , forall m. PrimMonad m => LinearInPlace m Double p
                        , forall m. PrimMonad m => LinearInPlace m Double s
                        )
                        => CLearn (ModelPick p s)

data ModelPick :: Type -> Type -> Type where
    MARIMA :: Sing p -> Sing d -> Sing q
           -> ModelPick (ARIMAp p q) (ARIMAs p d q)
    MFCRNN :: Sing h
           -> ModelPick (LRp h 1 :# LRp (1 + h) h) (R h)
    MLSTM  :: Sing h
           -> ModelPick (LRp h 1 :# LSTMp 1 h)     (R h :# R h)
    MGRU   :: Sing h
           -> ModelPick (LRp h 1 :# GRUp 1 h)      (R h)

data Process = PSin

modelLearn :: ModelPick p s -> Model ('Just p) ('Just s) Double Double
modelLearn = \case
    MARIMA (SNat :: Sing pp) (SNat :: Sing d) (SNat :: Sing q) ->
              arima @pp @d @q
    MFCRNN (SNat :: Sing h) ->
              funcD sumElements
           <~ fc
           <~ fcra logistic logistic
           <~ funcD konst
    MLSTM  (SNat :: Sing h) ->
              funcD sumElements
           <~ fc
           <~ lstm
           <~ funcD konst
    MGRU   (SNat :: Sing h) ->
              funcD sumElements
           <~ fc
           <~ gru
           <~ funcD konst

data Options = O { oMode     :: Mode
                 , oProcess  :: Process
                 , oNoise    :: Double
                 , oInterval :: Int
                 , oLookback :: Natural
                 }

processConduit
    :: PrimMonad m
    => Process
    -> MWC.Gen (PrimState m)
    -> ConduitT i Double m ()
processConduit = \case
    PSin -> \g -> void . flip (foldr (>=>) pure) (0, 1) . repeat $ \(t, v) -> do
      dv <- genContVar (normalDistr 0 0.025) g
      let v' = min 2 . max 0.5 $ v + dv
          t' = t + v'
      yield (sin (2 * pi * (1/25) * t))
      return (t', v')

noisyConduit
    :: PrimMonad m
    => Double
    -> MWC.Gen (PrimState m)
    -> ConduitT Double Double m ()
noisyConduit σ g = C.mapM $ \x ->
    (x + ) <$> genContVar (normalDistr 0 σ) g

-- type Context = ConduitT (Vector)

main :: IO ()
main = MWC.withSystemRandom @IO $ \g -> do
    O{..} <- execParser $ info (parseOpt <**> helper)
                            ( fullDesc
                           <> progDesc "Learning ARIMA"
                           <> header "backprop-learn-arima - backprop-learn demo"
                            )

    SomeNat (Proxy :: Proxy n) <- pure $ someNatVal oLookback
    Just Refl <- pure $ Proxy @1 `isLE` Proxy @n
    let generator = processConduit oProcess g
                 .| noisyConduit oNoise g

    case oMode of
      CSimulate lim ->
        runConduit $ generator
                  .| maybe (C.map id) C.take lim
                  .| C.mapM_ print
                  .| C.sinkNull

      CLearn (modelLearn->model) -> do
        let unrolled = trainState . unrollFinal @(SV.Vector n) $ model
        p0 <- initParamNormal unrolled 0.2 g

        let report n b = do
              liftIO $ printf "(Batch %d)\n" (b :: Int)
              t0 <- liftIO getCurrentTime
              C.drop (n - 1)
              mp <- mapM (liftIO . evaluate . force) =<< await
              t1 <- liftIO getCurrentTime
              case mp of
                Nothing -> liftIO $ putStrLn "Done!"
                Just p -> do
                  chnk <- lift . state $ (,[])
                  liftIO $ do
                    printf "Trained on %d points in %s.\n"
                      (length chnk)
                      (show (t1 `diffUTCTime` t0))
                    let trainScore = testModelAll absErrorTest unrolled (TJust p) chnk
                    printf "Training error:   %.8f\n" trainScore
                  report n (b + 1)

        flip evalStateT []
            . runConduit
            $ transPipe lift generator
           .| leadings
           .| skipSampling 0.05 g
           .| C.iterM (modify . (:))
           .| optoConduit def p0
                (adam def (modelGrad squaredError noReg unrolled))
           .| report oInterval 0
           .| C.sinkNull

parseOpt :: Parser Options
parseOpt = O <$> subparser ( command "simulate" (info (parseSim   <**> helper) (progDesc "Simulate ARIMA only"))
                          <> command "learn"    (info (parseLearn <**> helper) (progDesc "Simulate and learn model"))
                           )
             <*> option (maybeReader parseProcess)
                   ( long "process"
                  <> help "Process to learn"
                  <> showDefaultWith showProcess
                  <> value PSin
                  <> metavar "PROCESS"
                   )
             <*> option auto
                   ( short 'e'
                  <> help "Standard deviation of noise term"
                  <> showDefault
                  <> value 0.05
                  <> metavar "DOUBLE"
                   )
             <*> option auto
                   ( long "interval"
                  <> short 'i'
                  <> help "Report interval"
                  <> showDefault
                  <> value 5000
                  <> metavar "INT"
                   )
             <*> option auto
                   ( long "lookback"
                  <> short 'l'
                  <> help "Learn lookback"
                  <> showDefault
                  <> value 10
                  <> metavar "INT"
                   )
  where
    parseSim :: Parser Mode
    parseSim = CSimulate <$> optional (option auto
                                          ( short 'n'
                                         <> help "Number of items to generate (infinite, if none given)"
                                         <> metavar "INT"
                                          )
                                      )
    parseLearn :: Parser Mode
    parseLearn = subparser
        ( command "arima" (info parseARIMA (progDesc "Learn ARIMA(p,d,q) model"))
       <> command "fcrnn" (info parseFCRNN (progDesc "Learn Fully Connected RNN model"))
       <> command "lstm"  (info parseLSTM  (progDesc "Learn LSTM model"))
       <> command "gru"   (info parseGRU   (progDesc "Learn GRU model"))
        )
    parseARIMA :: Parser Mode
    parseARIMA = do
        p <- argument auto ( help "AR(p): Autoregressive lookback"
                          <> metavar "INT"
                           )
        d <- argument auto ( help "I(d): Differencing degree"
                          <> metavar "INT"
                           )
        q <- argument auto ( help "MA(q): Moving average lookback"
                          <> metavar "INT"
                           )
        pure $ withSomeSing @Nat p $ \sp@SNat ->
          withSomeSing @Nat d $ \sd@SNat ->
          withSomeSing @Nat q $ \sq@SNat ->
            CLearn (MARIMA sp sd sq)
    parseFCRNN :: Parser Mode
    parseFCRNN = do
        h <- argument auto ( help "Hidden layer size"
                          <> metavar "INT"
                          <> showDefault
                          <> value 10
                           )
        pure $ withSomeSing @Nat h $ \sh@SNat ->
          CLearn (MFCRNN sh)
    parseLSTM :: Parser Mode
    parseLSTM = do
        h <- argument auto ( help "Hidden layer size"
                          <> metavar "INT"
                          <> showDefault
                          <> value 10
                           )
        pure $ withSomeSing @Nat h $ \sh@SNat ->
          CLearn (MLSTM sh)
    parseGRU :: Parser Mode
    parseGRU = do
        h <- argument auto ( help "Hidden layer size"
                          <> metavar "INT"
                          <> showDefault
                          <> value 10
                           )
        pure $ withSomeSing @Nat h $ \sh@SNat ->
          CLearn (MGRU sh)

showProcess :: Process -> String
showProcess PSin = "sin"

parseProcess :: String -> Maybe Process
parseProcess s = case map toLower s of
    "sin" -> Just PSin
    _     -> Nothing

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
module AsteriskGaussianConvolution where

import           Control.Monad                  as M
import qualified Data.Array.Accelerate          as A
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.Repa                as R
import           Data.Binary
import           Data.Complex
import           Data.List                      as L
import           Data.Vector.Storable           as VS
import           Data.Vector.Unboxed            as VU
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.GreensFunction
import           FokkerPlanck.Histogram
import           Foreign.CUDA.Driver            as CUDA
import           FourierMethod.FourierSeries2D
import           FourierPinwheel
import           Image.IO
import           STC                            hiding (convolve)
import           System.Directory
import           System.Environment
import           System.FilePath
import           Utils.Array
import           Utils.List
import           Utils.Parallel
import           Utils.Time

main = do
  args@(deviceIDsStr:numPointsStr:deltaStr:thresholdStr:numPointsReconStr:deltaReconStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:tauStr:numR2FreqStr:periodR2Str:phiFreqsStr:rhoFreqsStr:thetaFreqsStr:scaleFreqsStr:initDistStr:initScaleStr:histFilePath:stdR2Str:stdStr:numBatchR2Str:numBatchR2FreqsStr:numBatchOriStr:batchSizeStr:sStr:radiusStr:numThreadStr:_) <-
    getArgs
  let deviceIDs = read deviceIDsStr :: [Int]
      numPoints = read numPointsStr :: Int
      delta = read deltaStr :: Double
      threshold = read thresholdStr :: Double
      numPointsRecon = read numPointsReconStr :: Int
      deltaRecon = read deltaReconStr :: Double
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      tau = read tauStr :: Double
      numR2Freq = read numR2FreqStr :: Int
      periodR2 = read periodR2Str :: Double
      phiFreq = read phiFreqsStr :: Int
      phiFreqs = L.map fromIntegral [-phiFreq .. phiFreq]
      rhoFreq = read rhoFreqsStr :: Int
      rhoFreqs = L.map fromIntegral [-rhoFreq .. rhoFreq]
      thetaFreq = read thetaFreqsStr :: Int
      thetaFreqs = L.map fromIntegral [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Int
      scaleFreqs = L.map fromIntegral [-scaleFreq .. scaleFreq]
      initScale = read initScaleStr :: Double
      initDist = read initDistStr :: [(Double, Double, Double, Double)]
      initPoints = L.map (\(x, y, t, s) -> Point x y t s) initDist
      initSource = [L.head initPoints]
      initSink = [L.last initPoints]
      numThread = read numThreadStr :: Int
      folderPath = "output/test/AsteriskGaussianConvolution"
      maxScale = read maxScaleStr :: Double
      halfLogPeriod = log maxScale
      stdR2 = read stdR2Str :: Double
      std = read stdStr :: Double
      numBatchR2 = read numBatchR2Str :: Int
      numBatchR2Freqs = read numBatchR2FreqsStr :: Int
      numBatchOri = read numBatchOriStr :: Int
      batchSize = read batchSizeStr :: Int
      s = read sStr :: Double
      radius = read radiusStr :: Double
      periodEnvelope = periodR2^2 / 2
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  initialise []
  devs <- M.mapM device deviceIDs
  ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
  ptxs <- M.mapM createTargetFromContext ctxs
  hist <-
    if flag
      then do
        printCurrentTime "Read coefficients..."
        decodeFile histFilePath
      else do
        printCurrentTime "Start computing coefficients..."
        sampleCartesian
          histFilePath
          folderPath
          ptxs
          numPoints
          periodEnvelope
          delta
          numOrientation
          initScale
          thetaSigma
          tau
          threshold
          s
          phiFreq
          rhoFreq
          thetaFreq
          scaleFreq
  printCurrentTime "Done"
  printCurrentTime "Start Convloution.."
  plan <-
    makePlan
      folderPath
      emptyPlan
      numPointsRecon
      numPointsRecon
      numR2Freq
      (2 * thetaFreq + 1)
      (2 * scaleFreq + 1)
      (2 * phiFreq + 1)
      (2 * rhoFreq + 1)
  let periodEnv = periodR2 * sqrt 2
      coefficients = getNormalizedHistogramArr hist
      harmonicsArray =
        createHarmonics
          numR2Freq
          phiFreq
          rhoFreq
          thetaFreq
          scaleFreq
          (-s)
          periodR2
          periodEnv
          coefficients
      asteriskGaussianVecs =
        asteriskGaussianFull
          numR2Freq
          thetaFreq
          scaleFreq
          (-s)
          periodR2
          periodEnv
          10
          10
      r2Freqs = getListFromNumber numR2Freq
      shiftFilter =
        fromListUnboxed
          (Z :. numR2Freq :. numR2Freq)
          [ L.foldl'
            (\b (Point x y theta scale) ->
               b +
               cis
                 (-(fromIntegral freqX * x + fromIntegral freqY * y) * 2 * pi /
                   periodR2))
            0
            initSource
          | freqY <- r2Freqs
          , freqX <- r2Freqs
          ]
      dftID = DFTPlanID DFT1DG [numR2Freq, numR2Freq] [0, 1]
      idftID = DFTPlanID IDFT1DG [numR2Freq, numR2Freq] [0, 1]
  shiftFilterF <-
    dftExecute plan dftID . VU.convert . toUnboxed . computeS . makeFilter2D $
    shiftFilter
  asteriskGaussianF <- dftExecuteBatchP plan dftID asteriskGaussianVecs
  filteredAsteriskGaussian <-
    fmap VS.concat .
    dftExecuteBatchP plan idftID . parMap rdeepseq (VS.zipWith (*) shiftFilterF) $
    asteriskGaussianF
  let numThetaFreq = 2 * thetaFreq + 1
      numRFreq = 2 * scaleFreq + 1
      sourceDist =
        FPArray
          numR2Freq
          numR2Freq
          (2 * scaleFreq + 1)
          (2 * thetaFreq + 1)
          (2 * rhoFreq + 1)
          (2 * phiFreq + 1) .
        parMap
          rdeepseq
          (\radialFreq ->
             VU.convert .
             toUnboxed .
             computeS .
             R.slice
               (fromUnboxed
                  (Z :. numRFreq :. numThetaFreq :. numR2Freq :. numR2Freq) .
                VS.convert $
                filteredAsteriskGaussian) $
             (Z :. radialFreq :. All :. All :. All)) $
        [0 .. 2 * scaleFreq]
  source <- convolve harmonicsArray sourceDist
  let sourceMat = A.transpose . toMatrixAcc $ source
  sourceR2' <-
    computeFourierSeriesR2StreamAcc
      ptxs
      (getFPArrayNumXFreq source)
      numPointsRecon
      (getFPArrayNumRFreq source * getFPArrayNumThetaFreq source)
      periodR2
      delta
      numBatchR2
      sourceMat
  let sourceR2 =
        R.reshape
          (Z :. numRFreq :. numThetaFreq :. numPointsRecon :. numPointsRecon)
          sourceR2'
  plotImageRepa (folderPath </> "Source.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numPointsRecon :. numPointsRecon) .
    VU.map sqrt .
    toUnboxed . sumS . sumS . R.map (\x -> (magnitude x) ** 2) . rotate4D2 $
    sourceR2

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 TemplateHaskell #-}
module Turtle2 (main) where
import Prelude hiding (dropWhile)
import Control.Applicative
import Control.Arrow
import Data.Complex
import AngleNum
import Ros.Node
import Ros.Topic (repeatM, force, dropWhile, metamorphM, yieldM)
import Ros.TopicUtil (everyNew, interruptible)
import Ros.Turtlesim.Pose
import Ros.Turtlesim.Velocity
import Ros.Logging
import System.IO (hFlush, stdout)

-- A type synonym for a 2D point.
type Point = Complex Float

-- A Topic of user-supplied waypoint trajectories.
getTraj :: Topic IO [Point]
getTraj = repeatM (do putStr "Enter waypoints: " >> hFlush stdout
                      $(logInfo "Waiting for new traj")
                      map (uncurry (:+)) . read <$> getLine)

-- Produce a new goal 'Point' every time a goal is reached.
destinations :: (Functor m, Monad m) =>
                Topic m Point -> Topic m Pose -> Topic m Point
destinations goals poses = metamorphM (start (p2v <$> poses)) goals
  where start t g = yieldM g (go g t)
        go g t g' = force (dropWhile (keepGoing g) t) >>= yieldM g' . go g'
        keepGoing goal pose = magnitude (goal - pose) > 1.5
        p2v (Pose x y _ _ _) = x :+ y

-- Compute linear distance to goal and bearing to goal
toGoal :: (Pose,Point) -> (Float, Angle Float)
toGoal (pos,goal) = (magnitude v, angle $ phase v)
  where v = goal - (x pos :+ y pos)

-- Steer based on a current pose estimate and distance from goal
steering :: Pose -> (Float, Angle Float) -> Velocity
steering pos (dpos, thetaDesired) = Velocity (min 2 dpos) angVel
  where thetaErr = toDegrees $ thetaDesired - angle (theta pos)
        angVel = signum thetaErr * min 2 (abs thetaErr)

navigate :: (Pose, Point) -> Velocity
navigate = uncurry ($) . (steering . fst &&& toGoal)

main = runNode "HaskellBTurtle" $
       do enableLogging (Just Warn)
          poses <- subscribe "/turtle1/pose"
          let goals = destinations (interruptible getTraj) poses
          advertise "/turtle1/command_velocity"
                    (navigate <$> everyNew poses goals)

max_stars_repo_name: rgleichman/roshask

max_issues_repo_name: rgleichman/roshask

max_forks_repo_name: rgleichman/roshask


text:
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Main (main) where

import Injection
import Projection

import Data.Complex (Complex ((:+)))
import Data.Dynamic (Dynamic)
import Data.Fixed (Fixed, E6)
import Data.Functor.Const (Const)
import Data.Functor.Identity (Identity)
import Data.List.NonEmpty (NonEmpty ((:|)))
import Data.Map (Map)
import qualified Data.Map as Map
import Data.Monoid (Dual)
import Data.Monoid (Product)
import Data.Monoid (Sum)
import Data.Monoid (All, Any)
import qualified Data.Monoid as Monoid (First, Last)
import Data.Ord (Down (..))
import Data.Ratio (Ratio, (%))
import Data.Semigroup (Max, Min)
import qualified Data.Semigroup as Semigroup (First, Last)
import Data.Set (Set)
import qualified Data.Set as Set
import Data.Text (Text)
import qualified Data.Text.Lazy as Lazy (Text)
import Numeric.Natural (Natural)
import Test.Hspec
import Test.QuickCheck hiding (Fixed)
import Test.QuickCheck.Instances ()

main :: IO ()
main = hspec $ do
    describe "instance Injection Integer Integer" $ do
        it "is resolvable" (resolveInjection @Integer @Integer)
    describe "instance Injection String String" $ do
        it "is resolvable" (resolveInjection @String @String)
    describe "instance Injection Text Text" $ do
        it "is resolvable" (resolveInjection @Text @Text)
    describe "instance Injection Lazy.Text Lazy.Text" $ do
        it "is resolvable" (resolveInjection @Lazy.Text @Lazy.Text)
    describe "instance Retraction Integer Dynamic" $ do
        it "is resolvable" (resolveRetraction @Integer @Dynamic)
        it "is the left inverse of inject" (lawLeftInverse @Integer @Dynamic)
    describe "instance Injection (Maybe Integer) [Integer]" $ do
        it "is resolvable" (resolveInjection @(Maybe Integer) @[Integer])
        it "outputs a value with the same length as the input" $ do
            let test = inject @(Maybe Integer) @[Integer]
            length (test Nothing) `shouldBe` 0
            length (test (Just 1)) `shouldBe` 1
        it "is injective" (lawInjective @(Maybe Integer) @[Integer])
    describe "instance Retraction (Maybe Integer) [Integer]" $ do
        it "is resolvable" (resolveRetraction @(Maybe Integer) @[Integer])
        it "is the left inverse of inject" (lawLeftInverse @(Maybe Integer) @[Integer])
    describe "instance Injection Natural Integer" $ do
        it "is resolvable" (resolveInjection @Natural @Integer)
        it "injects the value itself" $ do
            let test = inject @Natural @Integer
            test 0 `shouldBe` 0
            test 1 `shouldBe` 1
            test 2 `shouldBe` 2
        it "is injective" (lawInjective @Natural @Integer)
    describe "instance Retraction Natural Integer" $ do
        it "is resolvable" (resolveRetraction @Natural @Integer)
        it "is the left inverse of inject" (lawLeftInverse @Natural @Integer)
        it "is defined over non-negative integers" $ do
            retract @Natural @Integer 0 `shouldBe` Just 0
            retract @Natural @Integer 1 `shouldBe` Just 1
        it "is not defined over negative integers" $ do
            retract @Natural @Integer (-1) `shouldBe` Nothing
    describe "instance Injection Text String" $ do
        it "is resolvable" (resolveInjection @Text @String)
        it "is injective" (lawInjective @Text @String)
    describe "instance Injection Lazy.Text String" $ do
        it "is resolvable" (resolveInjection @Lazy.Text @String)
        it "is injective" (lawInjective @Lazy.Text @String)
    describe "instance Injection Text Lazy.Text" $ do
        it "is resolvable" (resolveInjection @Text @Lazy.Text)
        it "is injective" (lawInjective @Text @Lazy.Text)
    describe "instance Injection Lazy.Text Text" $ do
        it "is resolvable" (resolveInjection @Lazy.Text @Text)
        it "is injective" (lawInjective @Lazy.Text @Text)
    describe "instance Injection Integer (Fixed a)" $ do
        it "is resolvable" (resolveInjection @Integer @(Fixed E6))
        it "is injective" (lawInjective @Integer @(Fixed E6))
    describe "instance Retraction Integer (Fixed a)" $ do
        it "is resolvable" (resolveRetraction @Integer @(Fixed E6))
        it "is the left inverse of inject" (lawLeftInverse @Integer @(Fixed E6))
        it "is defined over integers" $ do
            retract @Integer @(Fixed E6) 0 `shouldBe` Just 0
            retract @Integer @(Fixed E6) 1 `shouldBe` Just 1
            retract @Integer @(Fixed E6) (-1) `shouldBe` Just (-1)
        it "is not defined over fractions" $ do
            retract @Integer @(Fixed E6) (1 / 2) `shouldBe` Nothing
    describe "instance Injection Integer (Const Integer String)" $ do
        it "is resolvable" (resolveInjection @Integer @(Const Integer String))
        it "is injective" (lawInjective @Integer @(Const Integer String))
    describe "instance Injection (Const Integer String) Integer" $ do
        it "is resolvable" (resolveInjection @(Const Integer String) @Integer)
        it "is injective" (lawInjective @(Const Integer String) @Integer)
    describe "instance Injection Integer (Ratio Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Ratio Integer))
        it "is injective" (lawInjective @Integer @(Ratio Integer))
    describe "instance Retraction Integer (Ratio Integer)" $ do
        it "is resolvable" (resolveRetraction @Integer @(Ratio Integer))
        it "is the left inverse of inject" (lawLeftInverse @Integer @(Ratio Integer))
        it "is defined over integers" $ do
            retract @Integer @(Ratio Integer) (0 % 1) `shouldBe` Just 0
            retract @Integer @(Ratio Integer) (1 % 1) `shouldBe` Just 1
            retract @Integer @(Ratio Integer) (1 % (-1)) `shouldBe` Just (-1)
        it "is not defined over fractions" $ do
            retract @Integer @(Ratio Integer) (1 % 2) `shouldBe` Nothing
    describe "instance Injection Double (Complex Double)" $ do
        it "is resolvable" (resolveInjection @Double @(Complex Double))
        it "is injective" (lawInjective @Double @(Complex Double))
    describe "instance Retraction Double (Complex Double)" $ do
        it "is resolvable" (resolveRetraction @Double @(Complex Double))
        it "is the left inverse of inject" (lawLeftInverse @Double @(Complex Double))
        it "is defined over real numbers" $ do
            retract @Double @(Complex Double) (0 :+ 0) `shouldBe` Just 0
            retract @Double @(Complex Double) (1 :+ 0) `shouldBe` Just 1
            retract @Double @(Complex Double) ((-1) :+ 0) `shouldBe` Just (-1)
        it "is not defined over imaginary and complex numbers" $ do
            retract @Double @(Complex Double) (0 :+ 1) `shouldBe` Nothing
            retract @Double @(Complex Double) (1 :+ 1) `shouldBe` Nothing
    describe "instance Injection Integer (Identity Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Identity Integer))
        it "is injective" (lawInjective @Integer @(Identity Integer))
    describe "instance Injection (Identity Integer) Integer" $ do
        it "is resolvable" (resolveInjection @(Identity Integer) @Integer)
        it "is injective" (lawInjective @(Identity Integer) @Integer)
    describe "instance Injection (NonEmpty Integer) [Integer]" $ do
        it "is resolvable" (resolveInjection @(NonEmpty Integer) @[Integer])
        it "is injective" (lawInjective @(NonEmpty Integer) @[Integer])
    describe "instance Retraction (NonEmpty Integer) [Integer]" $ do
        it "is resolvable" (resolveRetraction @(NonEmpty Integer) @[Integer])
        it "is the left inverse of inject" (lawLeftInverse @(NonEmpty Integer) @[Integer])
        it "is not defined on the empty list" $ do
            retract @(NonEmpty Integer) @[Integer] [] `shouldBe` Nothing
        it "is defined over non-empty lists" $ do
            retract @(NonEmpty Integer) @[Integer] [0] `shouldBe` Just (0 :| [])
            retract @(NonEmpty Integer) @[Integer] [0, 1] `shouldBe` Just (0 :| [1])
    describe "instance Injection Integer (Down Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Down Integer))
        it "is injective" (lawInjective @Integer @(Down Integer))
    describe "instance Injection (Down Integer) Integer" $ do
        it "is resolvable" (resolveInjection @(Down Integer) @Integer)
        it "is injective" (lawInjective @(Down Integer) @Integer)
    describe "instance Injection Integer (Product Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Product Integer))
        it "is injective" (lawInjective @Integer @(Product Integer))
    describe "instance Injection (Product Integer) Integer" $ do
        it "is resolvable" (resolveInjection @(Product Integer) @Integer)
        it "is injective" (lawInjective @(Product Integer) @Integer)
    describe "instance Injection Integer (Sum Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Sum Integer))
        it "is injective" (lawInjective @Integer @(Sum Integer))
    describe "instance Injection (Sum Integer) Integer" $ do
        it "is resolvable" (resolveInjection @(Sum Integer) @Integer)
        it "is injective" (lawInjective @(Sum Integer) @Integer)
    describe "instance Injection Integer (Dual Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Dual Integer))
        it "is injective" (lawInjective @Integer @(Dual Integer))
    describe "instance Injection (Dual Integer) Integer" $ do
        it "is resolvable" (resolveInjection @(Dual Integer) @Integer)
        it "is injective" (lawInjective @(Dual Integer) @Integer)
    describe "instance Injection Integer (Monoid.Last Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Monoid.Last Integer))
        it "is injective" (lawInjective @Integer @(Monoid.Last Integer))
    describe "instance Retraction Injection (Monoid.Last Integer)" $ do
        it "is resolvable" (resolveRetraction @Integer @(Monoid.Last Integer))
        it "is the left inverse of inject" (lawLeftInverse @Integer @(Monoid.Last Integer))
        it "is not defined over mempty" $ do
            retract @Integer @(Monoid.Last Integer) mempty `shouldBe` Nothing
    describe "instance Injection Integer (Monoid.First Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Monoid.First Integer))
        it "is injective" (lawInjective @Integer @(Monoid.First Integer))
    describe "instance Retraction Injection (Monoid.First Integer)" $ do
        it "is resolvable" (resolveRetraction @Integer @(Monoid.First Integer))
        it "is the left inverse of inject" (lawLeftInverse @Integer @(Monoid.First Integer))
        it "is not defined over mempty" $ do
            retract @Integer @(Monoid.First Integer) mempty `shouldBe` Nothing
    describe "instance Injection Integer (Semigroup.First Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Semigroup.First Integer))
        it "is injective" (lawInjective @Integer @(Semigroup.First Integer))
    describe "instance Injection (Semigroup.First Integer) Integer" $ do
        it "is resolvable" (resolveInjection @(Semigroup.First Integer) @Integer)
        it "is injective" (lawInjective @(Semigroup.First Integer) @Integer)
    describe "instance Injection Integer (Semigroup.Last Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Semigroup.Last Integer))
        it "is injective" (lawInjective @Integer @(Semigroup.Last Integer))
    describe "instance Injection (Semigroup.Last Integer) Integer" $ do
        it "is resolvable" (resolveInjection @(Semigroup.Last Integer) @Integer)
        it "is injective" (lawInjective @(Semigroup.Last Integer) @Integer)
    describe "instance Injection Integer (Max Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Max Integer))
        it "is injective" (lawInjective @Integer @(Max Integer))
    describe "instance Injection (Max Integer) Integer" $ do
        it "is resolvable" (resolveInjection @(Max Integer) @Integer)
        it "is injective" (lawInjective @(Max Integer) @Integer)
    describe "instance Injection Integer (Min Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(Min Integer))
        it "is injective" (lawInjective @Integer @(Min Integer))
    describe "instance Injection (Min Integer) Integer" $ do
        it "is resolvable" (resolveInjection @(Min Integer) @Integer)
        it "is injective" (lawInjective @(Min Integer) @Integer)
    describe "instance Injection Integer (String -> Integer)" $ do
        it "is resolvable" (resolveInjection @Integer @(String -> Integer))
    describe "instance Injection Bool Any" $ do
        it "is resolvable" (resolveInjection @Bool @Any)
        it "is injective" (lawInjective @Bool @Any)
    describe "instance Injection Any Bool" $ do
        it "is resolvable" (resolveInjection @Any @Bool)
        it "is injective" (lawInjective @Any @Bool)
    describe "instance Injection Bool All" $ do
        it "is resolvable" (resolveInjection @Bool @All)
        it "is injective" (lawInjective @Bool @All)
    describe "instance Injection All Bool" $ do
        it "is resolvable" (resolveInjection @All @Bool)
        it "is injective" (lawInjective @All @Bool)
    describe "instance Projection Integer Integer" $ do
        it "is resolvable" (resolveProjection @Integer @Integer)
        it "is surjective" (lawSurjective @Integer @Integer id)
    describe "instance Section Integer Integer" $ do
        it "is resolvable" (resolveSection @Integer @Integer)
        it "is the right inverse of project" (lawRightInverse @Integer @Integer)
    describe "instance Projection [(Integer, Integer)] (Map Integer Integer)" $ do
        it "is resolvable" (resolveProjection @[(Integer, Integer)] @(Map Integer Integer))
        it "is surjective" (lawSurjective @[(Integer, Integer)] @(Map Integer Integer) Map.toList)
    describe "instance Projection [Integer] (Set Integer)" $ do
        it "is resolvable" (resolveProjection @[Integer] @(Set Integer))
        it "is surjective" (lawSurjective @[Integer] @(Set Integer) Set.toList)

resolveInjection :: forall from into. Injection from into => Expectation
resolveInjection = seq (inject @from @into) return ()

resolveRetraction :: forall from into. Retraction from into => Expectation
resolveRetraction = seq (retract @from @into) return ()

resolveProjection :: forall whole part. Projection whole part => Expectation
resolveProjection = seq (project @whole @part) return ()

resolveSection :: forall whole part. Section whole part => Expectation
resolveSection = seq (update @whole @part) return ()

lawInjective
    :: forall from into
    .  Injection from into
    => Arbitrary from
    => (Eq from, Show from)
    => (Eq into, Show into)
    => Property
lawInjective = property $ \from1 from2 ->
    let into1 = inject @from @into from1
        into2 = inject @from @into from2
    in (from1 /= from2) ==> (into1 =/= into2)

lawLeftInverse
    :: forall from into
    .  Retraction from into
    => Arbitrary from
    => (Eq from, Show from)
    => Property
lawLeftInverse = property $ \from ->
    let into = inject @from @into from
        from' = retract @from @into into
    in Just from === from'

lawSurjective
    :: forall whole part
    .  Projection whole part
    => Arbitrary part
    => (Eq part, Show part)
    => (part -> whole)  -- ^ any total right inverse of @project@
    -> Property
lawSurjective rightInverse = property $ \part ->
    -- This test actually proves nothing: that a right inverse of 'project'
    -- merely exists proves that it is surjective. Here, we are just probing the
    -- putative right inverse.
    project (rightInverse part) === part

lawRightInverse
    :: forall whole part
    .  Section whole part
    => Arbitrary whole
    => Show whole
    => Arbitrary part
    => (Eq part, Show part)
    => Property
lawRightInverse = property $ \part whole ->
    project @whole @part (update whole part) === part

instance Arbitrary a => Arbitrary (Down a) where
    arbitrary = Down <$> arbitrary

max_stars_repo_name: ttuegel/injection

max_issues_repo_name: ttuegel/injection

max_forks_repo_name: ttuegel/injection


text:
module BuiltIn.Stats (builtInStats) where

import qualified Data.Map as M
import qualified Data.Vector as V
import           Numeric.SpecFunctions
import           Statistics.Distribution
import           Statistics.Distribution.Binomial
import           Statistics.Distribution.ChiSquared
import           Statistics.Distribution.Exponential
import           Statistics.Distribution.FDistribution
import           Statistics.Distribution.Gamma
import           Statistics.Distribution.Hypergeometric
import           Statistics.Distribution.StudentT
import qualified Statistics.Sample as SS

import           BuiltIn.Common
import           Type
import           Types

-- | Truncate a double back to a double
dTrunc :: Double -> Double
dTrunc x = fromIntegral (truncate x :: Integer)

biAveDev :: BuiltIn
biAveDev = BuiltIn
    { evalVal = aveDev
    , emitVal = const "AVEDEV"
    , typeSig = typeList typeNum :-> typeNum
    , argHelp = "[numbers]"
    , addHelp = Nothing
    }

aveDev :: [Value] -> Interpreter EvalError Value
aveDev [VList xs]
    | length xs > 255 = biErr "Too many arguments (>255)"
    | otherwise = do
        xs' <- vListToDbls xs
        let mn = SS.mean $ V.fromList xs'
        let var = fmap (\x -> abs (mn - x)) xs'
        return . VNum . SS.mean $ V.fromList var
aveDev _ = infError

biAvg :: BuiltIn
biAvg = BuiltIn
    { evalVal = avg
    , emitVal = const "AVERAGE"
    , typeSig = typeList typeNum :-> typeNum
    , argHelp = "[numbers]"
    , addHelp = Nothing
    }

avg :: [Value] -> Interpreter EvalError Value
avg [VList xs]
    | length xs > 255 = biErr "Too many arguments (>255)"
    | otherwise = do
        xs' <- vListToDbls xs
        return . VNum . SS.mean $ V.fromList xs'
avg _ = infError

biBinom :: BuiltIn
biBinom = BuiltIn
    { evalVal = binom
    , emitVal = const "BINOMDIST"
    , typeSig = typeNum :-> typeNum :-> typeNum :-> typeBool :-> typeNum
    , argHelp = "number_s trials probability_s cumulative"
    , addHelp = Nothing
    }

binom :: [Value] -> Interpreter EvalError Value
binom [VNum num, VNum trls, VNum prob, VBool cum]
    | prob < 0 || prob > 1 = biErr "Probability in binomdist must be between 0 and 1"
    | num  < 0 = biErr "Number of successes in binomdist cannot be negative"
    | num > trls = biErr "Successes > trials is not allowed in binomdist"
    | otherwise = if cum
        then return . VNum $ cumulative  dist (dTrunc num)
        else return . VNum $ probability dist (truncate num)
  where
    dist   = binomial (truncate trls) prob
binom _         = infError

biChiDist :: BuiltIn
biChiDist = BuiltIn
    { evalVal = chiDist
    , emitVal = const "CHIDIST"
    , typeSig = typeNum :-> typeNum :-> typeNum
    , argHelp = "x degrees_freedom"
    , addHelp = Nothing
    }

chiDist :: [Value] -> Interpreter EvalError Value
chiDist [VNum x, VNum df]
    | x < 0 = biErr "First argument to chidist must be a positive number"
    | df < 1  = biErr "Degrees of freedom for chidist must be 1 or greater"
    | df > 10**10 = biErr "Degrees of freedom for chidist too large"
    | otherwise = return . VNum $ complCumulative dist x
  where
    dist = chiSquared (truncate df)
chiDist _ = infError

{-
biChiInv :: BuiltIn
biChiInv = BuiltIn
    { evalVal = chiInv
    , emitVal = const "CHIINV""
    , typeSig = typeNum :-> typeNum :-> typeNum
    , argHelp = "probability degrees_freedom"
    , addHelp = Nothing
    }

chiInv :: [Value] -> Interpreter EvalError Value
chiInv [VNum prob, VNum df]
    | prob < 0 || prob > 1 = biErr "Probability in chiinv must be between 0 and 1"
    | df < 1  = biErr "Degrees of freedom for chiinv must be 1 or greater"
    | df > 10**10 = biErr "Degrees of freedom for chiinv too large"
    | otherwise = return . VNum $ cumulative dist prob
  where
    dist = chiSquared (truncate df)
chiInv _ = infError
-}

biConfd :: BuiltIn
biConfd = BuiltIn
    { evalVal = confd
    , emitVal = const "CONFIDENCE"
    , typeSig = typeNum :-> typeNum :-> typeNum :-> typeNum
    , argHelp = "alpha standard_dev size"
    , addHelp = Nothing
    }

confd :: [Value] -> Interpreter EvalError Value
confd [VNum alph, VNum std, VNum size]
    | alph <= 0 || alph >= 1 = biErr "Alpha for confidence must be between 0 and 1"
    | size < 1 = biErr "Sample size in confidence must be greater than one"
    | std <= 0 = biErr "Standard deviation for confidence must be >= 0"
    | otherwise = return $ VNum $ cv * se
  where
    se = std / sqrt (dTrunc size)
    df = size - 1
    cp = 1 - alph/2
    cv = quantile (studentT df) cp
confd _ = infError

biExpDist :: BuiltIn
biExpDist = BuiltIn
    { evalVal = expDist
    , emitVal = const "EXPONDIST"
    , typeSig = typeNum :-> typeNum :-> typeNum :-> typeNum
    , argHelp = "x lambda cumulative"
    , addHelp = Nothing
    }

expDist :: [Value] -> Interpreter EvalError Value
expDist [VNum x, VNum lam, VBool cum]
    | x < 0 = biErr "First paramter to expondist must be >= 0"
    | lam <= 0 = biErr "Lambda must be > 0 in expondist"
    | otherwise = if cum
        then return . VNum $ cumulative dist x
        else return . VNum $ density dist x
  where
    dist = exponential lam
expDist _ = infError

biFDist :: BuiltIn
biFDist = BuiltIn
    { evalVal = fDist
    , emitVal = const "FDIST"
    , typeSig = typeNum :-> typeNum :-> typeNum :-> typeNum
    , argHelp = "x degrees_freedom1 degrees_freedom2"
    , addHelp = Nothing
    }

fDist :: [Value] -> Interpreter EvalError Value
fDist [VNum x, VNum df1, VNum df2]
    | x < 0 = biErr "First argument to fdist must be positive"
    | df1 < 1 || df1 >= 10**10 = biErr "fdist numerator degrees of freedom out of range"
    | df2 < 1 || df2 >= 10**10 = biErr "fdist denominator degrees of freedom out of range"
    | otherwise = return . VNum $ complCumulative dist x
  where
    dist = fDistribution (truncate df1) (truncate df2)
fDist _ = infError

biFisher :: BuiltIn
biFisher = BuiltIn
    { evalVal = fisher
    , emitVal = const "FISHER"
    , typeSig = typeNum :-> typeNum
    , argHelp = "x"
    , addHelp = Nothing
    }

fisher :: [Value] -> Interpreter EvalError Value
fisher [VNum x]
    | x >= 1 || x <= -1 = biErr "fisher requires an argument between -1 and 1"
    | otherwise = return . VNum $ atanh x
fisher _ = infError

biFisherInv :: BuiltIn
biFisherInv = BuiltIn
    { evalVal = fisherinv
    , emitVal = const "FISHERINV"
    , typeSig = typeNum :-> typeNum
    , argHelp = "y"
    , addHelp = Nothing
    }

fisherinv :: [Value] -> Interpreter EvalError Value
fisherinv [VNum x] = return . VNum $ tanh x
fisherinv _ = infError

biGamDist :: BuiltIn
biGamDist = BuiltIn
    { evalVal = gamDist
    , emitVal = const "GAMMADIST"
    , typeSig = typeNum :-> typeNum :-> typeNum :-> typeBool :-> typeNum
    , argHelp = "x alpha beta cumulative"
    , addHelp = Nothing
    }

gamDist :: [Value] -> Interpreter EvalError Value
gamDist [VNum x, VNum alph, VNum beta, VBool cum]
    | x < 0 = biErr "gamma requires its first argument to be positive"
    | alph <= 0 = biErr "gammadist requires alpha > 0"
    | beta <= 0 = biErr "gammadist requires beta > 0"
    | otherwise = if cum
        then return . VNum $ cumulative dist x
        else return . VNum $ density dist x
  where
    dist = gammaDistr alph beta
gamDist _ = infError

biGamInv :: BuiltIn
biGamInv = BuiltIn
    { evalVal = gamInv
    , emitVal = const "GAMMAINV"
    , typeSig = typeNum :-> typeNum :-> typeNum :-> typeNum
    , argHelp = "probability alpha beta"
    , addHelp = Nothing
    }

gamInv :: [Value] -> Interpreter EvalError Value
gamInv [VNum p, VNum alph, VNum beta]
    | p < 0 || p > 1 = biErr "gammainv requires probability to be >= 0 and <= 1"
    | alph <= 0 = biErr "gammainv requires alpha > 0"
    | beta <= 0 = biErr "gammainv requires beta > 0"
    | otherwise = return . VNum $ quantile dist p
  where
    dist = gammaDistr alph beta
gamInv _ = infError

biGamLn :: BuiltIn
biGamLn = BuiltIn
    { evalVal = gamLn
    , emitVal = const "GAMMALN"
    , typeSig = typeNum :-> typeNum
    , argHelp = "x"
    , addHelp = Nothing
    }

gamLn :: [Value] -> Interpreter EvalError Value
gamLn [VNum x]
    | x <= 0 = biErr "gammaln requires its argument be > 0"
    | otherwise = return . VNum $ logGammaL x
gamLn _ = infError

biHGeoDist :: BuiltIn
biHGeoDist = BuiltIn
    { evalVal = hgeoDist
    , emitVal = const "HYPGEOMDIST"
    , typeSig = typeNum :-> typeNum :-> typeNum :-> typeNum :-> typeNum
    , argHelp = "sample_s number_sample population_s number_population"
    , addHelp = Nothing
    }

hgeoDist :: [Value] -> Interpreter EvalError Value
hgeoDist [VNum sam, VNum num, VNum pop, VNum psiz]
    | sam < 0 = biErr "hypgeomdist sample_s must be >=0"
    | sam > min num pop = biErr "hypgeomdist sample_s too large (> number_sample or population_s)"
    | num <= 0 = biErr "hypgeomdist number_sample must be > 0"
    | num > psiz = biErr "hypgeomdist number_sample too large (> number_population)"
    | pop <= 0 = biErr "hypgeomdist population_s must be > 0"
    | pop > psiz = biErr "hypgeomdist population_s too large (> number_population)"
    | psiz <= 0 = biErr "hypgeomdist number_population must be > 0"
    | otherwise = return . VNum $ probability dist $ truncate sam
  where
    dist = hypergeometric (truncate pop) (truncate psiz) (truncate num) -- psiz
hgeoDist _ = infError

builtInStats :: M.Map String BuiltIn
builtInStats = M.fromList
    [ ("avedev",      biAveDev)
    , ("average",     biAvg)
    , ("binomdist",   biBinom)
    , ("chidist",     biChiDist)
--  , ("chiinv",      biChiInv)
    , ("confidence",  biConfd)
    , ("expondist",   biExpDist)
    , ("fdist",       biFDist)
    , ("fisher",      biFisher)
    , ("fisherinv",   biFisherInv)
    , ("gammadist",   biGamDist)
    , ("gammainv",    biGamInv)
    , ("gammaln",     biGamLn)
    , ("hypgeomdist", biHGeoDist)
    ]

max_stars_repo_name: ahodgen/archer-calc

max_issues_repo_name: ahodgen/archer-calc

max_forks_repo_name: ahodgen/archer-calc


text:
-----------------------------------------------------------------------------
{- |
Module      :  Numeric.LinearAlgebra
Copyright   :  (c) Alberto Ruiz 2006-10
License     :  GPL-style

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional
Portability :  uses ffi

This module reexports all normally required functions for Linear Algebra applications.

It also provides instances of standard classes 'Show', 'Read', 'Eq',
'Num', 'Fractional', and 'Floating' for 'Vector' and 'Matrix'.
In arithmetic operations one-component vectors and matrices automatically
expand to match the dimensions of the other operand.

-}
-----------------------------------------------------------------------------
module Numeric.LinearAlgebra (
    module Numeric.Container,
    module Numeric.LinearAlgebra.Algorithms
) where

import Numeric.Container
import Numeric.LinearAlgebra.Algorithms
import Numeric.Matrix()
import Numeric.Vector()

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE FlexibleContexts #-}

module Main where

import Graphics.Plot

import Data.Packed.Vector
import Data.Packed.Matrix
import Numeric.Container
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.HMatrix (tr)

import System.Environment (getArgs)

-- TODO (a->b->c) -> Vector a -> Vector b -> Matrix c
outerOp :: (Element a, Element b) => (Matrix a -> Matrix b -> Matrix c) -> Vector a -> Vector b -> Matrix c
outerOp k x x' = k xp xq
    where xp = repmat (asColumn x) 1 (dim x')
          xq= repmat (asRow x') (dim x) 1

k_SE sigma l x x' = (sigma^2*) $ exp $ (-1/(2*l^2))*dx^2
    where dx = outerOp (-) x x'

k = k_SE 3 0.3

inferenceTrivial k x y ySigma2 domain = (mean, variance)
    where mean = k domain x <> invkxx <> y
          variance = k domain domain - (kdomainx <> invkxx <> k x domain)
          invkxx = inv $ k x x + diag ySigma2
          kdomainx = k domain x

inferenceCholesky k x y ySigma2 domain = (mean, variance)
    where mean = (tr kxdomain) <> alpha
          variance = k domain domain - (tr v <> v)
          alpha = tr l <\> (l <\> y)
          l = chol $ k x x + diag ySigma2
          v = l <\> kxdomain
          kxdomain = k x domain

plotProcessDistribution :: Vector Double -> Vector Double -> Vector Double -> IO ()
plotProcessDistribution domain m v = mplot [domain, m, m + v, m - v, m + 2 * v, m - 2 * v]

main :: IO ()
main = do
    datfile:[] <- getArgs
    dat <- fmap readMatrix $ readFile datfile
    let [x, y] = toColumns dat -- [x,y] or [x,f]
    let domain = linspace 128 (0.5, 10.5)

    let n = 1
    let (x',y') = (vjoin $ replicate n x, vjoin $ take n $ repeat y)

    let (m, v') = inferenceTrivial k x' y' (constant 1 $ dim y') domain
    let v = takeDiag v'

    plotProcessDistribution domain m v

max_stars_repo_name: Jim-Holmstroem/gaussianprocess

max_issues_repo_name: Jim-Holmstroem/gaussianprocess

max_forks_repo_name: Jim-Holmstroem/gaussianprocess


text:
module SimpleParser where

import SimpleDefs
import Control.Monad (liftM)
import Text.ParserCombinators.Parsec hiding (spaces)
import Numeric (readFloat, readHex, readOct)
import Data.Complex (Complex((:+)))
import Data.Char (toLower)
import Data.Ratio (Rational, (%))

symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=?>@^_~#"

spaces :: Parser ()
spaces = skipMany1 space

{-
parseChar :: Parser LispVal
parseChar = do string "#\\"
               c <- many1 letter
               return $ Char $ case map toLower c of
                   "space" -> ' '
                   "newline" -> '\n'
                   [x] -> x
-}

parseString :: Parser LispVal
parseString = do char '"'
                 x <- many $ parseEscape <|> (noneOf "\\\"")
                 char '"'
                 return $ String x

parseEscape :: Parser Char
parseEscape = do char '\\'
                 c <- oneOf "\\\"nrt"
                 return $ case c of
                     '\\' -> c
                     '\"' -> c
                     'n' -> '\n'
                     'r' -> '\r'
                     't' -> '\t'

parseAtom :: Parser LispVal
parseAtom = do first <- letter <|> symbol
               rest <- many $ letter <|> digit <|> symbol
               let atom = first : rest
               return $ case atom of
                 "#t" -> Bool True
                 "#f" -> Bool False
                 _ -> Atom atom

parseSimpleNumber :: Parser LispVal
parseSimpleNumber = many1 digit >>= return . Number . read

{-
parseRadixNumber :: Parser LispVal
parseRadixNumber = char '#' >>
                (   parseDecimal
                <|> parseOctal
                <|> parseHexadecimal
                <|> parseBinary
                )

parseDecimal, parseHexadecimal, parseBinary, parseOctal :: Parser LispVal
parseDecimal = do char 'd'
                  number <- many1 digit
                  return . Number . read $ number

parseHexadecimal = do char 'x'
                      number <- many $ digit <|> (oneOf "abcdef")
                      return . Number . fst . head . readHex $ number

parseBinary = do char 'b'
                 number <- many1 . oneOf $ "01"
                 return . Number . readBin $ number

parseOctal = do char 'o'
                number <- many1 . oneOf $ "01234567"
                return . Number . fst . head . readOct $ number

readBin :: String -> Integer
readBin bin = fst $ foldl (\(n, ix) x -> (n + (num x) * (2 ^ ix), ix - 1))
                          (0, (length bin) - 1)
                          bin where
                  num '0' = 0
                  num '1' = 1
-}
parseNumber :: Parser LispVal
parseNumber = parseSimpleNumber -- <|> parseRadixNumber

{-
parseFloat :: Parser LispVal
parseFloat = do before <- many1 digit
                char '.'
                after <- many1 digit
                return . Float . fst . head . readFloat $ before ++ '.' : after

parseRational :: Parser LispVal
parseRational = do before <- fmap read $ many1 digit
                   char '/'
                   after <- fmap read $ many1 digit
                   return . Rational $ (before % after)

parseComplex :: Parser LispVal
parseComplex = do real <- fmap toFloat $ parseFloat <|> parseSimpleNumber
                  char '+'
                  imag <- fmap toFloat $ parseFloat <|> parseSimpleNumber
                  char 'i'
                  return . Complex $ (real :+ imag) where
                      toFloat (Float x) = x
                      toFloat (Number x) = fromInteger x
-}

parseQuoted :: Parser LispVal
parseQuoted = do char '\''
                 x <- parseExpr
                 return $ List [Atom "quote", x]

parseList :: Parser LispVal
parseList = liftM List $ sepBy parseExpr spaces

parseDottedList :: Parser LispVal
parseDottedList = do head <- endBy parseExpr spaces
                     char '.'
                     many1 spaces
                     tail <- parseExpr
                     return $ DottedList head tail

parseExpr :: Parser LispVal
parseExpr = parseAtom
        -- <|> parseBool
        -- <|> parseChar
        <|> parseString
        -- <|> parseComplex
        -- <|> parseRational
        -- <|> parseFloat
        <|> parseNumber
        <|> parseQuoted
        <|> do char '('
               x <- (try parseList) <|> parseDottedList
               char ')'
               return x


text:
{-# LANGUAGE DataKinds, TypeApplications, RankNTypes, FlexibleContexts,
             PartialTypeSignatures, TypeFamilies, TypeOperators, ScopedTypeVariables #-}

import qualified Control.Foldl as L
import qualified Numeric.LinearAlgebra.Static as H
import Numeric.LinearAlgebra.Data
import GHC.TypeNats
import Data.List
import Numeric.Backprop
import System.Random
import System.Random.Shuffle
import Networks
import Optimierung
import MNIST
import Tuple
import Data.Proxy

main :: IO ()
main = testDigits

antwort :: KnownNat n => H.R n -> Int
antwort = maxIndex . H.extract

normalizeInput :: (KnownNat n, KnownNat m)
               => (H.R n, H.R m)
               -> (H.R n, H.R m)
normalizeInput (x, y) = (x / 255, y)

multilayerPerceptron :: forall i o. KN i o => Modell _ (H.R i) (H.R o)
multilayerPerceptron = sigmoidNoBias @30 <~ sigmoid

threeLayer :: forall i o n. (KN i o, KnownNat n) => Proxy n -> Modell _ (H.R i) (H.R o)
threeLayer _ = sigmoid @n <~ sigmoid

mP' :: (KN i o) => Modell _ (H.R i) (H.R o)
mP' = leakyReLU @30 <~ leakyReLU

emptyLn = putStr "\n"

randomOS :: (Fractional a, Random a) => a -> a -> IO a
randomOS x y = (* y) . subtract x <$> randomIO

trainableRBM s p0 = makeFoldAll (Minibatch (10, L2 (1.0, 0.000))) (rbmcd1 s) p0

testDigits :: IO ()
testDigits = do
    let importImage l i = fmap normalizeInput <$> unsafeGetBoth l i
    let inDownloads = (++) "/home/julian/Downloads/"
    let lpath  = inDownloads "train-labels-idx1-ubyte"
    let ipath  = inDownloads "train-images-idx3-ubyte"
    let tlpath = inDownloads "t10k-labels-idx1-ubyte"
    let tipath = inDownloads "t10k-images-idx3-ubyte"

    test  <- importImage tlpath tipath
    train <- importImage lpath  ipath
    trainData <- concat <$> (mapM shuffleM $ replicate 30 train)

    p0 <- randomIO

    let trained = L.fold (uberwachtesFold (Minibatch (10, SGD 1.5)) multilayerPerceptron
                          se' p0) trainData

    emptyLn
    print $ (/ 100) . sum . map (\(x, l) ->
        if ((== antwort l) . antwort . prediction trained $ x)
        then 1
        else 0) $ test

    print $ auswerten (antwort . evalBP2 multilayerPerceptron trained)
        (==) (map (\(x, a) -> (x, antwort a)) test)

    where pred p x i l = (antwort l == x) && ((== x) . antwort . prediction p $ i)
          prediction p i = evalBP2 multilayerPerceptron p i

testAND :: IO ()
testAND = do
    -- mische 100000 Samples
    daten <- shuffleM . take 100000 . cycle $ samps
    -- zufälliger Anfangswert, der resultierende Vektor von randomIO
    -- ist standard normal verteilt
    p0 <- (* 0.1) <$> randomIO
    -- Netz wird trainiert und optimiert mit SGD und der Lernrate 0.1
    let trainiert = L.fold (uberwachtesFold (SGD 0.1) linear2 se p0) daten
    -- Hilfsfunktion für Ausgabe
    let disp = printLogicOpWith (evalBP2 linear2 trainiert)

    emptyLn
    disp 0 0
    disp 1 0
    disp 0 1
    disp 1 1

    print $ auswerten (evalBP2 linear2 trainiert) (\x y -> round x == round y) daten

-- Die vier (Eingabe, Ziel) Paare, die gelernt werden
-- diese entsprechen dem logischen UND.
samps :: [(H.R 2, Double)]
samps = [(H.vec2 0 0, 0), (H.vec2 1 0, 0), (H.vec2 0 1, 0), (H.vec2 1 1, 1)]

printLogicOpWith :: (H.R 2 -> Double) -> Double -> Double -> IO ()
printLogicOpWith f a b = print $ disp (show . f $ H.vec2 a b)
    where disp ans = "[" ++ show a ++ ", " ++ show b  ++ "]: " ++ ans

auswerten :: Ord b
          => (a -> b)
          -> (b -> b -> Bool)
          -> [(a, b)]
          -> [Int]
auswerten f pred xs = map (foldr sumCorrect 0) . groupBy sndEq . sortBy comp $ xs
    where sndEq (_, b1) (_, b2) = b1 == b2
          comp  (_, b1) (_, b2) = compare b1 b2
          sumCorrect x n        = if (eval x) then n+1 else n
          eval (a, b)           = pred (f a) b

max_stars_repo_name: QQMBR/BLL-FP-NN

max_issues_repo_name: QQMBR/BLL-FP-NN

max_forks_repo_name: QQMBR/BLL-FP-NN


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE Strict #-}
{-# LANGUAGE FlexibleContexts #-}

module CPVO.IO.Reader.Ecalj.DOS (
  readPDOS
  )
  where

import CPVO.Numeric
import CPVO.IO
--
import Text.Printf as TP
import Numeric.LinearAlgebra
import qualified System.Process as SP
import System.IO (openTempFile,hClose)
-- ===============================================

getPDOS' :: Matrix Double -> String -> [Int] -> [String] -> IO (Matrix Double)
getPDOS' res _ _ []  = return res
getPDOS' res tmpf intAOs (nf:nfiles)  = do
  _ <- SP.system $ "mkdir -p temp; more +2 " ++ nf ++ " > " ++ tmpf
  aPDOS' <- fmap (\x -> sumRow $ (¿) x intAOs) $ loadMatrix tmpf
  getPDOS' (fromBlocks [[res,asColumn aPDOS']]) tmpf intAOs nfiles

-------------------------------------------------------------
-- Input Processing: read PDOS data
--
readPDOS :: Double -> String -> String -> [(Int, ([Char], (String, [Int])))]
         -> IO [(Matrix Double, (Int, (String, (Int, String))))]
readPDOS invStat tailer dir ctrlAtAOs = do
  putStrLn "========readPDOS"
  putStrLn $ show ctrlAtAOs
  sequence
    $ (\x ->  [f a | f <- (readPDOS' invStat dir tailer), a <- x]) ctrlAtAOs

readPDOS' :: Double -> String -> String -> [(Int, (String, (String, [Int])))
         -> IO (Matrix Double, (Int, (String, (Int, String))))]
readPDOS' invStat foldernya tailer = fmap (readOnePDOS foldernya tailer) $ flipBy invStat [1,2] -- spin 1 up n spin 2 down

--readPDOS :: String -> String -> Int -> ((String, String, [Int]), [(Int, String)]) -> IO (Int,String, Matrix Double)
--readPDOS ::(spin, (noAt, (symAt, (labelAt, PDOS :: Matrix Double))))
readOnePDOS :: String
            -> String -> Int -> (Int, (String, (String, [Int])))
            -> IO (Matrix Double, (Int, (String, (Int, String))))
readOnePDOS theFolder tailing spin (noAt,(symAt,(labelAt,intAOs))) = do
  let namaFao = theFolder ++ "/dos.isp" ++ show spin ++ ".site" ++ (TP.printf "%03d" noAt) ++ "." ++ tailing
  (tmpfile,h) <- openTempFile "temp" "aEDOS.suffix"
  hClose h
  _ <- SP.system $ "mkdir -p temp; sed '{1d}' " ++ namaFao ++ " > " ++ tmpfile
  aoE <- fmap (\x -> (¿) x [0]) $ loadMatrix $ tmpfile -- 0th column, Energy column
  let zeroE = asColumn $ konst 0 (rows aoE)                                      -- zero valued column instead of real column
  aPDOS <- fmap (dropColumns 1) $ getPDOS' zeroE tmpfile intAOs [namaFao]                                -- create sum of per atomic AOs (PDOS/atom)
  return $ (fromBlocks [[aoE, aPDOS]] , (spin, (hashSpaceText labelAt, (noAt, symAt))))
------------------------------------------------------------------

-- getPDOS :: String -> String -> Int -> ((String, String, [Int]), [(Int, String)]) -> IO (Int,String, Matrix Double)
-- getPDOS theFolder tailing spin (a@(namaAtom,jdAtom,intAOs),lsAtoms) = do
--   let namaFaos = map (\(x,_) -> theFolder ++ "/dos.isp" ++ show spin ++ ".site" ++ (TP.printf "%03d" x) ++ "." ++ tailing) lsAtoms
--   (tmpfile,h) <- openTempFile "temp" "aEDOS.suffix"
--   hClose h
--   _ <- inshell2text $ "mkdir -p temp; more +2 " ++ (head namaFaos) ++ " > " ++ tmpfile -- this is needed only to generate aoE, the real processing is in getPDOS'
--   aoE <- fmap (\x -> (¿) x [0]) $ loadMatrix tmpfile                             -- 0th column, Energy column
--   let zeroE = asColumn $ konst 0 (rows aoE)                                      -- zero valued column instead of real column
--   aPDOS <- fmap (dropColumns 1) $ getPDOS' zeroE tmpfile intAOs namaFaos                                -- create sum of per atomic AOs (PDOS/atom)
--   return $ (spin, hashSpaceText jdAtom, fromBlocks [[aoE, aPDOS]])
-- ------------------------------------------------------------------
--

max_stars_repo_name: hasanalrasyid/cpvoh

max_issues_repo_name: hasanalrasyid/cpvoh

max_forks_repo_name: hasanalrasyid/cpvoh


text:
module Mobius where

import Data.Complex

moebius :: (Floating a, Fractional a, Eq a) => a -> a -> a -> a -> a -> a
moebius a b c d z
    | a*d - b*c == 0 = undefined
    | otherwise = (a*z + b) / (c*z + d)

lambdamoebius :: (RealFloat a, Fractional a, Eq a) =>
    Complex a -> Complex a -> Complex a
lambdamoebius a z = (z - a) / (1 - ((conjugate a) * z))

circular :: (RealFloat a, Fractional a) => Complex a -> Complex a
circular z = möbius 0 (-1) 1 0 z

elliptic :: (Floating a, Fractional a, Eq a) => a -> a -> a
elliptic theta z = möbius (cos a) (-sin a) (sin a) (cos a) z --(cis theta) * z
    where a = theta

parabolic :: (Floating a, Fractional a) => a -> a -> a
parabolic a z = z + a -- translation

hyperbolic :: (Floating a, Fractional a) => a -> a -> a
hyperbolic theta z = (exp theta) * z

loxodromic :: (Floating a, Fractional a) => a -> a -> a
loxodromic a z = a * z -- S-shaped path

projection :: (Floating a, Fractional a) => a -> a
projection z = tanh

max_stars_repo_name: Vetii/Fractal

max_issues_repo_name: Vetii/Fractal

max_forks_repo_name: Vetii/Fractal


text:
{- |
Module      :  Numeric.GSL.Root
Copyright   :  (c) Alberto Ruiz 2009
License     :  GPL

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional
Portability :  uses ffi

Multidimensional root finding.

<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>

The example in the GSL manual:

@import Numeric.GSL
import Numeric.LinearAlgebra(format)
import Text.Printf(printf)

rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]

disp = putStrLn . format \"  \" (printf \"%.3f\")

main = do
    let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
    print sol
    disp path

\> main
[1.0,1.0]
 0.000  -10.000  -5.000  11.000  -1050.000
 1.000   -3.976  24.827   4.976     90.203
 2.000   -3.976  24.827   4.976     90.203
 3.000   -3.976  24.827   4.976     90.203
 4.000   -1.274  -5.680   2.274    -73.018
 5.000   -1.274  -5.680   2.274    -73.018
 6.000    0.249   0.298   0.751      2.359
 7.000    0.249   0.298   0.751      2.359
 8.000    1.000   0.878  -0.000     -1.218
 9.000    1.000   0.989  -0.000     -0.108
10.000    1.000   1.000   0.000      0.000
@

-}
-----------------------------------------------------------------------------

module Numeric.GSL.Root (
    uniRoot, UniRootMethod(..),
    uniRootJ, UniRootMethodJ(..),
    root, RootMethod(..),
    rootJ, RootMethodJ(..),
) where

import Data.Packed.Internal
import Data.Packed.Matrix
import Numeric.GSL.Internal
import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
import Foreign.C.Types
import System.IO.Unsafe(unsafePerformIO)

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

data UniRootMethod = Bisection
                   | FalsePos
                   | Brent
                   deriving (Enum, Eq, Show, Bounded)

uniRoot :: UniRootMethod
        -> Double
        -> Int
        -> (Double -> Double)
        -> Double
        -> Double
        -> (Double, Matrix Double)
uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit

uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
    fp <- mkDoublefun f
    rawpath <- createMIO maxit 4
                         (c_root m fp epsrel (fi maxit) xl xu)
                         "root"
    let it = round (rawpath @@> (maxit-1,0))
        path = takeRows it rawpath
        [sol] = toLists $ dropRows (it-1) path
    freeHaskellFunPtr fp
    return (sol !! 1, path)

foreign import ccall safe "root"
    c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM

-------------------------------------------------------------------------
data UniRootMethodJ = UNewton
                    | Secant
                    | Steffenson
                    deriving (Enum, Eq, Show, Bounded)

uniRootJ :: UniRootMethodJ
        -> Double
        -> Int
        -> (Double -> Double)
        -> (Double -> Double)
        -> Double
        -> (Double, Matrix Double)
uniRootJ method epsrel maxit fun dfun x = uniRootJGen (fi (fromEnum method)) fun
    dfun x epsrel maxit

uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
    fp <- mkDoublefun f
    dfp <- mkDoublefun df
    rawpath <- createMIO maxit 2
                         (c_rootj m fp dfp epsrel (fi maxit) x)
                         "rootj"
    let it = round (rawpath @@> (maxit-1,0))
        path = takeRows it rawpath
        [sol] = toLists $ dropRows (it-1) path
    freeHaskellFunPtr fp
    return (sol !! 1, path)

foreign import ccall safe "rootj"
    c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double)
            -> Double -> CInt -> Double -> TM

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

data RootMethod = Hybrids
                | Hybrid
                | DNewton
                | Broyden
                deriving (Enum,Eq,Show,Bounded)

-- | Nonlinear multidimensional root finding using algorithms that do not require
-- any derivative information to be supplied by the user.
-- Any derivatives needed are approximated by finite differences.
root :: RootMethod
     -> Double                     -- ^ maximum residual
     -> Int                        -- ^ maximum number of iterations allowed
     -> ([Double] -> [Double])     -- ^ function to minimize
     -> [Double]                   -- ^ starting point
     -> ([Double], Matrix Double)  -- ^ solution vector and optimization path

root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit

rootGen m f xi epsabs maxit = unsafePerformIO $ do
    let xiv = fromList xi
        n   = dim xiv
    fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
    rawpath <- vec xiv $ \xiv' ->
                   createMIO maxit (2*n+1)
                         (c_multiroot m fp epsabs (fi maxit) // xiv')
                         "multiroot"
    let it = round (rawpath @@> (maxit-1,0))
        path = takeRows it rawpath
        [sol] = toLists $ dropRows (it-1) path
    freeHaskellFunPtr fp
    return (take n $ drop 1 sol, path)

foreign import ccall safe "multiroot"
    c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM

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

data RootMethodJ = HybridsJ
                 | HybridJ
                 | Newton
                 | GNewton
                deriving (Enum,Eq,Show,Bounded)

-- | Nonlinear multidimensional root finding using both the function and its derivatives.
rootJ :: RootMethodJ
      -> Double                     -- ^ maximum residual
      -> Int                        -- ^ maximum number of iterations allowed
      -> ([Double] -> [Double])     -- ^ function to minimize
      -> ([Double] -> [[Double]])   -- ^ Jacobian
      -> [Double]                   -- ^ starting point
      -> ([Double], Matrix Double)  -- ^ solution vector and optimization path

rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun jac xinit epsabs maxit

rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
    let xiv = fromList xi
        n   = dim xiv
    fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
    jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
    rawpath <- vec xiv $ \xiv' ->
                   createMIO maxit (2*n+1)
                         (c_multirootj m fp jp epsabs (fi maxit) // xiv')
                         "multiroot"
    let it = round (rawpath @@> (maxit-1,0))
        path = takeRows it rawpath
        [sol] = toLists $ dropRows (it-1) path
    freeHaskellFunPtr fp
    freeHaskellFunPtr jp
    return (take n $ drop 1 sol, path)

foreign import ccall safe "multirootj"
    c_multirootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM

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

checkdim1 n v
    | dim v == n = v
    | otherwise = error $ "Error: "++ show n
                        ++ " components expected in the result of the function supplied to root"

checkdim2 n m
    | rows m == n && cols m == n = m
    | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
                        ++ " Jacobian expected in rootJ"


text:
{-# OPTIONS_GHC -Wall                     #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing  #-}
{-# OPTIONS_GHC -fno-warn-type-defaults   #-}
{-# OPTIONS_GHC -fno-warn-unused-do-bind  #-}
{-# OPTIONS_GHC -fno-warn-missing-methods #-}
{-# OPTIONS_GHC -fno-warn-orphans         #-}

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

-- | An extended Kalman filter. Note that this could be generalized
-- further. If you need a Kalman filter and have a state update matrix
-- \(A\) then you can just use @const bigA@ as the function which
-- returns the Jacobian.

module Kalman (
  extKalman
  )where

import GHC.TypeLits
import Numeric.LinearAlgebra.Static hiding ( create )
import Data.Maybe ( fromJust )

-- | The model for the extended Kalman filter is given by
--
-- \[
-- \begin{aligned}
-- \boldsymbol{x}_i &= \boldsymbol{a}(\boldsymbol{x}_{i-1}) + \boldsymbol{\psi}_{i-1} \\
-- \boldsymbol{y}_i &= {H}\boldsymbol{x}_i + \boldsymbol{\upsilon}_i
-- \end{aligned}
-- \]
--
-- where
--
-- * \(\boldsymbol{a}\)\ is some non-linear vector-valued state update
-- function.
--
-- * \(\boldsymbol{\psi}_{i}\) are independent identically normally
-- distributed random variables with mean 0 representing the fact that
-- the state update is noisy: \(\boldsymbol{\psi}_{i} \sim {\cal{N}}(0,Q)\).
--
-- * \({H}\) is a matrix which represents how we observe the
-- hidden state process.
--
-- * \(\boldsymbol{\upsilon}_i\) are independent identically normally
-- distributed random variables with mean 0 represent the fact that
-- the observations are noisy: \(\boldsymbol{\upsilon}_{i} \sim {\cal{N}}(0,R)\).
--
-- We assume the whole process starts at 0 with our belief of the state (aka the
-- prior state) being given by
-- \(\boldsymbol{x}_0 \sim {\cal{N}}(\boldsymbol{\mu}_0, \Sigma_0)\)
extKalman ::  forall m n .
              (KnownNat m, KnownNat n, (1 <=? n) ~ 'True, (1 <=? m) ~ 'True) =>
              R n              -- ^ Prior mean \(\boldsymbol{\mu}_0\)
              -> Sq n          -- ^ Prior variance \(\Sigma_0\)
              -> L m n         -- ^ Observation map \(H\)
              -> Sq m          -- ^ Observation noise \(R\)
              -> (R n -> R n)  -- ^ State update function \(\boldsymbol{a}\)
              -> (R n -> Sq n) -- ^ A function which returns the
                               -- Jacobian of the state update
                               -- function at a given point
                               -- \(\frac{\partial \boldsymbol{a}}{\partial \boldsymbol{x}}\)
              -> Sq n          -- ^ State noise \(Q\)
              -> [R m]         -- ^ List of observations \(\boldsymbol{y}_i\)
              -> [(R n, Sq n)] -- ^ List of posterior means and variances \((\hat{\boldsymbol{x}}_i, \hat{\Sigma}_i)\)
extKalman muPrior sigmaPrior bigH bigSigmaY
  littleA bigABuilder bigSigmaX ys = result
  where
    result = scanl update (muPrior, sigmaPrior) ys

    update :: (R n, Sq n) -> R m -> (R n, Sq n)
    update (xHatFlat, bigSigmaHatFlat) y =
      (xHatFlatNew, bigSigmaHatFlatNew)
      where

        v :: R m
        v = y - (bigH #> xHatFlat)

        bigS :: Sq m
        bigS = bigH <> bigSigmaHatFlat <> (tr bigH) + bigSigmaY

        bigK :: L n m
        bigK = bigSigmaHatFlat <> (tr bigH) <> (inv bigS)

        xHat :: R n
        xHat = xHatFlat + bigK #> v

        bigSigmaHat :: Sq n
        bigSigmaHat = bigSigmaHatFlat - bigK <> bigS <> (tr bigK)

        bigA :: Sq n
        bigA = bigABuilder xHat

        xHatFlatNew :: R n
        xHatFlatNew = littleA xHat

        bigSigmaHatFlatNew :: Sq n
        bigSigmaHatFlatNew = bigA <> bigSigmaHat <> (tr bigA) + bigSigmaX

inv :: (KnownNat n, (1 <=? n) ~ 'True) => Sq n -> Sq n
inv m = fromJust $ linSolve m eye

max_stars_repo_name: idontgetoutmuch/Kalman

max_issues_repo_name: f-o-a-m/Kalman

max_forks_repo_name: f-o-a-m/Kalman