text:
{-# LANGUAGE BangPatterns #-}
module Image.IO
  ( ImageRepa(..)
  , readImagePathList
  , imagePathSource
  , readImageRepa
  , readImageConduit
  , plotImageRepa
  , complexImageToColorImage
  , plotImageRepaComplex
  , plotImageRepaComplexGray
  ) where

import           Codec.Picture
import           Control.Monad.IO.Class       (liftIO)
import           Control.Monad.Trans.Resource
import           Data.Array.Repa              as R
import           Data.Complex
import           Data.Conduit                 as C
import           Data.Conduit.List            as CL
import           Data.List                    as L
import           Data.Vector.Unboxed          as VU
import           Data.Word
import           GHC.Float
import           Utils.Parallel

-- x = cols = width
-- y = rows = height
data ImageRepa a = ImageRepa
  { imageDepth   :: !Int
  , imageContent :: !(R.Array U DIM3 a) -- nf x cols x rows
  }

readImagePathList :: FilePath -> IO [String]
readImagePathList = fmap lines . readFile

imagePathSource :: FilePath -> ConduitT () FilePath (ResourceT IO) ()
imagePathSource filePath = do
  pathList <- liftIO $ readImagePathList filePath
  sourceList pathList

readImageRepa :: FilePath -> Bool -> IO (ImageRepa Double)
readImageRepa filePath isColor = do
  buffer <- liftIO $ readImage filePath
  case buffer of
    Left msg -> error msg
    Right dImg ->
      return $
      if isColor
        then case dImg of
               ImageY8 img ->
                 ImageRepa 8 . computeS $
                 fromFunction
                   (Z :. (1 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. _ :. x :. y) ->
                      fromIntegral $ pixelAt img x y :: Double)
               ImageY16 img ->
                 ImageRepa 16 . computeS $
                 fromFunction
                   (Z :. (1 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. _ :. x :. y) ->
                      fromIntegral $ pixelAt img x y :: Double)
               ImageYF img ->
                 ImageRepa 64 . computeS $
                 fromFunction
                   (Z :. (1 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. _ :. x :. y) ->
                      float2Double $ pixelAt img x y :: Double)
               ImageRGB8 img ->
                 ImageRepa 8 . computeS $
                 fromFunction
                   (Z :. (3 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. k :. x :. y) ->
                      let (PixelRGB8 r g b) = pixelAt img x y
                       in case k of
                            0 -> fromIntegral r
                            1 -> fromIntegral g
                            2 -> fromIntegral b
                            _ -> error "readImageConduit: dimension error.")
               ImageRGB16 img ->
                 ImageRepa 16 . computeS $
                 fromFunction
                   (Z :. (3 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. k :. x :. y) ->
                      let (PixelRGB16 r g b) = pixelAt img x y
                       in case k of
                            0 -> fromIntegral r
                            1 -> fromIntegral g
                            2 -> fromIntegral b
                            _ -> error "readImageConduit: dimension error.")
               ImageRGBF img ->
                 ImageRepa 64 . computeS $
                 fromFunction
                   (Z :. (3 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. k :. x :. y) ->
                      let (PixelRGBF r g b) = pixelAt img x y
                       in case k of
                            0 -> float2Double r
                            1 -> float2Double g
                            2 -> float2Double b
                            _ -> error "readImageConduit: dimension error.")
               img ->
                 let rgbImg = convertRGB8 img
                  in ImageRepa 8 . computeS $
                     fromFunction
                       (Z :. (3 :: Int) :. imageWidth rgbImg :.
                        imageHeight rgbImg)
                       (\(Z :. k :. x :. y) ->
                          let (PixelRGB8 r g b) = pixelAt rgbImg x y
                           in case k of
                                0 -> fromIntegral r
                                1 -> fromIntegral g
                                2 -> fromIntegral b
                                _ -> error "readImageConduit: dimension error.")
        else case dImg of
               ImageY8 img ->
                 ImageRepa 8 . computeS $
                 fromFunction
                   (Z :. (1 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. _ :. x :. y) ->
                      fromIntegral $ pixelAt img x y :: Double)
               ImageY16 img ->
                 ImageRepa 16 . computeS $
                 fromFunction
                   (Z :. (1 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. _ :. x :. y) ->
                      fromIntegral $ pixelAt img x y :: Double)
               ImageYF img ->
                 ImageRepa 64 . computeS $
                 fromFunction
                   (Z :. (1 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. _ :. x :. y) ->
                      float2Double $ pixelAt img x y :: Double)
               ImageRGB8 img ->
                 ImageRepa 8 . computeS $
                 fromFunction
                   (Z :. (1 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. _ :. x :. y) ->
                      let (PixelRGB8 r g b) = pixelAt img x y
                       in rgb2Gray
                            (fromIntegral (maxBound :: Word8))
                            (fromIntegral r)
                            (fromIntegral g)
                            (fromIntegral b))
               ImageRGB16 img ->
                 ImageRepa 16 . computeS $
                 fromFunction
                   (Z :. (1 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. _ :. x :. y) ->
                      let (PixelRGB16 r g b) = pixelAt img x y
                       in rgb2Gray
                            (fromIntegral (maxBound :: Word16))
                            (fromIntegral r)
                            (fromIntegral g)
                            (fromIntegral b))
               ImageRGBF img ->
                 ImageRepa 64 . computeS $
                 fromFunction
                   (Z :. (1 :: Int) :. imageWidth img :. imageHeight img)
                   (\(Z :. _ :. x :. y) ->
                      let (PixelRGBF r g b) = pixelAt img x y
                       in rgb2Gray
                            1
                            (float2Double r)
                            (float2Double g)
                            (float2Double b))
               img ->
                 let rgbImg = convertRGB8 img
                  in ImageRepa 8 . computeS $
                     fromFunction
                       (Z :. (1 :: Int) :. imageWidth rgbImg :.
                        imageHeight rgbImg)
                       (\(Z :. _ :. x :. y) ->
                          let (PixelRGB8 r g b) = pixelAt rgbImg x y
                           in rgb2Gray
                                (fromIntegral (maxBound :: Word8))
                                (fromIntegral r)
                                (fromIntegral g)
                                (fromIntegral b))

readImageConduit ::
     Bool -> ConduitT FilePath (ImageRepa Double) (ResourceT IO) ()
readImageConduit isColor =
  awaitForever
    (\filePath -> do
       img <- liftIO $ readImageRepa filePath isColor
       yield img)

{-# INLINE rgb2Gray #-}
rgb2Gray :: Double -> Double -> Double -> Double -> Double
rgb2Gray bound r g b
  | yLinear <= 0.0031308 = 12.92 * yLinear * bound
  | otherwise = (1.055 * (yLinear ** (1 / 2.4)) - 0.055) * bound
  where
    yLinear =
      0.2126 * func bound r + 0.7152 * func bound g + 0.0722 * func bound b

{-# INLINE func #-}
func :: Double -> Double -> Double
func bound x
  | y < 0.04045 = y / 12.92
  | otherwise = ((y + 0.055) / 1.055) ** 2.4
  where
    y = x / bound

{-# INLINE normalizeImageRepa #-}
normalizeImageRepa :: (RealFrac e, Eq e, Unbox e, Ord e) => ImageRepa e -> ImageRepa e
normalizeImageRepa i@(ImageRepa depth img)
  | maxV == minV =
    ImageRepa depth . computeS . R.map (\x -> fromIntegral (2 ^ depth - 1)) $
    img
  | otherwise =
    ImageRepa depth .
    computeS .
    R.map (\x -> (x - minV) / (maxV - minV) * fromIntegral (2 ^ depth - 1)) $
    img
  where
    !maxV = VU.maximum . toUnboxed $ img
    !minV = VU.minimum . toUnboxed $ img

{-# INLINE normalizeImageRepaPar #-}
normalizeImageRepaPar :: (RealFrac e, Eq e, Unbox e, Ord e) => ImageRepa e -> IO (ImageRepa e)
normalizeImageRepaPar i@(ImageRepa depth img) = do
  if maxV == minV
    then do
      arr <- computeP . R.map (\x -> fromIntegral (2 ^ depth - 1)) $img
      return . ImageRepa depth $ arr
    else do
      arr <-
        computeP .
        R.map (\x -> (x - minV) / (maxV - minV) * fromIntegral (2 ^ depth - 1)) $
        img
      return . ImageRepa depth $ arr
  where
    !maxV = VU.maximum . toUnboxed $ img
    !minV = VU.minimum . toUnboxed $ img

{-# INLINE plotImageRepa #-}
plotImageRepa :: (Unbox e, Ord e, RealFrac e) => FilePath -> ImageRepa e -> IO ()
plotImageRepa filePath img@(ImageRepa depth x) = do
  normalizedImg <- imageContent <$> normalizeImageRepaPar img
  let Z :. nfp' :. nxp' :. nyp' = extent x
      -- normalizedImg = imageContent . normalizeImageRepa $ img
      w =
        case nfp' of
          1
            -- ImageY8 $
            -- generateImage
            --   (\x y ->
            --      let v =
            --            fromIntegral . round $
            --            normalizedImg R.! (Z :. 0 :. x :. y)
            --       in v)
            --   nxp'
            --   nyp'
           ->
            ImageRGB8 $
            generateImage
              (\x y ->
                 let r =
                       fromIntegral . round $
                       normalizedImg R.! (Z :. 0 :. x :. y)
                     g =
                       fromIntegral . round $
                       normalizedImg R.! (Z :. 0 :. x :. y)
                     b =
                       fromIntegral . round $
                       normalizedImg R.! (Z :. 0 :. x :. y)
                  in PixelRGB8 r g b)
              nxp'
              nyp'
          3 ->
            ImageRGB8 $
            generateImage
              (\x y ->
                 let r =
                       fromIntegral . round $
                       normalizedImg R.! (Z :. 0 :. x :. y)
                     g =
                       fromIntegral . round $
                       normalizedImg R.! (Z :. 1 :. x :. y)
                     b =
                       fromIntegral . round $
                       normalizedImg R.! (Z :. 2 :. x :. y)
                  in PixelRGB8 r g b)
              nxp'
              nyp'
          _ ->
            error $
            "ImageRepa is neither a gray image nor a color image. There are " L.++
            show nfp' L.++
            " channels."
  savePngImage filePath w

-- Plot Complex Image
{-# INLINE complexScale #-}
complexScale :: (Ord e, Floating e, Unbox e) => ImageRepa (Complex e) -> e
complexScale (ImageRepa _ arr) =
  let (rVec, iVec) = VU.unzip . VU.map (\(a :+ b) -> (a, b)) . R.toUnboxed $ arr
      rMax = VU.maximum rVec
      rMin = VU.minimum rVec
      iMax = VU.maximum iVec
      iMin = VU.minimum iVec
   in 2 / ((max rMax iMax) - (min rMin iMin))

{-# INLINE complexImageToColorImage #-}
complexImageToColorImage ::
     (Floating e, Ord e, Unbox e) => ImageRepa (Complex e) -> ImageRepa e
complexImageToColorImage img@(ImageRepa depth arr) =
  let scale = complexScale img
      (Z :. _ :. cols :. rows) = extent arr
      (r, g, b) =
        VU.unzip3 .
        VU.map
          (\(r :+ i) ->
             let [x, y] = L.map (* scale) [r, i]
                 radius = sqrt (x * x + y * y)
                 a = 0.40824829046386301636 * x
                 b = 0.70710678118654752440 * y
                 d = 1 / (1 + radius * radius)
                 d' = 0.5 - radius * d
                 r' = 0.5 + (0.81649658092772603273 * x * d)
                 b' = 0.5 - (d * (a - b))
                 g' = 0.5 - (d * (a + b))
                 [red, grn, blu] = L.map (+ d') [r', g', b']
                 [red', grn', blu'] = L.map (flip (-) d') [r', g', b']
              in if radius < 1
                   then (red', grn', blu')
                   else (red, grn, blu)) .
        toUnboxed $
        arr
   in ImageRepa 8 . fromUnboxed (Z :. (3 :: Int) :. cols :. rows) . VU.concat $
      [r, g, b]

{-# INLINE normalizeComplexImage #-}
normalizeComplexImage ::
     ImageRepa (Complex Double) -> ImageRepa (Complex Double)
normalizeComplexImage (ImageRepa depth img) =
  let maxV = VU.maximum magVec
      minV = VU.minimum magVec
      maxVal = 2 ^ depth
      minVal = 0
      magVec = VU.map magnitude . toUnboxed $ img
   in if maxV == minV
        then (ImageRepa depth img)
        else ImageRepa depth .
             fromUnboxed (extent img) .
             VU.zipWith (\a b -> mkPolar b (phase a)) (toUnboxed img) .
             VU.map
               (\x -> (x - minV) / (maxV - minV) * (maxVal - minVal) + minVal) $
             magVec

{-# INLINE plotImageRepaComplex #-}
plotImageRepaComplex ::
     (Unbox e, Floating e, RealFrac e)
  => FilePath
  -> ImageRepa (Complex e)
  -> IO ()
plotImageRepaComplex filePath img =
  plotImageRepa filePath . complexImageToColorImage -- . normalizeComplexImage
  $
  img

{-# INLINE plotImageRepaComplexGray #-}
plotImageRepaComplexGray :: FilePath -> ImageRepa (Complex Double) -> IO ()
plotImageRepaComplexGray filePath img =
  let (ImageRepa depth arr) = complexImageToColorImage img
      (Z :. chs :. cols :. rows) = extent arr
  in plotImageRepa filePath .
     ImageRepa depth .
     computeS .
     extend (Z :. (1 :: Int) :. All :. All) .
     R.sumS .
     R.backpermute
       (Z :. cols :. rows :. chs)
       (\(Z :. c :. r :. k) -> (Z :. k :. c :. r)) $
     arr


text:
module Data.TestRational where

import           Prelude       hiding (Rational, gcd, negate, (*), (+), (-), (/), (^))

import           Domains.Field
import           Domains.Ring

import           Data.Complex
import           Data.Rational

import           TestUtil

fr :: Rational Int
fr = Rational 1 2

run :: IO ()
run = do
  putStrLn ("(1 / 2) = " ++ show fr)

ext: hs

max_stars_repo_name: pmilne/algebra

max_issues_repo_name: pmilne/algebra

max_forks_repo_name: pmilne/algebra


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Particle
-- Copyright   :  (c) 2016 FP Complete Corporation
-- License     :  MIT (see LICENSE)
-- Maintainer  :  dominic@steinitz.org
--
-- =The Theory
--
-- The particle filter, `runPF`, in this library is applicable to the
-- state space model given by
--
-- \[
-- \begin{aligned}
-- \boldsymbol{x}_i &= \boldsymbol{a}_i(\boldsymbol{x}_{i-1}) + \boldsymbol{\psi}_{i-1} \\
-- \boldsymbol{y}_i &= \boldsymbol{h}_i(\boldsymbol{x}_i) + \boldsymbol{\upsilon}_i
-- \end{aligned}
-- \]
--
-- where
--
-- * \(\boldsymbol{a_i}\)\ is some non-linear vector-valued possibly
-- time-varying state update function.
--
-- * \(\boldsymbol{\psi}_{i}\) are independent normally
-- distributed random variables with mean 0 representing the fact that
-- the state update is noisy: \(\boldsymbol{\psi}_{i} \sim {\cal{N}}(0,Q_i)\).
--
-- * \(\boldsymbol{h}_i\)\ is some non-linear vector-valued possibly time-varying
-- function describing how we observe the hidden state process.
--
-- * \(\boldsymbol{\upsilon}_i\) are independent normally
-- distributed random variables with mean 0 represent the fact that
-- the observations are noisy: \(\boldsymbol{\upsilon}_{i} \sim {\cal{N}}(0,R_i)\).
--
-- Clearly this could be generalised further; anyone wishing for such
-- a generalisation is encouraged to contribute to the library.
--
-- The smoother, `oneSmoothingPath` implements [Forward filtering /
-- backward
-- smoothing](https://en.wikipedia.org/wiki/Particle_filter#Backward_particle_smoothers)
-- and returns just one path from the particle filter; in most cases,
-- this will need to be run many times to provide good estimates of
-- the past. Note that `oneSmoothingPath` uses /all/ observation up
-- until the current time. This could be generalised to select only a
-- window of observations up to the current time. Again, contributions
-- to implement this generalisation are welcomed.
--
-- = An Extended Example
--
-- The equation of motion for a pendulum of unit length subject to
-- [Gaussian white
-- noise](https://en.wikipedia.org/wiki/White_noise#Mathematical_definitions)
-- is
--
-- \[
-- \frac{\mathrm{d}^2\alpha}{\mathrm{d}t^2} = -g\sin\alpha + w(t)
-- \]
--
-- We can discretize this via the usual [Euler method](https://en.wikipedia.org/wiki/Euler_method)
--
-- \[
-- \begin{bmatrix}
-- x_{1,i} \\
-- x_{2,i}
-- \end{bmatrix}
-- =
-- \begin{bmatrix}
-- x_{1,i-1} + x_{2,i-1}\Delta t \\
-- x_{2,i-1} - g\sin x_{1,i-1}\Delta t
-- \end{bmatrix}
-- +
-- \mathbf{q}_{i-1}
-- \]
--
-- where \(q_i \sim {\mathcal{N}}(0,Q)\) and
--
-- \[
-- Q
-- =
-- \begin{bmatrix}
-- \frac{q^c \Delta t^3}{3} & \frac{q^c \Delta t^2}{2} \\
-- \frac{q^c \Delta t^2}{2} & {q^c \Delta t}
-- \end{bmatrix}
-- \]
--
-- Assume that we can only measure the horizontal position of the
-- pendulum and further that this measurement is subject to error so that
--
-- \[
-- y_i = \sin x_i + r_i
-- \]
--
-- where \(r_i \sim {\mathcal{N}}(0,R)\).
--
--
-- First let's set the time step and the acceleration caused by
-- earth's gravity and the covariance matrices for the state and
-- observation.
--
--
-- > {-# LANGUAGE DataKinds #-}
-- >
-- > import Numeric.Particle
-- >
-- > deltaT, g :: Double
-- > deltaT = 0.01
-- > g      = 9.81
-- >
-- > bigQ :: Sym 2
-- > bigQ = sym $ matrix bigQl
-- >
-- > qc :: Double
-- > qc = 0.01
-- >
-- > bigQl :: [Double]
-- > bigQl = [ qc * deltaT^3 / 3, qc * deltaT^2 / 2,
-- >           qc * deltaT^2 / 2,       qc * deltaT
-- >         ]
-- >
-- > bigR :: Sym 1
-- > bigR  = sym $ matrix [0.2]
-- >
-- > data SystemState a = SystemState { x1  :: a, x2  :: a }
-- >   deriving Show
-- >
-- > newtype SystemObs a = SystemObs { y1  :: a }
-- >   deriving Show
-- >
--
-- We use the data generated using code made available with Simo
-- Särkkä's
-- [book](http://www.cambridge.org/gb/academic/subjects/statistics-probability/applied-probability-and-stochastic-networks/bayesian-filtering-and-smoothing)
-- which starts with an angle of 1.5 radians and an angular velocity
-- of 0.0 radians per second. We set our prior to have a mean of 1.6, not too far from the actual value.
--
-- > {-# LANGUAGE DataKinds #-}
-- >
-- > m0 :: PendulumState
-- > m0 = vector [1.6, 0]
-- >
-- > bigP :: Sym 2
-- > bigP = sym $ diag 0.1
-- >
-- > initParticles :: R.MonadRandom m =>
-- >                  m (Particles (SystemState Double))
-- > initParticles = V.replicateM nParticles $ do
-- >   r <- R.sample $ R.rvar (Normal m0 bigP)
-- >   let x1 = fst $ headTail r
-- >       x2 = fst $ headTail $ snd $ headTail r
-- >   return $ SystemState { x1 = x1, x2 = x2}
-- >
-- > nObs :: Int
-- > nObs = 200
-- >
-- > nParticles :: Int
-- > nParticles = 20
-- >
-- > (.+), (.*), (.-) :: (Num a) => V.Vector a -> V.Vector a -> V.Vector a
-- > (.+) = V.zipWith (+)
-- > (.*) = V.zipWith (*)
-- > (.-) = V.zipWith (-)
-- >
-- > stateUpdateP :: Particles (SystemState Double) ->
-- >                 Particles (SystemState Double)
-- > stateUpdateP xPrevs = V.zipWith SystemState x1s x2s
-- >   where
-- >     ix = V.length xPrevs
-- >
-- >     x1Prevs = V.map x1 xPrevs
-- >     x2Prevs = V.map x2 xPrevs
-- >
-- >     deltaTs = V.replicate ix deltaT
-- >     gs = V.replicate ix g
-- >     x1s = x1Prevs .+ (x2Prevs .* deltaTs)
-- >     x2s = x2Prevs .- (gs .* (V.map sin x1Prevs) .* deltaTs)
--
--
-- <<diagrams/src_Numeric_Particle_diagV.svg#diagram=diagV&height=600&width=500>>
--
-- === Code for Plotting
--
-- The full code for plotting the results:
--
-- > {-# LANGUAGE ExplicitForAll #-}
-- > {-# LANGUAGE TypeOperators  #-}
-- > {-# LANGUAGE DataKinds      #-}
-- >
-- > import qualified Graphics.Rendering.Chart as C
-- > import Graphics.Rendering.Chart.Backend.Diagrams
-- > import Data.Colour
-- > import Data.Colour.Names
-- > import Data.Default.Class
-- > import Control.Lens
-- >
-- > import Data.Csv
-- > import System.IO hiding ( hGetContents )
-- > import Data.ByteString.Lazy ( hGetContents )
-- > import qualified Data.Vector as V
-- >
-- > import           Data.Random.Distribution.Static.MultivariateNormal ( Normal(..) )
-- > import qualified Data.Random as R
-- > import           Data.Random.Source.PureMT ( pureMT )
-- > import           Control.Monad.State ( evalState, replicateM )
-- >
-- > import Data.List ( transpose )
-- >
-- > import GHC.TypeLits ( KnownNat )
-- >
-- > import Numeric.LinearAlgebra.Static
-- >
-- > chartEstimated :: String ->
-- >               [(Double, Double)] ->
-- >               [(Double, Double)] ->
-- >               [(Double, Double)] ->
-- >               C.Renderable ()
-- > chartEstimated title acts obs ests = C.toRenderable layout
-- >   where
-- >
-- >     actuals = C.plot_lines_values .~ [acts]
-- >             $ C.plot_lines_style  . C.line_color .~ opaque red
-- >             $ C.plot_lines_title .~ "Actual Trajectory"
-- >             $ C.plot_lines_style  . C.line_width .~ 1.0
-- >             $ def
-- >
-- >     measurements = C.plot_points_values .~ obs
-- >                  $ C.plot_points_style  . C.point_color .~ opaque blue
-- >                  $ C.plot_points_title .~ "Measurements"
-- >                  $ def
-- >
-- >     estimas = C.plot_lines_values .~ [ests]
-- >             $ C.plot_lines_style  . C.line_color .~ opaque black
-- >             $ C.plot_lines_title .~ "Inferred Trajectory"
-- >             $ C.plot_lines_style  . C.line_width .~ 1.0
-- >             $ def
-- >
-- >     layout = C.layout_title .~ title
-- >            $ C.layout_plots .~ [C.toPlot actuals, C.toPlot measurements, C.toPlot estimas]
-- >            $ C.layout_y_axis . C.laxis_title .~ "Angle / Horizontal Displacement"
-- >            $ C.layout_y_axis . C.laxis_override .~ C.axisGridHide
-- >            $ C.layout_x_axis . C.laxis_title .~ "Time"
-- >            $ C.layout_x_axis . C.laxis_override .~ C.axisGridHide
-- >            $ def
-- >
-- > stateUpdateNoisy :: R.MonadRandom m =>
-- >                     Sym 2 ->
-- >                     Particles (SystemState Double) ->
-- >                     m (Particles (SystemState Double))
-- > stateUpdateNoisy bigQ xPrevs = do
-- >   let xs = stateUpdateP xPrevs
-- >
-- >       x1s = V.map x1 xs
-- >       x2s = V.map x2 xs
-- >
-- >   let ix = V.length xPrevs
-- >   etas <- replicateM ix $ R.sample $ R.rvar (Normal 0.0 bigQ)
-- >
-- >   let eta1s, eta2s :: V.Vector Double
-- >       eta1s = V.fromList $ map (fst . headTail) etas
-- >       eta2s = V.fromList $ map (fst . headTail . snd . headTail) etas
-- >
-- >   return (V.zipWith SystemState (x1s .+ eta1s) (x2s .+ eta2s))
-- >
-- > obsUpdate :: Particles (SystemState Double) ->
-- >              Particles (SystemObs Double)
-- > obsUpdate xs = V.map (SystemObs . sin . x1) xs
-- >
-- > weight :: forall a n . KnownNat n =>
-- >           (a -> R n) ->
-- >           Sym n ->
-- >           a -> a -> Double
-- > weight f bigR obs obsNew = R.pdf (Normal (f obsNew) bigR) (f obs)
-- >
-- > runFilter :: Particles (SystemObs Double) -> V.Vector (Particles (SystemState Double))
-- > runFilter pendulumSamples = evalState action (pureMT 19)
-- >   where
-- >     action = do
-- >       xs <- initParticles
-- >       scanMapM
-- >         (runPF (stateUpdateNoisy bigQ) obsUpdate (weight f bigR))
-- >         return
-- >         xs
-- >         pendulumSamples
-- >
-- > testSmoothing :: Particles (SystemObs Double) -> Int -> [Double]
-- > testSmoothing ss n = V.toList $ evalState action (pureMT 23)
-- >   where
-- >     action = do
-- >       xss <- V.replicateM n $ oneSmoothingPath (stateUpdateNoisy bigQ) (weight h bigQ) nParticles (runFilter ss)
-- >       let yss = V.fromList $ map V.fromList $
-- >                 transpose $
-- >                 V.toList $ V.map (V.toList) $
-- >                 xss
-- >       return $ V.map (/ (fromIntegral n)) $ V.map V.sum $ V.map (V.map x1) yss
-- >
-- > type PendulumState = R 2
-- >
-- > f :: SystemObs Double -> R 1
-- > f = vector . pure . y1
-- >
-- > h :: SystemState Double -> R 2
-- > h u = vector [x1 u , x2 u]
-- >
-- > diagV = do
-- >   h <- openFile "data/matlabRNGs.csv" ReadMode
-- >   cs <- hGetContents h
-- >   let df = (decode NoHeader cs) :: Either String (V.Vector (Double, Double))
-- >   case df of
-- >     Left _ -> error "Whatever"
-- >     Right generatedSamples -> do
-- >       let preObs = V.take nObs $ V.map fst generatedSamples
-- >       let obs = V.toList preObs
-- >       let acts = V.toList $ V.take nObs $ V.map snd generatedSamples
-- >       let nus = take nObs (testSmoothing (V.map SystemObs preObs) 50)
-- >       denv <- defaultEnv C.vectorAlignmentFns 600 500
-- >       let charte = chartEstimated "Particle Smoother"
-- >                                   (zip [0,1..] acts)
-- >                                   (zip [0,1..] obs)
-- >                                   (zip [0,1..] nus)
-- >       return $ fst $ runBackend denv (C.render charte (600, 500))
--
-----------------------------------------------------------------------------

{-# OPTIONS_GHC -Wall                     #-}

{-# LANGUAGE BangPatterns                 #-}

module Numeric.Particle (
    runPF
  , oneSmoothingPath
  , Particles
  , Path
  ) where

import           Data.Random hiding ( StdNormal, Normal )
import           Control.Monad
import qualified Data.Vector as V
import           Data.Vector ( Vector )
import           Data.Bits ( shiftR )

type Particles a = Vector a -- ^ As an aid for the reader
type Path a      = Vector a -- ^ As an aid for the reader

runPF
  :: MonadRandom m
  => (Particles a -> m (Particles a)) -- ^ System evolution at a point
                                      -- \(\boldsymbol{a}_i(\boldsymbol{x})\)
  -> (Particles a -> Particles b)     -- ^ Measurement operator at a point \(\boldsymbol{h}_i\)
  -> (b -> b -> Double)               -- ^ Observation probability density function \(p(\boldsymbol{y}_k|\boldsymbol{x}^{(i)}_k)\)
  -> Particles a                      -- ^ Current estimate
                                      -- \(\big\{\big(w^{(i)}_k, \boldsymbol{x}^{(i)}_k\big) : i \in \{1,\ldots,N\}\big\}\)
                                      -- where \(N\) is the number of particles
  -> b                                -- ^ New measurement \(\boldsymbol{y}_i\)
  -> m (Particles a)                  -- ^ New estimate
                                      -- \(\big\{\big(w^{(i)}_{k+1}, \boldsymbol{x}^{(i)}_{k+1}\big) : i \in \{1,\ldots,N\}\big\}\)
                                      -- where \(N\) is the number of particles

runPF stateUpdate obsUpdate weight statePrevs obs = do
  stateNews <- stateUpdate statePrevs
  let obsNews = obsUpdate stateNews
  let weights       = V.map (weight obs) obsNews
      cumSumWeights = V.tail $ V.scanl (+) 0 weights
      totWeight     = V.last cumSumWeights
      nParticles    = V.length statePrevs
  vs <- V.replicateM nParticles (sample $ uniform 0.0 totWeight)
  let js = indices cumSumWeights vs
      stateTildes = V.map (stateNews V.!) js
  return stateTildes

oneSmoothingPath
  :: MonadRandom m
  => (Particles a -> m (Particles a)) -- ^ System evolution at a point
                                      -- \(\boldsymbol{a}_i(\boldsymbol{x})\)
  -> (a -> a -> Double)               -- ^ State probability density function \(p(\boldsymbol{x}^{(i)}_k|\boldsymbol{x}^{(i)}_{k-1})\)
  -> Int                              -- ^ Number of particles \(N\)
  -> (Vector (Particles a))           -- ^ Filtering estimates
                                      -- \(\bigg\{\big\{\big(w^{(i)}_{k+1}, \boldsymbol{x}^{(i)}_{k+1}\big) : i \in \{1,\ldots,N\}\big\} : k \in \{1,\ldots,T\}\bigg\}\)
                                      -- where \(T\) is the number of timesteps
  -> m (Path a)
                                      -- ^ A path for the smoothed particle \(\{\boldsymbol{x}_k : k \in T\}\) where \(T\) is the number of timesteps
oneSmoothingPath stateUpdate weight nParticles filterEstss = do
  let ys = filterEstss
  ix <- sample $ uniform 0 (nParticles - 1)
  let xn = (V.head ys) V.! ix
  scanMapM (oneSmoothingStep stateUpdate weight) return xn (V.tail ys)

oneSmoothingStep :: MonadRandom m =>
         (Particles a -> m (Particles a)) ->
         (a -> a -> Double) ->
         a ->
         Particles a ->
         m a
oneSmoothingStep stateUpdate
                 stateDensity
                 smoothingSampleAtiPlus1
                 filterSamplesAti = do it
  where
    it = do
      mus <- stateUpdate filterSamplesAti
      let weights =  V.map (stateDensity smoothingSampleAtiPlus1) mus
          cumSumWeights = V.tail $ V.scanl (+) 0 weights
          totWeight     = V.last cumSumWeights
      v <- sample $ uniform 0.0 totWeight
      let ix = binarySearch cumSumWeights v
          xnNew = filterSamplesAti V.! ix
      return $ xnNew

indices :: Vector Double -> Vector Double -> Vector Int
indices bs xs = V.map (binarySearch bs) xs

binarySearch :: (Ord a) =>
                Vector a -> a -> Int
binarySearch vec x = loop 0 (V.length vec - 1)
  where
    loop !l !u
      | u <= l    = l
      | otherwise = let e = vec V.! k in if x <= e then loop l k else loop (k+1) u
      where k = l + (u - l) `shiftR` 1

scanMapM :: Monad m => (s -> a -> m s) -> (s -> m b) -> s -> Vector a -> m (Vector b)
scanMapM f g !s0 !xs
  | V.null xs = do
    r <- g s0
    return $ V.singleton r
  | otherwise = do
    s <- f s0 (V.head xs)
    r <- g s0
    liftM (r `V.cons`) (scanMapM f g s (V.tail xs))

max_stars_repo_name: idontgetoutmuch/Kalman

max_issues_repo_name: f-o-a-m/Kalman

max_forks_repo_name: f-o-a-m/Kalman


text:
-- Exercisee/Ch2/Pt1/Ex7.hs
module Ex7 where

-- Exercise 7
-- Add data types and parsers to support the full numeric tower of Scheme
-- numeric types. Haskell has built-in types to represent many of these; check
-- the Prelude. For others, you can define compound types that represent eg. a
-- `Rational` as a numerator and denominator, or a `Complex` as a real and an
-- imaginary part (each itself a `Real`).

import Data.Char                            (digitToInt, isSpace)
import Data.Complex
import Data.Ratio
import Control.Monad
import Numeric
import System.Environment

import Text.ParserCombinators.Parsec hiding (spaces)

data LispVal = Atom String
             | List [LispVal]
             | DottedList [LispVal] LispVal -- (a b c ... . z)
             | Number LispNum
             | Character String
             | String String
             | Bool Bool
             deriving Show -- debug purposes

data LispNum = Complex (Complex Float)
             | Real Float
             | Rational (Ratio Integer)
             | Integer Integer
             deriving (Eq, Show)

parseCharacter :: Parser LispVal
parseCharacter = do
  string "#\\"
  character <- many $ satisfy $ not . isSpace
  return $ Character character

parseString :: Parser LispVal
parseString = do
  char '"'
  -- order is important here; need to try parsing an escape sequence first because
  -- otherwise we fail when we reach the second character of the escape
  -- sequence
  x <- many $ choice $ escChars ++ [nonQuote]
  char '"'
  return $ String x
  where
    nonQuote = noneOf "\""
    -- taken from here:
    -- https://en.wikipedia.org/wiki/Escape_sequences_in_C#Table_of_escape_sequences
    -- (excluded characters by oct/hex codes)
    escChars = [ char '\\' >> char x | x <- "abfnrtv\\'\"?" ]

parseBool :: Parser LispVal
parseBool = do
  try (char '#')
  v <- oneOf "tf"
  return $ case v of
    't' -> Bool True
    'f' -> Bool False

parseAtom :: Parser LispVal
parseAtom = do
              first <- letter <|> symbol
              rest <- many (letter <|> digit <|> symbol)
              return $ Atom (first:rest)

parseNumber :: Parser LispVal
parseNumber = do
  n <- (try parseBin) <|> (try parseOct) <|> (try parseHex) <|> (try parseFlt)
    <|> (try parseRat) <|> (try parseCmp) <|> parseDec
  return $ Number n
  where
    parseBin = do
      string "#b"
      binStr <- many1 $ oneOf "01"
      -- cribbed from http://stackoverflow.com/a/26961027/1893155
      let binVal = foldl (\acc x -> acc * 2 + digitToInt x) 0 binStr
      return $ Integer (toInteger binVal)
    parseOct = do
      string "#o"
      octStr <- many1 octDigit
      let octVal = fst $ (readOct octStr) !! 0
      return $ Integer octVal
    parseDec = parseDecNoPre <|> parseDecPre
    parseDecNoPre = do
      decStr <- many1 digit
      return $ (Integer . read) decStr
    parseDecPre = do
      string "#d"
      decStr <- many1 digit
      return $ (Integer . read) decStr
    parseHex = do
      string "#x"
      hexStr <- many1 hexDigit
      let hexVal = fst $ (readHex hexStr) !! 0
      return $ Integer hexVal
    parseFlt = do
      whole <- many1 digit
      string "."
      mantissa <- many1 digit
      let fltVal = fst $ (readFloat (whole ++ "." ++ mantissa) !! 0)
      return $ Real $ fltVal
    parseRat = do
      m <- many1 digit
      many space
      string "/"
      many space
      n <- many1 digit
      let (m', n') = (fst $ (readDec m) !! 0, fst $ (readDec n) !! 0)
      return $ Rational $ m' % n'
    parseCmp = do
      (Real a) <- (try parseDec) <|> parseFlt
      many space
      string "+"
      many space
      (Real b) <- (try parseDec) <|> parseFlt
      char 'i'
      return $ Complex $ a :+ b

-- TODO : I don't like all these `try`s -- better way to do this?
parseExpr :: Parser LispVal
parseExpr = (try parseNumber)
         <|> (try parseString)
         <|> (try parseCharacter)
         <|> (try parseBool)
         <|> (try parseAtom)

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

spaces :: Parser ()
spaces = skipMany1 space

readExpr :: String -> String
readExpr input = case parse parseExpr "lisp" input of
  Left err -> "No match: " ++ show err
  Right val -> "Found value"

main :: IO ()
main = do
  (expr:_) <- getArgs
  putStrLn (readExpr expr)

max_stars_repo_name: EFulmer/haskell-scheme-wikibook

max_issues_repo_name: EFulmer/haskell-scheme-wikibook

max_forks_repo_name: EFulmer/haskell-scheme-wikibook


text:
import           Control.Applicative            ( liftA2 )

import           Math.Combinat.Partitions.Set
import           Numeric.LinearAlgebra
import           Numeric.LinearAlgebra.Data

(∨) = min
(f .: g) x y = f $ g x y

-- Partitions [l]→[n] factoring through a 2-fibration [2l]↠[l]
isTwoPartition = all even . map length . fromSetPartition

-- l≤n: n

twoPartitions n l =
  filter isTwoPartition . setPartitionsWithKParts (n ∨ (l - n))

-- Constraint: lists must be finite (forcing `length` on both)
relationToMatrix :: (r -> c -> Bool) -> [r] -> [c] -> Matrix R
relationToMatrix p rs cs = (length rs >< length cs) $ liftA2 ((↯) .: p) rs cs
  where (↯) b = if b then 1 else 0

max_stars_repo_name: hseg/experiments

max_issues_repo_name: hseg/experiments

max_forks_repo_name: hseg/experiments


text:
{-# LANGUAGE GeneralizedNewtypeDeriving, BangPatterns, NoMonomorphismRestriction, ScopedTypeVariables, CPP #-}
-- vim: sts=2:sw=2:ai:et
--module Rolling where

import Data.Word
import Data.Int
import Data.Bits
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as UM
import Control.Monad.ST
import qualified Data.ByteString as B
import qualified Data.ByteString.Unsafe as B
import qualified Data.ByteString.Lazy as L
import Pipes
import Pipes.Lift
import qualified Pipes.Prelude as P
import Control.Monad.State.Strict
import Control.Exception (assert)
import qualified Test.QuickCheck as QC
import Data.List
import qualified Control.Lens as L
import Control.Lens.Operators
import Control.DeepSeq
import Criterion.Main
import Statistics.Sample.Histogram
-- import Debug.Trace

window :: Int
mask :: Word64

#if 0
window = 16
mask = 0xf
#else
window = 256
mask = 0x3fff
#endif

testMask :: Word64 -> Bool
testMask x = x .&. mask == mask
{-# INLINE testMask #-}

hash :: Word8 -> Word64
{-
hash x = v8
  where v = fromIntegral x
        vv = v `shiftL` 8 .|. v
        vvvv = vv `shiftL` 16 .|. vv
        v8 = vvvv `shiftL` 32 .|. vvvv
        -}
-- hash x = 31 * fromIntegral x
hash x = lut `U.unsafeIndex` fromIntegral x
{-# INLINE hash #-}

rhash :: Word8 -> Word8 -> Word64 -> Word64
rhash old new h = h `rotateL` 1 `xor` (hash old `rotateL` window) `xor` hash new

(+>) :: Word64 -> Word64 -> Word64
(+>) !o !n = (o `rotateL` window) `xor` n
{-# INLINE (+>) #-}

hashCombine :: Word64 -> Word64 -> Word64
hashCombine !x !y = x `rotateL` 1 `xor` y
{-# INLINE hashCombine #-}

type Data = B.ByteString -- S.Vector Word8

roll :: [Word64] -> [Word64]
roll hashed = tail $ scanl hashCombine h (zipWith (+>) hashed (drop window hashed))
  where h = foldl' hashCombine 0 (take window hashed)

prop_rolls :: QC.Property
prop_rolls = QC.forAll inputList $ \a ->
             QC.forAll (QC.suchThat inputList (\v -> length v > window)) $ \b ->
               last (roll b) == last (roll (a ++ b))

contiguous :: Data -> [Data]
contiguous xs = zipWith (\a b -> B.take (b - a) (B.drop a xs)) (0:boundaries) boundaries
  where hashed = map hash (B.unpack xs)
        rolled = roll hashed
        markers = findIndices testMask rolled
        boundaries = (++[B.length xs]) $ map (+(window+1)) markers

cprop_allInputIsOutput :: QC.Property
cprop_allInputIsOutput = QC.forAll inputData $ \xs -> B.concat (contiguous xs) == xs

cprop_prefix :: QC.Property
cprop_prefix = QC.forAll inputData $ \xs -> QC.forAll inputData $ \ys ->
  let a = contiguous (xs `B.append` ys); b = contiguous xs in init' b `isPrefixOf` a

cprop_suffix :: QC.Property
cprop_suffix = QC.forAll inputData $ \xs -> QC.forAll inputData $ \ys ->
  let a = contiguous (xs `B.append` ys); c = contiguous ys in tail' c `isSuffixOf` a

cprop_valid :: QC.Property
cprop_valid = cprop_allInputIsOutput QC..&. cprop_prefix QC..&. cprop_suffix

data HashState
  = HashState {
      _lastHash :: {-# UNPACK #-} !Word64
    , _lastWindow :: !Data
    }

initialState :: HashState
initialState = HashState 0 B.empty

lastHash :: L.Lens' HashState Word64
lastHash f (HashState h w) = fmap (\h' -> HashState h' w) (f h)

lastWindow :: L.Lens' HashState Data
lastWindow f (HashState h w) = fmap (HashState h) (f w)

data Output = Partial { getOutput :: Data } | Complete { getOutput :: Data }
  deriving (Show)

isComplete :: Output -> Bool
isComplete (Partial _) = False
isComplete (Complete _) = True

rollingBoundaries :: Data -> Data -> Word64 -> (Word64, U.Vector Int)
rollingBoundaries !old !new !h0 = runST $ do mv <- UM.new (B.length new)
                                             (h', len) <- runner mv
                                             v <- U.unsafeFreeze (UM.take len mv)
                                             return (h', v)
  where
    runner mv = go h0 0 0
      where
        go !h !iOut !iIn | iIn < B.length new =
                           do
                              let hi = hash (old `B.unsafeIndex` iIn) +> hash (new `B.unsafeIndex` iIn)
                                  h' = hashCombine h hi
                              iOut' <- case testMask h' of
                                         True -> do UM.unsafeWrite mv iOut (iIn + 1)
                                                    return (iOut + 1)
                                         False -> return iOut
                              go h' iOut' (iIn + 1)
                         | otherwise = return (h, iOut)

rollsplitP :: forall m. Monad m => Pipe Data Output (StateT HashState m) ()
rollsplitP =
    do w <- L.use lastWindow
       if B.length w >= window
        then await >>= warmedUpPhase
        else await >>= initialPhase
  where
    initialPhase !x =
      do
        w <- L.use lastWindow
        assert (B.length w < window) (return ())

        let n = window - B.length w
            (xi, xn) = B.splitAt n x
            w' = w `B.append` xi
        lastWindow .= w'

        h <- L.use lastHash
        lastHash .= foldl' hashCombine h (map hash (B.unpack $ B.take n x))

        yield (Partial xi)

        if B.length w' < window
        then await >>= initialPhase
        else warmedUpPhase xn

    warmedUpPhase :: Data -> Pipe Data Output (StateT HashState m) ()
    warmedUpPhase !x =
      do
        w <- L.use lastWindow
        assert (B.length w == window) (return ())

        w' <- hoist (L.zoom lastHash) $
          do
            let len = B.length x
            chow w (B.take window x)
            when (len > window) $ chow x (B.drop window x)
            return (B.drop len w `B.append` B.drop (len - window) x)
        lastWindow .= w'
        await >>= warmedUpPhase

    chow !old !new = assert (B.length old >= B.length new) $
      do
        h <- get
        let (newH, boundaries) = rollingBoundaries old new h
        let sliceAction a b = do
              yield (Complete (B.unsafeTake (b - a) (B.unsafeDrop a new)))
              return b
        n <- U.foldM sliceAction 0 boundaries
        when (n < B.length new) $ yield (Partial (B.drop n new))
        put newH
    {-# INLINE chow #-}

recombine :: Monad m => Int64 -> Int64 -> Producer Output m a -> Producer L.ByteString m a
recombine nmin nmax = loop L.empty
  where
    loop d p =
      case L.length d of
        n | n < nmax ->
            do e <- lift (next p)
               case e of
                 Left v -> yield d >> return v
                 Right (output, p') -> if isComplete output && L.length d' >= nmin
                                       then yield d' >> loop L.empty p'
                                       else loop d' p'
                   where d' = d `L.append` L.fromStrict (getOutput output)
          | otherwise ->
            do yield (L.take nmax d)
               loop (L.drop nmax d) p

rollsplit :: Monad m => Int64 -> Int64 -> Producer Data m () -> Producer L.ByteString m ()
rollsplit nmin nmax p = recombine nmin nmax $ evalStateP initialState $ hoist lift p >-> rollsplitP

rollsplitL :: [Data] -> [L.ByteString]
rollsplitL = filter (not . L.null) . rollsplitL'

rollsplitL' :: [Data] -> [L.ByteString]
rollsplitL' xs = P.toList $ rollsplit 0 (maxBound :: Int64) (each xs)

inputList :: QC.Arbitrary a => QC.Gen [a]
inputList = QC.sized $ \n -> do
  k <- QC.choose (0,n)
  a <- QC.choose (0,window)
  let len = max 0 $ (k-1)*window+a
  QC.vector len

inputData :: QC.Gen Data
inputData = fmap B.pack inputList

init' :: [a] -> [a]
init' xs = take (length xs - 1) xs

tail' :: [a] -> [a]
tail' xs = drop 1 xs

prop_allInputIsOutput :: QC.Property
prop_allInputIsOutput = QC.forAll (QC.listOf inputData) $ \xs -> L.concat (rollsplitL xs) == L.fromChunks xs

prop_inputSplit :: QC.Property
prop_inputSplit = QC.forAll (QC.listOf inputData) $ \xs -> rollsplitL xs == rollsplitL [B.concat xs]

prop_prefix :: QC.Property
prop_prefix = QC.forAll inputData $ \xs -> QC.forAll inputData $ \ys ->
  let a = rollsplitL [xs, ys]; b = rollsplitL [xs] in init' b `isPrefixOf` a

prop_suffix :: QC.Property
prop_suffix = QC.forAll inputData $ \xs -> QC.forAll inputData $ \ys ->
  let a = rollsplitL [xs, ys]; c = rollsplitL [ys] in tail' c `isSuffixOf` a

prop_concat :: QC.Property
prop_concat = prop_prefix QC..&. prop_suffix

prop_eq :: QC.Property
prop_eq = QC.forAll inputData $ \xs -> rollsplitL' [xs] == map L.fromStrict (contiguous xs)

prop_valid :: QC.Property
prop_valid = prop_allInputIsOutput QC..&. prop_inputSplit QC..&. prop_concat QC..&. prop_eq

qc :: IO ()
qc = QC.quickCheckWith (QC.stdArgs { QC.maxSuccess = 700 }) (cprop_valid QC..&. prop_valid)

test :: Int64 -> Int64 -> [Data] -> IO ()
test nmin nmax xs = runEffect $ for (rollsplit nmin nmax (each xs)) (lift . print)

test2 :: [Data] -> IO ()
test2 xs = runEffect $ for (evalStateP initialState $ each xs >-> rollsplitP) (lift . print)

lut :: U.Vector Word64
lut = U.fromList [
    0xfead5b707dc7705c, 0x377c1e06dc1e45cf, 0x0184179586d5ae76, 0xd23aa044f8193aa6,
    0xbd8ef5fcde7bd95e, 0x29a822b00a75ea90, 0x5ba03c1b2fdc2f86, 0x4f67d80bad410270,
    0xfcc4b6b0cb67bb75, 0x8f4359ea8777f5d0, 0x5a110ec6371430f5, 0xe15dae4e9709aa66,
    0xf8008efefbc29115, 0xd667f1ec0a5bcccf, 0xf56a2815cb29dbf8, 0xfc684381b1846cad,
    0x15be2b7ed5231683, 0x0c226b1f911732c5, 0xe449b5885fd1d930, 0xf6d86c27a242cf1d,
    0xc45cd857fdd32ca4, 0x36de9a99be8c3419, 0x8c34f2ee7050433f, 0x9fca9c3c61c05ce8,
    0xf3727fe070e026f1, 0x019172d2948847be, 0x547bcbfaaef62543, 0x97df2b9f6c70001b,
    0x2833bcc8513a3ae8, 0x63fb313206823b31, 0xe97ba036be16b286, 0xdde1143e9e8f7684,
    0xc4e0d8e6e491054c, 0x7c488b71d4cab5fd, 0x2e0986d1f8e2d362, 0x97978ef2435b5cb1,
    0x5506c42930dfc2a4, 0x3205d4c5886673cf, 0x22d5bfe8ba771cc5, 0x39503d4cb89afd6b,
    0x9fc1cad11e5b5755, 0xce967aae1fccb3c3, 0x0f164d18916d4043, 0x113893188dd7e602,
    0xedb9343bcf5f16ef, 0x0ce3b155cfc00cbf, 0x144c6c104492dff7, 0xe1ab67476c3c787f,
    0x06357095472efed7, 0x29573d1071bb1237, 0xe27e7ce26dcd54c6, 0x03a8623e84a0e89c,
    0x381194ea8cb4f8fe, 0x53cbe859f67ae71c, 0x85e08cd4432b31f5, 0x21f9758ed837239d,
    0x835aad64a22264b4, 0xf83ec94f1a0159f3, 0x016a63af376515a4, 0x8af0e82577e30b6e,
    0x4e3717bd197f8c8d, 0x20193e9699604cea, 0xfa8eb53310f65183, 0xb4980361f9a7e3b3,
    0xccb9771f04663244, 0xe975bed27cfd7a05, 0x46e9eb7dd3ffd802, 0xbe4b5af269bf4b64,
    0xf3ff155176f34e7e, 0xfc8e0e5c20814c3f, 0x3910b75cbe6677f5, 0x86fb5f8ef54ad45a,
    0xad430e08aa6aea72, 0x6a82d96b64a7f48c, 0x8aa8af38d787327e, 0x3a10990b93334a3d,
    0x010c38275c1fc3ac, 0x03eac81f78d29425, 0x356d5c2052ac066f, 0x9a14a98e51d4eafc,
    0x9c07059dac0831b5, 0x51c3aa247abb4a6c, 0x7f9d96a58c371a56, 0xe098fe430cedd615,
    0x960e3a1dbe524565, 0xf0401c427323f9e8, 0x2c22d55a6135331c, 0x8684f0837cb96bc9,
    0xb7ac495220cc7e9b, 0x4083bd84b86ad15e, 0xd7434c89c5464bfd, 0x6af7d49f8b658a9d,
    0x1b9878f861dd7d00, 0x6f5b7c3c3baaca36, 0xab6e6179f8b69d37, 0xda4c68dc98bcecd5,
    0x42cb0a636ba3fa38, 0x04f2814021d2f99a, 0x53a069c086a24b1d, 0xce39ccfebf0fd3cf,
    0x6da4576c5b938392, 0x73961a564d3684af, 0xe62079e996d934ab, 0x0fe43a5aded65708,
    0x8ea00b7c1ba4a6b4, 0xb3f325b0d1acc61a, 0xb7d68c724ee32277, 0x8838a5a56df55c2e,
    0x639e1a16eebc31dc, 0xc824476623ee8713, 0x6cfd35dd41d4a35f, 0xc0cf8dfb2eb67ca4,
    0x9b2db80c20bf8281, 0x8f43090206cbc3a5, 0xf018850ed3cce401, 0xdd9c873dbc24dacf,
    0x8c310f7792646b4a, 0xd64514da16af625c, 0x7de1ce3ec600b27d, 0x97742bfa39500785,
    0x37e57a0934a9e321, 0x0872c00e96058f5f, 0x9a9291f9452ff389, 0x1240421b117dc351,
    0x179b86fb16a829ed, 0xf773ae083f3cd908, 0xc32304664f8ce10f, 0x3bf48227de1ad114,
    0x1e051eb805a697ad, 0x70ef7848e5ff8346, 0xe7dc2e245cb0d7f3, 0xf9ccdaee9f48ca62,
    0x2ef082a2b4419ad3, 0x8a9287ed4ecd5750, 0x9f7b959d7c2ff1ba, 0x1b86fb1047eee0a9,
    0x342ac90ae166bfb4, 0x9c28c8b1a92ab370, 0x7248bbdc9faaf6c1, 0x1a3b7cdd85269cfa,
    0xcf3af94626a8c04f, 0x7a4ba5bef4f2856c, 0xa80dd551adc99bb8, 0xd7cef340d2e894ba,
    0x134c05ffd4928eba, 0x78b5ec19a589ff89, 0xd8d23f88cbd19231, 0xa95a65c91e98558c,
    0x0a996675f053aa80, 0xf35d31427e2dce98, 0x9985000bf5aefe7b, 0x73ac3f7cc6a7d6e0,
    0xc48e7406cc446de3, 0xf919aa74e167dbbe, 0x8b09e12fc48dbbcd, 0x33b5f1fbafc21ead,
    0xfbbf051f129c275d, 0x3239b2a8ee124a62, 0x8984d8d6b4a2ad6a, 0xe2293308df3acd2b,
    0x2529297459fa7fac, 0x36c72d90491e670f, 0xe7ca9d442ddae3b3, 0x8fa38fa2a012ce0d,
    0x5caac60bfad2c355, 0xa92158609a2e1285, 0x8a097df8fdb4c6cf, 0x1f1b6f8230e90ff0,
    0x8369b558bec7e569, 0x2e9d4f7aabcb210d, 0xafb7d767ac296c18, 0x2f5cd5c2515d81c1,
    0x45934465add6c676, 0xa5a45e774fd55fc3, 0xb4cb7e6e96d2a59d, 0xfdc6635b7ec562b0,
    0x1bcedbd2ea8b4bf9, 0x0ab5b24bb7ccf99d, 0x31aa59853645b444, 0x9149c77d16594ff5,
    0x77d9fee30ed21a67, 0x12982d192e88015b, 0xd45d1f788a6a1ef3, 0x332d24e0edf8ca0a,
    0x88c7695975b9d19b, 0x065f29d44b168060, 0xd2a2a94e28c99de3, 0x8235595da56d3cd7,
    0x4d58d451daa20b65, 0x36a1edff55e32372, 0xdf5c5c1e6b75c40b, 0x85c8a6bf77ecf2d3,
    0x7d72ca70ec2882f2, 0xfe286fdc7db952ca, 0x13e308d730603af4, 0x7803625ac832d270,
    0x920af03f31044bc6, 0x1fe4f2d207fecc03, 0x38537ef4067d59e9, 0x0eca7ef2100d20fc,
    0x0a6fe0e637366afc, 0x9d52574f1a85c9e5, 0xd0a0fdbd96ed75d2, 0xb6785ac2254c0d92,
    0x708122c30c97ff1f, 0x3e265597b73bc11e, 0xb3a870c781058225, 0xd93d905b304c4ec6,
    0x52ff676e213f6ccb, 0xde7503df2ee1aed6, 0xcc95a88b5c5fc05a, 0x351d157d93a36b21,
    0xf33702abc575fb08, 0xedcec8d4e602f76d, 0x210d869a4002a35a, 0xed6debb5a777ff11,
    0x4aa4ef159b20a5da, 0x1248e5993f479e53, 0x6afdd7979ed93e76, 0x0048f115bbf8fcc9,
    0x3c43dcc810f39b4c, 0xd596e217217ba38d, 0xf026411acc2ffc3f, 0xa9bfe1a48496f61c,
    0x3c5d0feeb403ffe6, 0xf3cc4252ff0a0e4b, 0x9287b23486699a23, 0xfcb9175b7c6ad7f1,
    0xd28a3ca6572757d0, 0x9679f32c0d2f01f8, 0x9aec761f1611e561, 0x00d429298b4544fc,
    0x08aa8292dcaff98f, 0x36c88c6cf0188da2, 0x674a34c8653e66b8, 0x0e6316f1731e02b8,
    0xd049743dfb30601e, 0xe99b74d4efa92d88, 0xa41f9c33729e01ad, 0xe0a2c2031c01b13d,
    0x04a396e2dfc4f61f, 0xdb2556b846e68964, 0xa5de683723f23704, 0x456f750ec820ea4f,
    0x47b81bb8dda336aa, 0x7c0b4cab5dffb32e, 0x020a4567a37ada75, 0x66325b8ce0bc3890,
    0x566e9b6d6a194d5e, 0x53efbb71787c4049, 0xee775a2e794219de, 0xa65c279ca51a6a46,
    0xaac56091fcec9c38, 0x0b34953b386ed0c4, 0xde46839785d1b945, 0x623a65d725974df2
  ]

stats :: [Int64] -> IO ()
stats xs = print (minimum xs) >> print (maximum xs) >> print (fromIntegral (sum xs) / fromIntegral (length xs) :: Double) >> print hist
  where
    hist = histogram_ 12 0 32 (U.map (logBase 2 . fromIntegral) $ U.fromList xs :: U.Vector Double) :: U.Vector Int

main :: IO ()
main = do
  benchRawData <- B.readFile "bench.dat"
  let benchData bs = force $ map (\p -> B.take bs (B.drop p benchRawData)) [0,bs..B.length benchRawData-1]
  stats $ map L.length (rollsplitL' $ benchData 4096)
  stats $ map L.length (rollsplitL' $ benchData 65536)
  defaultMain [
    let dat = benchData 4096 in dat `seq` bench "simple 4096" $ nf rollsplitL' dat,
    let dat = benchData 65536 in dat `seq` bench "simple 65536" $ nf rollsplitL' dat,
    let dat = benchData 1048576 in dat `seq` bench "simple 1048576" $ nf rollsplitL' dat,
    let dat = benchData 4096 in dat `seq` bench "clamped 4096" $ nf (P.toList . rollsplit 1024 16386 . each)  dat,
    let dat = benchData 65536 in dat `seq` bench "clamped 65536" $ nf (P.toList . rollsplit 1024 16386 . each) dat
    ]

max_stars_repo_name: aristidb/datastorage

max_issues_repo_name: aristidb/datastorage

max_forks_repo_name: aristidb/datastorage


text:
-- |
-- = Single-node reservoir
--
-- In this project, we exploit an established analogy between
-- spatially extended and delay systems (refs. 1-3). That makes possible
-- to employ DDEs as a reservoir substrate.
--
-- 1. Arecchi, F. T., et al. “Two-Dimensional Representation of
--    a Delayed Dynamical System.” Physical Review A, vol. 45, no. 7,
--    Jan. 1992, doi:10.1103/physreva.45.r4225.
-- 2. Appeltant, L., et al. “Information Processing Using a Single
--    Dynamical Node as Complex System.” Nature Communications, vol. 2,
--    2011, p. 468., doi:10.1038/ncomms1476.
-- 3. Virtual Chimera States for Delayed-Feedback Systems - Laurent Larger,
--    Bogdan Penkovsky, Yuri Maistrenko. Physical Review Letters - 08 / 2013.

{-# LANGUAGE BangPatterns #-}
module RC.NTC.Reservoir
  ( Reservoir (..)
  , genReservoir
  ) where

import           Numeric.LinearAlgebra
import qualified Data.Vector.Storable as V
import qualified Linear.V2 as V2
import qualified Numeric.DDE as DDE
import qualified Numeric.DDE.Model as DDEModel

import           RC.NTC.Types
import qualified RC.Helpers as H

-- | Substrate-specific low-level reservoir implementation
genReservoir :: DDEModel.RC -> Reservoir
genReservoir par@DDEModel.RC {
    DDEModel._filt = DDEModel.BandpassFiltering { DDEModel._tau = tau }
  } = Reservoir _r
  where
    _r sample = H.unflatten' nodes responseX
      where
        oversampling = 1 :: Int  -- No oversampling
        detuning = 1.0 :: Double  -- Delay detuning factor, 1 = no detuning
        nodes = rows sample
        delaySamples = round $ detuning * fromIntegral (oversampling * nodes)

        -- Matrix to timetrace
        trace1 = H.flatten' sample

        -- Duplicate the last element (DDE.integHeun2_2D consumes one extra input)
        trace = trace1 V.++ V.singleton (V.last trace1)

        -- Empirically chosen integration time step:
        -- twice faster than the system response time tau
        hStep = tau / 2

        (_, response) = DDE.integHeun2_2D [delaySamples] hStep (DDEModel.bandpassRhs par) (DDE.Input trace)
        responseX = V.map (\(V2.V2 x _) -> x) response

max_stars_repo_name: masterdezign/rc

max_issues_repo_name: masterdezign/rc

max_forks_repo_name: masterdezign/rc


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE TemplateHaskell   #-}
{-# LANGUAGE BangPatterns      #-}
{-| A simple, full dictionary, which maintains:

(1) A bijection of words and IDs (one to one);
(2) A count of every word
(3) Methods for keeping memory use under control using filtering

ID 0 is special in this implementation and represents the count and response of
the unknown word. This has some advantages:

(1) Remapping words after they are filtered is easy: map them to unknown, aka 0
(2) The lookup function is total (no undefined or Nothing to handle)
(3) It serves as a nice default (it looks like background noise)
(4) If you don't want that information, use (Vector.drop 1 vec) to eliminate it,
    or use ((Vector.//) vec [(0,0)]) to 0 it out.
|-}
module NLP.Albemarle.Dict
  ( -- * Using Dictionaries
  Dict(..)
  , counts
  , ids
  , countOf
  , idOf
  , select
  , selectMatrix
  , shift
  , size
  , asSparseMatrix
  , asSparseVector
  , asDenseVector

  -- * Creating Dictionaries
  , dictify
  , dictifyAllWords
  , dictifyFirstWords
  , dictifyDefault

  -- * Modifying Dictionaries
  , filterDict
  ) where
import ClassyPrelude
import Control.Monad.ST
import NLP.Albemarle
import qualified Data.Vector.Instances
import qualified Data.Vector.Generic as Vec
import Data.Vector.Generic ((//))
import qualified Data.Vector as BVec
import qualified Data.Vector.Unboxed as UVec
import qualified Data.Vector.Algorithms.Search as VSearch
import qualified Data.Vector.Algorithms.Intro as VSort
import qualified Data.List.Ordered as OList

import Data.Text (Text)
import qualified Data.Text as Text
import Data.HashSet (HashSet)
import qualified Data.HashSet as HashSet
import Data.HashMap.Strict (HashMap)
import qualified Data.HashMap.Strict as HashMap
import System.IO.Streams (InputStream, OutputStream)
import qualified System.IO.Streams as Streams
import Numeric.LinearAlgebra (Matrix)
import qualified Numeric.LinearAlgebra as HMatrix
import qualified Numeric.LinearAlgebra.Data as HMatrix
import Lens.Micro.TH
import Lens.Micro
import Data.Tuple (swap)
import Data.List.Split

data Dict = Dict {
  _counts :: !(HashMap Text Int),
  _ids :: !(BVec.Vector Text)
} deriving (Show, Eq)

makeLenses ''Dict

instance Monoid Dict where
  mempty = Dict {_counts=mempty, _ids=Vec.singleton ""}
  -- | Union two Dicts (a synonym for mappend). Dict's of huge corpora can be a
  --   memory hog, so consider using filterDict occasionally to manage that.
  mappend left right = Dict {
      _counts = HashMap.unionWith (+) (left^.counts) (right^.counts),
      _ids = Vec.fromList
        $ OList.nub
        $ Vec.toList
        $ Vec.modify VSort.sort (left^.ids <> right^.ids)
    }
instance Semigroup Dict where
  (<>) = mappend

sortNewV v = Vec.modify VSort.sort $ Vec.fromList v

-- | Filter the dictionary: remove stopwords and hepaxes.
filterDict :: Int   -- ^ the minimum word occurences (to remove rare words)
  -> Float  -- ^ the maximum proportion (0,1] of a word (to remove stopwords)
  -> Int    -- ^ the maximum number of unique words total (to cap memory)
  -> Dict   -- ^ the original dictionary
  -> Dict   -- ^ new dictionary
filterDict min_thresh max_fraction max_words dict = let
  max_thresh = floor
    $ (max_fraction *)
    $ fromIntegral
    $ sum
    $ HashMap.elems
    $ dict^.counts
  new_counts =
    HashMap.fromList
    -- This makes sure we preserve the unknown if we can, otherwise 0.
    $ cons ("", HashMap.lookupDefault 0 "" $ dict^.counts)
    $ take max_words
    $ sortOn (negate.snd)
    $ HashMap.toList
    $ HashMap.filter (\v -> min_thresh<=v && v<=max_thresh)
    $ dict^.counts
  in Dict {
    _counts = new_counts,
    _ids = sortNewV $ HashMap.keys new_counts
  }

-- | Get from a dictionary how many times a word was encountered.
--   What exactly this means depends on if you train the dictionary with
--   dictifyAllWords or dictifyFirstWords.
--   Unknown words return a 0. Keep in mind that on account of filterDict this
--   may not mean it doesn't exist in the corpus; it might have been trimmed
--   because it was too common, too rare, or there wasn't enough memory
countOf :: Dict -> Text -> Int
countOf dict word = fromMaybe 0 $ HashMap.lookup word $ dict^.counts

-- | Get the ID of a word from a Dictionary. Unknown words get a 0 (by design).
idOf :: Dict -> Text -> Int
idOf Dict{_ids=vec} word = runST $ do
  mvec <- Vec.thaw vec
  idx <- VSearch.binarySearch mvec word -- Why is this not pure?!
  return $! maybe 0 (\it-> if it == word then idx else 0) $ (Vec.!?) vec idx

-- | Make a sparse vector representing one document using this dict.
--
--   The length of the vector is the size of the vocabulary, plus one. Keep in
--   mind that the length is determined at runtime so pay attention that you do
--   not use vectors from different dictionaries together lest you cause
--   runtime errors.
asSparseVector :: Dict -> [Text] -> SparseVector
asSparseVector dict = SparseVector (Vec.length $ dict^.ids)
  . sort -- Strictly speaking this isn't always necessary but it is for LSA
  -- . cons (0,0.0) -- This shows we count from 0 (it just happens to be 0.0 there)
  . HashMap.toList
  . HashMap.fromListWith (+)
  . fmap (\w-> (idOf dict w, 1.0))

-- | Make a sparse matrix representing one document using this dict.
--
--   The width of the matrix is the size of the vocabulary, plus one. Keep in
--   mind that the length is determined at runtime so pay attention that you do
--   not use matrices from different dictionaries together lest you cause
--   runtime errors.
--
--   This is equavalent to but faster than repeated SparseVectors because it
--   builds a Hashmap first [O(v)] and then uses that [O(1)/lookup] rather than
--   doing no preparation [O(0)] and repeatedly using binary search [O(log v)]
--   Usually this is a 10-fold or better speedup.
asSparseMatrix :: Dict -> [[Text]] -> SparseMatrix
asSparseMatrix dict = let
  toid word = HashMap.lookupDefault 0 word           -- this evalues many times
    $ HashMap.fromList $ zip (v2l $ dict^.ids) [0..] -- this evaluates once
  tovec = SparseVector (Vec.length $ dict^.ids)
    . sort -- Strictly speaking this isn't always necessary but it is for LSA
    -- . cons (0,0.0) -- This shows we count from 0 (it just happens to be 0.0 there)
    . HashMap.toList
    . HashMap.fromListWith (+)
    . fmap (\w-> (toid w, 1.0))
  in SparseMatrix (Vec.length $ dict^.ids)
    . fmap tovec

-- | Make a dense vector representing one document using this dict.
--
--   The length of the vector is the size of the vocabulary, plus one. Keep in
--   mind that the length is determined at runtime so pay attention that you do
--   not use vectors from different dictionaries together lest you cause
--   runtime errors.
asDenseVector :: Dict -> [Text] -> DenseVector
asDenseVector dict =
  Vec.accum (+) (Vec.replicate (Vec.length $ dict^.ids) 0)
  . fmap (\w-> (idOf dict w, 1.0))

-- | O(n + v log v) Create a dictionary from a single document
--   (also consider the [[Text]] -> Dictionary alternatives, they may be faster)
dictify :: [Text] -> Dict
dictify wds = let
  cts = HashMap.fromListWith (+)
    $ cons ("", 0) -- This is to preserve the 0-is-nothing property
    $ (\x -> (x, 1)) <$> wds
  in Dict {
    _counts = cts,
    _ids = sortNewV $ HashMap.keys cts
  }

-- | Create a dictionary from many documents, where every instance of every word
--   counts. (As opposed to dictifyFirstWords.) Might be faster than FirstWords.
dictifyAllWords :: [[Text]] -> Dict
dictifyAllWords = dictify . mconcat

-- | Create a dictionary from many documents, counting only the first instance
--   of each word in a given document. Might be slower than AllWords.
dictifyFirstWords :: [[Text]] -> Dict
dictifyFirstWords = dictifyAllWords . fmap (HashSet.toList . HashSet.fromList)

-- | Use some sensible defaults for making large dictionaries:
-- handle 10k documents at a time, remove word with <= 2 instances in 10k docs,
-- remove words in >= 50% of documents, and limit to 1000000 unique words.
dictifyDefault :: [[Text]] -> Dict
dictifyDefault =
  foldl' (\l r -> filterDict 2 0.5 1000000 $ l <> r) mempty
  . fmap (dictifyAllWords)
  . chunksOf 10000

-- | O(1) Retrieve the number of unique words in the vocabulary
size :: Dict -> Int
size dict = HashMap.size $ dict^.counts

-- | O(v) Retrieve the total number of (non-unique) words processed
support :: Dict -> Int
support dict = HashMap.foldl' (+) 0 $ dict^.counts

-- Just for internal use. It's polymorphic, so don't use the eta reduction
v2l a = Vec.toList a
l2v a = Vec.fromList a

-- | What indices into the old matrix will generate an approximation of the new
--   matrix? This is intended for something like
-- > ((HMatrix.?) old_matrix $ select remap) == something like new_matrix
--
--   It's useful for merging termdoc-matrices, SVD's, or the like.
select :: Dict -> Dict -> [Int]
select left right = let
  l = zip [0..] $ v2l $ left^.ids
  r = v2l $ right^.ids
  in select' r l -- backward!

-- | Compute a variant on indirect sorted set intersection. Only used for `select`
select' :: Ord a
  => [a] -- ^ Sorted list X
  -> [(Int, a)] -- ^ Sorted list Y
  -> [Int] -- ^ if X[i] in Y: then Y[Z[i]] = X[i]; else Z[i] = 0
select' [] _  = []
select' ox [] = map (const 0) ox -- Must stay in lock step with X
select' ox@(x:xs) oy@((iy,y):ys) = case compare x y of
  LT -> 0  : select' xs oy
  EQ -> iy : select' xs ys
  GT ->      select' ox ys

-- | Rearrange rows in a matrix to match the order in a new matrix, after a
--   dictionary merge or filter event.
selectMatrix :: Dict -> Dict -> HMatrix.Matrix Double -> HMatrix.Matrix Double
selectMatrix left right old_matrix = (HMatrix.?) old_matrix $ select left right

-- | Given a Remap, convert the original word ID into the new ID
--   (might be more useful curried)
shift :: Dict -> Dict -> Int -> Int
shift left right i = HashMap.lookupDefault 0 i
  $ HashMap.fromList
  $ find'
  (zip [0..] (v2l $ left^.ids))
  (zip [0..] (v2l $ right^.ids))

-- | Compute an indirect sorted set intersection. Only used for `shift`
find' :: Ord a
  => [(Int, a)] -- ^ Left sorted list X
  -> [(Int, a)] -- ^ Right sorted list Y
  -> [(Int, Int)] -- ^ Indices into X, Y representing the sorted intersection
find' [] new  = []
find' orig [] = []
find' ox@((ix,x):xs) oy@((iy,y):ys) = case compare x y of
  LT -> find' xs oy
  EQ -> (ix, iy) : find' xs ys
  GT -> find' ox ys

max_stars_repo_name: SeanTater/albemarle

max_issues_repo_name: SeanTater/albemarle

max_forks_repo_name: SeanTater/albemarle