text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
module AI.Singularity.Examples.Mnist where
import GHC.TypeLits
import qualified Numeric.LinearAlgebra.Data as LD
import qualified Numeric.LinearAlgebra.Devel as LDev
import qualified Numeric.LinearAlgebra as LA
import qualified Control.Monad.State as S
import Control.Lens
import Data.Proxy
import Control.Monad.ST
import System.Random
import Control.Monad
import Control.Arrow
import qualified Data.Vector.Storable as VS

import AI.Singularity.Network

getLabels :: String -> IO [Float]
getLabels = map (fromIntegral . fromEnum) . drop 8 . BS.unpack <$< BS.readFile

getImages2 :: String -> IO [[Vector 28]]
getImages2 = (map . map) fromListV . splitByLen 28 . splitByLen 28 . map ((/255). fromIntegral . fromEnum) . drop 16 . BS.unpack <$< BS.readFile

getSet3 :: String -> IO [((Vector 400, Vector 384), Vector 10)]
getSet3 = map (toVec >>> first splitV) <$< getSet

getSet2 :: String -> IO [([Vector 28],Vector 10)]
getSet2 nm = zip <$> getImages2 ("mnist/"++nm++"-images-idx3"++rest) <*> (map wynik <$> getLabels ("mnist/"++nm++"-labels-idx1"++rest))
  where rest          = "-ubyte"
        wynik liczba  = Vec . VS.map (\x -> if x == liczba then 1 else 0) . VS.fromList $ [0..9]

getImages :: String -> IO [VS.Vector Float]
getImages = map VS.fromList . splitByLen (28*28) . map ((/255). fromIntegral . fromEnum) . drop 16 . BS.unpack <$< BS.readFile

-- getImagesP :: FilePath -> Producer PB.ByteString IO a
getImagesC fn = let s = do
                      x <- CB.take 28
                      if BS.length x == 0
                        then return ()
                        else yield (map ((/255) . fromIntegral . fromEnum) . BS.unpack $ x) >> s
                    v = do
                      x <- CL.take 28
                      if null x
                        then return ()
                        else yield x >> v
                in CB.sourceFile fn .| (CB.drop 16 >> s) .| CL.map (fromListV :: [Float] -> Vector 28) .| v

getLinearC fn = let s = do
                      x <- BS.unpack <$> CB.take (28*28)
                      if null x
                        then return ()
                        else do
                              yield . map ((/255). fromIntegral . fromEnum) $ x
                              s
                in CB.sourceFile fn .| (CB.drop 16 >> s) .| CL.map (fromListV :: [Float] -> Vector 784)

fix f = let x = f x in x

-- setConduit :: String -> Source
setConduit nm = zipSources (getImagesC im) (getLabelsC lab)
  where im   = "mnist/"++nm++"-images-idx3"++rest
        lab  = "mnist/"++nm++"-labels-idx1"++rest
        rest = "-ubyte"

setConduit2 nm = zipSources (getLinearC im) (getLabelsC lab)
  where im   = "mnist/"++nm++"-images-idx3"++rest
        lab  = "mnist/"++nm++"-labels-idx1"++rest
        rest = "-ubyte"

getLabelsC fn =
  CB.sourceFile fn .|
  (CB.drop 8 >> fix (\f -> CB.head >>= \case { Nothing -> return (); Just x -> yield x >> f} )) .| CL.map (fromEnum >>> wynik)
    where wynik n = Vec . VS.map (\x -> if x == n then 1 else 0) . VS.fromList $ [0..9] :: Vector 10

  --map (fromIntegral . fromEnum) . drop 8 . BS.unpack

getSet :: String -> IO [(VS.Vector Float, Float)]
getSet nm = zip <$> getImages ("mnist/"++nm++"-images-idx3"++rest) <*> getLabels ("mnist/"++nm++"-labels-idx1"++rest)
  where rest = "-ubyte"

printSet :: (VS.Vector Float, Float) -> IO ()
printSet (v,a) = do
  print a
  mapM_ (\(i,val) -> (if i `mod` 28 == 0 then putStr (printVal val ++ "\n") else putStr (printVal val) )) . zip [1..] . VS.toList $ v

printVal :: Float -> String
printVal x = if x < 0.00001 then "  " else "* "

toVec :: (VS.Vector Float, Float) -> (Vector 784, Vector 10)
toVec (!obraz, !liczba) = (Vec obraz, Vec wynik)
    where wynik = VS.map (\x -> if x == liczba then 1 else 0) . VS.fromList $ [0..9]

printVec :: Vector 784 -> String
printVec (Vec !v) =  concatMap (\(i,val) -> (if i `mod` 28 == 0 then printVal val ++ "\n" else printVal val )) . zip [1..] . VS.toList $ v

testVec :: (VS.Vector Float, Float) -> (Vector 784, Int)
testVec (!obrac, !liczba) = (Vec obrac, round liczba)

max_stars_repo_name: Antystenes/Memetic-Predictor

max_issues_repo_name: Antystenes/Memetic-Predictor

max_forks_repo_name: Antystenes/Memetic-Predictor


text:
-----------------------------------------------------------------------------
{- |
This module defines the monad of sampling functions. See Park, Pfenning and Thrun:
A probabilistic language based upon sampling functions, Principles of programming languages 2005

Sampling functions allow the composition of both discrete and continuous
probability distributions.

The implementation and interface are similar to those in the random-fu, monte-carlo
and monad-mersenne-random packages.

Example -- a biased coin:

@
data Throw = Head | Tail

throw bias = do
  b <- bernoulli bias
  return $ if b then Head else Tail

tenThrowsCrooked = replicateM 10 $ throw 0.3

countHeads = do
   throws <- tenThrowsCrooked
   return $ length [ () | Head <- throws]

main = do
 print =<< sampleIO tenThrowsCrooked
 print =<< eval ((\<4) \`fmap\` countHeads)
@

-}

{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}

module Math.Probably.Sampler where

import Control.Monad
import Control.Applicative
import qualified Math.Probably.PDF as PDF
import Numeric.LinearAlgebra hiding (find)
import System.Random.Mersenne.Pure64
import System.Environment
import Data.List
import Data.Maybe
import Data.Ord
import Control.Spoon

--import Debug.Trace

type Seed = PureMT

data Sampler a = Sam {unSam :: Seed -> (a, Seed) }
               | Samples [a]

instance Functor Sampler where
    fmap f (Sam sf) = Sam $ \rs -> let (x,rs') = sf rs in
                                   (f x, rs')
    fmap f (Samples xs) = Samples $ map f xs

instance Applicative Sampler where
    pure x = Sam (\rs-> (x, rs))
    (Sam sff) <*> (Sam sfx) = Sam $ \rs-> let (f ,rs') = sff rs
                                              (x, rs'') = sfx rs' in
                                          (f x, rs'')

instance Monad Sampler where
    return = pure
    (Sam sf) >>= f = Sam $ \rs-> let (x, rs'::Seed) = sf rs
                                     nextProb = f x
                                 in case nextProb of
                                      Sam g -> g rs'
                                      Samples xs -> primOneOf xs rs'
    (Samples xs) >>= f = Sam $ \rs-> let (x, rs'::Seed) = primOneOf xs rs
                                         nextProb = f x
                                     in case nextProb of
                                          Sam g -> g rs'
                                          Samples ys -> primOneOf ys rs'

-- | given a seed, return an infinite list of draws from sampling function
runSampler :: Seed -> Sampler a -> [a]
runSampler pmt s@(Sam sf)
   = let (x, pmt') = sf pmt
     in x : runSampler pmt' s
runSampler _ (Samples xs) = xs

-- | Get a seed
getSeedIO :: IO Seed
getSeedIO = do
   args <- getArgs
   case mapMaybe (stripPrefix "--seed=") args of
      [] ->  newPureMT
      sdStr:_ -> return $ pureMT $ read sdStr

-- | Return an infinite list of draws from sampling function in the IO monad
runSamplerIO :: Sampler a -> IO [a]
runSamplerIO s =
   fmap (`runSampler` s) $ getSeedIO

-- | Return a singe draw from sampling function
sampleIO :: Sampler a -> IO a
sampleIO s = head `fmap` runSamplerIO s

-- | Return a list of n draws from sampling function
sampleNIO :: Int -> Sampler a -> IO [a]
sampleNIO n s = take n `fmap` runSamplerIO s

-- | Estimate the probability that a hypothesis is true (in the IO monad)
eval :: Sampler Bool -> Sampler Double
eval s = do
  bs <- replicateM 1000 s
  return $ realToFrac (length (filter id bs)) / 1000

{-mu :: Vector Double
sigma :: Matrix Double

mystery = -1
mu = fromList [0,0,0]
sigma =  (3><3) [ 1,    1,       0,
                  1,    1,     mystery,
                  0,  mystery,   1]

samIt = sampleNIO 2 $ multiNormal mu sigma
 -}
-- | The joint distribution of two independent distributions
joint :: Sampler a -> Sampler b -> Sampler (a,b)
joint sf1 sf2 = liftM2 (,) sf1 sf2

-- | The joint distribution of two distributions where one depends on the other
jointConditional :: Sampler a -> (a-> Sampler b) -> Sampler (a,b)
jointConditional sf1 condsf
    = do x <- sf1
         y <- condsf x
         return (x,y)

--replicateM :: Monad m => Int -> m a -> m [a]
--replicateM n ma = forM [1..n] $ const ma

-- * Uniform distributions

-- | The unit interval U(0,1)
unitSample :: Sampler Double
unitSample = Sam randomDouble

-- | for x and y, the uniform distribution between x and y
uniform :: (Fractional a) => a -> a -> Sampler a
uniform a b = (\x->(realToFrac x)*(b-a)+a) `fmap` unitSample

-- * Normally distributed sampling function

--http://en.wikipedia.org/wiki/Box-Muller_transform
-- | The univariate gaussian (normal) distribution defined by mean and standard deviation
gauss :: (Floating b) => b -> b -> Sampler b
gauss m sd =
    do (u1,u2) <- (mapPair realToFrac) `fmap` joint unitSample unitSample
       return $ sqrt(-2*log(u1))*cos(2*pi*u2)*sd+m
  where mapPair f (x,y) = (f x, f y)

-- | Gaussians specialised for doubles
gaussD :: Double -> Double -> Sampler Double
gaussD m sd =
    do (u1,u2) <- joint unitSample unitSample
       return $ sqrt(-2*log(u1))*cos(2*pi*u2)*sd+m

gaussMany :: Floating b => [(b,b)] -> Sampler [b]
gaussMany means_sds = do gus <- gaussManyUnit (length means_sds)
                         return $ map f $ zip gus means_sds
    where f (gu, (mean, sd)) = gu*sd+mean

gaussManyD :: [(Double,Double)] -> Sampler [Double]
gaussManyD means_sds = do gus <- gaussManyUnitD (length means_sds)
                          return $ zipWith f gus means_sds
    where f gu (mean, sd) = gu*sd+mean

gaussManyUnit :: Floating b => Int -> Sampler [b]
gaussManyUnit 0 = return []
gaussManyUnit n | odd n = liftM2 (:) (gauss 0 1) (gaussManyUnit (n-1))
                | otherwise = do us <- forM [1..n] $ const $ unitSample
                                 return $ gaussTwoAtATime $ map realToFrac us
  where
    gaussTwoAtATime :: Floating a =>  [a] -> [a]
    gaussTwoAtATime (u1:u2:rest) = sqrt(-2*log(u1))*cos(2*pi*u2) : sqrt(-2*log(u1))*sin(2*pi*u2) : gaussTwoAtATime rest
    gaussTwoAtATime _ = []

gaussManyUnitD :: Int -> Sampler [Double]
gaussManyUnitD 0 = return []
gaussManyUnitD n | odd n = liftM2 (:) (gauss 0 1) (gaussManyUnit (n-1))
                | otherwise = do us <- forM [1..n] $ const $ unitSample
                                 return $ gaussTwoAtATimeD us
 where
   gaussTwoAtATimeD :: [Double] -> [Double]
   gaussTwoAtATimeD (u1:u2:rest) = sqrt(-2*log(u1))*cos(2*pi*u2) : sqrt(-2*log(u1))*sin(2*pi*u2) : gaussTwoAtATimeD rest
   gaussTwoAtATimeD _ = []

-- | Multivariate normal distribution
multiNormal :: Vector Double -> Matrix Double -> Sampler (Vector Double)
multiNormal mu sigma =
  let c =   cholSH sigma
      a = trans c
      k = dim mu
  in do z <- fromList `fmap` gaussManyUnitD k
--        return $ mu + (head $ toColumns $ a*asRow z)
        let c = asColumn z
        let r = asRow z
        return $ (mu + (head $ toColumns $ a `multiply` asColumn z))

multiNormalByChol :: Vector Double -> Matrix Double -> Sampler (Vector Double)
multiNormalByChol mu cholSigma =
  let a = trans $ cholSigma
      k = dim mu
  in do z <- fromList `fmap` gaussManyUnitD k
--        return $ mu + (head $ toColumns $ a*asRow z)
        let c = asColumn z
        let r = asRow z
        return $ (mu + (head $ toColumns $ a `multiply` asColumn z))

multiNormalIndep  :: Vector Double -> Vector Double -> Sampler (Vector Double)
multiNormalIndep vars mus = do
   let k = dim mus
   gs <- gaussManyUnitD k
   return $ fromList $ zipWith3 (\var mu g -> g*sqrt(var) + mu) (toList vars) (toList mus) gs

--http://en.wikipedia.org/wiki/Log-normal_distribution#Generating_log-normally-distributed_random_variates

-- | log-normal distribution <http://en.wikipedia.org/wiki/Log-normal_distribution>
logNormal :: (Floating b) => b -> b -> Sampler b
logNormal m sd =
    do n <- gauss 0 1
       return $ exp $ m + sd * n

-- * Other distribution

-- | Bernoulli distribution. Returns a Boolean that is 'True' with probability 'p'
bernoulli :: Double -> Sampler Bool
bernoulli p = (<p) `fmap` unitSample

discrete :: [(Double,a)] -> Sampler a
discrete weightedSamples =
   let sumWeights = sum $ map fst weightedSamples
       cummWeightedSamples = scanl (\(csum,_) (w,x) -> (csum+w,x)) (0,undefined) $ sortBy (comparing fst) weightedSamples
   in do u <- unitSample
         return . snd . fromJust $ find ((>=u*sumWeights) . fst) cummWeightedSamples

primOneOf :: [a] -> Seed -> (a, Seed)
primOneOf xs seed
   = let (u, nextSeed) = randomDouble seed
         idx = floor $ (realToFrac u)*(realToFrac $ length xs )
     in (xs !! idx, nextSeed)

oneOf :: [a] -> Sampler a
oneOf xs = do idx <- floor `fmap` uniform (0::Double) (realToFrac $ length xs )
              return $ xs !! idx

nOf :: Int -> [a] -> Sampler [a]
nOf n xs = sequence $ replicate n $ oneOf xs

{-main = do
  rnds <- take 1000 `fmap` runSamplerSysRan (oneOf [1,2,3])
  let diff = sort $ nub rnds
  print $ map (\x->(x, length $ filter (==x) rnds )) $ diff -}

-- | Bayesian inference from likelihood and prior using rejection sampling.
bayesRejection :: (PDF.PDF a) -> Double -> Sampler a -> Sampler a
bayesRejection p c q = bayes
    where bayes = do x <- q
                     u <- unitSample
                     if u < p x / c
                        then return x
                        else bayes

{-
expectation :: Fractional a =>  Int -> Sampler a -> IO a
expectation n sf =
    (mean . take n) `fmap` runSamplerIO sf

expectSD :: Floating a =>  Int -> Sampler a -> IO (a,a)
expectSD n sf =
    (meanSD . take n) `fmap` runSamplerIO sf
-}

--mapPair :: (a->b) -> (a,a) -> (b,b)

{-
--http://cgi.cse.unsw.edu.au/~dons/blog/2008/05/16#fast
mean :: Fractional a =>  [a] -> a
mean = go 0 0
        where
            -- go ::  -> Int -> [Double] -> Double
            go s n []     = s / fromIntegral n
            go !s !n (x:xs) = go (s+x) (n+1) xs

meanSD :: Floating a => [a] -> (a,a)
meanSD = go 0 0 0
    where go sq s n [] = let len = fromIntegral n in
                         (s/len, (recip len)*sqrt (len*sq-s*s))
          go !sq !s !n (x:xs) = go (sq+x*x) (s+x) (n+1) xs

test_mean = mean [0..1e8]
test_meanVar = meanSD [0..1e8]

main = do u <- expectSD 1000000 $ gauss 0 1
          print u

-}

--poisson ::  ::  Double -> [Double] -> IO Double
-- | Exponential distribution
expDist rate =  (\u-> negate $ (log(1-u))/rate) `fmap` unitSample

-- | Homogeneous poisson process defined by rate and duration
poissonMany :: Double -> Double -> Sampler [Double]
poissonMany rate tmax = aux 0
    where aux last = do
            next <- (+last) `fmap` expDist rate
            if next > tmax
               then return []
               else liftM2 (:) (return next) $ aux next

-- | binomial distribution
binomial :: Int -> Double -> Sampler Int
binomial n p = do
  bools <- forM [1..n] $ const $ fmap (<p) unitSample
  return $ length $ [t | t@True <- bools]

-- from random-fu
-- | Gamma distribution
gamma :: Double -> Double -> Sampler Double
gamma a b
     | a < 1
    = do
        u <- unitSample
        x <- gamma (1 + a) b
        return (x * u ** recip a)
    | otherwise
    = go
        where
            d = a - (1 / 3)
            c = recip (3 * sqrt d) -- (1 / 3) / sqrt d

            go = do
                x <- gaussD 0 1

                let cx = c * x
                    v = (1 + cx) ^ 3

                    x_2 = x * x
                    x_4 = x_2 * x_2

                if cx <= (-1)
                    then go
                    else do
                        u <- unitSample

                        if         u < 1 - 0.0331 * x_4
                            || log u < 0.5 * x_2  + d * (1 - v + log v)
                            then return (b * d * v)
                            else go

-- | inverse gamma distribution
invGamma :: Double -> Double -> Sampler Double
invGamma a b = recip `fmap` gamma a b

--http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution
--multiNormal :: Vector Double -> Matrix Double -> Sampler (Vector Double)

--http://www.xycoon.com/beta_randomnumbers.htm
-- | beta distribution
beta :: Int -> Int -> Sampler Double
beta a b =
    let gam n = do us <- forM [1..n] $ const unitSample
                   return $ log $ product us
    in do gama1 <- gamma (realToFrac a) 1
--          gama2 <- gamma (realToFrac a) 1

          gamb <- gamma (realToFrac b) 1
          return $ gama1/(gama1+gamb)

tbeta = sampleNIO 100 $ beta 1 1

ext: hs

max_stars_repo_name: glutamate/probably

max_issues_repo_name: glutamate/probably

max_forks_repo_name: glutamate/probably


text:
{-# LANGUAGE FlexibleContexts #-}
module FokkerPlanck.DomainChange where

import           Data.Array.Repa       as R
import           Data.Complex
import           Data.List             as L
import           FokkerPlanck.Types
import           Numeric.LinearAlgebra as NL
import           Types
import           Utils.Array

{-# INLINE r2z1t0Tor2s1t0 #-}
r2z1t0Tor2s1t0 :: Int -> [Double] -> R2Z1T0Array -> R2S1T0Array
r2z1t0Tor2s1t0 numOrientations freqs arr =
  let deltaTheta = 2 * pi / (fromIntegral numOrientations)
      (Z :. numFreqs :. numT0Freqs :. xLen :. yLen) = extent arr
      mat1 =
        (numOrientations >< numFreqs) . R.toList $
        R.traverse
          (fromListUnboxed (Z :. numFreqs) freqs)
          (const (Z :. numOrientations :. numFreqs))
          (\f (Z :. i :. j) ->
             exp (0 :+ (-1) * (deltaTheta * fromIntegral i) * f (Z :. j)))
      mat2 = (numFreqs >< (xLen * yLen * numT0Freqs)) . R.toList $ arr
   in fromListUnboxed (Z :. numOrientations :. numT0Freqs :. xLen :. yLen) .
      NL.toList . flatten $
      mat1 NL.<> mat2

{-# INLINE r2z1Tor2s1 #-}
r2z1Tor2s1 ::
     (Source s (Complex Double))
  => Int
  -> [Double]
  -> R.Array s DIM3 (Complex Double)
  -> R.Array U DIM3 (Complex Double)
r2z1Tor2s1 numOrientations freqs arr =
  let deltaTheta = 2 * pi / (fromIntegral numOrientations)
      (Z :. numFreqs :. xLen :. yLen) = extent arr
      mat1 =
        (numOrientations >< numFreqs) . R.toList $
        R.traverse
          (fromListUnboxed (Z :. numFreqs) freqs)
          (const (Z :. numOrientations :. numFreqs))
          (\f (Z :. i :. j) ->
             exp (0 :+ (-1) * (deltaTheta * fromIntegral i) * f (Z :. j)))
      mat2 = (numFreqs >< (xLen * yLen)) . R.toList $ arr
   in fromListUnboxed (Z :. numOrientations :. xLen :. yLen) .
      NL.toList . flatten $
      mat1 NL.<> mat2

{-# INLINE r2s1tor2z1 #-}
r2s1tor2z1 ::
     (Source s (Complex Double))
  => [Double]
  -> R.Array s DIM3 (Complex Double)
  -> R.Array U DIM3 (Complex Double)
r2s1tor2z1 freqs arr =
  let (Z :. numOrientations :. xLen :. yLen) = extent arr
      deltaTheta = 2 * pi / (fromIntegral numOrientations)
      numFreqs = L.length freqs
      mat1 =
        (numFreqs >< numOrientations) . R.toList $
        R.traverse
          (fromListUnboxed (Z :. numFreqs) freqs)
          (const (Z :. numFreqs :. numOrientations))
          (\f (Z :. i :. j) ->
             exp (0 :+ (deltaTheta * fromIntegral j) * f (Z :. i)))
      mat2 = (numOrientations >< (xLen * yLen)) . R.toList $ arr
   in fromListUnboxed (Z :. numFreqs :. xLen :. yLen) . NL.toList . flatten $
      mat1 NL.<> mat2

{-# INLINE r2z2t0s0Tor2s1rpt0s0 #-}
r2z2t0s0Tor2s1rpt0s0 ::
     Int -> [Double] -> Int -> [Double] -> R2Z2T0S0Array -> R2S1RPT0S0Array
r2z2t0s0Tor2s1rpt0s0 numOrientations thetaFreqs numScales scaleFreqs arr =
  let deltaTheta = 2 * pi / (fromIntegral numOrientations)
      (Z :. numThetaFreq :. numScaleFreq :. numTheta0Freq :. numScale0Freq :. xLen :. yLen) =
        extent arr
      mat1 =
        ((numOrientations * numScales) >< (numThetaFreq * numScaleFreq)) .
        R.toList $
        R.traverse3
          (fromFunction
             (Z :. numOrientations :. numScales :. numThetaFreq :. numScaleFreq)
             (const 0))
          (fromListUnboxed (Z :. numThetaFreq) thetaFreqs)
          (fromListUnboxed (Z :. numScaleFreq) scaleFreqs)
          (\a _ _ -> a)
          (\f1 f2 f3 (Z :. i :. j :. k :. l) ->
             exp
               (0 :+
                (-1) *
                ((deltaTheta * fromIntegral i) * f2 (Z :. k) +
                 (2 * pi * fromIntegral j / fromIntegral numScales) *
                 f3 (Z :. l))))
      mat2 =
        ((numThetaFreq * numScaleFreq) ><
         (xLen * yLen * numTheta0Freq * numScale0Freq)) .
        R.toList $
        arr
   in fromListUnboxed
        (Z :. numOrientations :. numScales :. numTheta0Freq :. numScale0Freq :.
         xLen :.
         yLen) .
      NL.toList . flatten $
      mat1 NL.<> mat2

{-# INLINE r2z2Tor2s1rp #-}
r2z2Tor2s1rp ::
     (Source s (Complex Double))
  => Int
  -> [Double]
  -> Int
  -> [Double]
  -> Double
  -> R.Array s DIM4 (Complex Double)
  -> R.Array U DIM4 (Complex Double)
r2z2Tor2s1rp numOrientations thetaFreqs numScales scaleFreqs maxScale arr =
  let deltaTheta = 2 * pi / (fromIntegral numOrientations)
      deltaScale = log maxScale / (fromIntegral numScales)
      (Z :. numThetaFreq :. numScaleFreq :. xLen :. yLen) = extent arr
      mat1 =
        ((numOrientations * numScales) >< (numThetaFreq * numScaleFreq)) .
        R.toList $
        R.traverse3
          (fromFunction
             (Z :. numOrientations :. numScales :. numThetaFreq :. numScaleFreq)
             (const 0))
          (fromListUnboxed (Z :. numThetaFreq) thetaFreqs)
          (fromListUnboxed (Z :. numScaleFreq) scaleFreqs)
          (\a _ _ -> a)
          (\f1 f2 f3 (Z :. i :. j :. k :. l) ->
             exp
               (0 :+
                (1) *
                ((deltaTheta * fromIntegral i) * f2 (Z :. k) +
                 (2 * pi * fromIntegral j * deltaScale / (log maxScale)) *
                 f3 (Z :. l))))
      mat2 = ((numThetaFreq * numScaleFreq) >< (xLen * yLen)) . R.toList $ arr
  in fromListUnboxed (Z :. numOrientations :. numScales :. xLen :. yLen) .
     NL.toList . flatten $
     mat1 NL.<> mat2

{-# INLINE r2z2Tor2s1rpP #-}
r2z2Tor2s1rpP ::
     (Source s (Complex Double))
  => Int
  -> [Double]
  -> Int
  -> [Double]
  -> Double
  -> R.Array s DIM4 (Complex Double)
  -> IO (R.Array U DIM4 (Complex Double))
r2z2Tor2s1rpP numOrientations thetaFreqs numScales scaleFreqs maxR arr = do
  let deltaTheta = 2 * pi / (fromIntegral numOrientations)
      deltaScale = log maxR / (fromIntegral numOrientations)
      (Z :. numThetaFreq :. numScaleFreq :. xLen :. yLen) = extent arr
      mat2 = ((numThetaFreq * numScaleFreq) >< (xLen * yLen)) . R.toList $ arr
  mat1 <-
    fmap
      (((numOrientations * numScales) >< (numThetaFreq * numScaleFreq)) .
       R.toList) .
    computeUnboxedP $
    R.traverse3
      (fromFunction
         (Z :. numOrientations :. numScales :. numThetaFreq :. numScaleFreq)
         (const 0))
      (fromListUnboxed (Z :. numThetaFreq) thetaFreqs)
      (fromListUnboxed (Z :. numScaleFreq) scaleFreqs)
      (\a _ _ -> a)
      (\f1 f2 f3 (Z :. i :. j :. k :. l) ->
         exp
           (0 :+
            (1) *
            ((deltaTheta * fromIntegral i) * f2 (Z :. k) +
             (2 * pi * fromIntegral j * deltaScale / (log maxR)) * f3 (Z :. l))))
  return .
    fromListUnboxed (Z :. numOrientations :. numScales :. xLen :. yLen) .
    NL.toList . flatten $
    mat1 NL.<> mat2

{-# INLINE r2z2t0s0Tor2s1rps1rp #-}
r2z2t0s0Tor2s1rps1rp ::
     Int -> [Double] -> [Double] -> Int -> [Double] -> [Double] -> Double -> R2Z2T0S0Array -> R.Array U DIM6 (Complex Double)
r2z2t0s0Tor2s1rps1rp numOrientations thetaFreqs theta0Freqs numScales scaleFreqs scale0Freqs maxScale arr =
  let deltaTheta = 2 * pi / (fromIntegral numOrientations)
      deltaScale =
        if maxScale == 0
          then 0
          else (2 * pi / (log maxScale)) / (fromIntegral numScales)
      (Z :. numThetaFreq :. numScaleFreq :. numTheta0Freq :. numScale0Freq :. xLen :. yLen) =
        extent arr
      mat1 =
        ((numOrientations * numScales * numOrientations * numScales) ><
         (numThetaFreq * numScaleFreq * numTheta0Freq * numScale0Freq)) .
        R.toList $
        R.traverse4
          (fromListUnboxed (Z :. numThetaFreq) thetaFreqs)
          (fromListUnboxed (Z :. numScaleFreq) scaleFreqs)
          (fromListUnboxed (Z :. numTheta0Freq) theta0Freqs)
          (fromListUnboxed (Z :. numScale0Freq) scale0Freqs)
          (\_ _ _ _ ->
             (Z :. numOrientations :. numScales :. numOrientations :. numScales :.
              numThetaFreq :.
              numScaleFreq :.
              numTheta0Freq :.
              numScale0Freq))
          (\ft fs ft0 fs0 (Z :. o :. s :. o0 :. s0 :. tf :. sf :. tf0 :. sf0) ->
             exp
               (0 :+
                (-1) *
                ((deltaTheta * fromIntegral o) * ft (Z :. tf) +
                 (deltaTheta * fromIntegral o0) * ft0 (Z :. tf0) +
                 (deltaScale * fromIntegral s) * fs (Z :. sf) +
                 (deltaScale * fromIntegral s0) * fs0 (Z :. sf0))))
      mat2 =
        ((numThetaFreq * numScaleFreq * numTheta0Freq * numScale0Freq) ><
         (xLen * yLen)) .
        R.toList $
        arr
   in fromListUnboxed
        (Z :. numOrientations :. numScales :. numOrientations :. numScales :.
         xLen :.
         yLen) .
      NL.toList . flatten $
      mat1 NL.<> mat2

{-# INLINE r2s1rps1rpTor2z2t0s0 #-}
r2s1rps1rpTor2z2t0s0 ::
     Int
  -> [Double]
  -> [Double]
  -> Int
  -> [Double]
  -> [Double]
  -> Double
  -> R.Array U DIM6 (Complex Double)
  -> R2Z2T0S0Array
r2s1rps1rpTor2z2t0s0 numOrientations thetaFreqs theta0Freqs numScales scaleFreqs scale0Freqs maxScale arr =
  let deltaTheta = 2 * pi / (fromIntegral numOrientations)
      deltaScale =
        if maxScale == 0
          then 0
          else (2 * pi / (log maxScale)) / (fromIntegral numScales)
      numThetaFreq = L.length thetaFreqs
      numScaleFreq = L.length scaleFreqs
      numTheta0Freq = L.length theta0Freqs
      numScale0Freq = L.length scale0Freqs
      (Z :. _ :. _ :. _ :. _ :. xLen :. yLen) =
        extent arr
      mat1 =
        ((numThetaFreq * numScaleFreq * numTheta0Freq * numScale0Freq) ><
         (numOrientations * numScales * numOrientations * numScales)) .
        R.toList $
        R.traverse4
          (fromListUnboxed (Z :. numThetaFreq) thetaFreqs)
          (fromListUnboxed (Z :. numScaleFreq) scaleFreqs)
          (fromListUnboxed (Z :. numTheta0Freq) theta0Freqs)
          (fromListUnboxed (Z :. numScale0Freq) scale0Freqs)
          (\_ _ _ _ ->
             (Z :. numThetaFreq :. numScaleFreq :. numTheta0Freq :.
              numScale0Freq :.
              numOrientations :.
              numScales :.
              numOrientations :.
              numScales))
          (\ft fs ft0 fs0 (Z :. tf :. sf :. tf0 :. sf0 :. o :. s :. o0 :. s0) ->
             exp
               (0 :+ (1) *
                ((deltaTheta * fromIntegral o) * ft (Z :. tf) +
                 (deltaTheta * fromIntegral o0) * ft0 (Z :. tf0) +
                 (deltaScale * fromIntegral s) * fs (Z :. sf) +
                 (deltaScale * fromIntegral s0) * fs0 (Z :. sf0))))
      mat2 =
        ((numOrientations * numScales * numOrientations * numScales) ><
         (xLen * yLen)) .
        R.toList $
        arr
   in fromListUnboxed
        (Z :. numThetaFreq :. numScaleFreq :. numTheta0Freq :. numScale0Freq :.
         xLen :.
         yLen) .
      NL.toList . flatten $
      mat1 NL.<> mat2

-- x y t s
{-# INLINE stsTor2s1rp #-}
stsTor2s1rp ::
     Int
  -> Double
  -> [Double]
  -> Int
  -> Double
  -> [Double]
  -> Int
  -> [Double]
  -> Int
  -> [Double]
  -> Double
  -> R.Array U DIM4 (Complex Double)
  -> IO (R.Array U DIM4 (Complex Double))
stsTor2s1rp numX xLen xFreqs numY yLen yFreqs numOrientations thetaFreqs numScales scaleFreqs maxScale arr = do
  let deltaX = xLen / fromIntegral numX
      deltaY = yLen / fromIntegral numY
      deltaTheta = 2 * pi / fromIntegral numOrientations
      deltaScale = (log maxScale) / fromIntegral numScales
      centerX = div numX 2
      centerY = div numY 2
      vec = NL.fromList . R.toList $ arr
  transformMat <-
    fmap
      (((numX * numY * numOrientations * numScales) ><
        (L.product . listOfShape . extent $ arr)) .
       R.toList) .
    computeUnboxedP .
    R.traverse4
      (fromListUnboxed (Z :. L.length xFreqs) xFreqs)
      (fromListUnboxed (Z :. L.length yFreqs) yFreqs)
      (fromListUnboxed (Z :. L.length thetaFreqs) thetaFreqs)
      (fromListUnboxed (Z :. L.length scaleFreqs) scaleFreqs)
      (\(Z :. numXFreq) (Z :. numYFreq) (Z :. numThetaFreq) (Z :. numScaleFreq) ->
         (Z :. numX :. numY :. numOrientations :. numScales :. numXFreq :.
          numYFreq :.
          numThetaFreq :.
          numScaleFreq)) $ \fx fy fo fs (Z :. x :. y :. o :. s :. xf :. yf :. orif :. sf) ->
      exp
        (0 :+
         (fo (Z :. orif) * fromIntegral o * deltaTheta +
          2 * pi *
          (fs (Z :. sf) * fromIntegral s * deltaScale / (log maxScale) +
           fx (Z :. xf) * fromIntegral ((x :: Int) - centerX) * deltaX / xLen +
           fy (Z :. yf) * fromIntegral ((y :: Int) - centerY) * deltaY / yLen)))
  return .
    fromListUnboxed (Z :. numX :. numY :. numOrientations :. numScales) .
    NL.toList $
    transformMat NL.#> vec

{-# INLINE stsTor2 #-}
stsTor2 ::
     Int
  -> Double
  -> [Double]
  -> Int
  -> Double
  -> [Double]
  -> R.Array U DIM4 (Complex Double)
  -> IO (R.Array U DIM2 (Complex Double))
stsTor2 numX xLen xFreqs numY yLen yFreqs arr = do
  let deltaX = xLen / fromIntegral numX
      deltaY = yLen / fromIntegral numY
      centerX = div numX 2
      centerY = div numY 2
  vec <-
    fmap (NL.fromList . R.toList) .
    sumP . sumS . R.map (\x -> (magnitude x) ^ 2 :+ 0) $
    arr
  transformMat <-
    fmap
      (((numX * numY) >< (L.product . L.drop 2 . listOfShape . extent $ arr)) .
       R.toList) .
    computeUnboxedP .
    R.traverse2
      (fromListUnboxed (Z :. L.length xFreqs) xFreqs)
      (fromListUnboxed (Z :. L.length yFreqs) yFreqs)
      (\(Z :. numXFreq) (Z :. numYFreq) ->
         (Z :. numX :. numY :. numXFreq :. numYFreq)) $ \fx fy (Z :. x :. y :. xf :. yf) ->
      exp
        (0 :+
         (2 * pi *
          (fx (Z :. xf) * fromIntegral (x - centerX) * deltaX / xLen +
           fy (Z :. yf) * fromIntegral (y - centerY) * deltaY / yLen)))
  return . fromListUnboxed (Z :. numX :. numY) . NL.toList $
    transformMat NL.#> vec

{-# INLINE stsTor2' #-}
stsTor2' ::
     Int
  -> Double
  -> [Double]
  -> Int
  -> Double
  -> [Double]
  -> R.Array U DIM4 (Complex Double)
  -> IO (R.Array U DIM2 Double)
stsTor2' numX xLen xFreqs numY yLen yFreqs arr = do
  let deltaX = xLen / fromIntegral numX
      deltaY = yLen / fromIntegral numY
      centerX = div numX 2
      centerY = div numY 2
      (Z :. numXF :. numYF :. numTF :. numSF) = extent arr
      mat = ((numXF * numYF) >< (numTF * numSF)) (R.toList arr)
  transformMat <-
    fmap
      (((numX * numY) >< (L.product . L.drop 2 . listOfShape . extent $ arr)) .
       R.toList) .
    computeUnboxedP .
    R.traverse2
      (fromListUnboxed (Z :. L.length xFreqs) xFreqs)
      (fromListUnboxed (Z :. L.length yFreqs) yFreqs)
      (\(Z :. numXFreq) (Z :. numYFreq) ->
         (Z :. numX :. numY :. numXFreq :. numYFreq)) $ \fx fy (Z :. x :. y :. xf :. yf) ->
      exp
        (0 :+
         (2 * pi *
          (fx (Z :. xf) * fromIntegral (x - centerX) * deltaX / xLen +
           fy (Z :. yf) * fromIntegral (y - centerY) * deltaY / yLen)))
  return .
    sumS .
    sumS .
    R.map (\x -> (magnitude x) ^ 2) .
    fromListUnboxed (Z :. numX :. numY :. numTF :. numSF) . NL.toList . flatten $
    transformMat NL.<> mat

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE TypeOperators #-}
---------------------------------
-- |
-- Module      :  HMM
-- Copyright   :  (C) Kaushik Chakraborty, 2019
-- License     :  Apache v2 (see the file LICENSE)
-- Maintainer  :  Kaushik Chakraborty <git@kaushikc.org>
-- Stability   :  experimental
--
---------------------------------

module HMM
  (
    HMM(..),
    -- ** Smart Constructors
    mkHMM_, mkHMM, mkHMM1, sensorDiagonal,
    -- * Utility Functions
    normalise, inverse, reverseV, extract1, (!**!), unitColumn,
    -- ** Converting to/from HMatrix matrices
    toHM, fromHM,
    -- * Inference Algorithms
    -- ** Forward-Backward
    forward, backward, forwardBackward,
    -- ** Fixed Lag Smoothing
    fixedLagSmoothing,
    -- * Sample HMM Models
    umbrellaHMM,
    -- * Miscellaneous
    -- ** Type Synonyms
    R, TransitionModel, SensorDiagonal, SensorModel, Message, Distribution,
    -- ** Data Types
    Persistent(..),
    -- ** Test Functions
    runFLSAlgo, runFBAlgo,
    -- ** Lenses
    prior, tModel, sDist, sModel
) where

import Text.Printf (printf)

import Papa
import GHC.TypeNats
import Data.Proxy (Proxy(..))

import Control.Monad.State (State, runState)
import Control.Monad.Trans.Maybe (MaybeT, runMaybeT)

import qualified Data.Vector as Vec

import Linear.V
import Linear

import qualified Numeric.LinearAlgebra.Static as HM
import qualified Numeric.LinearAlgebra as HMat
-----------------------------------------------------------------------------------

type R = Double
type TransitionModel (s :: Nat) a = V s (V s a)
type SensorDiagonal (s  :: Nat) a = V s (V s a)
type SensorModel (t :: Nat) (s :: Nat) a = V t (SensorDiagonal s a)
type Message (s :: Nat) a = V s (V 1 a)
type Distribution (s :: Nat) a = V s a

-- | Hidden Markov Models
data HMM (s :: Nat) (t :: Nat)  a b = HMM {
  -- | prior distribution which is used as initial forward message
  _prior :: Message s a,
  -- | transition model with @s@ states as @sxs@ matrix
  _tModel :: TransitionModel s a,
  -- | evidence value vector, each index maps to the corresponding evidence values
  _sDist :: V t b,
  -- | sensor model with @t@ evidence values having @sxs@ diagonal matrix capturing their ditributions for each state @s@
  _sModel :: SensorModel t s a
  } deriving (Show)

makeLenses ''HMM

mkHMM_ :: (KnownNat s, KnownNat t) => Maybe (Message s R) -- ^ Prior distribution on the initial state, @P(X0)@. If nothing then considered @(0.5,0.5)@
       -> TransitionModel s R -- ^ Transition Model as @sxs@ matrix, @s@ being the number of states
       -> V t b -- ^ evidence values vector map
       -> V t (V s R) -- ^ a @txs@ matrix for each evidence value @t@, a vector capturing conditional probabilities for each state @s@
       -> HMM s t R b
mkHMM_ mp xs evs = let p = fromMaybe (V $ Vec.replicate (dim xs) 0.5) mp
              in
              HMM p xs evs . (scaled <$>)

-- | A default HMM where both transition model states and evidence variables have boolean support
mkHMM :: KnownNat s => TransitionModel s R -> V 2 (V s R) -> HMM s 2 R Bool
mkHMM ts = mkHMM_ Nothing ts (V $ Vec.fromList [True, False])

-- | A default HMM with one state variable having boolean support
mkHMM1 :: V2 (V2 R) -> V2 (V2 R) -> HMM 2 2 R Bool
mkHMM1 ts ss = let f = (_V #) . over mapped toV
               in
                 uncurry mkHMM $ over both f (ts , ss)

-- | create a @sxs@ diagonal matrix with corresponding posterior probabilities from the sensor model of the @HMM@ for an input sensor value
sensorDiagonal :: (KnownNat s, KnownNat t, Eq b)
               => HMM s t a b
               -> b -- ^ sensor value
               -> Maybe (SensorDiagonal s a)
sensorDiagonal hmm e = findIndexOf (sDist.folded) (== e) hmm >>= (\i -> hmm ^? sModel . ix i)

-- | Multiply each number by a constant such that the sum is 1.0
normalise :: (Fractional a, Foldable f, Functor f) => f a -> f a
normalise xs | null xs = xs
             | otherwise =
                 let
                   s = sum xs
                 in
                   (/ s) <$> xs

-- | inverse a square matrix
inverse :: forall s. (KnownNat s) => V s (V s R) -> V s (V s R)
inverse  = fromHM . HM.inv . toHM

-- | reverse contents of a 'Linear.V' vector
reverseV :: forall s a. (KnownNat s) => V s a -> V s a
reverseV = V . Vec.fromList . Papa.reverse . toList

extract1 :: (KnownNat s) => V s (V 1 a) -> V s a
extract1 =  V . foldMap toVector

-- | Pointwise product of 2 @mxn@ 'V' vectors
infix 8 !**!
(!**!) :: (KnownNat m, KnownNat n) => V m (V n R) -> V m (V n R) -> V m (V n R)
as !**! bs = fromHM $ toHM as * toHM bs

-- | an unit vector having @s@ columns
unitColumn :: forall s. KnownNat s => Message s R
unitColumn = V $ Vec.replicate (fromIntegral $ natVal (Proxy :: Proxy s))  (toV $ V1 1.0)

-- | take a @sxt@ 'Linear.V.V' matrix and give corresponding 'Numeric.LinearAlgebra.Static.L' version
toHM :: (KnownNat s, KnownNat t) => V s (V t R) -> HM.L s t
toHM = HM.matrix . foldMap (Vec.toList . toVector)

-- | take a @sxt@ 'Numeric.LinearAlgebra.Static.L' matrix and give corresponding 'Linear.V.V' version
fromHM :: (KnownNat s, KnownNat t) => HM.L s t -> V s (V t R)
fromHM m = V $ V <$> Vec.fromList (fmap Vec.fromList $ HMat.toLists $ HM.extract m)

-- | Filtering message propagated forward
forward :: (KnownNat s) => HMM s t R b -> SensorDiagonal s R  -> Message s R -> Message s R
forward hmm ot f = normalise (ot !*! Linear.transpose (hmm ^. tModel) !*! f)

-- | Smoothing message propagated backward
backward :: (KnownNat s) => HMM s t R b -> SensorDiagonal s R -> Message s R -> Message s R
backward hmm ok1 bk2 = (hmm ^. tModel) !*! ok1 !*! bk2

-- | The forward–backward algorithm for smoothing: computing posterior prob- abilities of a sequence of states given a sequence of observations
forwardBackward :: forall s t u b . (KnownNat t, KnownNat s, KnownNat u, Eq b)
                => HMM s u R b -- ^ HMM model as a way to implement
                -> V t b -- ^ list of evidences for each time step
                -> V t (Distribution s R)
forwardBackward hmm evs = let
  -- reifying the number of evidences
  tNat = dim evs
  -- forward messages from time t .. 0
  fv :: V (t + 1) (Message s R)
  fv = V $ Vec.fromList $ foldl' (\m@(x:_) e -> forward hmm (sensorDiagonal hmm e ^?! _Just) x : m) [hmm ^. prior] evs
  -- getting rid of the prior message from the end of the list
  -- so now forward messages vector is from time t .. 1
  fv_0 :: V t (Message s R)
  fv_0 = V $ Vec.fromList $ fv ^.. taking tNat traversed

  --  backward messages from time 1 .. t
  bs :: V t (Message s R)
  bs = V $ Vec.fromList $ foldl' (\m@(x:_) e -> backward hmm (sensorDiagonal hmm e ^?! _Just) x : m) [unitColumn] $ reverseV evs
  -- reversing the backward messages
  -- so now the backward messages vector is from time t .. 1
  revBs :: V t (Message s R)
  revBs = reverseV bs
  in
    -- smoothing probabilities in reverse order of the list of evidences
    -- i.e. starting from time t .. 1
    liftA2 (\f b' -> extract1 $ normalise $ f !**! b') fv_0 revBs

-- | Persistent State
data Persistent (s :: Nat) a b d= Persistent {
  _t :: d, -- ^ current time
  _f_msg :: Message s a, -- ^ the forward message @P(Xt|e1:t)@
  _b :: V s (V s a), -- ^ the @d@-step backward transformation matrix
  _e_td_t :: Vec.Vector b -- ^ double-ended list of evidence from @t − d@ to @t@
  } deriving (Show)

makeLenses ''Persistent

-- | Initial persistent state where
--
--  * t = 1
--
--  * f_msg  = hmm.prior
--
--  * b = identity matrix
--
--  * e_td_t = empty vector
persistentInit :: forall s b d. (KnownNat s, Integral d) => Message s R -> Persistent s R b d
persistentInit p = Persistent { _t = 1, _f_msg = p, _b = identity, _e_td_t = Vec.empty}

-- | Online Algorithm for smoothing with a fixed time lag of @d@ steps
fixedLagSmoothing :: forall s u b d. (KnownNat s, KnownNat u, Eq b, Integral d)
                  => HMM s u R b -- ^ HMM model
                  -> d -- ^ length of lag
                  -> b -- ^ evidence at time @t@
                  -> MaybeT (State (Persistent s R b d)) (Distribution s R)
fixedLagSmoothing hmm d e = do
  e_td_t %= flip Vec.snoc e

  o_t <- uses e_td_t $ (^?! _Just) . sensorDiagonal hmm . Vec.last

  t' <- use t
  if t' > d then
    do
      e_td_t %= Vec.drop 1
      o_tmd <- uses e_td_t $ (^?! _Just) . sensorDiagonal hmm . Vec.head
      f_msg %= forward hmm o_tmd
      b %= \b' -> inverse o_tmd !*! inverse (hmm ^. tModel) !*! b' !*! (hmm ^. tModel) !*! o_t
      t += 1

      (f'', b'') <- liftA2 (,) (use f_msg) (use b)

      return $ extract1 $ normalise (f'' !**! (b'' !*! unitColumn))
  else
    do
      b %= \b' -> b' !*! (hmm ^. tModel) !*! o_t
      t += 1
      mzero

umbrellaHMM :: HMM 2 2 R Bool
umbrellaHMM = mkHMM1 (V2 (V2 0.7 0.3) (V2 0.3 0.7)) (V2 (V2 0.9 0.2) (V2 0.1 0.8))

runFLSAlgo  :: [Bool] -> Integer -> [Maybe (Distribution 2 String)]
runFLSAlgo bs d = let hmm = mkHMM1 (V2 (V2 0.7 0.3) (V2 0.3 0.7)) (V2 (V2 0.9 0.2) (V2 0.1 0.8))
                      initState = persistentInit (V $ Vec.fromList [V $ Vec.singleton 0.5, V $ Vec.singleton 0.5]) :: Persistent 2 R Bool Integer
                      algo = fixedLagSmoothing hmm d
                      a :: ([Maybe (Distribution 2 R)], Persistent 2 R Bool Integer)
                      a = foldl' (\(rs, s) x -> runState (runMaybeT $ algo x) s & _1 #%~ ((rs ++) . pure)) ([], initState) bs
                  in
                    (a ^. _1) & traverse.traverse.traverse %~ \(x :: R) -> printf "%.3f" x :: String

runFBAlgo :: [Bool] -> [Distribution 2 R]
runFBAlgo bs = reifyVectorNat (Vec.fromList bs) (toList . forwardBackward umbrellaHMM)

max_stars_repo_name: kaychaks/hmm

max_issues_repo_name: kaychaks/hmm

max_forks_repo_name: kaychaks/hmm


text:
module FFTSpec where

import qualified Clash.Prelude as Clash
import Clash.Prelude (Signal, Vec(..), BitVector, Index, Signed, Unsigned, SFixed, Bit, SNat(..),
                      simulate, simulate_lazy, listToVecTH, KnownNat, pack, unpack, (++#), mealy, mux, bundle, unbundle,
                      HiddenClockResetEnable)
import Test.Hspec
import Test.QuickCheck

import qualified Numeric.FFT as FFT
import qualified Data.Complex as C
import Clash.DSP.Complex (Complex, fromComplex, toComplex)
import Clash.DSP.FFT.Parallel

--FFT
spec = describe "Parallel FFTs" $ do
    specify "Recursive decimation in time equals known good implementation"      $ property prop_fftDITRec
    specify "Recursive decimation in frequency equals known good implementation" $ property prop_fftDIFRec
    specify "Iterative decimation in time equals known good implementation"      $ property prop_fftDITIter
    specify "Iterative decimtaion in frequency equals known good implementation" $ property prop_fftDIFIter

approxEqual :: Double -> Double -> Bool
approxEqual x y = abs (x - y) < 0.0001

approxEqualComplex (a C.:+ b) (c C.:+ d) = approxEqual a c && approxEqual b d

twiddles :: Vec 8 (Complex Double)
twiddles = $(listToVecTH (twiddleFactors 8))

prop_fftDITRec :: Vec 16 (C.Complex Double) -> Bool
prop_fftDITRec vec = and $ zipWith approxEqualComplex (map toComplex (Clash.toList (fftDITRec twiddles (Clash.map fromComplex vec)))) (FFT.fft (Clash.toList vec))

prop_fftDIFRec :: Vec 16 (C.Complex Double) -> Bool
prop_fftDIFRec vec = and $ zipWith approxEqualComplex (map toComplex (Clash.toList (fftDIFRec twiddles (Clash.map fromComplex vec)))) (FFT.fft (Clash.toList vec))

prop_fftDITIter :: Vec 16 (C.Complex Double) -> Bool
prop_fftDITIter vec = and $ zipWith approxEqualComplex (map toComplex (Clash.toList (fftDITIter twiddles (Clash.map fromComplex vec)))) (FFT.fft (Clash.toList vec))

prop_fftDIFIter :: Vec 16 (C.Complex Double) -> Bool
prop_fftDIFIter vec = and $ zipWith approxEqualComplex (map toComplex (Clash.toList (fftDIFIter twiddles (Clash.map fromComplex vec)))) (FFT.fft (Clash.toList vec))

max_stars_repo_name: adamwalker/clash-utils

max_issues_repo_name: adamwalker/clash-utils

max_forks_repo_name: adamwalker/clash-utils


text:
module Quadratic where

import Data.Complex

data Quadratic = Quadratic { a :: Double, b :: Double, c :: Double } deriving Show

type RootT = Complex Double

data Roots = Roots RootT RootT deriving Show

-- | Calculates roots of a polynomial and return set of roots
roots :: Quadratic -> Roots

-- Trivial, all constants are zero, error roots are not defined
roots (Quadratic 0 0 _) = error "Not a quadratic polynomial"

-- Is a polynomial of degree 1, x = -c / b
roots (Quadratic 0.0 b c) = let root = ( (-c) / b :+ 0)
                            in Roots root root

-- b^2 - 4ac = 0
roots (Quadratic a b c) =
  let discriminant = b * b - 4 * a * c
  in rootsInternal (Quadratic a b c) discriminant

rootsInternal :: Quadratic -> Double -> Roots
-- Discriminant is zero, roots are real
rootsInternal q d | d == 0 = let r = (-(b q) / 2.0 / (a q))
                                 root = r :+ 0
                             in Roots root root

-- Discriminant is negative, roots are complex
rootsInternal q d | d < 0 = Roots (realpart :+ complexpart) (realpart :+ (-complexpart))
  where plusd = -d
        twoa = 2.0 * (a q)
        complexpart = (sqrt plusd) / twoa
        realpart = - (b q) / twoa

-- discriminant is positive, all roots are real
rootsInternal q d = Roots (root1 :+ 0) (root2 :+ 0)
  where plusd = -d
        twoa = 2.0 * (a q)
        dpart = (sqrt plusd) / twoa
        prefix = - (b q) / twoa
        root1 = prefix + dpart
        root2 = prefix - dpart

max_stars_repo_name: subhashischatterjee/HuskellCookbook

max_issues_repo_name: subhashischatterjee/HuskellCookbook

max_forks_repo_name: subhashischatterjee/HuskellCookbook


text:
{-# LANGUAGE DataKinds, ExistentialQuantification, FlexibleContexts         #-}
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, TypeSynonymInstances #-}
{-# LANGUAGE UndecidableInstances                                           #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
module Numeric.Algebra.Ring.Noetherian (Noetherian) where
import qualified Data.Complex              as C
import           Data.Function
import           Data.Ord
import           Data.Ratio
import           Numeric.Algebra
import qualified Numeric.Algebra.Complex   as NA
import           Numeric.Algebra.Instances
import           Prelude                   hiding (negate, subtract, (*), (+),
                                            (-))
import qualified Prelude                   as P
-- | Noetherian ring (i.e. every ideal is finitely generated).
class (Commutative r, Ring r) => Noetherian r where

instance Noetherian Int where

instance Noetherian Integer where

instance (Commutative (NA.Complex r), Ring (NA.Complex r)) => Noetherian (NA.Complex r) where
instance (Commutative (C.Complex r), Ring (C.Complex r)) => Noetherian (C.Complex r) where
instance Integral n => Noetherian (Ratio n)

max_stars_repo_name: konn/algebra-extras

max_issues_repo_name: konn/algebra-extras

max_forks_repo_name: konn/algebra-extras


text:
{-# LANGUAGE DataKinds #-}

import Math.Tensor
import Math.Tensor.LorentzGenerator
import Math.Tensor.Examples.Gravity
import Math.Tensor.Examples.Gravity.DiffeoSymEqns
import Math.Tensor.Examples.Gravity.Schwarzschild

import Math.Tensor.Internal.LinearAlgebra

import Numeric.LinearAlgebra (rank)
import Numeric.LinearAlgebra.Data (toLists, fromLists, size, cmap, (===), Matrix, Z)
import Data.List (nub, nubBy, findIndex, sort)

import Data.Maybe (mapMaybe)

import Data.Int (Int64)
import Data.Ratio

import Control.Parallel.Strategies (parList, runEvalIO, rdeepseq)

import qualified Data.IntMap.Strict as I

main = do

  let ans0 = fromListT6 [((Empty, Empty, Empty, Empty, Empty, Empty), AnsVar $ I.singleton 1 (SField 1))] :: ATens 0 0 0 0 0 0 AnsVarR
--  let (eta4,eps4,ans4) = mkAnsatzTensorFastAbs 4 symList4 areaList4 :: (AnsatzForestEta, AnsatzForestEpsilon, ATens 1 0 0 0 0 0 AnsVarR)
  let ans4 = ZeroTensor :: ATens 1 0 0 0 0 0 AnsVarR
  let (eta6,eps6,ans6) = mkAnsatzTensorFastAbs 6 symList6 areaList6 :: (AnsatzForestEta, AnsatzForestEpsilon, ATens 1 0 1 0 0 0 AnsVarR)
  let (eta8,eps8,ans8) = mkAnsatzTensorFastAbs 8 symList8 areaList8 :: (AnsatzForestEta, AnsatzForestEpsilon, ATens 2 0 0 0 0 0 AnsVarR)
  let (eta10_1,eps10_1,ans10_1) = mkAnsatzTensorFastAbs 10 symList10_1 areaList10_1 :: (AnsatzForestEta, AnsatzForestEpsilon, ATens 2 0 0 0 2 0 AnsVarR)
  let (eta10_2,eps10_2,ans10_2) = mkAnsatzTensorFastAbs 10 symList10_2 areaList10_2 :: (AnsatzForestEta, AnsatzForestEpsilon, ATens 2 0 1 0 0 0 AnsVarR)

  let r0    = tensorRank6' ans0
  let r4    = tensorRank6' ans4
  let r6    = tensorRank6' ans6
  let r8    = tensorRank6' ans8
  let r10_1 = tensorRank6' ans10_1
  let r10_2 = tensorRank6' ans10_2

  let ans0'    = ans0
  let ans4'    = shiftLabels6 r0 ans4
  let ans6'    = shiftLabels6 (r0 + r4) ans6
  let ans8'    = shiftLabels6 (r0 + r4 + r6) ans8
  let ans10_1' = shiftLabels6 (r0 + r4 + r6 + r8) ans10_1
  let ans10_2' = shiftLabels6 (r0 + r4 + r6 + r8 + r10_1) ans10_2

  let ansaetze = ans0' &.&>
                 ans4' &.&>
                 ans6' &.&>
                 ans8' &.&>
                 ans10_1' &.&>
                 (singletonTList6 ans10_2')

  putStrLn $ "dimension of ansatz space   : " ++ show (tensorRank6 ansaetze)

  let two = SField (2 :: Rational)

  let e1 = eqn1 ans0' ans4'
  let e2 = eqn3 ans6'
  let e3 = eqn1A ans4' (two &. ans8')
  let e4 = eqn1AI ans6' ans10_2'
  let e5 = eqn2Aa ans6' (two &. ans10_1')
  let e6 = eqn3A ans6' ans10_2'

  let system = e1 &.&>
               e2 &.&>
               e3 &.&>
               e4 &.&>
               e5 &.&>
               (singletonTList6 e6)

  putStrLn $ "rank of system              : " ++ show (tensorRank6 system)

  let solution = redefineVarsSystem6 $ solveSystem6 system ansaetze

  case solution of
    (t1 `AppendTList6` (
     t2 `AppendTList6` (
     t3 `AppendTList6` (
     t4 `AppendTList6` (
     t5 `AppendTList6` (
     t6 `AppendTList6` (
     EmptyTList6)))))))
      -> do
          let Just ans0'' = tryAsATens ans0 t1
          let Just ans4'' = tryAsATens ans4 t2
          let Just ans6'' = tryAsATens ans6 t3
          let Just ans8'' = tryAsATens ans8 t4
          let Just ans10_1'' = tryAsATens ans10_1 t5
          let Just ans10_2'' = tryAsATens ans10_2 t6

          let e1' = removeZeros6 $ eqn1 ans0'' ans4''
          let e2' = removeZeros6 $ eqn3 ans6''
          let e3' = removeZeros6 $ eqn1A ans4'' (two &. ans8'')
          let e4' = removeZeros6 $ eqn1AI ans6'' ans10_2''
          let e5' = removeZeros6 $ eqn2Aa ans6'' (two &. ans10_1'')
          let e6' = removeZeros6 $ eqn3A ans6'' ans10_2''

          putStrLn $ "dimension of solution space : " ++ show (tensorRank6 solution)
          putStrLn ""
          putStrLn "Equations on solution space :"

          print e1'
          print e2'
          print e3'
          print e4'
          print e5'
          print e6'

          let kin' = ans10_2'' &- (contrATens3 (0,0) $ contrATens3 (1,1) $ interI2 &* ans10_1'')
          let kin = kin' &+ (tensorTrans1 (0,1) kin')

          putStrLn ""
          putStrLn $ "number of parameters in kinetic part of EOM : " ++ show (tensorRank6' kin)
          putStrLn $ "number of parameters in mass part of EOM    : " ++ show (tensorRank6' ans8'')

          let n1 = removeZeros6 $ (delta3A &* ans4'') &+ (contrATens1 (1,0) $ interArea &* ans4'') &+ (contrATens1 (0,0) $ contrATens1 (1,1) $ two &. flatArea &* interArea &* ans8'')

          let n2_1 = contrATens1 (0,0) $ contrATens1 (1,1) $ contrATens2 (0,0) $ flatArea &* interArea &* interJ2 &* ans10_2''
          let n2_2 = contrATens1 (0,0) $ contrATens1 (2,1) $ contrATens2 (0,0) $ flatArea &* interArea &* interJ2 &* ans10_2''
          let n2_3 = contrATens1 (0,0) $ contrATens1 (1,1) $ flatArea &* interArea &* ans10_1''

          let n2 = removeZeros6 $ cyclicSymATens5 [0,1,2] $ n2_1 &+ n2_2 &- (two &. n2_3)

          putStrLn ""
          putStrLn "Noether identity 1 :"
          print n1

          putStrLn ""
          putStrLn "Noether identity 2 :"
          print n2

  return ()


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleContexts #-}
-----------------------------------------------------------------------------
--
-- Module      :  AI.Network.RNN.Genetic
-- Copyright   :  (c) JP Moresmau
-- License     :  BSD3
--
-- Maintainer  :  JP Moresmau <jp@moresmau.fr>
-- Stability   :  experimental
-- Portability :
--
-- | Code for genetic evolution of a network
--
-----------------------------------------------------------------------------

module AI.Network.RNN.Genetic where

import Control.Monad
import Control.Monad.Random hiding (fromList)
import qualified Control.Monad.Random as MR
import Control.DeepSeq
import Data.IORef
import AI.GeneticAlgorithm.Simple
import AI.Network.RNN.Types
import AI.Network.RNN.Util
import Numeric.LinearAlgebra.HMatrix

import Debug.Trace

-- | Keep all relevant data together
data RNNData a b = (RNNEval a) => RNNData a !(TrainData b) !Double
   -- deriving (Show,Read,Eq)

-- | evaluate network
instance NFData (RNNData a b) where
    rnf (RNNData rnn _ _) = rnf rnn

-- | Error threshold (not currently used)
errorThreshold :: Double
errorThreshold = 0.1

-- | Chromosome instance
instance Chromosome (RNNData a b) where
    -- | cross network by taking half the values of one and half of the other
    crossover (RNNData rnn1 td l) (RNNData rnn2 _ _) = do
        rnns <- crossNetworkFull rnn1 rnn2
        --idx <- getRandomR (0, size rnn1 - 1)
        --let (v3,v4) = crossVector rnn1 rnn2 idx
        return $ map (\r->RNNData r td l) rnns -- [v3,v4]

    mutation (RNNData rnn1 td l) = do
        rnn2 <- mutateNetwork rnn1
        --rnn2 <- stdNormalMutate rnn1
        return $ RNNData rnn2 td l
    -- | fitness is "higher is better", inverse of cost
    fitness (RNNData rnn1 td _) = 1 / cost rnn1 td -- (fromVector sz rnn1)

-- | Cross 2 networks by taking the first half of one, the second half of the other (not currently used)
crossNetworkFull :: (RandomGen g,RNNEval a) =>  a -> a -> Rand g [a]
crossNetworkFull rnn1 rnn2 = do
    let v1 = toVector rnn1
        v2 = toVector rnn2
    idx <- getRandomR (0, size v1 - 1)
    let (v3,v4) = crossVector v1 v2 idx
        rnn3 = fromVector (rnnsize rnn1) v3
        rnn4 = fromVector (rnnsize rnn1) v4
    return $[rnn3,rnn4]

-- | cross 2 networks by taking half the values of one and half of the other
crossNetworkHalf :: (RandomGen g,RNNEval a) =>  a -> a -> Rand g [a]
crossNetworkHalf rnn1 rnn2 = do
    let v1 = toVector rnn1
        v2 = toVector rnn2
    (v3,v4) <- mixVector v1 v2 0.5
    let
        rnn3 = fromVector (rnnsize rnn1) v3
        rnn4 = fromVector (rnnsize rnn1) v4
    return $ force [rnn3,rnn4]

crossMutate :: (RandomGen g,RNNEval a) =>  a -> a -> Rand g [a]
crossMutate rnn1 rnn2 = do
    let v1 = toVector rnn1
        v2 = toVector rnn2
    v3 <- stdNormalMutate v1
    v4 <- stdNormalMutate v2
    let
        rnn3 = fromVector (rnnsize rnn1) v3
        rnn4 = fromVector (rnnsize rnn1) v4
    return [rnn3,rnn4]

-- | Cross networks by taking the average of their values (not currently used)
avgNetwork :: (RandomGen g,RNNEval a) =>  a -> a -> Rand g [a]
avgNetwork rnn1 rnn2 = do
    let v1 = toVector rnn1
        v2 = toVector rnn2
        v3 = fromList $ zipWith (\x y -> (x+y)/2) (toList v1) (toList v2)
        rnn3 = fromVector (rnnsize rnn1) v3
      --  rnn4 = fromVector (rnnsize rnn1) v4 of
      --      Right n -> n
      --      Left e  -> error e
    return [rnn3]

crossMatrixEq :: Matrix Double -> Matrix Double -> (Matrix Double,Matrix Double)
crossMatrixEq m1 m2 =
    let row1      = upperHalf $ rows m1
        col1      = upperHalf $ cols m1
        [[a1,a2],[a3,a4]] = toBlocksEvery row1 col1 m1
        [[b1,b2],[b3,b4]] = toBlocksEvery row1 col1 m2
    in (fromBlocks [[a1,b2],[a3,b4]],fromBlocks [[b1,a2],[b3,a4]])
    where
        upperHalf a = (if odd a then a+1 else a) `div` 2

crossVectorEq :: Vector Double -> Vector Double -> (Vector Double,Vector Double)
crossVectorEq v1 v2 = crossVector v1 v2 (size v1 `div` 2)

-- | Cross vector at the given point
crossVector :: Vector Double -> Vector Double -> Int -> (Vector Double,Vector Double)
crossVector v1 v2 idx =
--    let m1 = asRow v1
--        m2 = asRow v2
--        [[a1,a2]] = toBlocks [1] [idx,size v1-idx]  m1
--        [[b1,b2]] = toBlocks [1] [idx,size v1-idx] m2
--    in (head $ toRows $ fromBlocks [[a1,b2]],head $ toRows $ fromBlocks [[b1,a2]])
--    let sz2 = size v2 - idx - 1
--        [a1,a2] = takesV [idx,sz2] v1
--        [b1,b2] = takesV [idx,sz2] v2
--    in (vjoin [a1,b2],vjoin [a2,b1])
 let sz = size v1-idx
 in force $ (vjoin [subVector 0 idx v1,subVector idx sz v2],
                    vjoin [subVector 0 idx v2,subVector idx sz v1])

-- | Mix vector values given a probability to take a value from the first one and not the second
mixVector :: (Monad m,RandomGen g) =>  Vector Double -> Vector Double -> Double -> RandT g m (Vector Double,Vector Double)
mixVector v1 v2 prob = do
    -- rs <- replicateM (size v1) (getRandomR (0, 1))
    let (l1,l2) = swapL (toList v1) (toList v2) True
        -- unzip $ map swapR $ zip3 (toList v1) (toList v2) [0..]
    return (fromList l1,fromList l2)
    where
        swapR (a,b,idx) = if (idx `mod` 2) == 0 then (a,b) else (b,a)
        --if idx <= prob then (a,b) else (b,a):
        swapL :: [a] -> [a] -> Bool -> ([a],[a])
        swapL [] _ _ = ([],[])
        swapL (x:xs) (y:ys) True = let
            (rs1,rs2)=swapL xs ys False in (x:rs1,y:rs2)
        swapL (x:xs) (y:ys) False = let
            (rs1,rs2)=swapL xs ys True in (y:rs1,x:rs2)

-- | Mix vector values given a probability to take a value from the first one and not the second
mixList :: (Monad m,RandomGen g) =>  [Double] -> [Double] -> Double -> RandT g m ([Double],[Double])
mixList v1 v2 prob = do
    --rs <- getRandomRs (0,1)
    return $ unzip $ map swapR $ zip3 v1 v2 [0..]
    where
        swapR (a,b,idx) = if (idx `mod` 2) == 0 then (a,b) else (b,a)

stdNormalMutate :: (Monad m,RandomGen g) =>  Vector Double -> RandT g m (Vector Double)
stdNormalMutate v1 = do
    --rs <- replicateM (size v1) stdNormal
    -- let mx = maxElement v1
    -- let mn = minElement v1
    --let sz= trace ("v1:"++ show (mn,mx)) $ size v1
    seed <- getRandom
    let rs = randomVector seed Gaussian $ size v1
    let v2 = (rs) + v1
--    let nv2 = cmap (\x->
--                            if x < 0
--                                then 0
--                                else if x > 1
--                                    then 1
--                                    else x)
--         v2
    -- let mx2 = maxElement v2
    -- let mn2 = minElement v2
    return v2
    -- return $ trace ("v2:"++show (mn2,mx2)) v2

-- | network mutation: select randomly from several mutation algorithms
mutateNetwork :: (RandomGen g,RNNEval a) =>  a -> Rand g a
mutateNetwork rnn = do
    f <- MR.fromList [(stdNormalMutation,1),(pointMutation,1),(swapMutation,1),(insertMutation,1),(flipMutation,0.5)]
    f rnn

-- | Mutate values around a standard deviation
stdNormalMutation :: (RandomGen g,RNNEval a) =>  a -> Rand g a
stdNormalMutation rnn = do
    let v1 = toVector rnn
    v2 <- stdNormalMutate v1
    return $ force $ fromVector (rnnsize rnn) v2

-- | Mutate a single value of the vector, taking a random value
pointMutation :: (RandomGen g,RNNEval a) =>  a -> Rand g a
pointMutation rnn = do
    let v1 = toVector rnn
    idx <- getRandomR (0, size v1 - 1)
    dbl <- getRandomR (0,1)
    let v2 = accum v1 const [(idx,dbl)]
    return $ fromVector (rnnsize rnn) v2

-- | Swap two values in the vector
swapMutation :: (RandomGen g,RNNEval a) =>  a -> Rand g a
swapMutation rnn = do
    let v1 = toVector rnn
    idx1 <- getRandomR (0, size v1 - 1)
    idx2 <- getRandomR (0, size v1 - 1)
    let v2 = accum v1 const [(idx1,atIndex v1 idx2),(idx2,atIndex v1 idx1)]
    return $ fromVector (rnnsize rnn) v2

-- | Insert a random value at a random point in the vector, discarding the last value
insertMutation :: (RandomGen g,RNNEval a) =>  a -> Rand g a
insertMutation rnn = do
    let v1 = toVector rnn
    idx <- getRandomR (0, size v1 - 1)
    dbl <- getRandomR (0,1)
    let v2 = vjoin [subVector 0 idx v1,fromList [dbl],subVector idx (size v1 - idx -1) v1]
    return $ fromVector (rnnsize rnn) v2

flipMutation :: (RandomGen g,RNNEval a) =>  a -> Rand g a
flipMutation rnn = do
    let v1 = toVector rnn
    let v2 = cmap (\a->(-a)) v1
    return $ fromVector (rnnsize rnn) v2

-- | Stop function
-- we take a list of past fitnesses and stop when the fitness 5 generations ago was not worse that the current one
stopf :: IORef [Double] -> Int -> (RNNData a b,Double) -> Int -> IO Bool
stopf fitnessList maxGen (nd,fit) gen = do
    mfit <- atomicModifyIORef' fitnessList (\l->
        let l2 = take 5 (fit:l)
        in if length l2 == 5
            then (l2, Just $ last l2)
            else (l2, Nothing)
        )

    print $ "Fitness (" ++ show gen ++"): " ++ show fit
    let converged = case mfit of
                        Nothing -> False
                        Just f  -> f >= fit
    when converged $ print "Converged!"
    return ( gen >= maxGen)

-- | Build a random network data
buildNetworkData :: (RNNEval re,Monad m,RandomGen g) => TrainData b -> Size re -> FullSize re -> RandT g m (RNNData re b)
buildNetworkData td@(TrainData is _ _ _) dim fullSz = do
  --n <- randomNetwork dim fullSz
  let n = startNetwork dim fullSz
  return $ RNNData n td (fromIntegral $ length is)

max_stars_repo_name: JPMoresmau/rnn

max_issues_repo_name: JPMoresmau/rnn

max_forks_repo_name: JPMoresmau/rnn


text:
{-# LANGUAGE BangPatterns, ScopedTypeVariables, CPP #-}
{-# OPTIONS_GHC -fwarn-unused-imports #-}

-- The goal of this benchmark is to demonstrate kernels that have private
-- state which has to be read from memory on each invocation.  Further,
-- the state forms the majority of the working set -- it's larger than
-- the streaming input.

-- In theory this means that keeping a kernel on one core and
-- exploiting pipeline parallelism is better than following the data
-- through the stream graph in a depth first traversal.

--module Main(main) where
module Main where

import Control.Monad as C
import Control.Monad.Par
import Control.Monad.Par.Stream as S
import Control.Monad.Par.OpenList
import Control.DeepSeq
import Control.Exception
import Control.Parallel.Strategies as Strat

-- import Data.Array.Unboxed as U
import Data.Complex
import Data.Int
import Data.Word
import Data.List  (intersperse)
import Data.List.Split (chunk)

import Prelude as P
import System.Environment
import System.Exit
import System.CPUTime
import System.CPUTime.Rdtsc
import GHC.Conc as Conc
import GHC.IO (unsafePerformIO, unsafeDupablePerformIO, unsafeInterleaveIO)

import Debug.Trace
import Control.Monad.Par.Logging

import qualified Data.Vector.Unboxed as UV
import           Data.Vector.Unboxed hiding ((++))

-- Performs some (presently meaningless) computation on a state & a
-- window (stream element) to compute a new state and new window.
--
-- Assumes statesize is a multiple of bufsize:
statefulKern :: Vector Double -> Vector Double -> (Vector Double, Vector Double)
statefulKern state buf = (newstate, newelem)
 where
  -- We could probably test the memory behavior we're interested in
  -- better with inplace update here... but for now, this:
  newstate = UV.map (\d -> d/sum + 2) state
  newelem  = UV.map (+sum) buf

  sum = P.sum partialSums
  partialSums = [ sumslice (cutslice n) | n <- [0..factor-1] ]

  cutslice n    = UV.slice (n*blen) blen state
  sumslice slce = UV.sum (UV.zipWith (+) buf slce)

  factor       = slen `quot` blen
  slen = UV.length state
  blen = UV.length buf

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

monadpar_version (_,numfilters, bufsize, statecoef, numwins) = do
  putStrLn$ "Running monad-par version."

  let statesize = bufsize * statecoef
  results <- evaluate $ runPar$ do
       strm1 :: Stream (UV.Vector Double) <- S.generate numwins (\n -> UV.replicate bufsize 0)
       -- Make a pipeline of numfilters stages:
       let initstate = UV.generate statesize fromIntegral
       pipe_end <- C.foldM (\s _ -> streamScan statefulKern initstate s) strm1 [1..numfilters]

       sums  <- streamMap UV.sum pipe_end
#if 0
       return sums

  -- This is tricky, but two different consumers shouldn't prevent
  -- garbage collection.
  ls <- toListSpin results
--  Just (Cons h _) <- pollIVar results
  putStrLn$ "Sum of first window: " ++ show (P.head ls)
  forkIO$ measureRateList ls
  putStrLn$ "Final sum = "++ show (P.sum ls)
#else

       streamFold (+) 0 sums

  putStrLn$ "Final sum = "++ show results
#endif

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

sparks_version (_,numfilters, bufsize, statecoef, numwins) = do
  putStrLn$ "Running sparks version."

  -- Here we represent the stream as a plain list.
  let
      statesize = bufsize * statecoef
      strm1 :: [UV.Vector Double] = P.replicate numwins $ UV.replicate bufsize 0
      initstate = UV.generate statesize fromIntegral
      applyKern = scan statefulKern initstate

-- This one has the problem that it fully evaluates the stream for the
-- first kernel before moving on to the second:
--      strm_last = (parRepeatFun numfilters applyKern) strm1

      pipe_end = applyNKernels statefulKern numfilters initstate strm1

      sums = P.map UV.sum pipe_end
-- #define SERIAL
#ifndef SERIAL
	     `using` (Strat.parBuffer numCapabilities rwhnf)
#endif

  putStrLn$ "Sum of first window: "++ show (P.head sums)
  measureRateList (sums)
--  measureRateList (strm_last)
--  measureRateList (forceList strm_last)
  putStrLn$ "Final Sum = " ++ show (P.sum sums)

-- Make sure the cars of a list are evaluated before following each cdr:
forceList [] = []
forceList (h:t) = rnf h `seq` forceList t

-- A slightly different version of Data.List.scanl
scan :: (a -> b -> (a,c)) -> a -> [b] -> [c]
scan f q [] = []
scan f q (h:t) = h' : scan f q t
 where
  (q',h') = f q h

type StatefulKernel s a b = s -> a -> (s,b)

-- applyNKernels _ _ _ [] = []
-- applyNKernels :: NFData a => StatefulKernel s a a -> Int -> s -> [a] -> [a]
applyNKernels :: (NFData a, NFData s) => StatefulKernel s a a -> Int -> s -> [a] -> [a]
applyNKernels _    0 _    ls = ls
applyNKernels kern n init ls =
  applyNKernels kern (n-1) init (loop init ls)
 where
  tasklog = unsafeNewTaskSeries (nameFromValue (n,kern))

  loop st [] = []
  loop st (h:t) =
    let (st', x) =
#if 0
	           timePure tasklog$
#endif
		   kern st h in
#ifndef SERIAL
    rnf x `par`
#endif
     x : loop st' t

-- Compose two stateful kernels in parallel.
composeStatefulKernels :: (NFData b, NFData s1) =>
			  StatefulKernel s1 a b -> StatefulKernel s2 b c
		       -> StatefulKernel (s1,s2) a c
-- composeStatefulKernels (f1,f2) (s1,s2) x =
composeStatefulKernels f1 f2 (s1,s2) x =
    rnf pr1 `par` (newstate, snd pr2)
 where
  pr1 = f1 s1 x
  pr2 = f2 s2 (snd pr1)
  newstate = (fst pr1, fst pr2)

parRepeatFun n f =
--  P.foldr (.) id (P.replicate n f)
  P.foldr (.|| rdeepseq) id (P.replicate n f)

--------------------------------------------------------------------------------
-- Main script

default_version = "monad"
default_numfilters = 4
default_bufsize    = 256
default_statecoef  = 10   -- in MULTIPLES of bufsize
default_numwins    = 10 * 1000

main = do
  args <- getArgs
  arg_tup@(version,_,_,_,_) <-
       case args of
	 []          -> return (default_version, default_numfilters, default_bufsize, default_statecoef, default_numwins)
	 [a,b,c,d,e] -> return (a, read b, read c, read d, read e)
	 _         -> do
	               putStrLn$ "ERROR: Invalid arguments, must take 0 or 5 args."
		       putStrLn$ "  Expected args: (version='monad'|'sparks' #filters, bufsize, stateSizeMultiplier, #bufsToProcess)"
		       putStrLn$ "  Received args: "++ show args
		       exitFailure

  putStrLn$ "numCapabilities: "++ show numCapabilities
  putStrLn$ "  Frequency in measurable ticks:  "++ commaint oneSecond ++ "\n"

  case version of
    "monad"  -> monadpar_version arg_tup
    "sparks" -> sparks_version  arg_tup
    _        -> error$ "unknown version: "++version

  putStrLn$ "Finally, dumping all logs:"
  printAllLogs

-- It is not necessary to evaluate every element in the case of an unboxed vector.
instance NFData a => NFData (UV.Vector a) where
 rnf !vec = ()

print_ msg = trace msg $ return ()

-- work pop 1 peek N push 1
-- float->float filter
-- firFilter n coefs =
-- {

--     float sum = 0;
--     for (int i = 0; i < N; i++)
--       sum += peek(i) * COEFF[N-1-i];
--     pop();
--     push(sum);
--   }
-- }

{-

Here's what cachegrind says (on 4 core nehalem):

  $ valgrind --tool=cachegrind ./stream/disjoint_working_sets_pipeline monad 4 768 10 1000 +RTS -N4
   .....
      [measureRate] current rate: 58  Total elems&time 916  181,988,055,721
      [measureRate] Hit end of stream after 1000 elements.
     Final sum = 1.560518243231086e22
     ==21202==
     ==21202== I   refs:      7,111,462,273
     ==21202== I1  misses:          374,190
     ==21202== L2i misses:          298,364
     ==21202== I1  miss rate:          0.00%
     ==21202== L2i miss rate:          0.00%
     ==21202==
     ==21202== D   refs:      3,882,935,974  (3,542,949,529 rd   + 339,986,445 wr)
     ==21202== D1  misses:       14,606,684  (    9,824,455 rd   +   4,782,229 wr)
     ==21202== L2d misses:        6,774,479  (    2,088,565 rd   +   4,685,914 wr)
     ==21202== D1  miss rate:           0.3% (          0.2%     +         1.4%  )
     ==21202== L2d miss rate:           0.1% (          0.0%     +         1.3%  )
     ==21202==
     ==21202== L2 refs:          14,980,874  (   10,198,645 rd   +   4,782,229 wr)
     ==21202== L2 misses:         7,072,843  (    2,386,929 rd   +   4,685,914 wr)
     ==21202== L2 miss rate:            0.0% (          0.0%     +         1.3%  )

Sparks version:
     Final Sum = 1.560518243231086e22
     ==21226==
     ==21226== I   refs:      5,898,314,238
     ==21226== I1  misses:          291,271
     ==21226== L2i misses:          246,518
     ==21226== I1  miss rate:          0.00%
     ==21226== L2i miss rate:          0.00%
     ==21226==
     ==21226== D   refs:      3,264,359,909  (3,206,394,437 rd   + 57,965,472 wr)
     ==21226== D1  misses:       16,003,068  (   10,905,138 rd   +  5,097,930 wr)
     ==21226== L2d misses:        9,177,043  (    4,207,106 rd   +  4,969,937 wr)
     ==21226== D1  miss rate:           0.4% (          0.3%     +        8.7%  )
     ==21226== L2d miss rate:           0.2% (          0.1%     +        8.5%  )
     ==21226==
     ==21226== L2 refs:          16,294,339  (   11,196,409 rd   +  5,097,930 wr)
     ==21226== L2 misses:         9,423,561  (    4,453,624 rd   +  4,969,937 wr)
     ==21226== L2 miss rate:            0.1% (          0.0%     +        8.5%  )

 -}


text:
module STCR2Z2T0S0ReversalCornerPointSet where

import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           DFT.Plan
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           STC.PowerMethod
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoDecayStr:taoReversalStr:taoCornerStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:histFilePath:numIterationStr:writeSourceFlagStr:rStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      taoDecay = read taoDecayStr :: Double
      taoReversal = read taoReversalStr :: Double
      taoCorner = read taoCornerStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2Z2T0S0ReversalCornerPointSet"
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0ReversalCornerRadial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          taoDecay
          taoReversal
          taoCorner
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      radialArr
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  plan <- makeR2Z2T0S0Plan emptyPlan arrR2Z2T0S0
  let n = 30
      m = round $ (fromIntegral n) * (sqrt 2) / 2
      a = 10
      b = -10
      c = 10
      a' = round $ (fromIntegral a) * (sqrt 2) / 2
      b' = round $ (fromIntegral b) * (sqrt 2) / 2
      c' = round $ (fromIntegral c) * (sqrt 2) / 2
      r = read rStr :: Double
      rr = round r :: Int
      shift = 15
      numTheta = 45
      deltaTheta = 1 * pi / numTheta
      xs =
        -- [ R2S1RPPoint (rr, shift, 0, 1)
        -- , R2S1RPPoint (rr, -shift, 0, 1)
        -- , R2S1RPPoint (-rr, shift, 0, 1)
        -- , R2S1RPPoint (-rr, -shift, 0, 1)
        -- , R2S1RPPoint (shift, rr, 0, 1)
        -- , R2S1RPPoint (-shift, rr, 0, 1)
        -- , R2S1RPPoint (shift, -rr, 0, 1)
        -- , R2S1RPPoint (-shift, -rr, 0, 1)
        -- ]
        -- [R2S1RPPoint (i, -i, 0, 1) | i <- [-2,0 .. 20]] L.++
        [R2S1RPPoint (-i, -i, 0, 1) | i <- [-10,-8 .. 10]]
        -- ((L.map
        --     (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
        --     [ (r * cos (k * deltaTheta) + 0, r * sin (k * deltaTheta) + 0)
        --     | k <- [0 .. numTheta]
        --     ]) -- L.++
        --  -- (L.map
        --  --    (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
        --  --    [ ((r + 2) * cos (k * deltaTheta), (r + 2) * sin (k * deltaTheta))
        --  --    | k <- [0 .. numTheta - 1]
        --  --    ]) L.++
        --  -- (L.map
        --  --    (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
        --  --    [ ( (r + 5) * cos (k * deltaTheta)
        --  --      , (r + 5) * sin (k * deltaTheta))
        --  --    | k <- [0 .. numTheta - 1]
        --  --    ])
        --        -- L.++ [R2S1RPPoint (3, 5, 0, 1)]
        --  )
               -- ([R2S1RPPoint (i, i, 0, 1) | i <- [a',a' + b' .. c']] L.++
               --  [R2S1RPPoint (i, -i, 0, 1) | i <- [-a',-(a' + b') .. -c']] L.++
               --  [R2S1RPPoint (i, i, 0, 1) | i <- [-a',-(a' + b') .. -c']] L.++
               --  [R2S1RPPoint (i, -i, 0, 1) | i <- [a',a' + b' .. c']] L.++
               --  [R2S1RPPoint (i, 0, 0, 1) | i <- [a,a + b .. c]] L.++
               --  [R2S1RPPoint (0, i, 0, 1) | i <- [-a,-(a + b) .. -c]] L.++
               --  [R2S1RPPoint (i, 0, 0, 1) | i <- [-a,-(a + b) .. -c]] L.++
               --  [R2S1RPPoint (0, i, 0, 1) | i <- [a,a + b .. c]] L.++
               --  [R2S1RPPoint (3, 5, 0, 1)])
                -- [ R2S1RPPoint (n, 0, 0, 1)
                -- , R2S1RPPoint (0, n, 0, 1)
                -- , R2S1RPPoint (-n, 0, 0, 1)
                -- , R2S1RPPoint (0, -n, 0, 1)
                -- , R2S1RPPoint (m, m, 0, 1)
                -- , R2S1RPPoint (-m, m, 0, 1)
                -- , R2S1RPPoint (m, -m, 0, 1)
                -- , R2S1RPPoint (-m, -m, 0, 1)
                -- ]
  let bias = computeBiasR2T0S0 numPoint numPoint theta0Freqs scale0Freqs xs
      eigenVec =
        computeInitialEigenVectorR2T0S0
          numPoint
          numPoint
          theta0Freqs
          scale0Freqs
          thetaFreqs
          scaleFreqs
          xs
  powerMethodR2Z2T0S0
    plan
    folderPath
    numPoint
    numPoint
    numOrientation
    thetaFreqs
    theta0Freqs
    numScale
    scaleFreqs
    scale0Freqs
    arrR2Z2T0S0
    numIteration
    writeSourceFlag
    -- (printf "_%d" (round r :: Int))
    (printf
       "_%d_%d_%d_%d_%.2f_%.2f"
       (round maxScale :: Int)
       (round taoDecay :: Int)
       (round taoReversal :: Int)
       (round taoCorner :: Int)
       thetaSigma
       scaleSigma)
    0.5
    bias
    eigenVec


text:
-- |
-- Module      : SpriteClip.Types
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created September 7 2015

-- TODO | -
--        -

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
-- TODO: Bad idea to export all of these modules?
module SpriteClip.Types (module SpriteClip.Types,
                         module Data.Complex,
                         module T, module S, module M,
                         module Cairo,
                         module Gtk,
                         module Southpaw.Cartesian.Plane.BoundingBox) where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.Complex
import qualified Data.Set  as S
import qualified Data.Text as T
import qualified Data.Map  as M

import qualified Graphics.Rendering.Cairo as Cairo
import qualified Graphics.UI.Gtk          as Gtk

import Southpaw.Cartesian.Plane.BoundingBox

--------------------------------------------------------------------------------------------------------------------------------------------
-- Types
--------------------------------------------------------------------------------------------------------------------------------------------
-- Sprite sheets ---------------------------------------------------------------------------------------------------------------------------
-- |
-- type CutoutState = ([BoundingBox Double], Complex Double, Maybe (Complex Double))

-- |
-- TODO: Split up into logical, UI, graphics, etc. (or whatever grouping makes the most sense)
-- TODO: Rename '_cutouts' (?)
-- TODO: Coordinate transformations (scaling, panning)
data AppState = AppState { _cutouts :: [BoundingBox Double],                     --
                           _active  :: Maybe (BoundingBox Double, Complex Int),  -- The active cutout and marker (or Nothing)
                           _sheet   :: Cairo.Surface,                            -- The sprite sheet image
                           _mouse   :: Complex Double,                           -- The current position of the mouse
                           _keys    :: S.Set T.Text,                             -- The set of currently pressed keys
                           _pin     :: Maybe (Complex Double)                    --
                          --  _actionQueue :: []                                 --
                         }

-- |
-- data InputState = InputState { _mouse :: Complex Double,
--                                _keys  :: S.Set Key }

-- Input -----------------------------------------------------------------------------------------------------------------------------------
-- |
-- TODO: Rename (eg. 'Command') (?)
type Command = AppState -> IO AppState
-- data Action  = Action String (AppState -> IO AppState)

-- | Maps names to commands
-- type CommandPalette = M.Map String Command

-- |
-- TODO: Rename (?)
-- TODO: Figure out a better representation (?)
-- TODO: Composite types (?)
-- TODO: Mouse move/click/relative/scroll, key up/down/hold, drag'n'drop begin/end/motion/dropped,
-- type Key = String
type Event = S.Set String

-- |
data KeyBindings = KeyBindings (M.Map Event Command)

-- |
data Settings = Settings {}

-- | UI state machine
-- TODO: Rename
-- data InteractionState = Resizing Corner BoundingBox | Moving BoundingBox | Selected BoundingBox

-- Logical ---------------------------------------------------------------------------------------------------------------------------------
-- |
data Alignment f = Vertical f | Horizontal f deriving (Show, Eq)

max_stars_repo_name: SwiftsNamesake/SpriteClip

max_issues_repo_name: SwiftsNamesake/SpriteClip

max_forks_repo_name: SwiftsNamesake/SpriteClip


text:
-- | Markov chain weather change simulations
module Statistics.MarkovChain.MarkovChain
  ( steadyState1
  , steadyState2
  , steadyState3
  , steadyState4
  ) where

-- all modules required by the below
import Data.Complex                   -- base
import Data.Ord
import Numeric.LinearAlgebra.Data     -- hmatrix
import System.Random

import qualified Data.List as L
import qualified Numeric.LinearAlgebra.HMatrix as M

-- steady state of Markov chain -------------------------------------------------
-- we will write 4 approaches
-- 1. raise the state matrix to a high power, the limiting distribution is the
--    steady state
-- 2. solve the system of eqns: pi P = pi and sum_i=1^3 pi_i = 1
-- 3. using Perron-Frobenius theorem that eigenvector for the largest eval is
--    proportional to pi, so we find it and normalise it
-- 4. running a Monte Carlo simulation

-- The setup:
transition :: Matrix R
transition = (3><3)[0.5, 0.4, 0.1,
                    0.3, 0.2, 0.5,
                    0.5, 0.3, 0.2]

data Weather = Fine | Cloudy | Rain deriving (Eq, Show)

-- route 1:
steadyState1 :: Matrix R
steadyState1 = M.takeRows 1 $ foldl (M.<>) t [t | _ <- [1..n]]
  where
    n = 10^6
    t = transition
-- >>> disp 4 steadyState1
-- 1x3
-- 0.4375  0.3125  0.2500

-- route 2:
system :: Matrix R
system = takeRows 2 (tr (transition - ident 3)) === matrix 3 [1,1,1]

target :: Matrix R
target = (3><1) [0,0,1 :: R]

steadyState2 :: Maybe (Matrix R)
steadyState2 = M.linearSolve system target
-- >>> steadyState2
-- Just (3><1)
--  [ 0.4375
--  , 0.3125
--  ,   0.25 ]

-- route 3:
evals :: Vector (Complex Double)
evecs :: Matrix (Complex Double)
(evals, evecs) = M.eig (M.tr transition)

sizes :: [R]
sizes = map M.magnitude (toList evals)

maxid :: (Ord a, Num a, Enum a) => [a] -> (a, Int)
maxid xs = L.maximumBy (comparing fst) (zip xs [0..])

getEvec :: Matrix (Complex Double)
getEvec = M.dropColumns m $ M.takeColumns (m+1) evecs
  where
    m = snd $ maxid sizes

steadyState3 :: Matrix (Complex Double)
steadyState3 = getEvec / scalar (M.sumElements getEvec)
-- >>> disp 4 steadyState3
-- (3><1)
--  [ 0.43749999999999994 :+ 0.0
--  , 0.31250000000000006 :+ 0.0
--  , 0.24999999999999997 :+ 0.0 ]

-- route 4:
-- this is absolutely not the best way to do it at all, but its an intuitive
-- approach for those with little statistical training

-- replicate each rain type as many times it needs to be weighted as given,
-- this is in order to sample with the correct probability
genpop :: Matrix R -> Int -> [Weather]
genpop t s = replicate (round (10*t!s!0) :: Int) Fine ++
             replicate (round (10*t!s!1) :: Int) Cloudy ++
             replicate (round (10*t!s!2) :: Int) Rain

-- increment a count of each weather occurrence
increment :: [Int] -> Weather -> [Int]
increment [x,y,z] w = case w of
                        Fine   -> [x+1,   y,   z]
                        Cloudy -> [  x, y+1,   z]
                        Rain   -> [  x,   y, z+1]
increment xs _ = error $ "Needed a [x, y, z], got: " ++ show xs

wToInt :: Weather -> Int
wToInt Fine   = 0
wToInt Cloudy = 1
wToInt Rain   = 2

getNext :: Weather -> Int -> Weather
getNext w n = head . drop n . genpop transition . wToInt $ w

-- our pattern now becomes like:
-- getNext 3 (getNext 1 (getNext 4 (getNext 8 Fine)))

-- because we are now dealing with random numbers we have to move to IO town
steadyState4' :: Int -> IO [Double]
steadyState4' n = do
  g <- newStdGen
  let scores' = [fromIntegral x / fromIntegral n | x <- scores]
        where
          scores = L.foldl' increment [0,0,0] states
          states = init (L.scanl' getNext Fine steps)
          steps = take n $ randomRs (0, 9) g
  return scores'

steadyState4 :: IO [Double]
steadyState4 = steadyState4' (10^6)

max_stars_repo_name: karetsu/swh

max_issues_repo_name: karetsu/swh

max_forks_repo_name: karetsu/swh