text:
{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE TypeFamilies #-}

-----------------------------------------------------------------------------
{- |
Module      :  Numeric.LinearAlgebra.Algorithms
Copyright   :  (c) Alberto Ruiz 2006-9
License     :  GPL-style

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional
Portability :  uses ffi

High level generic interface to common matrix computations.

Specific functions for particular base types can also be explicitly
imported from "Numeric.LinearAlgebra.LAPACK".

-}
-----------------------------------------------------------------------------

module Numeric.LinearAlgebra.Algorithms (
-- * Supported types
    Field(),
-- * Linear Systems
    linearSolve,
    luSolve,
    cholSolve,
    linearSolveLS,
    linearSolveSVD,
    inv, pinv,
    det, invlndet,
    rank, rcond,
-- * Matrix factorizations
-- ** Singular value decomposition
    svd,
    fullSVD,
    thinSVD,
    compactSVD,
    singularValues,
    leftSV, rightSV,
-- ** Eigensystems
    eig, eigSH, eigSH',
    eigenvalues, eigenvaluesSH, eigenvaluesSH',
    geigSH',
-- ** QR
    qr, rq,
-- ** Cholesky
    chol, cholSH, mbCholSH,
-- ** Hessenberg
    hess,
-- ** Schur
    schur,
-- ** LU
    lu, luPacked,
-- * Matrix functions
    expm,
    sqrtm,
    matFunc,
-- * Nullspace
    nullspacePrec,
    nullVector,
    nullspaceSVD,
    orth,
-- * Norms
    Normed(..), NormType(..),
    relativeError,
-- * Misc
    eps, peps, i,
-- * Util
    haussholder,
    unpackQR, unpackHess,
    pinvTol,
    ranksv
) where

import Data.Packed.Internal hiding ((//))
import Data.Packed.Matrix
import Numeric.LinearAlgebra.LAPACK as LAPACK
import Data.List(foldl1')
import Data.Array
import Numeric.ContainerBoot

{- | Class used to define generic linear algebra computations for both real and complex matrices. Only double precision is supported in this version (we can
transform single precision objects using 'single' and 'double').

-}
class (Product t,
       Convert t,
       Container Vector t,
       Container Matrix t,
       Normed Matrix t,
       Normed Vector t) => Field t where
    svd'         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
    thinSVD'     :: Matrix t -> (Matrix t, Vector Double, Matrix t)
    sv'          :: Matrix t -> Vector Double
    luPacked'    :: Matrix t -> (Matrix t, [Int])
    luSolve'     :: (Matrix t, [Int]) -> Matrix t -> Matrix t
    linearSolve' :: Matrix t -> Matrix t -> Matrix t
    cholSolve'   :: Matrix t -> Matrix t -> Matrix t
    linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
    linearSolveLS'  :: Matrix t -> Matrix t -> Matrix t
    eig'         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
    eigSH''      :: Matrix t -> (Vector Double, Matrix t)
    eigOnly      :: Matrix t -> Vector (Complex Double)
    eigOnlySH    :: Matrix t -> Vector Double
    cholSH'      :: Matrix t -> Matrix t
    mbCholSH'    :: Matrix t -> Maybe (Matrix t)
    qr'          :: Matrix t -> (Matrix t, Matrix t)
    hess'        :: Matrix t -> (Matrix t, Matrix t)
    schur'       :: Matrix t -> (Matrix t, Matrix t)

instance Field Double where
    svd' = svdRd
    thinSVD' = thinSVDRd
    sv' = svR
    luPacked' = luR
    luSolve' (l_u,perm) = lusR l_u perm
    linearSolve' = linearSolveR                 -- (luSolve . luPacked) ??
    cholSolve' = cholSolveR
    linearSolveLS' = linearSolveLSR
    linearSolveSVD' = linearSolveSVDR Nothing
    eig' = eigR
    eigSH'' = eigS
    eigOnly = eigOnlyR
    eigOnlySH = eigOnlyS
    cholSH' = cholS
    mbCholSH' = mbCholS
    qr' = unpackQR . qrR
    hess' = unpackHess hessR
    schur' = schurR

instance Field (Complex Double) where
#ifdef NOZGESDD
    svd' = svdC
    thinSVD' = thinSVDC
#else
    svd' = svdCd
    thinSVD' = thinSVDCd
#endif
    sv' = svC
    luPacked' = luC
    luSolve' (l_u,perm) = lusC l_u perm
    linearSolve' = linearSolveC
    cholSolve' = cholSolveC
    linearSolveLS' = linearSolveLSC
    linearSolveSVD' = linearSolveSVDC Nothing
    eig' = eigC
    eigOnly = eigOnlyC
    eigSH'' = eigH
    eigOnlySH = eigOnlyH
    cholSH' = cholH
    mbCholSH' = mbCholH
    qr' = unpackQR . qrC
    hess' = unpackHess hessC
    schur' = schurC

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

square m = rows m == cols m

vertical m = rows m >= cols m

exactHermitian m = m `equal` ctrans m

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

-- | Full singular value decomposition.
svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
svd = {-# SCC "svd" #-} svd'

-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
--
-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@.
thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
thinSVD = {-# SCC "thinSVD" #-} thinSVD'

-- | Singular values only.
singularValues :: Field t => Matrix t -> Vector Double
singularValues = {-# SCC "singularValues" #-} sv'

-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
--
-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
fullSVD m = (u,d,v) where
    (u,s,v) = svd m
    d = diagRect 0 s r c
    r = rows m
    c = cols m

-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
compactSVD :: Field t  => Matrix t -> (Matrix t, Vector Double, Matrix t)
compactSVD m = (u', subVector 0 d s, v') where
    (u,s,v) = thinSVD m
    d = rankSVD (1*eps) m s `max` 1
    u' = takeColumns d u
    v' = takeColumns d v

-- | Singular values and all right singular vectors.
rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
          | otherwise  = let (_,s,v) = svd m     in (s,v)

-- | Singular values and all left singular vectors.
leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
leftSV m  | vertical m = let (u,s,_) = svd m     in (u,s)
          | otherwise  = let (u,s,_) = thinSVD m in (u,s)

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

-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
luPacked :: Field t => Matrix t -> (Matrix t, [Int])
luPacked = {-# SCC "luPacked" #-} luPacked'

-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
luSolve = {-# SCC "luSolve" #-} luSolve'

-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
linearSolve = {-# SCC "linearSolve" #-} linearSolve'

-- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'.
cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t
cholSolve = {-# SCC "cholSolve" #-} cholSolve'

-- | Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular valures which are less than 'eps' times the largest singular value.
linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD'

-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.
linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS'

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

-- | Eigenvalues and eigenvectors of a general square matrix.
--
-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
eig = {-# SCC "eig" #-} eig'

-- | Eigenvalues of a general square matrix.
eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
eigenvalues = {-# SCC "eigenvalues" #-} eigOnly

-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t)
eigSH' = {-# SCC "eigSH'" #-} eigSH''

-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
eigenvaluesSH' :: Field t => Matrix t -> Vector Double
eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH

-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.
--
-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
eigSH m | exactHermitian m = eigSH' m
        | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"

-- | Eigenvalues of a complex hermitian or real symmetric matrix.
eigenvaluesSH :: Field t => Matrix t -> Vector Double
eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m
                | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"

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

-- | QR factorization.
--
-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
qr :: Field t => Matrix t -> (Matrix t, Matrix t)
qr = {-# SCC "qr" #-} qr'

-- | RQ factorization.
--
-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
rq :: Field t => Matrix t -> (Matrix t, Matrix t)
rq m =  {-# SCC "rq" #-} (r,q) where
    (q',r') = qr $ trans $ rev1 m
    r = rev2 (trans r')
    q = rev2 (trans q')
    rev1 = flipud . fliprl
    rev2 = fliprl . flipud

-- | Hessenberg factorization.
--
-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary
-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).
hess        :: Field t => Matrix t -> (Matrix t, Matrix t)
hess = hess'

-- | Schur factorization.
--
-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary
-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
-- upper triangular in 2x2 blocks.
--
-- \"Anything that the Jordan decomposition can do, the Schur decomposition
-- can do better!\" (Van Loan)
schur       :: Field t => Matrix t -> (Matrix t, Matrix t)
schur = schur'

-- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
mbCholSH :: Field t => Matrix t -> Maybe (Matrix t)
mbCholSH = {-# SCC "mbCholSH" #-} mbCholSH'

-- | Similar to 'chol', without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
cholSH      :: Field t => Matrix t -> Matrix t
cholSH = {-# SCC "cholSH" #-} cholSH'

-- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
--
-- If @c = chol m@ then @c@ is upper triangular and @m == ctrans c \<> c@.
chol :: Field t => Matrix t ->  Matrix t
chol m | exactHermitian m = cholSH m
       | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"

-- | Joint computation of inverse and logarithm of determinant of a square matrix.
invlndet :: (Floating t, Field t)
         => Matrix t
         -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det))
invlndet m | square m = (im,(ladm,sdm))
           | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix"
  where
    lp@(lup,perm) = luPacked m
    s = signlp (rows m) perm
    dg = toList $ takeDiag $ lup
    ladm = sum $ map (log.abs) dg
    sdm = s* product (map signum dg)
    im = luSolve lp (ident (rows m))

-- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'.
det :: Field t => Matrix t -> t
det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup)
      | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix"
    where (lup,perm) = luPacked m
          s = signlp (rows m) perm

-- | Explicit LU factorization of a general matrix.
--
-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
lu = luFact . luPacked

-- | Inverse of a square matrix. See also 'invlndet'.
inv :: Field t => Matrix t -> Matrix t
inv m | square m = m `linearSolve` ident (rows m)
      | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix"

-- | Pseudoinverse of a general matrix.
pinv :: Field t => Matrix t -> Matrix t
pinv m = linearSolveSVD m (ident (rows m))

-- | Numeric rank of a matrix from the SVD decomposition.
rankSVD :: Element t
        => Double   -- ^ numeric zero (e.g. 1*'eps')
        -> Matrix t -- ^ input matrix m
        -> Vector Double -- ^ 'sv' of m
        -> Int      -- ^ rank of m
rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s)

-- | Numeric rank of a matrix from its singular values.
ranksv ::  Double   -- ^ numeric zero (e.g. 1*'eps')
        -> Int      -- ^ maximum dimension of the matrix
        -> [Double] -- ^ singular values
        -> Int      -- ^ rank of m
ranksv teps maxdim s = k where
    g = maximum s
    tol = fromIntegral maxdim * g * teps
    s' = filter (>tol) s
    k = if g > teps then length s' else 0

-- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave).
eps :: Double
eps =  2.22044604925031e-16

-- | 1 + 0.5*peps == 1,  1 + 0.6*peps /= 1
peps :: RealFloat x => x
peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x)

-- | The imaginary unit: @i = 0.0 :+ 1.0@
i :: Complex Double
i = 0:+1

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

-- | The nullspace of a matrix from its SVD decomposition.
nullspaceSVD :: Field t
             => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
                                  --   or Right \"theoretical\" matrix rank.
             -> Matrix t          -- ^ input matrix m
             -> (Vector Double, Matrix t) -- ^ 'rightSV' of m
             -> [Vector t]        -- ^ list of unitary vectors spanning the nullspace
nullspaceSVD hint a (s,v) = vs where
    tol = case hint of
        Left t -> t
        _      -> eps
    k = case hint of
        Right t -> t
        _       -> rankSVD tol a s
    vs = drop k $ toRows $ ctrans v

-- | The nullspace of a matrix. See also 'nullspaceSVD'.
nullspacePrec :: Field t
              => Double     -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps')
              -> Matrix t   -- ^ input matrix
              -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m)

-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision.
nullVector :: Field t => Matrix t -> Vector t
nullVector = last . nullspacePrec 1

orth :: Field t => Matrix t -> [Vector t]
-- ^ Return an orthonormal basis of the range space of a matrix
orth m = take r $ toColumns u
  where
    (u,s,_) = compactSVD m
    r = ranksv eps (max (rows m) (cols m)) (toList s)

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

{-  Pseudoinverse of a real matrix with the desired tolerance, expressed as a
multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').

@\> let m = 'fromLists' [[1,0,    0]
                    ,[0,1,    0]
                    ,[0,0,1e-10]]
\  --
\> 'pinv' m
1. 0.           0.
0. 1.           0.
0. 0. 10000000000.
\  --
\> pinvTol 1E8 m
1. 0. 0.
0. 1. 0.
0. 0. 1.@

-}
--pinvTol :: Double -> Matrix Double -> Matrix Double
pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
    (u,s,v) = thinSVDRd m
    sl@(g:_) = toList s
    s' = fromList . map rec $ sl
    rec x = if x < g*tol then 1 else 1/x
    tol = (fromIntegral (max r c) * g * t * eps)
    r = rows m
    c = cols m
    d = dim s
    u' = takeColumns d u
    v' = takeColumns d v

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

-- many thanks, quickcheck!

haussholder :: (Field a) => a -> Vector a -> Matrix a
haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w))
    where w = asColumn v

zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs)
              where xs = toList v

zt 0 v = v
zt k v = join [subVector 0 (dim v - k) v, konst 0 k]

unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
unpackQR (pq, tau) =  {-# SCC "unpackQR" #-} (q,r)
    where cs = toColumns pq
          m = rows pq
          n = cols pq
          mn = min m n
          r = fromColumns $ zipWith zt ([m-1, m-2 .. 1] ++ repeat 0) cs
          vs = zipWith zh [1..mn] cs
          hs = zipWith haussholder (toList tau) vs
          q = foldl1' mXm hs

unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
unpackHess hf m
    | rows m == 1 = ((1><1)[1],m)
    | otherwise = (uH . hf) m

uH (pq, tau) = (p,h)
    where cs = toColumns pq
          m = rows pq
          n = cols pq
          mn = min m n
          h = fromColumns $ zipWith zt ([m-2, m-3 .. 1] ++ repeat 0) cs
          vs = zipWith zh [2..mn] cs
          hs = zipWith haussholder (toList tau) vs
          p = foldl1' mXm hs

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

-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
rcond :: Field t => Matrix t -> Double
rcond m = last s / head s
    where s = toList (singularValues m)

-- | Number of linearly independent rows or columns.
rank :: Field t => Matrix t -> Int
rank m = rankSVD eps m (singularValues m)

{-
expm' m = case diagonalize (complex m) of
    Just (l,v) -> v `mXm` diag (exp l) `mXm` inv v
    Nothing -> error "Sorry, expm not yet implemented for non-diagonalizable matrices"
  where exp = vectorMapC Exp
-}

diagonalize m = if rank v == n
                    then Just (l,v)
                    else Nothing
    where n = rows m
          (l,v) = if exactHermitian m
                    then let (l',v') = eigSH m in (real l', v')
                    else eig m

-- | Generic matrix functions for diagonalizable matrices. For instance:
--
-- @logm = matFunc log@
--
matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
matFunc f m = case diagonalize m of
    Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v
    Nothing -> error "Sorry, matFunc requires a diagonalizable matrix"

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

golubeps :: Integer -> Integer -> Double
golubeps p q = a * fromIntegral b / fromIntegral c where
    a = 2^^(3-p-q)
    b = fact p * fact q
    c = fact (p+q) * fact (p+q+1)
    fact n = product [1..n]

epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]]

geps delta = head [ k | (k,g) <- epslist, g<delta]

{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
     based on a scaled Pade approximation.
-}
expm :: Field t => Matrix t -> Matrix t
expm = expGolub

expGolub :: ( Fractional t, Element t, Field t
            , Normed Matrix t
            , RealFrac (RealOf t)
            , Floating (RealOf t)
            ) => Matrix t -> Matrix t
expGolub m = iterate msq f !! j
    where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m
          a = m */ fromIntegral ((2::Int)^j)
          q = geps eps -- 7 steps
          eye = ident (rows m)
          work (k,c,x,n,d) = (k',c',x',n',d')
              where k' = k+1
                    c' = c * fromIntegral (q-k+1) / fromIntegral ((2*q-k+1)*k)
                    x' = a <> x
                    n' = n |+| (c' .* x')
                    d' = d |+| (((-1)^k * c') .* x')
          (_,_,_,nf,df) = iterate work (1,1,eye,eye,eye) !! q
          f = linearSolve df nf
          msq x = x <> x

          (<>) = multiply
          v */ x = scale (recip x) v
          (.*) = scale
          (|+|) = add

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

{- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia.
It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try @matFunc sqrt@.

@m = (2><2) [4,9
           ,0,4] :: Matrix Double@

@\>sqrtm m
(2><2)
 [ 2.0, 2.25
 , 0.0,  2.0 ]@
-}
sqrtm ::  Field t => Matrix t -> Matrix t
sqrtm = sqrtmInv

sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x))
    where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps   = a
                                | otherwise = fixedPoint (b:rest)
          fixedPoint _ = error "fixedpoint with impossible inputs"
          f (y,z) = (0.5 .* (y |+| inv z),
                     0.5 .* (inv y |+| z))
          (.*) = scale
          (|+|) = add
          (|-|) = sub

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

signlp r vals = foldl f 1 (zip [0..r-1] vals)
    where f s (a,b) | a /= b    = -s
                    | otherwise =  s

swap (arr,s) (a,b) | a /= b    = (arr // [(a, arr!b),(b,arr!a)],-s)
                   | otherwise = (arr,s)

fixPerm r vals = (fromColumns $ elems res, sign)
    where v = [0..r-1]
          s = toColumns (ident r)
          (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals)

triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]]
    where el p q = if q-p>=h then v else 1 - v

luFact (l_u,perm) | r <= c    = (l ,u ,p, s)
                  | otherwise = (l',u',p, s)
  where
    r = rows l_u
    c = cols l_u
    tu = triang r c 0 1
    tl = triang r c 0 0
    l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst 1 r) r r
    u = l_u |*| tu
    (p,s) = fixPerm r perm
    l' = (l_u |*| tl) |+| diagRect 0 (konst 1 c) r c
    u' = takeRows c (l_u |*| tu)
    (|+|) = add
    (|*|) = mul

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

data NormType = Infinity | PNorm1 | PNorm2 | Frobenius

class (RealFloat (RealOf t)) => Normed c t where
    pnorm :: NormType -> c t -> RealOf t

instance Normed Vector Double where
    pnorm PNorm1    = norm1
    pnorm PNorm2    = norm2
    pnorm Infinity  = normInf
    pnorm Frobenius = norm2

instance Normed Vector (Complex Double) where
    pnorm PNorm1    = norm1
    pnorm PNorm2    = norm2
    pnorm Infinity  = normInf
    pnorm Frobenius = pnorm PNorm2

instance Normed Vector Float where
    pnorm PNorm1    = norm1
    pnorm PNorm2    = norm2
    pnorm Infinity  = normInf
    pnorm Frobenius = pnorm PNorm2

instance Normed Vector (Complex Float) where
    pnorm PNorm1    = norm1
    pnorm PNorm2    = norm2
    pnorm Infinity  = normInf
    pnorm Frobenius = pnorm PNorm2

instance Normed Matrix Double where
    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns
    pnorm PNorm2    = (@>0) . singularValues
    pnorm Infinity  = pnorm PNorm1 . trans
    pnorm Frobenius = pnorm PNorm2 . flatten

instance Normed Matrix (Complex Double) where
    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns
    pnorm PNorm2    = (@>0) . singularValues
    pnorm Infinity  = pnorm PNorm1 . trans
    pnorm Frobenius = pnorm PNorm2 . flatten

instance Normed Matrix Float where
    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns
    pnorm PNorm2    = realToFrac . (@>0) . singularValues . double
    pnorm Infinity  = pnorm PNorm1 . trans
    pnorm Frobenius = pnorm PNorm2 . flatten

instance Normed Matrix (Complex Float) where
    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns
    pnorm PNorm2    = realToFrac . (@>0) . singularValues . double
    pnorm Infinity  = pnorm PNorm1 . trans
    pnorm Frobenius = pnorm PNorm2 . flatten

-- | Approximate number of common digits in the maximum element.
relativeError :: (Normed c t, Container c t) => c t -> c t -> Int
relativeError x y = dig (norm (x `sub` y) / norm x)
    where norm = pnorm Infinity
          dig r = round $ -logBase 10 (realToFrac r :: Double)

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

-- | Generalized symmetric positive definite eigensystem Av = lBv,
-- for A and B symmetric, B positive definite (conditions not checked).
geigSH' :: Field t
        => Matrix t -- ^ A
        -> Matrix t -- ^ B
        -> (Vector Double, Matrix t)
geigSH' a b = (l,v')
  where
    u = cholSH b
    iu = inv u
    c = ctrans iu <> a <> iu
    (l,v) = eigSH' c
    v' = iu <> v
    (<>) = mXm

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE StrictData #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE TypeFamilies #-}

module Main (main) where

import qualified Data.Vector as V
import           Data.Vector.Unboxed.Deriving (derivingUnbox)
import qualified Data.Vector.Generic as VG
import qualified Data.Vector.Generic.Mutable as VGM
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Unboxed as VU
import qualified Numeric.LinearAlgebra as HMat

newtype Vec = Vec { unV :: HMat.Vector Double }
  deriving (Eq, Ord, Show)

-- Does crash with GHC 8.6.5

dims :: Int
dims = 3

data instance VU.MVector s Vec = MVector_Vec (VU.MVector s Double)
data instance VU.Vector    Vec = Vector_Vec  (VU.Vector    Double)

instance VGM.MVector VU.MVector Vec where
  basicLength (MVector_Vec v) = VGM.basicLength v `div` dims
  basicUnsafeSlice start size (MVector_Vec v) = MVector_Vec (VGM.basicUnsafeSlice (start * dims) (size * dims) v)
  basicOverlaps (MVector_Vec v1) (MVector_Vec v2) = VGM.basicOverlaps v1 v2
  basicUnsafeNew size = MVector_Vec <$> VGM.basicUnsafeNew (size * dims)
  basicInitialize (MVector_Vec vec) = VGM.basicInitialize vec
  basicUnsafeRead vec idx = fmap (Vec . VG.convert) $ VU.freeze sliceVals
    where MVector_Vec sliceVals = VGM.unsafeSlice idx 1 vec
  basicUnsafeWrite vec idx (Vec val) = VG.unsafeCopy sliceVals (VG.convert val)
    where MVector_Vec sliceVals = VGM.unsafeSlice idx 1 vec
  basicUnsafeCopy (MVector_Vec target) (MVector_Vec source) = VGM.basicUnsafeCopy target source
  basicUnsafeMove (MVector_Vec target) (MVector_Vec source) = VGM.basicUnsafeMove target source

instance VG.Vector VU.Vector Vec where
  basicUnsafeFreeze (MVector_Vec mvec) = Vector_Vec <$> VG.basicUnsafeFreeze mvec
  basicUnsafeThaw (Vector_Vec vec) = MVector_Vec <$> VG.basicUnsafeThaw vec
  basicLength (Vector_Vec vec) = VG.basicLength vec `div` dims
  basicUnsafeSlice idx size (Vector_Vec vec) = Vector_Vec $ VG.basicUnsafeSlice (idx * dims) (size * dims) vec
  basicUnsafeIndexM (Vector_Vec vec) idx = return $! Vec $ VS.convert $! VG.unsafeSlice (idx * dims) dims vec
  basicUnsafeCopy (MVector_Vec mvec) (Vector_Vec vec) = VG.unsafeCopy mvec vec
  elemseq _ x y = x `seq` y

instance VU.Unbox Vec

-- Does not crash:

-- $(derivingUnbox "Vec"
--   [t| Vec -> (Double, Double, Double) |]
--   [| \(Vec hm) -> let [a, b, c] = HMat.toList hm in (a, b, c) |]
--   [| \(a, b, c) -> Vec (HMat.fromList [a, b, c]) |]
--  )

newtype NormalVec = NormalVec { getNormalVec :: Vec }
  deriving (Eq, Ord, Show)

$(derivingUnbox "NormalVec"
  [t| NormalVec -> Vec |]
  [| getNormalVec |]
  [| NormalVec |]
 )

unboxedVectorVecRoundtrip :: V.Vector Vec -> Bool
unboxedVectorVecRoundtrip v =
  let
    !indexVector = VU.fromList [0,1,2]
    !v' = VG.map (v VG.!) indexVector
  in
    VG.convert v' == V.map (v VG.!) (VG.convert indexVector)

unboxedVectorNormalVecRoundtrip :: V.Vector NormalVec -> Bool
unboxedVectorNormalVecRoundtrip v =
  let
    !indexVector = VU.fromList [0,1,2]
    !v' = VG.map (v VG.!) indexVector
  in
    VG.convert v' == V.map (v VG.!) (VG.convert indexVector)

main :: IO ()
main = do
  putStrLn "Unboxed.Vector (Vec World) roundtrip (does not crash in repro)"

  let v = Vec (HMat.fromList [1,2,3])

  print $ unboxedVectorVecRoundtrip $
    V.fromList [ v, v, v ]

  putStrLn "Unboxed.Vector (NormalVec World) roundtrip (CRASHES in repro)"

  print $ unboxedVectorNormalVecRoundtrip $
    V.fromList [ NormalVec v, NormalVec v, NormalVec v ]

ext: hs


text:
module Parser (renderImage, parseMNISTSet, Image, Label, MNISTSet, MNISTSample) where

-- import Codec.Compression.GZip (decompress)
import qualified Data.ByteString.Lazy as BS
import qualified Data.ByteString as B
import Data.Functor
-- import System.Random
import Data.Word8
import Data.Binary.Get
import Control.Monad
import Control.Monad.Except as E
import qualified Data.Vector.Storable as V
-- import Numeric.LinearAlgebra
-- import Numeric.LinearAlgebra.Devel

data IDXType = IDXUByte | IDXSByte | IDXShort | IDXInt | IDXFloat | IDXLong | IDXWTF deriving (Show)

-- IDX format header :: (datatype, dimensionality, [size of each dimension])
makeIDXType :: Word8 -> Maybe IDXType
makeIDXType 0x08 = Just IDXUByte
makeIDXType 0x09 = Just IDXSByte
makeIDXType 0x0B = Just IDXShort
makeIDXType 0x0C = Just IDXInt
makeIDXType 0x0D = Just IDXFloat
makeIDXType 0x0E = Just IDXLong
makeIDXType _    = Nothing

type IDXParser = E.ExceptT String Get

runIDXParser :: IDXParser a -> BS.ByteString -> Either String a
runIDXParser p bs = runGet (runExceptT p) bs

type Header = (IDXType, Int, [Int])

data IDXData = IDXData Header BS.ByteString

instance Show IDXData where
  show (IDXData (t, d, s:ss) _) = show s ++ " " ++ show d ++ "-dimensional samples of type " ++ show t

parseIDXHeader :: IDXParser Header
parseIDXHeader = do
  (0x00 : 0x00 : t : ds : []) <- lift $ replicateM 4 getWord8
  let ds' = fromIntegral ds
  case makeIDXType t of
    Just t' -> do
      sizes <- lift $ (replicateM ds' getWord32be)
      return (t', ds', fromIntegral <$> sizes)
    Nothing -> throwError "Invalid byte for datatype"

parseIDXData :: IDXParser IDXData
parseIDXData = do
  h <- parseIDXHeader
  bs <- lift getRemainingLazyByteString
  return $ IDXData h bs

type Image = V.Vector Double
type Label = V.Vector Double

type MNISTSample = (Image, Label)
-- MNISTSet = ((nrows, ncols), nsamples, mnistsamples)
type MNISTSet = ((Int, Int), Int, [MNISTSample])

-- maps a Word8 value into a number between 0 and 1
word8ToDouble :: Word8 -> Double
word8ToDouble w = (fromIntegral w) / 255.0

chunksOf :: V.Storable a => Int -> V.Vector a -> [V.Vector a]
chunksOf n v
  | V.length v == 0 = []
  | n < V.length v = let (c, c') = V.splitAt n v
                     in  c:chunksOf n c'
  | otherwise = [v]

-- Parse a sequence of bytes(Word8) into a vector of numbers(Double)
getImage :: (Int, Int) -> Get Image
getImage (nrows, ncols) = do
  bs <- getByteString (nrows*ncols)
  return $ V.generate (nrows*ncols) $ \i -> word8ToDouble (B.index bs i)

getImages :: Int -> (Int,Int) -> Get [Image]
getImages n (nrows, ncols) = do
  replicateM n (getImage (nrows, ncols))

parseMNISTLabels :: IDXParser (Int, [Label])
parseMNISTLabels = do
  idxdata <- parseIDXData
  case idxdata of
    IDXData (IDXUByte, 1, [nlabels]) lbs -> return $ (nlabels, fmap word8ToVector $ BS.unpack lbs)
    IDXData (_, 1, [nlabels]) lbs -> throwError "Wrong data type for mnist labels"
    IDXData (_, _, _) lbs -> throwError "Wrong data dimensionality for mnist labels set"

parseMNISTImages :: IDXParser (Int, (Int,Int), [Image])
parseMNISTImages = do
  idxdata <- parseIDXData
  case idxdata of
    IDXData (IDXUByte, 3, [nimages, nrows, ncols]) ibs ->
      let images = runGet (getImages nimages (nrows,ncols)) ibs
      in return (nimages, (nrows, ncols), images)
    IDXData (_, 3, [nimages, nrows, ncols]) ibs -> throwError "Wrong data type for mnist images"
    IDXData (_, _, _) ibs -> throwError "Wrong data dimensionality for mnist images set"

parseMNISTSet :: BS.ByteString -> BS.ByteString -> Either String MNISTSet
parseMNISTSet ibs lbs = do
  (nlabels, labels) <- runIDXParser parseMNISTLabels lbs
  (nimages, (nrows, ncols), images) <- runIDXParser parseMNISTImages ibs
  if nlabels == nimages
    then
    let samples = zip images labels
        nsamples = nlabels
    in return ((nrows, ncols), nsamples, samples)
    else throwError "Number of labels not equal to number of images"

renderPixel :: Double -> Char
renderPixel n = let s = " .:oO@" in s !! (floor (n*256) * length s `div` 256)

renderImage :: (Int, Int) -> Image -> String
renderImage (nrows, ncols) = unlines . (map V.toList) . (chunksOf ncols) . (V.map renderPixel)

-- converts a Label(a byte representing a digit) to a vector corresponding to the desired ouput of a network
word8ToVector :: Word8 -> V.Vector Double
word8ToVector l = V.generate 10 $ \i -> if i == (fromIntegral l) then 1 else 0

max_stars_repo_name: DrPyser/NeuralNetwork

max_issues_repo_name: DrPyser/NeuralNetwork

max_forks_repo_name: DrPyser/NeuralNetwork


text:
module Main where

import Numeric.LinearAlgebra

x1 = [  3,   4,   5,   6,   7,   8] :: [Double]
y1 = [150, 155, 150, 170, 160, 175] :: [Double]

x2 = [-0.2, 0.2, 1]     :: [Double]
y2 = [0.49, 0.64, 1.39] :: [Double]

x3 = [1,2,3,4]  :: [Double]
y3 = [6,5,7,10] :: [Double]

-- | Least square estimator
--
lse :: Int -> [Double] -> [Double] -> Matrix Double
lse deg x y = pinv (tx' <> x') <> tx' <> y'
    where

        tx' = trans x'

        x' :: Matrix Double
        x' = (lx><deg) (foldr (\x xs -> pow x (deg-1) ++ xs) [] x)

        y' :: Matrix Double
        y' = (ly><1) y

        lx = length x
        ly = length y

        pow x deg = go [] x deg
            where go :: [Double] -> Double -> Int -> [Double]
                  go acc _ 0   = 1:acc
                  go acc x deg = go (x**fromIntegral(deg):acc) x (deg-1)

-- | prettify for gnuplot
--
toF :: Matrix Double -> String
toF = concatMap (\(i,[x]) -> show (round' x) ++ "*" ++ "x**"++ (show i) ++ " + ") .  zip [0..] . toLists
  where round' = (/100) . fromInteger . round . (*100) :: Double -> Double

-- | Mean square estimator
mse :: Matrix Double -> [(Double, Double)] -> Double
mse m xys = nth . sum $ map (\(x,y) -> ((f x) - y)**2) xys
  where l = toLists m :: [[Double]]
        f x = sum $ zipWith (\(i, [w]) x -> w * x^i) (zip [0..] l) (repeat x)
        nth = (/ fromIntegral (length xys))

-- | data -> fitting -> gnuplot
plot :: Int -> [Double] -> [Double] -> String
plot deg x y = toF $ lse deg x y

max_stars_repo_name: cirquit/random-code

max_issues_repo_name: cirquit/Personal-Repository

max_forks_repo_name: cirquit/Personal-Repository


text:
-----------------------------------------------------------------------------
-- |
-- Module      : Berp.Base
-- Copyright   : (c) 2010 Bernie Pope
-- License     : BSD-style
-- Maintainer  : florbitous@gmail.com
-- Stability   : experimental
-- Portability : ghc
--
-- This module exports all the primitive functions which are needed by
-- the compiled programs. Avoid putting extraneous exports in this file
-- because it is imported by all compiled programs.
--
-----------------------------------------------------------------------------

module Berp.Base
   ( module Complex
   , int, none, string, true, false, def, lambda, (=:), stmt, ifThenElse, ret, pass, break
   , continue, while, whileElse, for, forElse, ifThen, (@@), tailCall, tuple, read, var
   , (%), (+), (-), (*), (.), (/), (==), (<), (>), (<=), (>=), and, or, not, klass, setattr, list, dictionary
   , subs, try, tryElse, tryFinally, tryElseFinally, except, exceptDefault, raise, reRaise, raiseFrom
   , pure, pureObject, yield, generator, returnGenerator, unaryMinus, unaryPlus, invert, runEval
   , interpretStmt, topVar, unpack, setitem, Pat (..), complex, set
   , importModule, readGlobal, writeGlobal, readLocal, writeLocal, run, runWithGlobals, importAll )
   where

import Prelude ()
import Data.Complex as Complex (Complex (..))
import Berp.Base.TopLevel (importModule, importAll, run, runWithGlobals)
import Berp.Base.Prims
   ( (=:), stmt, ifThenElse, ret, pass, break, continue, while, whileElse, for, forElse, ifThen,
     (@@), tailCall, read, var, setattr, subs, try, tryElse, tryFinally, tryElseFinally, except,
     exceptDefault, raise, reRaise, raiseFrom, yield, def, lambda, generator, returnGenerator, topVar, pure,
     pureObject, unpack, setitem, Pat (..), readGlobal, writeGlobal,
     readLocal, writeLocal )
import Berp.Base.Operators
   ((%), (+), (-), (*), (.), (/), (==), (<), (>), (<=), (>=), and, or, unaryMinus, unaryPlus, invert, not)
import Berp.Base.Monad (runEval, interpretStmt)
import Berp.Base.Class (klass)
import Berp.Base.StdTypes.Integer (int)
import Berp.Base.StdTypes.Tuple (tuple)
import Berp.Base.StdTypes.Bool (true, false)
import Berp.Base.StdTypes.String (string)
import Berp.Base.StdTypes.None (none)
import Berp.Base.StdTypes.List (list)
import Berp.Base.StdTypes.Dictionary (dictionary)
import Berp.Base.StdTypes.Set (set)
import Berp.Base.StdTypes.Complex (complex)

max_stars_repo_name: ppelleti/berp

max_issues_repo_name: ppelleti/berp

max_forks_repo_name: bjpop/berp


text:
module Main where

import           Control.Monad.Error (runErrorT)
import           Control.Monad.Mersenne.Random (evalRandom)
import           Numeric.IEEE (IEEE, epsilon)
import           System.Random.Mersenne.Pure64 (pureMT)
import qualified Data.Vector.Unboxed as V
import           Statistics.Test.ApproxRand
import           Test.HUnit (Assertion, assertBool, assertEqual)
import           Test.Framework
import           Test.Framework.Providers.HUnit

tests :: Test
tests = testGroup "Paired approximate randomization tests" $
  concat [statTests, randomizationTests]

main :: IO ()
main = defaultMain [ tests ]

-- Statistics tests

statTests :: [Test]
statTests = [meanDifferenceTest]

meanDifferenceTest :: Test
meanDifferenceTest =
  testIEEEEquality "mean difference robot competition"
    1.8 $ meanDifference cohenRobotsAlpha cohenRobotsBeta

-- Approximate andomization tests

randomizationTests :: [Test]
randomizationTests = [pairApproxExactTestScores]

pairApproxExactTestScores :: Test
pairApproxExactTestScores =
  testEquality  "number of extreme values robot competition"
    (Right 21) $ V.length `fmap` V.filter (>= 1.8) `fmap` scores
  where
    test = runErrorT $
      approxRandPairStats differenceMean 1024 cohenRobotsAlpha cohenRobotsBeta
    scores = evalRandom test $ pureMT 42

-- Helper functions

testEquality :: (Show a, Eq a) => String -> a -> a -> Test
testEquality msg a b = testCase msg $ assertEqual msg a b

testIEEEEquality :: IEEE a => String -> a -> a -> Test
testIEEEEquality msg a b = testCase msg $ assertEqualIEEE msg a b

assertEqualIEEE :: IEEE a => String -> a -> a -> Assertion
assertEqualIEEE msg a b = assertBool msg $ fracEq epsilon a b

-- Suggested in The Floating Point Guide: http://floating-point-gui.de/
fracEq :: (Fractional a, Ord a) => a -> a -> a -> Bool
fracEq eps a b
  | a == b     = True
  | a * b == 0 = diff < (eps * eps)
  | otherwise  = diff / (abs a + abs b) < eps
  where diff = abs (a - b)

-- Example from Cohen, 1995
cohenRobotsAlpha :: V.Vector Double
cohenRobotsAlpha = V.fromList [8,3,9,6,5,8,7,8,9,9]

cohenRobotsBeta :: V.Vector Double
cohenRobotsBeta  = V.fromList [7,0,9,4,5,9,8,3,4,5]

max_stars_repo_name: danieldk/approx-rand-test

max_issues_repo_name: danieldk/approx-rand-test

max_forks_repo_name: danieldk/approx-rand-test


text:
module Main
( main
) where

import Debug.Trace

import Codec.Wav (importFile)
import Data.Audio (Audio(..))

import Data.Array.Unboxed (elems, assocs)
import Data.Array.IArray (Array(..), listArray, array)
import Data.Complex (magnitude)
import Data.Int (Int64, Int32, Int16, Int8)

import Common (importantFrequencies, breakChunks, slice)
import SongId (ft, mapChunk2Freqs, chunkIdentifier)

type WordSize = Int32
chunkSize = 5000
file = ""

main :: IO ()
main = do
    input <- importFile file
    song_dft <- case input :: Either String (Audio WordSize) of
        Left err -> do
            putStrLn err
            return [(array (0,0) [])]
        Right a@(Audio _ _ samples) -> return (ft chunkSize samples)

    putStrLn $ show (length song_dft) ++ " " ++ show (length . elems $ head song_dft)

    let ys = map (mapChunk2Freqs . elems) song_dft

    let freqs = importantFrequencies 5 (8 * chunkSize `quot` 50 `quot` 2)
    let chunks = map (chunkIdentifier freqs) ys
    foldl (\acc x -> acc >> print x) (return ()) $ chunks


text:
{-|
Module: MachineLearning.NeuralNetwork.Topology
Description: Neural Network's Topology
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Neural Network's Topology
-}

module MachineLearning.NeuralNetwork.Topology
(
  Topology
  , LossFunc(..)
  , makeTopology
  , loss
  , propagateForward
  , propagateBackward
  , numberOutputs
  , initializeTheta
  , initializeThetaIO
  , initializeThetaM
  , flatten
  , unflatten
)

where

import Control.Monad (zipWithM)
import Data.List (foldl')
import qualified Control.Monad.Random as RndM
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA
import MachineLearning.Types (R, Vector, Matrix)
import MachineLearning.Utils (listOfTuplesToList)
import MachineLearning.NeuralNetwork.Layer (Layer(..), Cache(..))
import MachineLearning.NeuralNetwork.Regularization (Regularization, forwardReg, backwardReg)

-- | Loss function's type.
-- Takes x, weights and y.
type LossFunc = Matrix -> Matrix -> R

-- | Neural network topology has at least 2 elements: numver of input and number of outputs.
-- And sizes of hidden layers between 2 elements.
data Topology = Topology [(Int, Int)] [Layer] LossFunc

-- | Makes Neural Network's Topology.
-- Takes number of inputs, list of hidden layers, output layer and loss function.
makeTopology :: Int -> [Layer] -> Layer -> LossFunc -> Topology
makeTopology nInputs hiddenLayers outputLayer lossFunc =
  let layers = hiddenLayers ++ [outputLayer]
      layerSizes = nInputs : (map lUnits layers)
      sizes = getThetaSizes layerSizes
  in Topology sizes layers lossFunc

-- | Calculates loss for the given topology.
-- Takes topology, regularization, x, weights, y.
loss :: Topology -> Regularization -> Matrix -> [(Matrix, Matrix)] -> Matrix -> R
loss (Topology _ _ lf) reg x weights y =
  let lossValue = lf x y
      regValue = forwardReg reg weights
  in lossValue + regValue

-- | Implementation of forward propagation algorithm.
propagateForward :: Topology -> Matrix -> [(Matrix, Matrix)] -> (Matrix, [Cache])
propagateForward (Topology _ layers _) x thetaList = foldl' f (x, []) $ zip thetaList layers
  where f (a, cs) (theta, hl) =
          let (a', cache) = forwardPass hl a theta
          in (a', cache:cs)

-- | Makes one forward step for the given layer.
forwardPass :: Layer -> Matrix -> (Matrix, Matrix) -> (Matrix, Cache)
forwardPass layer a (b, w) = (a', Cache z a w)
  where z = lForward layer a b w
        a' = lActivation layer z

-- | Implementation of backward propagation algorithm.
propagateBackward :: Topology -> Regularization -> Matrix -> [Cache] -> Matrix -> [(Matrix, Matrix)]
propagateBackward (Topology _ layers _) reg scores (cache:cacheList) y = gradientList
  where cache' = Cache scores (cacheX cache) (cacheW cache)
        cacheList' = cache':cacheList
        gradientList = snd $ foldl' f (y, []) $ zip cacheList' $ reverse layers
        f (da, grads) (cache, hl) =
          let (da', db, dw) = backwardPass hl reg da cache
          in (da', (db, dw):grads)

-- | Makes one backward step for the given layer.
backwardPass :: Layer -> Regularization -> Matrix -> Cache -> (Matrix, Matrix, Matrix)
backwardPass layer reg da cache = (da', db, dw')
  where delta = lActivationGradient layer (cacheZ cache) da
        (da', db, dw) = lBackward layer delta cache
        dw' = dw + (backwardReg reg (cacheW cache))

-- | Returns number of outputs of the given topology.
numberOutputs :: Topology -> Int
numberOutputs (Topology nnt _ _) = fst $ last nnt

-- | Returns dimensions of weight matrices for given neural network topology
getThetaSizes :: [Int] -> [(Int, Int)]
getThetaSizes nn = zipWith (\r c -> (r, c)) (tail nn) nn

-- | Create and initialize weights vector with random values
-- for given neural network topology.
-- Takes a seed to initialize generator of random numbers as a first parameter.
initializeTheta :: Int -> Topology -> Vector
initializeTheta seed topology = RndM.evalRand (initializeThetaM topology) gen
  where gen = RndM.mkStdGen seed

-- | Create and initialize weights vector with random values
-- for given neural network topology inside IO Monad.
initializeThetaIO :: Topology -> IO Vector
initializeThetaIO = RndM.evalRandIO . initializeThetaM

-- | Create and initialize weights vector with random values
-- for given neural network topology inside RandomMonad.
initializeThetaM :: RndM.RandomGen g => Topology -> RndM.Rand g Vector
initializeThetaM topology = flatten <$> initializeThetaListM topology

-- | Create and initialize list of weights matrices with random values
-- for given neural network topology.
initializeThetaListM :: RndM.RandomGen g => Topology -> RndM.Rand g [(Matrix, Matrix)]
initializeThetaListM (Topology sizes layers _) = zipWithM lInitializeThetaM layers sizes

-- | Flatten list of matrices into vector.
flatten :: [(Matrix, Matrix)] -> Vector
flatten ms = V.concat $ map LA.flatten $ listOfTuplesToList ms

-- | Unflatten vector into list of matrices for given neural network topology.
unflatten :: Topology -> Vector -> [(Matrix, Matrix)]
unflatten (Topology sizes _ _) v =
  let offsets = reverse $ foldl' (\os (r, c) -> (r+r*c + head os):os) [0] (init sizes)
      ms = zipWith (\o (r, c) -> (LA.reshape r (slice o r), LA.reshape c (slice (o+r) (r*c)))) offsets sizes
      slice o n = V.slice o n v
  in ms

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE BangPatterns #-}
-----------------------------------------------------------------------------
--
-- Module      :  AI.Network.RNN.Util
-- Copyright   :  (c) JP Moresmau
-- License     :  BSD3
--
-- Maintainer  :  JP Moresmau <jp@moresmau.fr>
-- Stability   :  experimental
-- Portability :
--
-- | Utility functions
--
-----------------------------------------------------------------------------

module AI.Network.RNN.Util where

import Control.Parallel.Strategies
import Control.Monad.Random

import Numeric.LinearAlgebra.HMatrix

import Data.List
import Data.Ord
import qualified Data.Set as S

-- | Calculate the product of a matrix by a vector, with both represented by a list
listMProd :: (Num a) => [a] -> [a] -> [a]
listMProd mdt vdt = go mdt vdt 0
  where
    go [] _  !s = [s]
    go ls [] !s = s : go ls vdt 0
    go (y:ys) (x:xs) !ix = go ys xs (y*x+ix)
--    map (foldr (\(a,b) c->a*b+c) 0) $ go mdt vdt []
--  where
--    go [] _  s = [s]
--    go ls [] s = s : go ls vdt []
--    go (y:ys) (x:xs) acc = go ys xs ((y,x):acc)
--listMProd mdt vdt = let (_,l,lst)= foldl' go (vdt,[],0) mdt
--    in reverse $ lst:l
--    where
--        go ([],acc,ix) m =
--            let (x:xs)=vdt
--            in (xs,ix:acc,x*m)
--        go ((x:xs),acc,ix) m = (xs,acc,x*m+ix)
--listMProd mdt vdt = let (_,l,lst)= foldr go (rvdt,[],0) mdt
--    in lst:l
--    where
--        rvdt = reverse vdt
--        go m ([],acc,ix) =
--            let (x:xs)=rvdt
--            in (xs,ix:acc,x*m)
--        go m ((x:xs),acc,ix) = (xs,acc,x*m+ix)
--listMProd mdt vdt = map sum $ takes1 (length vdt) $ zipWith (*) mdt (cycle vdt)

-- | Calculate the mean of a list
calcMeanList :: (Fractional a) => [a] -> a
calcMeanList = uncurry (/) . foldr (\e (s,c) -> (e+s,c+1)) (0,0)

takes1 :: Int -> [a] -> [[a]]
takes1 _ [] = []
takes1 s xs = let
    (one,two) = splitAt s xs
    in one : takes1 s two

-- | Split a given list into a series of list whose length is given by the first argument
takes :: [Int] -> [a] -> [[a]]
takes [] _ = []
takes (_:sz) [] = [] : takes sz []
takes (s:sz) xs = let
    (one,two) = splitAt s xs
    in one : takes sz two

-- | It's a simple, differentiable sigmoid function.
sigmoid :: Floating a => a -> a
sigmoid x = 1 / (1 + exp (-x))
{-# INLINE sigmoid #-}

-- | Normalize a list of values between 0 and 1
normalize :: (Floating b,Ord b)=> [(a,b)] -> [(a,b)]
normalize xs =
    let
        ws = map snd xs
        mi = minimum ws
        mx = maximum ws
        df = mx-mi
     in  map (\(a,w)->(a,(w-mi)/df)) xs

-- | Parallel zipwith
parZipWith :: (a -> b -> c) -> [a] -> [b] -> [c]
parZipWith f xs ys = withStrategy (parList rseq) $ zipWith f xs ys

-- | Parallel zipWith3
parZipWith3 :: (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
parZipWith3 f xs ys zs = withStrategy (parList rseq) $ zipWith3 f xs ys zs

-- | Number of binary digits necessary to store the given number
binaryDigits :: (Num a, Ord a) => a -> a
binaryDigits a = go 1 1
    where go nb sz
            | nb == a = sz
            | nb > a  = sz - 1
            | otherwise = go (nb*2) (sz+1)

-- | Number of bits needed to store in a sparse manner (with 2 bits) the given number
sparseSize :: (Integral a, Ord a) => a -> a
sparseSize a = go 1
      where
        go nb =
            let pos = nb * (nb - 1) `div` 2
            in if pos >= a
                then nb
                else go (nb+1)

-- | Sparse representation: store all numbers from 0 to n using 2 bits
-- Returns the location of each two bits
sparse :: Int -> [(Int,Int)]
sparse n = let
    sz = sparseSize n
    in take n $ iterate (gen sz) (0,1)
    where gen sz (a,b) =
            if b==sz-1 then (a+1,a+2)
                       else (a,b+1)

-- | Standard Normal distribution
stdNormal :: (Monad m,RandomGen g) => RandT g m (Double,Double)
stdNormal = do
    (x1,x2,w1)<-gen
    let w  = sqrt( (-2 * log( w1 ) ) / w1 )
        y1 = x1 * w;
        y2 = x2 * w;
    return (y1,y2)
    where
        ranx = do
            r1 <- getRandomR (0,1)
            return (2 * r1 - 1)
        gen = do
            x1 <- ranx
            x2 <- ranx
            let w = x1 * x1 + x2 * x2
            if w>=1 then gen else return (x1,x2,w)

softmax :: [Double] -> [Double]
softmax is =
    let den = sum $ map exp is
    in map (( / den) . exp) is

softmaxV :: Vector Double -> Vector Double
softmaxV is =
    let den = sumElements $ cmap exp is
    in cmap (( / den) . exp) is

euclidian :: Vector Double -> Vector Double -> Double
euclidian v1 v2 = sumElements $ cmap (**2) (v1 - v2)

equilateralEncoding :: Int -> Matrix Double
equilateralEncoding n =
    let z1 = replicate (n-2) 0
        m = (n><(n-1)) (-1: z1 ++ (1 : z1) ++ repeat 0)
    in foldl' pass m [2..n-1]
    where
        r :: Double -> Double
        r k = -1 / k
        f :: Double -> Double
        f k = sqrt (k*k -1) / k
        pass :: Matrix Double -> Int -> Matrix Double
        pass m1 ik =
            let k = fromIntegral ik
                f1 = f k
                a0 = concatMap (\i->map (\j->((i,j),f1)) [0..ik-2]) [0..ik-1]
                m2 = accum m1 (*) a0
                r1 = r k
                a1 = map (\x->((x,ik-1),r1)) [0..ik-1] ++ [((ik,ik-1),1)]
                m3 = accum m2 const a1
            in m3

equilateralDecoding :: Matrix Double -> Vector Double -> Int
equilateralDecoding m os =
    let rs = zip (toRows m) [0..]
        es = map (\(vs,ix)->(euclidian vs os,ix)) rs
    in snd $ minimumBy (comparing fst) es

roundTo :: (Fractional a, Integral b, RealFrac r) =>
                 b -> r -> a
roundTo n f=  (fromInteger $ round $ f * (10^n)) / (10.0^^n)

ordNub :: (Ord a) => [a] -> [a]
ordNub l = go S.empty l
     where
       go _ []     = []
       go s (x:xs) = if x `S.member` s then go s xs
                                     else x : go (S.insert x s) xs

max_stars_repo_name: JPMoresmau/rnn

max_issues_repo_name: JPMoresmau/rnn

max_forks_repo_name: JPMoresmau/rnn


text:
{-# OPTIONS_GHC -Wall            #-}

{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE OverloadedLists     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeOperators       #-}

module Heat2D (main, a1, a2) where

import Data.Number.Symbolic
import Data.Proxy

import qualified Naperian as N
import qualified GHC.TypeLits as M
import           Data.Functor
import           System.IO

import           Numeric.Sundials
import           Numeric.LinearAlgebra

kx, ky :: Floating a => a
kx = 0.5
ky = 0.75

-- x mesh spacing
-- y mesh spacing
dx :: Floating a => a
dx = 1 / (fromIntegral nx - 1)

dy :: Floating a => a
dy = 1 / (fromIntegral ny - 1)

c1, c2 :: Floating a => a
c1 = kx/dx/dx
c2 = ky/dy/dy

preA2 :: forall b m n . (M.KnownNat m, M.KnownNat n, Num b) =>
        N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] b
preA2 = N.Prism $ N.Prism $ N.Prism $ N.Prism $ N.Scalar $
      N.viota @m <&> (\(N.Fin x) ->
      N.viota @n <&> (\(N.Fin w) ->
      N.viota @m <&> (\(N.Fin v) ->
      N.viota @n <&> (\(N.Fin u) ->
      (f m n x w v u)))))
        where
          m = fromIntegral $ M.natVal (undefined :: Proxy m)
          n = fromIntegral $ M.natVal (undefined :: Proxy n)
          f p q i j k l | i == 0               = 0
                        | j == 0               = 0
                        | i == p - 1           = 0
                        | j == q - 1           = 0
                        | k == i - 1 && l == j = 1
                        | k == i     && l == j = -2
                        | k == i + 1 && l == j = 1
                        | otherwise            = 0

a2Num :: forall a m n . (M.KnownNat m, M.KnownNat n, Floating a) =>
        N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] a
a2Num = N.binary (*) (N.Scalar c2) preA2

a2 :: forall a m n . (M.KnownNat m, M.KnownNat n, Floating a, Eq a) =>
          N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] (Sym a)
a2 = N.binary (*) (N.Scalar $ var "a") preA2

preA1 :: forall b m n . (M.KnownNat m, M.KnownNat n, Num b) =>
        N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] b
preA1 = N.Prism $ N.Prism $ N.Prism $ N.Prism $ N.Scalar $
      N.viota @m <&> (\(N.Fin x) ->
      N.viota @n <&> (\(N.Fin w) ->
      N.viota @m <&> (\(N.Fin v) ->
      N.viota @n <&> (\(N.Fin u) ->
      (f m n x w v u)))))
        where
          m = fromIntegral $ M.natVal (undefined :: Proxy m)
          n = fromIntegral $ M.natVal (undefined :: Proxy n)
          f :: Int -> Int -> Int -> Int -> Int -> Int -> b
          f p q i j k l | i == 0                   = 0
                        | j == 0                   = 0
                        | i == p - 1               = 0
                        | j == q - 1               = 0
                        | k == i     && l == j - 1 = 1
                        | k == i     && l == j     = -2
                        | k == i     && l == j + 1 = 1
                        | otherwise                = 0

a1Num :: forall a m n . (M.KnownNat m, M.KnownNat n, Floating a) =>
        N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] a
a1Num = N.binary (*) (N.Scalar c1) preA1

a1 :: forall a m n . (M.KnownNat m, M.KnownNat n, Floating a, Eq a) =>
           N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] (Sym a)
a1 = N.binary (*) (N.Scalar $ var "b") preA1

h :: forall m n a . (M.KnownNat m, M.KnownNat n, Floating a) =>
                     N.Hyper '[N.Vector m, N.Vector n] a
h = N.Prism $ N.Prism $ N.Scalar $
     N.viota @n <&> (\(N.Fin x) ->
     N.viota @m <&> (\(N.Fin w) ->
     sin (pi * (fromIntegral w) * dx)
             * sin (2 * pi * (fromIntegral x) * dy)))

-- spatial mesh size
nx, ny :: Int
nx = 30
ny = 60

bigA :: Matrix Double
bigA = fromLists $
       fmap (N.elements . N.Prism . N.Prism . N.Scalar) $
       N.elements $ N.crystal $ N.crystal $
       N.binary (+) (a2Num @Double @60 @30) (a1Num @Double @60 @30)

b :: Vector Double
b = fromList $ N.elements (h @30 @60 @Double)

t0, tf :: Double
t0 = 0.0
tf = 0.3

bigNt :: Int
bigNt = 20

dTout :: Double
dTout = (tf - t0) / (fromIntegral bigNt)

ts :: [Double]
ts = map (dTout *) $ map fromIntegral [1..bigNt]

sol :: Matrix Double
-- sol = odeSolveV SDIRK_5_3_4' Nothing 1.0e-5 1.0e-10 (const bigU') (assoc (nx * ny) 0.0 [] :: Vector Double) (fromList $ ts)
--   where
--     bigU' bigU = bigA #> bigU + b
sol = undefined

main :: IO ()
main = do
  h1 <- openFile "heat2e.000.txt" WriteMode
  mapM_ (hPutStrLn h1) $ map (concatMap (' ':)) $ map (map show) $ toLists sol
  hClose h1
  mapM_ (\i -> putStrLn $ show $ sqrt $ (sol!i) <.> (sol!i) / (fromIntegral nx) / (fromIntegral ny)) ([0 .. length ts - 1] :: [Int])

max_stars_repo_name: idontgetoutmuch/Diffusions

max_issues_repo_name: idontgetoutmuch/Diffusions

max_forks_repo_name: idontgetoutmuch/Diffusions


text:
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
module Main (main) where

import           Data.Complex
import qualified Data.Text           as T
import           Graphics.SvgTree
import           Linear.V2
import           Reanimate

main :: IO ()
main = reanimate $ pauseAtEnd 2
  fourierAnimation_

sWidth :: Double
sWidth = 0.02

piFourier :: Fourier
piFourier = mkFourier piPoints

piPoints :: [RPoint]
piPoints = lineToPoints 500 $
  toLineCommands $ extractPath $ scale 10 $ center $ latexAlign "\\pi"

fourierAnimation_ :: Animation
fourierAnimation_ = mkAnimation 50 $ \t ->
    let fLength = t
        circles = setFourierLength (fLength*maxLength) piFourier
        maxLength = sum $ map magnitude $ take 499 $ drop 1 $ fourierCoefficients piFourier
        phi = fromToS 0 15 t
    in mkGroup
    [ mkBackground "black"
    , drawCircles $ fourierCoefficients $ rotateFourier phi circles
    , withStrokeColor "green" $
      withFillOpacity 0 $
      mkLinePath $ mkFourierOutline circles
    , withFillColor "white" $
      translate (-screenWidth/16*7) (screenHeight/16*7) $
      latex $ T.pack $ "Circles: " ++ show (length $ fourierCoefficients circles)
    ]

newtype Fourier = Fourier {fourierCoefficients :: [Complex Double]}

pointAtFourier :: Fourier -> Complex Double
pointAtFourier = sum . fourierCoefficients

mkFourier :: [RPoint] -> Fourier
mkFourier points = Fourier $ findCoefficient 0 :
    concat [ [findCoefficient n, findCoefficient (-n)] | n <- [1..] ]
  where
    findCoefficient :: Int -> Complex Double
    findCoefficient n =
        sum [ toComplex point * exp (negate (fromIntegral n) * 2 *pi * i*t) * deltaT
            | (idx, point) <- zip [0::Int ..] points, let t = fromIntegral idx/nPoints ]
    i = 0 :+ 1
    toComplex (V2 x y) = x :+ y
    deltaT = recip nPoints
    nPoints = fromIntegral (length points)

-- setFourierCircles :: Double -> Fourier -> Fourier
-- setFourierCircles n _ | n < 1 = error "Invalid argument. Need at least one circle."
-- setFourierCircles n (Fourier coeffs) =
--     Fourier $ take iCircles coeffs ++ [coeffs!!iCircles * realToFrac fCircle]
--   where
--     (iCircles, fCircle) = divMod' n 1

setFourierLength :: Double -> Fourier -> Fourier
setFourierLength _ (Fourier []) = Fourier []
setFourierLength len0 (Fourier (first:lst)) = Fourier $ first : worker len0 lst
  where
    worker _len [] = []
    worker len (c:cs) =
      if magnitude c < len
        then c : worker (len - magnitude c) cs
        else [c * realToFrac (len / magnitude c)]

rotateFourier :: Double -> Fourier -> Fourier
rotateFourier phi (Fourier coeffs) =
    Fourier $ worker coeffs (0::Integer)
  where
    worker [] _ = []
    worker (x:rest) 0 = x : worker rest 1
    worker [left] n = worker [left,0] n
    worker (left:right:rest) n =
      let n' = fromIntegral n in
      left * exp (negate n' * 2 * pi * i * phi') :
      right * exp (n' * 2 * pi * i * phi') :
      worker rest (n+1)
    i = 0 :+ 1
    -- n = length coeffs `div` 2
    phi' = realToFrac phi

drawCircles :: [Complex Double] -> SVG
drawCircles circles = mkGroup
    [ worker circles
    , withStrokeWidth sWidth $
      withStrokeColor "white" $
      withStrokeLineJoin JoinRound $
      withFillOpacity 0 $
      mkLinePath [ (x, y) | x :+ y <- scanl (+) 0 circles ] ]
  where
    worker [] = None
    worker (x :+ y : rest) =
      let radius = sqrt(x*x+y*y) in
      mkGroup
      [ withStrokeWidth sWidth $
        withStrokeColor "dimgrey" $
        withFillOpacity 0 $
        mkCircle radius
      , translate x y $ worker rest ]

mkFourierOutline :: Fourier -> [(Double, Double)]
mkFourierOutline fourier =
    [ (x, y)
    | idx <- [0 .. granularity]
    , let x :+ y = pointAtFourier $ rotateFourier (idx/granularity) fourier
    ]
  where
    granularity = 500

max_stars_repo_name: jhrcek/reanimate

max_issues_repo_name: jhrcek/reanimate

max_forks_repo_name: jhrcek/reanimate


text:
{-# LANGUAGE Strict #-}
module Pinwheel.Asterisk where

import Data.Complex

{-# INLINE asteriskGaussianFunc #-}
asteriskGaussianFunc ::
     (RealFloat e) => e -> e -> Int -> e -> e -> e -> Complex e
asteriskGaussianFunc phi rho n a b period =
  (0 :+ (-1)) ^ n *
  exp
    ((-1) *
     (((fromIntegral n ^ 2) / (4 * b) + (pi * rho) ^ 2 / (a * period ^ 2)) :+
      (fromIntegral n * phi)))

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE BangPatterns #-}
module InteractionNetwork
  ( adjacencyMatrix
  , interactionNetwork
  , tcc
  , ctpl
  , tge
  , allPairsBFS
  , spectralRadius
  , spectralRadius'
  , centralities
  , oneHopCentralities
  , InteractionNetwork(..)
  ) where

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data (mkSparse, nCols)

import qualified Data.Sequence as S
import Data.Bits
import Data.Maybe

type Snapshot = AssocMatrix
data InteractionNetwork = INet { numNodes :: Int , graph :: [Snapshot] }

adjacencyMatrix :: Int -> [(Int,Int)] -> Snapshot
adjacencyMatrix numNodes interactions = [((n, n), 0)] ++ map (\(x, y) -> ((x,y),1)) interactions
  where n = numNodes - 1

interactionNetwork :: Int -> [[(Int, Int)]] -> InteractionNetwork
interactionNetwork size interactions = INet { numNodes = size, graph = map (adjacencyMatrix size) interactions }

allPairsBFS :: InteractionNetwork -> [[(Int, Maybe Int)]]
allPairsBFS network = [ bfs network u | u <- [0..(n-1)] ]
  where n = numNodes network

-- NOTE: the adjacency matrix is symmetric, symmetric matricies have real-valued eigenvalues.
spectralRadius :: InteractionNetwork -> Double
spectralRadius network = foldr (\s acc -> let rho = maxElement . cmap abs . eigenvaluesSH . trustSym $ toDense s in
                                   if rho > acc then rho else acc) 0.0 (graph network)

-- minimum sr over all timesteps
spectralRadius' :: InteractionNetwork -> Double
spectralRadius' network =
  minimum $ map (maxElement . cmap abs . eigenvaluesSH . trustSym . toDense) (graph network)

oneHopCommunicability :: InteractionNetwork -> Matrix Double
oneHopCommunicability network =
  communicability' g (ident n)
  where a = 0.1
        n = numNodes network
        g = graph network

        communicability' []       q = q
        communicability' (s:rest) q =
          -- by replacing (I - aA) with (I + aA) we can enforce only
          -- on hop per time step (apparently)
          let q' = q <> ((ident n) + (a*(toDense s))) in
            communicability' rest $ q' / scalar (norm_2 q')

-- | Compute the normalized communicability matrix.
communicability :: InteractionNetwork -> Matrix Double
communicability network =
  communicability' g (ident n)
  -- take a < 1 per Grindrod 2010
  -- a = 0.1 seems sufficient for teleportation, brownian, and crw models
  where a = 0.1 -- XXX: a < 1/max (ρ(A[k])) for k in 1..length g
        n = numNodes network
        g = graph network

        communicability' []       q = q
        communicability' (s:rest) q =
          let q' = q <> (inv $ (ident n) - (a*(toDense s))) in
            communicability' rest $ q' / scalar (norm_2 q')

centralities :: InteractionNetwork -> (Vector Double, Vector Double)
centralities network = (broadcastCentralities network q, receiveCentralities network q)
  where q = communicability network

oneHopCentralities :: InteractionNetwork -> (Vector Double, Vector Double)
oneHopCentralities network = (broadcastCentralities network q, receiveCentralities network q)
  where q = oneHopCommunicability network

broadcastCentralities :: InteractionNetwork -> Matrix Double -> Vector Double
broadcastCentralities network q =
  q #> (vector $ replicate (numNodes network) 1)

receiveCentralities :: InteractionNetwork -> Matrix Double -> Vector Double
receiveCentralities network q =
  (vector $ replicate (numNodes network) 1) <# q

-- | Compute the temporal correlation coefficient
--
-- Should be 1 if every snapshot is the same
--
-- Follows the equation from Clausset paper, accounting for vertices
-- that don't have any edges at all.
tcc :: InteractionNetwork -> Double
tcc net = (sum $ foldr (\xs acc -> zipWith (+) xs acc) (repeat 0.0) [ cit t t' | (t,t') <- zip g (tail g) ])
          / (fromIntegral $ (timesteps-1) * numNodes net)
  where timesteps = length (graph net)
        g         = graph net

        cit :: Snapshot -> Snapshot -> [Double]
        cit t t' = let rows  = toRows $ toDense t
                       rows' = toRows $ toDense t'
                   in
                     zipWith cij rows rows'
        cij :: Vector Double -> Vector Double -> Double
        cij r r' = if sr == 0 || sr' == 0 then 1 else  (r <.> r') / (sqrt $ (sr) * (sr'))
          where sr  = sumElements r
                sr' = sumElements r'

-- | Compute the characteristic temporal path length of the graph.
--
-- Any nodes that do not have a temporal path are excluded from the
-- average.
ctpl :: InteractionNetwork -> [[(Int, Maybe Int)]] -> Double
ctpl network bst = total / (fromIntegral . length $ pathLengths)
  where pathLengths = map (\(_, (Just x)) -> x) . filter (\(_, x) -> case x of
                                                                       Just 0 -> False
                                                                       Just x -> True
                                                                       Nothing -> False) $
                      concat bst
        total = fromIntegral $ sum pathLengths

-- | Temporal Global Efficiency
tge :: InteractionNetwork -> [[(Int, Maybe Int)]] -> Double
tge network bst = total / fromIntegral (n*(n-1))
  where n     = numNodes network
        total = sum $ map rpSum bst

        rpSum :: [(Int, Maybe Int)] -> Double
        rpSum = foldr (\(_, path) !acc -> case path of
                                           Just 0 -> acc -- ignore self paths
                                           Just x -> (1/fromIntegral x) + acc
                                           Nothing -> acc) 0.0

-- | Geeralized BFS for time varying graphs. Returns the length of the
-- shortest temporal path from a source node to all other node in the
-- network.
--
-- Any node `x' that is unreachable will apear in the result as (x,
-- Nothing).
bfs :: InteractionNetwork -> Int -> [(Int, Maybe Int)]
bfs net source = bfs' 1 (graph net) (bit source) (S.singleton source) [(source, Just 0)]
  where bfs' :: Int -> [Snapshot] -> Integer -> S.Seq Int -> [(Int, Maybe Int)] -> [(Int, Maybe Int)]
        bfs' t net@(s:future) !visited queue tpLengths =
          let front S.:< rest = S.viewl queue
              frontier   = explore visited s front -- get all adjacent nodes to front that have not already been visited
              visited'   = visited .|. (foldr (\x acc -> setBit acc x) zeroBits frontier)
              queue'     = rest S.>< (S.fromList $ frontier ++ [front])
              tpLengths' = tpLengths ++ map (\x -> (x, Just t)) frontier
          in
            if front /= source -- source acts as a sentinel, marking the "begining" of the queue
            then bfs' t net visited' queue' tpLengths'
            else bfs' (t+1) future visited' queue' tpLengths'
        bfs' _ [] visited _ tpLengths = [ (x, Nothing) | x <- [0..((numNodes net) - 1)], not $ testBit visited x ] ++ tpLengths

        -- | Get all nodes adjacent to n that have not beed discovered yet.
        explore :: Integer -> Snapshot -> Int -> [Int]
        explore visited s source = map (snd . fst) $ filter (\((x,y),z) -> (z /= 0) && (x == source) && (not $ testBit visited y)) s

max_stars_repo_name: wfvining/tvg-smallworld

max_issues_repo_name: wfvining/tvg-smallworld

max_forks_repo_name: wfvining/tvg-smallworld


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TemplateHaskell       #-}

module Lib where

import           Control.Monad.Except                      (MonadError,
                                                            runExceptT,
                                                            throwError)
import           Control.Monad.IO.Class                    (MonadIO, liftIO)
import           Data.ByteString                           (ByteString)
import           Data.FileEmbed                            (embedFile)
import           Data.Monoid                               ((<>))
import           Data.Scientific                           (Scientific,
                                                            toRealFloat)
import           Debug.Trace
import           GHC.TypeLits                              ()
import           Graphics.Plot                             (mplot)
import           Graphics.Rendering.Chart.Backend.Diagrams (toFile)
import           Graphics.Rendering.Chart.Easy             (def, layout_title,
                                                            line, plot, points,
                                                            re, (.=))
import           Numeric.GSL.Minimization                  (MinimizeMethod (NMSimplex),
                                                            MinimizeMethodD (SteepestDescent, VectorBFGS2),
                                                            minimizeD,
                                                            minimizeV,
                                                            minimizeVD)
import           Numeric.LinearAlgebra                     (Container, accum,
                                                            matFunc, tr, ( #> ),
                                                            (<.>))
import qualified Numeric.LinearAlgebra                     as Matrix
import           Numeric.LinearAlgebra.Data                (Matrix, Vector,
                                                            asColumn, asRow,
                                                            cols, cond,
                                                            dropColumns,
                                                            dropRows, rows,
                                                            takeColumns,
                                                            takeRows, (><),
                                                            (|||))
import qualified Numeric.LinearAlgebra.Data                as MData
import           Numeric.LinearAlgebra.Static              (L, dim, matrix, mul)
import qualified Numeric.LinearAlgebra.Static              as Static
import           Text.Megaparsec                           (Dec, ParseError,
                                                            Token, char,
                                                            newline, parse,
                                                            sepBy)
import           Text.Megaparsec.ByteString                (Parser)
import           Text.Megaparsec.Lexer                     (number)

data Error =
  ParseFailed (ParseError (Token ByteString) Dec)
  deriving (Show, Eq)

------------------------------------------------------------
-- Parsing.
rawDataParser :: Parser [[Scientific]]
rawDataParser =
  filter (not . null) <$> (number `sepBy` char ',') `sepBy` newline

toMatrix :: [[Scientific]] -> Matrix Double
toMatrix = Matrix.fromLists . fmap (fmap toRealFloat)

wineData :: ByteString
wineData = $(embedFile "data/wine.data")

wineParser :: Parser (Matrix Double)
wineParser = toMatrix <$> rawDataParser

loadData
  :: (MonadError Error m, MonadIO m)
  => m (Matrix Double)
loadData = either (throwError . ParseFailed) pure (parse wineParser "" wineData)

------------------------------------------------------------
-- Matrices
a :: L 3 2
a = matrix [1 .. 6]

b :: L 2 1
b = matrix [1 .. 2]

c :: L 3 1
c = mul a b

------------------------------------------------------------
-- Linear Regression
initialTheta :: Int -> Vector Double
initialTheta n = MData.vector $ replicate n 0

------------------------------------------------------------
-- Charts
signal :: [Double] -> [(Double, Double)]
signal xs =
  [(x, (sin (x * 3.14159 / 45) + 1) / 2 * sin (x * 3.14159 / 5)) | x <- xs]

writeChart :: IO ()
writeChart =
  toFile def "example.svg" $ do
    layout_title .= "Amplitude Modulation"
    plot (line "am" [signal [0,0.5 .. 800]])
    plot (points "am points" (signal [0,7 .. 800]))

------------------------------------------------------------
-- Main
main :: IO ()
main = do
  result <-
    runExceptT $ do
      dataset <- loadData
      process dataset
  print result

data Result = Result
  { crossValidationSet    :: Matrix Double
  , crossValidationResult :: Vector Double
  , finalTheta            :: Vector Double
  , finalHyp              :: Vector Double
  , finalCost             :: Double
  } deriving (Show)

process
  :: (MonadIO m, MonadError Error m)
  => Matrix Double -> m Result
process dataSet = do
  let is n v = cond v n 0 1 0
  let target = 1
  let scale :: Matrix Double =
        accum
          (MData.ident 13)
          (*)
          [((0, 0), 0.1), ((3, 3), 0.1), ((4, 4), 0.01), ((12, 12), 0.001)]
  let m = rows dataSet
  let f = cols dataSet
  let y = is target $ MData.flatten $ takeColumns 1 dataSet
  let x = Matrix.col (replicate m 1.0) ||| ((dropColumns 1 dataSet) <> scale)
  let cvSize = 10
  let trainingX = dropRows cvSize x
  let cvX = takeRows cvSize x
  let trainingY = MData.fromList $ drop cvSize $ MData.toList y
  let cvY = MData.fromList $ take cvSize $ MData.toList y
  let theta = initialTheta f
  let (finalTheta, path) =
        minimizeVD
          SteepestDescent -- VectorBFGS2
          10e-6
          50
          10e-6
          0.1
          (costFn x y)
          (gradFn x y)
          theta
  pure $
    Result
      cvX
      cvY
      finalTheta
      (hypothesis cvX finalTheta)
      (costFn cvX cvY finalTheta)

sigmoid
  :: Floating a
  => a -> a
sigmoid term = 1.0 / (1.0 + exp (-term))

hypothesis :: Matrix Double -> Vector Double -> Vector Double
hypothesis x theta = sigmoid (x #> theta)

costFn :: Matrix Double -> Vector Double -> Vector Double -> Double
costFn x y theta = traceShowId $ all
  where
    all :: Double
    all = vsum (leftTerm - rightTerm) / fromIntegral (rows x)
    h = hypothesis x theta
    leftTerm = (-y) * log h
    rightTerm = (1.0 - y) * log (1.0 - h)

gradFn :: Matrix Double -> Vector Double -> Vector Double -> Vector Double
gradFn x y theta = (tr x #> (hypothesis x theta - y)) / fromIntegral (rows x)

vsum :: Vector Double -> Double
vsum = sum . MData.toList

max_stars_repo_name: krisajenkins/sketch

max_issues_repo_name: krisajenkins/sketch

max_forks_repo_name: krisajenkins/sketch


text:
{-
  An attempt to implement a basic case of Cooley-Tukey FFT algorithm.
 -}
{-# LANGUAGE ViewPatterns #-}
module Main
  ( main
  ) where

import Control.DeepSeq
import Control.Monad
import Control.Monad.State.Strict
import Data.Bits
import Data.Complex
import Data.Time.Clock
import Statistics.Sample
import System.Random
import Text.Printf

import qualified Data.Vector as V
import qualified Data.Vector.Mutable as VM

type Cpx = Complex Double

-- e^((c pi i) * (k / N))
expAuxCommon :: Int -> Int -> Int -> Cpx
expAuxCommon c k n = cis theta
  where
    theta = pi * fromIntegral (c * k) / fromIntegral n

-- Computes e^((-2 pi i) * (k / N))
expFft :: Int -> Int -> Cpx
expFft = expAuxCommon (-2)

-- Computes e^((2 pi i) * (k / N))
expIfft :: Int -> Int -> Cpx
expIfft = expAuxCommon 2

splitByParity :: V.Vector Cpx -> (V.Vector Cpx, V.Vector Cpx)
splitByParity vs = (evens, odds)
  where
    l = V.length vs
    oddCount = (l + 1) `quot` 2
    evenCount = l - oddCount
    evens =
      V.fromListN evenCount $ (vs V.!) <$> [0,2..]
    odds =
      V.fromListN oddCount $ (vs V.!) <$> [1,3..]

-- TODO: For now we assume that vector length is a power of two,
-- this can be easily extended to any length, by pretending out-of-range
-- values are 0.
gDitFft :: (Int -> Int -> Cpx) -> V.Vector Cpx -> V.Vector Cpx
gDitFft expF = fix $
  \impl vs ->
    let l = V.length vs
    in if l <= 1
      then vs
      else V.create $ do
        let (impl -> es, impl -> os) = splitByParity vs
            hf = l `quot` 2
        vec <- VM.unsafeNew l
        forM_ [0 .. hf - 1] $ \k -> do
          let a = es V.! k
              b = expF k l * os V.! k
          VM.unsafeWrite vec k (a + b)
          VM.unsafeWrite vec (hf + k) (a - b)
        pure vec

ditFft :: V.Vector Cpx -> V.Vector Cpx
ditFft = gDitFft expFft . rightPadZeros

iditFft :: V.Vector Cpx -> V.Vector Cpx
iditFft (rightPadZeros -> vs) = V.map (/ fromIntegral l) (iditFftAux vs)
  where
    l = V.length vs
    iditFftAux :: V.Vector Cpx -> V.Vector Cpx
    iditFftAux = gDitFft expIfft

{-
  right padding zeros in the end of a matrix to make the length of the vector
  a power of 2.

  special case: returns empty vector when the input is empty.
 -}
rightPadZeros :: V.Vector Cpx -> V.Vector Cpx
rightPadZeros vs
    | l <= 1 = vs
    | isPowerOf2 = vs
    | otherwise =
        let closestPowOf2Gt =
              {-
                Find the closest power of two that is greater or equal to input v.
                only works when v > 1 and v is not a power of two.
               -}
              unsafeShiftL 1 (finiteBitSize l - countLeadingZeros l)
        in vs <> V.fromListN (closestPowOf2Gt - l) (repeat 0)
  where
    l = V.length vs
    isPowerOf2 =
      countLeadingZeros l + countTrailingZeros l + 1 == finiteBitSize l

evaluateOnVector :: V.Vector Cpx -> IO ()
evaluateOnVector vs = do
  let result = iditFft (ditFft vs)
      diffs :: V.Vector Double
      diffs = V.zipWith (\x y -> magnitude (x - y)) vs result
  printf "StdDev: %.9f\n" (stdDev diffs)

type M = StateT StdGen IO

genRandomR :: Random a => (a, a) -> M a
genRandomR = state . randomR

genCpxVector :: M (V.Vector Cpx)
genCpxVector =
  genRandomR (1024, 8192) >>= genCpxVectorOfSize

genCpxVectorOfSize :: Int -> M (V.Vector Cpx)
genCpxVectorOfSize l = do
  let rDouble = genRandomR (-100,100)
  xs <- replicateM l $ (:+) <$> rDouble <*> rDouble
  pure $ V.fromListN l xs

evalRandomVector :: M ()
evalRandomVector = do
  vs <- genCpxVector
  liftIO $ do
    printf "Vector length: %d\n" (V.length vs)
    evaluateOnVector vs

directConvolve :: V.Vector Cpx -> V.Vector Cpx -> V.Vector Cpx
directConvolve xs ys = V.create $ do
  let lX = V.length xs
      lY = V.length ys
      lZ = lX + lY - 1
  vec <- VM.replicate lZ 0
  forM_ [0 .. lX - 1] $ \i -> do
    let x = V.unsafeIndex xs i
    forM_ [0 .. lY - 1] $ \j -> do
      let y = V.unsafeIndex ys j
          k = i + j
      v <- VM.unsafeRead vec k
      VM.unsafeWrite vec k $! v + x * y
  pure vec

fftConvolve :: V.Vector Cpx -> V.Vector Cpx -> V.Vector Cpx
fftConvolve xsPre ysPre = V.fromListN lZ (V.toList (iditFft rs))
  where
    lX = V.length xsPre
    lY = V.length ysPre
    lZ = lX + lY - 1
    -- TOOD: admittedly this is a bit wasteful...
    xs = V.fromListN lZ (V.toList xsPre <> repeat 0)
    ys = V.fromListN lZ (V.toList ysPre <> repeat 0)
    xs' = ditFft xs
    ys' = ditFft ys
    rs = V.zipWith (*) xs' ys'

evalRandomConvolve :: M ()
evalRandomConvolve = do
  xs <- genCpxVector
  ys <- genCpxVector
  let lX = V.length xs
      lY = V.length ys
      lZ = lX + lY - 1
      zsDirect = directConvolve xs ys
      zsFft = fftConvolve xs ys
      diffs :: V.Vector Double
      diffs = V.zipWith (\x y -> magnitude (x - y)) zsDirect zsFft
  liftIO $ do
    putStrLn $ "Vector lengths (xs,ys,xs * ys): " <> show (lX, lY, lZ)
    _ <- forceAndMeasure "Direct" zsDirect
    _ <- forceAndMeasure "FFT" zsFft
    printf "StdDev: %.9f\n" (stdDev diffs)

forceAndMeasure :: NFData a => String -> a -> IO a
forceAndMeasure tag r = do
  tBefore <- getCurrentTime
  tAfter <- r `deepseq` getCurrentTime
  printf "Time (%s): %.6f\n" tag ((realToFrac $ diffUTCTime tAfter tBefore) :: Double)
  pure r

main :: IO ()
main = do
  let cs = (\x -> x*2 :+ (x*2 + 1)) <$> [0..19]
      vs = V.fromList cs
  evaluateOnVector vs
  g <- newStdGen
  evalStateT (evalRandomVector >> evalRandomConvolve) g

max_stars_repo_name: Javran/misc

max_issues_repo_name: Javran/misc

max_forks_repo_name: Javran/misc


text:
-- ex 1.15 (1.4, 1.17)

import qualified Data.Vector as V

import Statistics.Sample.Histogram (histogram)
import Statistics.Sample (mean, varianceUnbiased, stdDev, range)

-- Top 40 stocks on the Over-The-Counter marked:
samples :: V.Vector Double
samples = V.fromList
  [ 11.88,  6.27,  5.49,  4.81,  4.40,  3.78,  3.44,  3.11,  2.88,  2.68
  ,  7.99,  6.07,  5.26,  4.79,  4.05,  3.69,  3.36,  3.03,  2.74,  2.63
  ,  7.15,  5.98,  5.07,  4.55,  3.94,  3.62,  3.26,  2.99,  2.74,  2.62
  ,  7.13,  5.91,  4.94,  4.43,  3.93,  3.48,  3.20,  2.89,  2.69,  2.61
  ]

-- 1.2 histogram
hist :: (V.Vector Double, V.Vector Int)
hist = histogram 15 samples

main :: IO ()
main = do
  let (lowerBs, hs) = hist
  putStrLn $ "The lower bounds: " ++ show lowerBs
  putStrLn $ "The histogram: " ++ show hs

  let av = mean samples
      sd = stdDev samples
  putStrLn $ "The average value is: " ++ show av ++ " and the standard deviation is: " ++ show sd

  let f k =
       V.length . V.filter (\x -> (av-k*sd)<x && x< (av+k*sd)) $ samples
  putStrLn $ "Our sample size is: " ++ show (V.length samples)
  putStrLn $ "Within one sigma:   " ++ show (f 1)
  putStrLn $ "Within two sigma:   " ++ show (f 2)
  putStrLn $ "Within three sigma: " ++ show (f 3)

  putStrLn $ "range/4 is:           " ++ show (range samples / 4)
  putStrLn $ "The average value is: " ++ show av

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, ConstraintKinds, ExistentialQuantification, GADTs, RankNTypes, MultiParamTypeClasses   #-}
{-# LANGUAGE ImpredicativeTypes, RankNTypes, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : Count-Min Sketch
-- | Creator: Xiao Ling
-- | Created: 12/17/2015
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module CountMinSketch where

import Control.Monad.Random.Class
import Control.Monad.Random
import Control.Monad.State

import Data.Matrix
import Data.Conduit
import qualified Data.Conduit.List as Cl
import qualified Data.Matrix as M

import Core
import Utils
import Statistics

{-----------------------------------------------------------------------------
  Types
------------------------------------------------------------------------------}

type D        = Int              -- * num rows
type W        = Int              -- * num cols
type P        = Int              -- * some prime that is 2^n - 1
type Event    = Int              -- * some input event from universe:
                                 -- * {1, .., i, .., n}

type Sketch   = Matrix Int       -- * D x W matrix  of counters
type ABs      = Matrix Int       -- * 2 x D matrix  of a_j and b_j for j = 1..d
type ColIdx   = [Int]            -- * indices of counters to increment
type Hash     = Event -> ColIdx  -- * Maps event to d x 1 list of column indices

-- * Algorithm State
type AState    = (Sketch,Hash,(D,W,P))

-- * the algorithm uses `some m` equipped with state and randomness
type Some m    = (MonadRandom m, MonadState AState m)

-- * Algorithm outputes an `Estimator` that
-- * estimates the fequence of `event` i from  universe {1 ... i... n},
type Estimator = Event -> Int

{-----------------------------------------------------------------------------
    Count Min Sketch Conduit
------------------------------------------------------------------------------}

-- * Batch mode
sketch' :: EpsDelta -> Event -> Batch Event IO Estimator
sketch' t n = sketch t n `using` eval

-- * Streaming mode
sketch :: Some m => EpsDelta -> Event -> Streaming Event m ()
sketch t n = inits False t n $= awaitWith update'

-- * initalize algorithm given confidence/error `t` and maximum event `n`
inits :: Some m => Bool -> EpsDelta -> Event -> Conduit Event m Event
inits False t n = lift (inits' t n) >> inits True t n
inits _     _ _ = pass

-- * a concrete `eval`uation of Some `m`
-- * run with vacuous inital condition to be thrown away
eval :: StateT AState (Rand StdGen) () -> IO Estimator
eval m = fmap (\(s,h,_) -> \i -> freq i s h)
        . evalRandIO $ execStateT m $ (ones 1 1, const [], (0,0,0))

{-------------------------------------------------------------------------------
   CountMinSketch Primitives
--------------------------------------------------------------------------------}

-- * initalize algorithm with (w,d,p), zero counters, and hash function
inits' :: Some m => EpsDelta -> Event -> m ()
inits' ed n = do
    let t = newDwp ed n
    setParam t
    hash    <- newHash t
    let sketch = newSketch t
    setHash   hash
    setSketch sketch

-- * on event `i`, update the `sketch`
update' :: Some m => Event -> m ()
update' i = do
    pm      <- askParam
    sketch  <- getSketch
    hash    <- askHash
    let incr = toMask (hash i) pm
    setSketch $ incr .+. sketch

{-------------------------------------------------------------------------------
   Effectful Subroutines
--------------------------------------------------------------------------------}

-- * h_j (i) = (a_j x i + b_j mod p) mod w
-- * sumRow-wise  [ a_1 ... a_j ... a_d ]  .*.  [ i_1 ...  i_d ]     mod p  mod w
-- *              [ b_1 ... b_j ... b_d ]       [ 1   ...  1   ]
-- * result:      [ idx_1, ..., idx_d   ]
-- *               where the indices start at 1
newHash :: MonadRandom m => (D,W,P) -> m Hash
newHash t = do
    ab  <- newAbs t
    return $ toHash t ab

-- * intialize random cofficients a_j, b_j, for d x 2 matrix:
-- * [ a_1 ... a_j ... a_d ]
-- * [ b_1 ... b_j ... b_d ]
newAbs :: MonadRandom m => (D,W,P) -> m ABs
newAbs (d,_,p) = (M.fromList 2 d . take (2*d)) <$> getRandomRs (1,p)

{-----------------------------------------------------------------------------
    Pure Subroutines
------------------------------------------------------------------------------}

-- * given some event `i` from  universe {1 ... i... n}, a table of
-- * counters Sketch `s` and hash function `h`, compute the
-- * frequence of event `i`
freq :: Event -> Sketch -> Hash -> Int
freq i s h = minimum $ go s (h i) 1
    where go s []     _ = []
          go s (c:cs) d = (getElem d c s) : (go s cs $ succ d)

-- * compute d,w, and p given accuracy/confidence level and
-- * max value of event `n`
newDwp :: EpsDelta -> Event -> (D,W,P)
newDwp (ED e d) n = (round $ log2(1/d), round $ 2/e, round $ 2^n -1)

-- * initalize the sketch
newSketch :: (D,W,P) -> Sketch
newSketch (d,w,_) = zeros d w

-- * Note `is` maps some `i` from event universe onto d x 2 matrix:
-- * [ i_1 ... i_2 ]
-- * [ 1   ...  1  ]
toHash :: (D, W, P) -> ABs -> Hash
toHash (d, w, p) ab i = go $ sumR (ab .*. is) ... (\v -> v `mod` p `mod` w)
    where go = fmap (+1) . toList
          is = (i .+ (zeros 1 d)) <-> ones 1 d

-- * given hash function, construct a binary mask of size `w` x `d`
-- * to increment sketch
toMask :: ColIdx -> (D,W,P) -> Sketch
toMask ks (d,w,p) = go w ks
    where
        go w ks = let vec c w = zeros 1 (c-1) <|> single 1 <|> zeros 1 (w - c) in
            case ks of
                c:[] -> vec c w
                c:cs -> vec c w <-> go w cs

{-----------------------------------------------------------------------------
    Managing State
    TODO : redesign so this is not necessary
------------------------------------------------------------------------------}

askParam :: Some m => m (D,W,P)
askParam    = get >>= \(_,_,p) -> return p

setParam :: Some m => (D,W,P) -> m ()
setParam t  = get >>= \(s,h,_) -> put (s,h,t)

getSketch :: Some m => m Sketch
getSketch   = get >>= \(s,_,_) -> return s

setSketch :: Some m => Sketch -> m ()
setSketch s = get >>= \(_,h,t) -> put (s,h,t)

askHash :: Some m => m Hash
askHash     = get >>= \(_,h,_) -> return h

setHash :: Some m => Hash -> m ()
setHash h   = get >>= \(s,_,t) -> put (s,h,t)

{-----------------------------------------------------------------------------
    Find conduit lib function for these:
------------------------------------------------------------------------------}

-- * identity conduit, why can't I find this in the library?
pass :: Monad m => Conduit i m i
pass  = do
    mx <- await
    case mx of
        Nothing -> return ()
        Just x  -> yield x >> pass

awaitWith :: Monad m => (i -> m a) -> ConduitM i o m ()
awaitWith g = do
    mi <- await
    case mi of
        Nothing -> return ()
        Just i  -> lift (g i) >> awaitWith g


text:
module Strategy.Orientation where

import Prelude hiding (Right, Left)
import Numeric.LinearAlgebra ((<>))

import Cube (Color, Side(..), Cube)
import qualified Cube as Cube
import Rotation (MatrixRotation)
import qualified Rotation as Rotation
import RotationPath (CubeMutation)
import qualified RotationPath as RotationPath

{-|
  Orienting a cube so that other strategies have an oriented cube to work with.

  In the context of this module color of a side refers only to the middle field of a side.
|-}

sides :: Cube -> [(Side, Color)]
sides = map (\(s, _, _, c) -> (s, c)) . filter (\(_, x, y, _) -> x == 1 && y == 1)

isColorOnSide :: Color -> Side -> Cube -> Bool
isColorOnSide color side = any (\(s, c) -> (c == color && s == side)) . sides

-- | Returning (prefix, suffix) rotations to do something with a color on top
forColorOnTop :: Color -> Cube -> (MatrixRotation, MatrixRotation)
forColorOnTop color cube
  | isColorOnSide color Front cube = (Rotation.topToBack, Rotation.topToFront)
  | isColorOnSide color Right cube = (Rotation.topToLeft, Rotation.topToRight)
  | isColorOnSide color Back cube = (Rotation.topToFront, Rotation.topToBack)
  | isColorOnSide color Left cube = (Rotation.topToRight, Rotation.topToLeft)
  | isColorOnSide color Bottom cube =
    let r = Rotation.topToFront <> Rotation.topToFront
    in (r, r)
  | otherwise = (Rotation.identity, Rotation.identity)

putColorOnTop :: Color -> CubeMutation
putColorOnTop color cube =
  let (rotation, _) = forColorOnTop color cube
  in RotationPath.rotate rotation cube

withColorOnTop :: Color -> CubeMutation -> CubeMutation
withColorOnTop color λ cube =
  let (prefix, suffix) = forColorOnTop color cube
  in RotationPath.between prefix suffix λ cube

max_stars_repo_name: runjak/hRubiks

max_issues_repo_name: runjak/hRubiks

max_forks_repo_name: runjak/hRubiks


text:
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE UndecidableInstances #-}

-- | This module defines commonly useful functions that are not specifically related with vectors and matrices.
module FPNLA.Utils (
    iif,
    mapPair,
    splitSlice,

    Truncable(..),

    prevEnum,
    nextEnum,
    enumerate

) where

import Data.Complex (Complex ((:+)))

-- | Takes two functions @f@ and @g@, and a pair, and applies @f@ to the first element and @g@ to the second element of the pair.
mapPair :: (a -> b) -- ^ @f@
    -> (c -> d) -- ^ @g@
    -> (a, c) -- ^ Input pair @(a,b)@
    -> (b, d) -- ^ Output pair @(f a, g b)@
mapPair f g (a, c) = (f a,  g c)

-- | It's equivalent to @if cond then e1 else e2@, but defined as a function.
iif :: Bool -- ^ @cond@
    -> a -- ^ @e1@
    -> a -- ^ @e2@
    -> a -- ^ @if cond then e1 else e2@
iif cond expTrue expFalse = if cond then expTrue else expFalse

-- | Splits a list into pieces of size @n@.
-- The last piece could be smaller if the length of the list is not multiple of @n@.
splitSlice :: Int -- ^ @n@
    -> [a] -- ^ The input list
    -> [[a]] -- ^ The input list divided in pieces of size @n@
splitSlice 0 _ = error "n must be greater than 0"
splitSlice _ [] = []
splitSlice n ls = sl : splitSlice n sls
    where (sl, sls) = splitAt n ls

-- | This class is used to compare two float point numbers.
-- With float point numbers, two different computations that theoretically returns the same value in fact returns different values. Even if the difference is so small that is ignorable, if we test using equality the test fail.
-- For this, instead to check for equality we define that two float point numbers are the same if the distance between them is less than some epsilon value.
class Truncable a where
    -- | Check if the distance between two float point values is less than some epsilon value.
    inEpsilonRange :: Real e => e -- ^ Epsilon
        -> a -- ^ A float point value
        -> a -- ^ Another float point value
        -> Bool

-- | We provide a default instance of 'Truncable' for any data type that has 'Floating' and 'RealFrac' instances.
instance (Floating n, RealFrac n) => Truncable n where
    inEpsilonRange epsilon n1 n2 = expN1 == expN2 && abs(mantN1 - mantN2) <= eps
        where
            eps = fromRational $ toRational epsilon
            shift n d1 = n * (10 ** fromIntegral d1)
            getExponent r
                    | r == 0 = 0
                    | r < 0 = getExponent (abs r)
                    | r > 10 = 1 + getExponent (r / 10)
                    | r < 1 = (-1) + getExponent (r * 10)
                    | otherwise = 0
            expN1 :: Integer = getExponent n1
            expN2 :: Integer = getExponent n2
            mantN1 = shift n1 (-expN1)
            mantN2 = shift n2 (-expN2)

-- | We define a specific instance of 'Truncable' for the 'Complex' data type.
instance (RealFrac n, Truncable n) => Truncable (Complex n) where
    inEpsilonRange epsilon (r1:+i1) (r2:+i2) = inEpsilonRange epsilon r1 r2 && inEpsilonRange epsilon i1 i2

-- | Returns the predecessor of a enumerated value.
-- The operation is cyclic so, the predecessor of the first value is the last value.
prevEnum :: (Eq a, Bounded a, Enum a) => a -> a
prevEnum a | a == minBound = maxBound
           | otherwise = pred a

-- | Returns the successor of a enumerated value.
-- The operation is cyclic so, the successor of the last value is the first value.
nextEnum :: (Eq a, Bounded a, Enum a) => a -> a
nextEnum a | a == maxBound = minBound
           | otherwise = succ a

-- | Returls a list containing all the values of an enumerated.
enumerate :: forall a. (Enum a, Bounded a) => [a]
enumerate = map toEnum [fromEnum (minBound :: a) .. fromEnum (maxBound :: a)]

max_stars_repo_name: mauroblanco/fpnla

max_issues_repo_name: mauroblanco/fpnla

max_forks_repo_name: mauroblanco/fpnla