text:
-- {-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE DataKinds, GADTs, TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FunctionalDependencies #-}

module HBLAS.BLAS.Level1Spec (main, spec) where

import Data.Complex

import Numerical.HBLAS.MatrixTypes as Matrix
import Numerical.HBLAS.BLAS.Level1 as BLAS

import Test.Hspec

main :: IO ()
main = hspec spec

spec :: Spec
spec = do
  asumSpec
  axpySpec
  copySpec
  dotSpec
  sdotSpec
  dotcSpec
  dotuSpec
  nrm2Spec
  rotSpec
  rotgSpec
  rotmSpec
  rotmgSpec
  scalSpec
  swapSpec
  iamaxSpec
  iaminSpec

asumSpec :: Spec
asumSpec =
  context "?ASUM" $ do
    describe "SASUM" $ do
      it "vector of length 6 with incx 1" $ do
        vecTest1SASUM
      it "vector of length 12 with incx 2" $ do
        vecTest2SASUM

vecTest1SASUM :: IO ()
vecTest1SASUM = do
  vec <- Matrix.generateMutableDenseVector 6 (\idx -> [1 .. 6] !! idx)
  res <- BLAS.sasum vec
  res `shouldBe` (21.0 :: Float)

vecTest2SASUM :: IO ()
vecTest2SASUM = do
  vec <- Matrix.generateMutableDenseVectorWithStride 12 2 (\idx -> [1 .. 12] !! idx)
  res <- BLAS.sasum  vec
  res `shouldBe` (sum [1,3 .. 12] :: Float)

axpySpec :: Spec
axpySpec =
  context "?AXPY" $ do
    describe "SAXPY" $ do
      it "vectors of lengths 6 and 6 with both incx 1" $ do
        vecTest1SAXPY
      it "vectors of lengths 12 and 18 with incx 2 and 3" $ do
        vecTest2SAXPY

vecTest1SAXPY :: IO ()
vecTest1SAXPY = do
  input <- Matrix.generateMutableDenseVector 6 (\idx -> [1 .. 6] !! idx)
  output <- Matrix.generateMutableDenseVector 6 (\idx -> [2, 3, 4, 3, 5, 6] !! idx)
  BLAS.saxpy 6 (-1.0) input output
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector output
  resList `shouldBe` [1, 1, 1, -1, 0, 0]

vecTest2SAXPY :: IO ()
vecTest2SAXPY = do
  input <- Matrix.generateMutableDenseVectorWithStride 18 3 (\idx -> [1 .. 18] !! idx)
  output <- Matrix.generateMutableDenseVectorWithStride 12 2 (\idx -> [1 .. 12] !! idx)
  BLAS.saxpy 6 2.0 input output
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector output
  resList `shouldBe` [3, 2, 11, 4, 19, 6, 27, 8, 35, 10, 43, 12]

copySpec :: Spec
copySpec =
  context "?COPY" $ do
    describe "DCOPY" $ do
      it "vectors of lengths 6 and 6 with both incx 1" $ do
        vecTest1DCOPY
      it "vectors of lengths 6 and 9 with incx 2 and 3" $ do
        vecTest2DCOPY

vecTest1DCOPY :: IO ()
vecTest1DCOPY = do
  input <- Matrix.generateMutableDenseVector 6 (\idx -> [1 .. 6] !! idx)
  output <- Matrix.generateMutableDenseVector 6 (const 0.0)
  BLAS.dcopy 6 input output
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector output
  resList `shouldBe` [1, 2, 3, 4, 5, 6]

vecTest2DCOPY :: IO ()
vecTest2DCOPY = do
  input <- Matrix.generateMutableDenseVectorWithStride 6 2 (\idx -> [1 .. 6] !! idx)
  output <- Matrix.generateMutableDenseVectorWithStride 9 3 (const 0.0)
  BLAS.dcopy 3 input output
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector output
  resList `shouldBe` [1, 0, 0, 3, 0, 0, 5, 0, 0]

dotSpec :: Spec
dotSpec =
  context "?DOT" $ do
    describe "SDOT" $ do
      it "vectors of lengths 6 and 12 with incx 2 and 4" $ do
        vecTest1SDOT
    describe "DDOT" $ do
      it "vectors of lengths 12 and 6 with incx 2 and 1" $ do
        vecTest1DDOT

vecTest1SDOT :: IO ()
vecTest1SDOT = do
  left <- Matrix.generateMutableDenseVectorWithStride 6 2 (\idx -> [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] !! idx)
  right <- Matrix.generateMutableDenseVectorWithStride 12 4 (\idx -> [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0] !! idx)
  res <- sdot 3 left right
  res `shouldBe` ((1 + 15 + 45) :: Float)

vecTest1DDOT :: IO ()
vecTest1DDOT = do
  left <- Matrix.generateMutableDenseVectorWithStride 12 2 (((+) 1) . fromRational . toRational)
  right <- Matrix.generateMutableDenseVectorWithStride 6 1 (\idx -> [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] !! idx)
  res <- ddot 6 left right
  res `shouldBe` ((1 + 6 + 15 + 28 + 45 + 66) :: Double)

sdotSpec :: Spec
sdotSpec =
  context "?SDOT" $ do
    describe "SDSDOT" $ do
      it "vectors of lengths 6 and 12 with incx 2 and 4" $ do
        pending
    describe "DSDOT" $ do
      it "vectors of lengths 12 and 6 with incx 2 and 1" $ do
        vecTest1DSDOT

{-vecTest1SDSDOT :: IO ()-}
{-vecTest1SDSDOT = do-}
  {-left <- Matrix.generateMutableDenseVectorWithStride 6 2 (\idx -> ([1 .. 6] :: [Float]) !! idx) -- Float-}
  {-right <- Matrix.generateMutableDenseVectorWithStride 12 4 (\idx -> ([1 .. 12] :: [Float]) !! idx) -- Float-}
  {-res <- sdsdot 3 2.0 left right-}
  {-res `shouldBe` ((2 :: Float) + 1 + 15 + 45)-}

vecTest1DSDOT :: IO ()
vecTest1DSDOT = do
  left <- Matrix.generateMutableDenseVectorWithStride 12 2 (\idx -> ([1 .. 12] :: [Float]) !! idx)
  right <- Matrix.generateMutableDenseVectorWithStride 6 1 (\idx -> ([1 .. 6] :: [Float]) !! idx)
  res <- dsdot 6 left right
  res `shouldBe` ((1 + 6 + 15 + 28 + 45 + 66) :: Double)

dotcSpec :: Spec
dotcSpec =
  context "?DOTC" $ do
    describe "CDOTC" $ do
      it "vectors of length 6 and 9 with incx 2 and 3" $ do
        vecTest1CDOTC

vecTest1CDOTC :: IO ()
vecTest1CDOTC = do
  left <- Matrix.generateMutableDenseVectorWithStride 6 2 (\idx -> [2:+3, 1:+(-1), 1:+1, 1:+(-1), 1:+1, 1:+(-1)] !! idx)
  right <- Matrix.generateMutableDenseVectorWithStride 9 3 (\idx -> [1:+(-2), 1:+1, 1:+(-1), 1:+1, 1:+(-1), 1:+1, 1:+(-1), 1:+1, 1:+(-1)] !! idx)
  cdotc 3 left right `shouldReturn` ((-2):+(-9) :: Complex Float)

dotuSpec :: Spec
dotuSpec =
  context "?DOTU" $ do
    describe "CDOTU" $ do
      it "vectors of length 6 and 9 with incx of 2 and 3" $ do
        vecTest1CDOTU

vecTest1CDOTU :: IO ()
vecTest1CDOTU = do
  left <- Matrix.generateMutableDenseVectorWithStride 6 2 (\idx -> [1:+1, 1:+(-1), 1:+1, 1:+(-1), 1:+1, 1:+(-1)] !! idx)
  right <- Matrix.generateMutableDenseVectorWithStride 9 3 (\idx -> [1:+(-2), 1:+1, 1:+(-1), 1:+1, 1:+(-1), 1:+1, 1:+(-1), 1:+1, 1:+(-1)] !! idx)
  --res <- Matrix.generateMutableValue (1:+1)
  res <- cdotu 3 left right

  res `shouldBe` 5:+1

nrm2Spec :: Spec
nrm2Spec =
  context "?NRM2" $ do
    describe "SNRM2" $ do
      it "vector of length 6 with incx of 1" $ do
        vecTest1SNRM2
    describe "DZNRM2" $ do
      it "vector of length 8 with incx of 2" $ do
        vecTest1DZNRM2

vecTest1SNRM2 :: IO ()
vecTest1SNRM2 = do
  input <- Matrix.generateMutableDenseVector 6 (\idx -> [1.0, -2.0, 3.0, -4.0, 5.0, -6.0] !! idx)
  res <- snrm2 6 input
  True `shouldBe` 1e-6 > (abs $ res - (sqrt $ sum $ fmap (\x->x^2) [1, 2, 3, 4, 5, 6]))

vecTest1DZNRM2 :: IO ()
vecTest1DZNRM2 = do
  input <- Matrix.generateMutableDenseVectorWithStride 8 2 (\idx -> [1:+1, 1:+2, 2:+(-3), 2:+(-2), (-3):+1, (-3):+0, (-4):+2, (-4):+1] !! idx)
  res <- dznrm2 4 input
  True `shouldBe` 1e-12 > (abs $ res - (sqrt $ sum $ fmap (\x->x^2) [1, 1, 2, 3, 3, 1, 4, 2]))

rotSpec :: Spec
rotSpec =
  context "?ROT" $ do
    describe "SROT" $ do
      it "vectors of length 6 and 6 with incx of 2" $ do
        vecTest1SROT
    describe "DROT" $ do
      it "vectors of length 4 and 8 with incx of 1 and 2" $ do
        vecTest1DROT

vecTest1SROT :: IO ()
vecTest1SROT = do
  left <- Matrix.generateMutableDenseVectorWithStride 6 2 (\idx -> [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] !! idx)
  right <- Matrix.generateMutableDenseVectorWithStride 6 2 (\idx -> [6.0, 5.0, 4.0, 3.0, 2.0, 1.0] !! idx)
  srot 3 left right (-1) 2
  resLeft <- Matrix.mutableVectorToList $ _bufferMutDenseVector left
  resRight <- Matrix.mutableVectorToList $ _bufferMutDenseVector right
  resLeft `shouldBe` [11.0, 2.0, 5.0, 4.0, -1.0, 6.0]
  resRight `shouldBe` [-8.0, 5.0, -10.0, 3.0, -12.0, 1.0]

vecTest1DROT :: IO ()
vecTest1DROT = do
  left <- Matrix.generateMutableDenseVectorWithStride 4 1 (\idx -> [1, 2, 3, 4] !! idx)
  right <- Matrix.generateMutableDenseVectorWithStride 8 2 (\idx -> [8, 7, 6, 5, 4, 3, 2, 1] !! idx)
  drot 4 left right 0 (-2)
  resLeft <- Matrix.mutableVectorToList $ _bufferMutDenseVector left
  resRight <- Matrix.mutableVectorToList $ _bufferMutDenseVector right
  resLeft `shouldBe` [-16, -12, -8, -4]
  resRight `shouldBe` [2, 7, 4, 5, 6, 3, 8, 1]

rotgSpec :: Spec
rotgSpec =
  context "?ROTG" $ do
    describe "SROTG" $ do
      it "apply on 3 and 4" $ do
        pending
    describe "DROTG" $ do
      it "apply on 5.8 and 3.4" $ do
        pending

{-vecTest1SROTG :: IO ()-}
{-vecTest1SROTG = do-}
  {-a <- Matrix.generateMutableValue 3-}
  {-b <- Matrix.generateMutableValue 4-}
  {-c <- Matrix.generateMutableValue 0-}
  {-s <- Matrix.generateMutableValue 0-}
  {-srotg a b c s-}
  {-av <- Matrix.mutableValueToValue a-}
  {-bv <- Matrix.mutableValueToValue b-}
  {-cv <- Matrix.mutableValueToValue c-}
  {-sv <- Matrix.mutableValueToValue s-}
  {-av `shouldBe` 5-}
  {-True `shouldBe` 1e-6 > (abs $ bv - 1/0.6)-}
  {-cv `shouldBe` 0.6-}
  {-sv `shouldBe` 0.8-}

{-vecTest1DROTG :: IO ()-}
{-vecTest1DROTG = do-}
  {-a <- Matrix.generateMutableValue 5.8-}
  {-b <- Matrix.generateMutableValue 3.4-}
  {-c <- Matrix.generateMutableValue 0-}
  {-s <- Matrix.generateMutableValue 0-}
  {-drotg a b c s-}
  {-av <- Matrix.mutableValueToValue a-}
  {-bv <- Matrix.mutableValueToValue b-}
  {-cv <- Matrix.mutableValueToValue c-}
  {-sv <- Matrix.mutableValueToValue s-}
  {-True `shouldBe` 1e-12 > (abs $ av - sqrt(3.4^2 + 5.8^2))-}
  {-True `shouldBe` 1e-12 > (abs $ bv - 3.4 / sqrt(3.4^2 + 5.8^2))-}
  {-True `shouldBe` 1e-12 > (abs $ cv - 5.8 / sqrt(3.4^2 + 5.8^2))-}
  {-True `shouldBe` 1e-12 > (abs $ sv - 3.4 / sqrt(3.4^2 + 5.8^2))-}

rotmSpec :: Spec
rotmSpec =
  context "?ROTM" $ do
    it "vectors of length 4 and 8 with incx of 1 and 2, param starts with -1" $ do
      pending
    it "vectors of length 6 and 9 with incx of 2 and 3, param starts with 1" $ do
      pending

{-
vecTest1DROTM :: IO ()
vecTest1DROTM = do
  x <- Matrix.generateMutableDenseVectorWithStride 4 1 (\idx -> [1, 2, 3, 4] !! idx)
  y <- Matrix.generateMutableDenseVectorWithStride 8 2 (\idx -> [8, 7, 6, 5, 4, 3, 2, 1] !! idx)
  param <- Matrix.generateMutableDenseVector 5 (\idx -> [-1, 0, -1, 1, 0] !! idx)
  drotm 4 x y param
  resX <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
  resY <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resX `shouldBe` [8, 6, 4, 2]
  resY `shouldBe` [-1, 7, -2, 5, -3, 3, -4, 1]

vecTest1SROTM :: IO ()
vecTest1SROTM = do
  x <- Matrix.generateMutableDenseVectorWithStride 6 2 (\idx -> [1, 2, 3, 4, 5, 6] !! idx)
  y <- Matrix.generateMutableDenseVectorWithStride 9 3 (\idx -> [9, 8, 7, 6, 5, 4, 3, 2, 1] !! idx)
  param <- Matrix.generateMutableDenseVector 5 (\idx -> [1, 1, 2, -2, 1] !! idx)
  srotm 3 x y param
  resX <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
  resY <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resX `shouldBe` [10, 2, 9, 4, 8, 6]
  resY `shouldBe` [8, 8, 7, 3, 5, 4, -2, 2, 1]
-}

rotmgSpec :: Spec
rotmgSpec =
  context "?ROTMG" $ do
    describe "SROTMG" $ do
      it "todo" $ do
        pending

{-
vecTest1SROTMG :: IO ()
vecTest1SROTMG = do
  d1 <- Matrix.generateMutableValue 3
  d2 <- Matrix.generateMutableValue 6
  x <- Matrix.generateMutableValue 1
  let y = 1
  param <- Matrix.generateMutableDenseVector 5 (\idx -> [-1, 1, 1, -1, 1] !! idx)
  srotmg d1 d2 x y param
  paramR <- Matrix.mutableVectorToList $ _bufferMutDenseVector param
  updatedD1 <- Matrix.mutableValueToValue d1
  updatedD2 <- Matrix.mutableValueToValue d2
  updatedX <- Matrix.mutableValueToValue x
  paramR `shouldBe` [1, 0, 0.5, 0, 1]
  updatedD1 `shouldBe` 4
  updatedD2 `shouldBe` 2
  updatedX `shouldBe` 1.5
-}

scalSpec :: Spec
scalSpec =
  context "?SCAL" $ do
    describe "SSCAL" $ do
      it "vector of length 8 with incx 2" $ do
        vecTest1SSCAL
    describe "CSCAL" $ do
      it "vector of length 8 with incx 4" $ do
        vecTest1CSCAL
    describe "CSSCAL" $ do
      it "vector of length 8 with incx 1" $ do
        vecTest1CSSCAL

vecTest1SSCAL :: IO ()
vecTest1SSCAL = do
  x <- Matrix.generateMutableDenseVectorWithStride 8 2 (\idx -> [1, 2, 3, 4, 5, 6, 7, 8] !! idx)
  sscal 4 (-2) x
  xRes <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
  xRes `shouldBe` [-2, 2, -6, 4, -10, 6, -14, 8]

vecTest1CSCAL :: IO ()
vecTest1CSCAL = do
  x <- Matrix.generateMutableDenseVectorWithStride 8 4 (\idx -> [1:+1, 1:+2, 2:+(-3), 2:+(-2), (-3):+1, (-3):+0, (-4):+2, (-4):+1] !! idx)
  cscal 2 (2:+(-2)) x -- size 2?
  xRes <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
  xRes `shouldBe` [4:+0, 1:+2, 2:+(-3), 2:+(-2), (-4):+8, (-3):+0, (-4):+2, (-4):+1]

vecTest1CSSCAL :: IO ()
vecTest1CSSCAL = do
  x <- Matrix.generateMutableDenseVector 8 (\idx -> [1:+1, 1:+2, 2:+(-3), 2:+(-2), (-3):+1, (-3):+0, (-4):+2, (-4):+1] !! idx)
  csscal 8 (-2) x
  xRes <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
  xRes `shouldBe` [(-2):+(-2), (-2):+(-4), (-4):+6, (-4):+4, 6:+(-2), 6:+0, 8:+(-4), 8:+(-2)]

swapSpec :: Spec
swapSpec =
  context "?SWAP" $ do
    describe "SSWAP" $ do
      it "vectors of length 8 and 4 with incx 2 and 1" $ do
        vecTest1SSWAP
    describe "CSWAP" $ do
      it "vectors of length 9 and 6 with incx 3 and 2" $ do
        vecTest1CSWAP

vecTest1SSWAP :: IO ()
vecTest1SSWAP = do
  x <- Matrix.generateMutableDenseVectorWithStride 8 2 (\idx -> [1, 2, 3, 4, 5, 6, 7, 8] !! idx)
  y <- Matrix.generateMutableDenseVectorWithStride 4 1 (\idx -> [-1, -2, -3, -4] !! idx)
  sswap 4 x y
  xRes <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
  yRes <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  xRes `shouldBe` [-1, 2, -2, 4, -3, 6, -4, 8]
  yRes `shouldBe` [1, 3, 5, 7]

vecTest1CSWAP :: IO ()
vecTest1CSWAP = do
  x <- Matrix.generateMutableDenseVectorWithStride 9 3 (\idx -> [1:+1, 1:+2, 2:+(-3), 2:+(-2), (-3):+1, (-3):+0, (-4):+2, (-4):+1, 0:+9] !! idx)
  y <- Matrix.generateMutableDenseVectorWithStride 6 2 (\idx -> [1:+2, 1:+3, 3:+(-3), 2:+2, 3:+1, 3:+3] !! idx)
  cswap 3 x y
  xRes <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
  yRes <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  xRes `shouldBe` [1:+2, 1:+2, 2:+(-3), 3:+(-3), (-3):+1, (-3):+0, 3:+1, (-4):+1, 0:+9]
  yRes `shouldBe` [1:+1, 1:+3, 2:+(-2), 2:+2, (-4):+2, 3:+3]

iamaxSpec :: Spec
iamaxSpec =
  context "I?AMAX" $ do
    describe "ISAMAX" $ do
      it "vector of length 8 with incx 2" $ do
        vecTest1ISAMAX
    describe "ICAMAX" $ do
      it "vector of length 9 with incx 1" $ do
        vecTest1ICAMAX

vecTest1ISAMAX :: IO ()
vecTest1ISAMAX = do
  x <- Matrix.generateMutableDenseVectorWithStride 8 2 (\idx -> [1, 2, 3, 4, 5, 6, 7, 8] !! idx)
  idx <- isamax  x
  idx `shouldBe` 7

vecTest1ICAMAX :: IO ()
vecTest1ICAMAX = do
  x <- Matrix.generateMutableDenseVector 9 (\idx -> [1:+1, 1:+2, 2:+(-3), 2:+(-2), (-3):+1, (-3):+0, (-4):+2, (-4):+1, 0:+9] !! idx)
  idx <- icamax  x
  idx `shouldBe` 8

iaminSpec :: Spec
iaminSpec =
  context "I?AMIN" $ do
    describe "ISAMIN" $ do
      it "vector of length 8 with incx 2" $ do
        pending
    describe "ICAMIN" $ do
      it "vector of length 9 with incx 1" $ do
        pending

{-
vecTest1ISAMIN :: IO ()
vecTest1ISAMIN = do
  x <- Matrix.generateMutableDenseVector 8 (\idx -> [1, 2, 3, 4, -5, 6, 7, 8] !! idx)
  idx <- isamin 4 x 2
  idx `shouldBe` 2

vecTest1ICAMIN :: IO ()
vecTest1ICAMIN = do
  x <- Matrix.generateMutableDenseVector 9 (\idx -> [1:+2, 1:+2, (-2):+(-3), 2:+(-2), (-3):+1, (-2):+0, (-4):+2, (-4):+1, 0:+9] !! idx)
  idx <- icamin 9 x 1
  idx `shouldBe` 5
-}

max_stars_repo_name: schnecki/hblas

max_issues_repo_name: schnecki/hblas

max_forks_repo_name: schnecki/hblas


text:
module RC.Helpers
  ( addBiases
  , flatten'
  , unflatten'
  , randList
  , randMatrix
  , randSparse
  , hsigmoid
  ) where

import           System.Random
import           Data.List ( unfoldr )
import           Numeric.LinearAlgebra as LA

-- | Hard sigmoid
hsigmoid :: (Fractional a, Ord a)
         => (a, a, a)
         -- ^ Vertical scaling, width, offset
         -> a
         -> a
hsigmoid (β, width, offset) x = f
  where
    f | x < offset = 0.0
      | x < width = β * (x - offset)
      | otherwise = β * (width - offset)
{-# SPECIALISE hsigmoid :: (Double, Double, Double) -> Double -> Double #-}
{-# SPECIALISE hsigmoid :: (Float, Float, Float) -> Float -> Float #-}

-- | Prepend a row of ones
--
-- >>> addBiases $ (2><3) [20..26]
-- (3><3)
--  [  1.0,  1.0,  1.0
--  , 20.0, 21.0, 22.0
--  , 23.0, 24.0, 25.0 ]
addBiases :: LA.Matrix Double -> LA.Matrix Double
addBiases m = let no = LA.cols m
                  m' = LA.konst 1.0 (1, no)
              in m' LA.=== m

-- | Random matrix with elements in the range of [minVal; maxVal]
randMatrix
  :: StdGen
     -> (Int, Int)
     -- ^ Number of rows and columns
     -> (Double, Double)
     -- ^ Minimal and maximal values
     -> LA.Matrix Double
randMatrix seed (rows', cols') (minVal, maxVal) = (LA.reshape cols'. LA.vector) xs
  where
    xs = f <$> randList (rows' * cols') seed
    f x = (maxVal - minVal) * x + minVal

-- | Random sparse matrix
--
-- NB: at the moment, the matrix is stored in memory as
-- an ordinary (dense) matrix.
randSparse g (rows', cols') (minVal, maxVal) connectivity =
    LA.reshape cols' $ LA.vector xs
  where
    (g1, g2) = split g
    rlist = randList (rows' * cols')
    xs = zipWith f (rlist g1) (rlist g2)
    f lv rv | lv < connectivity = (maxVal - minVal) * rv + minVal
            | otherwise = 0.0

randList :: (Random a, Floating a) => Int -> StdGen -> [a]
randList n = take n. unfoldr (Just. random)
{-# SPECIALISE randList :: Int -> StdGen -> [Float] #-}
{-# SPECIALISE randList :: Int -> StdGen -> [Double] #-}

-- | Flatten by columns
flatten' :: Matrix Double -> Vector Double
flatten' = flatten. tr

-- | Reshape using columns instead of rows
unflatten' :: Int -> Vector Double -> Matrix Double
unflatten' nodes = tr. reshape nodes

ext: hs

max_stars_repo_name: masterdezign/rc

max_issues_repo_name: masterdezign/rc

max_forks_repo_name: masterdezign/rc


text:
{- |
Module      :  Numeric.GSL.ODE
Copyright   :  (c) Alberto Ruiz 2010
License     :  GPL

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

Solution of ordinary differential equation (ODE) initial value problems.

<http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>

A simple example:

@import Numeric.GSL
import Numeric.LinearAlgebra
import Graphics.Plot

xdot t [x,v] = [v, -0.95*x - 0.1*v]

ts = linspace 100 (0,20 :: Double)

sol = odeSolve xdot [10,0] ts

main = mplot (ts : toColumns sol)@

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

module Numeric.GSL.ODE (
    odeSolve, odeSolveV, ODEMethod(..), Jacobian
) where

import Data.Packed.Internal
import Numeric.GSL.Internal

import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
import Foreign.C.Types
import System.IO.Unsafe(unsafePerformIO)

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

type Jacobian = Double -> Vector Double -> Matrix Double

-- | Stepping functions
data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
               | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use the embedded methods.
               | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
               | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
               | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
               | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
               | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
               | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.
               | RK1imp Jacobian -- ^ Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method.
               | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12.
               | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems.

-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
odeSolve
    :: (Double -> [Double] -> [Double])        -- ^ xdot(t,x)
    -> [Double]        -- ^ initial conditions
    -> Vector Double   -- ^ desired solution times
    -> Matrix Double   -- ^ solution
odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
    where hi = (ts@>1 - ts@>0)/100
          epsAbs = 1.49012e-08
          epsRel = 1.49012e-08
          l2v f = \t -> fromList  . f t . toList

-- | Evolution of the system with adaptive step-size control.
odeSolveV
    :: ODEMethod
    -> Double -- ^ initial step size
    -> Double -- ^ absolute tolerance for the state vector
    -> Double -- ^ relative tolerance for the state vector
    -> (Double -> Vector Double -> Vector Double)   -- ^ xdot(t,x)
    -> Vector Double     -- ^ initial conditions
    -> Vector Double     -- ^ desired solution times
    -> Matrix Double     -- ^ solution
odeSolveV RK2 = odeSolveV' 0 Nothing
odeSolveV RK4 = odeSolveV' 1 Nothing
odeSolveV RKf45 = odeSolveV' 2 Nothing
odeSolveV RKck = odeSolveV' 3 Nothing
odeSolveV RK8pd = odeSolveV' 4 Nothing
odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac)
odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac)
odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac)
odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
odeSolveV MSAdams = odeSolveV' 9 Nothing
odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)

odeSolveV'
    :: CInt
    -> Maybe (Double -> Vector Double -> Matrix Double)   -- ^ optional jacobian
    -> Double -- ^ initial step size
    -> Double -- ^ absolute tolerance for the state vector
    -> Double -- ^ relative tolerance for the state vector
    -> (Double -> Vector Double -> Vector Double)   -- ^ xdot(t,x)
    -> Vector Double     -- ^ initial conditions
    -> Vector Double     -- ^ desired solution times
    -> Matrix Double     -- ^ solution
odeSolveV' method mbjac h epsAbs epsRel f  xiv ts = unsafePerformIO $ do
    let n   = dim xiv
    fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
    jp <- case mbjac of
        Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
        Nothing  -> return nullFunPtr
    sol <- vec xiv $ \xiv' ->
            vec (checkTimes ts) $ \ts' ->
             createMIO (dim ts) n
              (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
              "ode"
    freeHaskellFunPtr fp
    return sol

foreign import ccall safe "ode"
    ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM

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

checkdim1 n v
    | dim v == n = v
    | otherwise = error $ "Error: "++ show n
                        ++ " components expected in the result of the function supplied to odeSolve"

checkdim2 n m
    | rows m == n && cols m == n = m
    | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
                        ++ " Jacobian expected in odeSolve"

checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
              | otherwise = error "odeSolve requires increasing times"
    where ts' = toList ts

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{- MDP.hs       Code for CS 533, Mini-Project 2.
 -
 - Copyright © 2012 Loren B. Davis.  ALL RIGHTS RESERVED.  I hereby release
 - this code under the MIT License: (http://www.opensource.org/licenses/MIT).
 - Since this is an assigned exercise, remember your academic integrity!
 -}

module MDP {- ( MarkovDP, State, Action, policyIterate, showUpdates, test4x3,
             showT, decodeT, whenToPass, whenToPark, parkingProblem,
             ParkingState, encodeState, decodeState ) -}
  where
import Data.Array.Unboxed (UArray, (!), accumArray, assocs, listArray)
import Data.List (maximumBy, (\\))
import Data.Ord ()
import Numeric.LinearAlgebra (Matrix, Vector, build, linearSolveSVD, toLists)
-- import Numeric.LinearAlgebra hiding (i)
-- import Numeric.LinearAlgebra.Algorithms hiding (i)
import Prelude hiding(pi)
{- This program uses the hmatrix package (cabal install hmatrix) rather than
 - attempt to reinvent the wheel.  You might also need to install libgsl,
 - liblapack, and possibly MKL or ATLAS.
 -}

type State = Int
type Action = Int
type MarkovDP = ( Int,
                  Int,
                  UArray State Double,
                  UArray (State,Action,State) Double
                )
{- n = The number of states.  S = {0,...,n-1}
 - m = The number of actions. A = {0,...,m-1}
 - r = An array of size n representing the reward function.  r[s] is the
 -     reward associated with state s.
 - t = A three-dimensional array representing the transition function, such
 -     that T(s,a,s') = t ! (s,a,s').  Note that the values are between 0 and
 -     1, each row sums to 1, and that the matrices will typically be sparse.
 -}
type PIUpdate = [(State, Action)]
{- The update of each step of the policy iteration algorithm is the list
 - of pairs (s,a) such that the new policy will take action a in state s.
 -}

policyIterate :: MarkovDP ->
                 Double ->
                 (UArray State Double, UArray State Action, [String])
{- Parameters:
 - process  = The Markov decision process to solve, in the form
 -        (|S|,|A|,R,T)
 - gamma    = The discount parameter, between 0 and 1.
 -
 - Returns (v, pi, updates) where
 - v        = A size-n array representing the value function.  The expect-
 -        ed value of being in state s is v[s].
 - pi       = A size-n array representing the policy.  The action that pi
 -        will take in state s is pi[s].
 - updates  = A human-readable list of update strings, one per step of the
 -        iteration.
 -
 - This implementation always starts with a policy of taking the first action
 - in any situation.  This is not really correct (the algorithm in the text
 - says to try a "random" policy), but saves me from dealing with stateful
 - computation in this program.  To specify a different initial policy, call
 - policyIterate1.  It does not attempt any asynchronous policy iteration.
 -}
policyIterate (n, m, r, t) gamma = (v', pi'', translatePIUpdates updates'' )
  where
    ( v', pi'', updates'' ) = seq ((n,m,r,t), gamma, naive)
                                  (policyIterate1 naive [])
    naive = (initialPolicy n)
    policyIterate1 :: UArray State Action ->
                      [PIUpdate] ->
                      (UArray State Double, UArray State Action, [PIUpdate])
{- Parameters:
 - pi       = The initial policy to attempt to improve.
 - updates  = The previous iterations' updates.
 -
 - Returns (v, pi) where
 - v    = A size-n array representing the value function.  The expected value
 -     of being in state s is v[s].
 - pi'  = A size-n array representing the policy.  The action that pi' will
 -       take in state s is pi[s].
 - updates' = The list of updates for all iterations of the algorithm.  When
 -            there are no further updates, the algorithm terminates.
 -}
    policyIterate1 pi updates | diff == [] = (v, pi', updates')
                              | otherwise  = policyIterate1 pi' $! updates'
      where
        states :: [State]
        states = [0..n-1]
        actions :: [Action]
        actions = [0..m-1]
        v :: UArray State Double
{- The updated value function.  This is based on eq. 17.10 in (Russell &
 - Norvig 2010), rather than the modified-policy-iteration algorithm in fig.
 - 17.7.  It uses LAPACK to solve the system of linear equations in O(n^3)
 - time.
 -}
        v = let
              updateHelper :: Double -> Double -> Double
              updateHelper w z = gamma * (t!(i,pi!i,j)) -
                                           (if i == j then 1.0 else 0.0)
                where
                  i :: Int
                  i = truncate w
                  j :: Int
                  j = truncate z
              a :: Matrix Double
              a = build (n,n) updateHelper
              b :: Matrix Double
              b = build (n,1) (\y _ -> let i = truncate y :: Int in - (r!i))
              x :: Matrix Double
              x = linearSolveSVD a b
            in
              listArray (0,n-1) (concat . toLists $ x)
        pi' :: UArray Int Action
{- The updated policy pi' is the best action to take in any state (determined
 - by one-step lookahead).  Q.v. eq. 17.4 in (Russell & Norvig 2010).
 -}
        pi' = seq v (listArray (0, n-1)
                        (map (\i -> maximumBy (argmax i) actions)
                             states))
        argmax :: State -> Action -> Action -> Ordering
{- A helper function that compares the expected value of actions a1 and a2
 - from state s0.
 -}
        argmax s0 a1 a2 = let
                            u1 = sum .
                                 map (\s' -> t!(s0,a1,s') * v!s') $
                                 states
                            u2 = sum .
                                 map (\s' -> t!(s0,a2,s') * v!s') $
                                 states
                          in compare u1 u2
        diff :: PIUpdate
        diff = seq pi' ((assocs pi') \\ (assocs pi))
        updates' :: [PIUpdate]
        updates' = seq diff (updates ++ [diff])

initialPolicy :: Int -> UArray Int Action
{- A policy for a MDP with n states that always takes action zero.  This is
 - not the same as the "random" policy I'm supposed to use for the algorithm,
 - but is simpler to write in this language.
 -}
initialPolicy n = listArray (0, n-1) (take n (repeat 0))

initialVF :: Int -> UArray Int Double
{- A value function initialized to n zeroes.  Unlike the above, this is corr-
 - ect.  (But unused in policy iteration.)
 -}
initialVF n = listArray (0, n-1) (take n (repeat 0.0))

translatePIUpdates :: [PIUpdate] -> [String]
translatePIUpdates u =
  map (\i -> "Step " ++ show (i+1) ++ ": " ++ translatePIUpdate (u !! i))
      [0 .. (length u) - 1]

translatePIUpdate :: PIUpdate -> String
translatePIUpdate [] = "Done."
translatePIUpdate [(s,a)] = "take action " ++
                            show a ++
                            " in state " ++
                            show s ++
                            "."
translatePIUpdate ((s,a):xs) = "take action " ++
                               show a ++
                               " in state " ++
                               show s ++
                               "; " ++
                               translatePIUpdate xs

-- Display the returned update messages:
showUpdates :: (UArray State Double, UArray State Action, [String]) ->
               IO ()
showUpdates (_,_,u) = putStrLn . unlines $ u

-- A useful means of displaying only the nonzero entries of the sparse array
-- t:
showT :: MarkovDP -> [((State,Action,State),Double)]
showT (_,_,_,t) = filter ((> 0.0) . snd) (assocs t)

-- The parking-lot world:

{- The decoded values of a state of the parking-lot world.  Note that Exit
 - in this implementation is a terminal state that the park action automat-
 - ically reaches, not a third action.  Crash is another world state, always
 - and only reached by attempting to park in a filled spot.
 -
 - A i True is the state in which the agent is in front of row A, space i,
 - which is filled.  A i False is similar, but the spot is empty; likewise
 - for row B.  I take advantage of the symmetry between parking in row A and
 - in row B to define P i as the state of being parked i spaces from the
 - store.  This reduces the size of the state space from 8n to 5n+1. The
 - range of i is 1 to n.
 -}
data ParkingState = A Int Bool | B Int Bool | P Int | Exit | Crash
  deriving (Eq, Show)

encodeState :: Int -> ParkingState -> State
{- Converts a ParkingState into an index suitable for the PolicyIterate algo-
 - rithm.  The encoding depends on the number of rows, n.
 -}
encodeState _ Exit = 0
encodeState _ Crash = 1
encodeState _ (P i) = 1 + i
encodeState n (A i True) = 1 + n + i
encodeState n (A i False) = 1 + 2*n + i
encodeState n (B i True) = 1 + 3*n + i
encodeState n (B i False) = 1 + 4*n + i

decodeState :: Int -> State -> ParkingState
{- The inverse of the mapping above.
 -}
decodeState _ 0 = Exit
decodeState _ 1 = Crash
decodeState n m | m <= 1 + n = P (m-1)
                | m <= 1 + 2*n = A (m-n-1) True
                | m <= 1 + 3*n = A (m-2*n-1) False
                | m <= 1 + 4*n = B (m-3*n-1) True
                | m <= 1 + 5*n = B (m-4*n-1) False
                | otherwise = error $ "Invalid state " ++ show m ++ "."

-- Decode an association from t:
decodeT' :: Int ->
           ((State,Action,State),Double) ->
           ((ParkingState,String,ParkingState),Double)
decodeT' c ((s,a,s'),p) = ((decodeState c s,
                            if (a == 1) then "Park" else "Move",
                            decodeState c s'),
                           p)

-- Show the condensed, human-readable transition function.
decodeT :: Int ->
           [((State,Action,State),Double)] ->
           [((ParkingState,String,ParkingState),Double)]
decodeT c = map (decodeT' c)

-- Display the states in which the policy tells the agent to keep looking.
whenToPass :: Int ->
              (UArray State Double, UArray State Action, [String]) ->
              [ParkingState]
whenToPass c (_,pi,_) = map (decodeState c) .
                        filter (> 1+c) .
                        map fst .
                        filter ((==) 0 . snd) .
                        assocs $
                        pi

-- And when to park.
whenToPark :: Int ->
              (UArray State Double, UArray State Action, [String]) ->
              [ParkingState]
whenToPark c (_,pi,_) = map (decodeState c) .
                        filter (> 1+c) .
                        map fst .
                        filter ((==) 1 . snd) .
                        assocs $
                        pi

parkingProblem :: Int ->
                  Double ->
                  Double ->
                  Double ->
                  Double ->
                  Double ->
                  Double ->
                  MarkovDP
{- Generates a MDP representing a parking problem with the given parameters:
 -
 - c = The number of parking spaces per row.
 - alpha = The probability that a handicap spot will be empty.  Should be be-
 -         low, but close to, 1.0.
 - beta = The penalty for parking in a handicap spot.  Should be a large neg-
 -        ative number.
 - delta = A parameter representing the rate at which a shopper loses utility
 -         for distance.  I represent the utility of parking as inversely
 -         proportional to the square of the distance.  Should be a positive
 -         number.
 - iota = The rate of utility loss for driving around.  Should be a
 -        small positive number.
 - kappa = The penalty for crashing into another car.  Should be a larger
 -         negative number than beta.
 - lambda = A parameter determining the probability that a non-handicap spot
 -          will be empty.  I represent this as a Poisson-like distribution, as
 -          we can model the spaces as filling in from the closest on out, with
 -          cars entering as independent events and leaving after a certain
 -          period of time.  Therefore, this parameter should be a positive
 -          number, equal to the expected number of cars parked in the lot.
 -}
parkingProblem c alpha beta delta iota kappa lambda = let
  n = c*5 + 2
  r = listArray (0,n-1) ([0, kappa, beta] ++
                         [delta/(fromIntegral i) | i <- [2..c]] ++
                         (repeat (-iota)))
  t = accumArray (+) 0 ((0,0,0),(n-1,1,n-1)) (
-- Exit is a terminal state.
    [((0,0,0),1.0), ((0,1,0),1.0)] ++
-- Crashing or parking always obviates any further reward.
    [((i,0,0),1.0) | i <- [1..c+1]] ++ [((i,1,0),1.0) | i <- [1..c+1]] ++
-- Attempting to park in any occupied space always crashes the car.
    [((encodeState c (A i True),1,1),1.0) | i <- [1..c]] ++
    [((encodeState c (B i True),1,1),1.0) | i <- [1..c]] ++
-- However, parking in a free space succeeds.
    [((encodeState c (A i False),1,(encodeState c (P i))),1.0) | i <- [1..c]] ++
    [((encodeState c (B i False),1,(encodeState c (P i))),1.0) | i <- [1..c]] ++
-- An agent at A[1] can only move to B[1], which could be free or not.  The
-- probability that it is free is alpha.
    [((encodeState c (A 1 False),0,encodeState c (B 1 True)),1.0-alpha),
     ((encodeState c (A 1 True),0,encodeState c (B 1 True)),1.0-alpha),
     ((encodeState c (A 1 False),0,encodeState c (B 1 False)),alpha),
     ((encodeState c (A 1 True),0,encodeState c (B 1 False)),alpha)] ++
-- An agent at A[2] can only move to A[1], which again is free with probabil-
-- ity alpha.
    [((encodeState c (A 2 False),0,encodeState c (A 1 True)),1.0-alpha),
     ((encodeState c (A 2 True),0,encodeState c (A 1 True)),1.0-alpha),
     ((encodeState c (A 2 False),0,encodeState c (A 1 False)),alpha),
     ((encodeState c (A 2 True),0,encodeState c (A 1 False)),alpha)] ++
-- An agent at any other A[i] can only move to A[i-1].
    [((encodeState c (A i False),0,encodeState c (A (i-1) False)),quasiPoisson lambda (i-1)) | i <- [3..c]] ++
    [((encodeState c (A i True),0,encodeState c (A (i-1) False)),quasiPoisson lambda (i-1)) | i <- [3..c]] ++
    [((encodeState c (A i False),0,encodeState c (A (i-1) True)),1.0 - quasiPoisson lambda (i-1)) | i <- [3..c]] ++
    [((encodeState c (A i True),0,encodeState c (A (i-1) True)),1.0 - quasiPoisson lambda (i-1)) | i <- [3..c]] ++
-- An agent at B[i], where i < c, can only move to B[i+1].
    [((encodeState c (B i False),0,encodeState c (B (i+1) False)),quasiPoisson lambda (i+1)) | i <- [1..(c-1)]] ++
    [((encodeState c (B i True),0,encodeState c (B (i+1) False)),quasiPoisson lambda (i+1)) | i <- [1..(c-1)]] ++
    [((encodeState c (B i False),0,encodeState c (B (i+1) True)),1.0-quasiPoisson lambda (i+1)) | i <- [1..(c-1)]] ++
    [((encodeState c (B i True),0,encodeState c (B (i+1) True)),1.0-quasiPoisson lambda (i+1)) | i <- [1..(c-1)]] ++
-- Finally, an agent at B[c] can only move to A[c].
    [((encodeState c (B c False),0,encodeState c (A c True)),1.0-quasiPoisson lambda c),
     ((encodeState c (B c True),0,encodeState c (A c True)),1.0-quasiPoisson lambda c),
     ((encodeState c (B c False),0,encodeState c (A c False)),quasiPoisson lambda c),
     ((encodeState c (B c True),0,encodeState c (A c False)),quasiPoisson lambda c)]
       )
  in (n, 2, r, t)

poisson :: Double -> Int -> Double
-- The pmf of a Poisson distribution.
poisson lambda i = (exp (-lambda)) *
                   (lambda^i) /
                   (fromIntegral (factorial i))

poissonCMF :: Double -> Int -> Double
-- The memoized cmf of a Poisson distribution.
poissonCMF lambda = ((scanl1 (+) (map (poisson lambda) [0..])) !!)

quasiPoisson :: Double -> Int -> Double
{- Take into account that there are two rows to fill.  The parameter gamma
 - should be the expected number of parked cars; the parameter i should be
 - at least 2.
 -}
quasiPoisson lambda i = ((poissonCMF lambda (i*2-4)) +
                         (poissonCMF lambda (i*2-3)))
                        / 2.0

factorial :: Int -> Integer
-- A memoized factorial function.
factorial = ((scanl (*) 1 [1..]) !!)

-- Test data:

-- The 4x3 world in fig. 17.3 of (Russell & Norvig 2010):

test4x3 :: MarkovDP
test4x3 = ( 12, -- The 11 states in the example, plus a terminal state 0.
            4,  -- Left, right, up, down
            listArray (0,11)
                      [0,
                       -0.04, -0.04, -0.04, 1.0,
                       -0.04,        -0.04, -1.0,
                       -0.04, -0.04, -0.04, -0.04],
            accumArray (+) 0.0 ((0,0,0),(11,3,11)) [
 -- The synthetic terminal state:
 ((0,0,0),1.0),((0,1,0),1.0),((0,2,0),1.0),((0,3,0),1.0),
 -- Position (1,1):
 ((1,0,1),1.0),((1,1,2),1.0),((1,2,1),1.0),((1,3,5),1.0),
 -- Position (1,2):
 ((2,0,1),1.0),((2,1,3),1.0),((2,2,2),1.0),((2,3,2),1.0),
 -- Position (1,3):
 ((3,0,2),1.0),((3,1,4),1.0),((3,2,3),1.0),((3,3,6),1.0),
 -- Position (1,4):
 ((4,0,0),1.0),((4,1,0),1.0),((4,2,0),1.0),((4,3,0),1.0),
 -- Position (2,1):
 ((5,0,5),1.0),((5,1,5),1.0),((5,2,1),1.0),((5,3,8),1.0),
 -- Position (2,2) does not exist.
 -- Position (2,3):
 ((6,0,6),1.0),((6,1,7),1.0),((6,2,3),1.0),((6,3,10),1.0),
 -- Position (2,4):
 ((7,0,0),1.0),((7,1,0),1.0),((7,2,0),1.0),((7,3,0),1.0),
 -- Position (3,1):
 ((8,0,8),1.0),((8,1,9),1.0),((8,2,5),1.0),((8,3,8),1.0),
 -- Position (3,2):
 ((9,0,8),1.0),((9,1,10),1.0),((9,2,9),1.0),((9,3,9),1.0),
 -- Position (3,3):
 ((10,0,9),1.0),((10,1,11),1.0),((10,2,6),1.0),((10,3,10),1.0),
 -- Position (3,4):
 ((11,0,10),1.0),((11,1,11),1.0),((11,2,7),1.0),((11,3,11),1.0)
            ]
          )

-- The three-state world in problem 17.10 in (Russell & Norvig 2010):
test3 :: MarkovDP
test3 = ( 3, -- There are three states in the example.  I rename 3 to 0.
          2, -- and two actions, a and b.
          listArray (0,2) [0,-1,-2],
          accumArray (+) 0.0 ((0,0,0),(2,1,2)) [
 -- State 0 is a terminal state.
 ((0,0,0),1.0), ((0,1,0),1.0),
 -- In state 1, action a goes to state 2 (p=0.8) or state 1 (p=0.2)
 -- Action b goes to state 0 (p=0.1) or state 1 (p=0.9)
 ((1,0,2),0.8), ((1,0,1),0.2), ((1,1,0),0.1), ((1,1,1),0.9),
 -- In state 2, action a goes to state 1 (p=0.8) or state 2 (p=0.2)
 -- Action b goes to state 0 (p=0.1) or state 1 (p=0.9)
 ((2,0,1),0.8), ((2,0,2),0.2), ((2,1,0),0.1), ((2,1,2),0.9)
          ]
        )

-- The 101 x 3 world from exercise 17.8 in (Russell & Norvig 2010):
{- This example is useful because I’ve already solved the critical values of
 - gamma for it, and because it contains a larger number of states.
 -}
test101 :: MarkovDP
test101 = ( 204, -- The start state, an immediate reward state on top followed
                 -- by 100 penalties, an immediate penalty state on the bottom
                 -- followed by 100 rewards, and a synthetic terminal state.
            2,   -- Up or down at the start.  Both move the agent right any-
                 -- where else.
            listArray (0,203)
  ((0:50.0:(take 100 (repeat (-1.0))))++(-50.0:(take 100 (repeat 1.0)))++[0]),
            accumArray (+) 0.0 ((0,0,0),(203,1,203))
-- In the start state, up goes to state 1, down to state 102.
  ([((0,0,1),1.0), ((0,1,102),1.0)] ++
-- In states 1-100, any action moves the agent to the next consecutive state.
  [((i,0,i+1), 1.0) | i <- [1..100]] ++ [((i,1,i+1), 1.0) | i <- [1..100]] ++
-- State 101 always sends the agent to terminal state 203.
  [((101,0,203),1.0),((101,1,203),1.0)] ++
-- States 102-202 also send the agent to the next consecutive state:
  [((i,0,i+1), 1.0) | i <- [102..202]] ++
  [((i,1,i+1), 1.0) | i <- [102..202]] ++
-- Finally, state 203 is terminal:
  [((203,0,203),1.0),((203,1,203),1.0)])
          )

max_stars_repo_name: Davislor/markov

max_issues_repo_name: Davislor/markov

max_forks_repo_name: Davislor/markov


text:
{-
 - Module containing a few functions about rendering the fractals.
-}
module Rendering where

import Data.Complex
import Data.List

-- RENDERING ----
escape :: (RealFloat a, Fractional a) => [Complex a] -- list of samples
    -> (a -> a) -- Intensity mapping function
    -> a -- pixel brightness
escape samples mapFunc = sum $ zipWith (\x y -> mapFunc (magnitude (y - x)))
     samples (tail samples)

follow :: (RealFloat a, Fractional a) => [Complex a] -- list of samples
    -> (Complex a -> a) -- itensity function
    -> a
follow [] func = 0
follow samples func = func $ last samples

smooth :: (RealFloat a, Fractional a) => [Complex a]
    -> a
smooth samples = foldl' (+) 0 (map (exp . ((-1)*) . magnitude) samples)

-- to get all values from iterating Julia
sample :: (RealFloat a, Fractional a) => Int  -- nbSamples
    -> (Complex a -> Bool) -- stopping condition
    -> (Complex a -> Complex a) -- how to get next point in serie
    -> Complex a -- starting point of serie (z0)
    -> [Complex a] -- list of samples
sample nbsamples isValid funct start = ((takeWhile isValid).(take nbsamples).(iterate (funct))) start

sampleIFS :: (RealFloat a, Fractional a) => Int -- nbSamples
    -> (Complex a -> Bool) -- stop condition
    -> [(Complex a -> Complex a)] -- list of functions to apply
    -> Complex a -- starting point
    -> [Complex a]
sampleIFS nbSamples isValid functions start = (takeWhile isValid).(take nbSamples).
    (iterateList functions) $ start

iterateList :: [(a -> a)] -> a -> [a]
iterateList functions x = x : iterateList fs (f x)
    where (f:fs) = cycle functions

max_stars_repo_name: Vetii/Fractal

max_issues_repo_name: Vetii/Fractal

max_forks_repo_name: Vetii/Fractal


text:
module Math.IRT.Fisher
    ( FisherInfo(..)
    , fisherInfoObserved
    , fisherInfoExpected
    , l'
    , l''
    ) where

import Statistics.Distribution
import Math.IRT.Internal.Distribution

data FisherInfo = FisherInfo { items :: [Double]
                             , test  :: !Double
                             , sem   :: !Double
                             } deriving (Show)

-- |The observed Fisher Information
--  Applies a min/max prior to ensure a real-valued SEM
fisherInfoObserved :: (ContDistr d, DensityDeriv d) => Double -> [Bool] -> [d] -> FisherInfo
fisherInfoObserved theta resps params =
    let is = zipWith go resps params
        t  = sum is
        s  = sqrt (1 / t)
    in FisherInfo is t s
    where go u x = negate $ l'' u x theta

-- |The expected Fisher Information
fisherInfoExpected :: ContDistr d => Double -> [d] -> FisherInfo
fisherInfoExpected theta params =
    let is = map go params
        t  = sum is
        s  = sqrt (1 / t)
    in FisherInfo is t s
    where go x = let (cdf, ccdf, pdf) = pqDers x theta
                 in (square pdf) / (cdf * ccdf)

-- |The first derivative of `l`, whatever `l` is... Once again, the catIrt documentation disappoints.
-- According to catIrt, "u is the response, and x are the parameters."
-- We only implement the MLE route
l' :: ContDistr d => Bool -> d -> Double -> Double
l' u x theta =
    let (cdf, ccdf, pdf) = pqDers x theta
        denom = cdf * ccdf
        f k   = k * pdf / denom
    in case u of
         True  -> f ccdf
         False -> f cdf

-- |The second derivative of `l` (same one as in `l'`)
l'' :: (ContDistr d, DensityDeriv d) => Bool -> d -> Double -> Double
l'' u x theta =
    let (cdf, ccdf, pdf) = pqDers x theta
        pdf'             = densityDeriv x theta
    in case u of
         True  -> (((-1) / square cdf) * (square pdf))
                  + (1 / cdf * pdf')
         False -> negate $ ((1 / square ccdf) * square pdf)
                           + (1 / ccdf * pdf')

pqDers :: ContDistr d => d -> Double -> (Double, Double, Double)
pqDers x theta = let pComp = cumulative x theta
                     p'    = density x theta
                 in ( pComp
                    , (1 - pComp)
                    , p')

square :: Double -> Double
square = (^^ (2 :: Int))

max_stars_repo_name: argiopetech/irt

max_issues_repo_name: argiopetech/irt

max_forks_repo_name: argiopetech/irt


text:
-- |
-- Bidirectional version of "Data.Complex".
{-# LANGUAGE Safe, QuasiQuotes, TypeOperators #-}
module Data.Invertible.Complex
  ( complex
  , polar
  , conjugate
  ) where

import qualified Data.Complex as C

import Data.Invertible.Bijection
import Data.Invertible.Internal
import Data.Invertible.TH

-- |Convert between 'Complex' numbers and rectangular form.
complex :: (a, a) <-> C.Complex a
complex = [biCase|(r, i) <-> r C.:+ i|]

-- |Convert between complex numbers and 'C.polar' form.
polar :: RealFloat a => (a, a) <-> C.Complex a
polar = uncurry C.mkPolar :<->: C.polar

-- |The 'C.conjugate' of a complex number.
conjugate :: Num a => C.Complex a <-> C.Complex a
conjugate = involution C.conjugate

max_stars_repo_name: dylex/isomorphism

max_issues_repo_name: dylex/isomorphism

max_forks_repo_name: dylex/isomorphism


text:
module MachineLearning.NeuralNetwork.WeightInitializationTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus
import Control.Monad (replicateM)
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA
import qualified Numeric.Morpheus.Statistics as Stat
import qualified Control.Monad.Random as RndM
import MachineLearning.NeuralNetwork.WeightInitialization

assertDistribution eps x mu e = do
  let mean = Stat.mean x
  assertApproxEqual "mean" eps mu mean
  assertBool "maximum" $ (V.maximum x) <= e
  assertBool "minimum" $ (V.minimum x) >= (-e)

generateData algo n sz = do
   rndList <- replicateM n $ RndM.evalRandIO (algo sz)
   let (bs, ws) = unzip rndList
       b = LA.flatten $ LA.fromBlocks [bs]
       w = LA.flatten $ LA.fromBlocks [ws]
   return (b, w)

testWeightInitAlgo eps algo n sz mu e = do
  (b, w) <- generateData algo n sz
  assertDistribution eps b 0 0
  assertDistribution eps w mu e

heEps (r, c) = sqrt (2/(fromIntegral $ r + c))
nguyenEps (r, c) = (sqrt 6) / (sqrt . fromIntegral $ r + c)

sz = (7, 5)

tests = [ testGroup "flatten" [
            testCase "hu" $ testWeightInitAlgo 1e-1 he 100 sz 0 (heEps sz)
            , testCase "nguyen" $ testWeightInitAlgo 1e-1 nguyen 100 sz 0 (nguyenEps sz)
            ]
        ]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Vector (module Data.Complex, module Vector) where

import Data.Complex

type Vector = Complex Float
j, zero :: Vector
j = 0:+1
zero = 0:+0
normalize a = a / (sqrt $ a * conjugate a)

max_stars_repo_name: jarys/jampazen

max_issues_repo_name: jarys/jampazen

max_forks_repo_name: jarys/jampazen


text:
-- |
-- Defines Log-domain probability density functions

{-# LANGUAGE ViewPatterns, NoMonomorphismRestriction, FlexibleInstances #-}
module Math.Probably.PDF where

import qualified Math.Probably.Student as S
import Numeric.LinearAlgebra
import Control.Spoon
import Control.DeepSeq
import qualified Data.Vector.Storable as VS

-- | The type of probablility density functions
type PDF a = a->Double

instance Show (a->Double) where
    show f = error "Math.Prob.PDF: showing function"

instance Eq (a->Double) where
    f == g = error "Math.Prob.PDF: comparing functions!"

instance Num a => Num (a->Double) where
    f + g = \x-> f x + g x
    f * g = \x-> f x * g x
    f - g = \x-> f x - g x
    abs f = \x -> abs (f x)
    signum f = \x-> signum $ f x
    fromInteger = error $ "Math.Prob.PDF: frominteger on function"

-- | Uniform distributions
uniform :: (Real a, Fractional a) => a-> a-> PDF a
uniform from to = \x-> if x>=from && x <=to
                               then log $ realToFrac $ 1/(to-from )
                               else 0

--http://en.wikipedia.org/wiki/Normal_distribution
-- | Normal distribution
gauss :: (Real a, Floating a) => a-> a-> PDF a
gauss mean sd x = realToFrac $ negate $ (x-mean)**2/(2*sd*sd) + log(sd*sqrt(2*pi))

-- | Normal distribution by variance
normal :: (Real a, Floating a) => a-> a-> a -> a
normal = normalLogPdf

normalLogPdf :: (Real a, Floating a) => a-> a-> a -> a
normalLogPdf mean variance x
   = ((log 1)-(0.500*(log ((2.000*pi)*variance))))-(((x-mean)**2)/(2*variance))

-- | Normal distribution, specialised for Doubles
gaussD :: Double -> Double-> PDF Double
gaussD mean sd x = negate $ (x-mean)**2/(2*sd*sd) + log(sd*sqrt(2*pi))

gammafun = exp . S.gammaln

-- | <http://en.wikipedia.org/wiki/Gamma_distribution>
gammaD :: Double -> Double -> PDF Double
gammaD k theta x = log $ x**(k-1)*(exp(-x/theta)/(theta**k*(exp (S.gammaln k))))

invGammaD :: Double -> Double -> PDF Double
invGammaD a b x =log $ (b**a/gammafun a)*(1/x)**(a+1)*exp(-b/x)

logNormal m var x = negate $ log (x*sqrt(2*pi*var)) + (log x - m)**2 / (2*var)

logNormalD :: Double -> Double-> PDF Double
logNormalD = logNormal

-- | Beta distribution
beta a b x = log $ (recip $ S.beta a b) * x ** (a-1) * (1-x) ** (b-1)

-- | Binomial distribution. This is the Probability Mass Function, not PDF
binomial :: Int -> Double -> PDF Int
binomial n p k =
    let realk = realToFrac k
    in log (realToFrac (choose (toInteger n) (toInteger k))) + realk * log( p ) + (realToFrac $ n-k) * log( (1-p) )
 where --http://blog.plover.com/math/choose.html
   choose :: Integer -> Integer -> Integer
   choose n 0 = 1
   choose 0 k = 0
   choose (n) (k) |  k > n `div` 2 = choose n (n-k)
                  |  otherwise = (choose (n-1) (k-1)) * (n) `div` (k)

-- | the joint distribution of two independent distributions
zipPdfs :: Num a => PDF a -> PDF b -> PDF (a,b)
zipPdfs fx fy (x,y) = fx x + fy y

-- | the multiplication of two PDFS
mulPdf :: Num a => PDF a -> PDF a -> PDF a
mulPdf d1 d2 = \x -> (d1 x + d2 x)

--instance Num a => Num (PDF a) where

--instance NFData (Matrix Double)
--   where rnf mat = mapMatrix (\x-> x `seq` 1.0::Double) mat `seq` ()

-- | multivariate normal
multiNormal :: Vector Double -> Matrix Double -> PDF (Vector Double)
multiNormal mu sigma =
  let k = realToFrac $ dim mu
      invSigma =  inv sigma
      mat1 = head . head . toLists
  in \x-> log (recip ((2*pi)**(k/2) * sqrt(det sigma))) + (mat1 $ negate $ 0.5*(asRow $ x-mu) `multiply` invSigma `multiply` (asColumn $ x-mu) )

multiNormalByInv :: Double -> Matrix Double -> Vector Double -> PDF (Vector Double)
multiNormalByInv lndet invSigma mu =
  let k = realToFrac $ dim mu
      mat1 = head . head . toLists
  in \x-> log 1 - (k/2)*log (2*pi) - lndet/2 + (mat1 $ negate $ 0.5*(asRow $ x-mu) `multiply` invSigma `multiply` (asColumn $ x-mu) )

multiNormalByInvFixCov ::  Matrix Double -> Vector Double -> PDF (Vector Double)
multiNormalByInvFixCov invSigma mu =
  let k = realToFrac $ dim mu
      mat1 = head . head . toLists
  in \x-> (mat1 $ negate $ 0.5*(asRow $ x-mu) `multiply` invSigma `multiply` (asColumn $ x-mu) )

multiNormalIndep :: Vector Double -> Vector Double -> PDF (Vector Double)
multiNormalIndep vars mus xs= VS.sum $ VS.zipWith3 (\var mu x -> normalLogPdf mu var x) vars mus xs

{-mu1 = 2 |> [0, 0::Double]
sig1 = (2><2)[1::Double, 0,
              0, 1]

tst = multiNormal mu1 sig1 mu1
tsta = inv sig1
tstb = det sig1
tstc = let mu = mu1
           sigma = sig1
           x = mu1
           invSigma = inv sigma
       in (asRow $ x-mu) `multiply` invSigma `multiply` (asColumn $ x-mu)  -}

posdefify m =
   let (eigvals, eigvecM) = eigSH $ mkSym {- $ trace (show m) -}  m
       n = rows m
       eigValsVecs = map f $ zip (toList eigvals) (toColumns eigvecM)
       f (val,vec) = (abs val,vec)
       q = fromColumns $ map snd eigValsVecs
       bigLambda = diag $ fromList $ map fst eigValsVecs
   in mkSym $ q `multiply` bigLambda `multiply` inv q

mkSym m = buildMatrix (rows m) (cols m)$ \(i,j) ->if i>=j then m @@>(i,j)
                                                           else m @@>(j,i)

max_stars_repo_name: glutamate/probably

max_issues_repo_name: glutamate/probably

max_forks_repo_name: glutamate/probably


text:
module AutomatedReviewAnalyzer where

import Numeric.LinearAlgebra as LA
import Data.List.Split (splitOn)
import qualified Data.Map.Strict as M
import System.IO
import Data.Char (isAlpha, toLower)
import SentimentAnalysis

classify :: Matrix R -> Vector R -> R -> Vector R
classify x theta thata0 = cmap check (theta LA.<# LA.tr' x + theta0s)
    where
        theta0s = rows x LA.|> repeat 0
        check z = if z > 0 then 1.0 else -1.0

accuracy :: Vector R -> Vector R -> R
accuracy pred target = sum xs / fromIntegral (length xs)
    where
        xs = zipWith check (toList pred) (toList target)
        check a b = if a == b then 1 else 0

classifierAccuracy :: (Matrix R -> Vector R -> Int -> [Int] -> (Vector R, R))
                            -> Matrix R -> Matrix R -> Vector R -> Vector R -> Int -> [Int] -> (R, R)
classifierAccuracy classifier x x' y y' t zs =
    (accuracy (classify x theta theta0) y, accuracy (classify x' theta theta0) y')
    where
        (theta, theta0) = classifier x y t zs

classifierAccuracy' :: (Matrix R -> Vector R -> Int -> R -> [Int] -> (Vector R, R))
                            -> Matrix R -> Matrix R -> Vector R -> Vector R -> Int -> R -> [Int] -> (R, R)
classifierAccuracy' classifier x x' y y' t lambda zs =
    (accuracy (classify x theta theta0) y, accuracy (classify x' theta theta0) y')
    where
        (theta, theta0) = classifier x y t lambda zs

separate :: String -> String
separate = concatMap separate'
    where
        separate' c = if isAlpha c then [toLower c] else [' ', c, ' ']

addWord :: String -> M.Map String R -> M.Map String R
addWord w map = if M.member w map then map else M.insert w 1.0 map

countWord :: String -> M.Map String R -> M.Map String R
countWord w = M.insertWith (+) w 1

extractFeatures :: [String] -> M.Map String R -> Matrix R
extractFeatures xs d = (length xss >< M.size d) ys
    where
        xss = map (words . separate) xs
        ds = map (foldl (flip countWord) M.empty) xss
        ys = concatMap f $ zip ds (repeat d)
        f (d'', d') = map (g . flip M.lookup d'') (M.keys d')
        g Nothing = 0
        g (Just a) = a

main :: IO ()
main = do

    h1 <- openFile "reviews_train.tsv" ReadMode
    hSetEncoding h1 latin1
    trainDta <- hGetContents h1
    h2 <- openFile "reviews_val.tsv" ReadMode
    hSetEncoding h2 latin1
    valDta <- hGetContents h2
--    testDta <- readFile "reviews_test.tsv"
    orderDta <- readFile "4000.txt"

    let t = 20
    let lambda = 0.005
    let (labels, text) = unzip $ map ((\ i -> ((read . head) i, i !! 4)) . splitOn "\t") (drop 1 $ lines trainDta)
    let bagOfWords =  foldr addWord M.empty $ concatMap (words . separate) text
    let features = extractFeatures text bagOfWords
    let (labels', text') = unzip $ map ((\ i -> ((read . head) i, i !! 4)) . splitOn "\t") (drop 1 $ lines valDta)
    let features' = extractFeatures text' bagOfWords

    let indices = map read $ splitOn "," orderDta

    let (trainAccuracy, validationAccuracy) = classifierAccuracy perceptron features features' (vector labels) (vector labels') t indices
    let (trainAccuracy', validationAccuracy') = classifierAccuracy averagePerceptron features features' (vector labels) (vector labels') t indices
    let (trainAccuracy'', validationAccuracy'') = classifierAccuracy' pegasos features features' (vector labels) (vector labels') t lambda indices

    putStrLn "Classifier accuracy:"
--    print $ (takeColumns 10 . takeRows 4) features
    print $ "Training accuracy for perceptron: " ++ show trainAccuracy
    print $ "Validation accuracy for perceptron: " ++ show validationAccuracy
    print $ "Training accuracy for average perceptron: " ++ show trainAccuracy'
    print $ "Validation accuracy for average perceptron: " ++ show validationAccuracy'
    print $ "Training accuracy for pegasos: " ++ show trainAccuracy''
    print $ "Validation accuracy for pegasos: " ++ show validationAccuracy''


text:
module FxForHasql.Prelude
(
  module Exports,
)
where

-- base
-------------------------
import Control.Applicative as Exports
import Control.Arrow as Exports hiding (first, second)
import Control.Category as Exports
import Control.Concurrent as Exports
import Control.Exception as Exports
import Control.Monad as Exports hiding (fail, mapM_, sequence_, forM_, msum, mapM, sequence, forM)
import Control.Monad.IO.Class as Exports
import Control.Monad.Fail as Exports
import Control.Monad.Fix as Exports hiding (fix)
import Control.Monad.ST as Exports
import Data.Bifunctor as Exports
import Data.Bits as Exports
import Data.Bool as Exports
import Data.Char as Exports
import Data.Coerce as Exports
import Data.Complex as Exports
import Data.Data as Exports
import Data.Dynamic as Exports
import Data.Either as Exports
import Data.Fixed as Exports
import Data.Foldable as Exports
import Data.Function as Exports hiding (id, (.))
import Data.Functor as Exports
import Data.Functor.Compose as Exports
import Data.Functor.Contravariant as Exports
import Data.Int as Exports
import Data.IORef as Exports
import Data.Ix as Exports
import Data.List as Exports hiding (sortOn, isSubsequenceOf, uncons, concat, foldr, foldl1, maximum, minimum, product, sum, all, and, any, concatMap, elem, foldl, foldr1, notElem, or, find, maximumBy, minimumBy, mapAccumL, mapAccumR, foldl')
import Data.List.NonEmpty as Exports (NonEmpty(..))
import Data.Maybe as Exports
import Data.Monoid as Exports hiding (Alt)
import Data.Ord as Exports
import Data.Proxy as Exports
import Data.Ratio as Exports
import Data.STRef as Exports
import Data.String as Exports
import Data.Traversable as Exports
import Data.Tuple as Exports
import Data.Unique as Exports
import Data.Version as Exports
import Data.Void as Exports
import Data.Word as Exports
import Debug.Trace as Exports
import Foreign.ForeignPtr as Exports
import Foreign.Ptr as Exports
import Foreign.StablePtr as Exports
import Foreign.Storable as Exports
import GHC.Conc as Exports hiding (orElse, withMVar, threadWaitWriteSTM, threadWaitWrite, threadWaitReadSTM, threadWaitRead)
import GHC.Exts as Exports (lazy, inline, sortWith, groupWith)
import GHC.Generics as Exports (Generic)
import GHC.IO.Exception as Exports
import Numeric as Exports
import Prelude as Exports hiding (fail, concat, foldr, mapM_, sequence_, foldl1, maximum, minimum, product, sum, all, and, any, concatMap, elem, foldl, foldr1, notElem, or, mapM, sequence, id, (.))
import System.Environment as Exports
import System.Exit as Exports
import System.IO as Exports (Handle, hClose)
import System.IO.Error as Exports
import System.IO.Unsafe as Exports
import System.Mem as Exports
import System.Mem.StableName as Exports
import System.Timeout as Exports
import Text.ParserCombinators.ReadP as Exports (ReadP, ReadS, readP_to_S, readS_to_P)
import Text.ParserCombinators.ReadPrec as Exports (ReadPrec, readPrec_to_P, readP_to_Prec, readPrec_to_S, readS_to_Prec)
import Text.Printf as Exports (printf, hPrintf)
import Text.Read as Exports (Read(..), readMaybe, readEither)
import Unsafe.Coerce as Exports

-- text
-------------------------
import Data.Text as Exports (Text)

-- bytestring
-------------------------
import Data.ByteString as Exports (ByteString)

-- fx
-------------------------
import Fx as Exports

max_stars_repo_name: nikita-volkov/fx-for-hasql

max_issues_repo_name: nikita-volkov/fx-for-hasql

max_forks_repo_name: nikita-volkov/fx-for-hasql


text:
{-# LANGUAGE FlexibleInstances #-}
-----------------------------------------------------------------------------
--
-- Module      :  Math.Test
-- Copyright   :
-- License     :  AllRightsReserved
--
-- Maintainer  :
-- Stability   :
-- Portability :
--
-- |
--
-----------------------------------------------------------------------------

module Math.Test where

import Test.Framework (testGroup, Test)
import Test.Framework.Providers.HUnit
import Test.Framework.Providers.QuickCheck2 (testProperty)
import Test.HUnit hiding (Test)
import Test.QuickCheck
import Data.AEq (AEq,(~==))

import Numeric.LinearAlgebra (Matrix, Vector, Element, fromList, toRows, toList, flatten)
import Foreign.Storable (Storable)

import Math

mathSuite :: Test
mathSuite = testGroup "Auxillary math functions"
    [-- testCase "The mean of a vector is the vector itself" (testMean [(fromList [1.0]), (fromList [1.0])] (fromList [1.0]) ),
     -- testCase "should fail" (testMean [(fromList [1.0]), (fromList [2.0])] (fromList [1.0]) ),
     testProperty "The mean of a vector is the vector itself" meanProp1,
     testProperty "Mean addition" mean_addition,
     testProperty "Alternate ScatterMatrix" scatter_alternative,
     testProperty "Mean subtraction" mean_subtraction,
     testProperty "ScatterMatrix addition" scatter_addition,
     testProperty "ScatterMatrix subtraction" scatter_subtraction,
     testProperty "The addition of Summarystats" sstat_addition,
     testProperty "The subtraction of Summarystats" sstat_subtraction
     ]

testMean :: [Vector Double] -> Vector Double -> Assertion
testMean points expected = expected @=? meanv points

instance Arbitrary (Vector Double) where
    arbitrary = sized $ \n -> resize (n+1) $ fmap fromList arbitrary

instance Arbitrary X where
    arbitrary = sized $ \n -> sequence [fmap fromList (vector 3) | i<-[1..n]]

instance (AEq a, Storable a)=> AEq (Vector a) where
    x ~== y = (toList x) ~== (toList y)

instance AEq (Matrix Double) where
    x ~== y = (flatten x) ~== (flatten y)

scatter_alternative x = length x > 0 ==> scatterMatrix x ~== scatterMatrixAlt x
            where types = (x::X)

meanProp1 x = meanv [x] == x
        where types = (x::Vector Double)

mean_addition x k = k <= length x ==> add_means nx (meanv xs) ny (meanv ys) ~== meanv (xs++ys)
            where nx = length xs
                  ny = length ys
                  xs = take k x
                  ys = drop k x
                  types = (x::X, k::Int)

mean_subtraction x k = k <= length x ==> sub_means (nx+ny) (meanv (xs++ys)) ny (meanv ys) ~== meanv xs
            where nx = length xs
                  ny = length ys
                  xs = take k x
                  ys = drop k x
                  types = (x::X, k::Int)

scatter_addition x k = k <= length x ==> add_scatter nx sx mx ny sy my ~== scatterMatrix x
            where nx = length xs
                  ny = length ys
                  xs = take k x
                  ys = drop k x
                  mx = meanv xs
                  my = meanv ys
                  sx = scatterMatrix xs
                  sy = scatterMatrix ys
                  types = (x::X, k::Int)

scatter_subtraction x k = k>= 1 && k >= (quot (length x) 2)+1 && k <= length x ==> sub_scatter n s m ny sy my ~== scatterMatrix xs
            where xs = take k x
                  ys = drop k x
                  n = length x
                  ny = length ys
                  m = meanv x
                  my = meanv ys
                  s = scatterMatrix x
                  sy = scatterMatrix ys
                  types = (x::X, k::Int)

sstat_addition x k = k <= length x ==> (fromX xs) <+> (fromX ys) ~== (fromX x)
            where   xs = take k x
                    ys = drop k x
                    types = (x::X, k::Int)

sstat_subtraction x k = k < length x ==> (fromX x) <-> (fromX xs) ~== (fromX ys)
            where   xs = take k x
                    ys = drop k x
                    types = (x::X, k::Int)

validX xs ys = length xs > 1
               && length ys > 1
               && dimension xs > 1
               && dimension ys > 1
               && dimension xs == dimension ys

max_stars_repo_name: juliankopkalarsen/FpStats

max_issues_repo_name: juliankopkalarsen/FpStats

max_forks_repo_name: juliankopkalarsen/FpStats


text:
-- |
-- Module      : Data.Express.Name
-- Copyright   : (c) 2019-2021 Rudy Matela
-- License     : 3-Clause BSD  (see the file LICENSE)
-- Maintainer  : Rudy Matela <rudy@matela.com.br>
--
-- Defines the 'Name' type class.
module Data.Express.Name
  ( Name (..)
  , names
  , variableNamesFromTemplate
  )
where

import Data.Express.Utils.String

import Data.List

-- for instances:
import Data.Int
import Data.Word
import Data.Ratio
import Data.Complex
import Data.Char

-- |
-- If we were to come up with a variable name for the given type
-- what 'name' would it be?
--
-- An instance for a given type @ Ty @ is simply given by:
--
-- > instance Name Ty where name _ = "x"
--
-- Examples:
--
-- > > name (undefined :: Int)
-- > "x"
--
-- > > name (undefined :: Bool)
-- > "p"
--
-- > > name (undefined :: [Int])
-- > "xs"
--
-- This is then used to generate an infinite list of variable 'names':
--
-- > > names (undefined :: Int)
-- > ["x", "y", "z", "x'", "y'", "z'", "x''", "y''", "z''", ...]
--
-- > > names (undefined :: Bool)
-- > ["p", "q", "r", "p'", "q'", "r'", "p''", "q''", "r''", ...]
--
-- > > names (undefined :: [Int])
-- > ["xs", "ys", "zs", "xs'", "ys'", "zs'", "xs''", "ys''", ...]
class Name a where
  -- | /O(1)./
  --
  -- Returns a name for a variable of the given argument's type.
  --
  -- > > name (undefined :: Int)
  -- > "x"
  --
  -- > > name (undefined :: [Bool])
  -- > "ps"
  --
  -- > > name (undefined :: [Maybe Integer])
  -- > "mxs"
  --
  -- The default definition is:
  --
  -- > name _ = "x"
  name :: a -> String
  name _ = "x"

-- |
-- > name (undefined :: ()) = "u"
-- > names (undefined :: ()) = ["u", "v", "w", "u'", "v'", ...]
instance Name ()        where  name _  =  "u"

-- |
-- > name (undefined :: Bool) = "p"
-- > names (undefined :: Bool) = ["p", "q", "r", "p'", "q'", ...]
instance Name Bool      where  name _  =  "p"

-- |
-- > name (undefined :: Int) = "x"
-- > names (undefined :: Int) = ["x", "y", "z", "x'", "y'", ...]
instance Name Int       where  name _  =  "x"

-- |
-- > name (undefined :: Integer) = "x"
-- > names (undefined :: Integer) = ["x", "y", "z", "x'", ...]
instance Name Integer   where  name _  =  "x"

-- |
-- > name (undefined :: Char) = "c"
-- > names (undefined :: Char) = ["c", "d", "e", "c'", "d'", ...]
instance Name Char      where  name _  =  "c"

-- |
-- > name (undefined :: Ordering) = "o"
-- > names (undefined :: Ordering) = ["o", "p", "q", "o'", ...]
instance Name Ordering  where  name _  =  "o"

-- |
-- > name (undefined :: Rational) = "q"
-- > names (undefined :: Rational) = ["q", "r", "s", "q'", ...]
instance Name (Ratio a) where  name _  =  "q"

-- |
-- > name (undefined :: Complex) = "x"
-- > names (undefined :: Complex) = ["x", "y", "z", "x'", ...]
instance Name (Complex a)  where  name _  =  "x"

-- |
-- > name (undefined :: Float) = "x"
-- > names (undefined :: Float) = ["x", "y", "z", "x'", ...]
instance Name Float     where  name _  =  "x"

-- |
-- > name (undefined :: Double) = "x"
-- > names (undefined :: Double) = ["x", "y", "z", "x'", ...]
instance Name Double    where  name _  =  "x"

-- |
-- > names (undefined :: ()->()) = ["f", "g", "h", "f'", ...]
-- > names (undefined :: Int->Int) = ["f", "g", "h", ...]
instance Name (a -> b)  where  name _  =  "f"

-- |
-- > names (undefined :: Maybe Int) = ["mx", "mx1", "mx2", ...]
-- > nemes (undefined :: Maybe Bool) = ["mp", "mp1", "mp2", ...]
instance Name a => Name (Maybe a) where
  name mx  =  "m" ++ name x
    where
    Just x = mx

-- |
-- > names (undefined :: Either Int Int) = ["exy", "exy1", ...]
-- > names (undefined :: Either Int Bool) = ["exp", "exp1", ...]
instance (Name a, Name b) => Name (Either a b) where
  name exy  =  "e" ++ n ++ m
    where
    Left x  = exy
    Right y = exy
    n = name x
    m = head $ names y \\ [n]

-- |
-- > names (undefined :: (Int,Int)) = ["xy", "zw", "xy'", ...]
-- > names (undefined :: (Bool,Bool)) = ["pq", "rs", "pq'", ...]
instance (Name a, Name b) => Name (a,b) where
  name xy  =  n ++ m
    where
    (x,y)  =  xy
    n  =  name x
    m  =  head $ names y \\ [n]

-- |
-- > names (undefined :: (Int,Int,Int)) = ["xyz","uvw", ...]
-- > names (undefined :: (Int,Bool,Char)) = ["xpc", "xpc1", ...]
instance (Name a, Name b, Name c) => Name (a,b,c) where
  name xyz  =  n ++ m ++ o
    where
    (x,y,z)  =  xyz
    n  =  name x
    m  =  head $ names y \\ [n]
    o  =  head $ names z \\ [n,m]

-- |
-- > names (undefined :: ((),(),(),())) = ["uuuu", "uuuu1", ...]
-- > names (undefined :: (Int,Int,Int,Int)) = ["xxxx", ...]
instance (Name a, Name b, Name c, Name d) => Name (a,b,c,d) where
  name xyzw  =  name x ++ name y ++ name z ++ name w  where  (x,y,z,w)  =  xyzw

-- |
-- > names (undefined :: [Int]) = ["xs", "ys", "zs", "xs'", ...]
-- > names (undefined :: [Bool]) = ["ps", "qs", "rs", "ps'", ...]
instance Name a => Name [a] where
  name xs  =  name (head xs) ++ "s"

-- |
-- Returns na infinite list of variable names from the given type:
-- the result of 'variableNamesFromTemplate' after 'name'.
--
-- > > names (undefined :: Int)
-- > ["x", "y", "z", "x'", "y'", "z'", "x''", "y''", "z''", ...]
--
-- > > names (undefined :: Bool)
-- > ["p", "q", "r", "p'", "q'", "r'", "p''", "q''", "r''", ...]
--
-- > > names (undefined :: [Int])
-- > ["xs", "ys", "zs", "xs'", "ys'", "zs'", "xs''", "ys''", ...]
names :: Name a => a -> [String]
names  =  variableNamesFromTemplate . name

-- instances of further types and arities --

instance Name Word where  name _  =  "x"

instance Name Int8   where  name _  =  "x"
instance Name Int16  where  name _  =  "x"
instance Name Int32  where  name _  =  "x"
instance Name Int64  where  name _  =  "x"
instance Name Word8  where  name _  =  "x"
instance Name Word16 where  name _  =  "x"
instance Name Word32 where  name _  =  "x"
instance Name Word64 where  name _  =  "x"
instance Name GeneralCategory where  name _  =  "c"

instance (Name a, Name b, Name c, Name d, Name e) => Name (a,b,c,d,e) where
  name xyzwv  =  name x ++ name y ++ name z ++ name w ++ name v
    where  (x,y,z,w,v)  =  xyzwv

instance (Name a, Name b, Name c, Name d, Name e, Name f)
      => Name (a,b,c,d,e,f) where
  name xyzwvu  =  name x ++ name y ++ name z ++ name w ++ name v ++ name u
    where  (x,y,z,w,v,u)  =  xyzwvu

instance (Name a, Name b, Name c, Name d, Name e, Name f, Name g)
      => Name (a,b,c,d,e,f,g) where
  name xyzwvut  =  name x ++ name y ++ name z ++ name w
                ++ name v ++ name u ++ name t
    where  (x,y,z,w,v,u,t)  =  xyzwvut

instance (Name a, Name b, Name c, Name d, Name e, Name f, Name g, Name h)
      => Name (a,b,c,d,e,f,g,h) where
  name xyzwvuts  =  name x ++ name y ++ name z ++ name w
                 ++ name v ++ name u ++ name t ++ name s
    where  (x,y,z,w,v,u,t,s)  =  xyzwvuts

instance ( Name a, Name b, Name c, Name d
         , Name e, Name f, Name g, Name h
         , Name i)
      => Name (a,b,c,d,e,f,g,h,i) where
  name xyzwvutsr  =  name x ++ name y ++ name z ++ name w
                  ++ name v ++ name u ++ name t ++ name s
                  ++ name r
    where  (x,y,z,w,v,u,t,s,r)  =  xyzwvutsr

instance ( Name a, Name b, Name c, Name d
         , Name e, Name f, Name g, Name h
         , Name i, Name j )
      => Name (a,b,c,d,e,f,g,h,i,j) where
  name xyzwvutsrq  =  name x ++ name y ++ name z ++ name w
                   ++ name v ++ name u ++ name t ++ name s
                   ++ name r ++ name q
    where  (x,y,z,w,v,u,t,s,r,q)  =  xyzwvutsrq

instance ( Name a, Name b, Name c, Name d
         , Name e, Name f, Name g, Name h
         , Name i, Name j, Name k )
      => Name (a,b,c,d,e,f,g,h,i,j,k) where
  name xyzwvutsrqp  =  name x ++ name y ++ name z ++ name w
                    ++ name v ++ name u ++ name t ++ name s
                    ++ name r ++ name q ++ name p
    where  (x,y,z,w,v,u,t,s,r,q,p)  =  xyzwvutsrqp

instance ( Name a, Name b, Name c, Name d
         , Name e, Name f, Name g, Name h
         , Name i, Name j, Name k, Name l )
      => Name (a,b,c,d,e,f,g,h,i,j,k,l) where
  name xyzwvutsrqpo  =  name x ++ name y ++ name z ++ name w
                     ++ name v ++ name u ++ name t ++ name s
                     ++ name r ++ name q ++ name p ++ name o
    where  (x,y,z,w,v,u,t,s,r,q,p,o)  =  xyzwvutsrqpo

max_stars_repo_name: rudymatela/express

max_issues_repo_name: rudymatela/express

max_forks_repo_name: rudymatela/express


text:
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE RankNTypes #-}
{-# OPTIONS_GHC -Wincomplete-patterns #-}
{-# OPTIONS_GHC -Wno-incomplete-patterns #-}

module Num where

import AST
import Control.Monad.Except
import qualified Data.Complex as C

-- --------------------------------------------------------------------------
-- Pogger number
-- This module implements some conversion between pogger numbers.
-- ophan instance. Why not? Nobody will define another Num for PoggerNum
-- Somewhere else.
-- --------------------------------------------------------------------------

-- | convert a PoggerNum to Complex
poggerNumToComplex :: PoggerNum -> C.Complex Double
poggerNumToComplex (Integer n) = (realToFrac n) C.:+ 0.0
poggerNumToComplex (Real r) = r C.:+ 0.0
poggerNumToComplex (Rational a b) = (realToFrac a / realToFrac b) C.:+ 0.0
poggerNumToComplex (Complex a b) = a C.:+ b

poggerNumToDouble :: PoggerNum -> Double
poggerNumToDouble (Integer n) = fromIntegral n
poggerNumToDouble (Real r) = r
poggerNumToDouble (Rational a b) = fromIntegral a / fromIntegral b
poggerNumToDouble (Complex _ _) = error "Can't downgrade complex to real"

-- | Convert complex to pogger complex.
toPoggerComplex :: C.Complex Double -> PoggerNum
toPoggerComplex c = Complex (C.realPart c) (C.imagPart c)

-- | Generic binary operator that handles coersion rules.
poggerBinop :: (forall a. Num a => a -> a -> a) -> PoggerNum -> PoggerNum -> PoggerNum
poggerBinop op = mkop
  where
   mkop (Integer m) (Integer n) = Integer (m `op` n)
   mkop (Rational m n) (Rational m' n') = Rational (m `op` m') (n `op` n')
   mkop (Real r) (Real r') = Real (r `op` r')
   mkop (Complex m n) (Complex m' n') = Complex (m `op` m') (n `op` n')
   mkop (Rational a b) (Integer n) = simplify $ Rational (a `op` (n * b)) b
   mkop i@(Integer _) r@(Rational _ _) = r `op` i
   mkop m@(Complex _ _) n = toPoggerComplex (poggerNumToComplex m `op` poggerNumToComplex n)
   mkop m n@(Complex _ _) = m `op` n
   mkop m@(Real _) n = (Real . C.realPart) (poggerNumToComplex m `op` poggerNumToComplex n)
   mkop m n@(Real _) = m `op` n

{-# INLINE poggerBinop #-}

-- | absolute value of pogger num
-- >>> abs (Integer 2)
-- Integer 2
poggerUnary :: (forall a. Num a => a -> a) -> PoggerNum -> PoggerNum
poggerUnary op = \case
  Integer m -> Integer (op m)
  Rational m n -> Rational (op m) (op n)
  Real m -> Real (op m)
  c@(Complex _ _) ->
    let a = op (poggerNumToComplex c)
     in Complex (C.realPart a) (C.imagPart a)

instance Num PoggerNum where
  (+) = poggerBinop (+)
  {-# INLINE (+) #-}

  (*) = poggerBinop (*)
  {-# INLINE (*) #-}

  abs = poggerUnary abs
  {-# INLINE abs #-}

  negate = poggerUnary negate
  {-# INLINE negate #-}

  fromInteger = Integer
  {-# INLINE fromInteger #-}

  signum = poggerUnary signum
  {-# INLINE signum #-}

-- | division on pogger number
poggerNumDivide :: PoggerNum -> PoggerNum -> PoggerNum
poggerNumDivide (Real r) (Real r') = Real (r / r')
poggerNumDivide c@(Complex _ _) c'@(Complex _ _) =
  toPoggerComplex (poggerNumToComplex c / poggerNumToComplex c')
poggerNumDivide (Integer m) (Integer n) = (Integer . floor) ((fromIntegral m) / (fromIntegral n))
poggerNumDivide (Rational m n) (Rational m' n') = simplify (Rational (m * n') (n * m'))

-- >>> (abs . negate) (Rational 2 3 / Rational 5 3)
-- Rational 2 5
-- >>> Rational 2 3 / Rational 2 6
-- Integer 2
instance Fractional PoggerNum where
  (/) = poggerNumDivide
  {-# INLINE (/) #-}

  fromRational = Real . fromRational
  {-# INLINE fromRational #-}

-- | simply rational terms.
simplify :: PoggerNum -> PoggerNum
simplify (Rational a b) =
  let v = gcd a b
      r@(Rational m n) = Rational (a `div` v) (b `div` v)
   in if n == 1
        then Integer m
        else Rational m n
simplify _ = error "trying to simply non-rational"
{-# INLINE simplify #-}

-- | ord for pogger number
instance Ord PoggerNum where
  m <= n = poggerNumToDouble m <= poggerNumToDouble n
  {-# INLINE (<=) #-}

instance Real PoggerNum where
  toRational (Integer n) = fromIntegral n
  toRational (Real r) = toRational r
  toRational (Rational a b) = fromIntegral a / fromIntegral b
  toRational (Complex a _) = toRational a
  {-# INLINE toRational #-}

instance Enum PoggerNum where
  toEnum n = Integer (fromIntegral n)
  {-# INLINE toEnum #-}

  fromEnum (Integer n) = fromIntegral n
  fromEnum _ = (-1)
  {-# INLINE fromEnum #-}

-- | pogger
instance Integral PoggerNum where
  toInteger (Integer n) = n
  toInteger (Real r) = floor r
  toInteger (Rational a b) = floor (fromIntegral a / fromIntegral b)
  toInteger (Complex a _) = floor a
  {-# INLINE toInteger #-}

  quotRem (Integer a) (Integer b) =
    let (a', b') = (quotRem a b)
     in (Integer a', Integer b')
  quotRem a b = error ("no division algorithm defined for" ++ show a ++ show b)

-- | Safe wrappers
safeQuotRem :: PoggerNum -> PoggerNum -> ThrowsError (PoggerNum, PoggerNum)
safeQuotRem a@(Integer _) b@(Integer _) = return (quotRem a b)
safeQuotRem a _ = throwError (TypeMisMatch "can't mod on non integer values" (Number a))

safeMod :: PoggerNum -> PoggerNum -> ThrowsError PoggerNum
safeMod a b = safeQuotRem a b >>= return . snd

safeDiv :: PoggerNum -> PoggerNum -> ThrowsError PoggerNum
safeDiv a b = safeQuotRem a b >>= return . fst

max_stars_repo_name: ailrk/poggerscheme

max_issues_repo_name: ailrk/poggerscheme

max_forks_repo_name: ailrk/poggerscheme


text:
{-# LANGUAGE FlexibleContexts #-}

module Numeric.Trainee.Data.Csv (
	parseCsv, parseCsvFile
	) where

import Prelude.Unicode

import Control.DeepSeq
import Control.Exception (throwIO)
import qualified Data.Vector as V
import qualified Data.Text as T
import qualified Data.Text.IO as T
import Numeric.LinearAlgebra hiding (toList)

import Numeric.Trainee.Types
import Numeric.Trainee.Data.Common

parseCsv ∷ Bool → T.Text → CsvM (Sample (Vector a) (Vector a)) → Either String (Samples (Vector a) (Vector a))
parseCsv header' dat p = mapM (runCsvM p ∘ V.fromList ∘ T.split (≡ ',')) $ V.fromList $ dropHeader $ T.lines dat where
	dropHeader
		| header' = tail
		| otherwise = id

parseCsvFile ∷ Bool → FilePath → CsvM (Sample (Vector a) (Vector a)) → IO (Samples (Vector a) (Vector a))
parseCsvFile header' fpath p = do
	cts ← T.readFile fpath
	case parseCsv header' cts p of
		Left e → throwIO $ userError $ "Error parsing " ++ fpath ++ ": " ++ e
		Right rs → return $!! rs

max_stars_repo_name: mvoidex/trainee

max_issues_repo_name: mvoidex/trainee

max_forks_repo_name: mvoidex/trainee


text:
{-# LANGUAGE FlexibleInstances #-}

module HVX.SymbolicSubgradSpec (main, spec) where

import Test.Hspec
import Test.QuickCheck hiding ((><))
import Numeric.LinearAlgebra

import HVX.Primitives
import HVX.Internal.Matrix
import HVX.Internal.Primitives
import HVX.Internal.TestUtil
import HVX.Internal.SymbolicSubgrad

main :: IO ()
main = hspec spec

spec :: Spec
spec = describe "SymbolicSubgrad" $ do
  evaluateSpec
  jacobianWrtVarSpec

evaluateSpec :: Spec
evaluateSpec = describe "evaluate" $ do

  describe "log/exp" $
    it "log(exp(x)) = x" $ property $
      forAll smallVectors $
        \x -> fpequalsMat
          (evaluate (EVar "x") [("x", x)])
          (evaluate (EFun Log (EFun Exp (EVar "x"))) [("x", x)])

  describe "log_sum_exp" $
    it "log_sum_exp(x) = log(ones'*exp(x))" $ property $
      forAll smallVectors $
        \x -> let n = rows x in
          fpequalsMat
            (evaluate (EFun LogSumExp (EVar "x")) [("x", x)])
            (evaluate
              (EFun Log
                (EFun (Mul (konst 1.0 (1, n)))
                  (EFun Exp (EVar "x"))))
              [("x", x)])

  describe "max/min" $
    it "max(x) = -min(-x)" $ property $
      forAll vectors $
        \x -> fpequalsMat
          (evaluate (EFun Max (EVar "x")) [("x", x)])
          (evaluate (EFun Neg (EFun Min (EFun Neg (EVar "x")))) [("x", x)])

  describe "lpnorm" $ do
    it "lpnorm(x, 1) = sum(abs(x))" $ property $
      forAll vectors $
        \x -> fpequalsMat
          (evaluate
            (norm 1 (EVar "x"))
            [("x", x)])
          (evaluate
            (EFun (Mul (konst 1.0 (1, (rows x)))) (habs (EVar "x")))
            [("x", x)])
    it "lpnorm(x, 2)^2 = x^T * x" $ property $
      forAll vectors $
        \x -> fpequalsMat
          (tr x <> x)
          (evaluate
            (EFun (PowBaseP1InfEven 2) (norm 2 (EVar "x")))
            [("x", x)])

  describe "quadform" $
    it "quadform(I, x) = x^T x" $ property $
      forAll vectors $
        \x -> fpequalsMat
          (evaluate
            (EFun (Quadform (ident $ rows x)) (EVar "x"))
            [("x", x)])
          (tr x <> x)

  describe "pow (base variable)" $ do
    it "(x^a)^(1/a) = x (a non integral, a > 1)" $ property $
      forAll (pairsOf positiveVectors smallNonIntPowersGreaterThanOne) $
        \(x, p) -> fpequalsMat
          (evaluate (EVar "x") [("x", x)])
          (evaluate
            (EFun (PowBaseP1InfNotInt p) (EFun (PowBaseP01 (1.0/p)) (EVar "x")))
            [("x", x)])
    it "(x^a)^(1/a) = x (a even, a > 1)" $ property $
      forAll (pairsOf positiveVectors smallEvenPowersGreaterThanOne) $
        \(x, p) -> fpequalsMat
          (evaluate (EVar "x") [("x", x)])
          (evaluate
            (EFun (PowBaseP1InfEven (round p :: Integer)) (EFun (PowBaseP01 (1.0/p)) (EVar "x")))
            [("x", x)])

jacobianWrtVarSpec :: Spec
jacobianWrtVarSpec = describe "jacobianWrtVar" $ do

  it "works with respect to a variable that doesn't occur in Expr." $
    forAll (pairsOf vectors vectors) $
      \(x, y) -> fpequalsMat
        (jacobianWrtVar
          (EVar "x")
          [("x", x), ("y", y)]
          "y")
        (konst 0.0 ((rows x), (rows y)))

  describe "log/exp" $
    it "log(exp(x)) - x has zero jacobian" $ property $
      forAll smallVectors $
        \x -> fpequalsMat
          (jacobianWrtVar
            (EAdd (EFun Log (EFun Exp (EVar "x"))) (EFun Neg (EVar "x")))
            [("x", x)]
            "x")
          (zeroMat $ rows x)

  describe "log_sum_exp" $
    it "jacobian log_sum_exp(x) = jacobian log(ones'*exp(x))" $ property $
      forAll smallVectors $
        \x -> fpequalsMat
          (jacobianWrtVar (EFun LogSumExp (EVar "x")) [("x", x)] "x")
          (jacobianWrtVar
            (EFun Log
              (EFun (Mul (konst 1.0 (1, (rows x))))
                (EFun Exp (EVar "x"))))
            [("x", x)]
            "x")

  describe "max/min" $
    it "max(x) + min(-x) has zero jacobian" $ property $
      forAll vectors $
        \x -> fpequalsMat
          (jacobianWrtVar
            (EAdd (EFun Max (EVar "x")) (EFun Min (EFun Neg (EVar "x"))))
            [("x", x)]
            "x")
          (tr $ zeroVec $ rows x)

  describe "lpnorm" $ do
    it "jacobian lpnorm(x, 1) = jacobian sum(abs(x))" $ property $
      forAll vectors $
        \x -> fpequalsMat
          (jacobianWrtVar
            (norm 1 (EVar "x"))
            [("x", x)]
            "x")
          (jacobianWrtVar
            (EFun (Mul (konst 1.0 (1, (rows x)))) (habs (EVar "x")))
            [("x", x)]
            "x")
    it "jacobian lpnorm(0, 2) = 0" $ property $
      \(Positive n) -> fpequalsMat
        (jacobianWrtVar (norm 2 (EVar "x")) [("x", zeroVec n)] "x")
        (tr $ zeroVec n)

  describe "quadform" $
    it "jacobian quadform(I, x) = jacobian lpnorm(x, 2)^2" $ property $
      forAll vectors $
        \x -> fpequalsMat
          (jacobianWrtVar (EFun (Quadform (ident $ rows x)) (EVar "x")) [("x", x)] "x")
          (jacobianWrtVar (EFun (PowBaseP1InfEven 2) (norm 2 (EVar "x"))) [("x", x)] "x")

  describe "pow (base variable)" $ do
    it "(x^a)^(1/a) - x has zero jacobian (a non integral, a > 1)" $ property $
      forAll (pairsOf positiveVectors smallNonIntPowersGreaterThanOne) $
        \(x, p) -> fpequalsMat
          (jacobianWrtVar
            (EAdd
              (EFun (PowBaseP1InfNotInt p) (EFun (PowBaseP01 (1.0/p)) (EVar "x")))
              (EFun Neg (EVar "x")))
            [("x", x)]
            "x")
          (zeroMat $ rows x)
    it "(x^a)^(1/a) - x has zero jacobian (a even, a > 1)" $ property $
      forAll (pairsOf positiveVectors smallEvenPowersGreaterThanOne) $
        \(x, p) -> fpequalsMat
          (jacobianWrtVar
            (EAdd
              (EFun (PowBaseP1InfEven (round p :: Integer)) (EFun (PowBaseP01 (1.0/p)) (EVar "x")))
              (neg (EVar "x")))
            [("x", x)]
            "x")
          (zeroMat $ rows x)

max_stars_repo_name: chrisnc/hvx

max_issues_repo_name: chrisnc/hvx

max_forks_repo_name: chrisnc/hvx