text:
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE RankNTypes, GADTs, TypeFamilies, TypeOperators, ConstraintKinds #-}
{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-}
{-# LANGUAGE DataKinds, PolyKinds, ConstraintKinds #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE BangPatterns #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Language.Grappa.Distribution where

import Control.Monad
import Control.Monad.Trans
import GHC.Exts
import Data.Typeable
import qualified Data.Vector.Storable as V
--import Data.Number.LogFloat hiding (sum, product, log)

import Numeric.LinearAlgebra hiding (R, Uniform, (<>), (<#))
import qualified Numeric.LinearAlgebra as M

import qualified Numeric.Log as Log
import qualified Numeric.SpecFunctions as Gamma

import qualified Numeric.AD.Mode.Forward as ADF
import qualified Numeric.AD.Internal.Forward as ADF
import qualified Numeric.AD.Mode.Reverse as ADR
import qualified Numeric.AD.Internal.Reverse as ADR
import qualified Numeric.AD.Internal.Identity as ADR
import qualified Numeric.AD.Jacobian as ADR
import qualified Data.Reflection as ADR (Reifies)

-- | The R type alias
type R = Double

----------------------------------------------------------------------
-- * Supporting the Gamma and Beta Functions
----------------------------------------------------------------------

-- | Class for computing @ln (gamma x)@
class HasGamma a where
  logGamma :: a -> a
  digamma :: a -> a

instance HasGamma Double where
  logGamma = Gamma.logGamma
  digamma = Gamma.digamma

-- | The 'trigamma' function = the derivative of 'digamma' (which itself is the
-- derivative of 'logGamma'). This implementation was cribbed from
-- <https://stackoverflow.com/questions/2978979/how-to-improve-performance-of-this-numerical-computation-in-haskell/2983578 this>
-- Stack Overflow post
trigamma :: Double -> Double
trigamma x = go 0 (x + 5) $ computeP $ x + 6
  where
    go :: Int -> Double -> Double -> Double
    go !i !z !p
        | i >= 6    = p
        | otherwise = go (i+1) (z-1) (1 / (z*z) + p)
    invSq z = 1 / (z * z)
    computeP z = (((((5/66*p-1/30)*p+1/42)*p-1/30)*p+1/6)*p+1)/z+0.5*p
      where p = invSq z

----------------------------------------------------------------------
-- * Helper Definitions for Automatic Differentiation
----------------------------------------------------------------------

-- | Convert between two 'AD.Forward' types, using a conversion function which
-- should numerically be the identity, i.e., it should the representation of @r@
-- in type @a@ to the closest representable number to @r@ in @b@.
{-
unsafeConvertAD :: (a -> b) -> AD.Forward a -> AD.Forward b
unsafeConvertAD f (AD.Forward x y) =
  AD.Forward (f x) (f y)
unsafeConvertAD f (AD.Lift x) =
  AD.Lift (f x)
unsafeConvertAD _ AD.Zero = AD.Zero
-}

-- | Apply a function to an 'ADF.Forward' object that has a more efficient or
-- accurate version for the underlying type. For instance, log1p = log . (1+)
-- has a more accurate version in C code, but we still need to apply the above
-- mathematical definition to compute the derivative.
optimizedADOp :: Num a => (ADF.Forward a -> ADF.Forward a) -> (a -> a) ->
                 ADF.Forward a -> ADF.Forward a
optimizedADOp f f_opt x =
  ADF.bundle (f_opt $ ADF.primal x) (ADF.tangent $ f x)

-- Needed to use @'ADF.Forward' 'Double'@ Log
instance Log.Precise (ADF.Forward Double) where
  log1p = optimizedADOp Log.log1p (log . (1 +))
  expm1 = optimizedADOp Log.expm1 ((\x -> x - 1) . exp)

-- Needed to use @'ADR.Reverse' 'Double'@ Log
instance ADR.Reifies rs ADR.Tape => Log.Precise (ADR.Reverse rs Double) where
  log1p x = ADR.unary Log.log1p (ADR.Id (1 / (ADR.primal x + 1))) x
  expm1 x = ADR.unary Log.expm1 (ADR.Id ((ADR.primal x) * exp (ADR.primal x))) x

instance ADR.Reifies s ADR.Tape => HasGamma (ADR.Reverse s Double) where
  logGamma x =
    ADR.unary Gamma.logGamma (ADR.Id $ Gamma.digamma $ ADR.primal x) x
  digamma x =
    ADR.unary Gamma.digamma (ADR.Id $ trigamma $ ADR.primal x) x

----------------------------------------------------------------------
-- * Probabilities
----------------------------------------------------------------------

-- | We use the shorthand @Prob@ for the type of probabilities
newtype Prob = Prob { fromProb :: Log.Log Double }
             deriving (Num,Eq,Ord,Real,Fractional,Typeable)

instance Show Prob where
  show (Prob x) = show x

-- | Specialized method for testing if a 'Prob' is zero
probIsZero :: Prob -> Bool
probIsZero (Prob x) = x == 0

-- | Specialized method for testing if a 'Prob' is a NaN
probIsNaN :: Prob -> Bool
probIsNaN (Prob x) = isNaN $ Log.ln x

-- | Convert a (non-negative) real to a 'Prob'
rToProb :: R -> Prob
rToProb = Prob . Log.Exp . log

-- | Convert a real that is already in log-space to a 'Prob'
logRToProb :: R -> Prob
logRToProb = Prob . Log.Exp

-- | Convert a 'Prob' to a real (not in log space)
probToR :: Prob -> R
probToR = exp . Log.ln . fromProb

-- | Convert a 'Prob' to a real in log space
probToLogR :: Prob -> R
probToLogR = Log.ln . fromProb

-- | Lifting the 'Log.sum' method to 'Prob'
sumProb :: [Prob] -> Prob
sumProb ps = Prob $ Log.sum $ map fromProb ps

----------------------------------------------------------------------
-- * Matrices of reals and probabilities
----------------------------------------------------------------------

-- | A bundled up type for unboxed real-valued vectors
newtype RVector = RVector { unRVector :: Vector Double } deriving (Num,Eq,Show)

-- | Return the length of an 'RVector'
lengthV :: RVector -> Int
lengthV (RVector v) = V.length v

-- | Return the @i@th element of a 'RVector'
atV :: RVector -> Int -> Double
atV (RVector v) i = v V.! i

-- | Generate a probability vector of the given size
generateV :: Int -> (Int -> Double) -> RVector
generateV n f = RVector $ V.generate n f

-- | Map over an 'RVector'
mapV :: (R -> R) -> RVector -> RVector
mapV f (RVector v) = RVector (V.map f v)

-- | Fold over an 'RVector'
foldrV :: (R -> b -> b) -> b -> RVector -> b
foldrV f b (RVector v) = V.foldr' f b v

-- | Take the sum of an 'RVector'
sumV :: RVector -> Double
sumV (RVector v) = V.foldl' (+) 0 v

-- | Generate a vector from a list of its elements
fromListV :: [Double] -> RVector
fromListV = RVector . V.fromList

-- | Extract the elements of a vector as a list
toListV :: RVector -> [Double]
toListV = V.toList . unRVector

-- | The type of real-valued matrices
newtype RMatrix = RMatrix { unRMatrix :: Matrix Double } deriving (Num,Eq,Show)

-- | Return the number of rows of an 'RMatrix'
rowsM :: RMatrix -> Int
rowsM (RMatrix m) = rows m

-- | Return the number of columns of a probability matrix
colsM :: RMatrix -> Int
colsM (RMatrix m) = M.cols m

-- | Return the @(i,j)@th element of a 'RMatrix'
atM :: RMatrix -> Int -> Int -> Double
atM (RMatrix m) i j = atIndex m (i,j)

-- | Build a matrix from a list of its rows
fromRowsM :: [RVector] -> RMatrix
fromRowsM vs = RMatrix $ fromRows $ map unRVector vs

-- | Build a matrix of probabilities from a list of its columns
fromColsM :: [RVector] -> RMatrix
fromColsM vs = RMatrix $ fromColumns $ map unRVector vs

-- | Get the rows of a matrix as a list of vectors
toRowsM :: RMatrix -> [RVector]
toRowsM (RMatrix m) = map RVector $ toRows m

-- | Get the columns of a matrix as a list of vectors
toColsM :: RMatrix -> [RVector]
toColsM (RMatrix m) = map RVector $ toColumns m

-- | Matrix multiplication
mulM :: RMatrix -> RMatrix -> RMatrix
mulM (RMatrix m1) (RMatrix m2) = RMatrix (m1 M.<> m2)

-- | Matrix-vector multiplication
mulMV :: RMatrix -> RVector -> RVector
mulMV (RMatrix m) (RVector v) = RVector (m #> v)

-- | Vector-matrix multiplication
mulVM :: RVector -> RMatrix -> RVector
mulVM (RVector v) (RMatrix m) = RVector (v M.<# m)

-- | The type of (efficient storable) vectors of values in log space
newtype ProbVector = ProbVector { unProbVector :: Vector Double }
                   deriving (Eq, Show)

-- | Return the length of a 'ProbVector'
lengthPV :: ProbVector -> Int
lengthPV (ProbVector v) = V.length v

-- | Return the @i@th element of a 'ProbVector'
atPV :: ProbVector -> Int -> Prob
atPV (ProbVector v) i = Prob $ Log.Exp $ v V.! i

-- | Generate a probability vector of the given size
generatePV :: Int -> (Int -> Prob) -> ProbVector
generatePV n f = ProbVector $ V.generate n (Log.ln . fromProb . f)

-- | Generate a probability vector from a list of its elements
fromListPV :: [Prob] -> ProbVector
fromListPV = ProbVector . V.fromList . map probToLogR

-- | Extract the elements of a probability vector as a list
toListPV :: ProbVector -> [Prob]
toListPV = map logRToProb . V.toList . unProbVector

-- | Map over a 'ProbVector'
mapPV :: (Prob -> Prob) -> ProbVector -> ProbVector
mapPV f (ProbVector v) = ProbVector (V.map (probToLogR . f . logRToProb) v)

-- | Fold over a 'ProbVector'
foldrPV :: (Prob -> b -> b) -> b -> ProbVector -> b
foldrPV f b (ProbVector v) = V.foldr' (f . logRToProb) b v

-- | Apply a binary operator pointwise to two probability vectors
binOpPV :: (Prob -> Prob -> Prob) -> ProbVector -> ProbVector -> ProbVector
binOpPV f v1 v2 =
  generatePV (min (lengthPV v1) (lengthPV v2)) (\i -> f (atPV v1 i) (atPV v2 i))

instance Num ProbVector where
  (+) = binOpPV (+)
  (-) = binOpPV (-)
  (*) = binOpPV (*)
  abs = mapPV abs
  signum = mapPV signum
  fromInteger i = generatePV 1 (const $ fromInteger i)

-- | Take the sum of a 'ProbVector'. The algorithm for this is adapted from
-- 'Log.sum', though that function requires a 'Foldable' instance.
sumPV :: ProbVector -> Prob
sumPV (ProbVector v) | V.length v == 0 = Prob $ Log.Exp $ log 0
sumPV (ProbVector v) =
  let max_v = V.foldl1' max v in
  Prob $ Log.Exp $
  if isInfinite max_v then max_v else
    max_v + Log.log1p (V.foldl' (\r x -> r + Log.expm1 (x - max_v)) 0 v
                       + fromIntegral (V.length v - 1))

-- | The type of matrices of values in log space
newtype ProbMatrix = ProbMatrix (Matrix Double) deriving (Eq, Show)

-- | Return the number of rows of a probability matrix
rowsPM :: ProbMatrix -> Int
rowsPM (ProbMatrix m) = rows m

-- | Return the number of columns of a probability matrix
colsPM :: ProbMatrix -> Int
colsPM (ProbMatrix m) = M.cols m

-- | Return the @(i,j)@th element of a 'ProbMatrix'
atPM :: ProbMatrix -> Int -> Int -> Prob
atPM pm@(ProbMatrix m) i j =
  if i >= rowsPM pm then
    error ("atPM index out of bounds: row " ++ show i ++ " >= " ++ show (rowsPM pm))
    else
    if j >= colsPM pm then
      error ("atPM index out of bounds: column " ++ show j
             ++ " >= " ++ show (colsPM pm))
    else
      Prob $ Log.Exp $ atIndex m (i,j)

-- | Build a matrix of probabilities from a list of its rows
fromRowsPM :: [ProbVector] -> ProbMatrix
fromRowsPM vs = ProbMatrix $ fromRows $ map unProbVector vs

-- | Build a matrix of probabilities from a list of its columns
fromColsPM :: [ProbVector] -> ProbMatrix
fromColsPM vs = ProbMatrix $ fromColumns $ map unProbVector vs

-- | Get the rows of a matrix as a list of vectors
toRowsPM :: ProbMatrix -> [ProbVector]
toRowsPM (ProbMatrix m) = map ProbVector $ toRows m

-- | Get the columns of a matrix as a list of vectors
toColsPM :: ProbMatrix -> [ProbVector]
toColsPM (ProbMatrix m) = map ProbVector $ toColumns m

-- | Matrix-vector multiplication for probability matrices
mulPMV :: ProbMatrix -> ProbVector -> ProbVector
mulPMV m v =
  generatePV (rowsPM m) $ \i ->
  sumPV $ generatePV (min (colsPM m) (lengthPV v)) $ \j ->
  (atPM m i j) * (atPV v j)

-- | Vector-matrix multiplication
mulPVM :: ProbVector -> ProbMatrix -> ProbVector
mulPVM v m =
  generatePV (colsPM m) $ \j ->
  sumPV $ generatePV (min (rowsPM m) (lengthPV v)) $ \i ->
  (atPM m i j) * (atPV v i)

-- | Matrix-matrix multiplication for probability matrices
mulPM :: ProbMatrix -> ProbMatrix -> ProbMatrix
mulPM m1 (ProbMatrix m2) =
  ProbMatrix $ fromColumns $
  map (unProbVector . mulPMV m1 . ProbVector) $ toColumns m2

----------------------------------------------------------------------
-- * Distributions
----------------------------------------------------------------------

-- | This defines the support type of a distribution
type family Support (d :: *) :: *

-- | The base class for distributions, stating that we can evaluate the density
-- of any given distribution at a particular value.  For debugging purposes, we
-- also require distributions to be 'Show'able.
class PDFDist d where
  distDensity :: d -> Support d -> Prob

-- | This type class states that monad @m@ supports distribution type @d@ by
-- allowing distributions of type @d@ to be randomly sampled in @m@
class Monad m => SampleableIn m d where
  distSample :: d -> m (Support d)

instance (SampleableIn m d,
          Monad (t m), MonadTrans t) => SampleableIn (t m) d where
  distSample = lift . distSample

-- | This type family states that all distribution types in @ds@ satisfy @c@
type family DistsIn (ds :: [*]) (c :: * -> Constraint) :: Constraint where
  DistsIn '[] c = ()
  DistsIn (d ': ds) c = (c d, DistsIn ds c)

-- | This constraint says that the distribution type @d@ is continuous, i.e.,
-- isomorphic to the reals.
--
-- The class functions take @d@, not @Proxy d@, since e.g. we need
-- different implementations for @uniform 0 1 :: Uniform@ and @uniform
-- 1 2 :: Uniform@.
class Continuous d where
  -- | Note that 'toReal' need not be defined on all of @Support d@,
  -- and implementation should raise an error on out of bounds
  -- values. E.g., for @uniform 0 1@, we have @Support Uniform ~
  -- R@, but in fact 'toReal' is only defined on @[0,1]@. We
  -- don't use @Maybe R@ as the return type, since passing an out
  -- of bounds value to 'toReal' is probably a bug.
  toReal :: d -> Support d -> R
  fromReal :: d -> Double -> Support d

{- FIXME: this does not seem to be the derivative in the below, and is in fact
just fromReal in log-space. Figure this out or remove it!

  -- | The derivative of 'fromReal'.
  --
  -- In the multivariable case the derivative becomes the determinant
  -- of the Jacobian (matrix of partial derivatives).
  --
  -- We can factor this into a separate class if we end up with some
  -- continuous distributions for which we don't have Jacobian
  -- determinants.
  fromReal' :: d -> Double -> LogFloat
-}

-- | The PDF of a continuous dist transformed to be a dist on the
-- whole real line.
--
-- All of our single variable transforms are increasing, so I don't
-- think we have to worry about taking absolute values; indeed, the
-- @LogFloat@ is not even defined for negative values. (compare
-- formulas for single dimensional change of variables which uses
-- derivative and multiple dimension change of variables where the
-- *absolute value* of the Jacobian is used:
-- https://en.wikipedia.org/wiki/Integration_by_substitution).
{-
fromRealDistDensity :: (Continuous d, PDFDist d) => d -> Double -> Prob
fromRealDistDensity d x = distDensity d (fromReal d x) * Prob (fromReal' d x)
-}

----------------------------------------------------------------
-- * Helper functions for defining 'Continuous' instances
----------------------------------------------------------------

-- TODO(conathan): revisit these bijections once we have a better
-- understanding of any additional properties we want from them. We
-- plan to use these bijections to work in the reals and not have to
-- worry about bounds.
--
-- Potential trouble: consider the open interval @(lb,ub)@ bijection
-- below in @toRealLbUb@ and @fromRealLbUb@. The logit function @\p ->
-- log (p / (1 - p))@ is small for most of the reals, so if we do
-- anything where we choose random reals over a large range, nearly
-- all of them will cluster in the ends of underlying interval
-- (i.e. near @lb@ and @ub@). If we want more uniform choices in the
-- reals to map to more uniform choices in our underlying interval,
-- then we'll need a bijection that spreads the interval much wider.
--
-- Some bijections to consider (these need to be scaled to work with
-- @(lb,ub)@):
--
-- - @\x -> log (x / (1 - x))@:
--   https://www.wolframalpha.com/input/?i=plot+log+(x+/+(1+-+x))
--
-- - @\x -> (1 / (1 - x)) - (1 / x)@:
--   https://www.wolframalpha.com/input/?i=plot+(1/(1-x)+-+1/x)
--
-- - @\x -> tan x@:
--   https://www.wolframalpha.com/input/?i=plot+tan(x)
--
-- All of these can be made more uniform by scaling them by a
-- constant, e.g. @\x -> 100000 * tan x@.

{- FIXME: this is only used for fromReal'; see comments for it above
-- | The derivative of the identify function, returning @LogFloat@.
fromRealId' :: Double -> LogFloat
fromRealId' = logFloat
-}

-- | Biject the open interval @(lb,ub)@ onto the reals.
toRealLbUb :: Floating a => a -> a -> a -> a
toRealLbUb lb ub x = r
  where
    -- Biject @(lb,ub)@ onto @(0,1)@
    p = (x - lb) / (ub - lb)
    -- Biject @(0,1)@ onto the reals
    --r = log (p / (1 - p))
    r = (2*p - 1) / (p - p*p)

-- | Biject the reals on the open interval @(lb,ub)@
fromRealLbUb :: Floating a => a -> a -> a -> a
fromRealLbUb lb ub r = x
  where
    -- Inverse of @\p -> log (p / 1 - p)@, except that we have a special case
    -- for infinite values of exp(-r), to prevent us getting NaN values for AD
    -- exp_neg_r = exp (-r)
    -- p = if isInfinite exp_neg_r then 0 else 1 / (1 + exp_neg_r)
    p = (r - 2 + sqrt (r*r + 4)) / (2*r)
    -- Inverse of @\x -> (x - lb) / (ub - lb)@.
    x = p * (ub - lb) + lb

{-
-- | Derivative of 'fromRealLbUb'.
fromRealLbUb' :: R -> R -> R -> LogFloat
fromRealLbUb' lb ub r =
  -- The below is an expansion of
  -- @log ((ub - lb) * (exp (-r) / (1 + exp (-r))**2))@.
  --
  -- Here we use @log1p@ -- an implementation of @\p -> log (1 + p)@
  -- that's optimized for small @p@ -- from
  -- @Data.Numeric.LogFloat@. This implementation of @log1p@ will use
  -- an optimized C implementation if compiled with the @useffi@ flag,
  -- which is supposed to be enabled by default.
  --
  -- The @log-domain@ package that provides an native Haskell
  -- implementation of @log1p@, I think.
  logToLogFloat $ log (ub - lb) - r - 2 * log1p (exp (- r))
-}

-- | Biject the open interval @(0,\infty)@ onto the reals, using the piecewise
-- function that maps @x@ in @(0,1)@ to @log x@ and otherwise maps @x@ to @x-1@.
toRealZeroInftyPW :: Double -> Double
toRealZeroInftyPW x = if x < 1 then log x else x - 1

-- | Biject the reals on the open interval @(0,\infty)@, using the piecewise
-- function that maps negative @x@ to @exp x@ and otherwise maps @x@ to @x+1@.
fromRealZeroInftyPW :: Double -> Double
fromRealZeroInftyPW x = if x < 0 then exp x else x + 1

-- | Biject the open interval @(lb,\infty)@ onto the reals.
toRealLbInfty :: Double -> Double -> R
toRealLbInfty lb x = log (x - lb)

-- | Biject the reals onto the open interval @(lb,\infty)@.
fromRealLbInfty :: Double -> Double -> R
fromRealLbInfty lb r = exp r + lb

-- | Biject the reals onto the open interval @(lb,\infty)@ in log-space
fromRealLbInfty' :: Double -> Double -> Prob
fromRealLbInfty' lb r = logRToProb lb + logRToProb r

-- | Biject the open interval @(-\infty,ub)@ onto the reals.
toRealInftyUb :: Double -> Double -> R
toRealInftyUb ub x = - log (ub - x)

-- | Biject the reals onto the open interval @(-\infty,ub)@.
fromRealInftyUb :: Double -> Double -> R
fromRealInftyUb ub r = - (exp (- r) - ub)

-- | Biject the reals onto the open interval @(lb,\infty)@ in log-space
fromRealInftyUb' :: Double -> Double -> Prob
fromRealInftyUb' _ub _r =
  error "FIXME: fromRealInftyUb' is not right!"
  -- logRToProb ub + logRToProb (- r)

----------------------------------------------------------------------
-- * Example distribution types
----------------------------------------------------------------------

-- | Evaluate the density of a normal distribution at value @r@. This is
-- "unchecked" in that it assumes @sigma@ is non-negative.
normalDensityUnchecked :: Floating r => r -> r -> r -> Log.Log r
normalDensityUnchecked mu sigma x =
  Log.Exp $ - 0.5 * sq ((x - mu) / sigma) - log sigma - halfLogTwoPi
  where
    sq n = n * n
    halfLogTwoPi = 0.5 * log (2 * pi)

-- | Evaluate the density of a normal distribution at value @r@
normalDensityGen :: (Ord r, Floating r, Show r) => r -> r -> r -> Log.Log r
normalDensityGen mu sigma x =
  if sigma <= 0 then
    error $ "normalDensity: sigma is <= 0: sigma = " ++ show sigma
  else
    normalDensityUnchecked mu sigma x

-- | Evaluate the density of a normal distribution at value @r@
normalDensity :: Double -> Double -> Double -> Prob
normalDensity mu sigma x = Prob $ normalDensityGen mu sigma x

data StdNormal = StdNormal deriving Show
type instance Support StdNormal = R
instance PDFDist StdNormal where
  distDensity StdNormal x = logRToProb $ - 0.5 * sq x - halfLogTwoPi
    where
    sq n = n * n
    halfLogTwoPi = 0.5 * log (2 * pi)

instance Continuous StdNormal where
  toReal _ = id
  fromReal _ = id
  --fromReal' _ = fromRealId'

-- | The normal distribution
data Normal = Normal R R deriving Show
type instance Support Normal = R

instance PDFDist Normal where
  distDensity (Normal mu sigma) r = normalDensity mu sigma r

instance Continuous Normal where
  toReal _ = id
  fromReal _ = id
  --fromReal' _ = fromRealId'

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

-- | Multivariate normal distribution.
--
-- Our multivariate normal of @k@ dimensions is parameterized by the
-- mean vectur @mu = [mu1,...,muk]@ and the *Cholesky factor* @c@ of
-- the covariance matrix @Sigma@. I.e., for multivariate normal with
-- mean @mu@ and covariance matrix @Sigma@, we use
--
-- > MVNormal mu c
--
-- for
--
-- > transpose c * c = Sigma
--
-- The requirements for @c@ to be well formed are that
--
-- - @c@ is upper triangular.
--
-- - all diagonal entries of @c@ are positive.
--
-- There are several reasons to prefer the Cholesky factor @c@ to the
-- covariance matrix @Sigma@:
--
-- - computing the PDF is more efficient using @c@.
--
-- - the sampling algorithm uses @c@, so we need to compute it anyway.
--
-- - all upper triangular matrices with positive diagonal entries are
--   valid Cholesky factors, and a matrix is positive definite iff it
--   has a Cholesky factorization, and the Cholesky factorization is
--   unique. So, it's easy to uniformly sample the space of
--   multivariate normals by uniformly sampling the space upper
--   triangular matrices with positive diagonal entries. We care about
--   sampling multivariate normals e.g. when using multivariate
--   normals as a variational family in Black Box Variation Inference.
--
-- There is a more general definition of multivariate normal which
-- only requires non-negative entries on the diagonal of @c@, where
-- zero values correspond to zero variance random variables,
-- i.e. constants, but we don't support those here (the PDFs would be
-- Dirac deltas).
--
-- A covariance matrix @Sigma@ can be converted into a Cholesky factor
-- @c@ using the @Numeric.LinearAlgebra.chol@ function:
--
-- > c = chol Sigma
--
-- The covariance matrix @Sigma@ must be symmetric and positive
-- definite; you can use @Numeric.LinearAlgebra.trustsym@ to construct
-- it.
data MVNormal = MVNormal (Vector Double) (Matrix Double) deriving Show

type instance Support MVNormal = Vector Double

instance PDFDist MVNormal where
  distDensity (MVNormal mu c) x =
    mvnDensityMuC mu c x

instance (Monad m, SampleableIn m Normal) => SampleableIn m MVNormal where
  -- https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution
  distSample (MVNormal mu c) = do
    z <- V.fromList <$> (replicateM (size mu) $ distSample (Normal 0 1))
    -- Our Cholesky factorization is @Sigma = transpose c * c@ and
    -- Wikipedia has @Sigma = c * transpose c@, so we need to
    -- transpose here.
    return $ mu + (tr c #> z)

-- | Multivariate normal density in terms of Cholesky factor.
mvnDensityMuC :: Vector Double -> Matrix Double -> Vector Double -> Prob
mvnDensityMuC mu c x =
  if not isValidCholeskyFactor
  then error "distDensity: Cholesky factor is malformed!"
  -- https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Likelihood_function.
  else logRToProb $ -0.5 * (logAbsDet + prod + k * log2pi)
  where
    log2pi = log (2 * pi)
    k = fromIntegral $ size x
    -- We have
    --
    -- > Sigma = transpose c * c
    --
    -- and so
    --
    -- > det Sigma = det c * det c .
    --
    -- We know @det c > 0@ since the diag is positive.
    logAbsDet = 2 * (log $ V.product (takeDiag c))
    -- Here @cholSolve@ solves a linear system for the matrix @Sigma@
    -- for which @c@ is the Cholesky factor. I.e. if
    --
    -- > u == cholSolve m v
    --
    -- then
    --
    -- > transpose c * c * v == u .
    prod = (x - mu) <.> flatten (cholSolve c (asColumn $ x - mu))

    isValidCholeskyFactor = all (> 0) (V.toList $ takeDiag c)

-- | Multivariate normal density in terms of covariance matrix @Sigma@.
--
-- Can use this to test the other implementation above; see
-- 'testMvDensity' below.
mvnDensityMuSigma :: Vector Double -> Matrix Double -> Vector Double -> Prob
mvnDensityMuSigma mu sigma x =
  if signDet < 0
  then error "distDensity: matrix is not positive definite!"
  -- https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Likelihood_function.
  else logRToProb $ -0.5 * (logAbsDet + prod + k * log2pi)
  where
    log2pi = log (2 * pi)
    k = fromIntegral $ size x
    (sigmaInv, (logAbsDet, signDet)) = invlndet sigma
    prod = (x - mu) <.> (sigmaInv #> (x - mu))

-- | Check that the two mv normal density functions agree.
--
-- They agree if this function returns a small number.
testMvnDensity :: Vector Double -> Matrix Double -> Vector Double -> Double
testMvnDensity mu c x =
   abs $ (probToR $ mvnDensityMuSigma mu (tr c <> c) x) -
         (probToR $ mvnDensityMuC mu c x)

-- | Should return a smalllllll number.
testMvnDensityExample :: Double
testMvnDensityExample = sum [ testMvnDensity mu c x | x <- xs ]
  where
    mu = V.fromList [1,2]
    c = chol (trustSym $ (2><2) [1,0.5,0.5,1])
    xs = [ V.fromList [x0,x1] | x0 <- [-5,-4.5..5] , x1 <- [-5,-4.5..5] ]

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

-- | The uniform distribution
data Uniform = Uniform Double R deriving Show
type instance Support Uniform = R

uniformDensityGen :: (Ord r, RealFloat r, Log.Precise r) =>
                     r -> r -> r -> Log.Log r
uniformDensityGen lb ub r
  | lb >= ub =
    error "Malformed uniform distribution! Lower bound >= upper bound!"
  | r >= lb && r < ub =
    -- The expression 1 / (ub - lb) in log-space
      Log.Exp $ negate $ log $ ub - lb
  | otherwise = 0

uniformDensity :: Double -> Double -> Double -> Prob
uniformDensity lb ub r
  | lb >= ub =
    error "Malformed uniform distribution! Lower bound >= upper bound!"
  | r >= lb && r < ub = 1 / rToProb (ub - lb)
  | otherwise = 0

instance PDFDist Uniform where
  distDensity (Uniform lb ub) = uniformDensity lb ub

instance Continuous Uniform where
  toReal (Uniform lb ub) = toRealLbUb lb ub
  fromReal (Uniform lb ub) = fromRealLbUb lb ub
  --fromReal' (Uniform lb ub) = fromRealLbUb' lb ub

data StudentT = StudentT Double deriving Show
type instance Support StudentT = R

instance Continuous StudentT where
  toReal _ = id
  fromReal _ = id
  --fromReal' _ = fromRealId'

data Cauchy = Cauchy deriving Show
type instance Support Cauchy = R
instance PDFDist Cauchy where
  distDensity Cauchy x = logRToProb $ - log (1 + sq x) - log pi
    where sq n = n * n

instance Continuous Cauchy where
  toReal _ = id
  fromReal _ = id
  --fromReal' _ = fromRealId'

data For d = forall t. For [t] (t -> d)
type instance Support (For d) = [Support d]
instance PDFDist d => PDFDist (For d) where
  distDensity (For xs f) ys =
    product [distDensity (f x) y | (x,y) <- zip xs ys]

data StdUniform = StdUniform deriving Show
type instance Support StdUniform = R

instance PDFDist StdUniform where
  distDensity _ = uniformDensity 0 1

instance Continuous StdUniform where
  toReal _ = toRealLbUb 0 1
  fromReal _ = fromRealLbUb 0 1
  --fromReal' _ = fromRealLbUb' 0 1

-- | The categorical distribution, that picks a natural number between 1 and @n@
-- from a list of probabilisties of each number. Note that the probabilities
-- need not be normalized, e.g., the list @[1,1,1]@ is treated the same as
-- @[1/3,1/3,1/3]@, i.e., that each of 0, 1, and 2 has equal probability.
data Categorical = Categorical [Prob] deriving Show
type instance Support Categorical = Int

instance PDFDist Categorical where
  distDensity (Categorical ws) n =
    let total_w = sum ws in
    if n >= length ws then 0 else ws!!n / total_w

data Exponential = Exponential R deriving Show
type instance Support Exponential = R

exponentialDensityUnchecked :: (RealFloat a, Log.Precise a) =>
                               a -> a -> Log.Log a
exponentialDensityUnchecked rate x = Log.Exp $ log rate - rate * x

exponentialDensity :: (Ord a, RealFloat a, Log.Precise a) => a -> a -> Log.Log a
exponentialDensity rate x =
  if rate > 0 then exponentialDensityUnchecked rate x else 0

instance PDFDist Exponential where
  distDensity (Exponential rate) = Prob . exponentialDensity rate

data Beta = Beta R R deriving Show
type instance Support Beta = R

betaDensityUnchecked :: (HasGamma a, Floating a) => a -> a -> a -> Log.Log a
betaDensityUnchecked alpha beta x =
  Log.Exp $
  ((alpha - 1) * log x) + ((beta - 1) * log (1 - x)) -
  (logGamma alpha + logGamma alpha - logGamma (alpha + beta))

betaDensity :: (Ord a, HasGamma a, RealFloat a, Log.Precise a) =>
               a -> a -> a -> Log.Log a
betaDensity alpha beta x =
  if x > 0 && x < 1 then betaDensityUnchecked alpha beta x else 0

-- | Compute the PDF of the Beta distribution where the argument is already in
-- log space
betaDensityLog :: (Ord a, HasGamma a, RealFloat a, Log.Precise a) =>
                  a -> a -> Log.Log a -> Log.Log a
betaDensityLog alpha beta x =
  if x > 0 && x < 1 then
    Log.Exp $
    ((alpha - 1) * Log.ln x) + ((beta - 1) * Log.ln (1 - x)) -
    (logGamma alpha + logGamma alpha - logGamma (alpha + beta))
  else 0

instance PDFDist Beta where
  distDensity (Beta alpha beta) = Prob . betaDensity alpha beta

data Gamma = Gamma R R deriving Show
type instance Support Gamma = R

-- | Calculate the PDF of the gamma distribution, given shape parameter @k@ and
-- scale parameter @theta@, but without the @1 / (gamma k * theta ** k)@
-- normalization constant
relativeGammaDensity :: Floating a => a -> a -> a -> Log.Log a
relativeGammaDensity k theta x =
  Log.Exp $ (k - 1) * log x - (x / theta)

-- | Calculate the PDF of the gamma distribution, assuming @x > 0@, given shape
-- parameter @k@ and scale parameter @theta@
gammaDensityUnchecked :: (HasGamma a, Floating a) => a -> a -> a -> Log.Log a
gammaDensityUnchecked k theta x =
  Log.Exp (Log.ln (relativeGammaDensity k theta x) - logGamma k - k * log theta)

-- | Calculate the PDF of the gamma distribution, given shape parameter @k@ and
-- scale parameter @theta@
gammaDensity :: (Ord a, HasGamma a, RealFloat a, Log.Precise a) =>
                a -> a -> a -> Log.Log a
gammaDensity k theta x =
  if x <= 0 then 0 else gammaDensityUnchecked k theta x

instance PDFDist Gamma where
  distDensity (Gamma k theta) = Prob . gammaDensity k theta

data Dirichlet = Dirichlet [R] deriving Show
type instance Support Dirichlet = [R]

-- | The log of the multivariate beta function
logMVBeta :: (HasGamma a, Floating a) => [a] -> Log.Log a
logMVBeta alphas =
  -- Gamma (a_1) * ... * Gamma (a_n) / Gamma (a_1 + ... + a_n) in log space
  Log.Exp $ sum (map logGamma alphas) - logGamma (sum alphas)

-- | Calculate the density of the Dirichlet distribution at any type that
-- supports the gamma function
dirichletDensity :: (HasGamma a, Floating a) => [a] -> [a] -> Log.Log a
dirichletDensity alphas xs =
  -- x_1 ** (alpha_1 - 1) * ... * x_n ** (alpha_n - 1) / Beta (alphas)
  Log.Exp $
  sum (zipWith (\alpha x -> (alpha - 1) * log x) alphas xs) -
  Log.ln (logMVBeta alphas)

-- | Calculate the density of the Dirichlet distribution over log space
dirichletDensityLog :: (HasGamma a, Floating a) => [a] -> [Log.Log a] ->
                       Log.Log a
dirichletDensityLog alphas xs =
  -- x_1 ** (alpha_1 - 1) * ... * x_n ** (alpha_n - 1) / Beta (alphas)
  Log.Exp $
  sum (zipWith (\alpha x -> (alpha - 1) * Log.ln x) alphas xs) -
  Log.ln (logMVBeta alphas)

-- | The log of the multivariate beta function
logMVBetaV :: RVector -> Log.Log Double
logMVBetaV (RVector alphas) =
  -- Gamma (a_1) * ... * Gamma (a_n) / Gamma (a_1 + ... + a_n) in log space
  Log.Exp $ V.foldl' (\r alpha -> r + logGamma alpha) 0 alphas

-- | Calculate the density of the Dirichlet distribution over an 'RVector'
dirichletDensityV :: RVector -> RVector -> Prob
dirichletDensityV (RVector alphas) (RVector xs) =
  -- x_1 ** (alpha_1 - 1) * ... * x_n ** (alpha_n - 1) / Beta (alphas)
  Prob $ Log.Exp $
  V.ifoldl' (\r i alpha -> r + (alpha - 1) * log (xs!i)) 0 alphas -
  Log.ln (logMVBetaV $ RVector alphas)

-- | Calculate the density of the Dirichlet distribution over a 'ProbVector'
dirichletDensityPV :: RVector -> ProbVector -> Prob
dirichletDensityPV (RVector alphas) (ProbVector xs) =
  -- x_1 ** (alpha_1 - 1) * ... * x_n ** (alpha_n - 1) / Beta (alphas)
  Prob $ Log.Exp $
  V.ifoldl' (\r i alpha -> r + (alpha - 1) * xs!i) 0 alphas -
  Log.ln (logMVBetaV $ RVector alphas)

instance PDFDist Dirichlet where
  distDensity (Dirichlet alphas) xs = Prob $ dirichletDensity alphas xs

data Bernoulli = Bernoulli Double deriving Show

data Binomial = Binomial Int Double deriving Show

data Wishart = Wishart (Matrix Double) Double deriving Show
type instance Support Wishart = Matrix Double
instance PDFDist Wishart where
  distDensity (Wishart _scaleChol _df) _x = undefined

-- | A dirac delta distribution, that has a 100 percent chance of returning a
-- given constant value
data Dirac a = Dirac a deriving Show
type instance Support (Dirac a) = a

instance Eq a => PDFDist (Dirac a) where
  -- FIXME: this density is not really correct for R; it should really be
  -- infinite for continuous functions and 1 for discrete ones
  distDensity (Dirac a) x = if a == x then 1 else 0

-- | A technically invalid distribution, that always returns the value @x@ but
-- always has a score of @1@
data DontCare a = DontCare a deriving Show
type instance Support (DontCare a) = a

instance PDFDist (DontCare a) where
  distDensity (DontCare _) _ = 1

----------------------------------------------------------------------
-- * Continuity
----------------------------------------------------------------------

class Continuous (d f) => ReprContinuous d f
instance Continuous (d f) => ReprContinuous d f

-- | This type class states that monad @m@ supports distribution type @d@ by
-- allowing distributions of type @d@ to be randomly sampled in @m@
class Monad m => ReprSampleableIn m d f where
  distSampleRepr :: d f -> m (f (ReprSupport d))

class DiffContinuous d f where
  toRealF   :: d f -> f (ReprSupport d) -> f R
  fromRealF :: d f -> f R -> f (ReprSupport d)

class DiffPDF d f where
  distDensityF :: d f -> f (ReprSupport d) -> f Prob

{-
instance DiffContinuous ReprNormal AD.Forward where
  toRealF _ d = d
  fromRealF _ d = d

instance DiffPDF ReprNormal AD.Forward where
  distDensityF (ReprNormal mu sigma) x =
    if sigma <= 0
      then error $ "distDensityAD Normal: sigma is <= 0: sigma = " ++ show sigma
           -- FIXME: instead of using exp below, figure out how to do an
           -- adToProb that takes a real in log space and preserves
      else adToProb $ exp $ -0.5 * sq((x - mu) / sigma) - log sigma - halfLogTwoPi
      where sq n = n * n
            halfLogTwoPi = 0.5 * log (2 * pi)

instance DiffContinuous ReprUniform AD.Forward where
  toRealF (ReprUniform lb ub) =
    toRealLbUbNum lb ub
  fromRealF (ReprUniform lb ub) =
    fromRealLbUbNum lb ub

instance DiffPDF ReprUniform AD.Forward where
  distDensityF (ReprUniform lb ub) _
    | lb >= ub = error "Malformed uniform distribution! Lower bound >= upper bound!"
    | otherwise =
      (1 / adToProb (lb - ub))
-}

----------------------------------------------------------------------
-- Entropy
----------------------------------------------------------------------

-- | FIXME: document this!
class Entropy d where
  entropy :: d -> Double

instance Entropy Normal where
  entropy (Normal _mu sigma) = k + log sigma
    where
    k = 0.5 * (1 + log (2 * pi))

----------------------------------------------------------------------
-- * The ReprDistribution stuff
----------------------------------------------------------------------

-- | This defines the support type of a distribution
type family ReprSupport (d :: (* -> *) -> *) :: *

data ReprNormal f = ReprNormal (f R) (f R)
type instance ReprSupport ReprNormal = R

data ReprUniform f = ReprUniform (f R) (f R)
type instance ReprSupport ReprUniform = R

data ReprCauchy f = ReprCauchy
type instance ReprSupport ReprCauchy = R

data ReprCategorical f = ReprCategorical [f Prob]
type instance ReprSupport ReprCategorical = Int

-- FIXME: this uses different vector and matrix types than MVNormal
data ReprMVNormal f = ReprMVNormal (f RMatrix) (f RMatrix)
type instance ReprSupport ReprMVNormal = RMatrix

max_stars_repo_name: GaloisInc/grappa

max_issues_repo_name: GaloisInc/grappa

max_forks_repo_name: GaloisInc/grappa


text:
module MachineLearning.PCATest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx

import qualified Numeric.LinearAlgebra as LA

import MachineLearning.PCA

getDimReducerSmokeTest =
  let nFeatures = 4
      nExamples = 7
      m = LA.matrix nFeatures [1 .. fromIntegral $ nFeatures*nExamples]
      m10 = m * 10
      (reduceDims, retainedVariance, mReduced) = getDimReducer m 2
      m10Reduced = reduceDims m10
  in do
    assertEqual "dimension equality (getDimReducer)" (LA.cols mReduced) 2
    assertEqual "dimension equality (reduceDims)" (LA.cols m10Reduced) 2
    assertApproxEqual "retained variance" 1e-10 1 retainedVariance

getDimReducer_rvSmokeTest rv=
  let nFeatures = 4
      nExamples = 7
      m = LA.matrix nFeatures [1 .. fromIntegral $ nFeatures*nExamples]
      m10 = m * 10
      (reduceDims, k, mReduced) = getDimReducer_rv m rv
      m10Reduced = reduceDims m10
  in do
    assertEqual "dimension equality (getDimReducer_rv)" (LA.cols mReduced) k
    assertEqual "dimension equality (reduceDims_rv)" (LA.cols m10Reduced) k
    assertBool "reduced number of dimensions" $ k <= nFeatures

tests = [ testGroup "smoke test" [
            testCase "getDimReducer" getDimReducerSmokeTest
            , testCase "getDimReducer_rv, rv = 0" $ getDimReducer_rvSmokeTest 0
            , testCase "getDimReducer_rv, rv = 0.5" $ getDimReducer_rvSmokeTest 0.5
            , testCase "getDimReducer_rv, rv = 1" $ getDimReducer_rvSmokeTest 1
            , testCase "getDimReducer_rv, rv = 2" $ getDimReducer_rvSmokeTest 2
            ]
        ]

ext: hs

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Day3 (solvePart1, solvePart2) where

import           Data.Complex
import           Data.Map     (Map, (!))
import qualified Data.Map     as Map

createMovementList :: [Complex Double] -> [Int] -> [Complex Double]
createMovementList (f:s:rest) (x:xs) = replicate x f ++ replicate x s ++ createMovementList rest xs

movementList :: [Complex Double]
movementList = createMovementList (iterate (*(0 :+ 1)) (1 :+ 0)) [1..]

positions :: [Complex Double]
positions = scanl (+) (0 :+ 0) movementList

manhattenDistance :: Complex Double -> Int
manhattenDistance c = round (sum (fmap abs c))

solvePart1 :: String -> Int
solvePart1 input = let
        n = read input
    in manhattenDistance (positions !! (n-1))

solvePart2 :: String ->  Int
solvePart2 input = let
        n = read input
    in populateMap n

populateMap :: Int -> Int
populateMap n = let
        foobar (p:ps) m = let
                v = eval p m
            in if v <= n
                then foobar ps (Map.insert p v m)
                else v
    in foobar positions Map.empty

eval :: Complex Double -> Map (Complex Double) Int -> Int
eval c m
    | Map.null m = 1
    | otherwise  = sum (map snd (neighbours c m))

neighbours :: Complex Double -> Map (Complex Double) Int -> [(Complex Double, Int)]
neighbours (a :+ b) m = let
        member k = Map.member k m
        entry k = (k, m ! k)
        neighbourhood = [(a+x) :+ (b+y) | x <- [-1, 0, 1], y <- [-1, 0, 1], x /= 0 || y /= 0]
    in map entry (filter member neighbourhood)

instance (Ord a) => Ord (Complex a) where
   (a :+ b) `compare` (c :+ d)
        | a > c = GT
        | a < c = LT
        | b > d = GT
        | b < d = LT
        | b == d = EQ


text:
module Math.Probably.Student where

import Numeric.LinearAlgebra
import Math.Probably.FoldingStats
import Text.Printf

--http://www.haskell.org/haskellwiki/Gamma_and_Beta_function
--cof :: [Double]
cof = [76.18009172947146,-86.50532032941677,24.01409824083091,
       -1.231739572450155,0.001208650973866179,-0.000005395239384953]

--ser :: Double
ser = 1.000000000190015

--gammaln :: Double -> Double
gammaln xx = let tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5)
                 ser' = foldl (+) ser $ map (\(y,c) -> c/(xx+y)) $ zip [1..] cof
             in -tmp' + log(2.5066282746310005 * ser' / xx)

beta z w = exp (gammaln z + gammaln w - gammaln (z+w))

fac :: (Enum a, Num a) => a -> a
fac n = product [1..n]

ixbeta :: (Enum a, Floating a) => a -> a -> a -> a
ixbeta x a b = let top = fac $ a+b-1
                   down j = fac j * fac (a+b-1-j)
               in sum $ map (\j->(top/down j)*(x**j)*(1-x)**(a+b-1-j)) [a..a+b-1]

studentIntegral :: (Enum a, Floating a) => a -> a -> a
studentIntegral t v = 1-ixbeta (v/(v+t*t)) (v/2) (1/2)

tTerms :: Vector Double
tTerms = fromList $ map tTermUnmemo [1..100]

tTermUnmemo :: Int -> Double
tTermUnmemo nu = gammaln ((realToFrac nu+1)/2) - log(realToFrac nu*pi)/2 - gammaln (realToFrac nu/2)

tTerm1 :: Int -> Double
tTerm1 df | df <= 100 = tTerms@>df
          | otherwise = tTermUnmemo df

tDist :: Int -> Double -> Double
tDist df t = tTerm1 df - (realToFrac df +1/2) * log (1+(t*t)/(realToFrac df))

tDist3 :: Double -> Double -> Int -> Double -> Double
tDist3 mean prec df x
    =   tTerm1 df
      + log(prec)/2
      - (realToFrac df +1/2) * log (1+(prec*xMinusMu*xMinusMu)/(realToFrac df))
 where xMinusMu = x-mean

oneSampleT :: Floating b => b -> Fold b b
oneSampleT v0 = fmap (\(mean,sd,n)-> (mean - v0)/(sd/(sqrt n))) meanSDNF

pairedSampleT :: Fold (Double, Double) Double
pairedSampleT  = before (fmap (\(mean,sd,n)-> (mean)/(sd/(sqrt n))) meanSDNF)
                        (uncurry (-))

data TTestResult = TTestResult { degFreedom :: Int,
                                 tValue :: Double,
                                 pValue :: Double }

ppTTestRes :: TTestResult -> String
ppTTestRes (TTestResult df tval pval) = printf "paired t(%d)=%.3g, p=%.3g" df tval pval

pairedTTest :: [(Double,Double)] -> TTestResult
pairedTTest vls =
    let tval = runStat pairedSampleT vls
        df =  length vls - 1
        pval =  (1-) $ studentIntegral (tval) (realToFrac df)
    in TTestResult df tval pval

oneSampleTTest :: [Double] -> TTestResult
oneSampleTTest vls =
  let tval = runStat (oneSampleT 0) vls
      df =  length vls - 1
      pval =  (1-) $ studentIntegral (tval) (realToFrac df)
  in TTestResult df tval pval

max_stars_repo_name: glutamate/probably-baysig

max_issues_repo_name: silky/probably-baysig

max_forks_repo_name: silky/probably-baysig


text:
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FunctionalDependencies #-}
module VAD(
      vad
    ) where

--  let framed = frameWithWinAndOverlap 256 hann 10 s
--      sp = mapS (fft . U.map (:+ 0)) framed
--

import Prelude hiding(splitAt,(:),foldr1,tail,(++))
import Internal
import Signal
import qualified Data.Vector.Unboxed as U
import Data.Vector.Unboxed((!),Unbox(..))
import Data.List.Stream
import Common
import Windows
import Data.Complex
import Transform
import qualified Trace as F
import Fixed
import GHC.TypeLits
import SpecialInt
import Data.Int

ltseF :: Sample a
      => Int
      -> Signal (U.Vector a)
      -> Signal (U.Vector a)
ltseF n (Signal s) =
    let (before,remaining) = splitAt n s
        _lste b (h:t) =
            let (future,tl) = splitAt n t
                valueMax = U.zipWith max (foldr1 (U.zipWith max) (h:future))  (foldl1' (U.zipWith max) b)
            in
            valueMax:(_lste (tail before ++ [h]) t)
    in
    Signal (_lste before remaining)

lstdD :: Int
      -> U.Vector Double
      -> U.Vector Double
      -> Double
lstdD winSize noiseEnergy lste =
        let d a b | b == 0 = fromDouble 1e12
                  | otherwise = a / b
        in
        10 * log(U.sum (U.zipWith (d) lste noiseEnergy) / fromIntegral winSize) / log 10

theTotal :: (SingI r, SingI s, SingI n)
         => Fixed Int32 n s r
         -> U.Vector (Fixed Int32 n s r)
         -> U.Vector (Fixed Int32 n s r)
         -> Fixed Int40 n s r
theTotal m lste noiseEnergy =
           let  d a b | b == 0 = convert m
                      | otherwise = convert a / convert b
           in
           U.sum (U.zipWith (d) lste noiseEnergy)
{-# INLINE [0] theTotal #-}

lstdF :: (SingI n, SingI s, SingI r)
      => Int
      -> U.Vector (Fixed Int32 n s r)
      -> U.Vector (Fixed Int32 n s r)
      -> (Fixed Int32 n s r)
lstdF winSize noiseEnergy lste = 10 * log(convert $ theTotal maxBound lste noiseEnergy/ fromIntegral winSize) / log 10

getDecisionF :: (SingI n, SingI s, SingI r, SingI (n + n))
             => Int
             -> U.Vector (Fixed Int32 (n + n) s r)
             -> [(U.Vector (Fixed Int32 (n + n) s r), U.Vector (Fixed Int32 (n + n) s r) )]
             -> [Fixed Int16 n s r]
getDecisionF winSize energy ((c,currentE):r) =
    let tv = fromDouble 31.0
        tn = fromDouble 27.0
        l = lstdF winSize energy c
        result | l >= tv = 1 : getDecisionF winSize energy r
               | l <= tn = 0 : getDecisionF winSize currentE r
               | otherwise = 0 : getDecisionF winSize energy r
    in
    result
getDecisionF winSize energy [] = 0:getDecisionF winSize energy []

getDecisionD :: Int
             -> U.Vector Double
             -> [(U.Vector Double, U.Vector Double)]
             -> [Double]
getDecisionD winSize energy ((c,currentE):r) =
    let tv = fromDouble 31.0
        tn = fromDouble 27.0
        l = lstdD winSize energy c
        result | l >= tv = 1 : getDecisionD winSize energy r
               | l <= tn = 0 : getDecisionD winSize currentE r
               | otherwise = 0 : getDecisionD winSize energy r
    in
    result
getDecisionD winSize energy [] = 0:getDecisionD winSize energy []

bandEnergy :: Complex Double -> Double
bandEnergy (x :+ y) = x*x + y*y

bandEnergyF :: (SingI n, SingI s, SingI r, SingI (n + n)) => Complex (Fixed Int16 n s r) -> Fixed Int32 (n + n) s r
bandEnergyF (x :+ y) = amul x x + amul y y

class VAD a where
    vad :: (Sample a, FFT a)
        => Sampled Time a
        -> Sampled Time a

instance (SingI n, SingI s, SingI r, SingI (n + n)) => VAD (Fixed Int16 n s r) where
    vad s =
        let winSize = 256
            overlap = 20
            n = 2
            framed = frameWithWinAndOverlap winSize overlap hann s
            energy = mapS (U.map bandEnergyF . fft . U.map (:+ 0)) (getSignal framed)
            noiseEnergy0 = U.generate winSize (const (fromDouble 0.00001))
            lt = ltseF n energy
        in
        --mapS (lstd winSize noiseEnergy0) lt
        Sampled (period framed) (onSamples (getDecisionF winSize noiseEnergy0) (zipS lt (dropS n energy)))

instance VAD Double where
    vad s =
        let winSize = 256
            overlap = 20
            n = 2
            framed = frameWithWinAndOverlap winSize overlap hann s
            energy = mapS (U.map bandEnergy . fft . U.map (:+ 0)) (getSignal framed)
            noiseEnergy0 = U.generate winSize (const (fromDouble 0.00001))
            lt = ltseF n energy
        in
        --mapS (lstd winSize noiseEnergy0) lt
        Sampled (period framed) (onSamples (getDecisionD winSize noiseEnergy0) (zipS lt (dropS n energy)))

max_stars_repo_name: alpheccar/Signal

max_issues_repo_name: cpehle/Signal

max_forks_repo_name: cpehle/Signal


text:
{-#Language TupleSections, MultiWayIf#-}
module Stratification where

-- Внешние импорты
import Data.Matrix (Matrix(..))

import Data
import Point
import Data.Array as A
import Data.Map as M
import Data.Complex
import SymbolicImage
import Jac (jac, e)

str :: F Point -> Imagination -> Int -> Stratification
str fp im i = iterate
    (stepStratification fp)
    (toStratification im) !! i

-- Один шаг с подразбиением.
stepStratification :: F Point-> F Stratification
stepStratification fp (Stratification d1 d2 ls) = Stratification
    d1 (d2*2) $ concatMap celling $ cyclics $
    (shiftCeil d1 d2 fp #=) <$> ls
  where

{-
  where
    -- считаем граф
    !gr     = (shiftCeil d1 d2 fp ##=) <$> ls
    -- индех вершин
    !index  = formArray gr
    -- обратный индех
    !indexR = formMap index

    unpack !gr = (index A.!) <$> gr
    pack !gr = do
        (a, b) <- gr
        let a' = fromJust $ a`M.lookup`indexR
        return (a', a', MB.mapMaybe (`M.lookup`indexR) b)

(##=) !f !c = (c, f c)

-}

-- вспомогательные структуры данных
-- массив с ячейками п индексам.
formArray :: [(Ceil3, [Ceil3])] -> Array Int Ceil3
formArray gr = array (1, length a) (zip [1..] a)
  where a = plat gr

-- находим индекс по ячейке.
formMap :: Array Int Ceil3 -> Map Ceil3 Int
formMap  arr = M.fromList $ (\(a,b) -> (b,a)) <$> A.assocs arr

plat :: Ord c => [(c, [c])] -> [c]
plat a = myNub.concat $ (\(a, b) -> a:b) <$> a

toStratification :: Imagination -> Stratification
toStratification (Imagination d c) = Stratification d 1
    $ concatMap (\c -> [Ceil3 c 1, Ceil3 c 0]) c

-- подразбиение на более мелкие уровни
celling :: Ceil3 -> [Ceil3]
celling (Ceil3 c i) = [Ceil3 c (i*2), Ceil3 c (i*2-1)]

-- Тут что-то не то.
shiftCeil :: Diameter -> Diameter -> F Point -> Ceil3 -> [Ceil3]
shiftCeil d1 d2 fp (Ceil3 c l) = myNub $ do
    let Point x y = fromCeil d1 c
    func d1 d2 fp
        <$> mySequence x d1
        <*> mySequence y d1
        <*> mySequence (rad d2 l) (pi/2/d2)

func :: Double -> Double -> F Point -> Double -> Double -> Double -> Ceil3
func d1 d2 fp x y t = Ceil3 (ceil d1 fp x y) (int d2 (arcOf fp x y t))

arcOf :: F Point -> Double -> Double -> Double -> Double
arcOf fp x y t = arc $ toPoint $ jacobian fp x y * e t

ceil :: Double -> F Point -> Double -> Double -> Ceil
ceil d fp x y = toCeil d $ fp (Point x y)

jacobian :: F Point -> Double -> Double -> Matrix Double
jacobian fp x y = jac fp $ Point x y

arc :: Point -> Double
arc (Point x y) = if
    | x >  0    ->       phase $          x :+ y
    | x == 0    -> abs $ phase $          x :+ y
    | otherwise ->       phase $ negate $ x :+ y

int :: Double -> Double -> Int
int d r = ceiling $ r/pi*2*d

rad :: Double -> Int ->Double
rad d r = toEnum r * pi/2/d

max_stars_repo_name: vojiranto/dsa

max_issues_repo_name: vojiranto/dsa

max_forks_repo_name: vojiranto/dsa


text:
-- | Functionality for graphing 1-dimensional vectors.
module System.Console.Ansigraph.Internal.Horizontal (
    displayRV
  , displayCV
  , displayPV
  , renderPV
  , renderRV
  , renderCV
) where

import System.Console.Ansigraph.Internal.Core

import Data.Complex
import Control.Monad.IO.Class (MonadIO)
-- for GHC <= 7.8
import Control.Applicative

---- Graphing Infrastructure  ----

barChars = "█▇▆▅▄▃▂▁ "

-- These values delineate regions for rounding to the nearest 1/8.
barVals :: [Double]
barVals = (+ 0.0625) . (/8) <$> [7,6..0]
     -- = [15/16, 13/16, 11/16, 9/16, 7/16, 5/16, 3/16, 1/16]

{- forward and reverse versions of unicode bar selection
   for positive and negative graph regions respectively -}

bars, barsR :: [(Double,Char)]
bars  = zip barVals barChars

barsR = zip barVals (reverse barChars)

selectBar, selectBarR :: Double -> Char
selectBar x = let l = filter (\p -> fst p < x) bars in
  case l of
       []     -> ' '
       (p:_) -> snd p

selectBarR x = let l = filter (\p -> fst p < x) barsR in
  case l of
       []     -> '█'
       (p:_) -> snd p

-- | Simple vector to String rendering that assumes positive input. Yields String of Unicode chars
--   representing graph bars varying in units of 1/8. The IO 'display' functions are preferable
--   for most use cases.
renderPV :: [Double] -> String
renderPV xs = let mx = maximum (filter (>= 0) $ 0:xs) in
              (selectBar . (/mx)) <$> xs

-- | Simple real vector rendering as a pair of strings. The IO 'display' functions are
--   preferable for most use cases.
renderRV :: [Double] -> (String,String)
renderRV l = let rp = l
                 rm = negate <$> rp
                 mx = maximum $ rp ++ rm
  in (selectBar  . (/mx) <$> rp,
      selectBarR . (/mx) <$> rm)

-- | Simple complex vector rendering as a pair of strings. The IO 'display' functions are
--   preferable for most use cases.
renderCV :: [Complex Double] -> (String,String,String,String)
renderCV l = let rp = realPart <$> l
                 rm = negate   <$> rp
                 ip = imagPart <$> l
                 im = negate   <$> ip
                 mx = maximum $ rp ++ rm ++ ip ++ im
  in (selectBar  . (/mx) <$> rp,
      selectBarR . (/mx) <$> rm,
      selectBar  . (/mx) <$> ip,
      selectBarR . (/mx) <$> im)

-- | ANSI based display for positive real vectors. Primarily invoked via 'graph', 'graphWith',
--   'animate', 'animateWith'.
displayPV :: MonadIO m => GraphSettings -> [Double] -> m ()
displayPV s l = let (rp,_) = renderRV l
                    rcol   = realColors s in colorStrLn rcol rp

-- | ANSI based display for real vectors. Primarily invoked via 'graph', 'graphWith',
--   'animate', 'animateWith'.
displayRV :: MonadIO m => GraphSettings -> [Double] -> m ()
displayRV s l = let (rp,rm) = renderRV l
                    rcol    = realColors s
  in do colorStrLn rcol          rp
        colorStrLn (invert rcol) rm

-- | ANSI based display for complex vectors. Primarily invoked via 'graph', 'graphWith',
--   'animate', 'animateWith'.
displayCV :: MonadIO m => GraphSettings -> [Complex Double] -> m ()
displayCV s l = let (rp,rm,ip,im) = renderCV l
                    (rcol,icol)   = colorSets s
  in do colorStrLn rcol          rp
        colorStrLn (invert rcol) rm
        colorStrLn icol          ip
        colorStrLn (invert icol) im

max_stars_repo_name: fieldstrength/ansigraph

max_issues_repo_name: fieldstrength/ansigraph

max_forks_repo_name: fieldstrength/ansigraph


text:
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE DataKinds #-}
{-# Language TypeOperators, TypeApplications, ScopedTypeVariables #-}

module ConCat.Circ where

import Prelude hiding ((.), id, curry, uncurry, const)
import ConCat.Circuit ( CircuitM, (:>)
  , Bus(..),Comp(..),Input,Output, Ty(..), busTy, Source(..), Template(..)
  , GenBuses(..), GS, genBusesRep', tyRep
  -- , delayCRep, unDelayName
  , namedC, constC -- , constS
  , genUnflattenB'
  , SourceToBuses(..), CompS(..), simpleComp
  , mkGraph
  , Attr
  , graphDot, writeDot, displayDot,Name,Report,Graph
  -- , simpleComp
  , tagged
  , systemSuccess
  -- For AbsTy
  , BusesM, abstB,abstC,reprC,Buses(..)
  , OkCAR
  -- , Ty(..)
  )
import GHC.Generics (unComp1, (:.:)(..), Generic)

import ConCat.Misc
import qualified ConCat.RunCircuit as RC
import ConCat.Rep
import ConCat.Category

import ConCat.AltCat ()
import ConCat.Rebox () -- necessary for reboxing rules to fire
import ConCat.Orphans ()
import ConCat.Free.VectorSpace
import ConCat.Free.Affine
import ConCat.Free.LinearRow
import ConCat.Additive

import Data.Vector.Sized (Vector)
import Data.Constraint
import Data.Complex
import Data.Key

-- This is what a morphism looks like in LCirc
data RLC a = Res a
  | Cap a
  | Ind a
  deriving (Eq, Show, Functor, Generic)

instance (HasRep a) => HasRep (RLC a) where
  type Rep (RLC a) = (Maybe a, Maybe a, Maybe a)
  abst (Just r, Nothing, Nothing) = Res r
  abst (Nothing, Just c, Nothing) = Cap c
  abst (Nothing, Nothing, Just i) = Ind i
  repr (Res r) = (Just r, Nothing, Nothing)
  repr (Cap c) = (Nothing, Just c, Nothing)
  repr (Ind i) = (Nothing, Nothing, Just i)

instance (HasV R a) => HasV R (Maybe a) where
  type V R (Maybe a) = Maybe :.: V R a
  toV = Comp1 . fmap toV
  unV = fmap unV . unComp1

instance (HasV R a) => HasV R (Complex a)

instance (HasRep a, HasV R a) => HasV R (RLC a)

instance (GenBuses a) => GenBuses (RLC a) where
  genBuses' = genBusesRep'
  ty = tyRep @(RLC a)
  unflattenB' = genUnflattenB'

newtype VI a = VI (a, a) deriving (Eq, Ord, Generic)

instance HasRep a => HasRep (VI a) where
  type Rep (VI a) = (a, a)
  repr (VI (v, i)) = (v, i)
  abst = VI

instance (GenBuses a) => GenBuses (VI a) where
  genBuses' = genBusesRep'
  ty = tyRep @(VI a)
  unflattenB' = genUnflattenB'

instance (HasRep a, HasV R a) => HasV R (VI a)

--elem :: VI' :> RLC
--elem = namedC "elem"

type VI' = VI R

type LC a = (VI a, RLC a)

-- laplace transform all the tings

--class Circ i b where
--  codup :: (Num b) => (LC b :* LC b) -> (LC b) -> (LC b)
--  dup :: (Num b) => (LC b) -> ((LC b) :* (LC b))
--  unit :: (Num b) => () -> (LC b)
--  counit :: (Num b) => (LC b) -> ()

blackbox :: Fractional a => (VI a) -> (RLC a) -> (VI a)
blackbox (VI (v,i)) (Res r) = VI (v + r*i, i)
blackbox (VI (v,i)) (Cap c) = VI (v + (i / c), i)
blackbox (VI (v, i)) (Ind l) = VI (v + (l * i), i)

--bb :: (Fractional b) => Circ a -> VI b
--bb (Codup (lc1, lc2) lc3) = (uncurry blackbox $ lc1)

--blackbox :: (RLC :> RLC) -> Affine R R R
--blackbox = undefined

newtype Blackbox a = Blackbox (L R a R) deriving (Generic)

instance (HasRep a) => HasRep (Blackbox a) where
  type Rep (Blackbox a) = a
  abst = undefined $ Blackbox
  repr (Blackbox a) = undefined $ a

instance (GenBuses a) => GenBuses (Blackbox a) where
  genBuses' = genBusesRep'
  ty = tyRep @(Blackbox a)
  unflattenB' = genUnflattenB'

instance OkFunctor (:>) (Blackbox) where
  okFunctor = Entail (Sub Dict)

max_stars_repo_name: eco-energy/LCirc

max_issues_repo_name: eco-energy/LCirc

max_forks_repo_name: eco-energy/LCirc


text:
module Main where

import qualified Statistics.Matrix as M
import qualified Statistics.Matrix.Fast as F
import qualified Statistics.Matrix.Fast.Algorithms as FA
import           Statistics.Matrix (Matrix (..))
import qualified Statistics.Matrix.Algorithms as A

import qualified Data.Vector.Unboxed         as U
import           Data.Vector.Unboxed         (Vector)

import qualified System.Random.MWC as Mwc

import qualified Weigh         as W

n :: Int
n = 100

testVector :: IO (Vector Double)
testVector = do
    gen <- Mwc.create
    Mwc.uniformVector gen (n*n)

testMatrix :: IO Matrix
testMatrix = do
    vec <- testVector
    return $ Matrix n n vec

weight :: Vector Double -> Matrix -> Matrix -> IO ()
weight v a b = do
    let
      v2 = U.take n v
    W.mainWith (do
        W.func "norm" M.norm v2
        W.func "Fast.norm" F.norm v2

        W.func "multiplyV" (M.multiplyV a) (v2)
        W.func "Fast.multiplyV" (F.multiplyV a) (v2)
        W.func "transpose" M.transpose a
        W.func "Fast.transpose" F.transpose a
        W.func "ident" M.ident n
        W.func "diag" M.diag v2

        W.func "multiply" (M.multiply a) b
        W.func "Fast.multiply" (F.multiply a) b

        W.func "qr" A.qr a
        W.func "Fast.qr" FA.qr a

        )

main :: IO ()
main = do
    v <- testVector
    a <- testMatrix
    b <- testMatrix

    putStrLn "---Benchmarking memory consumption---"
    weight v a b

max_stars_repo_name: arvindd/dh-core

max_issues_repo_name: arvindd/dh-core

max_forks_repo_name: arvindd/dh-core


text:
-- Modulo destinado a todas as funcoes relativas a execucao da Rede.

module Execution
(execute,
 feedforward) where

import InputOutput
import Types
import Numeric.LinearAlgebra.HMatrix

type Image = [Double]
type Sample = (Int, Image)

-- execute, por retornar IO, pode interagir chamando outras IO Actions como
-- leitura dos weights, biases e input da Rede, por exemplo.
execute :: IO String -- mero esqueleto da funcao de execucao
execute = do
        network <- readIn
        let image = [1.0 .. 784.0]
            executedNetwork = feedforward image network
        do definitiveAnswer executedNetwork

-- recebe a imagem, a network e computa os calculos,
-- retornando a nova data com os valores de ativacao
-- e zeta do hidden e output alterados.
feedforward :: Image -> Data -> Data
feedforward image network = let zH = add ((wHidden network) #> (fromList image)) (bHidden network)
                                aH = fromList $ sigV (toList zH)
                                zO = add ((wOutput network) #> aH) (bOutput network)
                                aO = fromList $ sigV (toList zO)
                            in Data (wHidden network) (bHidden network) aH zH (wOutput network) (bOutput network) aO zO


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE LambdaCase, TupleSections, TypeOperators #-}
{-# LANGUAGE RankNTypes, MultiParamTypeClasses, FlexibleInstances, FlexibleContexts, ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies, FunctionalDependencies, LiberalTypeSynonyms #-}
{-# LANGUAGE AllowAmbiguousTypes #-}

-- | Predefined references for commonly used data structures and reference generators.
--
-- When defining lenses one should use the more general types. For instance 'Lens' instead of the more strict 'Lens''. This way references with different @m1@ and @m2@ monads can be combined if there is a monad @m'@ for @MMorph m1 m'@ and @MMorph m2 m'@.
module Control.Reference.Predefined where

import Control.Reference.Representation
import Control.Reference.Types
import Control.Reference.Generators
import Control.Reference.Combinators
import Control.Reference.Operators

import Control.Instances.Morph
import Control.Monad
import Data.Ratio
import qualified Data.Text as Text
import Data.Complex
import Control.Monad.State
import Control.Concurrent.MVar hiding (modifyMVarMasked_)
import Control.Concurrent.Chan
import Data.IORef
import Data.Either.Combinators
import Data.STRef
import System.Directory
import System.FilePath
import System.IO
import System.IO.Error
import qualified Control.Exception as Ex

-- * Trivial references

-- | An identical lens. Accesses the context.
--
-- > self & a = a & self = a
self :: Lens a b a b
self = reference return (const . return) id

-- | An empty reference that do not traverse anything
--
-- > emptyRef &+& a = a &+& emptyRef = a
--
-- > a & emptyRef = emptyRef & a = emptyRef
emptyRef :: Simple RefPlus s a
emptyRef = reference (const mzero) (const return) (const return)

-- * References for simple data structures

-- | An indexed lens for accessing points a function
atArg :: Eq a => a -> Simple Lens (a -> b) b
atArg a = lens ($ a) (\b f -> \x -> if a == x then b else f x)

-- | A partial lens to access the value that may not exist
just :: Prism (Maybe a) (Maybe b) a b
just = prism Just Just (maybe (Left Nothing) Right) id

-- | A partial lens to access the right option of an 'Either'
right :: Prism (Either a b) (Either a c) b c
right = prism Right Right (either (Left . Left) Right) rightToMaybe

-- | A partial lens to access the left option of an 'Either'
left :: Prism (Either a c) (Either b c) a b
left = prism Left Left (either Right (Left . Right)) leftToMaybe

-- | Access the value that is in the left or right state of an 'Either'
anyway :: Lens (Either a a) (Either b b) a b
anyway = reference (either return return)
                   (\b -> return . mapBoth (const b) (const b))
                   (\f -> either (f >=> return . Left) (f >=> return . Right))

-- | References both elements of a tuple
both :: Traversal (a,a) (b,b) a b
both = reference (\(x,y)    -> morph [x,y])
                 (\v        -> return . const (v,v))
                 (\f (x,y)  -> (,) <$> f x <*> f y)

-- | References the head of a list
atHead :: Simple Lens [a] (Maybe a)
atHead = lens (\case [] -> Nothing; x:_ -> Just x)
              (\case Nothing -> drop 1;
                     Just v  -> \case []    -> [v]
                                      _:xs  -> v:xs)

-- | References the element at the head of the list
headElem :: Simple Partial [a] a
headElem = atHead & just

-- | References the tail of a list
_tail :: Simple Partial [a] [a]
_tail = simplePartial (\case [] -> Nothing; x:xs -> Just (xs,(x:)))

-- | References a suffix of a list
dropped :: Int -> Simple Partial [a] [a]
dropped 0 = self
dropped i = _tail & dropped (i-1)

-- | Views a list as an optinal pair
view :: Iso [a] [b] (Maybe (a,[a])) (Maybe (b,[b]))
view = iso' to to from from
  where to :: [x] -> Maybe (x,[x])
        to [] = Nothing
        to (x:xs) = Just (x,xs)
        from :: Maybe (x,[x]) -> [x]
        from Nothing = []
        from (Just (x,xs)) = x:xs

-- | An isomorphism between the list and text representation of a string
text :: Simple Iso String Text.Text
text = iso Text.pack Text.unpack

-- | Accesses the reversed version of a list
--
-- > 'turn' reversed == reversed
reversed :: Iso [a] [b] [a] [b]
reversed = iso' reverse reverse reverse reverse

-- | Accesses the numerator of a ratio
_numerator :: Integral a => Simple Lens (Ratio a) a
_numerator = lens numerator (\num' r -> num' % denominator r)

-- | Accesses the denominator of a ratio
_denominator :: Integral a => Simple Lens (Ratio a) a
_denominator = lens denominator (\denom' r -> numerator r % denom')

-- | Accesses the real part of a complex number
_realPart :: RealFloat a => Simple Lens (Complex a) a
_realPart = lens realPart (\real' c -> real' :+ imagPart c)

-- | Accesses the imaginary part of a complex number
_imagPart :: RealFloat a => Simple Lens (Complex a) a
_imagPart = lens imagPart (\imag' c -> realPart c :+ imag')

-- | Accesses the polar representation of a complex number
_polar :: RealFloat a => Simple Lens (Complex a) (a, a)
_polar = iso polar (uncurry mkPolar)

-- * Stateful references

-- | A dummy object to interact with the user through the console.
data Console = Console

-- | Interacts with a line of text on the console. Values set are printed, getting
-- is reading from the console.
consoleLine :: Simple IOLens Console String
consoleLine
  = reference (const (morph getLine))
              (\str -> const (morph (putStrLn str) >> return Console))
              (\f -> const (morph getLine >>= f
                                           >>= morph . putStrLn
                                           >> return Console))

-- | Reference to the contents of the file. Not thread-safe.
--
-- An empty file's content is @Just ""@ while a non-existent file's is @Nothing@
--
-- Creates a temporary file to store the result.
fileContent :: Simple IOLens FilePath (Maybe String)
fileContent
  = reference (\fp -> morph (getFileCont fp))
              (\cont fp -> morph (setFileCont fp cont) >> return fp)
              (\trf fp -> morph (getFile fp) >>= \hcnt -> trf (fmap snd hcnt)
                             >>= morph . updateFileCont fp (fmap fst hcnt) >> return fp)
  where getFileCont :: FilePath -> IO (Maybe String)
        getFileCont fp = (Just <$> readFile fp)
                           `Ex.catch` \e -> if isDoesNotExistError e then return Nothing
                                                                     else Ex.throw e
        getFile :: FilePath -> IO (Maybe (Handle, String))
        getFile fp = do h <- (Just <$> openFile fp ReadMode)
                               `Ex.catch` \e -> if isDoesNotExistError e then return Nothing
                                                                         else Ex.throw e
                        case h of
                          Just handle -> do cont <- hGetContents handle
                                              `Ex.catch` \e -> hClose handle >> Ex.throw (e :: Ex.SomeException)
                                            return $ Just (handle, cont)
                          Nothing -> return Nothing

        setFileCont :: FilePath -> Maybe String -> IO ()
        setFileCont fp Nothing = removeFile fp
        setFileCont fp (Just cont) = writeFile fp cont

        updateFileCont :: FilePath -> Maybe Handle -> Maybe String -> IO ()
        updateFileCont fp h Nothing = (just !| hClose) h >> removeFile fp
        updateFileCont fp h (Just cont)
          = Ex.bracket (openTempFile (takeDirectory fp) (takeFileName fp))
                       (\(tfp,th) -> hClose th >> (just !| hClose) h >> removeFile tfp)
                       (\(_,th) -> hPutStr th cont >> (just !| hClose) h >> hSeek th AbsoluteSeek 0
                                      >> hGetContents th >>= writeFile fp)

-- | Access a value inside an MVar.
-- Setting is not atomic. If there is two supplier that may set the accessed
-- value, one may block and can corrupt the following updates.
--
-- Reads and updates are done in sequence, always using consistent data.
mvar :: Simple IOLens (MVar a) a
mvar = rawReference
         (\f mv -> pullBack $ withMVarMasked mv (sink . f))
         (\newVal mv -> morph $ do empty <- isEmptyMVar mv
                                   if empty then putMVar mv newVal
                                            else swapMVar mv newVal >> return ()
                                   return mv)
         (\trf mv -> modifyMVarMasked_ mv trf >> return mv)
         unusableOp unusableOp unusableOp

-- | Generalized version of 'Control.Concurrent.MVar.modifyMVarMasked_'.
modifyMVarMasked_ :: (Monad m, Morph IO m, MorphControl IO m) => MVar a -> (a -> m a) -> m ()
modifyMVarMasked_ m io =
  mask_ $ do
    a  <- morph (takeMVar m)
    a' <- io a `onException` morph (putMVar m a)
    morph (putMVar m a')

-- | Generalized version of 'Control.Exception.mask_'.
mask_ :: (MorphControl IO m) => m a -> m a
mask_ = pullBack . Ex.mask_ . sink

-- | Generalized version of 'Control.Exception.onException'.
onException :: (MorphControl IO m) => m a -> m b -> m a
onException a b = pullBack $ Ex.onException (sink a) (sink b)

chan :: Simple IOLens (Chan a) a
chan = reference (morph . readChan)
                 (\a ch -> morph (writeChan ch a) >> return ch)
                 (\trf ch -> morph (readChan ch) >>= trf
                               >>= morph . writeChan ch >> return ch)

-- | Access the value of an IORef.
ioref :: Simple IOLens (IORef a) a
ioref = reference (morph . readIORef)
                  (\v ior -> morph (atomicWriteIORef ior v) >> return ior)
                  (\trf ior -> morph (readIORef ior)
                                 >>= trf >>= morph . writeIORef ior >> return ior)

-- | Access the state inside a state monad (from any context).
state :: forall s m a . Monad m => Simple (StateLens s m) a s
state = reference (morph . const get') (\a s -> morph (put' a) >> return s)
                  (\trf s -> (morph get' >>= trf >> return s))
  where put' = put :: s -> StateT s m ()
        get' = get :: StateT s m s

-- | Access the value inside an 'STRef'
stRef :: Simple (STLens s) (STRef s a) a
stRef = reference (morph . readSTRef)
                  (\newVal ref -> morph $ writeSTRef ref newVal >> return ref)
                  (\trf ref -> morph (readSTRef ref) >>= trf
                                 >>= morph . writeSTRef ref >> return ref)

-- | Filters an indexed reference based on the index
whereOf :: (RefMonads w r, MonadPlus r)
        => (i -> Bool) -> (IndexedReference i w r MU MU s s a a) -> (IndexedReference i w r MU MU s s a a)
whereOf p iref i | p i       = iref i
                 | otherwise = emptyRef

max_stars_repo_name: mingmingrr/references

max_issues_repo_name: mingmingrr/references

max_forks_repo_name: mingmingrr/references


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

import           Control.Monad                  as M
import qualified Data.Array.Accelerate          as A
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.Repa                as R
import           Data.Complex
import           Data.List                      as L
import           Data.Vector.Storable           as VS
import           Data.Vector.Unboxed            as VU
import           DFT.Plan
import           Filter.Utils
import           Foreign.CUDA.Driver            as CUDA
import           FourierMethod.BlockCudaMatrix
import           FourierMethod.FourierSeries2D
import           Image.IO
import           Pinwheel.Base
import           Pinwheel.FourierSeries2D
import           System.Directory
import           System.Environment
import           System.FilePath
import           System.Random
import           Utils.SimpsonRule
import Utils.Array
import Numeric.GSL.Special.Bessel
import           Utils.Distribution

main = do
  args@(deviceIDsStr:numPointsStr:numR2FreqStr:deltaStr:deltaReconStr:periodR2Str:angularFreqStr:radialFreqStr:sigmaStr:periodEnvelopeStr:stdStr:numBatchStr:_) <-
    getArgs
  let deviceIDs = read deviceIDsStr :: [Int]
      numPoints = read numPointsStr :: Int
      numR2Freqs = read numR2FreqStr :: Int
      delta = read deltaStr :: Double
      deltaRecon = read deltaReconStr :: Double
      periodR2 = read periodR2Str :: Double
      angularFreq = read angularFreqStr :: Int
      radialFreq = read radialFreqStr :: Int
      sigma = read sigmaStr :: Double
      periodEnvelope = read periodEnvelopeStr :: Double
      std = read stdStr :: Double
      numBatch = read numBatchStr :: Int
      folderPath = "output/test/DFTPinwheel"
  removePathForcibly folderPath
  createDirectoryIfMissing True folderPath
  initialise []
  devs <- M.mapM device deviceIDs
  ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
  ptxs <- M.mapM createTargetFromContext ctxs
  let center = div numPoints 2
      pinwheel =
        computeUnboxedS . R.fromFunction (Z :. numPoints :. numPoints) $ \(Z :. x :. y) ->
          fourierMellinInvPeriod
            sigma
            (periodEnvelope * sqrt 2)
            angularFreq
            radialFreq
            ( deltaRecon * fromIntegral (x - center)
            , deltaRecon * fromIntegral (y - center))
      pinwheelMat =
        CuMat (numPoints ^ 2) 1 . CuVecHost . VU.convert . toUnboxed $ pinwheel
      centerFreq = div numR2Freqs 2
      -- pinwheelFreq =
      --   computeUnboxedS . R.fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. x :. y) ->
      --     fourierMellin
      --       sigma
      --       angularFreq
      --       radialFreq
      --       ( delta * fromIntegral (x - centerFreq)
      --       , delta * fromIntegral (y - centerFreq))
      -- pinwheelFreqMat =
      --   CuMat 1 (numR2Freqs ^ 2) . CuVecHost . VU.convert . toUnboxed $ -- . computeS . weightedArray
      --   pinwheelFreq
  -- invHarmonics <-
  --   createInverseHarmonicMatriesGPU ptxs 1 numPoints numR2Freqs periodR2 delta
  -- pinwheelCoef <-
  --   computeFourierCoefficientsR2Stream
  --     deviceIDs
  --     ptxs
  --     numR2Freqs
  --     numPoints
  --     periodR2
  --     delta
  --     1
  --     numBatch
  --     [pinwheelMat]
  -- pinwheelFreqSeries <-
  --   computeFourierSeriesR2
  --     deviceIDs
  --     numR2Freqs
  --     numPoints
  --     periodR2
  --     invHarmonics
  --     [pinwheelFreqMat]
  plotImageRepaComplex (folderPath </> "Pinwheel.png") .
    ImageRepa 8 .
    computeS . extend (Z :. (1 :: Int) :. All :. All) . R.traverse pinwheel id $ \f idx@(Z :. i :. j) ->
    let x = fromIntegral $ i - center
        y = fromIntegral $ j - center
        r = sqrt $ x ^ 2 + y ^ 2
       -- if r <= 6
       --   then 0
       --   else
     in f idx
  -- plotImageRepaComplex (folderPath </> "PinwheelTest.png") .
  --   ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --   pinwheelTest
  -- plotImageRepaComplex (folderPath </> "PinwheelCoefficients.png") .
  --   ImageRepa 8 . computeS . R.traverse pinwheelCoef id $ \f idx@(Z :. _ :. i :. j) ->
  --   let x = fromIntegral $ i - div numR2Freqs 2
  --       y = fromIntegral $ j - div numR2Freqs 2
  --       r = sqrt $ x ^ 2 + y ^ 2
  --   in if r <= 0
  --        then 0
  --        else f idx
  -- plotImageRepaComplex (folderPath </> "PinwheelFrequency.png") .
  --   ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --   pinwheelFreq
  -- plotImageRepaComplex (folderPath </> "PinwheelFrequencySeries.png") .
  --   ImageRepa 8 $
  --   pinwheelFreqSeries
  initVec <-
    (VS.fromList . L.map (\x -> x :+ 0)) <$>
    M.replicateM (numPoints ^ 2) randomIO
  lock <- getFFTWLock
  plan <-
    fst <$> idft1dGPlan lock emptyPlan [numPoints, numPoints] [0, 1] initVec
  -- dftPinwheel <-
  --   dftExecute plan (DFTPlanID IDFT1DG [numPoints, numPoints] [0, 1]) .
  --   VU.convert . toUnboxed . computeUnboxedS . makeFilter2D $
  --   pinwheel
  -- plotImageRepaComplex (folderPath </> "DFTPinwheel.png") .
  --   ImageRepa 8 .
  --   computeUnboxedS .
  --   makeFilter2D .
  --   fromUnboxed (Z :. (1 :: Int) :. numPoints :. numPoints) . VS.convert $
  --   dftPinwheel
  let centerR2Freq = div numR2Freqs 2
      gaussian2D =
        computeUnboxedS . fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i' :. j') ->
          let i = i' - centerR2Freq
              j = j' - centerR2Freq
              rho = (sqrt . fromIntegral $ i ^ 2 + j ^ 2) / periodR2 * 2*pi
           in (-- rho ^ (abs angularFreq) *
               exp
                 (pi * fromIntegral (i ^ 2 + j ^ 2) /
                  ((-1) * periodR2 ^ 2 * std ^ 2)) /
               (2 * pi * std ^ 2)) :+
              0
      -- pinwheelFreqAnaticalNonGaussian =
      --   centerHollowArray numR2Freqs . computeUnboxedS $
      --   analyticalFourierCoefficients2
      --     numR2Freqs
      --     delta
      --     angularFreq
      --     radialFreq
      --     (-sigma) --(sigma - 1)
      --     periodR2
      --     (periodEnvelope * sqrt 2)
      pinwheelFreqAnatical =
        centerHollowArray numR2Freqs .
        computeUnboxedS
        -- . R.zipWith (*) gaussian2D
        $
        analyticalFourierCoefficients2
          numR2Freqs
          delta
          angularFreq
          radialFreq
          (-sigma) --(sigma - 1)
          periodR2
          (periodEnvelope * sqrt 2)
      -- pinwheelFreqAnaticalMat =
      --   A.use .
      --   A.fromList (A.Z A.:. (numR2Freqs ^ 2) A.:. (1 :: Int)) . R.toList $
      --   pinwheelFreqAnatical
      -- pinwheelFreqAnaticalInner =
      --   centerHollowArray numR2Freqs . computeUnboxedS $
      --   analyticalFourierCoefficients2'
      --     numR2Freqs
      --     (1 / fromIntegral (div numR2Freqs 2))
      --     angularFreq
      --     radialFreq
      --     (-sigma)
      --     periodR2
      --     (periodEnvelope * sqrt 2)
      -- pinwheelFreqAnaticalInnerMat =
      --   A.use .
      --   A.fromList (A.Z A.:. (numR2Freqs ^ 2) A.:. (1 :: Int)) . R.toList $
      --   pinwheelFreqAnaticalInner
      -- pinwheelFreqSeriesAnatical =
      --   computeUnboxedS $
      --   analyticalFourierSeries1
      --     numPoints
      --     1 -- delta
      --     angularFreq
      --     radialFreq
      --     (sigma - 1)
      --     periodR2
  plotImageRepaComplex (folderPath </> "PinwheelCoefficientsAnalytical.png") .
    ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
    pinwheelFreqAnatical
  -- plotImageRepaComplex
  --   (folderPath </> "PinwheelCoefficientsAnalyticalNonGaussian.png") .
  --   ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --   pinwheelFreqAnaticalNonGaussian
  -- plotImageRepaComplex
  --   (folderPath </> "PinwheelCoefficientsAnalyticalInner.png") .
  --   ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --   pinwheelFreqAnaticalInner
  -- plotImageRepaComplex (folderPath </> "PinwheelFrequencySeriesAnalytical.png") .
  --   ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --   pinwheelFreqSeriesAnatical
  -- pinwheelFreqSeriesStream <-
  --   computeFourierSeriesR2Stream
  --     deviceIDs
  --     ptxs
  --     numR2Freqs
  --     numPoints
  --     periodR2
  --     delta
  --     numBatch
  --     [transposeCuMat pinwheelFreqMat]
  -- plotImageRepaComplex (folderPath </> "PinwheelFrequencySeriesStream.png") .
  --   ImageRepa 8 . computeS . R.traverse pinwheelFreqSeriesStream id $ \f idx@(Z :. _ :. i :. j) ->
  --   let x = fromIntegral $ i - div numPoints 2
  --       y = fromIntegral $ j - div numPoints 2
  --       r = sqrt $ x ^ 2 + y ^ 2
  --   in if r <= 0
  --        then 0
  --        else f idx
  -- pinwheelCoefficientsAnalyticalSeriesStream <-
  --   computeFourierSeriesR2StreamAcc
  --     ptxs
  --     numR2Freqs
  --     numPoints
  --     1
  --     periodR2
  --     deltaRecon
  --     numBatch
  --     pinwheelFreqAnaticalMat
  --     -- [ transposeCuMat . createCuMat $
  --     --   [VU.convert . toUnboxed $ pinwheelFreqAnatical]
  --     -- ]
  -- plotImageRepaComplex
  --   (folderPath </> "PinwheelCoefficientsAnalyticalSeries.png") .
  --   ImageRepa 8 .
  --   computeS .
  --   extend (Z :. (1 :: Int) :. All :. All) .
  --   R.traverse (sumS . rotate3D $ pinwheelCoefficientsAnalyticalSeriesStream) id $ \f idx@(Z :. i :. j) ->
  --   let x = fromIntegral $ i - center
  --       y = fromIntegral $ j - center
  --       r = sqrt $ x ^ 2 + y ^ 2
  --      -- if r <= 6
  --      --   then 0
  --      --   else
  --    in f idx
  pinwheelCoefficientsAnalyticalIDFT <-
    dftExecute plan (DFTPlanID IDFT1DG [numPoints, numPoints] [0, 1]) .
    VU.convert . toUnboxed . computeS . makeFilter2D $
    pinwheelFreqAnatical
  plotImageRepaComplex (folderPath </> "PinwheelCozefficientsIDFT.png") .
    ImageRepa 8 .
    computeS .
    makeFilter2DInverse .
    fromUnboxed (Z :. (1 :: Int) :. numPoints :. numPoints) . VS.convert $
    pinwheelCoefficientsAnalyticalIDFT
  -- pinwheelCoefficientsAnalyticalInnerSeriesStream <-
  --   computeFourierSeriesR2StreamAcc'
  --     ptxs
  --     numR2Freqs
  --     numPoints
  --     1
  --     periodR2
  --     1 -- (1 / fromIntegral (div numR2Freqs 2))
  --     numBatch
  --     pinwheelFreqAnaticalMat
  -- plotImageRepaComplex
  --   (folderPath </> "PinwheelCoefficientsAnalyticalInnverSeries.png") .
  --   ImageRepa 8 .
  --   computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
  --   pinwheelCoefficientsAnalyticalInnerSeriesStream
  -- plotImageRepaComplex
  --   (folderPath </> "PinwheelCoefficientsAnalyticalSumSeries.png") .
  --   ImageRepa 8 .
  --   computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
  --   R.zipWith
  --     (+)
  --     pinwheelCoefficientsAnalyticalSeriesStream
  --     pinwheelCoefficientsAnalyticalInnerSeriesStream
  --   computeS . R.traverse pinwheelCoefficientsAnalyticalSeriesStream id $ \f idx@(Z :. _ :. i :. j) ->
  --   let x = fromIntegral $ i - div numPoints 2
  --       y = fromIntegral $ j - div numPoints 2
  --       r = sqrt $ x ^ 2 + y ^ 2
  --   in if r <= 0
  --        then 0
  --        else f idx
  -- -- pinwheelCoefficientsSeriesStream <-
  --   computeFourierSeriesR2Stream
  --     deviceIDs
  --     ptxs
  --     numR2Freqs
  --     numPoints
  --     periodR2
  --     delta
  --     numBatch
  --     [transposeCuMat . createCuMat $ [VU.convert . toUnboxed $ pinwheelCoef]]
  -- plotImageRepaComplex (folderPath </> "PinwheelCoefficientsSeriesStream.png") .
  --   ImageRepa 8 . computeS . R.traverse pinwheelCoefficientsSeriesStream id $ \f idx@(Z :. _ :. i :. j) ->
  --   let x = fromIntegral $ i - div numPoints 2
  --       y = fromIntegral $ j - div numPoints 2
  --       r = sqrt $ x ^ 2 + y ^ 2
  --   in if r <= 0
  --        then 0
  --        else f idx
  -- let besselFilter =
  --       computeUnboxedS . fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i :. j) ->
  --         let xFreq = fromIntegral $ i - centerFreq
  --             yFreq = fromIntegral $ j - centerFreq
  --             rho = sqrt $ xFreq ^ 2 + yFreq ^ 2
  --         in if rho == 0
  --              then 1
  --              else (bessel_J1 (2 * pi * rho / periodR2)) / rho :+ 0
  -- plotImageRepaComplex (folderPath </> "BesselFilter.png") .
  --   ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --   besselFilter
  -- initVec1 <-
  --   (VS.fromList . L.map (\x -> x :+ 0)) <$>
  --   M.replicateM (numR2Freqs ^ 2) randomIO
  -- lock <- getFFTWLock
  -- plan <-
  --   fst <$>
  --   (dft1dGPlan lock emptyPlan [numR2Freqs, numR2Freqs] [0, 1] initVec1 >>= \(plan, vec) ->
  --      idft1dGPlan lock plan [numR2Freqs, numR2Freqs] [0, 1] vec)
  -- let besselVec =
  --       VU.convert . toUnboxed . computeS . makeFilter2D $ besselFilter
  --     pinwheelFreqAnaticalVec = VU.convert . toUnboxed $ pinwheelFreqAnatical
  -- besselVecF <-
  --   dftExecute plan (DFTPlanID DFT1DG [numR2Freqs, numR2Freqs] [0, 1]) besselVec
  -- pinwheelFreqAnaticalVecF <-
  --   dftExecute
  --     plan
  --     (DFTPlanID DFT1DG [numR2Freqs, numR2Freqs] [0, 1])
  --     pinwheelFreqAnaticalVec
  -- pinwheelFreqAnatical1 <-
  --   fmap (fromUnboxed (Z :. numR2Freqs :. numR2Freqs) . VS.convert) .
  --   dftExecute plan (DFTPlanID IDFT1DG [numR2Freqs, numR2Freqs] [0, 1]) $
  --   VS.zipWith (*) besselVecF pinwheelFreqAnaticalVecF
  -- plotImageRepaComplex
  --   (folderPath </> "BesselFilterConveledPinwheelFreqAnatical.png") .
  --   ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --   pinwheelFreqAnatical1
  -- let pinwheelFreqAnatical2 =
  --       R.zipWith (-) pinwheelFreqAnatical pinwheelFreqAnatical1
  --     pinwheelFreqAnaticalMat2 =
  --       A.use .
  --       A.fromList (A.Z A.:. (numR2Freqs ^ 2) A.:. (1 :: Int)) . R.toList $
  --       pinwheelFreqAnatical2
  -- plotImageRepaComplex (folderPath </> "pinwheelFreqAnatical2.png") .
  --   ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --   pinwheelFreqAnatical2
  -- pinwheelCoefficientsAnalyticalSeriesStream2 <-
  --   computeFourierSeriesR2StreamAcc
  --     ptxs
  --     numR2Freqs
  --     numPoints
  --     1
  --     periodR2
  --     delta
  --     numBatch
  --     pinwheelFreqAnaticalMat2
  -- plotImageRepaComplex
  --   (folderPath </> "PinwheelCoefficientsAnalyticalSeries2.png") .
  --   ImageRepa 8 .
  --   computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
  --   pinwheelCoefficientsAnalyticalSeriesStream2
  -- dftPinwheel <-
  --   fmap (fromUnboxed (Z :. numR2Freqs :. numR2Freqs) . VS.convert) .
  --   dftExecute plan (DFTPlanID DFT1DG [numR2Freqs, numR2Freqs] [0, 1]) .
  --   VU.convert . toUnboxed . computeS . makeFilter2D $
  --   pinwheel
  -- let zeroArr =
  --       fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i :. j) ->
  --         if (sqrt . fromIntegral $ (i - centerFreq) ^ 2 + (j - centerFreq) ^ 2) <
  --            1
  --           then 0
  --           else 1
  -- plotImageRepaComplex (folderPath </> "DFTPinwheel.png") .
  --   ImageRepa 8 .
  --   computeS .
  --   extend (Z :. (1 :: Int) :. All :. All) .
  --   R.zipWith (*) zeroArr . makeFilter2D $
  --   dftPinwheel
  -- let normFunc arr =
  --       let m = VU.maximum . toUnboxed . computeS . R.map magnitude $ arr
  --       in R.map (/ (m :+ 0)) arr
  -- plotImageRepaComplex (folderPath </> "Diff.png") .
  --   ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --   R.zipWith
  --     (\a b ->
  --        if magnitude b == 0
  --          then 0
  --          else a / b)
  --     (normFunc . R.zipWith (*) zeroArr . makeFilter2D $ dftPinwheel)
  --     (normFunc . R.zipWith (*) zeroArr $ pinwheelFreqAnatical)

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:
module TransitionMatrix
( makeMatrix
, randomMatrix
, randomMatrixNullDiag
, mkRandomTransitionMatrix
, mkRandomNullDiagTransitionMatrix
, mkBlockTransitionMatrix
, toTransitionMatrix
, toNLMatrix
, fromNLMatrix
, eig
, estimatePAbs
, TransitionMatrix (..)
, TransitionMatrixGenerator (..)
) where

import Data.Matrix
import qualified Data.Vector as V
import System.Random
import Control.Monad (forM)
import qualified Numeric.LinearAlgebra as NL

import MC

newtype TransitionMatrix = TransitionMatrix { toMatrix :: Matrix Double }
                           deriving (Show, Eq)

data TransitionMatrixGenerator =
    RandomMG Int Double |
    RandomNullDiagMG Int Double |
    BlockMG Int Double Double deriving (Eq, Ord, Show)

-- | Produce a transition matrix
makeMatrix :: TransitionMatrixGenerator -- ^ a generator
              -> MC TransitionMatrix    -- ^ a TransitionMatrix, in the MC monad
makeMatrix (RandomMG size absorptionProbability) = mkRandomTransitionMatrix size absorptionProbability
makeMatrix (RandomNullDiagMG size absorptionProbability) = mkRandomNullDiagTransitionMatrix size absorptionProbability
makeMatrix (BlockMG size absorptionProbability coupling) = mkBlockTransitionMatrix size absorptionProbability coupling

{- |
   Produce a random transition matrix
-}
mkRandomTransitionMatrix :: Int           -- ^ the matrix size
                         -> Double        -- ^ the absorption probability
                         -> MC TransitionMatrix
mkRandomTransitionMatrix size pAbs = do
    m <- randomMatrix size
    return $ toTransitionMatrix pAbs m

{- |
   Produce a random transition matrix with null elements on the diagonal
-}
mkRandomNullDiagTransitionMatrix :: Int           -- ^ the matrix size
                                 -> Double        -- ^ the absorption probability
                                 -> MC TransitionMatrix
mkRandomNullDiagTransitionMatrix size pAbs = do
    m <- randomMatrixNullDiag size
    return $ toTransitionMatrix pAbs m

{- |
   Produce a block transition matrix
-}
mkBlockTransitionMatrix :: Int           -- ^ the matrix size
                        -> Double        -- ^ the absorption probability
                        -> Double        -- ^ the coupling
                        -> MC TransitionMatrix
mkBlockTransitionMatrix size pAbs coupling = do
    let blockSize1 = size `div` 2
    let blockSize2 = size - blockSize1
    m1 <- randomMatrix blockSize1
    m4 <- randomMatrix blockSize2
    let fromCoupling = blockSize1 `div` 2
    let toCoupling = blockSize2 `div` 2
    let m2 = zero blockSize2 blockSize1
    let m3 = matrix blockSize1 blockSize2 (\(i,j) -> if j==fromCoupling && i==toCoupling then coupling else 0.0)
    let m = joinBlocks (m1, m2, m3, m4)
    return $ toTransitionMatrix pAbs m

-- | Return a random square matrix of given size, with elements in the [0,1[
--   range
randomMatrix :: Int                 -- ^ size of the matrix
             -> MC (Matrix Double)  -- ^ the matrix, in the MC monad because it's random
randomMatrix n = do
    lss <- forM [0..n-1] $ \ _ -> randomList n
    return $ fromLists lss

-- | Return a random square matrix of given size, with elements in the [0,1[
--   range and null diagonal elements
randomMatrixNullDiag :: Int                 -- ^ size of the matrix
                     -> MC (Matrix Double)  -- ^ the matrix, in the MC monad because it's random
randomMatrixNullDiag n = do
    lss <- forM [0..n-1] $ \ k -> do
        ls <- randomList (n-1)
        let (left, right) = splitAt k ls
        let ls' = left ++ [0.0] ++ right
        return ls'
    return $ fromLists lss

-- | Generate a random list of size n
randomList :: (Random a, Fractional a)
           => Int   -- ^ the list size
           -> MC [a]
randomList 0 = return []
randomList n = do
    x <- uniform
    xs <- randomList (n-1)
    return (x:xs)

--normalizeRow :: Double -> Int -> Matrix Double -> Matrix Double
--normalizeRow norm i m = mapRow (\ _ x -> x*norm/tot) i m
--    where tot = sum $ getRow i m

normalizeCol :: Double -> Int -> Matrix Double -> Matrix Double
normalizeCol norm i m = mapCol (\ _ x -> x*norm/tot) i m
    where tot = sum $ getCol i m

{- |
   Rescale the elements of the given matrix by multiplying by the
   absorption probability
-}
toTransitionMatrix :: Double            -- ^ the absorption probability
                   -> Matrix Double     -- ^ the input matrix
                   -> TransitionMatrix
toTransitionMatrix pAbs m = let nc   = ncols m
                                norm = 1.0 - pAbs
                             in TransitionMatrix $ foldr (normalizeCol norm) m [1..nc]

-- functions for conversion to/from Numeric.LinearAlgebra

toNLMatrix :: NL.Element a => Matrix a -> NL.Matrix a
toNLMatrix = NL.fromLists . toLists

fromNLMatrix :: NL.Element a => NL.Matrix a -> Matrix a
fromNLMatrix = fromLists . NL.toLists

eig :: (NL.Element a, NL.Field a) => Matrix a -> (V.Vector (NL.Complex Double), Matrix (NL.Complex Double))
eig m = let (vnl, mnl) = NL.eig $ toNLMatrix m
            v' = V.fromList $ NL.toList vnl
            m' = fromNLMatrix mnl
         in (v', m')

avg :: (Num a, Fractional a, Foldable t) => t a -> a
avg v = sum v / fromIntegral (length v)

estimatePAbs :: TransitionMatrix -> Double
estimatePAbs (TransitionMatrix m) = 1.0 - fromIntegral (nrows m) * avg (toList m)

max_stars_repo_name: arekfu/markov-chains

max_issues_repo_name: arekfu/markov-chains

max_forks_repo_name: arekfu/markov-chains


text:
#!/usr/bin/env runhaskell
{-# LANGUAGE OverloadedStrings, NoImplicitPrelude #-}
import ClassyPrelude hiding ((<>)) -- this <> is for semigroups (usually better)
import Data.Monoid ((<>)) -- ... but HMatrix only has Monoid
import Numeric.LinearAlgebra (tr) -- transpose
import qualified Numeric.LinearAlgebra as HMatrix
import qualified Data.Text.Format as Format

main = do
  -- Load the term-document matrix
  putStrLn "Loading"
  matrix <- HMatrix.loadMatrix "termdoc.txt"
  (u, sigma, vt) <- stochasticTruncatedSVD 50 2 matrix

  print $ HMatrix.size u
  print $ HMatrix.size sigma
  print $ HMatrix.size vt

-- Stochastic SVD. See Halko, Martinsson, and Tropp, 2010 for an explanation
stochasticTruncatedSVD :: Int -> Int -> HMatrix.Matrix Double -> IO (HMatrix.Matrix Double, HMatrix.Vector Double, HMatrix.Matrix Double)
stochasticTruncatedSVD top_vectors num_iterations original = do
  let (m, n) = HMatrix.size original
  let k = top_vectors + 10
  -- Stage A
  putStrLn "Stage A\n\tPower Iteration"
  omega <- HMatrix.randn n k
  let
    y 0 = original <> omega
    y n = seq (y (n-1)) $ original <> tr original <> y (n-1)
  bigY <- return $! y num_iterations -- Force evaluation here
  putStrLn "\tOrthogonalizing"
  q <- return $! HMatrix.orth bigY

  -- Stage B
  putStrLn "Stage B"
  let b = tr q <> original
      (uhat, sigma, vt) = HMatrix.thinSVD b
      u = q <> uhat

  return $! (HMatrix.takeColumns top_vectors u,
    HMatrix.subVector 0 top_vectors sigma,
    HMatrix.takeRows top_vectors $ tr vt)
  --putStrLn "Checking"
  --let ahat = u <> HMatrix.diag sigma <> tr vt

  --Format.print "MSE of approximation {}\n" [HMatrix.sumElements ((a-ahat)^2) / (fromIntegral $ m*n)]

max_stars_repo_name: SeanTater/albemarle

max_issues_repo_name: SeanTater/albemarle

max_forks_repo_name: SeanTater/albemarle