text:
{-# LANGUAGE BangPatterns #-}
module MnistData (loadDataWrapper) where

import           Codec.Compression.GZip (decompress)
import qualified Data.ByteString.Lazy   as BL
import qualified Numeric.LinearAlgebra  as H

loadDataWrapper :: IO ([(H.Vector H.R, H.Vector H.R)], [(H.Vector H.R, Int)], [(H.Vector H.R, Int)])
loadDataWrapper = do
  (_, _, testImage) <- loadImage "data/t10k-images-idx3-ubyte.gz"
  (_, testLabel) <- loadLabel "data/t10k-labels-idx1-ubyte.gz"
  (_, _, trainImage) <- loadImage "data/train-images-idx3-ubyte.gz"
  (_, trainLabel) <- loadLabel "data/train-labels-idx1-ubyte.gz"
  let !testData = zip testImage testLabel
  let (!trainData', !validateDate) = splitAt 50000 $ zip trainImage trainLabel
  let trainData = map (\(a, b)-> (a, labelToVec b)) trainData'
  return (trainData, validateDate, testData)

labelToVec :: Int -> H.Vector H.R
labelToVec n = H.fromList $ replicate n 0 ++ [1.0] ++ replicate (9-n) 0

loadLabel :: FilePath -> IO (Int, [Int])
loadLabel fp = do
  content <- decompress <$> BL.readFile fp
  let (hd, dat) = BL.splitAt 8 content
  let (magic, size) = BL.splitAt 4 hd
  case BL.unpack magic of
    [0, 0, 8, 1] -> return (readInt size, take (readInt size) (parseData dat))
    _            -> return (0, [])

loadImage :: FilePath -> IO (Int, (Int, Int), [H.Vector H.R])
loadImage fp = do
  content <- decompress <$> BL.readFile fp
  let (hd, dat) = BL.splitAt 16 content
  let (magic, tl) = BL.splitAt 4 hd
  let (size, xy) =  BL.splitAt 4 tl
  let (rowW, colW) = BL.splitAt 4 xy
  let (row, col) = (readInt rowW, readInt colW)
  let images = map (H.fromList . map ((/256.0). fromIntegral)) (every (row * col) $ BL.unpack dat)
  case BL.unpack magic of
    [0, 0, 8, 3] -> return (readInt size, (row, col), images)
    _            -> return (0, (0, 0), images)

parseData :: BL.ByteString -> [Int]
parseData = map fromIntegral . BL.unpack

readInt :: Integral a => BL.ByteString -> a
readInt bs = foldl (\a b -> a * 256 + b) 0 $ map fromIntegral $ BL.unpack bs

every :: Int -> [a] -> [[a]]
every n xs = hd : every n tl
  where
    (hd, tl) = splitAt n xs

max_stars_repo_name: fujiisat/nnadl-haskell-study

max_issues_repo_name: fujiisat/nnadl-haskell-study

max_forks_repo_name: fujiisat/nnadl-haskell-study


text:
module SceneRandom where

import Numeric.Vector
import Numeric.Scalar (fromScalar, scalar)
import Control.Monad.State
import System.Random
import Types
import VecUtil

randDouble :: State SceneContext Double
randDouble = do
  ss <- get
  let (v, gen) = random (ss_getGen ss) :: (Double, StdGen)
  put $ ss { ss_getGen = gen }
  return v

randRange :: Double -> Double -> State SceneContext Double
randRange l h = randDouble >>= (\v -> return $ v * (h - l) + l)

randNeg1To1 :: State SceneContext Double
randNeg1To1 = randRange (-1) 1

randVecSphere :: State SceneContext Vec3d
randVecSphere = do
  x <- randNeg1To1
  y <- randNeg1To1
  z <- randNeg1To1
  return $ normalized $ vec3 x y z

randVecHemisphere :: Vec3d -> State SceneContext Vec3d
randVecHemisphere normal = do
  let (tangent, bitangent) = tangentBitangent normal

  theta <- randRange 0 (2 * pi)
  phi <- randRange 0 (pi / 2)

  let t = fromScalar $ scalar $ sin theta * cos phi
  let b = fromScalar $ scalar $ cos theta * cos phi
  let n = fromScalar $ scalar $ sin phi
  return $ tangent * t + bitangent * b + normal * n

randVecDisc :: Vec3d -> Vec3d -> Double -> State SceneContext Vec3d
randVecDisc normal center radius = do
  let (tangent, bitangent) = tangentBitangent normal
  theta <- randRange 0 (2 * pi)
  let t = fromScalar $ scalar $ sin theta * radius
  let b = fromScalar $ scalar $ cos theta * radius
  return $ tangent * t + bitangent * b + center

-- https://math.stackexchange.com/a/182936
randVecCone :: Double -> Vec3d -> State SceneContext Vec3d
randVecCone angle vector = do
  let (tangent, bitangent) = tangentBitangent vector
  z <- randRange (cos angle) 1
  phi <- randRange 0 (2 * pi)
  let theta = acos z
  let t = fromScalar $ scalar $ cos phi
  let b = fromScalar $ scalar $ sin phi
  let v = fromScalar $ scalar $ cos theta
  let sinTheta = fromScalar $ scalar $ sin theta
  return $ sinTheta * (t * tangent + b * tangent) + v * vector

randColor :: State SceneContext Vec3d
randColor = do
  r <- randDouble
  g <- randDouble
  b <- randDouble
  return $ vec3 r g b

ext: hs


text:
module PCA.GPExample where

import Universum hiding (transpose, Vector)

import Data.Vector.Unboxed as V (fromList)

import Data.Array.Repa
import Numeric.LinearAlgebra.Repa hiding (Matrix, Vector)

import PCA.GaussianProcess
import PCA.Types
import PCA.Util

testList :: [Double]
testList = [ -5.0
           , -4.79591837
           , -4.59183673
           , -4.3877551
           , -4.18367347
           , -3.97959184
           , -3.7755102
           , -3.57142857
           , -3.36734694
           , -3.16326531
           , -2.95918367
           , -2.75510204
           , -2.55102041
           , -2.34693878
           , -2.14285714
           , -1.93877551
           , -1.73469388
           , -1.53061224
           , -1.32653061
           , -1.12244898
           , -0.91836735
           , -0.71428571
           , -0.51020408
           , -0.30612245
           , -0.10204082
           ,  0.10204082
           ,  0.30612245
           ,  0.51020408
           ,  0.71428571
           ,  0.91836735
           ,  1.12244898
           ,  1.32653061
           ,  1.53061224
           ,  1.73469388
           ,  1.93877551
           ,  2.14285714
           ,  2.34693878
           ,  2.55102041
           ,  2.75510204
           ,  2.95918367
           ,  3.16326531
           ,  3.36734694
           ,  3.57142857
           ,  3.7755102
           ,  3.97959184
           ,  4.18367347
           ,  4.3877551
           ,  4.59183673
           ,  4.79591837
           ,  5.0
           ]

xTest :: Vector D Double
xTest = delay $ fromUnboxed (Z :. 50) (V.fromList testList)

inputObserve :: InputObservations Double
inputObserve = InputObservations xTest

computeSD
    :: Shape sh
    => Array D sh Double
    -> Array U sh Double
computeSD = computeS @D @_ @_ @U

sqDist
    :: Vector D Double
    -> Vector D Double
    -> Matrix D Double
sqDist vec1 vec2 = matrixSquare vec1 ++^ (transposeMatrix matrix2')
     where
         matrix2' = matrixSquare vec2 --^
                      ((delay . smap ((*) 2)) $
                      matrix1 `mulS` matrix2)
         matrixSquare matrix = smap flipExp matrix
         flipExp = (flip (^)) 2
         matrix1 = toMatrix' $ vec1
         matrix2 = transposeMatrix . toMatrix' $ vec2

kernelFunction
    :: Vector D Double
    -> Vector D Double
    -> Matrix D Double
kernelFunction matrix1 matrix2 =
    smap fun $ sqDist matrix1 matrix2
    where
        fun = exp . ((*) (-0.5)) . ((*) (1/0.1))

xTrainList :: [Double]
xTrainList =
    [ -4
    , -3
    , -2
    , -1
    , 1
    ]

xTrain :: Vector D Double
xTrain = delay $ fromUnboxed (Z :. 5) (V.fromList xTrainList)

yTrain :: Vector D Double
yTrain = smap sin xTrain

testTrainingData :: GPTrainingData Double
testTrainingData = GPTrainingData xTrain yTrain

k, kS, l, lK :: Matrix D Double
k = kernelFunction xTrain xTrain
kS = kernelFunction xTrain xTest
l = cholSH (k +^ (delay . smap (* 0.0005) $ identD 50))
lK = fromMaybe (error "izviniti") (delay <$> linearSolveS l kS)

testGaussianProcess :: GaussianProcess Double
testGaussianProcess = GaussianProcess kernelFunction

max_stars_repo_name: serokell/PCA

max_issues_repo_name: serokell/PCA

max_forks_repo_name: serokell/PCA


text:
{-|
Module      : Qubism.Algebra
Description : Typeclasses and utilities for Vector and Hilbert spaces.
Copyright   : (c) Keith Pearson, 2018
License     : MIT
Maintainer  : keith@qubitrot.org
-}

module Qubism.Algebra where

import Data.Complex
import Data.Monoid

type C = Complex Double

-- | A vector space over the complex numbers
class Eq v => VectorSpace v where
  -- | The additive identity
  zero :: v
  -- | Scalar multiplication
  (.:) :: C -> v -> v
  -- | Vector addition
  (+:) :: v -> v -> v
  -- | Vector subtraction
  (-:) :: v -> v -> v
  -- | Inverse element
  neg  :: v -> v
  a -: b = a +: (neg  b)

-- | A Hilbert space over the complex numbers
class VectorSpace v => HilbertSpace v where
  -- | Sesquilinear inner product
  (<.>) :: v -> v -> C
  -- | Norm induced by the inner product
  norm  :: v -> Double
  norm a = realPart $ a <.> a

class VectorSpace v => Algebra v where
  -- | Bilinear product
  (*:) :: v -> v -> v

commutator :: Algebra v  => v -> v -> v
commutator a b = (a *: b) -: (b *: a)

anticommutator :: Algebra v => v -> v -> v
anticommutator a b = (a *: b) +: (b *: a)

max_stars_repo_name: qubitrot/qubism

max_issues_repo_name: qubitrot/qubism

max_forks_repo_name: qubitrot/qubism


text:
{-# OPTIONS -Wall #-}

module Parsing
  ( parseExpr
  , readExpr
  , readExprList
  ) where

import Text.ParserCombinators.Parsec hiding (spaces)

import Control.Monad.Except
import Data.Array.Base (listArray)
import Data.Char (toLower, toUpper)
import Data.Complex
import Data.Ratio
import Numeric

import Types

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

spaces :: Parser ()
spaces = skipMany1 space

parseAtom :: Parser LispVal
parseAtom = do
  first <- letter <|> symbol
  rest <- many (letter <|> digit <|> symbol)
  let atom = first : rest
  return $
    case atom of
      "#t" -> Bool True
      "#f" -> Bool False
      _ -> Atom atom

parseList :: Parser LispVal
parseList = List <$> sepBy parseExpr spaces

parseDottedList :: Parser LispVal
parseDottedList = do
  x <- endBy parseExpr spaces
  y <- char '.' >> spaces >> parseExpr
  return $ DottedList x y

parseAnyList :: Parser LispVal
parseAnyList = do
  _ <- char '('
  x <- try parseList <|> parseDottedList
  _ <- char ')'
  return x

parseQuoted :: Parser LispVal
parseQuoted = do
  _ <- char '\''
  x <- parseExpr
  return $ List [Atom "quote", x]

parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
  _ <- char '`'
  x <- parseExpr
  return $ List [Atom "quasiquote", x]

parseUnquoted :: Parser LispVal
parseUnquoted = do
  _ <- char ','
  x <- parseExpr
  return $ List [Atom "unquote", x]

parseDecimal :: Parser LispVal
parseDecimal = parseDecimal1 <|> parseDecimal2

parseDecimal1 :: Parser LispVal
parseDecimal1 = Number . read <$> many1 digit

parseDecimal2 :: Parser LispVal
parseDecimal2 = do
  _ <- try $ string "#d"
  x <- many1 digit
  (return . Number . read) x

parseOctal :: Parser LispVal
parseOctal = do
  _ <- try $ string "#o"
  x <- many1 digit
  (return . Number . fst . (!! 0) . readOct) x

parseHexadecimal :: Parser LispVal
parseHexadecimal = do
  _ <- try $ string "#x"
  x <- many1 digit
  (return . Number . fst . (!! 0) . readHex) x

parseBinary :: Parser LispVal
parseBinary = do
  _ <- try $ string "#b"
  x <- many1 (oneOf "01")
  (return . Number . readBin) x

readBin :: String -> Integer
readBin = readBin' 0

readBin' :: Integer -> String -> Integer
readBin' c "" = c
readBin' c (x:xs) = readBin' cc xs
  where
    cc =
      2 * c +
      if x == '0'
        then 0
        else 1

parseNumber :: Parser LispVal
parseNumber = parseDecimal <|> parseHexadecimal <|> parseOctal <|> parseBinary

escape :: Parser String
escape = do
  d <- char '\\'
  c <- oneOf "\\\"0nrvtbf"
  return [d, c]

nonEscape :: Parser Char
nonEscape = noneOf "\\\"\0\n\r\v\t\b\f"

character :: Parser String
character = fmap return nonEscape <|> escape

parseString :: Parser LispVal
parseString = do
  _ <- char '"'
  xs <- many character
  _ <- char '"'
  return $ String $ concat xs

parseBool :: Parser LispVal
parseBool = do
  _ <- char '#'
  (char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False))

caseInsensitiveChar :: Char -> Parser Char
caseInsensitiveChar c = char (toLower c) <|> char (toUpper c)

caseInsensitiveString :: String -> Parser String
caseInsensitiveString s = try (mapM caseInsensitiveChar s) <?> "\"" ++ s ++ "\""

parseCharacter :: Parser LispVal
parseCharacter = do
  _ <- try $ string "#\\"
  val <-
    try (caseInsensitiveString "newline" <|> caseInsensitiveString "space") <|> do
      x <- anyChar
      notFollowedBy alphaNum
      return [x]
  return $
    Character $
    case val of
      "newline" -> '\n'
      "space" -> ' '
      _ -> head val

parseFloat :: Parser LispVal
parseFloat = do
  x <- many1 digit
  _ <- char '.'
  y <- many1 digit
  (return . Float . fst . head . readFloat) (x ++ "." ++ y)

parseRatio :: Parser LispVal
parseRatio = do
  x <- many1 digit
  _ <- char '/'
  y <- many1 digit
  return $ Ratio (read x % read y)

toDouble :: LispVal -> Double
toDouble (Float f) = realToFrac f
toDouble (Number n) = fromIntegral n
toDouble _ = undefined

parseComplex :: Parser LispVal
parseComplex = do
  x <- try parseFloat <|> parseDecimal
  _ <- char '+'
  y <- try parseFloat <|> parseDecimal
  _ <- char 'i'
  return $ Complex (toDouble x :+ toDouble y)

parseVector :: Parser LispVal
parseVector = do
  _ <- string "#("
  x <- parseVector'
  _ <- char ')'
  return x

parseVector' :: Parser LispVal
parseVector' = do
  values <- sepBy parseExpr spaces
  return $ Vector (listArray (0, length values - 1) values)

parseExpr :: Parser LispVal
parseExpr =
  parseAtom <|> parseString <|> try parseFloat <|> try parseRatio <|> try parseComplex <|> try parseNumber <|>
  parseQuoted <|>
  parseQuasiQuoted <|>
  parseUnquoted <|>
  try parseVector <|>
  parseAnyList <|>
  try parseBool <|>
  try parseCharacter

readOrThrow :: Parser a -> String -> ThrowsError a
readOrThrow parser input =
  case parse parser "lisp" input of
    Left err -> throwError $ Parser err
    Right x -> return x

readExpr :: String -> ThrowsError LispVal
readExpr = readOrThrow parseExpr

readExprList :: String -> ThrowsError [LispVal]
readExprList = readOrThrow (endBy parseExpr spaces)

max_stars_repo_name: willGuimont/write-yourself-a-scheme

max_issues_repo_name: willGuimont/write-yourself-a-scheme

max_forks_repo_name: willGuimont/write-yourself-a-scheme


text:
{-# LANGUAGE FlexibleContexts #-}

module AI.Util.Matrix where

import Control.Monad.Random hiding (fromList)
import Foreign.Storable (Storable)
import Numeric.LinearAlgebra

-- |Return the size of a matrix as a 2-tuple.
size :: Matrix a -> (Int,Int)
size x = (rows x, cols x)

-- |Concatenate matrices horizontally.
horzcat :: Element a => [Matrix a] -> Matrix a
horzcat = fromBlocks . return

-- |Concatenate matrices vertically.
vertcat :: Element a => [Matrix a] -> Matrix a
vertcat = fromBlocks . map return

-- |Add a column of ones to a matrix.
addOnes :: Matrix Double -> Matrix Double
addOnes x = fromBlocks [[1, x]]

-- |Create a row matrix.
row :: [Double] -> Matrix Double
row = asRow . fromList

-- |Create a column matrix.
column :: [Double] -> Matrix Double
column = asColumn . fromList

--------------------------
-- Functions on Vectors --
--------------------------

takeVector :: Storable a => Int -> Vector a -> Vector a
takeVector n v = subVector 0 n v

dropVector :: Storable a => Int -> Vector a -> Vector a
dropVector n v = subVector n (dim v - n) v

sumVector :: (Num a, Storable a) => Vector a -> a
sumVector xs = foldVector (+) 0 xs

prodVector :: (Num a, Storable a) => Vector a -> a
prodVector xs = foldVector (*) 1 xs

---------------------------
-- Functions on Matrices --
---------------------------

mapRows :: Element a => (Vector a -> b) -> Matrix a -> [b]
mapRows f m = map f (toRows m)

mapCols :: Element a => (Vector a -> b) -> Matrix a -> [b]
mapCols f m = map f (toColumns m)

eachRow :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
eachRow f = fromRows . mapRows f

eachCol :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
eachCol f = fromColumns . mapCols f

sumRows :: (Element a, Num (Vector a)) => Matrix a -> Vector a
sumRows m = sum $ toRows m

sumCols :: (Element a, Num (Vector a)) => Matrix a -> Vector a
sumCols m = sum $ toColumns m

sumMatrix :: Matrix Double -> Double
sumMatrix = sumVector . sum . toRows

------------------------
-- Subset Referencing --
------------------------

subRefVec :: Storable a => Vector a -> [Int] -> Vector a
subRefVec v is = fromList $ map (v@>) is

subRefRows :: Element a => Matrix a -> [Int] -> Matrix a
subRefRows m is = fromRows $ map (r!!) is where r = toRows m

subRefCols :: Element a => Matrix a -> [Int] -> Matrix a
subRefCols m is = fromColumns $ map (c!!) is where c = toColumns m

max_stars_repo_name: cagix/aima-haskell


text:
#!/usr/bin/env stack
-- stack --resolver lts-15.04 runghc --package reanimate
module Main where

import           Codec.Picture.Types
import           Control.Lens
import           Control.Monad
import           Data.List
import           Data.Maybe
import qualified Data.Text                     as T
import           Data.Tuple
import qualified Data.Vector                   as V
import           Linear.V2
import           Linear.Vector
import qualified Numeric.LinearAlgebra         as Matrix
import           Numeric.LinearAlgebra.HMatrix (Matrix, linearSolve, toLists,
                                                (><))
import           Reanimate
import           Reanimate.Math.Common         (barycentricCoords, isBetween,
                                                rayIntersect)

type Points = V.Vector (V2 Double)
type Edges = [(Int, Int, Int)]
data Mesh = Mesh { meshPoints :: Points, meshEdges :: Edges }
data MeshPair = MeshPair Points Points Edges
-- The points in a RelMesh are:
--   relMeshStatic ++ x where Ax = B
data RelMesh = RelMesh
  { relMeshStatic :: Points
  , relMeshEdges  :: Edges
  , relMeshA      :: Matrix Double
  , relMeshB      :: Matrix Double
  }
data RelMeshPair = RelMeshPair Points Edges (Matrix Double) (Matrix Double) (Matrix Double) (Matrix Double)
-- Linear interpolation on RelMesh gives smooth morph.
-- solveMesh :: RelMesh -> Mesh
-- mkRelative :: Mesh -> RelMesh
-- mkRelativePair :: MeshPair -> RelMeshPair
-- triangulate :: Polygon -> Polygon -> MeshPair ?
-- embed :: MeshPair -> MeshPair
-- compatible :: Mesh -> Mesh -> Maybe MeshPair
-- linearInterpolate :: MeshPair -> Double -> Mesh
-- convexInterpolate :: RelMeshPair -> Double -> RelMesh

main :: IO ()
main = reanimate morphAnimation

morphAnimation :: Animation
morphAnimation = playThenReverseA $ pauseAround 1 1 $ addStatic (mkBackground "black") $
  signalA (curveS 2) $ mkAnimation 5 $ \t -> lowerTransformations $ scale 3 $ pathify $ center $
  mkGroup
  [ translate (-1) 0 $ drawTrig (linearInterpolate meshPair t)
  , translate 1 0 $ drawTrig $ solveMesh $ convexInterpolate relPair t
  ]
  where
    relPair = mkRelativePair meshPair
    meshPair = fromJust $ compatible example1 example2

mkLineP :: P -> P -> SVG
mkLineP (V2 x1 y1) (V2 x2 y2) = mkLine (x1,y1) (x2,y2)

-- FIXME: Check that the triangles are all anticlockwise.
-- FIXME: Check that the edges connect all the points.
-- FIXME: Check that the edgse leave no gaps.
compatible :: Mesh -> Mesh -> Maybe MeshPair
compatible a b =
  if meshEdges a == meshEdges b && V.length (meshPoints a) == V.length (meshPoints b)
    then Just $ MeshPair (meshPoints a) (meshPoints b) (meshEdges a)
    else Nothing

linearInterpolate :: MeshPair -> Double -> Mesh
linearInterpolate (MeshPair aP bP edges) t = Mesh
    { meshPoints = V.zipWith (lerp (1-t)) aP bP
    , meshEdges = edges }

example1 :: Mesh
example1 = Mesh points edges
  where
    points = V.fromList
      [ V2 1 0
      , V2 (-1/2) (sqrt 3 / 2)
      , V2 (-1/2) (-sqrt 3 / 2)
      , 3 * points V.! 0 ^/ 4
      , 3 * points V.! 1 ^/ 4
      , 3 * points V.! 2 ^/ 4
      , points V.! 0 ^/ 2
      , points V.! 1 ^/ 2
      , points V.! 2 ^/ 2
      ]
    edges =
      [ (1,5,4), (1,2,5), (2,6,5), (2,3,6), (3,4,6), (3,1,4)
      , (4,8,7), (4,5,8), (5,9,8), (5,6,9), (6,7,9), (6,4,7), (7,8,9)]

example2 :: Mesh
example2 = Mesh points edges
  where
    points = V.fromList
      [ V2 1 0
      , V2 (-1/2) (sqrt 3 / 2)
      , V2 (-1/2) (-sqrt 3 / 2)
      , 3 * points V.! 2 ^/ 4
      , 3 * points V.! 0 ^/ 4
      , 3 * points V.! 1 ^/ 4
      , points V.! 1  ^/ 2
      , points V.! 2 ^/ 2
      , points V.! 0 ^/ 2
      ]
    edges = meshEdges example1

drawTrig (Mesh points gs) = withStrokeColor "grey" $ withFillColor "white" $
  withStrokeWidth (defaultStrokeWidth/2) $ mkGroup
  [ mkGroup
    [ mkGroup
      [ mkLine (ax, ay) (bx, by)
      , mkLine (bx, by) (cx, cy)
      , mkLine (cx, cy) (ax, ay)
      ]
    | (a, b, c) <- gs
    , let V2 ax ay = points V.! (a-1)
          V2 bx by = points V.! (b-1)
          V2 cx cy = points V.! (c-1)
    ]
  , mkGroup $ concat
    [ [ colored v $ translate ax ay $ mkCircle circleRadius
      , withStrokeWidth 0 $
        withStrokeColor "white" $ withFillColor "black" $ mkGroup
        [ translate ax ay $ ppNum v ]
      ]
    | v <- nub $ concat [ [a,b,c] | (a,b,c) <- gs]
    , let V2 ax ay = points V.! (v-1)
    ]]
  where
    colored n =
      let c = promotePixel $ turbo (fromIntegral n / fromIntegral (length gs-1))
      in withStrokeColorPixel c . withFillColorPixel c
    ppNum n = cachedNumbers !! n
     --scaleToHeight (circR*1.5) $ center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}"

cachedNumbers =
  [ scaleToHeight (circleRadius*1.5) $ center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}"
  | n <- [0 .. ] ]

circleRadius :: Double
circleRadius = 0.05

-- Anticlockwise. No duplicate vertices. length >= 3
type Polygon = [V2 Double]
type P = V2 Double

-- T = (U, G)
-- G = [Polygon]
-- U = nub $ concat G

findStarNeighbours :: Eq a => [(a,a,a)] -> a -> [(a, a)]
findStarNeighbours allTrig self =
  [ (b,c)
  | (a,b,c) <- allTrig
  , self == a
  ] ++
  [ (c,a)
  | (a,b,c) <- allTrig
  , self == b
  ] ++
  [ (a,b)
  | (a,b,c) <- allTrig
  , self == c
  ]

isInterior :: Eq a => [(a, a)] -> Bool
isInterior = isJust . getExteriorPoly

getExteriorPoly :: Eq a => [(a, a)] -> Maybe [a]
getExteriorPoly [] = Nothing
getExteriorPoly ((a,b):rest) = worker [a] a b rest
  where
    worker acc start this [] = do
      guard (start == this)
      return (reverse acc)
    worker acc start this xs =
      case lookup this xs of
        Just next -> worker (this:acc) start next (delete (this,next) xs)
        Nothing   ->
          case lookup this (map swap xs) of
            Just next -> worker (this:acc) start next (delete (next, this) xs)
            Nothing   -> Nothing

convexInterpolate :: RelMeshPair -> Double -> RelMesh
convexInterpolate (RelMeshPair static edges leftM leftB rightM rightB) t =
  RelMesh
  { relMeshStatic = static
  , relMeshEdges  = edges
  , relMeshA      = Matrix.scale (1-t) leftM +
                    Matrix.scale t rightM
  , relMeshB      = Matrix.scale (1-t) leftB +
                    Matrix.scale t rightB
  }

solveMesh :: RelMesh -> Mesh
solveMesh (RelMesh static edges m b) =
  case linearSolve m b of
    Nothing -> error "Failed to solve mesh"
    Just ret ->
      Mesh (static <> V.fromList (worker (toLists ret))) edges
  where
    worker []             = []
    worker ([x]:[y]:rest) = V2 x y : worker rest
    worker _              = error "invalid result"

mkRelative :: Mesh -> RelMesh
mkRelative (Mesh points edges) = RelMesh (V.fromList exteriorPoints) edges mM bM
  where
    mM = (s><s) (concat m)
    bM = (s><1) b
    (s,exterior, (m, b)) = toParameters points edges
    exteriorPoints =
      [ points V.! (i-1)
      | i <- exterior
      ]

mkRelativePair :: MeshPair -> RelMeshPair
mkRelativePair (MeshPair p1 p2 edges) =
  let RelMesh static _ leftM leftB = mkRelative (Mesh p1 edges)
      RelMesh _ _ rightM rightB = mkRelative (Mesh p2 edges)
  in RelMeshPair static edges leftM leftB rightM rightB

toParameters points groups = (length interior*2,exterior,unzip $ concat
  [ let lst = [(if i == j then -1 else t)
              | j <- interior
              , let t = fromMaybe 0 $ lookup (i,j) lam_ij_cache
              ]
        pos = negate $ sum
           [ pj ^* t
           | j <- exterior
           , let t = fromMaybe 0 $ lookup (i,j) lam_ij_cache
                 pj = points V.! (j-1)
           ]
    in [ (dupX lst, pos ^. _x)
       , (dupY lst, pos ^. _y)]
  | i <- interior ])
  where
    lam_ij_cache = lam_ij points groups
    dupX []     = []
    dupX (x:xs) = x:0:dupX xs
    dupY []     = []
    dupY (x:xs) = 0:x:dupY xs
    (interior, exterior) =
      partition (isInterior . findStarNeighbours groups) [1 .. length points]

lam_ij points groups =
  [ ((i, j), t)
  | i <- [1..V.length points]
  , (j, t) <- lam_j points groups i
  ]
lam_j points groups p =
  [ (nP, sum [ t | (j,k,t) <- mu, j == nP ] / fromIntegral (length nPoints))
  | let n = findStarNeighbours groups p
        nPoints = fromMaybe [] $ getExteriorPoly n
        mu = calcMu points groups p
  , nP <- nPoints ]
calcMu points groups p = concat
    [ [ (nP, nP, t1)
      , (a, nP, t2)
      , (b, nP, t3) ]
      -- (nP, a, b)
    | {-p <- [1..length points]-}
      let selfVert = points V.! (p-1)
          n = findStarNeighbours groups p
          nPoints :: [Int]
          nPoints = fromMaybe [] $ getExteriorPoly n
    , (i, nP) <- zip [1 .. ] nPoints
    , let vert = points V.! (nP-1)
    , let line = (vert, selfVert)
    , let (a,b,aP,bP) = head $
            [ (a,b,aP,bP)
            | (a,b) <- n
            , let aP = points V.! (a-1)
                  bP = points V.! (b-1)
                  segment = (points V.! (a-1), points V.! (b-1))
            , case rayIntersect line segment of
                Nothing -> False
                Just u  -> isBetween u segment
            , a /= nP
            , b /= nP ]
    -- , b == (nPoints ++ nPoints) !! i
    , let (t1,t2,t3) = barycentricCoords vert aP bP selfVert
    ]

max_stars_repo_name: TristanCacqueray/reanimate

max_issues_repo_name: TristanCacqueray/reanimate

max_forks_repo_name: TristanCacqueray/reanimate


text:
module Main where

import MNIST.Prelude
import Control.Monad.Random
import qualified Numeric.LinearAlgebra.Static as SA

import Grenade

type MNIST
  = Network
    '[ FullyConnected 784 300, Logit
     , FullyConnected 300 100, Logit
     , FullyConnected 100   9, Logit
     ]
    '[ 'D1 784, 'D1 300
     , 'D1 300, 'D1 100
     , 'D1 100, 'D1 9
     , 'D1 9
     ]

randomNet :: MonadRandom m => m MNIST
randomNet = randomNetwork

type Input  = S ('D1 784)
type Output = S ('D1   9)
type DataSet = [(Input, Output)]

trainAll :: LearningParameters -> MNIST -> DataSet -> MNIST
trainAll rate net0 dataset = foldl' trainOne net0 dataset
  where
  trainOne :: MNIST -> (Input, Output) -> MNIST
  trainOne !network (i,o) = train rate network i o

netTrain :: MNIST -> LearningParameters -> Int -> IO MNIST
netTrain net0 rate n = do
  inps <- replicateM n randomTrainingData
  let outs = map randomTrainingLabels inps
  return $ trainAll rate net0 (zip inps outs)
  where
    randomTrainingData :: IO Input
    randomTrainingData = do
      s <- getRandom
      return . S1D $ SA.randomVector s SA.Uniform * 2 - 1

    randomTrainingLabels :: Input -> Output
    randomTrainingLabels (S1D v) = S1D . fromIntegral $ fromEnum isTrue
      where
        isTrue :: Bool
        isTrue = (v `inCircle` (fromRational    0.33, 0.33))
              || (v `inCircle` (fromRational (-0.33), 0.33))

    inCircle :: KnownNat n => R n -> (R n, Double) -> Bool
    inCircle v (o, r) = SA.norm_2 (v - o) <= r

netScore :: MNIST -> IO ()
netScore network = putStrLn . unlines $
  (fmap.fmap) (showNorm . go) testIns

  where
    testIns :: [[(Double, Double)]]
    testIns = [ [ (x,y)  | x <- [0..50] ] | y <- [0..20] ]

    go :: (Double, Double) -> Output
    go (x,y) = runNet network (S1D $ SA.vector [x / 25 - 1, y / 10 - 1])

    showNorm :: Output -> Char
    showNorm (S1D r) = render $ SA.mean r

    render :: Double -> Char
    render n | n <= 0.2  = ' '
             | n <= 0.4  = '.'
             | n <= 0.6  = '-'
             | n <= 0.8  = '='
             | otherwise = '#'

main :: IO ()
main = do
  net0 <- randomNet
  net  <- netTrain net0 params examples
  netScore net
  where
    examples :: Int
    examples = 10000

    params :: LearningParameters
    params = LearningParameters
      { learningRate = 0.01
      , learningMomentum = 0.9
      , learningRegulariser = 0.0005
      }

max_stars_repo_name: stites/haskell-mnist-benchmarks

max_issues_repo_name: stites/haskell-mnist-benchmarks

max_forks_repo_name: stites/haskell-mnist-benchmarks


text:
module Phi.Window
  ( Window, mkWindow
  , length
  , push
  , phi
  ) where

import           Data.Vector.Unboxed            (Vector)
import qualified Data.Vector.Unboxed            as V
import           Prelude                        hiding (length)
import           Statistics.Distribution        (complCumulative)
import           Statistics.Distribution.Normal (normalFromSample)

data Window =
  Window
  { sample :: !(Vector Double)
  , size   :: !Int
  }

mkWindow :: Int -> Window
mkWindow = Window V.empty

push :: Double -> Window -> Window
push d window = window {sample = sample'}
  where sample' =
          V.take (size window)
                 (V.cons d (sample window))

phi :: Window -> Double -> Double
phi w = negate . logBase 10 . pLater w
  where pLater :: Window -> Double -> Double
        pLater = complCumulative . normalFromSample . sample

length :: Window -> Int
length = V.length . sample

max_stars_repo_name: reinh/hs-phi

max_issues_repo_name: reinh/hs-phi

max_forks_repo_name: reinh/hs-phi


text:
{-# LANGUAGE ConstraintKinds        #-}
{-# LANGUAGE DataKinds              #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE KindSignatures         #-}
{-# LANGUAGE RecordWildCards        #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE TypeApplications       #-}
{-# LANGUAGE TypeOperators          #-}

module Numeric.Optimization.TwoPhase.Tableau.Build
    ( IsBuilder
    , mkTable
    ) where

import           Prelude hiding                             ((++))

import qualified Data.List as List
import           Data.Maybe                                 (fromJust)
import qualified Data.Vector.Sized as Vec
import           Data.Vector.Storable.Sized                 (Vector, (++))
import qualified Data.Vector.Storable.Sized as SVec
import           GHC.TypeLits
import           Numeric.LinearAlgebra.Static               (L)
import qualified Numeric.LinearAlgebra.Static as LS
import qualified Numeric.LinearAlgebra.Static.Vector as LS

import           Numeric.Optimization.Problem
import           Numeric.Optimization.TwoPhase.Types

type IsBuilder v s a c =
    ( KnownNat v
    , KnownNat s
    , KnownNat a
    , KnownNat c
    )

data Builder (v :: Nat) (s :: Nat) (a :: Nat) (c :: Nat) = Builder
    { slackSupply       :: [Vector s Double]
    , artificialSupply  :: [Vector a Double]
    , rowsBuilt         :: [Vector (Cols 'PhaseI v s a) Double]
    , adjustmentsBuilt  :: [Vector (Cols 'PhaseI v s a) Double]
    }

mkBuilder :: (IsBuilder v s a c) => Builder v s a c
mkBuilder = Builder slackI artificialI [] []
  where
    slackI      = reverse . fmap LS.rVec . LS.toRows $ LS.eye
    artificialI = reverse . fmap LS.rVec . LS.toRows $ LS.eye

buildPhaseI
    :: forall v s a c.
        (IsBuilder v s a c)
    => Builder v s a c
    -> Builder v s a c
buildPhaseI acc = acc
    { rowsBuilt = row : rowsBuilt acc }
  where
    row = SVec.snoc (os ++ xs ++ ss ++ as) 0
    os  = SVec.fromTuple (1, 0)
    xs  = SVec.replicate @v 0
    ss  = SVec.replicate @s 0
    as  = SVec.replicate @a (-1)

buildPhaseII
    :: forall v s a c.
        (IsBuilder v s a c)
    => Builder v s a c
    -> Coeffs v
    -> Builder v s a c
buildPhaseII acc xs' = acc
    { rowsBuilt = row : rowsBuilt acc }
  where
    row = SVec.snoc (os ++ xs ++ ss ++ as) 0
    os  = SVec.fromTuple (0, 1)
    xs  = SVec.map negate xs'
    ss  = SVec.replicate @s 0
    as  = SVec.replicate @a 0

buildObjective :: (IsBuilder v s a c)
    => Builder v s a c -> Coeffs v -> Builder v s a c
buildObjective acc = buildPhaseI . buildPhaseII acc

buildLEQ
    :: forall v s a c.
        (IsBuilder v s a c)
    => Builder v s a c
    -> Coeffs v
    -> Double
    -> Builder v s a c
buildLEQ acc xs rhs = acc
    { slackSupply   = tail (slackSupply acc)
    , rowsBuilt     = row : rowsBuilt acc
    }
  where
    row = SVec.snoc (os ++ xs ++ ss ++ as) rhs
    os  = SVec.fromTuple (0, 0)
    ss  = head (slackSupply acc)
    as  = SVec.replicate @a 0

buildGEQ
    :: forall v s a c.
        (IsBuilder v s a c)
    => Builder v s a c
    -> Coeffs v
    -> Double
    -> Builder v s a c
buildGEQ acc xs' rhs' = acc
    { slackSupply       = tail (slackSupply acc)
    , artificialSupply  = tail (artificialSupply acc)
    , rowsBuilt         = row : rowsBuilt acc
    }
  where
    row = SVec.snoc (os ++ xs ++ ss ++ as) rhs
    os  = SVec.fromTuple (0, 0)
    xs  = SVec.map negate xs'
    ss  = head (slackSupply acc)
    as  = head (artificialSupply acc)
    rhs = negate rhs'

buildEQU
    :: forall v s a c.
        (IsBuilder v s a c)
    => Builder v s a c
    -> Coeffs v
    -> Double
    -> Builder v s a c
buildEQU acc xs rhs = acc
    { artificialSupply  = tail (artificialSupply acc)
    , rowsBuilt         = row : rowsBuilt acc
    , adjustmentsBuilt  = row : adjustmentsBuilt acc
    }
  where
    row = SVec.snoc (os ++ xs ++ ss ++ as) rhs
    os  = SVec.fromTuple (0, 0)
    ss  = SVec.replicate @s 0
    as  = head (artificialSupply acc)

buildConstraint :: (IsBuilder v s a c)
    => Builder v s a c -> Constraint v x y -> Builder v s a c
buildConstraint acc (xs :< y) = buildLEQ acc xs y
buildConstraint acc (xs :> y) = buildGEQ acc xs y
buildConstraint acc (xs := y) = buildEQU acc xs y

toMatrix :: forall v s a c. (IsBuilder v s a c)
    => Builder v s a c -> L (Rows 'PhaseI c) (Cols 'PhaseI v s a)
toMatrix acc =
    LS.rowsL . fromJust . Vec.fromList $ fmap LS.vecR adjusted
  where
    rows    = rowsBuilt acc
    adjusts = adjustmentsBuilt acc

    adjusted = case adjusts of
        [] -> rows
        xs ->
            let oldObjP1  = head rows
                newObjP1  = List.foldl1' (SVec.zipWith (+)) (oldObjP1 : xs)
            in  newObjP1 : tail rows

mkTable :: forall d v s a c. (IsBuilder v s a c)
    => Problem d v s a c -> L (Rows 'PhaseI c) (Cols 'PhaseI v s a)
mkTable = toMatrix . go mkBuilder
  where
    go :: Builder v s a c -> Problem d v x y z -> Builder v s a c
    go acc (Maximize xs)    = buildObjective acc xs
    go acc (Minimize xs)    = buildObjective acc xs
    go acc (SuchThat p c)   = go (buildConstraint acc c) p

max_stars_repo_name: alex-mckenna/optimum

max_issues_repo_name: alex-mckenna/optimum

max_forks_repo_name: alex-mckenna/optimum


text:
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleInstances #-}

module Math.Matrix
       ( DimVector (..)
       , DimMatrix (..)
       , GPConstraint
       , withMat
       , toMatrix
       , mulM
       , transposeM
       , mapMM
       , mapDiagonalM
       , invSM
       , substractMeanM
       , (+^^)
       , (-^^)
       , (*^^)
       , cholM
       , foldAllSM
       , sumAllSM
       , trace2SM
       , toDimMatrix
       , detSM
       , toMatrixM'
       , zipWithDim
       , zipWithArray
       , identM
       , delayMatrix
       , linearSolveM
       , vectorLength
       , randomMatrixD
       , matrixRowsNum
       , matrixColsNum
       , pinvSM
       , mulPM
       , eigSHM
       )
       where

import Universum hiding (Vector, transpose, map, natVal,
                         zipWith)

import           GHC.TypeLits hiding (someNatVal)

import           Data.Array.Repa
import           Data.Random.Normal (normals)
import           Data.Vector.Unboxed.Base (Unbox)
import           Data.Vinyl.TypeLevel (AllConstrained)
--import           Numeric.Dimensions
import           Numeric.LinearAlgebra.Repa hiding (Matrix, Vector)
import           System.Random (Random, RandomGen)

import PCA.Types
import PCA.Util hiding (randomMatrixD, zipWithArray)

newtype DimVector r (n :: Nat) a
  = DimVector { runVector :: Vector r a }
  deriving Eq

newtype DimMatrix r (y :: Nat) (x :: Nat) a
  = DimMatrix { getInternal :: Matrix r a}

instance Functor (DimVector D n) where
  fmap f (DimVector vector) = DimVector (smap f vector)

withMat
  :: Matrix D a
  -> (forall x y. (KnownNat x, KnownNat y) => DimMatrix D x y a -> k)
  -> k
withMat m f =
  let (Z :. x :. y) = extent m
  in
  case someNatVal (fromIntegral x) of
    SomeNat (Proxy :: Proxy m) -> case someNatVal (fromIntegral y) of
      SomeNat (Proxy :: Proxy n) -> f (DimMatrix @_ @m @n m)

mulM
  :: forall y1 x1 y2 x2 a.
  ( AllConstrained KnownNat [x1, x2, y1, y2]
  , Numeric a
  , x1 ~ y2
  )
  => DimMatrix D y1 x1 a
  -> DimMatrix D y2 x2 a
  -> DimMatrix D y1 x2 a
mulM (DimMatrix m1) (DimMatrix m2) = DimMatrix $ delay $ m1 `mulS` m2

transposeM
  :: (KnownNat y, KnownNat x)
  => DimMatrix D y x a
  -> DimMatrix D x y a
transposeM (DimMatrix m) = DimMatrix $ transpose m

mapMM
  ::
  ( KnownNat y
  , KnownNat x
  , Unbox a
  , Unbox b
  )
  => (a -> b)
  -> DimMatrix D y x a
  -> DimMatrix D y x b
mapMM f (DimMatrix m) =  DimMatrix $ map f m

mapDiagonalM
  ::
  ( KnownNat y
  , KnownNat x
  , Unbox a
  )
  => (a -> a)
  -> DimMatrix D y x a
  -> DimMatrix D y x a
mapDiagonalM f (DimMatrix m) = DimMatrix $ mapDiagonal f m

invSM
  ::
  ( KnownNat y
  , KnownNat x
  , Field a
  , Numeric a
  , y ~ x
  )
  => DimMatrix D y x a
  -> DimMatrix D y x a
invSM (DimMatrix m) = DimMatrix $ delay $ invS m

substractMeanM
  ::
  ( KnownNat y
  , KnownNat x
  )
  => DimMatrix D y x Double
  -> DimMatrix D y x Double
substractMeanM (DimMatrix m) = DimMatrix $ substractMean m

infixl 6 +^^, -^^
infixl 7 *^^

(+^^)
  :: forall y1 x1 y2 x2 a.
  ( AllConstrained KnownNat [x1, x2, y1, y2]
  , x1 ~ x2
  , y1 ~ y2
  , Num a
  )
  => DimMatrix D y1 x1 a
  -> DimMatrix D y2 x2 a
  -> DimMatrix D y2 x2 a
(+^^) (DimMatrix m1) (DimMatrix m2) = DimMatrix $ m1 +^ m2

(-^^)
  :: forall y1 x1 y2 x2 a.
  ( AllConstrained KnownNat [x1, x2, y1, y2]
  , x1 ~ x2
  , y1 ~ y2
  , Num a
  )
  => DimMatrix D y1 x1 a
  -> DimMatrix D y2 x2 a
  -> DimMatrix D y2 x2 a
(-^^) (DimMatrix m1) (DimMatrix m2) = DimMatrix $ m1 -^ m2

(*^^)
  :: forall y1 x1 y2 x2 a.
  ( AllConstrained KnownNat [x1, x2, y1, y2]
  , x1 ~ x2
  , y1 ~ y2
  , Num a
  )
  => DimMatrix D y1 x1 a
  -> DimMatrix D y2 x2 a
  -> DimMatrix D y2 x2 a
(*^^) (DimMatrix m1) (DimMatrix m2) = DimMatrix $ m1 *^ m2

cholM
  ::
  ( KnownNat y
  , KnownNat x
  , Field a
  , y ~ x
  )
  => DimMatrix D y x a
  -> DimMatrix D y x a
cholM (DimMatrix m) = DimMatrix $ delay $ chol $ trustSym $ computeS m

sumAllSM
  ::
  ( KnownNat y
  , KnownNat x
  , Num a)
  => DimMatrix D y x a
  -> a
sumAllSM (DimMatrix m) = sumAllS m

foldAllSM
  :: (KnownNat y, KnownNat x)
  => (Double -> Double -> Double)
  -> Double
  -> DimMatrix D y x Double
  -> Double
foldAllSM f initValue (DimMatrix m) = foldAllS f initValue m

detSM
  :: (KnownNat y, KnownNat x)
  => DimMatrix D y x Double
  -> Double
detSM (DimMatrix m) = detS m

trace2SM
  :: (KnownNat x, KnownNat y)
  => DimMatrix D x y Double
  -> Double
trace2SM (DimMatrix m) = trace2S $ computeS m

toDimMatrix
  ::
  ( Source r a
  , KnownNat m
  , KnownNat n
  )
  => DimVector r m a
  -> Int
  -> DimMatrix D m n a
toDimMatrix (DimVector arr) desiredSize =
  DimMatrix (toMatrix arr desiredSize)

toMatrixM'
  :: (Source r a, KnownNat n, KnownNat m)
  => DimVector r n a
  -> DimMatrix D n m a
toMatrixM' (DimVector arr) =
  DimMatrix $ fromFunction (Z :. dimension :. 1) generator
  where
    dimension = size . extent $ arr
    generator (Z :. rows :. _) = linearIndex arr rows

zipWithDim
  ::
  ( Source r1 a
  , Source r2 b
  , KnownNat m
  , KnownNat n
  )
  => (a -> b -> c)
  -> DimMatrix r1 m n a
  -> DimMatrix r2 m n b
  -> DimMatrix D m n c
zipWithDim f (DimMatrix mat1) (DimMatrix mat2) =
  DimMatrix $ zipWith f mat1 mat2

zipWithArray
  ::
  ( KnownNat n
  , KnownNat m
  )
  => (a -> b -> c)
  -> DimVector D n a
  -> DimMatrix D n m b
  -> DimMatrix D n m c
zipWithArray f (DimVector array1) (DimMatrix array2) =
  DimMatrix $ zipWith f (toMatrix array1 n) array2
  where
    (Z :. n :. _) = extent array2

identM
  :: forall m n a.
  ( KnownNat m
  , KnownNat n
  , GPConstraint a
  , m ~ n
  )
  => DimMatrix D m n a
identM =
  let dim = fromEnum $
            natVal @m @Proxy Proxy in DimMatrix $
            identD dim

delayMatrix
  ::
  ( KnownNat m
  , KnownNat n
  , Source r a
  )
  => DimMatrix r m n a
  -> DimMatrix D m n a
delayMatrix (DimMatrix matrix) = DimMatrix (delay matrix)

linearSolveM
  ::
  ( Field a
  , AllConstrained KnownNat '[m, n]
  )
  => DimMatrix D m m a
  -> DimMatrix D m n a
  -> Maybe (DimMatrix D m n a)
linearSolveM (DimMatrix mat1) (DimMatrix mat2) =
  case linearSolveS mat1 mat2 of
    Nothing  -> Nothing
    Just sol -> Just . delayMatrix . DimMatrix $ sol

vectorLength
  :: forall r m a. KnownNat m
  => DimVector r m a
  -> Int
vectorLength _ = fromEnum $ natVal (Proxy @m)

matrixRowsNum
  :: forall r m n a. KnownNat m
  => DimMatrix r m n a
  -> Int
matrixRowsNum _ = fromEnum $ natVal (Proxy @m)

matrixColsNum
  :: forall r m n a. KnownNat n
  => DimMatrix r m n a
  -> Int
matrixColsNum _ = fromEnum $ natVal (Proxy @n)

randomMatrixD
  :: forall a g m n.
  ( RandomGen g
  , Random a
  , Unbox a
  , Floating a
  , KnownNat m
  , KnownNat n
  )
  => g
  -> (Int, Int)
  -> DimMatrix D m n a
randomMatrixD gen (rows, cols) =
  let randomList = take (rows * cols) (normals gen) in
    DimMatrix . delay $ fromListUnboxed (Z :. rows :. cols) randomList

type GPConstraint a =
  ( Field a
  , Random a
  , Unbox a
  , Floating a
  , Eq a
  )
{-
hermToMatrixM
  :: Herm a
  -> DimMatrix D m n a
hermToMatrixM = undefined
-}
pinvSM
  ::
  ( KnownNat y
  , KnownNat x
  , Field a
  , Numeric a
  , y ~ x
  )
  => DimMatrix D y x a
  -> DimMatrix D y x a
pinvSM (DimMatrix m) = DimMatrix . delay $ pinvS m

mulPM
  :: forall y1 x1 y2 x2 a.
  ( AllConstrained KnownNat [x1, x2, y1, y2]
  , Numeric a
  , x1 ~ y2
  )
  => DimMatrix D y1 x1 a
  -> DimMatrix D y2 x2 a
  -> DimMatrix D y1 x2 a
mulPM (DimMatrix m) (DimMatrix n) =
  DimMatrix . delay . runIdentity $ m `mulP` n

eigSHM
  :: (Field a, Numeric a)
  => Herm a
  -> (DimVector D m Double, DimMatrix D m m a)
eigSHM hermM =
  (DimVector $ delay $ fst out, DimMatrix $ delay $ snd out)
  where
    out = eigSH hermM


text:
module Statistics.Information.Continuous.Entropy where

import Data.Matrix
import Numeric.SpecFunctions
import Statistics.Information.Utils.KNN
import Statistics.Information.Utils.Random
import System.Random

entropy :: RandomGen g => g -> Int -> Int -> Matrix Double -> Double
entropy g k base xs | k < nrows xs =
  let noised = toLists $ noisy g xs
      distances = knnDistances chebyshevDistance noised k in
    (const + (fromIntegral d) * mean (map log distances)) /
    (log $ fromIntegral base)
  where
    d = ncols xs
    n = nrows xs
    const = digamma (fromIntegral n) - digamma (fromIntegral k) +
            (fromIntegral d) * log 2
    mean xs = sum xs / (fromIntegral $ length xs)

centropy :: RandomGen g => g -> Int -> Int -> Matrix Double -> Matrix Double ->
            Double
centropy g k base xs ys = joint - hy where
  joint = entropy g1 k base (xs <|> ys)
  hy = entropy g2 k base ys
  (g1, g2) = split g

max_stars_repo_name: eligottlieb/Shannon

max_issues_repo_name: eligottlieb/Shannon

max_forks_repo_name: eligottlieb/Shannon


text:
{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE Strict #-}

module Main where

import qualified Data.ByteString as BS
import qualified Data.Vector as V
import           Data.Vector (Vector)

import           Control.Monad (replicateM)
import           Control.Arrow ((&&&))
-- import           System.Random

import           Linear.Trace

import           Data.Vector.Split (chunksOf)
import           Data.List (nub)

import           Linear.Vector

import           Net
import           Dense --(softplus, sigmoid, reLU)
import           Backprop (Net, computeNetState)
import           Parser
import           Utils

import           System.Random.MWC
import           Statistics.Distribution
import           Statistics.Distribution.Normal

import           Control.DeepSeq (force)
import           Control.Exception (evaluate)

imageSize :: Int
imageSize = 28*28

stepSize :: Double
stepSize = 0.1
-- stepSize = 0.15

genWeight :: GenIO -> Int -> IO Double
genWeight gen numWeights = genContVar (normalDistr 0 (1/sqrt (fromIntegral numWeights))) gen

genBias :: GenIO -> Int -> IO Double
genBias gen _numWeights = genContVar (normalDistr 0 1) gen

actFn :: Floating a => a -> a
actFn = sigmoid

boolToDouble :: Bool -> Double
boolToDouble False = 0
boolToDouble True  = 1

instance Trace Vector where
  diagonal m =
    V.fromList $ map (\i -> (m V.! i) V.! i) [0..size-1]
    where
      size = V.length m

classify :: Vector Double -> Vector Double
classify v =
  let idx = V.maxIndex v
  in fmap boolToDouble (oneHotEncode (fromIntegral idx))

main :: IO ()
main = do
  -- setStdGen (mkStdGen 200)
  -- setStdGen (mkStdGen 203)
  genIO <- createSystemRandom

  trainLabels0 <- parseLabels <$> BS.readFile "images/train/train-labels-idx1-ubyte"
  trainImages <- parseImages <$> BS.readFile "images/train/train-images-idx3-ubyte"

  testLabels0 <- parseLabels <$> BS.readFile "images/test/t10k-labels-idx1-ubyte"
  testImages <- parseImages <$> BS.readFile "images/test/t10k-images-idx3-ubyte"

  let trainLabels = fmap (fmap boolToDouble) trainLabels0
      testLabels  = fmap (fmap boolToDouble) testLabels0

  let testSampleSize   = 1000
      sampleTestImages = V.take testSampleSize testImages
      sampleTestLabels = V.take testSampleSize testLabels

  let trainLabelsAndImages
        = V.zip trainImages trainLabels

  initialNet <-
    initNet V.replicateM V.fromList imageSize [16, 16, 10] actFn (genWeight genIO) (genBias genIO)
      :: IO (Net Vector Double)

  minibatches <- evaluate $ force $ chunksOf 10 (V.take 700 trainLabelsAndImages)

  putStrLn "Minibatches created."

  let trainedNet =
        -- train 50 stepSize initialNet (map V.fromList $ chunksOf 10 (take 700 trainLabelsAndImages))
        train 50 stepSize initialNet minibatches

        -- train 500 stepSize actFn initialNet [V.fromList $ take 3 trainLabelsAndImages]
        -- train 10 stepSize actFn initialNet (take 2 (map V.fromList $ chunksOf 200 trainLabelsAndImages))

  putStrLn "Measuring training set accuracy..."
  putStr "Training accuracy: "
  putStr (show (netTestAccuracy classify trainedNet (V.take 1000 trainImages) (V.take 1000 trainLabels)*100))
  putStrLn "%"

  putStrLn "Measuring testing set accuracy..."
  putStr "Test accuracy: "
  putStr (show (netTestAccuracy classify trainedNet sampleTestImages sampleTestLabels*100))
  putStrLn "%"


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TupleSections       #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE TypeOperators       #-}

module Main where
import           Debug.Trace
import           Control.Monad
import qualified Data.ByteString              as B
import qualified Data.ByteString.Lazy         as L
import           Data.List
import           Data.Serialize               (decodeLazy, encodeLazy)
import qualified Data.Vector.Storable         as V
import           Grenade
import           Numeric.LinearAlgebra        (maxIndex)
import qualified Numeric.LinearAlgebra.Static as S
import           System.Random
import           Graphics.Image              as I hiding (map)
import           Graphics.Image.Interface    as I hiding (map)
import           Graphics.Image.Processing   as I hiding (map)
import qualified Numeric.LinearAlgebra as NLA
import qualified Numeric.LinearAlgebra.Static as NLA
import qualified Numeric.LinearAlgebra.Data as NLA

readF::FilePath->IO (S ('D3 100 100 3))
readF path =do
     Right x<- (fmap (resize Bilinear Edge (100,100)))<$>readImage path::IO (Either String (Image VS RGB Double))
     let Just res=fromStorable $V.concatMap (\(PixelRGB r g b) -> V.fromList [r,g,b]) $ toVector x
     return res

scons a !b = V.cons a b

elemAt = (V.!)

toImg::S ('D3 100 100 3)->Image VS RGB Double
toImg x = fromVector (100,100) $groupRGB $ext x
          where groupRGB v= if V.null v then V.empty else (PixelRGB (v `elemAt` 0) (v `elemAt` 1) (v `elemAt` 2)) `scons` groupRGB (V.drop 3 v)
                ext::S ('D3 100 100 3)->V.Vector Double
                ext (S3D m) = NLA.flatten$NLA.extract m

 type FE
  = Network
    '[ Convolution 3 24 5 5 1 1,  Relu, Pooling 2 2 2 2
     , Reshape
     , FullyConnected 55296 1024, Logit
     , FullyConnected 1024 55296, Relu
     , Reshape
     , Deconvolution 24 3 6 6 2 2, Logit
     ]
    '[ 'D3 100 100 3, 'D3 96 96 24, 'D3 96 96 24, 'D3 48 48 24
     , 'D1 55296
     , 'D1 1024, 'D1 1024
     , 'D1 55296, 'D1 55296
     , 'D3 48 48 24
     , 'D3 100 100 3,'D3 100 100 3]

trainOne::LearningParameters->FE->S ('D3 100 100 3)->FE
trainOne rate !net i=
  train rate net i i

trainT::FE->Int->S ('D3 100 100 3)->FE
trainT !net x img=
    let lp = LearningParameters 0.0001 0.9 0.0005
    in  foldl' (\n i ->
      trace (show i) $trainOne (lp{ learningRate = learningRate lp *(0.999 ^ i)}) n img)
      net [0..x]

main' = do
  i<-readF "image/n11669921_50878.JPEG"
  writeImage "newimg.jpg" $ toImg i

main = do
  i<-readF "image/n11669921_50878.JPEG"
  x<-Prelude.read<$>getLine
  net<-randomNetwork::IO FE
  let net'=trainT net x i
  writeImage "newimg.jpg" $ toImg $runNet net' i


text:
{-# LANGUAGE RecordWildCards #-}
module Statistics.BBVI.StepSize
  ( KucP(..)
  , rhoKuc
  , defaultKucP
  , Rho
  )
where

import qualified Data.Vector                   as V
import           Statistics.BBVI.Propagator     ( DistCell(..)
                                                , Gradient
                                                , Memory
                                                )

-- | step size is a vector the same length as n parameters in
-- the distribution
type Rho = V.Vector Double

-- | calculates a step size following Kucukelbir et al 2017
rhoKuc
  :: KucP -- ^ parameters
  -> Gradient -- ^ gradient at time t
  -> DistCell a -- ^ distribution cell (needed for memory data
                -- e.g. \( s_{t-1} \) in Kucukelbir et al 2017)
  -> (Memory, Rho) -- ^ change in memory, stepsize
rhoKuc KucP {..} gra (Node time memory _dist) =
  ( deltaM
  , V.zipWith
    (\ds s ->
      eta
        *  (fromIntegral time)
        ** (negate 0.5 + eps)
        *  (1.0 / (tau + sqrt (s + ds)))
    )
    deltaM
    memory
  )
 where
  deltaM = V.zipWith (\g s -> alpha * g ^ (2 :: Int) - alpha * s) gra memory
rhoKuc _kp _gra (U{}) = error "called step size cell on an update cell!"

-- | parameters for 'rhoKuc'
data KucP = KucP
  { alpha :: Double
  , eta :: Double
  , tau :: Double
  , eps :: Double
  } deriving (Show, Eq, Ord, Read)

-- | default parameters for 'rhoKuc.' eta is what you probably want to
-- tune: kucukelbir et al try 0.01 0.1 1 10 100
defaultKucP :: KucP
defaultKucP = KucP { alpha = 0.1, eta = 0.1, tau = 1.0, eps = 1e-16 }

max_stars_repo_name: massma/propagator-bbvi

max_issues_repo_name: massma/propagator-bbvi

max_forks_repo_name: massma/propagator-bbvi


text:
module LU ( luFactor, luFactorM ) where

import Control.Monad( zipWithM_ )
import Control.Monad.ST( ST, runST )

import Numeric.LinearAlgebra
import qualified Numeric.LinearAlgebra.Matrix as M
import qualified Numeric.LinearAlgebra.Vector as V

luFactor :: (BLAS3 e) => Matrix e -> Either Int (Matrix e, [Int])
luFactor a = runST $ do
    ma <- M.newCopy a
    luFactorM ma >>=
        either (return . Left) (\pivots -> do
            a' <- M.unsafeFreeze ma
            return $ Right (a',pivots)
        )

{-
 - Recursive LU factorization with row pivoting.  Takes a matrix
 - A and factors it as P A = L U, where P is a permutation matrix,
 - L is a lower triangular matrix with ones along the diagonal, and
 - U is an upper triangular matrix.  On successful return, the values of
 - L and U are stored in A, and a list of the row swaps are returned.
 - On failure, the index of the failing column is returned.
 -}
luFactorM :: (BLAS3 e) => STMatrix s e -> ST s (Either Int [Int])
luFactorM a = do
    (m,n) <- M.getDim a
    let mn = min m n
        nleft = mn `div` 2

    case undefined of
        _ | mn > 1 ->
            M.withSplitColsAtM nleft a $ \a_1 a_2 ->
            M.withSplitRowsAtM nleft a_1 $ \a11 a21 ->
            M.withSplitRowsAtM nleft a_2 $ \a12 a22 ->
            luFactorM a_1 >>=
                either (return . Left) (\pivots -> do
                    zipWithM_ (M.swapRows a_2) [ 0.. ] pivots
                    M.triSolvMatrixM_ LeftSide NoTrans (Tri Lower Unit a11) a12
                    M.addMulMatrixWithScalesM_ (-1) NoTrans a21 NoTrans a12 1 a22
                    luFactorM a22 >>=
                        either (return . Left . (nleft+)) (\pivots' -> do
                            zipWithM_ (M.swapRows a21) [ 0.. ] pivots'
                            return $ Right (pivots ++ map (nleft+) pivots')
                        )
                )

        _ | mn == 1 ->
            M.withColM a 0 $ \x ->
                V.getWhichMaxAbs x >>= \(i,e) ->
                    if (e /= 0)
                        then do
                            V.scaleM_ (1/e) x
                            V.read x 0 >>= V.write x i
                            V.write x 0 e
                            return $ Right [i]
                        else
                            return $ Left 0

        _ | otherwise ->
            return $ Right []

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.Poisson
-- Copyright : (c) 2009, 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The Poisson distribution.  This is the discrete probability
-- distribution of a number of events occurring in a fixed interval if
-- these events occur with a known average rate, and occur
-- independently from each other within that interval.

module Statistics.Distribution.Poisson
    (
      PoissonDistribution
    -- * Constructors
    , poisson
    , poissonE
    -- * Accessors
    , poissonLambda
    -- * References
    -- $references
    ) where

import Control.Applicative
import Data.Data            (Data, Typeable)
import GHC.Generics         (Generic)
import Numeric.SpecFunctions (incompleteGamma,logFactorial)
import Numeric.MathFunctions.Constants (m_neg_inf)

import qualified Statistics.Distribution as D
import qualified Statistics.Distribution.Poisson.Internal as I
import Statistics.Internal

newtype PoissonDistribution = PD {
      poissonLambda :: Double
    } deriving (Eq, Typeable, Data, Generic)

instance Show PoissonDistribution where
  showsPrec i (PD l) = defaultShow1 "poisson" l i
instance Read PoissonDistribution where
  readPrec = defaultReadPrecM1 "poisson" poissonE

instance D.Distribution PoissonDistribution where
    cumulative (PD lambda) x
      | x < 0        = 0
      | isInfinite x = 1
      | isNaN      x = error "Statistics.Distribution.Poisson.cumulative: NaN input"
      | otherwise    = 1 - incompleteGamma (fromIntegral (floor x + 1 :: Int)) lambda

instance D.DiscreteDistr PoissonDistribution where
    probability (PD lambda) x = I.probability lambda (fromIntegral x)
    logProbability (PD lambda) i
      | i < 0     = m_neg_inf
      | otherwise = log lambda * fromIntegral i - logFactorial i - lambda

instance D.Variance PoissonDistribution where
    variance = poissonLambda

instance D.Mean PoissonDistribution where
    mean = poissonLambda

instance D.MaybeMean PoissonDistribution where
    maybeMean = Just . D.mean

instance D.MaybeVariance PoissonDistribution where
    maybeStdDev   = Just . D.stdDev

instance D.Entropy PoissonDistribution where
  entropy (PD lambda) = I.poissonEntropy lambda

instance D.MaybeEntropy PoissonDistribution where
  maybeEntropy = Just . D.entropy

-- | Create Poisson distribution.
poisson :: Double -> PoissonDistribution
poisson l = maybe (error $ errMsg l) id $ poissonE l

-- | Create Poisson distribution.
poissonE :: Double -> Maybe PoissonDistribution
poissonE l
  | l >=  0   = Just (PD l)
  | otherwise = Nothing

errMsg :: Double -> String
errMsg l = "Statistics.Distribution.Poisson.poisson: lambda must be non-negative. Got "
        ++ show l

-- $references
--
-- * Loader, C. (2000) Fast and Accurate Computation of Binomial
--   Probabilities. <http://projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf>
-- * Adell, J., Lekuona, A., and Yu, Y. (2010) Sharp Bounds on the
--   Entropy of the Poisson Law and Related Quantities
--   <http://arxiv.org/pdf/1001.2897.pdf>

max_stars_repo_name: vmchale/statistics

max_issues_repo_name: vmchale/statistics

max_forks_repo_name: vmchale/statistics


text:
{-# Language RecordWildCards #-}
-- |
-- Module    : GRN.EM
-- Copyright : (c) 2011 Jason Knight
-- License   : BSD3
--
-- Maintainer  : jason@jasonknight.us
-- Stability   : experimental
-- Portability : portable
--
-- Expectation Maximization Algorithm at work
--

module GRN.EM where

import GRN.StateTransition
import GRN.Render
import GRN.Sparse
import GRN.Types
import GRN.Utils

import qualified Data.Map as M
import Data.List
import Data.Char
import Data.Maybe
import Control.Monad
import Text.Printf
import System.Console.ParseArgs

import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Generic as G

import Control.Parallel.Strategies
import Data.Vector.Strategies

import Data.Graph.Inductive
import Data.Graph.Analysis.Algorithms.Common (componentsOf)

import Statistics.Distribution
import Statistics.Distribution.Normal

{-import Numeric.LBFGSB -}

import Graphics.Gnuplot.Simple

smoothingVar = 0.2 -- Variance of smoothing gaussian

ssaMap :: [Gene] -> SSA -> SSAMap
ssaMap glist s = M.fromList $ zip glist (U.toList s)

-- | Returns a map Gene->Double of any gene that has a measurement recorded.
measureMap :: ParseData -> MeasureMap
measureMap p = M.map fromJust $ M.fromList $ filter (isJust.snd) $ map (\(x,gi)-> (x,measurement gi)) (M.toList p)

-- | Returns a list of Uncertain locations.
uncertainVals :: KmapSet -> [UncertainSpot]
uncertainVals ks = filter filtFun $ concatMap explode (M.elems ks)
  where explode (Kmap g upstream x) = map (\(pos,entry)->(g,pos,entry)) (M.toList x)

filtFun :: (Gene,[Int],Kentry) -> Bool
filtFun (g,pos,x) = case x of
  C _ -> True
  X -> True
  _ -> False

fillKentries :: [UncertainSpot] -> KmapSet -> Theta -> KmapSet
fillKentries spots ks mm = foldr fillKentry ks (zip spots (U.toList mm))

fillKentry :: (UncertainSpot,Double) -> KmapSet -> KmapSet
fillKentry ((g,pos,_),val) ks = M.adjust fill g ks
  where fill (Kmap g up x) = Kmap g up (M.insert pos (V 0 val) x)

-- truncated Norm:: Mu  -> Var    -> Query
truncatedNorm :: Double -> Double -> Double -> Double
truncatedNorm mu var pt = 1/var * std ((pt-mu)/var) / denom
  where std = density standard
        cstd = cumulative standard
        denom = cstd ((1-mu)/var) - cstd ((0-mu)/var)

-- Given a theta vector, ssa vector, and a set of measurements, multiply the
-- measurement lookups convolved with a measurement gaussian against each other
-- to get a joint lookup.
msLookup :: SSAMap -> MeasureMap -> Double
msLookup ssam mm = product $ M.elems $ M.mapWithKey draw mm
  where draw gene meas = truncatedNorm (ssam M.! gene) smoothingVar meas

runAndSplit :: Args String -> KmapSet -> Masks -> Attractors
runAndSplit args ks masks = map ((\(d,s)->(d, ssaMap (M.keys ks) (ssdToSSA s))).normalizeSSD) raw
  where base = simulateCSC args csc (U.replicate n (1/(fromIntegral n)))
        dok = kmapToDOK ks 0
        csc = dokToCSC dok
        n = 2^(M.size ks)
        raw = map (filtSSD base) masks
        filtSSD base mask = zero `U.update` inds
          where zero = U.replicate n 0.0
                inds = U.map (\i->(i,base U.! i)) mask

runAndSplit2 :: Args String -> KmapSet -> Masks -> Attractors
runAndSplit2 args ks masks = map ((\(d,g)->(d, ssaMap (M.keys ks) (genSSA g))).normalizeGraph) raw
  where ssd = simulateDOKUnif args.kmapToDOK ks $ 0
        raw = componentsOf $ convertProbsVG ssd stripG
        stripG = last . take 10.iterate stripTransNodes $ kmapToStateGraph ks

-- Given an unfilled KMS and a new theta set, fill the kmap, then generate an
-- Attractors list
generateNewAttractors :: Args String -> KmapSet -> Masks -> [UncertainSpot] -> Theta -> Attractors
generateNewAttractors args ks masks us theta = runAndSplit args newks masks
  where newks = fillKentries us ks theta

expectationRun :: EMData -> Theta -> Double
expectationRun EMData{..} new = emNormSum * (G.sum (V.generate (length emAttractors) eachAttractor))
  where newAgivenTheta = generateNewAttractors emArgs emUnfilledKMS emCompNodes emUncertainLocs new
        eachAttractor i = msLookup oldssam emMeasurements * w * log (newAgivenO * msLookup newssam emMeasurements)
          where (w,oldssam) = emAttractors !! i
                (newAgivenO, newssam) = newAgivenTheta !! i

normSum :: Attractors -> MeasureMap -> Double
normSum atracs mm = (-1) / (sum $ map (\(w,issam)-> w * msLookup issam mm) atracs)

posteriorProb :: EMData -> Theta -> U.Vector Double
posteriorProb EMData{..} new = G.convert $ G.map (*(1/normSum)) (V.generate (length emAttractors) eachAttractor)
  where normSum = sum $ map (\(w,issam)-> w * msLookup issam emMeasurements ) newAgivenTheta
        newAgivenTheta = generateNewAttractors emArgs emUnfilledKMS emCompNodes emUncertainLocs new
        eachAttractor i = msLookup newssam emMeasurements * newAgivenO
          where (newAgivenO, newssam) = newAgivenTheta !! i

{-fullEM :: EMData -> Int -> Theta -> [(Theta,EMData)]-}
{-fullEM emd 0 oldtheta = [(oldtheta,emd)]-}
{-fullEM emd maxiter oldtheta -}
  {-| thetaClose = (oldtheta,emd) : [(newtheta,newemd)]-}
  {-| otherwise = (oldtheta,emd) : fullEM newemd (maxiter-1) newtheta-}
  {-where newemd = emd { emAttractors=newAtracs, emNormSum=newNormSum } -}
        {-newAtracs = generateNewAttractors (emArgs emd) (emUnfilledKMS emd) -}
                        {-(emCompNodes emd) (emUncertainLocs emd) newtheta-}
        {-newNormSum = normSum newAtracs (emMeasurements emd)-}
        {-newtheta = optimize (expectationRun emd) oldtheta-}
        {-thetaClose = (<1e-5).U.sum.U.map (**2).U.zipWith (-) newtheta $ oldtheta-}

--data EMData = EMData {
--      emAttractors :: Attractors
--    , emUnfilledKMS :: KmapSet
--    , emComps :: [ColoredStateGraph]
--    , emCompNodes :: Masks
--    , emNormSum :: Double
--    , emArgs :: Args String
--    , emMeasurements :: MeasureMap
--    , emUncertainLocs :: [UncertainSpot]
--    }

{-emRun :: Args String -> ParseData -> IO () -}
{-emRun args p = do-}
  {-let-}
    {-gen = gotArg args "generate"-}
    {-ks = buildKmaps p -}
    {-ms = measureMap p-}
    {-uvals = uncertainVals ks-}

    {-ng = M.size ks-}
    {-n = length uvals-}
    {-printList :: [U.Vector Double] -> String-}
    {-printList x = concat $ intersperse "\n" (map (concatMap (printf "%7.3f").G.toList) $ x)-}

    {---startTheta = U.fromList [0.379,0.5,1.0]-}
    {-startTheta = U.replicate n 0.5-}

    {-newks = fillKentries uvals ks startTheta-}
    {-stripG = componentsOf . last . take 10.iterate stripTransNodes $ kmapToStateGraph newks-}
    {-stripNodes = map (U.fromList.sort.nodes) stripG-}
    {-startAtracs = runAndSplit args newks stripNodes-}
    {-norms = normSum startAtracs ms-}

    {-emdat = EMData startAtracs ks stripG stripNodes norms args ms uvals-}

    {-mapt = map (U.fromList.map fromIntegral.dec2bin n) [0..2^n-1]-}
    {---mapt = map (U.replicate n) [0,0.1,0.9,1.0]-}
    {---mapt = map (U.fromList) [replicate n 0.0, replicate n 0.001]-}
    {---prVal = map (expectationRun emdat) mapt-}
    {-prVal = map (posteriorProb emdat) mapt -}
    {-prVal2 = map (expectationRun emdat) mapt -}

    {-elOpt = optimize (expectationRun emdat) (U.replicate n 0.5) -}
    {-path = fullEM emdat 10 startTheta-}

    {-xs = clamped [0,0.05..1]-}

    {-clamped xs = map (\x-> if (x>1.0) then 1.0 else (if (x<0.0) then 0.0 else x)) xs-}

    {-printSeg :: String -> String -> IO ()-}
    {-printSeg title st = do-}
          {-let l = 2 + length title-}
              {-upperTitle = map toUpper title-}
              {-width = 80-}
              {-half = (width-l) `div` 2-}
              {-stars = replicate half '*'-}
          {-printf "\n%s %s %s\n\n%s\n\n%s\n\n" stars upperTitle stars st (replicate width '*')-}

  --print n

  --forM_ mapt $ \t -> do
  --  let path = fullEM emdat 10 t
  {-printSeg "theta path" (printList.map (\(t,e)->t) $ path)-}
  {-printSeg "posterior probs" (printList.map (\(t,e)->posteriorProb e t) $ path)-}

  {---printSeg "attractors" (concatMap (printf "%7s".map intToDigit.dec2bin ng.G.head) stripNodes)-}
  {-printSeg "attractors" (concatMap (printf "%7d".G.head) stripNodes)-}
  {-print $ M.keys ks-}

  {---printf "%d Evaluations\n" ((length xs)^2)-}
  {-plotFunc3d [XLabel "X", YLabel "Y", XRange (0,1), YRange (0,1)] [] xs xs -}
     {-(\x y -> (-1)*(expectationRun emdat (U.fromList [0.5,x,y])))-}

  {---putStrLn "Expectations:"-}
  {---mapM_ (printf "%7.3f") prVal2-}
  {---putStrLn ""-}

  {---print elOpt-}

  {---putStrLn "Posterior Probs:"-}
  {---mapM_ (sequence.printList.U.fromList.sort.U.toList) prVal-}
  {---putStrLn ""-}
  {---mapM_ (printf "%5.2f") (map U.sum prVal)-}
  {---putStrLn ""-}

  {---putStrLn ""-}
  {---putStrLn ""-}
  {---print $ generateNewAttractors args ks stripNodes uvals (last mapt)-}
  {----}
  {---putStrLn ""-}
  {---putStrLn "Expectations:"-}
  {---mapM_ (printf "%7.3f\n") prVal2-}
  {---putStrLn ""-}

  {---putStrLn ""-}
  {---putStrLn ""-}
  {---putStrLn "**** Optimal Posterior Probs ****"-}
  {---mapM_ (sequence.printList) [posteriorProb emdat elOpt]-}
  {---mapM_ (printf "%7.3d".(U.head)) stripNodes-}
  {---putStrLn ""-}
  {---putStrLn "*********************************"-}
  {----}
  {---putStrLn ""-}
  {---printf "Optimal Expectation: %5.2f\n" (expectationRun emdat elOpt)-}
  {---putStrLn ""-}
  {---putStrLn "**** Optimal Theta ****"-}
  {---mapM_ (sequence.printList) [elOpt]-}

  {-printSeg "uncertain locs" (show uvals)-}
  {-printSeg "measurements" (show ms)-}

  {-when gen $ do-}
    {-let newks2 = fillKentries uvals ks (fst.last $ path)-}
        {-finalGraph = simulate args $ kmapToStateGraph newks2-}
        {-stripG2 = componentsOf . last . take 10.iterate stripTransNodes $ finalGraph-}
        {-stripNodes2 = map (U.fromList.sort.nodes) stripG2-}

    {-printSeg "final attractors" (concatMap (printf "%7d".G.head) stripNodes2)-}
    {-drawStateGraph finalGraph args-}

    {---let newks2 = fillKentries uvals ks (fst.head $ path)-}
    {---    interGraph = kmapToStateGraph newks2-}
    {---    ssd = simulateDOK args $ kmapToDOK newks2-}

    {----- WATCH OUT TODO I just found a huge bug: StateTransition.simulate does-}
    {----- not respect edge weights when moving probabilities around-}

max_stars_repo_name: binarybana/grn-pathways

max_issues_repo_name: binarybana/grn-pathways

max_forks_repo_name: binarybana/grn-pathways


text:
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE RankNTypes #-}
module Numeric.Tensor where

import qualified Prelude
import Control.Category
import Data.Kind(Type)
import Dependent.Size
import Data.Array.Accelerate((+),(*))
import qualified Data.Array.Accelerate as Acc
import Numeric.Vector.Sized

class Semigroup g where
  (<>) :: g -> g -> g

-- Associativity:
-- (a <> b) <> c = a <> (b <> c)

class Semigroup m => Monoid m where
  zero :: m

-- Identity:
-- a <> zero = a = zero <> a

class Monoid g => Group g where
    negate :: g -> g

-- Additive inverse
-- a <> negate a = zero = negate a <> a

-- Multiplication
class Group r => Ring r where
  one :: r
  (.*.) :: r -> r -> r
-- one .*. x = x = x .*. one

class Ring k => Field k where
  recip :: k -> k

-- x .*. recip x = one = recip x .*. x

class (Field k, Group v) => Module k v where
  (*.) :: k -> v -> v
  (.*) :: v -> k -> v

newtype Vect k v w = Vect {linear :: v -> w}

--
-- instance ShapeSize size => Semigroup (Sized size) where
--   -- (<>) :: Sized size -> Sized size -> Sized size
--   x <> y = zipWith (+) x y
--
-- instance ShapeSize size => Monoid (Sized size) where
--   -- zero :: Sized size
--   zero = generate (\_shape -> 0)
--
-- instance ShapeSize size => Group (Sized size) where
--   negate v = map Acc.negate v
--
instance Semigroup (Acc.Exp Acc.Double) where
  a <> b = a + b

instance Monoid (Acc.Exp Acc.Double) where
  zero = 0

instance Group (Acc.Exp Acc.Double) where
  negate = Acc.negate

instance Ring (Acc.Exp Acc.Double) where
  one = 1
  a .*. b = a * b

instance Field (Acc.Exp Acc.Double) where
  recip x = Acc.recip x

-- instance (ShapeSize size) => Module (Acc.Exp Acc.Double) (Sized size) where
--    r *. v = map (* r) v
--    v .* r = map (* r) v

instance Field k => Category (Vect k) where
  id = Vect id
  Vect f . Vect g = Vect (f . g)

class (Field f, Group v)  => VectorSpace f v where
  distribute :: f -> v -> v

-- instance ShapeSize size => VectorSpace (Acc.Exp Double) (Sized size) where
--   type (Sized size1 ⊗ Sized size2) = Sized (TensorSize size1 size2)

max_stars_repo_name: mixed-signals/mixed-signals

max_issues_repo_name: mixed-signals/mixed-signals

max_forks_repo_name: mixed-signals/mixed-signals


text:
-- stack bench --ba "--output bench.html"

import Criterion.Main
import Numeric.LinearAlgebra
import Linear.Matrix
import Linear.V4

-- the function to benchmark
fib m | m < 0     = error "negative!"
      | otherwise = go m
  where
    go 0 = 0
    go 1 = 1
    go n = go (n-1) + go (n-2)

hmatrixMult :: Double -> Matrix R
hmatrixMult m =
  let m1 = matrix 4 [m..(m + 15)] * ident 4
      m2 = (4><4) [m..] + row [1, 2, 3, 4]
  in m1 Numeric.LinearAlgebra.<> m2

linearMult :: Double -> V4 (V4 Double)
linearMult m =
  let m1 = V4 (V4 (m + 1) m m m) (V4 m (m + 6) m m) (V4 m m (m + 11) m) (V4 m m m (m + 16))
      m2 = (V4 (V4 (m + 1) (m + 3) (m + 5) (m + 7))
               (V4 (m + 5) (m + 7) (m + 9) (m + 11))
               (V4 (m + 9) (m + 11) (m + 13) (m + 15))
               (V4 (m + 13) (m + 15) (m + 17) (m + 19)))
  in m1 !*! m2

--data Row4 = Row4 [Double]
--data Matrix44 = Matrix44 [Row4]

matMult :: [[Double]] -> [[Double]] -> [[Double]]
matMult a b =
  let r0 = getRow a 0
      r1 = getRow a 1
      r2 = getRow a 2
      r3 = getRow a 3
      c0 = getColumn b 0
      c1 = getColumn b 1
      c2 = getColumn b 2
      c3 = getColumn b 3
  in [ [dotProduct r0 c0, dotProduct r0 c1, dotProduct r0 c2, dotProduct r0 c3]
     , [dotProduct r1 c0, dotProduct r1 c1, dotProduct r1 c2, dotProduct r1 c3]
     , [dotProduct r2 c0, dotProduct r2 c1, dotProduct r2 c2, dotProduct r2 c3]
     , [dotProduct r3 c0, dotProduct r3 c1, dotProduct r3 c2, dotProduct r3 c3]
     ]

getRow :: [[Double]] -> Int -> [Double]
getRow m i = m !! i

getColumn :: [[Double]] -> Int -> [Double]
getColumn m i = [ getRow m 0 !! i
                , getRow m 1 !! i
                , getRow m 2 !! i
                , getRow m 3 !! i
                ]

dotProduct :: [Double] -> [Double] -> Double
dotProduct row col = sum $ zipWith (*) row col

-- Matrix44 [(Row4 [1, 2, 3, 4]), (Row4 [1,2,3,4]), (Row4 [1,2,3,5]), (Row4 [1,2,3,6])]

naiveMult :: Double -> [[Double]]
naiveMult m =
  let m1 = [ [m + 1, m, m, m]
           , [m, m + 6, m, m]
           , [m, m, m + 11, m]
           , [m, m, m, m + 16]
           ]
      m2 = [ [m + 1, m + 3, m + 5, m + 7]
           , [m + 5, m + 7, m + 9, m + 11]
           , [m + 9, m + 11, m + 13, m + 15]
           , [m + 13, m + 15, m + 17, m + 19]
           ]
  in m1 `matMult` m2

-- Our benchmark harness.
main :: IO()
main = defaultMain [ bgroup "hmatrix" [ bench "mult 1" $ nf hmatrixMult 1.0
                                      ]
                   , bgroup "linear" [ bench "mult 1" $ nf linearMult 1.0
                                     ]
                   , bgroup "naive" [ bench "mult 1" $ nf naiveMult 1.0
                                     ]
                   ]

max_stars_repo_name: DavidAntliff/RTC-Haskell

max_issues_repo_name: DavidAntliff/RTC-Haskell

max_forks_repo_name: DavidAntliff/RTC-Haskell


text:
{-|
Module      : ControlSystems.DynamicSystems.Conversion
Description : Converts continuous time systems to discrete time
Copyright   : (c) Ryan Orendorff, 2020
License     : BSD3
Stability   : experimental
-}
module ControlSystems.DynamicSystems.Conversions (
  c2d
)

where

import Numeric.LinearAlgebra

-- | Convert a continuous time linear differential system to a discrete form,
-- using state space models.
c2d :: Double                         -- ^ The sampling time
    -> Matrix Double                  -- ^ The state transition matrix, aka A
    -> Matrix Double                  -- ^ The state input matrix, aka B
    -> (Matrix Double, Matrix Double) -- ^ The discretized A & B system
c2d t a b = (a_d, b_d)
  where
    -- First we have to extract the sizes of all the matrices at run time.
    (m_a, n_a) = (rows a, cols a)
    (_  , n_b) = (rows b, cols b)

    -- Then manually make the correct block. This must be a square matrix
    -- because the matrix exponential function (`expm`) expects to only have
    -- square inputs (given that it is a series expansion of matrix to higher
    -- powers). We need to form this matrix
    -- ⌈ A B ⌉
    -- ⌊ 0 0 ⌋
    block = a ||| b
              ===
              konst 0 (m_a, n_a + n_b) -- I left in a bug here! ;-D

    -- Calculate the discrete time block matrix, which is the following
    -- ⌈ A_d B_d ⌉
    -- ⌊  0   I  ⌋
    exp_block = expm (scale t block)

    -- Finally we can extract the correct submatrices assuming we pull out the
    -- right pieces.
    a_d = subMatrix (0, 0) (m_a, n_a) exp_block
    b_d = subMatrix (0, n_a) (m_a, n_b) exp_block

max_stars_repo_name: ryanorendorff/convex

max_issues_repo_name: ryanorendorff/convex

max_forks_repo_name: ryanorendorff/convex


text:
{-# LANGUAGE ScopedTypeVariables, GADTs #-}
module Main where

import Control.Monad
import Data.Complex
import Data.Int
import Data.List
import Data.Mat
import qualified Data.Map
import Data.Maybe
import System.Environment
import System.Exit
import System.IO

-- A universal datatype for MAT data.
data MatVal =
    VDouble  Double
  | VInt64   Int64
  | VBool    Bool
  | VChar    Char
  | VComplex (Complex MatVal)
  | VStruct  (Data.Map.Map String MatVal)
  | VArray   Int [MatVal]
  deriving (Eq, Show)

-- | Transform from the types into the universal representation.
toMatVal :: MatData -> MatVal
toMatVal (MatData MatDouble d)               = VDouble d
toMatVal (MatData MatInt64  i)               = VInt64  i
toMatVal (MatData MatBool   b)               = VBool   b
toMatVal (MatData MatChar   c)               = VChar   c
toMatVal (MatData MatCell   d)               = toMatVal d
toMatVal (MatData (MatStruct (Known fs)) ds) = VStruct (Data.Map.fromList (zip fs (map toMatVal ds)))
toMatVal (MatData (MatArray t (Known d)) xs) = VArray d (map (toMatVal . MatData t) xs)
toMatVal (MatData (MatComplex t) (x :+ y))   = matZipWith (\ r i -> VComplex (r :+ i)) (toMatVal (MatData t x)) (toMatVal (MatData t y))

-- | Push down complex values.
matZipWith :: (MatVal -> MatVal -> MatVal) -> MatVal -> MatVal -> MatVal
matZipWith op (VStruct xs) (VStruct ys)   = VStruct (Data.Map.unionWith (matZipWith op) xs ys)
matZipWith op (VArray d xs) (VArray _ ys) = VArray d (zipWith (matZipWith op) xs ys)
matZipWith op x y = op x y

main :: IO ()
main = do
  (fp:_) <- getArgs
  -- Read an entire MAT file
  xs :: [(String, MatData)] <- liftM Data.Map.toList (readMatFile fp)
  putStrLn "Variables contained in MAT file:"
  -- Print all the contents
  mapM_ (\ (n, d) -> putStr (n ++ " ->\n  " ++ show (toMatVal d) ++ "\n")) xs
  exitSuccess

max_stars_repo_name: glutamate/matio

max_issues_repo_name: glutamate/matio

max_forks_repo_name: glutamate/matio


text:
-- | Post-mortem tool main.

import           Universum

import           Data.List (last)
import qualified Data.Map.Strict as M
import           System.FilePath
import           System.IO (hPutStrLn)
import           Text.Printf (hPrintf)

import           Pos.Util.Util (histogram)

import           Options
import           Statistics
import           Types

main :: IO ()
main = parseOptions >>= \case
    Overview sampleProb logDirs -> do
        showLogDirs logDirs
        err $ "sample probability: " ++ show sampleProb
        err ""
        xs <- forM logDirs $ flip processLogDirOverview sampleProb
        chart xs "times.png"
        err "wrote times chart"
    Focus txHash logDir         -> do
        err $ "transaction hash: " ++ show txHash
        err $ "logs directory: " ++ show logDir
        let focusFile = ("focus_" ++ extractName logDir ++ "_" ++ toString txHash) <.> "csv"
        runJSONFold logDir (focusF txHash) >>= focusToCSV focusFile
        err $ "wrote result to " ++ show focusFile
    TxRelay logDirs             -> do
        showLogDirs logDirs
        err ""
        for_ logDirs processLogDirTxRelay
    Throughput txWindow waitWindow logDirs   -> do
        showLogDirs logDirs
        err $ "tx window: " ++ show txWindow
        err $ "wait window: " ++ show waitWindow
        err ""
        for_ logDirs $ processLogDirThroughput txWindow waitWindow

showLogDirs :: [FilePath] -> IO ()
showLogDirs logDirs = do
    err "log directories: "
    for_ logDirs $ \d -> err $ " - " ++ show d

processLogDirOverview :: FilePath -> Double -> IO (String, Map TxHash (Maybe Timestamp))
processLogDirOverview logDir sampleProb = do
    err $ "processing log directory " ++ show logDir ++ " ..."

    (rc, g, mp, cr, ft) <-
        runJSONFold logDir $ (,,,,) <$> receivedCreatedF
                                    <*> graphF
                                    <*> memPoolF
                                    <*> txCntInChainF
                                    <*> txFateF
    let total    = M.size rc
        included = sort $ mapMaybe snd $ M.toList rc
        lost     = total - length included
    err $ "total number of received transactions: " ++ show total
    err $ "included in blockchain: " ++ show (length included)
    err $ "lost transactions: " ++ show lost
    err $ "blocks in chain: " ++ show (length cr)
    err $ "mem pool entries: " ++ show (length mp)

    let dirName   = extractName logDir
    let graphFile = getName "graph" dirName "png"
    b <- writeGraph graphFile g
    when b $ err $ "wrote graph png to " ++ show graphFile

    let csvFile = getName "csv" dirName "csv"
    txCntInChainMemPoolToCSV csvFile sampleProb cr mp
    err $ "wrote csv file to " ++ show csvFile
    let reportFile = getName "report" dirName "txt"
    void (reportTxFate reportFile ft)
    err $ "wrote report file to " ++ show reportFile
    err $ "processing log directory " ++ show logDir ++ " done"
    err ""
    return (dirName, rc)

processLogDirTxRelay :: FilePath -> IO ()
processLogDirTxRelay logDir = do
    err $ "processing log directory " ++ show logDir ++ " ..."
    m <- runJSONFold logDir txReceivedF
    err $ "total number of received transactions: " ++ show (M.size m)
    let xs         = sortBy (compare `on` \(_, ys) -> negate $ length ys) $ M.toList m
        dirName    = extractName logDir
        relaysFile = getName "relays" dirName "txt"
    withFile relaysFile WriteMode $ \h ->
        for_ xs $ \(tx, ys) ->
            hPutStrLn h $ toString tx ++ ": " ++ show ys
    err $ "wrote relay data to " ++ show relaysFile
    let hist     = M.toList $ histogram $ map (length . snd) xs
        histFile = getName "histogram" dirName "txt"
    withFile histFile WriteMode $ \h ->
        for_ hist $ uncurry $ hPrintf h "%4d:%6d\n"
    err $ "wrote histogram to " ++ show histFile

processLogDirThroughput :: Double -> Double -> FilePath -> IO ()
processLogDirThroughput txWindow waitWindow logDir = do
    err $ "processing log directory " ++ show logDir ++ " ..."
    (xs, ys) <- runJSONFold logDir $ (,) <$> txCntInChainF <*> memPoolF
    err $ "chain length: " ++ show (length xs) ++ " block(s)"
    err $ show (length ys) ++ " mem pool event(s)"
    let svgFile = getName "throughput" (extractName logDir) "svg"
    throughput svgFile txWindow waitWindow 1000 xs ys
    err $ "wrote chart to " ++ show svgFile

err :: String -> IO ()
err = hPutStrLn stderr

getName :: FilePath -> String -> String -> FilePath
getName template name ext =
    let dir   = takeDirectory template
        base  = takeBaseName template
        base' =  base ++ "_" ++ name
    in  dir </> base' <.> ext

extractName :: FilePath -> String
extractName = last . splitDirectories


text:
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE BangPatterns #-}

module Numeric.MCMC.RiemannianLangevin (
            MarkovChain(..)
          , Options, createOptions
          , runChain, localMean, perturb
          ) where

import Numeric.MCMC.Langevin hiding (Options, runChain)
import Numeric.LinearAlgebra
import Control.Arrow
import Control.Monad
import Control.Monad.Primitive
import Control.Monad.Reader
import System.Random.MWC
import System.Random.MWC.Distributions

-- | Options for the chain.  The target (expected to be a log density), the
--   right matrix product of its curvature and gradient, the inverse Fisher metric tensor,
--   its (Cholesky) square root, and the step size.
data Options = Options {
    _target              :: [Double] -> Double        -- Target (log density)
  , _curvatureXgradient  :: [Double] -> [Double]      -- Curvature right-multiplied by gradient
  , _invFisherMetric     :: [Double] -> Matrix Double -- Inverse Fisher metric tensor
  , _sqrtInvFisherMetric :: [Double] -> Matrix Double -- Square root of the tensor
  , _eps                 :: {-# UNPACK #-} !Double    -- Step size
  }

-- | A result with this type has a view of the chain options.
type ViewsOptions = ReaderT Options

-- | Construct Options (data constructor not exported).
-- FIXME sqrtm is not converging; alternatives?
createOptions :: ([Double] -> Double)     -- Target (log density)
              -> ([Double] -> [Double])   -- Gradient
              -> ([Double] -> [[Double]]) -- Hessian
              -> Double                   -- Step size
              -> Options                  -- Options
createOptions t g h =
    Options t curvatureXgradient invFisherMetric sqrtInvFisherMetric
  where curvatureXgradient xs =
            let mat = invFisherMetric xs <> fromColumns [fromList (g xs)]
            in  concat . toLists . trans $ mat
        invFisherMetric       = pinv . fromRows . map (fromList . map (* (-1))) . h
        sqrtInvFisherMetric x = let z = schur (invFisherMetric x)
                                in  fst z <> sqrtm (snd z) <> trans (fst z)
{-# INLINE createOptions #-}

-- | Non-isotropic Gaussian density.
nonIsoGauss :: [Double] -> [Double] -> Matrix Double -> Double
nonIsoGauss xs mu sig = exp val
  where val = -0.5*p*log (2*pi) - 0.5*ldet - 0.5*
            (trans (xsM `sub` muM) <> invSig <> (xsM `sub` muM)) @@> (0, 0)
        (xsM, muM) = (\f (a, b) -> (f a, f b)) (\l -> fromColumns [fromList l]) (xs, mu)
        p                   = fromIntegral $ cols sig
        (invSig, (ldet, _)) = invlndet sig
{-# INLINE nonIsoGauss #-}

-- | Mean function for the discretized Riemannian Langevin diffusion.
localMean :: Monad m
          => [Double]                -- Current state
          -> ViewsOptions m [Double] -- Localized mean of proposal distribution
localMean t = do
    Options _ c _ _ e <- ask
    return $! zipWith (+) t (map (* (0.5 * e^(2 :: Int))) (c t))
{-# INLINE localMean #-}

-- | Perturb the state, creating a new proposal.
perturb :: PrimMonad m
        => [Double]                   -- Current state
        -> Gen (PrimState m)          -- MWC PRNG
        -> ViewsOptions m [Double]    -- Resulting perturbation.
perturb t g = do
    Options _ _ _ s _ <- ask
    zs <- replicateM (length t) (lift $ standard g)
    t0 <- localMean t
    let adjustedBrownianMotion = s t <> fromColumns [fromList zs]
        abmList = concat . toLists . trans $ adjustedBrownianMotion
        perturbedState = zipWith (+) t0 t1
        t1 = map (* eps) abmList
    return $! perturbedState
{-# INLINE perturb #-}

-- | Perform a Metropolis accept/reject step.
metropolisStep :: PrimMonad m
               => MarkovChain                -- Current state
               -> Gen (PrimState m)          -- MWC PRNG
               -> ViewsOptions m MarkovChain -- New state
metropolisStep state g = do
    Options target _ iF _ e <- ask
    let (t0, nacc) = (theta &&& accepts) state
    zc           <- lift $ uniformR (0, 1) g
    proposal     <- perturb t0 g
    t0Mean       <- localMean t0
    proposalMean <- localMean proposal
    let mc = if   zc < acceptProb
             then (proposal, 1)
             else (t0,       0)

        acceptProb = if isNaN val then 1 else val where val = arRatio

        arRatio = exp . min 0 $
            target proposal + log (nonIsoGauss proposal t0Mean ((e^(2::Int)) `scale` iF t0))
          - target t0       - log (nonIsoGauss t0 proposalMean ((e^(2::Int)) `scale` iF proposal))

    return $! MarkovChain (fst mc) (nacc + snd mc)
{-# INLINE metropolisStep #-}

-- | Diffuse through states.
runChain :: Options         -- Options of the Markov chain.
         -> Int             -- Number of epochs to iterate the chain.
         -> Int             -- Print every nth iteration
         -> MarkovChain     -- Initial state of the Markov chain.
         -> Gen RealWorld   -- MWC PRNG
         -> IO MarkovChain  -- End state of the Markov chain, wrapped in IO.
runChain = go
  where go o n t !c g | n == 0 = return c
                      | n `rem` t /= 0 = do
                            r <- runReaderT (metropolisStep c g) o
                            go o (n - 1) t r g
                      | otherwise = do
                            r <- runReaderT (metropolisStep c g) o
                            print r
                            go o (n - 1) t r g
{-# INLINE runChain #-}

max_stars_repo_name: jtobin/lazy-langevin

max_issues_repo_name: jtobin/lazy-langevin

max_forks_repo_name: jtobin/lazy-langevin


text:
module ParseNumber where

import Data.Char (digitToInt)
import Data.Complex
import Data.Ratio
import DataTypes (LispVal(..))

import Numeric (readDec, readHex, readInt, readOct)
import ParserCombinators
       (Parser(..), (<|>), char, digit, failure, hexDigit, many1,
        octDigit, oneOf)

readBinary :: ReadS Integer
readBinary = readInt 2 (`elem` "01") digitToInt

parseNumber :: Parser LispVal
parseNumber = parseComplex <|> parseRational <|> parseFloat <|> parseInteger

--------------
-- Integer
--------------
parseInteger :: Parser LispVal
parseInteger =
  parseIntegerDecimal <|> do
    _ <- char '#'
    base <- oneOf "bdox"
    parseIntegerBase base

parseIntegerBase :: Char -> Parser LispVal
parseIntegerBase base =
  case base of
    'd' -> parseIntegerDecimal
    'o' -> parseIntegerOctal
    'x' -> parseIntegerHex
    'b' -> parseIntegerBinary
    _ -> failure "Bad integer format"

parseIntegerHelper :: ReadS Integer -> Parser Char -> Parser LispVal
parseIntegerHelper reader parser =
  (char '-' >> (parseIntegerHelper' negate)) <|>
  (char '+' >> (parseIntegerHelper' id)) <|>
  (parseIntegerHelper' id)
  where
    parseIntegerHelper' :: (Integer -> Integer) -> Parser LispVal
    parseIntegerHelper' op = Integer . op . fst . head . reader <$> many1 parser

parseIntegerDecimal :: Parser LispVal
parseIntegerDecimal = parseIntegerHelper readDec digit

parseIntegerOctal :: Parser LispVal
parseIntegerOctal = parseIntegerHelper readOct octDigit

parseIntegerHex :: Parser LispVal
parseIntegerHex = parseIntegerHelper readHex hexDigit

parseIntegerBinary :: Parser LispVal
parseIntegerBinary = parseIntegerHelper readBinary (oneOf "01")

--------------
-- Float
--------------
parseFloat :: Parser LispVal
parseFloat =
  parseFloatDecimal <|> do
    _ <- char '#'
    base <- oneOf "bdox"
    parseFloatBase base

parseFloatBase :: Char -> Parser LispVal
parseFloatBase base =
  case base of
    'd' -> parseFloatDecimal
    'o' -> parseFloatOctal
    'x' -> parseFloatHex
    'b' -> parseFloatBinary
    _ -> failure "Bad float format"

parseFloatDecimal :: Parser LispVal
parseFloatDecimal = parseFloatHelper 10 digit readDec

parseFloatOctal :: Parser LispVal
parseFloatOctal = parseFloatHelper 8 octDigit readOct

parseFloatBinary :: Parser LispVal
parseFloatBinary = parseFloatHelper 2 (oneOf "01") readBinary

parseFloatHex :: Parser LispVal
parseFloatHex = parseFloatHelper 16 hexDigit readHex

parseFloatHelper :: Int -> Parser Char -> (ReadS Integer) -> Parser LispVal
parseFloatHelper base p reader =
  (char '-' >> (parseFloat' negate)) <|> (char '+' >> parseFloat' id) <|>
  (parseFloat' id)
  where
    helper :: Integer -> Integer -> Double
    helper whole decimal =
      if decimal == 0
        then fromIntegral whole
        else let d = fromIntegral decimal :: Double
                 w = fromIntegral whole :: Double
                 b = fromIntegral base :: Double
                 e = logBase b d
                 f = fromIntegral $ floor e
                 g = d / (b ** (f + 1))
             in w + g
    parseFloat' :: (Double -> Double) -> Parser LispVal
    parseFloat' op = do
      w <- many1 p
      _ <- char '.'
      d <- many1 p
      let whole = fst . head $ reader w
          decimal = fst . head $ reader d
      return $ Float $ (op (helper whole decimal))

--------------
-- Complex
--------------
parseComplex :: Parser LispVal
parseComplex =
  parseComplexDecimal <|> do
    _ <- char '#'
    base <- oneOf "bdox"
    parseComplexBase base

parseComplexBase :: Char -> Parser LispVal
parseComplexBase base =
  case base of
    'd' -> parseComplexDecimal
    'o' -> parseComplexOctal
    'x' -> parseComplexHex
    'b' -> parseComplexBinary
    _ -> failure "Bad complex format"

parseComplexDecimal :: Parser LispVal
parseComplexDecimal =
  parseComplexHelper parseIntegerDecimal parseFloatDecimal parseRationalDecimal

parseComplexOctal :: Parser LispVal
parseComplexOctal =
  parseComplexHelper parseIntegerOctal parseFloatOctal parseRationalOctal

parseComplexHex :: Parser LispVal
parseComplexHex =
  parseComplexHelper parseIntegerHex parseFloatHex parseRationalHex

parseComplexBinary :: Parser LispVal
parseComplexBinary =
  parseComplexHelper parseIntegerBinary parseFloatBinary parseRationalBinary

parseComplexHelper ::
     Parser LispVal -> Parser LispVal -> Parser LispVal -> Parser LispVal
parseComplexHelper pn pf pr = do
  real <- fmap toDouble (pr <|> pf <|> pn)
  imaginary <- fmap toDouble (pr <|> pf <|> pn)
  _ <- char 'i'
  return $ Complex (real :+ imaginary)
  where
    toDouble (Float x) = x
    toDouble (Integer x) = fromIntegral x
    toDouble (Rational x) = fromRational x

--------------
-- Rational
--------------
parseRational :: Parser LispVal
parseRational =
  parseRationalDecimal <|> do
    _ <- char '#'
    base <- oneOf "bdox"
    parseRationalBase base

parseRationalBase :: Char -> Parser LispVal
parseRationalBase base =
  case base of
    'd' -> parseRationalDecimal
    'o' -> parseRationalOctal
    'x' -> parseRationalHex
    'b' -> parseRationalBinary
    _ -> failure "Bad rational format"

parseRationalHelper :: Parser LispVal -> Parser LispVal
parseRationalHelper p = do
  num <- fmap toInt p
  _ <- char '/'
  denom <- fmap toInt p
  return $ Rational (num % denom)
  where
    toInt (Integer x) = x

parseRationalDecimal :: Parser LispVal
parseRationalDecimal = parseRationalHelper parseIntegerDecimal

parseRationalOctal :: Parser LispVal
parseRationalOctal = parseRationalHelper parseIntegerOctal

parseRationalHex :: Parser LispVal
parseRationalHex = parseRationalHelper parseIntegerHex

parseRationalBinary :: Parser LispVal
parseRationalBinary = parseRationalHelper parseIntegerBinary

max_stars_repo_name: ebenpack/hascheme

max_issues_repo_name: ebenpack/hascheme

max_forks_repo_name: ebenpack/hascheme


text:
module Language.Hakaru.Util.Coda where

import Statistics.Autocorrelation
import qualified Data.Packed.Vector as V
import qualified Data.Vector.Generic as G

effectiveSampleSize :: [Double] -> Double
effectiveSampleSize samples = n / (1 + 2*(G.sum rho))
  where n = fromIntegral (V.dim vec)
        vec = V.fromList samples
        cov = autocovariance vec
        rho = G.map (/ G.head cov) cov

max_stars_repo_name: zaxtax/hakaru-old

max_issues_repo_name: zaxtax/hakaru-old

max_forks_repo_name: zaxtax/hakaru-old


text:
module SpaceState.Init(lostLife, initState, die)
where

import Data.List
import Control.Monad
import Control.Monad.State as State
import Prelude hiding (catch)

import Graphics.Rendering.OpenGL as OpenGL
import Graphics.UI.SDL as SDL
import qualified Data.Edison.Assoc.StandardMap as M

import OpenGLUtils
import Statistics
import Entity
import AObject
import TextScreen
import Mission

import SpaceState.Game
import SpaceState.City
import SpaceState.Common

lostLife :: String -> String -> StateT SpaceState IO Bool
lostLife s1 s2 = do
  modify $ modLives pred
  modify $ modPlCash $ const 0 -- so that no points are given for cash
  state <- State.get
  if lives state <= 0
    then do
      gameOver s1 gameoverText
    else do
      pressKeyScreen (liftIO $ makeTextScreen (30, 550) [(gamefont state,
                         Color4 1.0 0.2 0.2 1.0, s1 ++ "\n\n" ++ s2 ++ "\n\nPress ENTER to continue")]
                         (return ())) SDLK_RETURN
      initState
      return False

initState :: StateT SpaceState IO ()
initState = do
  modify $ modAObjects $ setupBarycenters
  lc <- getRandomPlanet
  modify $ modTri $ modifyPosition (const $ (getPosition lc *+* (glVector3UnitX *** (2 * AObject.size lc))))
  modify $ modPlCash $ const startCash
  modify $ modPlHoldspace $ const maxHold
  modify $ modPlCargo $ const M.empty
  modify $ modPlHealth $ const startPlHealth
  modify $ modPlMissions $ const emptyMissionMap
  modify $ modAvailMission $ const Nothing
  gotoCity lc
  releaseKeys

getRandomPlanet :: StateT SpaceState IO AObject
getRandomPlanet = do
  state <- State.get
  n <- liftIO $ chooseIO (aobjects state)
  if (aobjName n == "Star")
    then getRandomPlanet
    else return n

die :: StateT SpaceState IO (Maybe Int)
die = do
  state <- State.get
  return $ Just $ finalPoints state


text:
module NeuralNet.Layer (
  NeuronLayer (NeuronLayer),
  ForwardPropStep (),
  forwardPropA,
  forwardPropZ,
  createInputForwardPropStep,
  createForwardPropStep,
  layerW,
  layerB,
  layerActivation,
  layerNumInputs,
  layerForward,
  layerForwardSet,
  updateLayerParams
) where

import NeuralNet.Activation
import NeuralNet.Matrix
import Numeric.LinearAlgebra

data NeuronLayer = NeuronLayer Activation (Matrix Double) (Matrix Double)
  deriving (Show, Eq)

data ForwardPropStep = ForwardPropStep (Matrix Double) (Matrix Double)
  deriving (Show, Eq)

forwardPropA :: ForwardPropStep -> Matrix Double
forwardPropA (ForwardPropStep a _) = a

forwardPropZ :: ForwardPropStep -> Matrix Double
forwardPropZ (ForwardPropStep _ z) = z

createInputForwardPropStep :: Matrix Double -> ForwardPropStep
createInputForwardPropStep a = createForwardPropStep a (zero (rows a) (cols a))

createForwardPropStep :: Matrix Double -> Matrix Double -> ForwardPropStep
createForwardPropStep a z
  | rows a /= rows z = error "a and z need same num rows"
  | cols a /= cols z = error "a and z need same num cols"
  | otherwise          = ForwardPropStep a z

layerForward :: NeuronLayer -> [Double] -> [Double]
layerForward l inputs = concat (toLists (takeColumns 1 (forwardPropA (layerForwardSet l x))))
  where x = col inputs

layerForwardSet :: NeuronLayer -> Matrix Double -> ForwardPropStep
layerForwardSet l x
  | expInputs /= n = error ("Layer " ++ show l ++ " expected " ++ show expInputs ++ " inputs but got " ++ show n)
  | otherwise      = createForwardPropStep a z
    where
      expInputs = layerNumInputs l
      n = rows x
      w = layerW l
      b = layerB l
      z = w <> x + b
      a = forward (layerActivation l) z

layerW :: NeuronLayer -> Matrix Double
layerW (NeuronLayer _ w _) = w

layerNumInputs :: NeuronLayer -> Int
layerNumInputs = cols . layerW

layerB :: NeuronLayer -> Matrix Double
layerB (NeuronLayer _ _ b) = b

layerActivation :: NeuronLayer -> Activation
layerActivation (NeuronLayer a _ _) = a

updateLayerParams :: NeuronLayer -> Matrix Double -> Matrix Double -> NeuronLayer
updateLayerParams l newW newB
  | dims (layerW l) /= dims newW = error "Ws don't align"
  | dims (layerB l) /= dims newB = error "bs don't align"
  | otherwise                    = NeuronLayer (layerActivation l) newW newB

max_stars_repo_name: danielholmes/neural-net

max_issues_repo_name: danielholmes/neural-net

max_forks_repo_name: danielholmes/neural-net