text:
{-# LANGUAGE FlexibleContexts #-}
module Features where

import Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra.HMatrix (eigSH')
import Codec.Compression.GZip (decompress)
import qualified Data.ByteString.Lazy as BS
import DigitRecognition
import Data.Bifunctor (Bifunctor(bimap))

mean :: Matrix R -> Vector R
mean x = tr x #> konst k r
    where
        r = rows x
        k = 1 / fromIntegral r

centreData :: Matrix R -> Matrix R
centreData x = x - fromRows (replicate (rows x) (mean x))

principalComponents :: Matrix R -> Matrix R
principalComponents x = vs
    where
        x' = centreData x
        cov = tr x' LA.<> x'
        (_, vs) = eigSH' cov

projectOntoPC :: Matrix R -> Matrix R -> Int -> Matrix R
projectOntoPC x pcs n = cmap (\x -> if x < 1.0e-16 && x > -1.0e-16 then 0 else x) $ centreData x LA.<> (pcs ¿ [0 .. n - 1])

triUpper :: Int -> Int -> [(Int, Int)]
triUpper m offset = [(a,b) | a <- [0.. m - 1], b <- [offset .. m - 1], a + offset <= b]

twoTerms :: [R] -> [R]
twoTerms xs = concat $ twoTerms' (drop 1 xs) xs
    where
        twoTerms' [] _ = []
        twoTerms' xs ys = twoTerms'' (zip xs ys) ++ twoTerms' (drop 1 xs) ys
        twoTerms'' xs = map (\(x, y) -> [sqrt 6 * x * y]) xs

threeTerms :: Int -> Int -> Int -> [(Int, Int)]
threeTerms m i offset = threeTerms' m i offset
    where
        threeTerms' 1 _ _ = []
        threeTerms' m i offset = fmap (bimap (offset +) (+i)) (triUpper m 1) ++ threeTerms' (m - 1) (i+1) (offset + m)

squareTerms :: [R] -> [R]
squareTerms xs = concat $ squareTerms' (drop 1 xs) xs
    where
        squareTerms' [] _ = []
        squareTerms' xs ys = squareTerms'' (zip xs ys) ++ squareTerms' (drop 1 xs) ys
        squareTerms'' xs = map (\(x, y) -> [sqrt 3 * x ** 2 * y, sqrt 3 * x * y ** 2]) xs

crossTerms :: Vector R -> [R]
crossTerms v = coords
    where
        symm = v `outer` v
        pairs = vjoin $ zipWith pairs' (toRows symm) [1..]
        triples = pairs `outer` v
        pairs' vs i = subVector i (size vs - i) vs
        coords = map ((sqrt 6 * ) . (triples `atIndex`)) $ threeTerms (cols triples - 1) 1 0

cubicFeatures :: Matrix R -> Matrix R
cubicFeatures x = cubeTerms ||| (n >< squareCols) squareTerms' ||| (n >< linearCols) linearTerms' ||| (n >< crossCols) crossTerms'
    where
        n = rows x
        x' = x ||| (rows x >< 1) (replicate n 1.0)
        cubeTerms = cmap (**3) x'
        squareTerms' = concatMap (squareTerms . toList) (toRows x')
        linearTerms' = concatMap (twoTerms . toList) (toRows x)
        crossTerms' = concatMap crossTerms (toRows x)
        squareCols = length squareTerms' `div` n
        linearCols = length linearTerms' `div` n
        crossCols = length crossTerms' `div` n

main :: IO ()
main = do

    tr <- decompress <$> BS.readFile "train-images-idx3-ubyte.gz"
    trl <- decompress <$> BS.readFile "train-labels-idx1-ubyte.gz"
    te <- decompress <$> BS.readFile "t10k-images-idx3-ubyte.gz"
    tel <- decompress <$> BS.readFile "t10k-labels-idx1-ubyte.gz"

    let trainData = buildData 60000 tr
    let trainLabels = buildLabels trl
    let testData = buildData 10000 te
    let testLabels = buildLabels tel
    let temp = 1
    let alpha = 0.3
    let lambda = 0.0001
    let k = 10
    let iter = 150

--    let nComponents = 18
    let pcs = principalComponents trainData
--    let trainData' = projectOntoPC trainData pcs nComponents
--    let testData' = projectOntoPC testData pcs nComponents

--    let (theta, cs) = softmaxRegression trainData' trainLabels temp alpha lambda k iter
--    print $ computeTestError testData' testLabels theta temp

    let nCubic = 10
    let trainData'' = projectOntoPC trainData pcs nCubic
    let testData'' = projectOntoPC testData pcs nCubic

    let trainCube = cubicFeatures trainData''
    let testCube = cubicFeatures testData''

    let (theta', _) = softmaxRegression trainCube trainLabels temp alpha lambda k iter
    print $ computeTestError testCube testLabels theta' temp


text:
module HashedExpression.Utils where

import Data.Array
import Data.Complex (Complex (..))
import qualified Data.Complex as Complex
import Data.Function ((&))
import qualified Data.IntMap.Strict as IM
import Data.List (foldl')
import Data.List.Split (splitOn)
import Data.Map (Map, fromList)
import qualified Data.Map as Map
import Data.Maybe
import qualified Data.Set as Set
import qualified Data.Text as T
import Data.Time (diffUTCTime, getCurrentTime)
import GHC.IO.Unsafe (unsafePerformIO)
import GHC.Stack (HasCallStack)
import HashedExpression.Internal.Base
import HashedExpression.Internal.Node
import Prelude hiding ((^))
import qualified Prelude

-- | Forward pipe operator
(|>) :: a -> (a -> b) -> b
(|>) = (&)

infixl 1 |>

-- |
measureTime :: IO a -> IO ()
measureTime action = do
  beforeTime <- getCurrentTime
  action
  afterTime <- getCurrentTime
  putStrLn $ "Took " ++ show (diffUTCTime afterTime beforeTime) ++ " seconds"

-- | Check if all elements of the list is equal
allEqual :: (Eq a) => [a] -> Bool
allEqual [] = True
allEqual [x] = True
allEqual (x : y : xs) = x == y && allEqual (y : xs)

fromR :: Double -> Complex Double
fromR x = x :+ 0

-- |
pullConstant :: ExpressionMap -> NodeID -> Maybe (Shape, Double)
pullConstant mp n
  | (shape, R, Const c) <- retrieveNode n mp = Just (shape, c)
  | otherwise = Nothing

-- |
pullConstants :: ExpressionMap -> [NodeID] -> Maybe (Shape, [Double])
pullConstants mp ns
  | xs@(x : _) <- mapMaybe (pullConstant mp) ns = Just (fst x, map snd xs)
  | otherwise = Nothing

-- |
isZero :: ExpressionMap -> NodeID -> Bool
isZero mp nId
  | Const 0 <- retrieveOp nId mp = True
  | RealImag arg1 arg2 <- retrieveOp nId mp,
    Const 0 <- retrieveOp arg1 mp,
    Const 0 <- retrieveOp arg2 mp =
    True
  | otherwise = False

-- |
isOne :: ExpressionMap -> NodeID -> Bool
isOne mp nId
  | Const 1 <- retrieveOp nId mp = True
  | RealImag arg1 arg2 <- retrieveOp nId mp,
    Const 1 <- retrieveOp arg1 mp,
    Const 0 <- retrieveOp arg2 mp =
    True
  | otherwise = False

-- |
isConstant :: ExpressionMap -> NodeID -> Bool
isConstant mp nId
  | Const _ <- retrieveOp nId mp = True
  | otherwise = False

-- |
pullSumOperands :: ExpressionMap -> NodeID -> [NodeID]
pullSumOperands mp nId
  | Sum operands <- retrieveOp nId mp = operands
  | otherwise = [nId]

-- |
pullProdOperands :: ExpressionMap -> NodeID -> [NodeID]
pullProdOperands mp nId
  | Mul operands <- retrieveOp nId mp = operands
  | otherwise = [nId]

toRange :: DimSelector -> Int -> (Int, Int, Int)
toRange (At i) size = (i, i, 1)
toRange (Range start end step) size = (start, if end >= start then end else end + size, step)

mkIndices :: DimSelector -> Int -> [Int]
mkIndices (At i) size = [i]
mkIndices (Range start end step) size =
  map (`mod` size) [start, start + step .. if end >= start then end else end + size]

breakSub :: Int -> [a] -> [([a], [a], [a])]
breakSub len [] = []
breakSub len xs | len > length xs = []
breakSub len xs@(y : ys) =
  ([], take len xs, drop len xs) :
  map (\(prefix, sub, suffix) -> (y : prefix, sub, suffix)) (breakSub len ys)

coercible :: Shape -> Shape -> Bool
coercible sml target =
  any
    ( \(prefix, sub, suffix) ->
        sub == sml && all (== 1) prefix && all (== 1) suffix
    )
    $ breakSub (length sml) target

-------------------------------------------------------------------------------
zipMp :: forall a b c. Ord a => Map a b -> Map a c -> Map a (b, c)
zipMp mp1 mp2 = foldl' f Map.empty $ Map.keys mp1
  where
    f acc k = case (Map.lookup k mp1, Map.lookup k mp2) of
      (Just v1, Just v2) -> Map.insert k (v1, v2) acc
      _ -> acc

-------------------------------------------------------------------------------
zipMp3 :: forall a b c d. Ord a => Map a b -> Map a c -> Map a d -> Map a (b, c, d)
zipMp3 mp1 mp2 mp3 = foldl' f Map.empty $ Map.keys mp1
  where
    f acc k = case (Map.lookup k mp1, Map.lookup k mp2, Map.lookup k mp3) of
      (Just v1, Just v2, Just v3) -> Map.insert k (v1, v2, v3) acc
      _ -> acc

instance PowerOp Double Int where
  (^) x y = x Prelude.^ y

instance ComplexRealOp Double (Complex Double) where
  (+:) = (:+)
  xRe = Complex.realPart
  xIm = Complex.imagPart
  conjugate = Complex.conjugate

-------------------------------------------------------------------------------
class ToText a where
  tt :: a -> T.Text
  ttq :: a -> T.Text
  ttq s = "\"" <> tt s <> "\""

instance ToText Double where
  tt s = T.pack $ show s

instance ToText Int where
  tt s = T.pack $ show s

instance ToText String where
  tt = T.pack

instance ToText T.Text where
  tt = id

ext: hs

max_stars_repo_name: McMasterU/HashedExpression

max_issues_repo_name: McMasterU/HashedExpression

max_forks_repo_name: McMasterU/HashedExpression


text:
{-# LANGUAGE ScopedTypeVariables #-}

module Jakway.Blackjack.Random.Checks
(
testDeckEvenDistribution,
testRNGDistribution
)
where

import qualified Statistics.Test.Types as Stats
import qualified Statistics.Test.ChiSquared as Stats
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as UM
import System.Random
import Jakway.Blackjack.AI
import Jakway.Blackjack.CardOps
import Jakway.Blackjack.Game
import Data.Hashable
import Control.Monad (foldM)
import Data.List (genericReplicate, genericTake)

deck_value_range :: Int
deck_value_range = 104 -- ^ 52 cards plus the visibility flag

-- |record a new observation, updating the U.Vector in place
deckToObservation :: U.Vector Int -> Deck -> U.Vector Int
deckToObservation vec deck = U.modify (\v -> UM.write v index newCount) vec
        where vecRange = U.length vec
              index = (`mod` vecRange) . hash $ (deck !! 13) -- ^ we've arbitrarily chosen the position to take a card from
              newCount = (vec U.! index) + 1

vecIncrement :: U.Vector Int -> Int -> U.Vector Int
vecIncrement v pos = U.modify (\mv -> UM.write mv pos newCount) v
        where newCount = (v U.! pos) + 1

newDeckIO :: IO Deck
newDeckIO = getStdGen >>= return . infiniteShuffledDeck

observeDeck :: Integer -> AI -> [AI] -> IO (U.Vector Int)
observeDeck numSamples dealerAI playerAIs =
        foldM (\vec _ -> newDeckIO >>= (\d -> (return $ evalGameKeepDeck dealerAI playerAIs d) >>= (\maybeDeck ->
                                case maybeDeck of Nothing -> return vec
                                                  -- | modify the list of observations and return it
                                                  (Just (resDeck,_)) -> return $ deckToObservation vec resDeck))) samplesVec [1..numSamples]

        where samplesVec = (U.fromList (replicate deck_value_range 0)) :: U.Vector Int

testDeckEvenDistribution :: Double -> Integer -> AI -> [AI] -> IO Stats.TestResult
testDeckEvenDistribution pvalue numSamples dealerAI playerAIs = do
        observations <- observeDeck numSamples dealerAI playerAIs
        return $ Stats.chi2test pvalue additionalDF (evenDistribution observations)
        where additionalDF = 0

testRNGDistribution :: Double -> Integer -> Integer -> Integer -> AI -> [AI] -> IO Stats.TestResult
testRNGDistribution pvalue numSamples numRNGSamples rngMaxRange dealerAI playerAIs = getStdGen >>= (\gen -> observeDeck numSamples dealerAI playerAIs >>= \samples -> return $ Stats.chi2test pvalue 0 $ rngDistribution numRNGSamples rngMaxRange gen samples)

--instead of testing against an even distribution, test against the
--distribution of the default StdGen

rngDistribution :: (RandomGen g) => Integer -> Integer -> g -> U.Vector Int -> U.Vector (Int, Double)
rngDistribution numRngObservations maxr gen observed = U.zip observed percents
        where -- how many observations to get the standard RNG
              -- distribution?
              startingVec = (U.fromList (genericReplicate maxr 0)) :: U.Vector Int
              rngObservations = genericTake numRngObservations $ randomRs (0, maxr) gen -- ^ lower bound on random numbers is 0
              -- |increment the count for this observation
              rngObservationCounts = foldr (\thisObservation vec -> vecIncrement vec thisObservation) startingVec $ map fromIntegral rngObservations
              -- |the distribution as a percent for each bin (i.e. a list
              -- of doubles)
              percents :: U.Vector Double
              percents = U.map (\o -> (fromIntegral o) / (fromIntegral numRngObservations)) rngObservationCounts

-- |Takes as input a set of bins with the number of observations in each
-- returns the set of bins you would expect given an even distribution
-- across all bins zipped with the observed values
evenDistribution :: U.Vector Int -> U.Vector (Int, Double)
evenDistribution observed = U.zip observed expectedObservations
        where numBins = U.length observed
              totalObservations = fromIntegral . U.sum $ observed
              expectedObservations = U.fromList $ replicate numBins (totalObservations / (fromIntegral numBins))

max_stars_repo_name: tjakway/haskell-blackjack

max_issues_repo_name: tjakway/haskell-blackjack

max_forks_repo_name: tjakway/haskell-blackjack


text:
module Main where

import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning.Types as T
import qualified MachineLearning as ML
import qualified MachineLearning.Classification.OneVsAll as OVA
import qualified MachineLearning.PCA as PCA

processFeatures :: T.Matrix -> T.Matrix
processFeatures = ML.addBiasDimension . (ML.mapFeatures 2)

calcAccuracy :: T.Matrix -> T.Vector -> [T.Vector] -> Double
calcAccuracy x y thetas = OVA.calcAccuracy y yPredicted
  where yPredicted = OVA.predict x thetas

main = do
  -- Step 1. Data loading.
  -- Step 1.1 Training Data loading.
  (x, y) <- pure ML.splitToXY <*> LA.loadMatrix "digits_classification/optdigits.tra"
  -- Step 1.1 Testing Data loading.
  (xTest, yTest) <- pure ML.splitToXY <*> LA.loadMatrix "digits_classification/optdigits.tes"

  -- Step 2. Outputs and features preprocessing.
  let numLabels = 10
      x' = processFeatures x
  -- Step 3. Dimensionality Reduction using PCA.
      (reduceDims, reducedDimensions, x1) = PCA.getDimReducer_rv x' 0.99999
      initialTheta = LA.konst 0 (LA.cols x1)
      initialThetas = replicate numLabels initialTheta
  -- Step 4. Learning.
      (thetas, optPath) = OVA.learn (OVA.BFGS2 0.1 0.5) 0.0001 30 (OVA.L2 30) numLabels x1 y initialThetas
  -- Step 5. Prediction and checking accuracy
      accuracyTrain = calcAccuracy x1 y thetas
      accuracyTest = calcAccuracy (reduceDims $ processFeatures xTest) yTest thetas

  -- Step 6. Printing results.
  putStrLn "\n=== Logistic Regression with PCA (Digits Recognition) ==="

  putStrLn $ "\nNumber of iterations to learn: " ++ show (map LA.rows optPath)

  putStrLn $ "\nReduced from " ++ show (LA.cols x') ++ " to " ++ show (LA.cols x1) ++ " dimensions"

  putStrLn $ "\nAccuracy on train set (%): " ++ show (accuracyTrain*100)
  putStrLn $ "Accuracy on test set (%): " ++ show (accuracyTest*100)

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Scene.Material.Diffuse where

import Numeric.LinearAlgebra

import Scene.Base
import Math.Color
import Math.Ray
import Scene.Monad

import Control.Applicative
import Control.Monad.Random

doWhileM :: Monad m => (a -> Bool) -> m a -> m a
doWhileM cond prog = prog >>= (\x -> if cond x then return x else doWhileM cond prog)

diffuse :: (Color -> Color) -> Material
diffuse absorb (Ray p normal) (Ray q incidence) = do
    s_unit <- doWhileM (\s -> s <.> s < 1) $ vec3 <$> getRandomR (-1,1) <*> getRandomR (-1,1) <*> getRandomR (-1,1)
    let s = scale (norm_2 normal) s_unit + normal
    ahead <- trace (Ray p s)
    return $! absorb ahead

max_stars_repo_name: danrocag/raytracer

max_issues_repo_name: danrocag/raytracer

max_forks_repo_name: danrocag/raytracer


text:
import Data.Bool
import Data.Complex (Complex((:+)), magnitude)

mandelbrot
  :: RealFloat a
  => Complex a -> Complex a
mandelbrot a = iterate ((a +) . (^ 2)) 0 !! 50

main :: IO ()
main =
  mapM_
    putStrLn
    [ [ bool ' ' '*' (2 > magnitude (mandelbrot (x :+ y)))
      | x <- [-2,-1.9685 .. 0.5] ]
    | y <- [1,0.95 .. -1] ]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
-- |
-- Module      :  DirichletSpec
-- Description :  Unit tests for DirichletSpec
-- Copyright   :  (c) Dominik Schrempf, 2021
-- License     :  GPL-3.0-or-later
--
-- Maintainer  :  dominik.schrempf@gmail.com
-- Stability   :  unstable
-- Portability :  portable
--
-- Creation date: Tue Oct 20 10:40:06 2020.
module Statistics.Distribution.DirichletSpec
  ( spec,
  )
where

import Control.Monad
import Data.Either
import qualified Data.Vector.Unboxed as V
import Numeric.Log hiding (sum)
import Statistics.Distribution.Dirichlet
import Test.Hspec
import System.Random.MWC

eps :: Double
eps = 1e-14

dd3 :: DirichletDistribution
dd3 = either error id $ dirichletDistribution $ V.fromList [0.2, 10, 20]

alphas10 :: V.Vector Double
alphas10 = V.fromList [0.1, 0.2, 0.3, 0.4, 0.5, 16, 17, 18, 19, 20]

dd10 :: DirichletDistribution
dd10 = either error id $ dirichletDistribution alphas10

ddSym :: Int -> Double -> DirichletDistribution
ddSym n a = either error id $ dirichletDistribution $ V.replicate n a

-- Extract means.
xbar :: Int -> [V.Vector Double] -> Double
xbar n xss = sum xs / fromIntegral (length xs)
  where xs = map (V.! n) xss

spec :: Spec
spec = do
  describe "dirichletDistribution" $ do
    it "only works for valid parameter vectors" $ do
      dirichletDistribution (V.fromList []) `shouldSatisfy` isLeft
      dirichletDistribution (V.fromList [1]) `shouldSatisfy` isLeft
      dirichletDistribution (V.fromList [-3, 3, 192]) `shouldSatisfy` isLeft
      dirichletDistribution (V.fromList [0, 3, 192]) `shouldSatisfy` isLeft
  describe "dirichletDensity" $ do
    it "is correct for some test cases" $ do
      let ddSym2 = ddSym 2 0.5
      let r = dirichletDensity ddSym2 (V.fromList [0.2, 0.8])
      abs (exp (ln r) - 0.7957747154594766) `shouldSatisfy` (< eps)
      let rBound = dirichletDensity ddSym2 (V.fromList [0, 1])
      rBound `shouldBe` 0
      let rOutBound1 = dirichletDensity ddSym2 (V.fromList [0, 1.1])
      rOutBound1 `shouldBe` 0
      let rOutBound2 = dirichletDensity ddSym2 (V.fromList [-0.1, 0.9])
      rOutBound2 `shouldBe` 0
      let r3 = dirichletDensity dd3 (V.fromList [0.3, 0.3, 0.4])
      abs (exp (ln r3) - 0.0001217825570884453) `shouldSatisfy` (< eps)
      let rWrongDim = dirichletDensity dd3 (V.fromList [0.3, 0.7])
      rWrongDim `shouldBe` 0
  describe "dirichletSample" $ do
    it "returns valid value vectors with expected mean" $ do
      g <- create
      let ddSym10 = ddSym 10 10
      xs <- replicateM 1000 (dirichletSample ddSym10 g)
      map V.length xs `shouldBe` replicate 1000 10
      -- print [ xbar i xs | i <- [0..9]]
      [ abs (xbar i xs - 0.1) > 0.01 | i <- [0..9]] `shouldBe` replicate 10 False
      xs' <- replicateM 1000 (dirichletSample dd10 g)
      map V.length xs' `shouldBe` replicate 1000 10
      -- print [ xbar i xs | i <- [0..9]]
      let aSum = V.sum alphas10
      [ abs (xbar i xs' - (alphas10 V.! i / aSum)) > 0.01 | i <- [0..9]] `shouldBe` replicate 10 False

max_stars_repo_name: dschrempf/dirichlet

max_issues_repo_name: dschrempf/dirichlet

max_forks_repo_name: dschrempf/dirichlet


text:
module Main (main) where

import Data.Complex
import Data.HDF5
import qualified Data.Vector.Storable as V
import System.Directory
import Test.Hspec
import Prelude hiding (withFile)

anyH5Exception :: Selector H5Exception
anyH5Exception = const True

main :: IO ()
main = hspec $ do
  describe "Files" $ do
    it "opens an existing file" $ do
      doesFileExist "test/reading_test.h5" `shouldReturn` True
      _ <- withFile "test/reading_test.h5" ReadOnly (\_ -> return ())
      return ()
    it "throws an exceptions when file does not exist" $ do
      withFile "test/not_there.h5" ReadOnly (\_ -> return ()) `shouldThrow` anyH5Exception
    it "truncates a file if it already exists" $ do
      copyFile "test/reading_test.h5" "test/truncate_test.h5"
      _ <- withFile "test/truncate_test.h5" WriteTruncate (\_ -> return ())
      isSmaller <- (<) <$> getFileSize "test/truncate_test.h5" <*> getFileSize "test/reading_test.h5"
      isSmaller `shouldBe` True
      removeFile "test/truncate_test.h5"
    it "appends to a file if it already exists" $ do
      copyFile "test/reading_test.h5" "test/append_test.h5"
      _ <- withFile "test/append_test.h5" WriteAppend (\_ -> return ())
      isSame <- (==) <$> getFileSize "test/append_test.h5" <*> getFileSize "test/reading_test.h5"
      isSame `shouldBe` True
      removeFile "test/append_test.h5"
    it "creates a file if it does not exist" $ do
      doesFileExist "test/a_non_existant_file.h5" `shouldReturn` False
      _ <- withFile "test/a_non_existant_file.h5" WriteTruncate (\_ -> return ())
      doesFileExist "test/a_non_existant_file.h5" `shouldReturn` True
      removeFile "test/a_non_existant_file.h5"
  describe "Datasets" $ do
    it "writes strided datasets" $ do
      (matrix :: TemporaryStridedMatrix Float) <-
        withFile' "test/strided_test_file.h5" WriteTruncate $ \g -> do
          writeDataset g "A" $ TemporaryStridedMatrix (2, 3) 4 (V.fromList [(1 :: Float) .. 8])
          openDataset @Text g "A" >>= readDataset
      matrix `shouldBe` (TemporaryStridedMatrix (2, 3) 3 (V.fromList [1, 2, 3, 5, 6, 7]))
      removeFile "test/strided_test_file.h5"
    it "writes complex datasets" $ do
      xs <-
        withFile' "test/complex_test_file.h5" WriteTruncate $ \g -> do
          writeDataset g "A" $ [((1 :: Double) :+ 2), ((-5) :+ 0.1)]
          openDataset @Text g "A" >>= readDataset
      xs `shouldBe` [((1 :: Double) :+ 2), ((-5) :+ 0.1)]
      removeFile "test/complex_test_file.h5"

-- removeFile "test/strided_test_file.h5"


text:
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE OverloadedStrings #-}

module Main where

import Data.List (dropWhileEnd, intersperse)

import System.Random.MWC.Probability (Gen)
import qualified System.Random.MWC.Probability as MWC
import Control.Monad.Trans.State.Strict (execStateT)

import Statistics.Distribution.Normal
import Statistics.Distribution.Uniform

import Conduit
import System.Environment (getArgs)

import Data.Model
import Numeric.MCMC

import Data.Map.Strict (Map)
import qualified Data.Map.Strict as M

import qualified Data.ByteString as BS (readFile)
import Data.Aeson

-- ~from https://hackage.haskell.org/package/declarative-0.2.1/docs/src/Numeric-MCMC.html#mcmc
-- A Markov chain driven by an arbitrary transition operator.
-- now using Conduit instead of Pipes
chain
  :: PrimMonad m
  => Transition m b
  -> b
  -> Gen (PrimState m)
  -> Producer m b
chain transition = loop where
  loop state prng = do
    next <- lift (MWC.sample (execStateT transition state) prng)
    yield next
    loop next prng

takeWhileEnd :: (a -> Bool) -> [a] -> [a]
takeWhileEnd f xs = case span f xs of
                         (ys, [])    -> ys
                         (_, (_:ys)) -> takeWhileEnd f ys

takeEveryC :: Monad m => Int -> Conduit a m a
takeEveryC n = do dropC (n-1)
                  x <- await
                  case x of
                       Just y -> yield y >> takeEveryC n
                       Nothing -> return ()

showList' :: [String] -> String
showList' = mconcat . intersperse ","

regions :: [RegName]
regions = [ "lvJ_1tag_0addtag_unblind_SR"
          , "lvJ_2tag_0addtag_unblind_SR"
          , "lvJ_1tag_1paddtag_unblind_SR"
          , "lvJ_1tag_0addtag_unblind_lowMH"
          , "lvJ_2tag_0addtag_unblind_lowMH"
          , "lvJ_1tag_1paddtag_unblind_lowMH"
          , "lvJ_1tag_0addtag_unblind_highMH"
          , "lvJ_2tag_0addtag_unblind_highMH"
          , "lvJ_1tag_1paddtag_unblind_highMH"
          ]

processes :: [ProcName]
processes  = [ "TTbar"
             , "Wb"
             , "Wc"
             , "Wl"
             , "stopWt"
             , "stopt"
             , "HVTWHlvqq2000"
             ]

shapeNPs :: [String]
shapeNPs = [ "EG_RESOLUTION_ALL__1up"
           , "EG_SCALE_ALL__1up"
           , "EL_EFF_ID_TotalCorrUncertainty__1up"
           , "EL_EFF_Iso_TotalCorrUncertainty__1up"
           , "EL_EFF_Reco_TotalCorrUncertainty__1up"
           , "EL_EFF_Trigger_TotalCorrUncertainty__1up"
           -- , "EXTRAP__1up"
           , "FATJET_JER__1up"
           , "FATJET_JMR__1up"
           , "FT_EFF_Eigen_B_0__1up"
           , "FT_EFF_Eigen_B_1__1up"
           , "FT_EFF_Eigen_B_2__1up"
           , "FT_EFF_Eigen_C_0__1up"
           , "FT_EFF_Eigen_C_1__1up"
           , "FT_EFF_Eigen_C_2__1up"
           , "FT_EFF_Eigen_C_3__1up"
           , "FT_EFF_Eigen_Light_0__1up"
           , "FT_EFF_Eigen_Light_1__1up"
           , "FT_EFF_Eigen_Light_2__1up"
           , "FT_EFF_Eigen_Light_3__1up"
           , "FT_EFF_Eigen_Light_4__1up"
           , "FT_EFF_extrapolation__1up"
           , "FT_EFF_extrapolation_from_charm__1up"
           , "JET_GroupedNP_1__1up"
           , "JET_GroupedNP_2__1up"
           , "JET_GroupedNP_3__1up"
           , "JET_JER_SINGLE_NP__1up"
           , "JET_Rtrk_Baseline__1up"
           , "JET_Rtrk_Modelling__1up"
           , "JET_Rtrk_Tracking__1up"
           , "MET_SoftTrk_ResoPara__1up"
           , "MET_SoftTrk_ResoPerp__1up"
           , "MET_SoftTrk_Scale__1up"
           -- , "MODEL_TTbar_Herwig__1up"
           -- , "MODEL_TTbar_aMcAtNlo__1up"
           -- , "MODEL_TTbar_rad__1up"
           , "MUONS_ID__1up"
           , "MUONS_MS__1up"
           , "MUONS_SCALE__1up"
           , "MUON_EFF_STAT__1up"
           , "MUON_EFF_SYS__1up"
           , "MUON_EFF_TrigStatUncertainty__1up"
           , "MUON_EFF_TrigSystUncertainty__1up"
           , "MUON_ISO_STAT__1up"
           , "MUON_ISO_SYS__1up"
           , "MUON_TTVA_STAT__1up"
           , "MUON_TTVA_SYS__1up"
           ]

main :: IO ()
main = do infiles <- getArgs
          let f Nothing = error "parse error."
              f (Just m) = fmap (rebinH 10 . snd) m

          let g k = let (procN, regN) = break (== '_') . takeWhileEnd (/= '/') . dropWhileEnd (/= '_') $ k
                    in  (ProcName procN, RegName . init . tail $ regN)

          fhs <- traverse (\fn -> (fn,) . decodeStrict <$> BS.readFile fn) infiles
                :: IO [(String, Maybe (Map String (String, Hist)))]

          let hs' = M.fromList $ map (\(ss, m) -> (g ss, f m)) fhs

          let regH regname procname varname = hs' M.! (procname, regname) M.! varname
          let procHs procname varname = M.fromList [(regn, regH regn procname varname) | regn <- regions]
          let varHs varname = M.fromList [(procn, procHs procn varname) | procn <- processes]

          let systHs = M.fromList $ map (\v -> (v, varHs v)) shapeNPs

          let nomHs = M.fromList $ map (\p -> (p, procHs p "nominal")) processes :: Prediction
          let dataHs = M.fromList $ map (\r -> (r, fmap round $ regH r "data" "nominal")) regions :: Dataset
          -- let dataHs = expectedData nomHs

          let sysDiffs = fmap (M.intersectionWith (M.intersectionWith (flip subH)) nomHs) systHs

          let shapeSysts = map (\(n, p) -> shapeParam standard n p) $ M.toList sysDiffs
          let normSysts = [ procNormParam (uniformDistr (-100.0) 100.0) 1.0 "HVTWHlvqq2000"
                          , procNormParam standard 0.3 "TTbar"
                          , procNormParam standard 0.3 "Wb"
                          , procNormParam standard 0.2 "Wc"
                          , procNormParam standard 0.1 "Wl"
                          ]

          let systs = normSysts ++ shapeSysts

          putStrLn . showList' . ("likelihood" :) $ map mpName systs

          -- start every NP at 0.0.
          let initial = map (const 0.0) systs

          let llh = modelLLH dataHs nomHs systs :: [Double] -> Double

          let t = Target llh Nothing
          let c = Chain t (llh initial) initial Nothing

          -- TODO
          -- I think the signal normalization needs a different step
          -- size...
          let trans = metropolis 0.1
          -- let trans = slice 0.1

          withSystemRandom . asGenIO $
                \gen -> chain trans c gen
                     =$ (dropC 5000 >> takeEveryC 20 =$ takeC 10000)
                     $$ mapM_C (\(Chain _ lhood xs _) -> putStrLn . showList' . fmap show $ lhood:xs)

max_stars_repo_name: cspollard/HMCMC

max_issues_repo_name: cspollard/HMCMC

max_forks_repo_name: cspollard/HMCMC


text:
module StructureKit.Prelude
  ( module Exports,
    guarded,
    mapCoercible,
    mapCoercible2,
    traverseCoercible,
    showBin,
    showBinFinite,

    -- * --
    BVec,
    UVec,
  )
where

import Control.Applicative as Exports hiding (WrappedArrow (..))
import Control.Arrow as Exports hiding (first, second)
import Control.Category as Exports
import Control.Concurrent as Exports
import Control.DeepSeq as Exports
import Control.Exception as Exports
import Control.Monad as Exports hiding (fail, forM, forM_, mapM, mapM_, msum, sequence, sequence_)
import Control.Monad.Fail as Exports
import Control.Monad.Fix as Exports hiding (fix)
import Control.Monad.IO.Class as Exports
import Control.Monad.Primitive as Exports
import Control.Monad.ST as Exports
import Control.Monad.Trans.Class as Exports
import Control.Monad.Trans.Cont as Exports hiding (callCC, shift)
import Control.Monad.Trans.Except as Exports (Except, ExceptT (ExceptT), catchE, except, mapExcept, mapExceptT, runExcept, runExceptT, throwE, withExcept, withExceptT)
import Control.Monad.Trans.Maybe as Exports
import Control.Monad.Trans.Reader as Exports (Reader, ReaderT (ReaderT), ask, mapReader, mapReaderT, runReader, runReaderT, withReader, withReaderT)
import Control.Monad.Trans.State.Strict as Exports (State, StateT (StateT), evalState, evalStateT, execState, execStateT, mapState, mapStateT, runState, runStateT, withState, withStateT)
import Control.Monad.Trans.Writer.Strict as Exports (Writer, WriterT (..), execWriter, execWriterT, mapWriter, mapWriterT, runWriter)
import Data.Bifunctor as Exports
import Data.Bits as Exports
import Data.Bool as Exports
import Data.ByteString as Exports (ByteString)
import Data.Char as Exports
import Data.Coerce as Exports
import Data.Complex as Exports
import Data.Data as Exports
import Data.Dynamic as Exports
import Data.Either as Exports
import Data.Fixed as Exports
import Data.Foldable as Exports hiding (toList)
import Data.Function as Exports hiding (id, (.))
import Data.Functor as Exports
import Data.Functor.Compose as Exports
import Data.Functor.Identity as Exports
import Data.Group as Exports
import Data.Hashable as Exports (Hashable (hash, hashWithSalt))
import Data.IORef as Exports
import Data.Int as Exports
import Data.IntMap.Strict as Exports (IntMap)
import Data.IntSet as Exports (IntSet)
import Data.Ix as Exports
import Data.List as Exports hiding (all, and, any, concat, concatMap, elem, filter, find, foldl, foldl', foldl1, foldr, foldr1, isSubsequenceOf, map, mapAccumL, mapAccumR, maximum, maximumBy, minimum, minimumBy, notElem, or, product, singleton, sortOn, sum, uncons)
import Data.List.NonEmpty as Exports (NonEmpty (..))
import Data.Maybe as Exports hiding (catMaybes, mapMaybe)
import Data.Monoid as Exports hiding (Alt)
import Data.Ord as Exports
import Data.Primitive as Exports hiding (alignment, sizeOf)
import Data.Proxy as Exports
import Data.Ratio as Exports
import Data.STRef as Exports
import Data.Semigroup as Exports hiding (First (..), Last (..))
import Data.String as Exports
import Data.String.ToString as Exports (ToString (..))
import Data.Text as Exports (Text)
import Data.Traversable as Exports
import Data.Tuple as Exports
import Data.Unique as Exports
import qualified Data.Vector as BVec
import qualified Data.Vector.Unboxed as UVec
import Data.Version as Exports
import Data.Void as Exports
import Data.Word as Exports
import Debug.Trace as Exports
import DeferredFolds.Unfoldr as Exports (Unfoldr (..))
import Deque.Strict as Exports (Deque)
import Foreign.ForeignPtr as Exports
import Foreign.Ptr as Exports
import Foreign.StablePtr as Exports
import Foreign.Storable as Exports
import GHC.Conc as Exports hiding (orElse, threadWaitRead, threadWaitReadSTM, threadWaitWrite, threadWaitWriteSTM, withMVar)
import GHC.Exts as Exports (IsList (..), groupWith, inline, lazy, sortWith)
import GHC.Generics as Exports (Generic)
import GHC.IO.Exception as Exports
import GHC.OverloadedLabels as Exports
import Numeric as Exports hiding (showBin)
import System.Environment as Exports
import System.Exit as Exports
import System.IO as Exports (Handle, hClose)
import System.IO.Error as Exports
import System.IO.Unsafe as Exports
import System.Mem as Exports
import System.Mem.StableName as Exports
import System.Timeout as Exports
import Text.ParserCombinators.ReadP as Exports (ReadP, readP_to_S, readS_to_P)
import Text.ParserCombinators.ReadPrec as Exports (ReadPrec, readP_to_Prec, readPrec_to_P, readPrec_to_S, readS_to_Prec)
import Text.Printf as Exports (hPrintf, printf)
import Text.Read as Exports (Read (..), readEither, readMaybe)
import Unsafe.Coerce as Exports
import Witherable as Exports
import Prelude as Exports hiding (all, and, any, concat, concatMap, elem, fail, filter, foldl, foldl1, foldr, foldr1, id, map, mapM, mapM_, maximum, minimum, notElem, or, product, sequence, sequence_, sum, (.))

guarded :: Alternative f => (a -> Bool) -> a -> f a
guarded pred x = if pred x then pure x else empty

mapCoercible :: Coercible a b => (a -> a) -> b -> b
mapCoercible mapper =
  coerce . mapper . coerce

mapCoercible2 :: Coercible f g => (g a -> g b) -> f a -> f b
mapCoercible2 mapper =
  coerce . mapper . coerce

traverseCoercible :: (Coercible a b, Functor f) => (a -> f a) -> b -> f b
traverseCoercible f =
  fmap coerce . f . coerce

showBin :: (Integral a, Show a) => a -> String
showBin a = showIntAtBase 2 intToDigit a ""

showBinFinite :: (Integral a, Show a, FiniteBits a) => a -> String
showBinFinite a =
  let str = showBin a
      len = length str
      rem = finiteBitSize a - len
      pad = replicate rem '0'
   in pad <> str

type BVec = BVec.Vector

type UVec = UVec.Vector

max_stars_repo_name: nikita-volkov/structure-kit

max_issues_repo_name: nikita-volkov/structure-kit

max_forks_repo_name: nikita-volkov/structure-kit


text:
{-# LANGUAGE AllowAmbiguousTypes    #-}
{-# LANGUAGE BangPatterns           #-}
{-# LANGUAGE DefaultSignatures      #-}
{-# LANGUAGE DeriveGeneric          #-}
{-# LANGUAGE DerivingVia            #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE KindSignatures         #-}
{-# LANGUAGE LambdaCase             #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE NoStarIsType           #-}
{-# LANGUAGE RankNTypes             #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE TupleSections          #-}
{-# LANGUAGE TypeApplications       #-}
{-# LANGUAGE TypeFamilyDependencies #-}
{-# LANGUAGE TypeInType             #-}
{-# LANGUAGE TypeOperators          #-}
{-# LANGUAGE UndecidableInstances   #-}
{-# OPTIONS_GHC -fno-warn-orphans   #-}

-- |
-- Module      : Numeric.Opto.Ref
-- Copyright   : (c) Justin Le 2019
-- License     : BSD3
--
-- Maintainer  : justin@jle.im
-- Stability   : experimental
-- Portability : non-portable
--
-- Abstract over different types for mutable references of values.
module Numeric.Opto.Ref (
    MR(..), ML(..)
  ) where

import           Control.Monad.Primitive
import           Data.Complex
import           Data.Mutable
import           GHC.TypeNats
import qualified Data.Vector.Generic                 as VG
import qualified Data.Vector.Generic.Mutable.Sized   as SMVG
import qualified Data.Vector.Generic.Sized           as SVG
import qualified Data.Vector.Storable                as VS
import qualified Numeric.LinearAlgebra               as HU
import qualified Numeric.LinearAlgebra.Devel         as HU
import qualified Numeric.LinearAlgebra.Static        as H
import qualified Numeric.LinearAlgebra.Static.Vector as H

instance VG.Vector v a => Mutable s (SVG.Vector v n a) where
    type Ref s (SVG.Vector v n a) = SVG.MVector (VG.Mutable v) n s a
    thawRef         = SVG.thaw
    freezeRef       = SVG.freeze
    copyRef         = SVG.copy
    moveRef         = SMVG.move
    cloneRef        = SMVG.clone
    unsafeThawRef   = SVG.unsafeThaw
    unsafeFreezeRef = SVG.unsafeFreeze

instance HU.Element a => Mutable s (HU.Matrix a) where
    type Ref s (HU.Matrix a) = HU.STMatrix s a
    thawRef x   = stToPrim $ HU.thawMatrix x
    freezeRef v = stToPrim $ HU.freezeMatrix v
    copyRef v x = stToPrim $ HU.setMatrix v 0 0 x
    moveRef         = undefined
    cloneRef        = undefined
    unsafeThawRef   = undefined
    unsafeFreezeRef = undefined

-- | Mutable ref for hmatrix's statically sized vector types, 'H.R' and
-- 'H.C'.
newtype MR s n a = MR { getMR :: SVG.MVector VS.MVector n s a }

instance KnownNat n => Mutable s (H.R n) where
    type Ref s (H.R n) = MR s n Double

    thawRef = fmap MR . thawRef . H.rVec
    freezeRef = fmap H.vecR . freezeRef . getMR
    copyRef (MR v) x = copyRef v (H.rVec x)
    moveRef         = undefined
    cloneRef        = undefined
    unsafeThawRef   = undefined
    unsafeFreezeRef = undefined

instance KnownNat n => Mutable s (H.C n) where
    type Ref s (H.C n) = MR s n (Complex Double)

    thawRef = fmap MR . thawRef . H.cVec
    freezeRef = fmap H.vecC . freezeRef . getMR
    copyRef (MR v) x = copyRef v (H.cVec x)
    moveRef         = undefined
    cloneRef        = undefined
    unsafeThawRef   = undefined
    unsafeFreezeRef = undefined

-- | Mutable ref for hmatrix's statically sized matrix types, 'H.L' and
-- 'H.M'.
newtype ML s n k a = ML { getML :: SVG.MVector VS.MVector (n * k) s a }

instance (KnownNat n, KnownNat k) => Mutable s (H.L n k) where
    type Ref s (H.L n k) = ML s n k Double

    thawRef = fmap ML . thawRef . H.lVec
    freezeRef = fmap H.vecL . freezeRef . getML
    copyRef (ML v) x = copyRef v (H.lVec x)
    moveRef         = undefined
    cloneRef        = undefined
    unsafeThawRef   = undefined
    unsafeFreezeRef = undefined

instance (KnownNat n, KnownNat k) => Mutable s (H.M n k) where
    type Ref s (H.M n k) = ML s n k (Complex Double)

    thawRef = fmap ML . thawRef . H.mVec
    freezeRef = fmap H.vecM . freezeRef . getML
    copyRef (ML v) x = copyRef v (H.mVec x)
    moveRef         = undefined
    cloneRef        = undefined
    unsafeThawRef   = undefined
    unsafeFreezeRef = undefined

max_stars_repo_name: mstksg/opto

max_issues_repo_name: mstksg/opto

max_forks_repo_name: mstksg/opto


text:
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric, FlexibleContexts #-}
-- |
-- Module    : Statistics.Sample.KernelDensity.Simple
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Kernel density estimation code, providing non-parametric ways to
-- estimate the probability density function of a sample.
--
-- The techniques used by functions in this module are relatively
-- fast, but they generally give inferior results to the KDE function
-- in the main 'Statistics.KernelDensity' module (due to the
-- oversmoothing documented for 'bandwidth' below).

module Statistics.Sample.KernelDensity.Simple
    {-# DEPRECATED "Use Statistics.Sample.KernelDensity instead." #-}
    (
    -- * Simple entry points
      epanechnikovPDF
    , gaussianPDF
    -- * Building blocks
    -- These functions may be useful if you need to construct a kernel
    -- density function estimator other than the ones provided in this
    -- module.

    -- ** Choosing points from a sample
    , Points(..)
    , choosePoints
    -- ** Bandwidth estimation
    , Bandwidth
    , bandwidth
    , epanechnikovBW
    , gaussianBW
    -- ** Kernels
    , Kernel
    , epanechnikovKernel
    , gaussianKernel
    -- ** Low-level estimation
    , estimatePDF
    , simplePDF
    -- * References
    -- $references
    ) where

import Data.Data (Data, Typeable)
import GHC.Generics (Generic)
import Numeric.MathFunctions.Constants (m_1_sqrt_2, m_2_sqrt_pi)
import Prelude hiding (sum)
import Statistics.Function (minMax)
import Statistics.Sample (stdDev)
import Statistics.Sample.Internal (sum)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U

-- | Points from the range of a 'Sample'.
newtype Points = Points {
      fromPoints :: U.Vector Double
    } deriving (Eq, Read, Show, Typeable, Data, Generic)

-- | Bandwidth estimator for an Epanechnikov kernel.
epanechnikovBW :: Double -> Bandwidth
epanechnikovBW n = (80 / (n * m_2_sqrt_pi)) ** 0.2

-- | Bandwidth estimator for a Gaussian kernel.
gaussianBW :: Double -> Bandwidth
gaussianBW n = (4 / (n * 3)) ** 0.2

-- | The width of the convolution kernel used.
type Bandwidth = Double

-- | Compute the optimal bandwidth from the observed data for the
-- given kernel.
--
-- This function uses an estimate based on the standard deviation of a
-- sample (due to Deheuvels), which performs reasonably well for
-- unimodal distributions but leads to oversmoothing for more complex
-- ones.
bandwidth :: G.Vector v Double =>
             (Double -> Bandwidth)
          -> v Double
          -> Bandwidth
bandwidth kern values = stdDev values * kern (fromIntegral $ G.length values)

-- | Choose a uniform range of points at which to estimate a sample's
-- probability density function.
--
-- If you are using a Gaussian kernel, multiply the sample's bandwidth
-- by 3 before passing it to this function.
--
-- If this function is passed an empty vector, it returns values of
-- positive and negative infinity.
choosePoints :: G.Vector v Double =>
                Int             -- ^ Number of points to select, /n/
             -> Double          -- ^ Sample bandwidth, /h/
             -> v Double        -- ^ Input data
             -> Points
choosePoints n h sample = Points . U.map f $ U.enumFromTo 0 n'
  where lo     = a - h
        hi     = z + h
        (a, z) = minMax sample
        d      = (hi - lo) / fromIntegral n'
        f i    = lo + fromIntegral i * d
        n'     = n - 1

-- | The convolution kernel.  Its parameters are as follows:
--
-- * Scaling factor, 1\//nh/
--
-- * Bandwidth, /h/
--
-- * A point at which to sample the input, /p/
--
-- * One sample value, /v/
type Kernel =  Double
            -> Double
            -> Double
            -> Double
            -> Double

-- | Epanechnikov kernel for probability density function estimation.
epanechnikovKernel :: Kernel
epanechnikovKernel f h p v
    | abs u <= 1 = f * (1 - u * u)
    | otherwise  = 0
    where u = (v - p) / (h * 0.75)

-- | Gaussian kernel for probability density function estimation.
gaussianKernel :: Kernel
gaussianKernel f h p v = exp (-0.5 * u * u) * g
    where u = (v - p) / h
          g = f * 0.5 * m_2_sqrt_pi * m_1_sqrt_2

-- | Kernel density estimator, providing a non-parametric way of
-- estimating the PDF of a random variable.
estimatePDF :: G.Vector v Double =>
               Kernel           -- ^ Kernel function
            -> Bandwidth        -- ^ Bandwidth, /h/
            -> v Double         -- ^ Sample data
            -> Points           -- ^ Points at which to estimate
            -> U.Vector Double
estimatePDF kernel h sample
    | n < 2     = errorShort "estimatePDF"
    | otherwise = U.map k . fromPoints
  where
    k p = sum . G.map (kernel f h p) $ sample
    f   = 1 / (h * fromIntegral n)
    n   = G.length sample
{-# INLINE estimatePDF #-}

-- | A helper for creating a simple kernel density estimation function
-- with automatically chosen bandwidth and estimation points.
simplePDF :: G.Vector v Double =>
             (Double -> Double) -- ^ Bandwidth function
          -> Kernel             -- ^ Kernel function
          -> Double             -- ^ Bandwidth scaling factor (3 for a Gaussian kernel, 1 for all others)
          -> Int                -- ^ Number of points at which to estimate
          -> v Double           -- ^ sample data
          -> (Points, U.Vector Double)
simplePDF fbw fpdf k numPoints sample =
    (points, estimatePDF fpdf bw sample points)
  where points = choosePoints numPoints (bw*k) sample
        bw     = bandwidth fbw sample
{-# INLINE simplePDF #-}

-- | Simple Epanechnikov kernel density estimator.  Returns the
-- uniformly spaced points from the sample range at which the density
-- function was estimated, and the estimates at those points.
epanechnikovPDF :: G.Vector v Double =>
                   Int          -- ^ Number of points at which to estimate
                -> v Double     -- ^ Data sample
                -> (Points, U.Vector Double)
epanechnikovPDF = simplePDF epanechnikovBW epanechnikovKernel 1

-- | Simple Gaussian kernel density estimator.  Returns the uniformly
-- spaced points from the sample range at which the density function
-- was estimated, and the estimates at those points.
gaussianPDF :: G.Vector v Double =>
               Int              -- ^ Number of points at which to estimate
            -> v Double         -- ^ Data sample
            -> (Points, U.Vector Double)
gaussianPDF = simplePDF gaussianBW gaussianKernel 3

errorShort :: String -> a
errorShort func = error ("Statistics.KernelDensity." ++ func ++
                        ": at least two points required")

-- $references
--
-- * Deheuvels, P. (1977) Estimation non paramétrique de la densité
--   par histogrammes
--   généralisés. Mhttp://archive.numdam.org/article/RSA_1977__25_3_5_0.pdf>

max_stars_repo_name: vmchale/statistics

max_issues_repo_name: vmchale/statistics

max_forks_repo_name: vmchale/statistics


text:
{-# LANGUAGE Rank2Types #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Matrix.Tri
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Triangular views of matrices.
--

module Numeric.LinearAlgebra.Matrix.Tri (
    -- * Immutable interface

    -- ** Vector multiplication
    triMulVector,

    -- ** Matrix multiplication
    triMulMatrix,
    triMulMatrixWithScale,

    -- ** Vector solving
    triSolvVector,

    -- ** Matrix solving
    triSolvMatrix,
    triSolvMatrixWithScale,

    -- * Mutable interface
    triCreate,

    -- ** Vector multiplication
    triMulVectorM_,

    -- ** Matrix multiplication
    triMulMatrixM_,
    triMulMatrixWithScaleM_,

    -- ** Vector solving
    triSolvVectorM_,

    -- ** Matrix solving
    triSolvMatrixM_,
    triSolvMatrixWithScaleM_,

    ) where

import Control.Monad( when )
import Control.Monad.ST( ST, runST, unsafeIOToST )
import Text.Printf( printf )

import Numeric.LinearAlgebra.Vector( Vector, STVector )
import qualified Numeric.LinearAlgebra.Vector as V
import Numeric.LinearAlgebra.Matrix.Base( Matrix )
import Numeric.LinearAlgebra.Matrix.STBase( STMatrix, RMatrix )
import qualified Numeric.LinearAlgebra.Matrix.STBase as M
import Numeric.LinearAlgebra.Types
import qualified Foreign.BLAS as BLAS

-- | A safe way to create and work with a mutable Tri Matrix before returning
-- an immutable one for later perusal.
triCreate :: (Storable e)
           => (forall s. ST s (Tri (STMatrix s) e))
           -> Tri Matrix e
triCreate mt = runST $ do
    (Tri u d ma) <- mt
    a <- M.unsafeFreeze ma
    return $ Tri u d a

-- | @triMulVector trans a x@ returns @op(a) * x@, where @op(a)@ is
-- determined by @trans@.
triMulVector :: (BLAS2 e)
             => Trans
             -> Tri Matrix e
             -> Vector e
             -> Vector e
triMulVector trans a x =
    V.create $ do
        x' <- V.newCopy x
        triMulVectorM_ trans a x'
        return x'

-- | @triMulMatrix side a b@
-- returns @alpha * op(a) * b@ when @side@ is @LeftSide@ and
-- @alpha * b * op(a)@ when @side@ is @RightSide@.  Operation
-- @op(a)@ is determined by @trans@.
triMulMatrix :: (BLAS3 e)
              => Side
              -> Trans -> Tri Matrix e
              -> Matrix e
              -> Matrix e
triMulMatrix side trans a b =
    M.create $ do
        b' <- M.newCopy b
        triMulMatrixM_ side trans a b'
        return b'

-- | @triMulMatrixWithScale alpha side trans a b@
-- returns @alpha * op(a) * b@ when @side@ is @LeftSide@ and
-- @alpha * b * op(a)@ when @side@ is @RightSide@.  Operation
-- @op(a)@ is determined by @trans@.
triMulMatrixWithScale :: (BLAS3 e)
                       => e
                       -> Side
                       -> Trans -> Tri Matrix e
                       -> Matrix e
                       -> Matrix e
triMulMatrixWithScale alpha side trans a b =
    M.create $ do
        b' <- M.newCopy b
        triMulMatrixWithScaleM_ alpha side trans a b'
        return b'

-- | @triMulVectorM_ a x@ sets @x := op(a) * x@, where @op(a)@ is determined
-- by @trans@.
triMulVectorM_ :: (RMatrix m, BLAS2 e)
               => Trans -> Tri m e
               -> STVector s e
               -> ST s ()
triMulVectorM_ trans (Tri uplo diag a) x = do
    (ma,na) <- M.getDim a
    nx <- V.getDim x
    let n = nx

    when (ma /= na) $ error $
        printf ("triMulVectorM_"
                ++ " _"
                ++ " (Tri _ _ <matrix with dim (%d,%d)>)"
                ++ " _"
                ++ ": matrix is not square")
               ma na

    when ((not . and) [ (ma,na) == (n,n)
                      , nx == n
                      ]) $ error $
        printf ("triMulVectorM_"
                ++ " _"
                ++ " (Tri _ _ <matrix with dim (%d,%d)>)"
                ++ " <vector with dim %d>"
                ++ ": dimension mismatch")
               ma na
               nx

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        V.unsafeWith x $ \px ->
            BLAS.trmv uplo trans diag n pa lda px 1

-- | @triMulMatrixM_ side trans a b@
-- sets @b := op(a) * b@ when @side@ is @LeftSide@ and
-- @b := b * op(a)@ when @side@ is @RightSide@.  Operation
-- @op(a)@ is determined by @trans@.
triMulMatrixM_ :: (RMatrix m, BLAS3 e)
               => Side
               -> Trans -> Tri m e
               -> STMatrix s e
               -> ST s ()
triMulMatrixM_ = triMulMatrixWithScaleM_ 1

-- | @triMulMatrixWithScaleM_ alpha side trans a b@
-- sets @b := alpha * op(a) * b@ when @side@ is @LeftSide@ and
-- @b := alpha * b * op(a)@ when @side@ is @RightSide@.  Operation
-- @op(a)@ is determined by @trans@.
triMulMatrixWithScaleM_ :: (RMatrix m, BLAS3 e)
                         => e
                         -> Side
                         -> Trans -> Tri m e
                         -> STMatrix s e
                         -> ST s ()
triMulMatrixWithScaleM_ alpha side trans (Tri uplo diag a) b = do
    (ma,na) <- M.getDim a
    (mb,nb) <- M.getDim b
    let (m,n) = (mb,nb)

    when (ma /= na) $ error $
        printf ("triMulMatrixWithScaleM_"
                ++ " _"
                ++ " _"
                ++ " _"
                ++ " (Tri _ _ <matrix with dim (%d,%d)>)"
                ++ " _"
                ++ ": matrix is not square")
               ma na

    when ((not . and) [ case side of LeftSide  -> (ma,na) == (m,m)
                                     RightSide -> (ma,na) == (n,n)
                      , (mb, nb ) == (m,n)
                      ]) $ error $
        printf ("triMulMatrixWithScaleM_"
                ++ " _"
                ++ " %s"
                ++ " _"
                ++ " (Tri _ _ <matrix with dim (%d,%d)>)"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch")
               (show side)
               ma na
               mb nb

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        M.unsafeWith b $ \pb ldb ->
            BLAS.trmm side uplo trans diag m n alpha pa lda pb ldb

-- | @triSolvVector trans a x@ returns @op(a) \\ x@, where @op(a)@ is
-- determined by @trans@.
triSolvVector :: (BLAS2 e)
             => Trans
             -> Tri Matrix e
             -> Vector e
             -> Vector e
triSolvVector trans a x =
    V.create $ do
        x' <- V.newCopy x
        triSolvVectorM_ trans a x'
        return x'

-- | @triSolvMatrix side a b@
-- returns @alpha * op(a) \\ b@ when @side@ is @LeftSide@ and
-- @alpha * b * op(a)@ when @side@ is @RightSide@.  Operation
-- @op(a)@ is determined by @trans@.
triSolvMatrix :: (BLAS3 e)
              => Side
              -> Trans -> Tri Matrix e
              -> Matrix e
              -> Matrix e
triSolvMatrix side trans a b =
    M.create $ do
        b' <- M.newCopy b
        triSolvMatrixM_ side trans a b'
        return b'

-- | @triSolvMatrixWithScale alpha side trans a b@
-- returns @alpha * op(a) \\ b@ when @side@ is @LeftSide@ and
-- @alpha * b * op(a)@ when @side@ is @RightSide@.  Operation
-- @op(a)@ is determined by @trans@.
triSolvMatrixWithScale :: (BLAS3 e)
                       => e
                       -> Side
                       -> Trans -> Tri Matrix e
                       -> Matrix e
                       -> Matrix e
triSolvMatrixWithScale alpha side trans a b =
    M.create $ do
        b' <- M.newCopy b
        triSolvMatrixWithScaleM_ alpha side trans a b'
        return b'

-- | @triSolvVectorM_ a x@ sets @x := op(a) \\ x@, where @op(a)@ is determined
-- by @trans@.
triSolvVectorM_ :: (RMatrix m, BLAS2 e)
               => Trans -> Tri m e
               -> STVector s e
               -> ST s ()
triSolvVectorM_ trans (Tri uplo diag a) x = do
    (ma,na) <- M.getDim a
    nx <- V.getDim x
    let n = nx

    when (ma /= na) $ error $
        printf ("triSolvVectorM_"
                ++ " _"
                ++ " (Tri _ _ <matrix with dim (%d,%d)>)"
                ++ " _"
                ++ ": matrix is not square")
               ma na

    when ((not . and) [ (ma,na) == (n,n)
                      , nx == n
                      ]) $ error $
        printf ("triSolvVectorM_"
                ++ " _"
                ++ " (Tri _ _ <matrix with dim (%d,%d)>)"
                ++ " <vector with dim %d>"
                ++ ": dimension mismatch")
               ma na
               nx

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        V.unsafeWith x $ \px ->
            BLAS.trsv uplo trans diag n pa lda px 1

-- | @triSolvMatrixM_ side trans a b@
-- sets @b := op(a) \\ b@ when @side@ is @LeftSide@ and
-- @b := b * op(a)@ when @side@ is @RightSide@.  Operation
-- @op(a)@ is determined by @trans@.
triSolvMatrixM_ :: (RMatrix m, BLAS3 e)
               => Side
               -> Trans -> Tri m e
               -> STMatrix s e
               -> ST s ()
triSolvMatrixM_ = triSolvMatrixWithScaleM_ 1

-- | @triSolvMatrixWithScaleM_ alpha side trans a b@
-- sets @b := alpha * op(a) \\ b@ when @side@ is @LeftSide@ and
-- @b := alpha * b * op(a)@ when @side@ is @RightSide@.  Operation
-- @op(a)@ is determined by @trans@.
triSolvMatrixWithScaleM_ :: (RMatrix m, BLAS3 e)
                         => e
                         -> Side
                         -> Trans -> Tri m e
                         -> STMatrix s e
                         -> ST s ()
triSolvMatrixWithScaleM_ alpha side trans (Tri uplo diag a) b = do
    (ma,na) <- M.getDim a
    (mb,nb) <- M.getDim b
    let (m,n) = (mb,nb)

    when (ma /= na) $ error $
        printf ("triSolvMatrixWithScaleM_"
                ++ " _"
                ++ " _"
                ++ " _"
                ++ " (Tri _ _ <matrix with dim (%d,%d)>)"
                ++ " _"
                ++ ": matrix is not square")
               ma na

    when ((not . and) [ case side of LeftSide  -> (ma,na) == (m,m)
                                     RightSide -> (ma,na) == (n,n)
                      , (mb, nb ) == (m,n)
                      ]) $ error $
        printf ("triSolvMatrixWithScaleM_"
                ++ " _"
                ++ " %s"
                ++ " _"
                ++ " (Tri _ _ <matrix with dim (%d,%d)>)"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch")
               (show side)
               ma na
               mb nb

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        M.unsafeWith b $ \pb ldb ->
            BLAS.trsm side uplo trans diag m n alpha pa lda pb ldb

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas