text:
{-# LANGUAGE BangPatterns           #-}
{-# LANGUAGE DataKinds              #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE LambdaCase             #-}
{-# LANGUAGE PolyKinds              #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE TupleSections          #-}
{-# LANGUAGE TypeApplications       #-}
{-# LANGUAGE TypeFamilyDependencies #-}

module Backprop.Learn.Model (
    module M, Backprop(..)
  -- * Running and Grad
  , runModel, runModelStoch, runModelStateless, runModelStochStateless
  , gradModel, gradModelStoch
  -- * Work with parameters
  , initParam, initParamNormal
  , encodeParam, decodeParam, decodeParamOrFail, saveParam, loadParam, loadParamOrFail
  -- * Iterated runners
  , iterateModel, iterateModelM, iterateModelStoch
  , scanModel, scanModelStoch
  -- * No final state
  , iterateModel_, iterateModelM_, iterateModelStoch_
  , scanModel_, scanModelStoch_
  -- * "Prime" runners
  , primeModel, primeModelStoch, selfPrime, selfPrimeM
  ) where

import           Backprop.Learn.Initialize
import           Backprop.Learn.Model.Combinator  as M
import           Backprop.Learn.Model.Function    as M
import           Backprop.Learn.Model.Neural      as M
import           Backprop.Learn.Model.Neural.LSTM as M
import           Backprop.Learn.Model.Parameter   as M
import           Backprop.Learn.Model.Regression  as M
import           Backprop.Learn.Model.State       as M
import           Backprop.Learn.Model.Stochastic  as M
import           Backprop.Learn.Model.Types       as M
import           Control.Monad.Primitive
import           Control.Monad.ST
import           Control.Monad.Trans.State
import           Data.Bifunctor
import           Data.Foldable
import           Data.Functor.Identity
import           Data.Type.Functor.Product
import           Data.Type.Tuple
import           Data.Word
import           Numeric.Backprop
import           Statistics.Distribution
import qualified Data.Binary                      as Bi
import qualified Data.ByteString.Lazy             as BSL
import qualified Data.Vector.Unboxed              as VU
import qualified System.Random.MWC                as MWC

-- TODO: this can be more efficient by breaking out into separate functions
runModel
    :: forall p s a b. (AllConstrainedProd Backprop s, Backprop b)
    => Model p s a b
    -> TMaybe p
    -> a
    -> TMaybe s
    -> (b, TMaybe s)
runModel f mp x ms = evalBP0 go
  where
    go :: forall z. Reifies z W => BVar z (b, TMaybe s)
    go = case ms' of
        PNothing    -> T2 y $ auto PNothing
        PJust s' -> T2 y $ isoVar (PJust . TF) (getTF . fromPJust) s'
      where
        y :: BVar z b
        ms' :: PMaybe (BVar z) s
        (y, ms') = runLearn f (mapProd (auto . getTF) mp)
                              (auto x)
                              (mapProd (auto . getTF) ms)

-- TODO: this can be more efficient by breaking out into separate functions
runModelStoch
    :: forall p s a b m. (AllConstrainedProd Backprop s, Backprop b, PrimMonad m)
    => Model p s a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> a
    -> TMaybe s
    -> m (b, TMaybe s)
runModelStoch f g mp x ms = do
    -- TODO: is this the best way to handle this?
    seed <- MWC.uniformVector @_ @Word32 @VU.Vector g 2
    pure $ evalBP0 (runST (go seed))
  where
    go  :: forall q z. Reifies z W
        => VU.Vector Word32
        -> ST q (BVar z (b, TMaybe s))
    go seed = do
      g' <- MWC.initialize seed
      (y :: BVar z b, ms') <- runLearnStoch f g'
        (mapProd (auto . getTF) mp)
        (auto x)
        (mapProd (auto . getTF) ms)
      pure $ case ms' of
        PNothing    -> T2 y $ auto PNothing
        PJust s' -> T2 y $ isoVar (PJust . TF) (getTF . fromPJust) s'

runModelStateless
    :: Model p 'Nothing a b
    -> TMaybe p
    -> a
    -> b
runModelStateless f = \case
    PNothing       -> evalBP  (runLearnStateless f PNothing  )
    PJust (TF p) -> evalBP2 (runLearnStateless f . PJust) p

runModelStochStateless
    :: PrimMonad m
    => Model p 'Nothing a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> a
    -> m b
runModelStochStateless f g mp x = do
    seed <- MWC.uniformVector @_ @Word32 @VU.Vector g 2
    pure $ case mp of
      PNothing       -> evalBP  (\x' -> runST $ do
          g' <- MWC.initialize seed
          runLearnStochStateless f g' PNothing x'
        ) x
      PJust (TF p) -> evalBP2 (\p' x' -> runST $ do
          g' <- MWC.initialize seed
          runLearnStochStateless f g' (PJust p') x'
        ) p x

gradModel
    :: (Backprop a, Backprop b, AllConstrainedProd Backprop p)
    => Model p 'Nothing a b
    -> TMaybe p
    -> a
    -> (TMaybe p, a)
gradModel f = \case
    PNothing       ->       (PNothing,)    . gradBP  (runLearnStateless f   PNothing)
    PJust (TF p) -> first (PJust . TF) . gradBP2 (runLearnStateless f . PJust) p

gradModelStoch
    :: (Backprop a, Backprop b, AllConstrainedProd Backprop p, PrimMonad m)
    => Model p 'Nothing a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> a
    -> m (TMaybe p, a)
gradModelStoch f g mp x = do
    seed <- MWC.uniformVector @_ @Word32 @VU.Vector g 2
    pure $ case mp of
      PNothing       ->       (PNothing,)    $ gradBP  (\x' -> runST $ do
          g' <- MWC.initialize seed
          runLearnStochStateless f g' PNothing x'
        ) x
      PJust (TF p) -> first (PJust . TF) $ gradBP2 (\p' x' -> runST $ do
          g' <- MWC.initialize seed
          runLearnStochStateless f g' (PJust p') x'
        ) p x

iterateModel
    :: (Backprop b, AllConstrainedProd Backprop s)
    => (b -> a)         -- ^ loop
    -> Int              -- ^ num times
    -> Model p s a b
    -> TMaybe p
    -> a
    -> TMaybe s
    -> ([b], TMaybe s)
iterateModel l n f p x = runIdentity . iterateModelM (Identity . l) n f p x

iterateModel_
    :: (Backprop b, AllConstrainedProd Backprop s)
    => (b -> a)         -- ^ loop
    -> Model p s a b
    -> TMaybe p
    -> a
    -> TMaybe s
    -> [b]
iterateModel_ l f p = go
  where
    go !x !s = y : go (l y) s'
      where
        (y, s') = runModel f p x s

selfPrime
    :: (Backprop b, AllConstrainedProd Backprop s)
    => (b -> a)         -- ^ loop
    -> Model p s a b
    -> TMaybe p
    -> a
    -> TMaybe s
    -> [TMaybe s]
selfPrime l f p = go
  where
    go !x !s = s' : go (l y) s'
      where
        (y, s') = runModel f p x s

iterateModelM
    :: (Backprop b, AllConstrainedProd Backprop s, Monad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> Model p s a b
    -> TMaybe p
    -> a
    -> TMaybe s
    -> m ([b], TMaybe s)
iterateModelM l n f p = go 0
  where
    go !i !x !s
      | i <= n    = do
          let (y, s') = runModel f p x s
          (ys, s'') <- flip (go (i + 1)) s' =<< l y
          pure (y : ys, s'')
      | otherwise = pure ([], s)

iterateModelM_
    :: (Backprop b, AllConstrainedProd Backprop s, Monad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> Model p s a b
    -> TMaybe p
    -> a
    -> TMaybe s
    -> m [b]
iterateModelM_ l n f p x = fmap fst . iterateModelM l n f p x

selfPrimeM
    :: (Backprop b, AllConstrainedProd Backprop s, Monad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> Model p s a b
    -> TMaybe p
    -> a
    -> TMaybe s
    -> m (TMaybe s)
selfPrimeM l n f p x = fmap snd . iterateModelM l n f p x

iterateModelStoch
    :: (Backprop b, AllConstrainedProd Backprop s, PrimMonad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> Model p s a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> a
    -> TMaybe s
    -> m ([b], TMaybe s)
iterateModelStoch l n f g p = go 0
  where
    go !i !x !s
      | i <= n    = do
          (y , s' ) <- runModelStoch f g p x s
          (ys, s'') <- flip (go (i + 1)) s' =<< l y
          pure (y : ys, s'')
      | otherwise = pure ([], s)

iterateModelStoch_
    :: (Backprop b, AllConstrainedProd Backprop s, PrimMonad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> Model p s a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> a
    -> TMaybe s
    -> m [b]
iterateModelStoch_ l n f g p x = fmap fst . iterateModelStoch l n f g p x

scanModel
    :: (Traversable t, Backprop b, AllConstrainedProd Backprop s)
    => Model p s a b
    -> TMaybe p
    -> t a
    -> TMaybe s
    -> (t b, TMaybe s)
scanModel f p = runState . traverse (state . runModel f p)

scanModel_
    :: (Traversable t, Backprop b, AllConstrainedProd Backprop s)
    => Model p s a b
    -> TMaybe p
    -> t a
    -> TMaybe s
    -> t b
scanModel_ f p xs = fst . scanModel f p xs

primeModel
    :: (Foldable t, Backprop b, AllConstrainedProd Backprop s)
    => Model p s a b
    -> TMaybe p
    -> t a
    -> TMaybe s
    -> TMaybe s
primeModel f p = execState . traverse_ (state . runModel f p)

scanModelStoch
    :: (Traversable t, Backprop b, AllConstrainedProd Backprop s, PrimMonad m)
    => Model p s a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> t a
    -> TMaybe s
    -> m (t b, TMaybe s)
scanModelStoch f g p = runStateT . traverse (StateT . runModelStoch f g p)

scanModelStoch_
    :: (Traversable t, Backprop b, AllConstrainedProd Backprop s, PrimMonad m)
    => Model p s a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> t a
    -> TMaybe s
    -> m (t b)
scanModelStoch_ f g p xs = fmap fst . scanModelStoch f g p xs

primeModelStoch
    :: (Foldable t, Backprop b, AllConstrainedProd Backprop s, PrimMonad m)
    => Model p s a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> t a
    -> TMaybe s
    -> m (TMaybe s)
primeModelStoch f g p = execStateT . traverse_ (StateT . runModelStoch f g p)

initParam
    :: (Initialize p, ContGen d, PrimMonad m)
    => model ('Just p) s a b                    -- ^ ignored
    -> d
    -> MWC.Gen (PrimState m)
    -> m p
initParam _ = initialize

initParamNormal
    :: (Initialize p, PrimMonad m)
    => model ('Just p) s a b                    -- ^ ignored
    -> Double
    -> MWC.Gen (PrimState m)
    -> m p
initParamNormal _ = initializeNormal

encodeParam
    :: Bi.Binary p
    => model ('Just p) s a b                    -- ^ ignored
    -> p
    -> BSL.ByteString
encodeParam _ = Bi.encode

decodeParam
    :: Bi.Binary p
    => model ('Just p) s a b                    -- ^ ignored
    -> BSL.ByteString
    -> p
decodeParam _ = Bi.decode

decodeParamOrFail
    :: Bi.Binary p
    => model ('Just p) s a b                    -- ^ ignored
    -> BSL.ByteString
    -> Either String p
decodeParamOrFail _ = bimap thrd thrd . Bi.decodeOrFail

saveParam
    :: Bi.Binary p
    => model ('Just p) s a b                    -- ^ ignored
    -> FilePath
    -> p
    -> IO ()
saveParam p fp = BSL.writeFile fp . encodeParam p

loadParam
    :: Bi.Binary p
    => model ('Just p) s a b                    -- ^ ignored
    -> FilePath
    -> IO p
loadParam p fp = decodeParam p <$> BSL.readFile fp

loadParamOrFail
    :: Bi.Binary p
    => model ('Just p) s a b                    -- ^ ignored
    -> FilePath
    -> IO (Either String p)
loadParamOrFail p fp = decodeParamOrFail p <$> BSL.readFile fp

thrd :: (a,b,c) -> c
thrd (_,_,z) = z

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
module Main where

import           Control.Arrow                         (first, (&&&))

import           Criterion
import           Criterion.IO                          (readRecords)
import           Criterion.Main                        (defaultConfig,
                                                        defaultMainWith)
import           Criterion.Types
import           Statistics.Resampling.Bootstrap       (estPoint)

import qualified Data.ByteString.Lazy                  as BL
import           Data.Csv                              as Csv
import           Data.List                             (isPrefixOf, sortBy, zip4)
import           Data.Ord                              (comparing)
import qualified Data.Vector.Unboxed                   as V

import           Test.QuickCheck

import           Graphics.Gnuplot.Advanced
import qualified Graphics.Gnuplot.Frame                as Frame
import qualified Graphics.Gnuplot.Frame.Option         as Opt
import qualified Graphics.Gnuplot.Frame.OptionSet      as Opts
import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D
import qualified Graphics.Gnuplot.LineSpecification    as LineSpec
import qualified Graphics.Gnuplot.Plot.TwoDimensional  as Plot2D
import qualified Graphics.Gnuplot.Terminal.PNG         as PNG

import           Quickselect

range :: [Int]
range = [0, 1000 .. 50000]

testData :: [(Int, IO (V.Vector Int))]
testData = map f range
  where f n = (n, generate $ (V.fromList <$> vectorOf n arbitrary))

rawOut = "../data/kth-element-haskell.raw"
csvOut = "../data/kth-element-haskell.csv"
pngOut = "../data/kth-element-haskell.png"

config :: Config
config = defaultConfig { resamples = 3, rawDataFile = Just rawOut }

main :: IO ()
main = do
  defaultMainWith config [
    bgroup "main" [
        bgroup "find" $ map (run quickselect) testData
        , bgroup "cheat" $ map (run sortselect) testData
        ]
    ]
  (Right records) <- readRecords rawOut
  let records' = map (\(Analysed r) -> r) records
  let find' = extractOutcomes (filterReports "main/find/" records')
      cheat' = extractOutcomes (filterReports "main/cheat/" records')
      results = zip3 find' cheat' range
  BL.writeFile csvOut (encode results)
  plot (PNG.cons pngOut) resultPlot
  return ()
  where run f (n, xs) = env xs $ \ ~(xs') -> bench (show n) $ nf (f n) xs'
        filterReports name = filter ((name `isPrefixOf`) . reportName)
        extractOutcomes = map (estPoint . anMean . reportAnalysis)

resultPlot :: Frame.T (Graph2D.T Int Double)
resultPlot =
  let lineSpec title =
        Graph2D.lineSpec $
        LineSpec.title title $
        LineSpec.deflt
      lineSpec1 = lineSpec "find"
      lineSpec2 = lineSpec "cheat"
      frameOpts =
        Opts.xLabel "input size" $
        Opts.yLabel "execution time (seconds)" $
        Opts.title "Performance of Kth Element" $
        Opts.add (Opt.custom "datafile separator" "") ["\",\""] $
        Opts.deflt
      path1 = fmap lineSpec1 $ Plot2D.pathFromFile Graph2D.lines csvOut 3 1
      path2 = fmap lineSpec2 $ Plot2D.pathFromFile Graph2D.lines csvOut 3 2
  in Frame.cons frameOpts $ mconcat [path1, path2]

ext: hs

max_stars_repo_name: mooreniemi/BAM

max_issues_repo_name: mooreniemi/BAM

max_forks_repo_name: mooreniemi/BAM


text:
{-# LANGUAGE NegativeLiterals#-}

{-# LANGUAGE NoMonomorphismRestriction, ExtendedDefaultRules#-}
module DocTest.Flat.Instances.Base where
import qualified DocTest
import Test.Tasty(TestTree,testGroup)
import Flat.Instances.Base
import Flat.Instances.Test
import Data.Fixed
import Data.Int
import Data.Complex(Complex(..))
import Numeric.Natural
import Data.Word
import Data.Ratio
import Flat.Run
import Data.Monoid
import qualified Data.List.NonEmpty as B
test = tstBits
y = 33

tests :: IO TestTree
tests = testGroup "Flat.Instances.Base" <$> sequence [  DocTest.test "src/Flat/Instances/Base.hs:104" "[ExpectedLine [LineChunk \"True\"]]" (DocTest.asPrint( let w = Just (11::Word8); a = Alt w <> Alt (Just 24) in tst a == tst w )),  DocTest.test "src/Flat/Instances/Base.hs:107" "[ExpectedLine [LineChunk \"True\"]]" (DocTest.asPrint( let w = Just (11::Word8); a = Alt Nothing <> Alt w in tst a == tst w )),  DocTest.test "src/Flat/Instances/Base.hs:145" "[ExpectedLine [LineChunk \"(True,0,\\\"\\\")\"]]" (DocTest.asPrint( test () )),  DocTest.test "src/Flat/Instances/Base.hs:158" "[ExpectedLine [LineChunk \"(True,1,\\\"0\\\")\"]]" (DocTest.asPrint( test False )),  DocTest.test "src/Flat/Instances/Base.hs:161" "[ExpectedLine [LineChunk \"(True,1,\\\"1\\\")\"]]" (DocTest.asPrint( test True )),  DocTest.test "src/Flat/Instances/Base.hs:176" "[ExpectedLine [LineChunk \"(True,8,\\\"01100001\\\")\"]]" (DocTest.asPrint( test 'a' )),  DocTest.test "src/Flat/Instances/Base.hs:181" "[ExpectedLine [LineChunk \"(True,16,\\\"11001000 00000001\\\")\"]]" (DocTest.asPrint( test 'È' )),  DocTest.test "src/Flat/Instances/Base.hs:184" "[ExpectedLine [LineChunk \"(True,24,\\\"10001101 10011100 00000001\\\")\"]]" (DocTest.asPrint( test '不' )),  DocTest.test "src/Flat/Instances/Base.hs:200" "[ExpectedLine [LineChunk \"(True,1,\\\"0\\\")\"]]" (DocTest.asPrint( test (Nothing::Maybe Bool) )),  DocTest.test "src/Flat/Instances/Base.hs:203" "[ExpectedLine [LineChunk \"(True,2,\\\"10\\\")\"]]" (DocTest.asPrint( test (Just False::Maybe Bool) )),  DocTest.test "src/Flat/Instances/Base.hs:209" "[ExpectedLine [LineChunk \"(True,2,\\\"00\\\")\"]]" (DocTest.asPrint( test (Left False::Either Bool ()) )),  DocTest.test "src/Flat/Instances/Base.hs:212" "[ExpectedLine [LineChunk \"(True,1,\\\"1\\\")\"]]" (DocTest.asPrint( test (Right ()::Either Bool ()) )),  DocTest.test "src/Flat/Instances/Base.hs:218" "[ExpectedLine [LineChunk \"(True,16,\\\"11110110 00000001\\\")\"]]" (DocTest.asPrint( test (MkFixed 123 :: Fixed E0) )),  DocTest.test "src/Flat/Instances/Base.hs:221" "[ExpectedLine [LineChunk \"True\"]]" (DocTest.asPrint( test (MkFixed 123 :: Fixed E0) == test (MkFixed 123 :: Fixed E2) )),  DocTest.test "src/Flat/Instances/Base.hs:234" "[ExpectedLine [LineChunk \"(True,8,\\\"00000000\\\")\"]]" (DocTest.asPrint( test (0::Word8) )),  DocTest.test "src/Flat/Instances/Base.hs:237" "[ExpectedLine [LineChunk \"(True,8,\\\"11111111\\\")\"]]" (DocTest.asPrint( test (255::Word8) )),  DocTest.test "src/Flat/Instances/Base.hs:281" "[ExpectedLine [LineChunk \"(True,8,\\\"01111111\\\")\"]]" (DocTest.asPrint( test (127::Word) )),  DocTest.test "src/Flat/Instances/Base.hs:286" "[ExpectedLine [LineChunk \"(True,16,\\\"11111110 00000001\\\")\"]]" (DocTest.asPrint( test (254::Word) )),  DocTest.test "src/Flat/Instances/Base.hs:291" "[ExpectedLine [LineChunk \"(True,24,\\\"10000000 10000000 00000010\\\")\"]]" (DocTest.asPrint( test (32768::Word32) )),  DocTest.test "src/Flat/Instances/Base.hs:296" "[ExpectedLine [LineChunk \"True\"]]" (DocTest.asPrint( all (test (3::Word) ==) [test (3::Word16),test (3::Word32),test (3::Word64)] )),  DocTest.test "src/Flat/Instances/Base.hs:309" "[ExpectedLine [LineChunk \"(True,8,\\\"00000000\\\")\"]]" (DocTest.asPrint( test (0::Natural) )),  DocTest.test "src/Flat/Instances/Base.hs:312" "[ExpectedLine [LineChunk \"(True,144,\\\"10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 00000010\\\")\"]]" (DocTest.asPrint( test (2^120::Natural) )),  DocTest.test "src/Flat/Instances/Base.hs:362" "[ExpectedLine [LineChunk \"(True,8,\\\"00000000\\\")\"]]" (DocTest.asPrint( test (0::Int) )),  DocTest.test "src/Flat/Instances/Base.hs:365" "[ExpectedLine [LineChunk \"(True,8,\\\"00000001\\\")\"]]" (DocTest.asPrint( test (-1::Int) )),  DocTest.test "src/Flat/Instances/Base.hs:368" "[ExpectedLine [LineChunk \"(True,8,\\\"00000010\\\")\"]]" (DocTest.asPrint( test (1::Int) )),  DocTest.test "src/Flat/Instances/Base.hs:371" "[ExpectedLine [LineChunk \"(True,8,\\\"00000011\\\")\"]]" (DocTest.asPrint( test (-2::Int) )),  DocTest.test "src/Flat/Instances/Base.hs:374" "[ExpectedLine [LineChunk \"(True,8,\\\"00000100\\\")\"]]" (DocTest.asPrint( test (2::Int) )),  DocTest.test "src/Flat/Instances/Base.hs:387" "[ExpectedLine [LineChunk \"(True,8,\\\"00000000\\\")\"]]" (DocTest.asPrint( test (0::Integer) )),  DocTest.test "src/Flat/Instances/Base.hs:390" "[ExpectedLine [LineChunk \"(True,8,\\\"00000001\\\")\"]]" (DocTest.asPrint( test (-1::Integer) )),  DocTest.test "src/Flat/Instances/Base.hs:393" "[ExpectedLine [LineChunk \"(True,8,\\\"00000010\\\")\"]]" (DocTest.asPrint( test (1::Integer) )),  DocTest.test "src/Flat/Instances/Base.hs:396" "[ExpectedLine [LineChunk \"(True,8,\\\"00011111\\\")\"]]" (DocTest.asPrint( test (-(2^4)::Integer) )),  DocTest.test "src/Flat/Instances/Base.hs:399" "[ExpectedLine [LineChunk \"(True,8,\\\"00100000\\\")\"]]" (DocTest.asPrint( test (2^4::Integer) )),  DocTest.test "src/Flat/Instances/Base.hs:402" "[ExpectedLine [LineChunk \"(True,144,\\\"11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 00000011\\\")\"]]" (DocTest.asPrint( test (-(2^120)::Integer) )),  DocTest.test "src/Flat/Instances/Base.hs:405" "[ExpectedLine [LineChunk \"(True,144,\\\"10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 00000100\\\")\"]]" (DocTest.asPrint( test (2^120::Integer) )),  DocTest.test "src/Flat/Instances/Base.hs:416" "[ExpectedLine [LineChunk \"(True,8,\\\"00000000\\\")\"]]" (DocTest.asPrint( test (0::Int8) )),  DocTest.test "src/Flat/Instances/Base.hs:419" "[ExpectedLine [LineChunk \"(True,8,\\\"11111110\\\")\"]]" (DocTest.asPrint( test (127::Int8) )),  DocTest.test "src/Flat/Instances/Base.hs:422" "[ExpectedLine [LineChunk \"(True,8,\\\"11111111\\\")\"]]" (DocTest.asPrint( test (-128::Int8) )),  DocTest.test "src/Flat/Instances/Base.hs:433" "[ExpectedLine [LineChunk \"(True,8,\\\"00000000\\\")\"]]" (DocTest.asPrint( test (0::Int16) )),  DocTest.test "src/Flat/Instances/Base.hs:436" "[ExpectedLine [LineChunk \"(True,8,\\\"00000010\\\")\"]]" (DocTest.asPrint( test (1::Int16) )),  DocTest.test "src/Flat/Instances/Base.hs:439" "[ExpectedLine [LineChunk \"(True,8,\\\"00000001\\\")\"]]" (DocTest.asPrint( test (-1::Int16) )),  DocTest.test "src/Flat/Instances/Base.hs:442" "[ExpectedLine [LineChunk \"(True,24,\\\"11111111 11111111 00000011\\\")\"]]" (DocTest.asPrint( test (minBound::Int16) )),  DocTest.test "src/Flat/Instances/Base.hs:447" "[ExpectedLine [LineChunk \"(True,24,\\\"11111110 11111111 00000011\\\")\"]]" (DocTest.asPrint( test (maxBound::Int16) )),  DocTest.test "src/Flat/Instances/Base.hs:460" "[ExpectedLine [LineChunk \"(True,8,\\\"00000000\\\")\"]]" (DocTest.asPrint( test (0::Int32) )),  DocTest.test "src/Flat/Instances/Base.hs:463" "[ExpectedLine [LineChunk \"(True,40,\\\"11111111 11111111 11111111 11111111 00001111\\\")\"]]" (DocTest.asPrint( test (minBound::Int32) )),  DocTest.test "src/Flat/Instances/Base.hs:466" "[ExpectedLine [LineChunk \"(True,40,\\\"11111110 11111111 11111111 11111111 00001111\\\")\"]]" (DocTest.asPrint( test (maxBound::Int32) )),  DocTest.test "src/Flat/Instances/Base.hs:477" "[ExpectedLine [LineChunk \"(True,8,\\\"00000000\\\")\"]]" (DocTest.asPrint( test (0::Int64) )),  DocTest.test "src/Flat/Instances/Base.hs:480" "[ExpectedLine [LineChunk \"(True,80,\\\"11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 00000001\\\")\"]]" (DocTest.asPrint( test (minBound::Int64) )),  DocTest.test "src/Flat/Instances/Base.hs:483" "[ExpectedLine [LineChunk \"(True,80,\\\"11111110 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 00000001\\\")\"]]" (DocTest.asPrint( test (maxBound::Int64) )),  DocTest.test "src/Flat/Instances/Base.hs:502" "[ExpectedLine [LineChunk \"(True,32,\\\"00000000 00000000 00000000 00000000\\\")\"]]" (DocTest.asPrint( test (0::Float) )),  DocTest.test "src/Flat/Instances/Base.hs:505" "[ExpectedLine [LineChunk \"(True,32,\\\"00000000 00000000 00000000 00000001\\\")\"]]" (DocTest.asPrint( test (1.4012984643E-45::Float) )),  DocTest.test "src/Flat/Instances/Base.hs:508" "[ExpectedLine [LineChunk \"(True,32,\\\"00000000 01111111 11111111 11111111\\\")\"]]" (DocTest.asPrint( test (1.1754942107E-38::Float) )),  DocTest.test "src/Flat/Instances/Base.hs:535" "[ExpectedLine [LineChunk \"(True,16,\\\"00000100 00000010\\\")\"]]" (DocTest.asPrint( test (4 :+ 2 :: Complex Word8) )),  DocTest.test "src/Flat/Instances/Base.hs:543" "[ExpectedLine [LineChunk \"(True,16,\\\"00000011 00000100\\\")\"]]" (DocTest.asPrint( test (3%4::Ratio Word8) )),  DocTest.test "src/Flat/Instances/Base.hs:555" "[ExpectedLine [LineChunk \"(True,1,\\\"0\\\")\"]]" (DocTest.asPrint( test ([]::[Bool]) )),  DocTest.test "src/Flat/Instances/Base.hs:558" "[ExpectedLine [LineChunk \"(True,5,\\\"10100\\\")\"]]" (DocTest.asPrint( test [False,False] )),  DocTest.test "src/Flat/Instances/Base.hs:592" "[ExpectedLine [LineChunk \"(True,2,\\\"10\\\")\"]]" (DocTest.asPrint( test (B.fromList [True]) )),  DocTest.test "src/Flat/Instances/Base.hs:595" "[ExpectedLine [LineChunk \"(True,4,\\\"0100\\\")\"]]" (DocTest.asPrint( test (B.fromList [False,False]) )),  DocTest.test "src/Flat/Instances/Base.hs:605" "[ExpectedLine [LineChunk \"(True,1,\\\"0\\\")\"]]" (DocTest.asPrint( test (False,()) )),  DocTest.test "src/Flat/Instances/Base.hs:608" "[ExpectedLine [LineChunk \"(True,0,\\\"\\\")\"]]" (DocTest.asPrint( test ((),()) )),  DocTest.test "src/Flat/Instances/Base.hs:613" "[ExpectedLine [LineChunk \"(True,7,\\\"0111011\\\")\"]]" (DocTest.asPrint( test (False,True,True,True,False,True,True) ))]

max_stars_repo_name: lagunoff/flat

max_issues_repo_name: tittoassini/flat

max_forks_repo_name: tittoassini/flat


text:
{-# LANGUAGE FlexibleContexts, RankNTypes #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
{-# OPTIONS_GHC -fno-warn-type-defaults      #-}
module Main where
import Control.Monad.ST (runST)
import qualified Data.Vector.Unboxed as U
import Control.Monad (replicateM)
import System.Random.MWC
import Statistics.Distribution
import Statistics.Distribution.Uniform
import Statistics.Distribution.Normal
import Statistics.Distribution.Exponential
import Graphics.Rendering.Chart
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Diagrams
import Data.Default.Class
import Control.Lens

plotNumbers xs t = toRenderable $ do
  layout_title .= t
  plot $ fmap histToPlot $ liftEC $ do
    plot_hist_bins .= 15
    plot_hist_values .= xs
    plot_hist_norm_func .= const id

-- lineare Kongruenzmethode:
createRandom a b m = random
  where random x = let y = ((a * x) + b) `mod` m in (y:random y)

main = do
    putStrLn $ show $ take 100 $ createRandom 4 1 9 5

    -- erzeuge gleichverteilten Zufallszahlengenerator
    -- dieser Generator erzeugt immer die gleiche Folge an Zufallszahlen!
    -- Das ist hilfreich um reproduzierbare Simulationen zu haben
    -- Um unterschiedliche Zuffalszahlen, bei jedem Lauf zu erhalten,
    -- sollte create mit createSystemRandom ersetzt werden
    mwc <- create
    n_uniform <- replicateM 1000 $ genContVar (uniformDistr 0 5) mwc
    renderableToFile def "histogramm_uniform.svg" (plotNumbers n_uniform "Standardverteilung mit a = 0 und b = 5")
    n_exp <- replicateM 1000 $ genContVar (exponential 0.6) mwc
    renderableToFile def "histogramm_exp.svg" (plotNumbers n_exp "Exponentialverteilung mit λ = 0.6")
    n_normal <- replicateM 1000 $ genContVar (normalDistr 0.9 1) mwc
    renderableToFile def "histogramm_normal.svg" (plotNumbers n_normal "Gaussverteilung mit μ = 0.9 und σ = 1")

max_stars_repo_name: Zortaniac/simulation

max_issues_repo_name: Zortaniac/simulation

max_forks_repo_name: Zortaniac/simulation


text:
import Data.Vector
import Statistics.Sample

import Statistics.Distribution.Normal
import Statistics.Distribution.Poisson
import qualified Statistics.Distribution as S

s1 :: Vector Double
s1 = fromList [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

s2 :: PoissonDistribution
s2 = poisson 2.5

s3 :: NormalDistribution
s3 = normalDistr mean stdDev
  where
    mean   = 1
    stdDev = 1

descriptive :: IO ()
descriptive = do
  print $ range s1
  -- 9.0
  print $ mean s1
  -- 5.5
  print $ stdDev s1
  -- 3.0276503540974917
  print $ variance s1
  -- 8.25
  print $ harmonicMean s1
  -- 3.414171521474055
  print $ geometricMean s1
  -- 4.5287286881167645

discrete :: IO ()
discrete = do
  print $ S.cumulative s2 0
  -- 8.208499862389884e-2
  print $ S.mean s2
  -- 2.5
  print $ S.variance s2
  -- 2.5
  print $ S.stdDev s2
  -- 1.5811388300841898

continuous :: IO ()
continuous = do
  print $ S.cumulative s3 0
  -- 0.15865525393145707
  print $ S.quantile s3 0.5
  -- 1.0
  print $ S.density s3 0
  -- 0.24197072451914334
  print $ S.mean s3
  -- 1.0
  print $ S.variance s3
  -- 1.0
  print $ S.stdDev s3
  -- 1.0

max_stars_repo_name: aavogt/wiwinwlh

max_issues_repo_name: aavogt/wiwinwlh

max_forks_repo_name: aavogt/wiwinwlh


text:
{-# LANGUAGE DeriveDataTypeable, FlexibleContexts, Rank2Types #-}
{-# OPTIONS_HADDOCK hide #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Matrix.STBase
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : asinerimental
--

module Numeric.LinearAlgebra.Matrix.STBase
    where

import Control.Monad( forM_, when )
import Control.Monad.ST( ST, RealWorld, runST, unsafeInterleaveST,
    unsafeIOToST )
import Data.Maybe( fromMaybe )
import Data.Typeable( Typeable )
import Foreign( Ptr, advancePtr, peek, peekElemOff, pokeElemOff,
    mallocForeignPtrArray )
import Text.Printf( printf )
import Unsafe.Coerce( unsafeCoerce )

import Numeric.LinearAlgebra.Types
import qualified Foreign.BLAS as BLAS
import Numeric.LinearAlgebra.Matrix.Base hiding ( unsafeWith,
    unsafeToForeignPtr, unsafeFromForeignPtr, )
import qualified Numeric.LinearAlgebra.Matrix.Base as M
import Numeric.LinearAlgebra.Vector( STVector, RVector )
import qualified Numeric.LinearAlgebra.Vector as V

-- | Mutable dense matrices in the 'ST' monad.
newtype STMatrix s e = STMatrix { unSTMatrix :: Matrix e }
    deriving (Typeable)

-- | Mutable dense matrices in the 'IO' monad.
type IOMatrix = STMatrix RealWorld

-- | A safe way to create and work with a mutable matrix before returning
-- an immutable matrix for later perusal. This function avoids copying
-- the matrix before returning it - it uses 'unsafeFreeze' internally,
-- but this wrapper is a safe interface to that function.
create :: (Storable e) => (forall s . ST s (STMatrix s e)) -> Matrix e
create mx = runST $ mx >>= unsafeFreeze
{-# INLINE create #-}

-- | Converts a mutable matrix to an immutable one by taking a complete
-- copy of it.
freeze :: (RMatrix m, Storable e) => m e -> ST s (Matrix e)
freeze a = do
    a' <- newCopy a
    unsafeFreeze a'
{-# INLINE freeze #-}

-- | Read-only matrices
class RMatrix m where
    -- | Get the dimensions of the matrix (number of rows and columns).
    getDim :: (Storable e) => m e -> ST s (Int,Int)

    -- | Same as 'withCol' but does not range-check index.
    unsafeWithCol :: (Storable e)
                  => m e
                  -> Int
                  -> (forall v. RVector v => v e -> ST s a)
                  -> ST s a

    -- | Perform an action with a list of views of the matrix columns.
    withCols :: (Storable e)
                 => m e
                 -> (forall v . RVector v => [v e] -> ST s a)
                 -> ST s a

    -- | Same as 'withSlice' but does not range-check index.
    unsafeWithSlice :: (Storable e)
                    => (Int,Int)
                    -> (Int,Int)
                    -> m e
                    -> (forall m'. RMatrix m' => m' e -> ST s a)
                    -> ST s a

    -- | Possibly view a matrix as a vector and perform an action on the
    -- view.  This only succeeds if the matrix is stored contiguously in
    -- memory, i.e. if the matrix contains a single column or the \"lda\"
    -- of the matrix is equal to the number of rows.
    maybeWithVector :: (Storable e)
                    => m e
                    -> (forall v . RVector v => v e -> ST s a)
                    -> Maybe (ST s a)

    -- | Converts a read-only matrix into an immutable matrix. This simply
    -- casts the matrix from one type to the other without copying.
    -- Note that because the matrix is possibly not copied, any subsequent
    -- modifications made to the read-only version of the matrix may be shared
    -- with the immutable version. It is safe to use, therefore, if the
    -- read-only version is never modified after the freeze operation.
    unsafeFreeze :: (Storable e) => m e -> ST s (Matrix e)

    -- | Unsafe cast from a read-only matrix to a mutable matrix.
    unsafeThaw :: (Storable e)
               => m e -> ST s (STMatrix s e)

    -- | Execute an 'IO' action with a pointer to the first element in the
    -- matrix and the leading dimension (lda).
    unsafeWith :: (Storable e) => m e -> (Ptr e -> Int -> IO a) -> IO a

instance RMatrix Matrix where
    getDim = return . dim
    {-# INLINE getDim #-}
    unsafeWithCol a j f = f (unsafeCol a j)
    {-# INLINE unsafeWithCol #-}
    withCols a f = f (cols a)
    {-# INLINE withCols #-}
    unsafeWithSlice ij mn a f = f (unsafeSlice ij mn a)
    {-# INLINE unsafeWithSlice #-}
    maybeWithVector a f | isContig a = Just $ f (toVector a)
                        | otherwise = Nothing
    {-# INLINE maybeWithVector #-}
    unsafeWith = M.unsafeWith
    {-# INLINE unsafeWith #-}
    unsafeFreeze = return
    {-# INLINE unsafeFreeze #-}
    unsafeThaw = return . STMatrix
    {-# INLINE unsafeThaw #-}

instance RMatrix (STMatrix s) where
    getDim = return . dim . unSTMatrix
    {-# INLINE getDim #-}
    unsafeWithCol = unsafeWithCol . unSTMatrix
    {-# INLINE unsafeWithCol #-}
    withCols = withCols . unSTMatrix
    {-# INLINE withCols #-}
    unsafeWithSlice ij mn = unsafeWithSlice ij mn . unSTMatrix
    {-# INLINE unsafeWithSlice #-}
    maybeWithVector = maybeWithVector . unSTMatrix
    {-# INLINE maybeWithVector #-}
    unsafeWith = unsafeWith . unSTMatrix
    {-# INLINE unsafeWith #-}
    unsafeFreeze = return . unSTMatrix
    {-# INLINE unsafeFreeze #-}
    unsafeThaw v = return $ cast v
      where
        cast :: STMatrix s e -> STMatrix s' e
        cast = unsafeCoerce
    {-# INLINE unsafeThaw #-}

-- | Perform an action with a view of a mutable matrix column
-- (no index checking).
unsafeWithColM :: (Storable e)
               => STMatrix s e
               -> Int
               -> (STVector s e -> ST s a)
               -> ST s a
unsafeWithColM a j f =
    unsafeWithCol a j $ \c -> do
        mc <- V.unsafeThaw c
        f mc
{-# INLINE unsafeWithColM #-}

-- | Perform an action with a list of views of the mutable matrix columns. See
-- also 'withCols'.
withColsM :: (Storable e)
               => STMatrix s e
               -> ([STVector s e] -> ST s a)
               -> ST s a
withColsM a f =
    withCols a $ \cs -> do
        mcs <- thawVecs cs
        f mcs
  where
    thawVecs [] = return []
    thawVecs (c:cs) = unsafeInterleaveST $ do
        mc <- V.unsafeThaw c
        mcs <- thawVecs cs
        return $ mc:mcs
{-# INLINE withColsM #-}

-- | Possibly view a matrix as a vector and perform an action on the
-- view.  This succeeds when the matrix is stored contiguously in memory,
-- i.e. if the matrix contains a single column or the \"lda\" of the matrix
-- is equal to the number of rows.  See also 'maybeWithVector'.
maybeWithVectorM :: (Storable e)
                 => STMatrix s e
                 -> (STVector s e -> ST s a)
                 -> Maybe (ST s a)
maybeWithVectorM a f =
    maybeWithVector a $ \v -> do
        mv <- V.unsafeThaw v
        f mv
{-# INLINE maybeWithVectorM #-}

-- | View a vector as a matrix of the given shape and pass it to
-- the specified function.
withFromVector :: (RVector v, Storable e)
                   => (Int,Int)
                   -> v e
                   -> (forall m . RMatrix m => m e -> ST s a)
                   -> ST s a
withFromVector mn@(m,n) v f = do
    nv <- V.getDim v
    when (nv /= m*n) $ error $
        printf ("withFromVector (%d,%d) <vector with dim %d>:"
                ++ " dimension mismatch") m n nv
    iv <- V.unsafeFreeze v
    f $ fromVector mn iv
{-# INLINE withFromVector #-}

-- | View a mutable vector as a mutable matrix of the given shape and pass it
-- to the specified function.
withFromVectorM :: (Storable e)
                     => (Int,Int)
                     -> STVector s e
                     -> (STMatrix s e -> ST s a)
                     -> ST s a
withFromVectorM mn@(m,n) v f = do
    nv <- V.getDim v
    when (nv /= m*n) $ error $
        printf ("withFromVectorM (%d,%d) <vector with dim %d>:"
                ++ " dimension mismatch") m n nv
    withFromVector mn v $ \a -> do
        ma <- unsafeThaw a
        f ma
{-# INLINE withFromVectorM #-}

-- | View a vector as a matrix with one column and pass it to
-- the specified function.
withFromCol :: (RVector v, Storable e)
                => v e
                -> (forall m . RMatrix m => m e -> ST s a)
                -> ST s a
withFromCol v f = do
    m <- V.getDim v
    withFromVector (m,1) v f
{-# INLINE withFromCol #-}

-- | View a mutable vector as a mutable matrix with one column and pass it to
-- the specified function.
withFromColM :: (Storable e)
                  => STVector s e
                  -> (STMatrix s e -> ST s a)
                  -> ST s a
withFromColM v f = do
    m <- V.getDim v
    withFromVectorM (m, 1) v f
{-# INLINE withFromColM #-}

-- | View a vector as a matrix with one row and pass it to
-- the specified function.
withFromRow :: (RVector v, Storable e)
                => v e
                -> (forall m . RMatrix m => m e -> ST s a)
                -> ST s a
withFromRow v f = do
    n <- V.getDim v
    withFromVector (1,n) v f
{-# INLINE withFromRow #-}

-- | View a mutable vector as a mutable matrix with one row and pass it to
-- the specified function.
withFromRowM :: (Storable e)
                  => STVector s e
                  -> (STMatrix s e -> ST s a)
                  -> ST s a
withFromRowM v f = do
    n <- V.getDim v
    withFromVectorM (1,n) v f
{-# INLINE withFromRowM #-}

-- | Perform an action with a view of a matrix column.
withCol :: (RMatrix m, Storable e)
            => m e
            -> Int
            -> (forall v . RVector v => v e -> ST s a)
            -> ST s a
withCol a j f = do
    (m,n) <- getDim a
    when (j < 0 || j >= n) $ error $
        printf ("withCol <matrix with dim (%d,%d)> %d:"
                ++ " index out of range") m n j

    unsafeWithCol a j f
{-# INLINE withCol #-}

-- | Like 'withCol', but perform the action with a mutable view.
withColM :: (Storable e)
         => STMatrix s e
         -> Int
         -> (STVector s e -> ST s a)
         -> ST s a
withColM a j f = do
    (m,n) <- getDim a
    when (j < 0 || j >= n) $ error $
        printf ("withColM <matrix with dim (%d,%d)> %d:"
                ++ " index out of range") m n j

    unsafeWithColM a j f
{-# INLINE withColM #-}

-- | Create a new matrix of given shape, but do not initialize the elements.
new_ :: (Storable e) => (Int,Int) -> ST s (STMatrix s e)
new_ (m,n)
    | m < 0 || n < 0 = error $
        printf "new_ (%d,%d): invalid dimensions" m n
    | otherwise = unsafeIOToST $ do
        f <- mallocForeignPtrArray (m*n)
        return $ STMatrix $ M.unsafeFromForeignPtr f 0 (m,n) (max 1 m)

-- | Create a matrix with every element initialized to the same value.
new :: (Storable e) => (Int,Int) -> e -> ST s (STMatrix s e)
new (m,n) e = do
    a <- new_ (m,n)
    setElems a $ replicate (m*n) e
    return a

-- | Creates a new matrix by copying another one.
newCopy :: (RMatrix m, Storable e) => m e -> ST s (STMatrix s e)
newCopy a = do
    mn <- getDim a
    b <- new_ mn
    unsafeCopyTo b a
    return b

-- | @copyTo dst src@ replaces the values in @dst@ with those in
-- source.  The operands must be the same shape.
copyTo :: (RMatrix m, Storable e) => STMatrix s e -> m e -> ST s ()
copyTo = checkOp2 "copyTo" unsafeCopyTo
{-# INLINE copyTo #-}

-- | Same as 'copyTo' but does not range-check indices.
unsafeCopyTo :: (RMatrix m, Storable e) => STMatrix s e -> m e -> ST s ()
unsafeCopyTo = vectorOp2 V.unsafeCopyTo
{-# INLINE unsafeCopyTo #-}

-- | Get the indices of the elements in the matrix, in column-major order.
getIndices :: (RMatrix m, Storable e) => m e -> ST s [(Int,Int)]
getIndices a = do
    (m,n) <- getDim a
    return $ [ (i,j) | j <- [ 0..n-1 ], i <- [ 0..m-1 ] ]

-- | Lazily get the elements of the matrix, in column-major order.
getElems :: (RMatrix m, Storable e) => m e -> ST s [e]
getElems a = case maybeWithVector a V.getElems of
    Just es -> es
    Nothing -> withCols a $ \xs ->
                   concat `fmap` mapM V.getElems xs

-- | Get the elements of the matrix, in column-major order.
getElems' :: (RMatrix m, Storable e) => m e -> ST s [e]
getElems' a = case maybeWithVector a V.getElems' of
    Just es -> es
    Nothing -> withCols a $ \xs ->
                   concat `fmap` mapM V.getElems' xs

-- | Lazily get the association list of the matrix, in column-major order.
getAssocs :: (RMatrix m, Storable e) => m e -> ST s [((Int,Int),e)]
getAssocs a = do
    is <- getIndices a
    es <- getElems a
    return $ zip is es

-- | Get the association list of the matrix, in column-major order.
getAssocs' :: (RMatrix m, Storable e) => m e -> ST s [((Int,Int),e)]
getAssocs' a = do
    is <- getIndices a
    es <- getElems' a
    return $ zip is es

-- | Set all of the values of the matrix from the elements in the list,
-- in column-major order.
setElems :: (Storable e) => STMatrix s e -> [e] -> ST s ()
setElems a es =
    case maybeWithVectorM a (`V.setElems` es) of
        Just st  -> st
        Nothing -> do
            (m,n) <- getDim a
            go m n 0 es
  where
    go _ n j [] | j == n = return ()
    go m n j [] | j < n = error $
        printf ("setElems <matrix with dim (%d,%d>"
                ++ "<list with length %d>: not enough elements)") m n (j*m)
    go m n j es' =
        let (es1', es2') = splitAt m es'
        in do
            withColM a j (`V.setElems` es1')
            go m n (j+1) es2'

-- | Set the given values in the matrix.  If an index is repeated twice,
-- the value is implementation-defined.
setAssocs :: (Storable e) => STMatrix s e -> [((Int,Int),e)] -> ST s ()
setAssocs a ies =
    sequence_ [ write a i e | (i,e) <- ies ]

-- | Same as 'setAssocs' but does not range-check indices.
unsafeSetAssocs :: (Storable e) => STMatrix s e -> [((Int,Int),e)] -> ST s ()
unsafeSetAssocs a ies =
    sequence_ [ unsafeWrite a i e | (i,e) <- ies ]

-- | Set the specified row of the matrix to the given vector.
setRow :: (RVector v, Storable e)
       => STMatrix s e -> Int -> v e -> ST s ()
setRow a i x = do
    (m,n) <- getDim a
    nx <- V.getDim x

    when (i < 0 || i >= m) $ error $
        printf ("setRow <matrix with dim (%d,%d)> %d:"
                ++ " index out of range") m n i
    when (nx /= n) $ error $
        printf ("setRow <matrix with dim (%d,%d)> _"
                ++ " <vector with dim %d>:"
                ++ " dimension mismatch") m n nx

    unsafeSetRow a i x
{-# INLINE setRow #-}

-- | Same as 'setRow' but does not range-check index or check
-- vector dimension.
unsafeSetRow :: (RVector v, Storable e)
             => STMatrix s e -> Int -> v e -> ST s ()
unsafeSetRow a i x = do
    jes <- V.getAssocs x
    sequence_ [ unsafeWrite a (i,j) e | (j,e) <- jes ]
{-# INLINE unsafeSetRow #-}

-- | Exchange corresponding elements in the given rows.
swapRows :: (BLAS1 e)
         => STMatrix s e -> Int -> Int -> ST s ()
swapRows a i1 i2 = do
    (m,n) <- getDim a
    when (i1 < 0 || i1 >= m || i2 < 0 || i2 >= m) $ error $
        printf ("swapRows <matrix with dim (%d,%d)> %d %d"
                ++ ": index out of range") m n i1 i2
    unsafeSwapRows a i1 i2

-- | Same as 'swapRows' but does not range-check indices.
unsafeSwapRows :: (BLAS1 e)
               => STMatrix s e -> Int -> Int -> ST s ()
unsafeSwapRows a i1 i2 = when (i1 /= i2) $ do
    (_,n) <- getDim a
    unsafeIOToST $
        unsafeWith a $ \pa lda ->
            let px = pa `advancePtr` i1
                py = pa `advancePtr` i2
                incx = lda
                incy = lda
            in
                BLAS.swap n px incx py incy

-- | Exchange corresponding elements in the given columns.
swapCols :: (BLAS1 e)
         => STMatrix s e -> Int -> Int -> ST s ()
swapCols a j1 j2 = do
    (m,n) <- getDim a
    when (j1 < 0 || j1 >= n || j2 < 0 || j2 >= n) $ error $
        printf ("swapCols <matrix with dim (%d,%d)> %d %d"
                ++ ": index out of range") m n j1 j2
    unsafeSwapCols a j1 j2

-- | Same as 'swapCols' but does not range-check indices.
unsafeSwapCols :: (BLAS1 e)
               => STMatrix s e -> Int -> Int -> ST s ()
unsafeSwapCols a j1 j2 = when (j1 /= j2) $ do
    (m,_) <- getDim a
    unsafeIOToST $
        unsafeWith a $ \pa lda ->
            let px = pa `advancePtr` (j1*lda)
                py = pa `advancePtr` (j2*lda)
                incx = 1
                incy = 1
            in
                BLAS.swap m px incx py incy

-- | Copy the specified row of the matrix to the vector.
rowTo :: (RMatrix m, Storable e)
         => STVector s e -> m e -> Int -> ST s ()
rowTo x a i = do
    (m,n) <- getDim a
    nx <- V.getDim x
    when (i < 0 || i >= m) $ error $
        printf ("rowTo"
               ++ " _"
               ++ " <matrix with dim (%d,%d)>"
               ++ " %d:"
               ++ ": index out of range"
               ) m n i
    when (nx /= n) $ error $
        printf ("rowTo"
                ++ " <vector with dim %d>"
                ++ " <matrix with dim (%d,%d)>"
                ++ " _"
                ++ ": dimension mismatch") nx m n

    unsafeRowTo x a i
{-# INLINE rowTo #-}

-- | Same as 'rowTo' but does not range-check index or check dimension.
unsafeRowTo :: (RMatrix m, Storable e)
            =>  STVector s e -> m e -> Int ->ST s ()
unsafeRowTo x a i = do
    (_,n) <- getDim a
    forM_ [ 0..n-1 ] $ \j -> do
        e <- unsafeRead a (i,j)
        V.unsafeWrite x j e
{-# INLINE unsafeRowTo #-}

-- | Set the diagonal of the matrix to the given vector.
setDiag :: (RVector v, Storable e)
        => STMatrix s e -> v e -> ST s ()
setDiag a x = do
    (m,n) <- getDim a
    nx <- V.getDim x
    let mn = min m n

    when (nx /= mn) $ error $
        printf ("setRow <matrix with dim (%d,%d)>"
                ++ " <vector with dim %d>:"
                ++ " dimension mismatch") m n nx

    unsafeSetDiag a x
{-# INLINE setDiag #-}

-- | Same as 'setDiag' but does not range-check index or check dimension.
unsafeSetDiag :: (RVector v, Storable e)
              => STMatrix s e -> v e -> ST s ()
unsafeSetDiag a x = do
    ies <- V.getAssocs x
    sequence_ [ unsafeWrite a (i,i) e | (i,e) <- ies ]
{-# INLINE unsafeSetDiag #-}

-- | Copy the diagonal of the matrix to the vector.
diagTo :: (RMatrix m, Storable e)
       => STVector s e -> m e -> ST s ()
diagTo x a = do
    nx <- V.getDim x
    (m,n) <- getDim a
    let mn = min m n

    when (nx /= mn) $ error $
        printf ("diagTo"
                ++ " <vector with dim %d>"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch") nx m n

    unsafeDiagTo x a
{-# INLINE diagTo #-}

-- | Same as 'diagTo' but does not range-check index or check dimensions.
unsafeDiagTo :: (RMatrix m, Storable e)
             => STVector s e -> m e -> ST s ()
unsafeDiagTo x a = do
    (m,n) <- getDim a
    let mn = min m n
    forM_ [ 0..mn-1 ] $ \i -> do
        e <- unsafeRead a (i,i)
        V.unsafeWrite x i e
{-# INLINE unsafeDiagTo #-}

-- | Get the element stored at the given index.
read :: (RMatrix m, Storable e) => m e -> (Int,Int) -> ST s e
read a (i,j) = do
    (m,n) <- getDim a
    when (i < 0 || i >= m || j < 0 || j >= n) $ error $
        printf ("read <matrix with dim (%d,%d)> (%d,%d):"
                ++ " index out of range") m n i j
    unsafeRead a (i,j)
{-# INLINE read #-}

-- | Same as 'read' but does not range-check index.
unsafeRead :: (RMatrix m, Storable e) => m e -> (Int,Int) -> ST s e
unsafeRead a (i,j) = unsafeIOToST $
    unsafeWith a $ \p lda ->
        peekElemOff p (i + j * lda)
{-# INLINE unsafeRead #-}

-- | Set the element stored at the given index.
write :: (Storable e)
      => STMatrix s e -> (Int,Int) -> e -> ST s ()
write a (i,j) e = do
    (m,n) <- getDim a
    when (i < 0 || i >= m || j < 0 || j >= n) $ error $
        printf ("write <matrix with dim (%d,%d)> (%d,%d):"
                ++ " index out of range") m n i j
    unsafeWrite a (i,j) e

{-# INLINE write #-}

-- | Same as 'write' but does not range-check index.
unsafeWrite :: (Storable e)
            => STMatrix s e -> (Int,Int) -> e -> ST s ()
unsafeWrite a (i,j) e = unsafeIOToST $
    unsafeWith a $ \p lda ->
        pokeElemOff p (i + j * lda) e
{-# INLINE unsafeWrite #-}

-- | Modify the element stored at the given index.
modify :: (Storable e)
       => STMatrix s e -> (Int,Int) -> (e -> e) -> ST s ()
modify a (i,j) f = do
    (m,n) <- getDim a
    when (i < 0 || i >= m || j < 0 || j >= n) $ error $
        printf ("modify <matrix with dim (%d,%d)> (%d,%d):"
                ++ " index out of range") m n i j
    unsafeModify a (i,j) f
{-# INLINE modify #-}

-- | Same as 'modify' but does not range-check index.
unsafeModify :: (Storable e)
             => STMatrix s e -> (Int,Int) -> (e -> e) -> ST s ()
unsafeModify a (i,j) f = unsafeIOToST $
    unsafeWith a $ \p lda ->
        let o = i + j * lda
        in do
            e <- peekElemOff p o
            pokeElemOff p o $ f e
{-# INLINE unsafeModify #-}

-- | @mapTo dst f src@ replaces @dst@ elementwise with @f(src)@.
mapTo :: (RMatrix m, Storable e, Storable f)
      => STMatrix s f
      -> (e -> f)
      -> m e
      -> ST s ()
mapTo dst f src = (checkOp2 "mapTo _" $ \z x -> unsafeMapTo z f x) dst src
{-# INLINE mapTo #-}

-- | Same as 'mapTo' but does not check dimensions.
unsafeMapTo :: (RMatrix m, Storable e, Storable f)
            => STMatrix s f
            -> (e -> f)
            -> m e
            -> ST s ()
unsafeMapTo dst f src =
    fromMaybe colwise $ maybeWithVectorM dst $ \vdst ->
    fromMaybe colwise $ maybeWithVector  src $ \vsrc ->
        V.unsafeMapTo vdst f vsrc
  where
    colwise = withColsM dst $ \zs ->
              withCols   src $ \xs ->
                  sequence_ [ V.unsafeMapTo z f x
                            | (z,x) <- zip zs xs
                            ]

-- | @zipWithTo dst f x y@ replaces @dst@ elementwise with @f(x, y)@.
zipWithTo :: (RMatrix m1, RMatrix m2, Storable e1, Storable e2, Storable f)
          => STMatrix s f
          -> (e1 -> e2 -> f)
          -> m1 e1
          -> m2 e2
          -> ST s ()
zipWithTo dst f x y =
    (checkOp3 "zipWithTo _" $ \dst1 x1 y1 -> unsafeZipWithTo dst1 f x1 y1)
        dst x y
{-# INLINE zipWithTo #-}

-- | Same as 'zipWithTo' but does not check dimensions.
unsafeZipWithTo :: (RMatrix m1, RMatrix m2, Storable e1, Storable e2, Storable f)
                => STMatrix s f
                -> (e1 -> e2 -> f)
                -> m1 e1
                -> m2 e2
                -> ST s ()
unsafeZipWithTo dst f x y =
    fromMaybe colwise $ maybeWithVectorM dst $ \vdst ->
    fromMaybe colwise $ maybeWithVector    x $ \vx ->
    fromMaybe colwise $ maybeWithVector    y $ \vy ->
        V.unsafeZipWithTo vdst f vx vy
  where
    colwise = withColsM dst $ \vdsts ->
              withCols   x $ \vxs ->
              withCols   y $ \vys ->

                  sequence_ [ V.unsafeZipWithTo vdst f vx vy
                            | (vdst,vx,vy) <- zip3 vdsts vxs vys
                            ]

-- | Set every element in the matrix to a default value.  For
-- standard numeric types (including 'Double', 'Complex Double', and 'Int'),
-- the default value is '0'.
clear :: (Storable e) => STMatrix s e -> ST s ()
clear a = fromMaybe colwise $ maybeWithVectorM a V.clear
  where
    colwise = withColsM a $ mapM_ V.clear

-- | @withSlice (i,j) (m,n) a@ performs an action with a view of the
-- submatrix of @a@ starting at index @(i,j)@ and having dimension @(m,n)@.
withSlice :: (RMatrix m, Storable e)
          => (Int,Int)
          -> (Int,Int)
          -> m e
          -> (forall m'. RMatrix m' => m' e -> ST s a)
          -> ST s a
withSlice ij mn a f = do
    ia <- unsafeFreeze a
    f $ slice ij mn ia

-- | Like 'withSlice', but perform the action with a mutable view.
withSliceM :: (Storable e)
           => (Int,Int)
           -> (Int,Int)
           -> STMatrix s e
           -> (STMatrix s e -> ST s a)
           -> ST s a
withSliceM ij mn a f =
    withSlice ij mn a $ \a' -> do
        ma <- unsafeThaw a'
        f ma

-- | Perform an action with a view gotten from taking the given number of
-- rows from the start of the matrix.
withTakeRows :: (RMatrix m, Storable e)
             => Int
             -> m e
             -> (forall m'. RMatrix m' => m' e -> ST s a)
             -> ST s a
withTakeRows i a f = do
    ia <- unsafeFreeze a
    f $ takeRows i ia

-- | Like 'withTakeRows', but perform the action with a mutable view.
withTakeRowsM :: (Storable e)
              => Int
              -> STMatrix s e
              -> (STMatrix s e -> ST s a)
              -> ST s a
withTakeRowsM i a f =
    withTakeRows i a $ \a' -> do
        ma <- unsafeThaw a'
        f ma

-- | Perform an action with a view gotten from dropping the given number of
-- rows from the start of the matrix.
withDropRows :: (RMatrix m, Storable e)
             => Int
             -> m e
             -> (forall m'. RMatrix m' => m' e -> ST s a)
             -> ST s a
withDropRows n a f = do
    ia <- unsafeFreeze a
    f $ dropRows n ia

-- | Like 'withDropRows', but perform the action with a mutable view.
withDropRowsM :: (Storable e)
              => Int
              -> STMatrix s e
              -> (STMatrix s e -> ST s a)
              -> ST s a
withDropRowsM i a f =
    withDropRows i a $ \a' -> do
        ma <- unsafeThaw a'
        f ma

-- | Perform an action with views from splitting the matrix rows at the given
-- index.
withSplitRowsAt :: (RMatrix m, Storable e)
                => Int
                -> m e
                -> (forall m1 m2. (RMatrix m1, RMatrix m2) => m1 e -> m2 e -> ST s a)
                -> ST s a
withSplitRowsAt i a f = do
    ia <- unsafeFreeze a
    uncurry f $ splitRowsAt i ia

-- | Like 'withSplitRowsAt', but perform the action with a mutable view.
withSplitRowsAtM :: (Storable e)
                 => Int
                 -> STMatrix s e
                 -> (STMatrix s e -> STMatrix s e -> ST s a)
                 -> ST s a
withSplitRowsAtM i a f =
    withSplitRowsAt i a $ \a1' a2' -> do
        ma1 <- unsafeThaw a1'
        ma2 <- unsafeThaw a2'
        f ma1 ma2

-- | Perform an action with a view gotten from taking the given number of
-- columns from the start of the matrix.
withTakeCols :: (RMatrix m, Storable e)
             => Int
             -> m e
             -> (forall m'. RMatrix m' => m' e -> ST s a)
             -> ST s a
withTakeCols i a f = do
    ia <- unsafeFreeze a
    f $ takeCols i ia

-- | Like 'withTakeCols', but perform the action with a mutable view.
withTakeColsM :: (Storable e)
              => Int
              -> STMatrix s e
              -> (STMatrix s e -> ST s a)
              -> ST s a
withTakeColsM i a f =
    withTakeCols i a $ \a' -> do
        ma <- unsafeThaw a'
        f ma

-- | Perform an action with a view gotten from dropping the given number of
-- columns from the start of the matrix.
withDropCols :: (RMatrix m, Storable e)
             => Int
             -> m e
             -> (forall m'. RMatrix m' => m' e -> ST s a)
             -> ST s a
withDropCols n a f = do
    ia <- unsafeFreeze a
    f $ dropCols n ia

-- | Like 'withDropCols', but perform the action with a mutable view.
withDropColsM :: (Storable e)
              => Int
              -> STMatrix s e
              -> (STMatrix s e -> ST s a)
              -> ST s a
withDropColsM i a f =
    withDropCols i a $ \a' -> do
        ma <- unsafeThaw a'
        f ma

-- | Perform an action with views from splitting the matrix columns at the given
-- index.
withSplitColsAt :: (RMatrix m, Storable e)
                => Int
                -> m e
                -> (forall m1 m2. (RMatrix m1, RMatrix m2) => m1 e -> m2 e -> ST s a)
                -> ST s a
withSplitColsAt i a f = do
    ia <- unsafeFreeze a
    uncurry f $ splitColsAt i ia

-- | Like 'withSplitColsAt', but perform the action with mutable views.
withSplitColsAtM :: (Storable e)
                 => Int
                 -> STMatrix s e
                 -> (STMatrix s e -> STMatrix s e -> ST s a)
                 -> ST s a
withSplitColsAtM i a f =
    withSplitColsAt i a $ \a1' a2' -> do
        ma1 <- unsafeThaw a1'
        ma2 <- unsafeThaw a2'
        f ma1 ma2

-- | Add a vector to the diagonal of a matrix.
shiftDiagM_ :: (RVector v, BLAS1 e)
              => v e -> STMatrix s e -> ST s ()
shiftDiagM_ s a = do
    (m,n) <- getDim a
    ns <- V.getDim s
    let mn = min m n

    when (ns /= mn) $ error $
        printf ("shiftDiagM_"
                ++ " <vector with dim %d>"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch")
                ns
                m n

    shiftDiagWithScaleM_ 1 s a

-- | Add a scaled vector to the diagonal of a matrix.
shiftDiagWithScaleM_ :: (RVector v, BLAS1 e)
                       => e -> v e -> STMatrix s e -> ST s ()
shiftDiagWithScaleM_ e s a = do
    (m,n) <- getDim a
    ns <- V.getDim s
    let mn = min m n

    when (ns /= mn) $ error $
        printf ("shiftDiagWithScaleM_"
                ++ " _"
                ++ " <vector with dim %d>"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch")
                ns
                m n

    unsafeIOToST $
        V.unsafeWith s $ \ps ->
        unsafeWith a $ \pa lda ->
            BLAS.axpy mn e ps 1 pa (lda+1)

-- | Add two matrices.
addTo :: (RMatrix m1, RMatrix m2, VNum e)
      =>  STMatrix s e -> m1 e -> m2 e -> ST s ()
addTo = checkOp3 "addTo" $ vectorOp3 V.addTo

-- | Subtract two matrices.
subTo :: (RMatrix m1, RMatrix m2, VNum e)
      => STMatrix s e -> m1 e -> m2 e -> ST s ()
subTo = checkOp3 "subTo" $ vectorOp3 V.subTo

-- | Conjugate the entries of a matrix.
conjugateTo :: (RMatrix m, VNum e)
            => STMatrix s e -> m e -> ST s ()
conjugateTo = checkOp2 "conjugateTo" $
    vectorOp2 V.conjugateTo

-- | Negate the entries of a matrix.
negateTo :: (RMatrix m, VNum e)
         => STMatrix s e -> m e -> ST s ()
negateTo = checkOp2 "negateTo" $
    vectorOp2 V.negateTo

-- | Scale the entries of a matrix by the given value.
scaleM_ :: (BLAS1 e)
          => e -> STMatrix s e -> ST s ()
scaleM_ e = vectorOp (V.scaleM_ e)

-- | @addWithScaleM_ alpha x y@ sets @y := alpha * x + y@.
addWithScaleM_ :: (RMatrix m, BLAS1 e)
               => e -> m e -> STMatrix s e -> ST s ()
addWithScaleM_ e = checkOp2 "addWithScaleM_" $
    unsafeAddWithScaleM_ e

unsafeAddWithScaleM_ :: (RMatrix m, BLAS1 e)
                     => e -> m e -> STMatrix s e -> ST s ()
unsafeAddWithScaleM_ alpha x y =
    fromMaybe colwise $ maybeWithVector  x $ \vx ->
    fromMaybe colwise $ maybeWithVectorM y $ \vy ->
        V.unsafeAddWithScaleM_ alpha vx vy
  where
    colwise = withCols   x $ \vxs ->
              withColsM y $ \vys ->
                  sequence_ [ V.unsafeAddWithScaleM_ alpha vx vy
                            | (vx,vy) <- zip vxs vys ]

-- | Scale the rows of a matrix; @scaleRowsM_ s a@ sets
-- @a := diag(s) * a@.
scaleRowsM_ :: (RVector v, BLAS1 e)
              => v e -> STMatrix s e -> ST s ()
scaleRowsM_  s a = do
    (m,n) <- getDim a
    ns <- V.getDim s
    when (ns /= m) $ error $
        printf ("scaleRowsM_"
                ++ " <vector with dim %d>"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch")
                ns
                m n

    unsafeIOToST $
        V.unsafeWith s $ \ps ->
        unsafeWith a   $ \pa lda ->
            go m n lda pa ps 0
  where
    go m n lda pa ps i | i == m    = return ()
                       | otherwise = do
                             e <- peek ps
                             BLAS.scal n e pa lda
                             go m n lda (pa `advancePtr` 1)
                                        (ps `advancePtr` 1)
                                        (i+1)

-- | Scale the columns of a matrix; @scaleColBysM_ s a@ sets
-- @a := a * diag(s)@.
scaleColsM_ :: (RVector v, BLAS1 e)
            => v e -> STMatrix s e -> ST s ()
scaleColsM_ s a = do
    (m,n) <- getDim a
    ns <- V.getDim s
    when (ns /= n) $ error $
        printf ("scaleColsM_"
                ++ " <vector with dim %d>"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch")
                ns
                m n

    es <- V.getElems s
    withColsM a $ \xs ->
        sequence_ [ V.scaleM_ e x
                  | (e,x) <- zip es xs
                  ]

-- | @rank1UpdateM_ alpha x y a@ sets @a := alpha * x * y^H + a@.
rank1UpdateM_ :: (RVector v1, RVector v2, BLAS2 e)
              => e -> v1 e -> v2 e -> STMatrix s e -> ST s ()
rank1UpdateM_ alpha x y a = do
    (m,n) <- getDim a
    nx <- V.getDim x
    ny <- V.getDim y

    when (nx /= m || ny /= n) $ error $
        printf ("rank1UpdateTo"
                ++ " _"
                ++ " <vector with dim %d>"
                ++ " <vector with dim %d>"
                ++ ": dimension mismatch"
                ++ "<matrix with dim (%d,%d)>"                )
                nx
                ny
                m n

    unsafeIOToST $
        V.unsafeWith x $ \px ->
        V.unsafeWith y $ \py ->
        unsafeWith a $ \pa lda ->
            BLAS.gerc m n alpha px 1 py 1 pa lda

-- | @transTo dst a@ sets @dst := trans(a)@.
transTo :: (RMatrix m, BLAS1 e)
        => STMatrix s e
        -> m e
        -> ST s ()
transTo a' a = do
    (ma,na) <- getDim a
    (ma',na') <- getDim a'
    let (m,n) = (ma,na)

    when ((ma,na) /= (na',ma')) $ error $
        printf ( "transTo"
               ++ " <matrix with dim (%d,%d)>"
               ++ " <matrix with dim (%d,%d)>"
               ++ ": dimension mismatch"
               )
               ma' na'
               ma na

    unsafeIOToST $
        unsafeWith a' $ \pa' lda' ->
        unsafeWith a $ \pa lda -> let
            go j px py | j == n = return ()
                       | otherwise = do
                           BLAS.copy m px 1 py lda'
                           go (j+1) (px `advancePtr` lda) (py `advancePtr` 1)
            in go 0 pa pa'

-- | @conjTransTo dst a@ sets @dst := conjugate(trans(a))@.
conjTransTo :: (RMatrix m, BLAS1 e)
            => STMatrix s e
            -> m e
            -> ST s ()
conjTransTo a' a = do
    transTo a' a
    conjugateTo a' a'

-- | @mulVectorTo dst transa a x@
-- sets @dst := op(a) * x@, where @op(a)@ is determined by @transa@.
mulVectorTo :: (RMatrix m, RVector v, BLAS2 e)
            => STVector s e
            -> Trans -> m e
            -> v e
            -> ST s ()
mulVectorTo dst = mulVectorWithScaleTo dst 1

-- | @mulVectorWithScaleTo dst alpha transa a x@
-- sets @dst := alpha * op(a) * x@, where @op(a)@ is determined by @transa@.
mulVectorWithScaleTo :: (RMatrix m, RVector v, BLAS2 e)
                     => STVector s e
                     -> e
                     -> Trans -> m e
                     -> v e
                     -> ST s ()
mulVectorWithScaleTo dst alpha t a x =
    addMulVectorWithScalesM_ alpha t a x 0 dst

-- | @addMulVectorWithScalesM_ alpha transa a x beta y@
-- sets @y := alpha * op(a) * x + beta * y@, where @op(a)@ is
-- determined by @transa@.
addMulVectorWithScalesM_ :: (RMatrix m, RVector v, BLAS2 e)
                         => e
                         -> Trans -> m e
                         -> v e
                         -> e
                         -> STVector s e
                         -> ST s ()
addMulVectorWithScalesM_ alpha transa a x beta y = do
    (ma,na) <- getDim a
    nx <- V.getDim x
    ny <- V.getDim y
    let (m,n) = (ny,nx)

    when ((not . and) [ case transa of NoTrans -> (ma,na) == (m,n)
                                       _       -> (ma,na) == (n,m)
                      , nx == n
                      , ny == m
                      ]) $ error $
        printf ("addMulVectorWithScalesTo"
                ++ " _"
                ++ " %s"
                ++ " <matrix with dim (%d,%d)>"
                ++ " <vector with dim %d>"
                ++ " _"
                ++ " <vector with dim %d>"
                ++ ": dimension mismatch")
               (show transa)
               ma na
               nx
               ny

    unsafeIOToST $
        unsafeWith a $ \pa lda ->
        V.unsafeWith x $ \px ->
        V.unsafeWith y $ \py ->
            if n == 0
                then BLAS.scal m beta py 1
                else BLAS.gemv transa ma na alpha pa lda px 1 beta py 1

-- | @mulMatrixTo dst transa a transb b@
-- sets @dst := op(a) * op(b)@, where @op(a)@ and @op(b)@ are determined
-- by @transa@ and @transb@.
mulMatrixTo :: (RMatrix m1, RMatrix m2, BLAS3 e)
            => STMatrix s e
            -> Trans -> m1 e
            -> Trans -> m2 e
            -> ST s ()
mulMatrixTo dst = mulMatrixWithScaleTo dst 1

-- | @mulMatrixWithScaleTo alpha transa a transb b c@
-- sets @c := alpha * op(a) * op(b)@, where @op(a)@ and @op(b)@ are determined
-- by @transa@ and @transb@.
mulMatrixWithScaleTo :: (RMatrix m1, RMatrix m2, BLAS3 e)
                     => STMatrix s e
                     -> e
                     -> Trans -> m1 e
                     -> Trans -> m2 e
                     -> ST s ()
mulMatrixWithScaleTo dst alpha ta a tb b =
    addMulMatrixWithScalesM_ alpha ta a tb b 0 dst

-- | @addMulMatrixWithScalesM_ alpha transa a transb b beta c@
-- sets @c := alpha * op(a) * op(b) + beta * c@, where @op(a)@ and
-- @op(b)@ are determined by @transa@ and @transb@.
addMulMatrixWithScalesM_ :: (RMatrix m1, RMatrix m2, BLAS3 e)
                         => e
                         -> Trans -> m1 e
                         -> Trans -> m2 e
                         -> e
                         -> STMatrix s e
                         -> ST s ()
addMulMatrixWithScalesM_ alpha transa a transb b beta c = do
    (ma,na) <- getDim a
    (mb,nb) <- getDim b
    (mc,nc) <- getDim c
    let (m,n) = (mc,nc)
        k = case transa of NoTrans -> na
                           _       -> ma

    when ((not . and) [ case transa of NoTrans -> (ma,na) == (m,k)
                                       _       -> (ma,na) == (k,m)
                      , case transb of NoTrans -> (mb,nb) == (k,n)
                                       _       -> (mb,nb) == (n,k)
                      , (mc, nc) == (m,n)
                      ]) $ error $
        printf ("addMulMatrixWithScalesM_"
                ++ " _"
                ++ " %s <matrix with dim (%d,%d)>"
                ++ " %s <matrix with dim (%d,%d)>"
                ++ " _"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch")
               (show transa) ma na
               (show transb) mb nb
               mc nc

    unsafeIOToST $
        unsafeWith a $ \pa lda ->
        unsafeWith b $ \pb ldb ->
        unsafeWith c $ \pc ldc ->
            BLAS.gemm transa transb m n k alpha pa lda pb ldb beta pc ldc

checkOp2 :: (RMatrix x, RMatrix y, Storable e, Storable f)
         => String
         -> (x e -> y f -> ST s a)
         -> x e
         -> y f
         -> ST s a
checkOp2 str f x y = do
    (m1,n1) <- getDim x
    (m2,n2) <- getDim y
    when ((m1,n1) /= (m2,n2)) $ error $
        printf ("%s <matrix with dim (%d,%d)> <matrix with dim (%d,%d)>:"
                ++ " dimension mismatch") str m1 n1 m2 n2
    f x y
{-# INLINE checkOp2 #-}

checkOp3 :: (RMatrix x, RMatrix y, RMatrix z, Storable e, Storable f, Storable g)
         => String
         -> (x e -> y f -> z g -> ST s a)
         -> x e
         -> y f
         -> z g
         -> ST s a
checkOp3 str f x y z = do
    (m1,n1) <- getDim x
    (m2,n2) <- getDim y
    (m3,n3) <- getDim z
    when((m1,n1) /= (m2,n2) || (m1,n1) /= (m3,n3)) $ error $
        printf ("%s <matrix with dim (%d,%d)> <matrix with dim (%d,%d)>:"
                ++ " <matrix with dim (%d,%d)> dimension mismatch")
               str m1 n1 m2 n2 m3 n3
    f x y z
{-# INLINE checkOp3 #-}

vectorOp :: (Storable e)
         => (STVector s e -> ST s ())
         -> STMatrix s e -> ST s ()
vectorOp f x =
    fromMaybe colwise $ maybeWithVectorM x $ \vx -> f vx
  where
    colwise = withColsM x $ \vxs ->
                  sequence_ [ f vx | vx <- vxs ]

vectorOp2 :: (RMatrix m, Storable e, Storable f)
          => (forall v . RVector v => STVector s f -> v e -> ST s ())
          -> STMatrix s f -> m e -> ST s ()
vectorOp2 f dst x =
    fromMaybe colwise $ maybeWithVectorM dst $ \vdst ->
    fromMaybe colwise $ maybeWithVector    x $ \vx ->
        f vdst vx
  where
    colwise = withColsM dst $ \vdsts ->
              withCols   x   $ \vxs ->
                  sequence_ [ f vdst vx | (vdst,vx) <- zip vdsts vxs ]
{-# INLINE vectorOp2 #-}

vectorOp3 :: (RMatrix m1, RMatrix m2, Storable e1, Storable e2, Storable f)
          => (forall v1 v2 . (RVector v1, RVector v2) =>
                  STVector s f -> v1 e1 -> v2 e2 -> ST s ())
          -> STMatrix s f -> m1 e1 -> m2 e2 -> ST s ()
vectorOp3 f dst x y =
    fromMaybe colwise $ maybeWithVectorM dst $ \vdst ->
    fromMaybe colwise $ maybeWithVector    x $ \vx ->
    fromMaybe colwise $ maybeWithVector    y $ \vy ->
        f vdst vx vy
  where
    colwise = withColsM dst $ \vdsts ->
              withCols   x   $ \vxs ->
              withCols   y   $ \vys ->
                  sequence_ [ f vdst vx vy
                            | (vdst,vx,vy) <- zip3 vdsts vxs vys ]
{-# INLINE vectorOp3 #-}

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
module PostgreSQL.Binary.Prelude
(
  module Exports,
  LazyByteString,
  ByteStringBuilder,
  LazyText,
  TextBuilder,
  mapLeft,
  joinMap,
)
where

-- base
-------------------------
import Control.Applicative as Exports
import Control.Arrow as Exports hiding (first, second)
import Control.Category as Exports
import Control.Concurrent as Exports
import Control.Exception as Exports
import Control.Monad as Exports hiding (fail, mapM_, sequence_, forM_, msum, mapM, sequence, forM)
import Control.Monad.IO.Class as Exports
import Control.Monad.Fail as Exports
import Control.Monad.Fix as Exports hiding (fix)
import Control.Monad.ST as Exports
import Data.Bifunctor as Exports
import Data.Bits as Exports
import Data.Bool as Exports
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
import Data.Function as Exports hiding (id, (.))
import Data.Functor as Exports
import Data.Functor.Identity as Exports
import Data.Int as Exports
import Data.IORef as Exports
import Data.Ix as Exports
import Data.List as Exports hiding (sortOn, isSubsequenceOf, uncons, concat, foldr, foldl1, maximum, minimum, product, sum, all, and, any, concatMap, elem, foldl, foldr1, notElem, or, find, maximumBy, minimumBy, mapAccumL, mapAccumR, foldl')
import Data.List.NonEmpty as Exports (NonEmpty(..))
import Data.Maybe as Exports
import Data.Monoid as Exports hiding (Last(..), First(..), (<>))
import Data.Ord as Exports
import Data.Proxy as Exports
import Data.Ratio as Exports
import Data.Semigroup as Exports
import Data.STRef as Exports
import Data.String as Exports
import Data.Traversable as Exports
import Data.Tuple as Exports
import Data.Unique as Exports
import Data.Version as Exports
import Data.Void as Exports
import Data.Word as Exports
import Debug.Trace as Exports
import Foreign.ForeignPtr as Exports
import Foreign.Ptr as Exports
import Foreign.StablePtr as Exports
import Foreign.Storable as Exports hiding (sizeOf, alignment)
import GHC.Conc as Exports hiding (orElse, withMVar, threadWaitWriteSTM, threadWaitWrite, threadWaitReadSTM, threadWaitRead)
import GHC.Exts as Exports (lazy, inline, sortWith, groupWith, IsList(fromList, Item))
import GHC.Generics as Exports (Generic, Generic1)
import GHC.IO.Exception as Exports
import Numeric as Exports
import Prelude as Exports hiding (fail, concat, foldr, mapM_, sequence_, foldl1, maximum, minimum, product, sum, all, and, any, concatMap, elem, foldl, foldr1, notElem, or, mapM, sequence, id, (.))
import System.Environment as Exports
import System.Exit as Exports
import System.IO as Exports
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, ReadS, readP_to_S, readS_to_P)
import Text.ParserCombinators.ReadPrec as Exports (ReadPrec, readPrec_to_P, readP_to_Prec, readPrec_to_S, readS_to_Prec)
import Text.Printf as Exports (printf, hPrintf)
import Text.Read as Exports (Read(..), readMaybe, readEither)
import Unsafe.Coerce as Exports

-- transformers
-------------------------
import Control.Monad.Trans.State.Strict as Exports hiding (liftCallCC, liftCatch)
import Control.Monad.Trans.Reader as Exports hiding (liftCallCC, liftCatch)
import Control.Monad.Trans.Class as Exports
import Data.Functor.Identity as Exports

-- bytestring
-------------------------
import Data.ByteString as Exports (ByteString)

-- text
-------------------------
import Data.Text as Exports (Text)

-- vector
-------------------------
import Data.Vector as Exports (Vector)

-- scientific
-------------------------
import Data.Scientific as Exports (Scientific)

-- uuid
-------------------------
import Data.UUID as Exports (UUID)

-- time
-------------------------
import Data.Time as Exports

-- unordered-containers
-------------------------
import Data.HashMap.Strict as Exports (HashMap)

-- containers
-------------------------
import Data.Map.Strict as Exports (Map)

-- custom
-------------------------
import qualified Data.ByteString.Lazy
import qualified Data.ByteString.Builder
import qualified Data.Text.Lazy
import qualified Data.Text.Lazy.Builder

type LazyByteString =
  Data.ByteString.Lazy.ByteString

type ByteStringBuilder =
  Data.ByteString.Builder.Builder

type LazyText =
  Data.Text.Lazy.Text

type TextBuilder =
  Data.Text.Lazy.Builder.Builder

{-# INLINE mapLeft #-}
mapLeft :: (a -> b) -> Either a x -> Either b x
mapLeft f =
  either (Left . f) Right

joinMap :: Monad m => (a -> m b) -> m a -> m b
joinMap f =
  join . liftM f

max_stars_repo_name: nikita-volkov/postgresql-binary

max_issues_repo_name: nikita-volkov/postgresql-binary

max_forks_repo_name: nikita-volkov/postgresql-binary


text:
module Raytracer.Light where

import Numeric.LinearAlgebra.Data (Vector, (|>), toList)
import Numeric.LinearAlgebra ((<.>))
import Raytracer.Geometry (Pos, rayTo, Intersectable, Ray(Ray))
import Raytracer.Camera (fire_ray)
import Data.Maybe (isJust)

import Codec.Picture (PixelRGB8(PixelRGB8), colorMap)

-- This has the position of the light, the cutoff range of the light, and the colour
data Light a = Light Pos Double a deriving (Show)

--computeLighting :: (Intersectable a) => a b -> [Light b] -> Pos -> [b]
computeLighting mesh lights pos = fmap attenuate $ filter (not . visible) $ fmap makeRay lights
  where
  makeRay light@(Light lPos _ _) = (pos `rayTo` lPos, light)
  visible (ray, _) = isJust $ fire_ray mesh $ ray

attenuate (Ray ray _, Light _ cutoff a) = if dist > cutoff then PixelRGB8 0 0 0 else blend
  where
  dist = sqrt $ ray <.> ray
  blend = colorMap (\x -> round $ (fromIntegral x) * (cutoff - dist) / cutoff) a


text:
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE MultiWayIf #-}

module Isumi.Math.FFT.Internal.Bruun
  ( fftDirect
  , cruN
  , fftBruun
  , fftBruunC
  )
where

import           Data.Bits                          (shiftL)
import           Data.Complex
import           Data.HashMap.Strict                (HashMap)
import qualified Data.HashMap.Strict                as HashMap
import           Data.List                          (concatMap)
import           Data.Maybe
import qualified Data.Vector                        as V
import qualified Data.Vector.Unboxed                as UV
import           Math.NumberTheory.Logarithms       (integerLog2)

import           Isumi.Math.FFT.Internal.Polynomial
import           Isumi.Math.FFT.Internal.Utils

-- | Computes fft using definition.
fftDirect :: UV.Vector (Complex Double) -> UV.Vector (Complex Double)
fftDirect xs = UV.generate len genEntry
 where
  len = UV.length xs
  genEntry k = UV.sum . flip UV.imap xs $ \n xn ->
    (crus UV.! k)^n * xn
  crus = cruN len

{-|
  Computes the complex roots of unity
-}
cruN :: Int -> UV.Vector (Complex Double)
cruN n = UV.generate n $ \k ->
  exp (realToComplex k * (-2) * pi * compI / fromIntegral n)

-- | FFT with Bruun's algorithm, for real input
fftBruun :: UV.Vector Double -> Maybe (UV.Vector (Complex Double))
fftBruun = fftBruunG RealPolynomial

-- | FFT with Bruun's alogrithm, for complex input
fftBruunC :: UV.Vector (Complex Double) -> Maybe (UV.Vector (Complex Double))
fftBruunC = fftBruunG ComplexPolynomial

fftBruunG :: (UV.Unbox a, Floating a, Eq a)
          => (Polynomial a -> RealOrComplexPolynomial Double)
          -> UV.Vector a
          -> Maybe (UV.Vector (Complex Double))
fftBruunG f xs =
  if | (not . isPowerOf2) n -> Nothing
     | UV.all (== 0) xs -> Just $ UV.replicate n 0
     | otherwise ->
       let r0 = f . fromJust . polyFromVecRep . UV.reverse $ xs
           d0 = BruunDivisorMinus n
        in Just $ fftBruun' (cruN n) (d0, r0)
  where
  n = UV.length xs

fftBruun' :: UV.Vector (Complex Double) -- ^cru array
          -> DivisorRemainder -- ^initial divisor and remainder
          -> UV.Vector (Complex Double)
fftBruun' cruValues dr =
  UV.generate n (\i -> cruToFFT HashMap.! (cruValues' V.! i))
  where
  n = UV.length cruValues
  cruToFFT :: HashMap (Integer, Integer) (Complex Double)
  cruToFFT = HashMap.fromList
           . fmap (\(d, r) -> (normalizeCRU . finalDivisorToCRU $ d,
               finalRemainderToValue r))
           . factorDRAll $ dr
  cruValues' = fmap normalizeCRU (UV.convert cruValues)
    :: V.Vector (Integer, Integer)

normalizeCRU :: Complex Double -> (Integer, Integer)
normalizeCRU (r :+ i) = (toInt r, toInt i)
  where
  toInt x = round (x * 1000000)

finalDivisorToCRU :: BruunDivisor Double -> Complex Double
finalDivisorToCRU (BruunDivisorComplex c) = (-1) * c
finalDivisorToCRU (BruunDivisorMinus n)   = if n == 1 then 1 else error ""
finalDivisorToCRU (BruunDivisorPlus n c)  =
  if n == 1 && c == 0
     then (-1)
     else error ""

finalRemainderToValue :: RealOrComplexPolynomial Double -> Complex Double
finalRemainderToValue (ComplexPolynomial p) =
  case p of
    PolyVecRep values ->
      if UV.length values == 1
         then values UV.! 0
         else error $ "Final stage remainder MUST have a length of 1. Actual: "
                      ++ show (UV.length values)
    PolyExpCoefRep _ ->
      error "Final stage remainder MUST have exp-coef representation"
finalRemainderToValue (RealPolynomial p) =
  finalRemainderToValue . ComplexPolynomial . mapPoly realToComplex $ p

-- | Recursively divisor-remainder pairs until they can not be factored
factorDRAll :: DivisorRemainder -> [DivisorRemainder]
factorDRAll dr = go [dr]
  where
  go :: [DivisorRemainder] -> [DivisorRemainder]
  go xs = let ys = concatMap (extractListFromMaybe . factorDR) xs
           in if length ys /= 2 * length xs then xs else go ys

type DivisorRemainder = (BruunDivisor Double, RealOrComplexPolynomial Double)

-- | Factor a divisor-remainder pair into TWO sub divisor remainder pairs.
factorDR :: DivisorRemainder -> Maybe [DivisorRemainder]
factorDR (d, r) =
  if divDegree d == 1
     then Nothing
     else fmap (fromJust . flip updateDR r) . pairToList <$> factorDivisor d

updateDR :: BruunDivisor Double -- ^new divisor
         -> RealOrComplexPolynomial Double -- ^old remainder
         -> Maybe DivisorRemainder -- ^(new divisor, new remainder)
updateDR d r = (d,) <$> r `modPolyG` divToPoly d

pairToList :: (a, a) -> [a]
pairToList (x, y) = [x, y]

isPowerOf2 :: Int -> Bool
isPowerOf2 x = (1 `shiftL` integerLog2 (fromIntegral x)) == x

extractListFromMaybe :: Maybe [a] -> [a]
extractListFromMaybe Nothing   = []
extractListFromMaybe (Just xs) = xs

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun


text:
{-# language TypeFamilies #-}
module Bases where

import qualified Data.Vector.Generic as G
import Data.Complex

haarV, haarU :: (G.Vector v a, a ~ Complex Double) => Int -> v a
haarV dim = G.fromList $ [1 / sqrt 2, -1 / sqrt 2] ++ take (dim - 2) (repeat 0)
haarU dim = G.fromList $ [1 / sqrt 2,  1 / sqrt 2] ++ take (dim - 2) (repeat 0)

max_stars_repo_name: RubenAstudillo/proyecto-lena

max_issues_repo_name: RubenAstudillo/proyecto-lena

max_forks_repo_name: RubenAstudillo/proyecto-lena


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

module Examples.MultiTargetTracking where

import Control.Arrow (returnA)
import Control.Monad (forM)
import Control.Monad.Bayes.Class (MonadSample, MonadInfer)

import qualified Data.Map as M
import Data.Maybe (fromJust)

import Numeric.LinearAlgebra.Static hiding ((<>))

import Inference (zdsparticles)
import DelayedSampling (DelayedSample, DelayedInfer, Result (..))
import qualified SymbolicDistr as DS
import DSProg (DeepForce (..), Expr' (..), Expr, marginal, zdeepForce, deepForce', forgetE)
import Distributions
import Util.ZStream (ZStream)
import qualified Util.ZStream as ZS
import Metaprob ((~~), Gen, lift, (|->), obs)
import qualified Metaprob as MP
import Examples.Demo (Sampler, Delayed, Weighted)

data ObsType = Clutter | NewTrack | Track Int
  deriving (Eq, Show, Ord)

filterM :: Monad m => (k -> v -> m Bool) -> M.Map k v -> m (M.Map k v)
filterM f = M.traverseMaybeWithKey $ \k v -> do
  b <- f k v
  return (if b then Just v else Nothing)

type TrackMap = M.Map Int (Expr (R 4))

tdiff, birthRate, deathRate, clutterLambda, newTrackLambda, pd, survivalProb :: Double
tdiff = 1
birthRate = 0.1
deathRate = 0.02
clutterLambda = 3
newTrackLambda = birthRate * tdiff
pd = 0.8
survivalProb = exp (- tdiff * deathRate)

clutterDistr :: DS.Distr (Expr (R 2))
clutterDistr = DS.mvNormal (Const (0 :: R 2)) (10 * sym eye)

newTrackD :: DelayedSample m => m (Expr (R 4))
newTrackD = DS.sample (DS.mvNormal (Const mu) cov)
  where
  mu = 0 :: R 4
  cov = sym (((5 * eye :: Sq 2) ||| (0 :: Sq 2))
            ===
            ((0 :: Sq 2) ||| (0.1 * eye :: Sq 2)))

trackMotion :: DelayedSample m => Double -> Expr (R 4) -> m (Expr (R 4))
trackMotion tdiff track = DS.sample (DS.mvNormal (MVMul (Const motionMatrix) track) motionCov)
  where
  motionMatrix :: Sq 4
  motionMatrix = eye + konst tdiff *
    (((0 :: Sq 2) ||| (eye :: Sq 2))
    ===
    ((- 1 / 100 * eye :: Sq 2) ||| (-0.1 * eye :: Sq 2)))
  motionCov :: Sym 4
  motionCov = sym $ konst tdiff * (((0.01 * eye :: Sq 2) ||| (0 :: Sq 2))
                                  ===
                                  ((0 :: Sq 2) ||| (0.1 * eye :: Sq 2)))

trackMeasurement :: Expr (R 4) -> DS.Distr (Expr (R 2))
trackMeasurement posvel = DS.mvNormal (MVMul (Const posFromPosVel) posvel) (sym eye)
  where
  posFromPosVel :: L 2 4
  posFromPosVel = (eye :: Sq 2) ||| (0 :: Sq 2)

updateWithAssocs :: DelayedSample m => Int -> TrackMap -> [ObsType] -> m ((TrackMap, Int), [DS.Distr (Expr (R 2))])
updateWithAssocs nextTrackID updatedOldTracks assocs = do
  (obsDists, newTrackPVs) <- fmap mconcat . forM assocs $ \k -> case k of
    Clutter -> return ([clutterDistr], [])
    NewTrack -> do
      newTrackPV <- newTrackD
      return ([trackMeasurement newTrackPV], [newTrackPV])
    Track i -> return ([trackMeasurement (updatedOldTracks M.! i)], [])
  coasted <- filterM (\_ _ -> sample (bernoulli ((1 - pd) / (1 - survivalProb * pd)))) notObserved
  return ((mconcat (observed : coasted : zipWith M.singleton [nextTrackID..] newTrackPVs),
           nextTrackID + length newTrackPVs), obsDists)
  where
  (observed, notObserved) = M.partitionWithKey (\i _ -> Track i `elem` assocs) updatedOldTracks

sampleStep :: DelayedSample m => (TrackMap, Int) -> Gen m ((TrackMap, Int), [Expr (R 2)])
sampleStep (allOldTracks, nextTrackID) = do
  updatedOldTracks <- lift (mapM (trackMotion tdiff) allOldTracks)
  assocs <- "assocs" ~~ MP.prim (shuffleWithRepeats' countDistrs)
  (newTrackInfo, obsDists) <- lift (updateWithAssocs nextTrackID updatedOldTracks assocs)
  observations <- "observations" ~~ MP.isequence (map MP.dsPrim obsDists)
  return (newTrackInfo, observations)
  where
  countDistrs = M.singleton Clutter (poisson clutterLambda)
             <> M.singleton NewTrack (poisson newTrackLambda)
             <> M.mapKeys Track (fmap (\_ -> bernoulli01 (survivalProb * pd)) allOldTracks)

observeStep :: DelayedInfer m => (TrackMap, Int) -> [R 2] -> m ((TrackMap, Int), [Expr (R 2)])
observeStep (allOldTracks, nextTrackID) observations = do
  updatedOldTracks <- mapM (trackMotion tdiff) allOldTracks
  -- assocs <- proposeAssocs observations updatedOldTracks
  MP.observingWithProposal
    ("observations" |-> MP.trList observations) (sampleStep (allOldTracks, nextTrackID))
    "assocs" (proposeAssocs 0 observations updatedOldTracks)

proposeAssocs :: DelayedInfer m => Int -> [R 2] -> TrackMap -> Gen m [ObsType]
proposeAssocs j (obs : observations) remainingTracks = do
  newTrack <- lift $ newTrackD
  let distrs = M.singleton Clutter (clutterLambda, clutterDistr)
        <> M.singleton NewTrack (newTrackLambda, trackMeasurement newTrack)
        <> M.mapKeys Track (fmap (\t -> (survivalProb * pd, trackMeasurement t)) remainingTracks)
  likes <- lift $ mapM (\(intensity, d) -> (\ll -> intensity * exp ll) <$> DS.score d obs) distrs
  let assocDistr = let probs = fmap (/ sum likes) likes in categoricalM probs
  i <- ("assoc" ++ show j) ~~ MP.prim assocDistr
  let remainingTracks' = case i of
        Track k -> M.delete k remainingTracks
        _ -> remainingTracks
  (i :) <$> proposeAssocs (j + 1) observations remainingTracks'
proposeAssocs _ [] remainingTracks = return []

generateGroundTruth :: DelayedSample m => ZStream m () (TrackMap, [Expr (R 2)])
generateGroundTruth = ZS.fromStep stepf initState where
  initState = (mempty, 0)
  stepf state () = do
    (newState@(tracks', _), observations) <- MP.sim (sampleStep state)
    return (newState, (tracks', observations))

processObservations :: DelayedInfer m => ZStream m [R 2] (M.Map Int (Result (R 4)))
processObservations = ZS.fromStep stepf initState where
  initState = (mempty, 0)
  stepf state obs = do
    (newState@(tracks', _), _) <- observeStep state obs
    marginalTracks <- mapM (fmap fromJust . marginal) tracks'
    return (newState, marginalTracks)

runMTTPF :: Int -> ZStream Sampler () ([(Int, R 4)], [[(Int, Result (R 4))]], [R 2])
runMTTPF numParticles = proc () -> do
  (groundTruth, obs) <- zdeepForce generateGroundTruth -< ()
  particles <- zdsparticles numParticles processObservations -< obs
  returnA -< (M.assocs groundTruth, map M.assocs particles, obs)

max_stars_repo_name: psg-mit/probzelus-haskell

max_issues_repo_name: psg-mit/probzelus-haskell

max_forks_repo_name: psg-mit/probzelus-haskell


text:
import Numeric.LinearAlgebra
import ActivationFunction
import Mnist
import SampleWeight

batchSize = 100

predict :: SampleWeight -> Vector R -> Vector R
predict ([w1,w2,w3],[b1,b2,b3]) x =
    softMax' . (\x'' -> sumInput x'' w3 b3) . sigmoid . (\x' -> sumInput x' w2 b2) . sigmoid $ sumInput x w1 b1

sumInput :: Vector R -> Weight -> Bias -> Vector R
sumInput x w b = (x <# w) + b

maxIndexPredict :: SampleWeight -> Vector R -> Double
maxIndexPredict sw x = fromIntegral . maxIndex $ predict sw x

take' :: Indexable c t => Int -> Int -> c -> [t]
take' n1 n2 x
    | n1 >= n2  = []
    | otherwise = (x ! n1) : take' (n1+1) n2 x

increment :: [Double] -> [Double] -> Double
increment ps l = fromIntegral . length . filter id $ zipWith (==) ps l

countAccuracy :: Double -> Int -> SampleWeight -> DataSet -> Double
countAccuracy a n sw ds@(i,l)
    | n <= 0    = a
    | otherwise =
      if maxIndexPredict sw (i ! (n-1)) == l ! (n-1)
        then countAccuracy (a+1) (n-1) sw ds
        else countAccuracy a (n-1) sw ds

-- For batch
countAccuracy' :: Double -> Int -> SampleWeight -> DataSet -> Double
countAccuracy' a n sw ds@(i,l)
    | n <= 0    = a
    | otherwise = countAccuracy' (a+cnt) (n-batchSize) sw ds
        where ps = maxIndexPredict sw <$> take' (n-batchSize) n i
              ls = take' (n-batchSize) n l
              cnt = increment ps ls

main = do
    [_, ds] <- loadMnist True
    sw <- loadSW
    let r = rows $ fst ds
        cnt = countAccuracy' 0 r sw ds

    putStrLn $ "Accuracy: " ++ show (cnt / fromIntegral r)


text:
module HVX.DcpTests.NoAffineNonmonConvexNondec where

import Numeric.LinearAlgebra

import HVX

main :: IO ()
main = do
  let x = EVar "x"
      _ = (EConst $ (1><2) [1,2]) *~ hexp x
  return ()

max_stars_repo_name: wellposed/hvx

max_issues_repo_name: wellposed/hvx

max_forks_repo_name: chrisnc/hvx


text:
-- | Interface to the FFTW routine
--
-- Ubuntu installation of C library:
--
--     sudo apt-get install libfftw3-dev
--
-- Basic usage:
-- <http://www.fftw.org/fftw3_doc/Complex-One_002dDimensional-DFTs.html#Complex-One_002dDimensional-DFTs>
--
-- Linker flags:
--
--     cc file.c -lm -lfftw3

module FFTW where

import qualified Prelude as P

import qualified Data.Complex as Complex

import qualified Test.QuickCheck as QC
import qualified Test.QuickCheck.Monadic as QC

import Language.Embedded.Expression (varExp)
  -- `varExp` can be used to make named constants in expressions

import Feldspar.Run
import Feldspar.Data.Vector

import DFT
import FFT_bench (printTime_def)

-- | Wrapper for the FFTW routine
--
-- This wrapper is mainly used for testing. It's not suitable for real code
-- because it constructs a plan every time it's called, and it silently
-- allocates an array for the output.
fftw :: DManifest (Complex Double) -> Run (DManifest (Complex Double))
fftw inp = do
    addInclude "<fftw3.h>"
    out  <- newArr (length inp)
    plan <- newObject "fftw_plan" False
    callProcAssign plan "fftw_plan_dft_1d"
      [ valArg (length inp)
      , iarrArg inp
      , arrArg out
      , valArg (varExp "FFTW_FORWARD"  :: Data Word32)
      , valArg (varExp "FFTW_ESTIMATE" :: Data Word32)
      ]
    callProc "fftw_execute" [objArg plan]
    freezeArr out

almostEq a b
    =    Complex.magnitude d P.< 1e-7
    P.&& Complex.phase d     P.< 1e-7
  where
    d = abs (a-b)

a ~= b = P.and $ P.zipWith almostEq a b

prop_fft_dft dft' fft' = QC.monadicIO $ do
    n   :: Int              <- QC.pick $ QC.choose (2,5)
    inp :: [Complex Double] <- QC.pick $ QC.vector (2 P.^ n)
    outd <- QC.run $ dft' inp
    outf <- QC.run $ fft' inp
    QC.assert (outd ~= outf)

-- | Compare 'fftw' against 'dft'
testFFTW =
    marshalledM (return . dft) $ \dft'  ->
    marshalledM fftw           $ \fftw' ->
      QC.quickCheck $ prop_fft_dft dft' fftw'
  where
    marshalledM = marshalled' def def
      { externalFlagsPre  = ["-Wno-incompatible-pointer-types"]
      , externalFlagsPost = ["-lm -lfftw3"]
      }

sizeOf_fftw_complex :: Data Length
sizeOf_fftw_complex = 16
  -- Checked on an x86_64 system

-- | Measure the time for 100 runs of 'fftw' (excluding initialization) for
-- arrays of the given size
benchmark n = do
    addInclude "<stdio.h>"
    addInclude "<string.h>"
    addInclude "<time.h>"
    addInclude "<fftw3.h>"

    addDefinition printTime_def

    inp  <- newObject "fftw_complex" True
    out  <- newObject "fftw_complex" True
    plan <- newObject "fftw_plan" False

    callProcAssign inp "fftw_malloc" [valArg (n*sizeOf_fftw_complex)]
    callProcAssign out "fftw_malloc" [valArg (n*sizeOf_fftw_complex)]
    callProc "memset"
      [ objArg inp
      , valArg (0 :: Data Index)
      , valArg (n*sizeOf_fftw_complex)
      ]

    callProcAssign plan "fftw_plan_dft_1d"
      [ valArg (n :: Data Word32)
      , objArg inp
      , objArg out
      , valArg (varExp "FFTW_FORWARD"  :: Data Word32)
      , valArg (varExp "FFTW_ESTIMATE" :: Data Word32)
          -- Change to `FFTW_MEASURE` to enable tuning
      ]

    start <- newObject "clock_t" False
    end   <- newObject "clock_t" False
    callProcAssign start "clock" []

    for (0,1,Excl 100) $ \(_ :: Data Index) ->
      callProc "fftw_execute" [objArg plan]

    callProcAssign end "clock" []
    callProc "printTime" [objArg start, objArg end]

    callProc "fftw_destroy_plan" [objArg plan]
    callProc "fftw_free" [objArg inp]
    callProc "fftw_free" [objArg out]

runBenchmark n = runCompiled'
    def
    def {externalFlagsPre = ["-O3"], externalFlagsPost = ["-lm","-lfftw3"]}
    (benchmark n)

max_stars_repo_name: Abhiroop/mu-feldspar

max_issues_repo_name: Abhiroop/mu-feldspar

max_forks_repo_name: Abhiroop/mu-feldspar


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Vector.ST
-- Copyright  : Copyright (c) , Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Mutable vectors in the ST monad.

module Numeric.LinearAlgebra.Vector.ST (
    -- * Mutable vectors
    STVector,
    IOVector,
    create,

    -- * Read-only vectors
    RVector(..),

    -- * Creating new vectors
    new_,
    new,

    -- * Copying vectors
    newCopy,
    copyTo,
    swap,

    -- * Reading and writing vector elements
    read,
    write,
    modify,
    getIndices,
    getElems,
    getElems',
    getAssocs,
    getAssocs',
    setElems,
    setAssocs,
    clear,

    -- * List-like operations
    mapTo,
    zipWithTo,

    -- * Vector linear algebra
    getSumAbs,
    getNorm2,
    getWhichMaxAbs,
    getDot,
    scaleM_,
    addWithScaleM_,
    kroneckerTo,

    -- * Vector views
    withSlice,
    withDrop,
    withTake,
    withSplitAt,

    withSliceM,
    withDropM,
    withTakeM,
    withSplitAtM,

    -- * Vector math operations
    -- ** Num
    addTo,
    subTo,
    mulTo,
    negateTo,
    conjugateTo,
    absTo,
    signumTo,

    -- ** Fractional
    divTo,
    recipTo,

    -- ** Floating
    sqrtTo,
    expTo,
    logTo,
    powTo,
    sinTo,
    cosTo,
    tanTo,
    asinTo,
    acosTo,
    atanTo,
    sinhTo,
    coshTo,
    tanhTo,
    asinhTo,
    acoshTo,
    atanhTo,

    -- * Conversions between mutable and immutable vectors
    freeze,

    -- * Unsafe operations
    unsafeCopyTo,
    unsafeSwap,
    unsafeRead,
    unsafeWrite,
    unsafeModify,
    unsafeMapTo,
    unsafeZipWithTo,
    unsafeAddWithScaleM_,
    unsafeGetDot,

    ) where

import Prelude()
import Numeric.LinearAlgebra.Vector.STBase


text:
{-# LANGUAGE CPP #-}
-- |
-- Module    : Numeric.SpecFunctions
-- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Special functions and factorials.
module Numeric.SpecFunctions (
    -- * Error function
    erf
  , erfc
  , invErf
  , invErfc
    -- * Gamma function
  , logGamma
  , logGammaL
  , incompleteGamma
  , invIncompleteGamma
  , digamma
    -- * Beta function
  , logBeta
  , incompleteBeta
  , incompleteBeta_
  , invIncompleteBeta
    -- * Sinc
  , sinc
    -- * Logarithm
    -- $log1p
  , log1p
  , log1pmx
  , log2
    -- * Exponent
  , expm1
    -- * Factorial
  , factorial
  , logFactorial
  , stirlingError
    -- * Combinatorics
  , choose
  , logChoose
    -- * References
    -- $references
  ) where

import Numeric.SpecFunctions.Internal

-- $log1p
--
-- Base starting from @4.9.0@ (GHC 8.0) provides 'log1p' and 'expm1'
-- as method of class 'Floating'. In this case we simply reexport
-- these function. Otherwise we provide our own with more restrictive
-- signature @Double → Double@.

-- $references
--
-- * Bernardo, J. (1976) Algorithm AS 103: Psi (digamma)
--   function. /Journal of the Royal Statistical Society. Series C
--   (Applied Statistics)/ 25(3):315-317.
--   <http://www.jstor.org/stable/2347257>
--
-- * Cran, G.W., Martin, K.J., Thomas, G.E. (1977) Remark AS R19
--   and Algorithm AS 109: A Remark on Algorithms: AS 63: The
--   Incomplete Beta Integral AS 64: Inverse of the Incomplete Beta
--   Function Ratio. /Journal of the Royal Statistical Society. Series
--   C (Applied Statistics)/ Vol. 26, No. 1 (1977), pp. 111-114
--   <http://www.jstor.org/pss/2346887>
--
-- * Lanczos, C. (1964) A precision approximation of the gamma
--   function.  /SIAM Journal on Numerical Analysis B/
--   1:86–96. <http://www.jstor.org/stable/2949767>
--
-- * Loader, C. (2000) Fast and Accurate Computation of Binomial
--   Probabilities. <http://projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf>
--
-- * Macleod, A.J. (1989) Algorithm AS 245: A robust and reliable
--   algorithm for the logarithm of the gamma function.
--   /Journal of the Royal Statistical Society, Series C (Applied Statistics)/
--   38(2):397–402. <http://www.jstor.org/stable/2348078>
--
-- * Majumder, K.L., Bhattacharjee, G.P. (1973) Algorithm AS 63: The
--   Incomplete Beta Integral. /Journal of the Royal Statistical
--   Society. Series C (Applied Statistics)/ Vol. 22, No. 3 (1973),
--   pp. 409-411. <http://www.jstor.org/pss/2346797>
--
-- * Majumder, K.L., Bhattacharjee, G.P. (1973) Algorithm AS 64:
--   Inverse of the Incomplete Beta Function Ratio. /Journal of the
--   Royal Statistical Society. Series C (Applied Statistics)/
--   Vol. 22, No. 3 (1973), pp. 411-414
--   <http://www.jstor.org/pss/2346798>
--
-- * Temme, N.M. (1992) Asymptotic inversion of the incomplete beta
--   function. /Journal of Computational and Applied Mathematics
--   41(1992) 145-157.
--
-- * Temme, N.M. (1994) A set of algorithms for the incomplete gamma
--   functions. /Probability in the Engineering and Informational
--   Sciences/, 8, 1994, 291-307. Printed in the U.S.A.

max_stars_repo_name: Shimuuar/math-functions

max_issues_repo_name: Shimuuar/math-functions

max_forks_repo_name: Shimuuar/math-functions