text:
{-
   Copyright 2016, Dominic Orchard, Andrew Rice, Mistral Contrastin, Matthew Danish

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
-}

{-
  Units of measure extension to Fortran: backend
-}

{-# LANGUAGE TupleSections #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ForeignFunctionInterface #-}

module Camfort.Specification.Units.InferenceBackendFlint where

import           Control.Monad
import           Data.List (partition)
-- import           Debug.Trace (trace, traceM, traceShowM)
import           Foreign
import           Foreign.C.Types
import           Numeric.LinearAlgebra (atIndex, rows, cols)
import qualified Numeric.LinearAlgebra as H
import           System.IO.Unsafe (unsafePerformIO)

foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_init" fmpz_mat_init :: Ptr FMPZMat -> CLong -> CLong -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_set" fmpz_mat_set :: Ptr FMPZMat -> Ptr FMPZMat -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_entry" fmpz_mat_entry :: Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
foreign import ccall unsafe "flint/fmpz.h fmpz_set_si" fmpz_set_si :: Ptr CLong -> CLong -> IO ()
foreign import ccall unsafe "flint/fmpz.h fmpz_get_si" fmpz_get_si :: Ptr CLong -> IO CLong
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_clear" fmpz_mat_clear :: Ptr FMPZMat -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_print_pretty" fmpz_mat_print_pretty :: Ptr FMPZMat -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_mul" fmpz_mat_mul :: Ptr FMPZMat -> Ptr FMPZMat -> Ptr FMPZMat -> IO ()

-- fmpz_mat_window_init(fmpz_mat_t window, const fmpz_mat_t mat, slong r1, slong c1, slong r2, slong c2)
--
-- Initializes the matrix window to be an r2 - r1 by c2 - c1 submatrix
-- of mat whose (0,0) entry is the (r1, c1) entry of mat. The memory
-- for the elements of window is shared with mat.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_window_init" fmpz_mat_window_init :: Ptr FMPZMat -> Ptr FMPZMat -> CLong -> CLong -> CLong -> CLong -> IO ()

-- Frees the window (leaving underlying matrix alone).
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_window_clear" fmpz_mat_window_clear :: Ptr FMPZMat -> IO ()

-- r <- fmp_mat_rref B den A
--
-- Uses fraction-free Gauss-Jordan elimination to set (B, den) to the
-- reduced row echelon form of A and returns the rank of A. Aliasing
-- of A and B is allowed. r is rank of A.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_rref" fmpz_mat_rref :: Ptr FMPZMat -> Ptr CLong -> Ptr FMPZMat -> IO CLong

-- r <- fmp_mat_inv B den A
--
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_inv" fmpz_mat_inv :: Ptr FMPZMat -> Ptr CLong -> Ptr FMPZMat -> IO CLong

-- fmpz_mat_hnf H A
--
-- H is the Hermite Normal Form of A.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_hnf" fmpz_mat_hnf :: Ptr FMPZMat -> Ptr FMPZMat -> IO ()

foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_rank" fmpz_mat_rank :: Ptr FMPZMat -> IO CLong

data FMPZMat

instance Storable FMPZMat where
  sizeOf _ = 4 * sizeOf (undefined :: CLong)
  alignment _ = alignment (undefined :: CLong)
  peek _ = undefined
  poke _ = undefined

-- testFlint :: IO ()
-- testFlint = do
--   traceM "***********************8 testFlint 8***********************"
--   alloca $ \ a -> do
--     alloca $ \ b -> do
--       alloca $ \ den -> do
--         let n = 10
--         fmpz_mat_init a n n
--         fmpz_mat_init b n n
--         forM_ [0..n-1] $ \ i -> do
--           forM_ [0..n-1] $ \ j -> do
--             e <- fmpz_mat_entry a i j
--             fmpz_set_si e (2 * i + j)

--         -- fmpz_mat_mul b a a
--         r <- fmpz_mat_rref b den a
--         fmpz_mat_print_pretty a
--         fmpz_mat_print_pretty b
--         d <- peek den
--         traceM $ "r = " ++ show r ++ " den = " ++ show d
--         fmpz_mat_hnf b a
--         fmpz_mat_print_pretty b
--         fmpz_mat_clear a
--         fmpz_mat_clear b

rref :: H.Matrix Double -> (H.Matrix Double, Int, Int)
rref m = unsafePerformIO $ do
  alloca $ \ outputM -> do
    alloca $ \ inputM -> do
      alloca $ \ den -> do
        let numRows = fromIntegral $ rows m
        let numCols = fromIntegral $ cols m
        fmpz_mat_init outputM numRows numCols
        fmpz_mat_init inputM numRows numCols
        forM_ [0..numRows-1] $ \ i -> do
          forM_ [0..numCols-1] $ \ j -> do
            e <- fmpz_mat_entry inputM i j
            fmpz_set_si e (floor (m `atIndex` (fromIntegral i, fromIntegral j)))
        r <- fmpz_mat_rref outputM den inputM
        d <- peek den

        -- DEBUG:
        -- fmpz_mat_print_pretty outputM
        -- traceM $ "r = " ++ show r ++ " den = " ++ show d
        --

        lists <- forM [0..numRows-1] $ \ i -> do
          forM [0..numCols-1] $ \ j -> do
            e <- fmpz_mat_entry outputM i j
            fromIntegral `fmap` fmpz_get_si e
        let m' = H.fromLists lists
        fmpz_mat_clear inputM
        fmpz_mat_clear outputM
        return (m', fromIntegral d, fromIntegral r)

hnf :: H.Matrix Double -> H.Matrix Double
hnf m = unsafePerformIO $ do
  alloca $ \ outputM -> do
    alloca $ \ inputM -> do
      let numRows = fromIntegral $ rows m
      let numCols = fromIntegral $ cols m
      fmpz_mat_init outputM numRows numCols
      fmpz_mat_init inputM numRows numCols
      forM_ [0..numRows-1] $ \ i -> do
        forM_ [0..numCols-1] $ \ j -> do
          e <- fmpz_mat_entry inputM i j
          fmpz_set_si e (floor (m `atIndex` (fromIntegral i, fromIntegral j)))
      fmpz_mat_hnf outputM inputM
      r <- fmpz_mat_rank outputM

      -- DEBUG:
      -- fmpz_mat_print_pretty outputM
      -- traceM $ "rank = " ++ show r
      --

      lists <- forM [0..fromIntegral r-1] $ \ i -> do
        forM [0..numCols-1] $ \ j -> do
          e <- fmpz_mat_entry outputM i j
          fromIntegral `fmap` fmpz_get_si e
      let m' = H.fromLists lists
      fmpz_mat_clear inputM
      fmpz_mat_clear outputM
      return m'

inv :: H.Matrix Double -> Maybe (H.Matrix Double, Int)
inv m = unsafePerformIO $ do
  alloca $ \ outputM -> do
    alloca $ \ inputM -> do
      alloca $ \ den -> do
        let numRows = fromIntegral $ rows m
        let numCols = fromIntegral $ cols m
        fmpz_mat_init outputM numRows numCols
        fmpz_mat_init inputM numRows numCols
        forM_ [0..numRows-1] $ \ i -> do
          forM_ [0..numCols-1] $ \ j -> do
            e <- fmpz_mat_entry inputM i j
            fmpz_set_si e (floor (m `atIndex` (fromIntegral i, fromIntegral j)))
        r <- fmpz_mat_inv outputM den inputM
        d <- peek den

        -- DEBUG:
        fmpz_mat_print_pretty outputM
        -- traceM $ "r = " ++ show r ++ " den = " ++ show d
        --

        lists <- forM [0..numRows-1] $ \ i -> do
          forM [0..numCols-1] $ \ j -> do
            e <- fmpz_mat_entry outputM i j
            fromIntegral `fmap` fmpz_get_si e
        let m' = H.fromLists lists
        fmpz_mat_clear inputM
        fmpz_mat_clear outputM
        if r == 1 then
          return $ Just (m', fromIntegral d)
        else
          return Nothing

withMatrix :: H.Matrix Double -> ((CLong, CLong, Ptr FMPZMat) -> IO b) -> IO b
withMatrix m f = do
  alloca $ \ outputM -> do
    let numRows = fromIntegral $ rows m
    let numCols = fromIntegral $ cols m
    fmpz_mat_init outputM numRows numCols
    forM_ [0 .. numRows - 1] $ \ i ->
      forM_ [0 .. numCols - 1] $ \ j -> do
        e <- fmpz_mat_entry outputM i j
        fmpz_set_si e (floor (m `atIndex` (fromIntegral i, fromIntegral j)) :: CLong)
    x <- f (numRows, numCols, outputM)
    fmpz_mat_clear outputM
    return x

withBlankMatrix :: CLong -> CLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix numRows numCols f = do
  alloca $ \ outputM -> do
    fmpz_mat_init outputM (fromIntegral numRows) (fromIntegral numCols)
    x <- f outputM
    fmpz_mat_clear outputM
    return x

withWindow :: Ptr FMPZMat -> CLong -> CLong -> CLong -> CLong -> (Ptr FMPZMat -> IO b) -> IO b
withWindow underM r1 c1 r2 c2 f = do
  alloca $ \ window -> do
    fmpz_mat_window_init window underM r1 c1 r2 c2
    x <- f window
    fmpz_mat_window_clear window
    return x

flintToHMatrix :: CLong -> CLong -> Ptr FMPZMat -> IO (H.Matrix Double)
flintToHMatrix numRows numCols flintM = do
  lists <- forM [0..numRows-1] $ \ i -> do
    forM [0..numCols-1] $ \ j -> do
      e <- fmpz_mat_entry flintM i j
      fromIntegral `fmap` fmpz_get_si e
  return $ H.fromLists lists

pokeM :: Ptr FMPZMat -> CLong -> CLong -> CLong -> IO ()
pokeM flintM i j v = do
  e <- fmpz_mat_entry flintM i j
  fmpz_set_si e v

peekM :: Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM flintM i j = do
  e <- fmpz_mat_entry flintM i j
  fmpz_get_si e

copyMatrix :: Ptr FMPZMat -> Ptr FMPZMat -> CLong -> CLong -> CLong -> CLong -> IO ()
copyMatrix m1 m2 r1 c1 r2 c2 =
  forM_ [r1 .. r2-1] $ \ i ->
    forM_ [c1 .. c2-1] $ \ j ->
      peekM m2 i j >>= pokeM m1 i j

-- copyMatrix' m1 m2 r1 c1 r2 c2 =
--   withWindow m2 r1 c1 r2 c2 $ \ w2 ->
--     withWindow m1 r1 c1 r2 c2 $ \ w1 -> do
--       fmpz_mat_set w1 w2

--------------------------------------------------
-- 'normalising' Hermite Normal Form, for lack of a better name
--
-- The problem with HNF is that it will happily return a matrix where
-- the leading-coefficient diagonal contains integers > 1.
--
-- In some cases the leading-coefficient does not divide some of the
-- remaining numbers in the row.
--
-- This corresponds to the case where one of the solutions would be
-- fractional if we were dealing with rational matrices.
--
-- But we don't want fractional solutions. We need to bump the matrix
-- so that this situation does not arise.
--
-- This tends to happen due to implicit polymorphism.
--
-- We do this by adding another column that copies the column with the
-- problematic leading-coefficient.
--
-- We then set up a new row expressing the constraint that the old
-- column is equal to the new column, 1:1.
--
-- This prevents a 'fractional' solution.
--
-- Along the way we look for leading-coefficients > 1 that do divide
-- their entire row. Then we simply apply elementary row scaling to
-- the row so that the leading-coefficients is 1.
--
-- The result is a matrix where the leading-coefficients of all the
-- rows that matter are 1. It's not quite a 'reduced' form though
-- because it can be bigger than the original matrix.
--
-- As a result we also return a list of columns that were cloned this
-- way.
--
-- The resulting matrix can have non-zeroes above the
-- leading-coefficients in the same column, so I apply some elementary
-- row operations to fix that and bring things into RREF.
--
-- Running time is computation of Hermite Normal Form twice, plus
-- construction of a slightly larger matrix (possibly), plus scanning
-- through the matrix for leading-coefficients, and then again to
-- divide them out sometimes.

normHNF :: H.Matrix Double -> (H.Matrix Double, [Int])
normHNF m
  | (rows m, cols m) == (1, 1) = (H.ident (if H.atIndex m (0, 0) /= 0 then 1 else 0), [])
  | otherwise = (m', indices)
  where
    numCols = cols m
    indexLookup j | j < numCols = j
                  | otherwise = indices !! (j `mod` numCols)
    indices = map indexLookup $ concatMap snd rs1
    (rs1, (m',[]):_) = break (null . snd) results
    results = tail $ iterate (normHNF' . fst) (m, [])

normHNF' :: H.Matrix Double -> (H.Matrix Double, [Int])
normHNF' m = fmap (map fromIntegral) . unsafePerformIO $ withMatrix m normhnf

normhnf :: (CLong, CLong, Ptr FMPZMat) -> IO (H.Matrix Double, [CLong])
normhnf (numRows, numCols, inputM) = do
  withBlankMatrix numRows numCols $ \ outputM -> do
    fmpz_mat_hnf outputM inputM
    rank <- fmpz_mat_rank outputM
    -- scale the rows so that it is a RREF, returning the indices
    -- where leading-coefficients had to be scaled to 1.
    indices <- elemrowscale outputM rank numCols
    -- HNF allows non-zeroes above the leading-coefficients that
    -- were greater than 1, so also fix that to bring into RREF.
    elemrowadds outputM rank numCols indices

    -- column indices of leading co-efficients > 1
    lcoefs <- filter ((> 1) . head . snd) <$> forM [0 .. rank-1] (\ i -> do
      cs <- zip [0..] `fmap` sequence [ peekM outputM i j | j <- [0 .. numCols-1] ]
      let (_, (j, lcoef):rest) = span ((== 0) . snd) cs
      return ((i, j), lcoef:map snd rest))

    -- split the identified rows into two categories: those that can
    -- be divided out by the leading co-efficient, and those that
    -- cannot.
    let (multCands, consCands) = partition (\ (_, lcoef:rest) -> all ((== 0) . (`rem` lcoef)) rest) lcoefs

    -- apply elementary row scaling
    forM_ multCands $ \ ((i, j), lcoef:_) ->
      forM_ [j..numCols - 1] $ \ j' -> do
        x <- peekM outputM i j'
        pokeM outputM i j' $ x `div` lcoef

    -- identify columns that need additional constraints generated and
    -- their associated value d, which is the value of the non-leading
    -- co-efficient divided by the GCD of the row.
    let genColCons ((_, j), lcoef:rest) = (j, minNLcoef `div` gcd lcoef minNLcoef)
          where
            restABS   = map abs rest
            minNLcoef = minimum (filter (/= 0) restABS)
        genColCons _ = error "normhnf: genColCons: impossible: missing leading co-efficient"

    let consCols = map genColCons consCands

    -- generate operations that poke 1.0 and -d into columns of the
    -- given row (matrix parameter supplied later); d is the value of
    -- the non-leading co-efficient that isn't divisible by the
    -- leading-coefficient, divided by the GCD of the row.
    let ops = [ \ flintM i -> do pokeM flintM i j 1
                                 pokeM flintM i (numCols + k) (-d)
              | ((j, d), k) <- zip consCols [0..] ]
    let numOps = fromIntegral (length ops)
    let numRows' = numOps + rank
    let numCols' = numOps + numCols

    withBlankMatrix numRows' numCols' $ \ outputM' -> do
      -- create larger matrix containing outputM as submatrix
      copyMatrix outputM' outputM 0 0 rank numCols

      -- apply the operations on the extra space to write the
      -- additional rows and columns providing the additional
      -- constraints
      forM_ (zip [rank..] ops) $ \ (i, op) -> op outputM' i

      -- re-run HNF
      withBlankMatrix numRows' numCols' $ \ outputM'' -> do
        fmpz_mat_hnf outputM'' outputM'
        rank' <- fmpz_mat_rank outputM''
        -- scale the rows so that it is a RREF, returning the indices
        -- where leading-coefficients had to be scaled to 1.
        indices' <- elemrowscale outputM'' rank' numCols'
        -- HNF allows non-zeroes above the leading-coefficients that
        -- were greater than 1, so also fix that to bring into RREF.
        elemrowadds outputM'' rank' numCols' indices'
        -- convert back to HMatrix form
        h <- flintToHMatrix rank' numCols' outputM''
        return (h, map fst consCols)

-- find leading-coefficients that are greater than 1 and scale those
-- rows accordingly to reach RREF
--
-- precondition: matrix outputM is in HNF
elemrowscale :: Ptr FMPZMat -> CLong -> CLong -> IO [(CLong, CLong)]
elemrowscale outputM rank numCols = do
  -- column indices of leading co-efficients > 1
  lcoefs <- filter ((> 1) . head . snd) <$> forM [0 .. rank-1] (\ i -> do
    cs <- zip [0..] `fmap` sequence [ peekM outputM i j | j <- [0 .. numCols-1] ]
    let (_, (j, lcoef):rest) = span ((== 0) . snd) cs
    return ((i, j), lcoef:map snd rest))

  let multCands = filter (\ (_, lcoef:rest) -> all ((== 0) . (`rem` lcoef)) rest) lcoefs

  -- apply elementary row scaling
  forM multCands $ \ ((i, j), lcoef:_) -> do
    forM_ [j..numCols - 1] $ \ j' -> do
      x <- peekM outputM i j'
      pokeM outputM i j' $ x `div` lcoef
    return (i, j)

-- use indices to guide elementary row additions to reach RREF
--
-- precondition: matrix outputM is in HNF save for work done by
-- elemrowscale, and indices is a list of coordinates where we have
-- just scaled the leading-coefficient to 1 and now we must look to
-- see if there are any non-zeroes in the column above the
-- leading-coefficient, because that is allowed by HNF.
elemrowadds :: Ptr FMPZMat -> CLong -> CLong -> [(CLong, CLong)] -> IO ()
elemrowadds outputM _ numCols indices = do
  -- look for non-zero members of the columns above the
  -- leading-coefficient and wipe them out.
  forM_ indices $ \ (lcI, lcJ) -> do
    let j = lcJ
    forM_ [0..lcI-1] $ \ i -> do
      -- (i, j) ranges over the elements of the column above leading
      -- co-efficient (lcI, lcJ).
      x <- peekM outputM i j
      if x == 0 then
        pure () -- nothing to do at row i
      else do
        -- (i, j) is non-zero and row i must be cancelled x times
        let sf = x -- scaling factor = non-zero magnitude
        forM_ [lcJ..numCols-1] $ \ j' -> do
          -- (i,   j') ranges over the row where we discovered a non-zero
          -- (lcI, j') ranges over the row with the leading co-efficient
          x1 <- peekM outputM i j'
          x2 <- peekM outputM lcI j'
          -- add lower row scaled by sf to upper row
          pokeM outputM i j' (x1 - x2 * sf)

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

-- m1 :: H.Matrix Double
-- m1 = (8><6)
--  [ 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 0.0, 0.0, 0.0,  1.0, -4.0,  0.0
--  , 0.0, 0.0, 0.0,  0.0,  4.0, -3.0 ]

-- m2 :: H.Matrix Double
-- m2 = (8><6)
--  [ 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 0.0, 0.0, 0.0,  1.0, -6.0,  0.0
--  , 0.0, 0.0, 0.0,  0.0,  6.0, -4.0 ]

max_stars_repo_name: apthorpe/camfort

max_issues_repo_name: apthorpe/camfort

max_forks_repo_name: apthorpe/camfort


text:
module AND where

import Numeric.LinearAlgebra
import AI.HNN.FF.Network

samples :: Samples Double
samples = [ (fromList [0,0], fromList [1,0])
          , (fromList [0,1], fromList [1,0])
          , (fromList [1,0], fromList [1,0])
          , (fromList [1,1], fromList [0,1]) ]

run :: IO ()
run = do
    n <- createNetwork 2 [2] 2
    mapM_ (print . output n tanh . fst) samples
    putStrLn "----------------"
    let n' = trainNTimes 1000 0.8 tanh tanh' n samples
    mapM_ (print . output n' tanh . fst) samples

ext: hs

max_stars_repo_name: yohad/networks

max_issues_repo_name: yohad/networks

max_forks_repo_name: yohad/networks


text:
import Data.Complex

-- Polynomial type
type Poly a = [Complex a]

__fftAux :: (RealFloat a) => a -> Poly a -> Poly a
__fftAux _ [] = []
__fftAux _ [x] = [x]
__fftAux c poly = uncurry (++) res
  where
    n = length poly
    ls = __fftAux c [ x | (x, i) <- zip poly [1..],  odd i ]
    rs = __fftAux c [ x | (x, i) <- zip poly [1..], even i ]
    res = combine ls rs 1
      where
        ang = c * 2 * pi / fromIntegral n
        wn = cos ang :+ sin ang
        combine [] [] _ = ([], [])
        combine (x:xs) (y:ys) w = (x + w * y : rxs, x - w * y : rys)
          where
            (rxs, rys) = combine xs ys (w * wn)

__closestPowerOf2 :: Int -> Int
__closestPowerOf2 n = compute 1
  where
    compute x
      | x >= n = x
      | otherwise = compute $ 2 * x

__fftWrapper :: (RealFloat a) => a -> Poly a -> Poly a
__fftWrapper c poly
  | oddf (length poly) == 1 = __fftAux c poly
  | otherwise = error "FFT: size is not a power of 2!"
    where
      oddf x
        | odd x = x
        | otherwise = oddf $ x `div` 2

fft :: (RealFloat a) => Poly a -> Poly a
fft = __fftWrapper 1

ifft :: (RealFloat a) => Poly a -> Poly a
ifft poly = map (/n) $ __fftWrapper (-1) poly
  where
    n = fromIntegral $ length poly

multiply :: (RealFloat a) => Poly a -> Poly a -> Poly a
multiply ps qs = ifft $ zipWith (*) (fft $ pad ps) (fft $ pad qs)
  where
    size = __closestPowerOf2 $ length ps + length qs - 1
    pad xs = xs ++ replicate (size - length xs) 0

multiplyReal :: (RealFloat a) => [a] -> [a] -> [a]
multiplyReal ps qs = map realPart $ multiply (conv ps) (conv qs)
  where
    conv = map (:+ 0)

multiplyIntegral :: (Integral a) => [a] -> [a] -> [a]
multiplyIntegral ps qs = map round $ multiplyReal (conv ps) (conv qs)
  where
    conv = map fromIntegral

max_stars_repo_name: anurudhp/CPHaskell

max_issues_repo_name: anurudhp/CPHaskell

max_forks_repo_name: anurudhp/CPHaskell


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE KindSignatures #-}

#if __GLASGOW_HASKELL__ >= 706
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE ExplicitNamespaces #-}
#endif

#if __GLASGOW_HASKELL__ >= 806
{-# LANGUAGE NoStarIsType #-}
#endif

{-# OPTIONS_GHC -fno-warn-deprecations #-}

{-|
Exports modules that Data.Orphans needs. Because Data.Orphans uses several modules
that only need to be in scope for certain versions of GHC, exporting all of the modules
separately eliminates the need to use CPP pragmas for GHC-version-specific imports.
This makes it much easier to be -Wall-compliant.

Note that this module does not export any modules that could introduce name clashes.
-}
module Data.Orphans.Prelude
#if MIN_VERSION_base(4,16,0)
    () where
#else
    ( module OrphansPrelude
    , realPart
    , imagPart
    , Typeable1
    , Typeable2
    , Typeable3
    , Typeable4
    , Typeable5
    , Typeable6
    , Typeable7
    ) where

import Control.Applicative as OrphansPrelude
import Control.Arrow as OrphansPrelude hiding ((+++), first, loop, second)
import Control.Category as OrphansPrelude hiding ((.), id)
import Control.Concurrent.QSem as OrphansPrelude
import Control.Monad as OrphansPrelude hiding (mapM, sequence)
import Control.Monad.Fix as OrphansPrelude

import Data.Bits as OrphansPrelude
import Data.Char as OrphansPrelude
import Data.Complex as OrphansPrelude (Complex(..))
import Data.Data as OrphansPrelude (Data(..), Constr, DataType, mkConstr, mkDataType)
import Data.Fixed as OrphansPrelude
import Data.Int as OrphansPrelude
import Data.List as OrphansPrelude (genericLength)
import Data.Monoid as OrphansPrelude
  hiding ( Any(..)
# if MIN_VERSION_base(4,5,0)
         , (<>)
# endif
         )
import Data.String as OrphansPrelude (IsString)
import Data.Typeable as OrphansPrelude (Typeable)
import Data.Version as OrphansPrelude
import Data.Word as OrphansPrelude

import Foreign.C.Error as OrphansPrelude
import Foreign.C.Types as OrphansPrelude
import Foreign.Marshal.Pool as OrphansPrelude
import Foreign.Ptr as OrphansPrelude
import Foreign.Storable as OrphansPrelude

import GHC.Arr as OrphansPrelude (Ix(..))
import GHC.Base as OrphansPrelude
import GHC.Conc as OrphansPrelude
import GHC.Desugar as OrphansPrelude (AnnotationWrapper)
import GHC.ForeignPtr as OrphansPrelude
import GHC.IO.Buffer as OrphansPrelude
import GHC.IO.BufferedIO as OrphansPrelude (BufferedIO)
import GHC.IO.Device as OrphansPrelude (IODevice, IODeviceType(..), RawIO)
import GHC.IO.Encoding as OrphansPrelude
import GHC.IO.Handle as OrphansPrelude
import GHC.IO.Handle.Types as OrphansPrelude
import GHC.Real as OrphansPrelude (Ratio(..), (%))
import GHC.ST as OrphansPrelude

import Numeric as OrphansPrelude (showHex)

import System.Console.GetOpt as OrphansPrelude
import System.IO as OrphansPrelude
import System.Posix.Internals as OrphansPrelude
import System.Posix.Types as OrphansPrelude

import Text.Printf as OrphansPrelude
import Text.Read as OrphansPrelude

import Unsafe.Coerce as OrphansPrelude (unsafeCoerce)

# if defined(mingw32_HOST_OS)
import GHC.IO.Encoding.CodePage.Table as OrphansPrelude
# endif

# if MIN_VERSION_base(4,4,0)
import Control.Monad.Zip as OrphansPrelude
import GHC.Fingerprint as OrphansPrelude
import GHC.IO.Encoding.Failure as OrphansPrelude

#  if !defined(mingw32_HOST_OS) && !defined(__GHCJS__)
import GHC.Event as OrphansPrelude
#   if !(MIN_VERSION_base(4,7,0))
      hiding (step)
#   endif
#  endif
# endif

# if MIN_VERSION_base(4,5,0)
import GHC.Stack as OrphansPrelude
import GHC.Stats as OrphansPrelude
# endif

# if MIN_VERSION_base(4,6,0)
import Data.Ord as OrphansPrelude (Down(..))
import GHC.GHCi as OrphansPrelude
import GHC.TypeLits as OrphansPrelude hiding (type (*))
# else
import GHC.Exts as OrphansPrelude (Down(..))
# endif

# if MIN_VERSION_base(4,6,0) && !(MIN_VERSION_base(4,8,2))
import GHC.IP as OrphansPrelude
# endif

# if MIN_VERSION_base(4,7,0)
import Data.Proxy as OrphansPrelude
import Data.Type.Coercion as OrphansPrelude (Coercion, TestCoercion)
import Data.Type.Equality as OrphansPrelude ((:~:)(..), TestEquality(..))
import GHC.Exts as OrphansPrelude (IsList(..))
import Text.Read.Lex as OrphansPrelude (Number)
# else
import Control.Concurrent.SampleVar as OrphansPrelude
# endif

# if MIN_VERSION_base(4,8,0)
import Data.Bifunctor as OrphansPrelude
import Data.Functor.Identity as OrphansPrelude
#endif

# if MIN_VERSION_base(4,9,0)
import Data.Functor.Classes as OrphansPrelude
import Data.Functor.Compose as OrphansPrelude
import Data.Semigroup as OrphansPrelude (stimesMonoid)
# endif

# if MIN_VERSION_base(4,9,0) && !(MIN_VERSION_base(4,12,0))
import Data.List.NonEmpty as OrphansPrelude (NonEmpty(..))
import Data.Semigroup as OrphansPrelude (Semigroup(..))
# endif

# if !(MIN_VERSION_base(4,11,0))
import Data.Typeable ( Typeable1, Typeable2, Typeable3, Typeable4
                     , Typeable5, Typeable6, Typeable7 )
# endif

# if MIN_VERSION_base(4,4,0)
realPart, imagPart :: Complex a -> a
# else
realPart, imagPart :: RealFloat a => Complex a -> a
# endif
realPart (x :+ _) = x
imagPart (_ :+ y) = y

# if MIN_VERSION_base(4,11,0)
type Typeable1 (a :: Type -> Type)                                                 = Typeable a
type Typeable2 (a :: Type -> Type -> Type)                                         = Typeable a
type Typeable3 (a :: Type -> Type -> Type -> Type)                                 = Typeable a
type Typeable4 (a :: Type -> Type -> Type -> Type -> Type)                         = Typeable a
type Typeable5 (a :: Type -> Type -> Type -> Type -> Type -> Type)                 = Typeable a
type Typeable6 (a :: Type -> Type -> Type -> Type -> Type -> Type -> Type)         = Typeable a
type Typeable7 (a :: Type -> Type -> Type -> Type -> Type -> Type -> Type -> Type) = Typeable a
# endif
#endif

max_stars_repo_name: haskell-compat/base-orphans

max_issues_repo_name: haskell-compat/base-orphans

max_forks_repo_name: haskell-compat/base-orphans


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}

module Examples.DelayedSamplingAbstract where

import GHC.TypeNats (KnownNat)
import Numeric.LinearAlgebra.Static (L, R, Sym, konst, sym, eye)

import Control.Monad.Bayes.Sampler (SamplerIO)
import Util.MStream (MStream)
import qualified Util.MStream as MS

import qualified Distributions as D
import DSProg
import SymbolicDistr

kalman :: Bool -> Double -> Double -> [Double] -> MStream M (Expr Double) (Expr Double)
kalman forgetb m b yobs = do
  x1 <- MS.lift $ sample (normal 0 1)
  MS.lift $ observe (normal x1 1) (head yobs)
  MS.yield x1
  loop 2 x1 (tail yobs)
  where
  forgett x = if forgetb then MS.lift (forgetE x) else pure ()
  loop :: Integer -> Expr Double -> [Double] -> MStream M (Expr Double) (Expr Double)
  loop t xpredt [] = pure xpredt
  loop t xpredt (y : ys) = do
    xt <- MS.lift $ sample (normal (Const m * xpredt + Const b) 1)
    MS.lift $ observe (normal xt 1) y
    MS.yield xt
    forgett xpredt
    loop (t + 1) xt ys

exampleKalman :: SamplerIO (Double, Maybe (Result Double))
exampleKalman = evalM (justRun (kalman True 1 1 [1..]))

kalmanMV :: forall n. KnownNat n => Bool -> L n n -> R n -> [R n] -> MStream M (Expr (R n)) (Expr (R n))
kalmanMV forgetb m b yobs = do
  x1 <- MS.lift $ sample (mvNormal (Const (konst 0)) (sym eye))
  MS.lift $ observe (mvNormal x1 (sym eye)) (head yobs)
  MS.yield x1
  loop 2 x1 (tail yobs)
  where
  forgett x = if forgetb then MS.lift (forgetE x) else pure ()
  loop :: Integer -> Expr (R n) -> [R n] -> MStream M (Expr (R n)) (Expr (R n))
  loop t xpredt [] = pure xpredt
  loop t xpredt (y : ys) = do
    xt <- MS.lift $ sample (mvNormal (MVMul (Const m) xpredt + Const b) (sym eye))
    MS.lift $ observe (mvNormal xt (sym eye)) y
    MS.yield xt
    forgett xpredt
    loop (t + 1) xt ys

exampleKalmanMV :: SamplerIO (Double, Maybe (Result (R 3)))
exampleKalmanMV = evalM (justRun (kalmanMV True eye (konst 1 :: R 3) (map konst [1..])))

betaBernoulli :: Double -> Double -> [Bool] -> MStream M (Expr Double) (Expr Double)
betaBernoulli a b xs = do
  p <- MS.lift $ sample (betaSD a b)
  loop p xs
  where
  loop :: Expr Double -> [Bool] -> MStream M (Expr Double) (Expr Double)
  loop p [] = pure p
  loop p (y : ys) = do
    MS.lift $ observe (bernoulliSD p) y
    MS.yield p
    loop p ys

exampleBetaBernoulli :: SamplerIO (Double, Maybe (Result Double))
exampleBetaBernoulli = evalM (justRun (betaBernoulli 1 1 (cycle [True,  True, False])))

yieldValue :: Expr p -> MStream M (Expr p) ()
yieldValue e = do
  x <- MS.lift $ force e
  MS.yield (Const x)

delay_triplet :: Double -> MStream M (Expr Double) (Expr Double)
delay_triplet zobs = do
  x <- MS.lift $ sample (normal 0 1)
  y <- MS.lift $ sample (normal x 1)
  MS.lift $ observe (normal y 1) zobs
  MS.yield y
  MS.yield x -- XXX: UNSOUND! NEED TO FIX GET_MARGINAL to force here!!!
  yieldValue y
  yieldValue x
  pure x

delay_iid :: [Double] -> MStream M (Expr Double) (Expr Double)
delay_iid yobs = do
  x <- MS.lift $ sample (normal 0 1)
  loop x yobs
  where
  loop x [] = yieldValue x >> pure x
  loop x (y : ys) = do
    MS.lift $ observe (normal x 1) y
    MS.yield x
    loop x ys

louis :: M ()
louis = do
  b <- sample $ bernoulliSD (Const 0.2)
  o <- sample $ normal (Const 0.2) 1
  b' <- force b
  if b' then D.factor (-0.1) else pure ()

max_stars_repo_name: psg-mit/probzelus-haskell

max_issues_repo_name: psg-mit/probzelus-haskell

max_forks_repo_name: psg-mit/probzelus-haskell


text:
module NeuralNet.TrainSpec (trainSpec) where

import Test.Hspec
import Numeric.LinearAlgebra
import NeuralNet.Net
import NeuralNet.Activation
import NeuralNet.Example
import NeuralNet.Layer
import NeuralNet.Train

trainSpec :: SpecWith ()
trainSpec =
  describe "NeuralNet.Train" $ do
    describe "nnBackward" $ do
      it "calculates correctly for simple example" $
        let
          layer1W = [-1.31386475, 0.88462238, 0.88131804, 1.70957306
                    , 0.05003364, -0.40467741, -0.54535995, -1.54647732
                    , 0.98236743, -1.10106763, -1.18504653, -0.2056499]
          layer1B = [1.48614836, 0.23671627, -1.02378514]
          layer2W = [-1.02387576, 1.12397796, -0.13191423]
          layer2B = [-1.62328545]
          nn = buildNNFromList (4, [LayerDefinition ReLU 3, LayerDefinition Sigmoid 1]) (layer1W ++ layer1B ++ layer2W ++ layer2B)
          examples = createExampleSet [ ([0.09649747, -0.2773882, -0.08274148, -0.04381817], 1)
                                      , ([-1.8634927, -0.35475898, -0.62700068, -0.47721803], 0)]

          inputs = exampleSetX examples
          layer1As = fromLists [[1.97611078, -1.24412333], [-0.62641691, -0.80376609], [-2.41908317, -0.92379202]]
          layer1Zs = fromLists [[-0.7129932, 0.62524497], [-0.16051336, -0.76883635], [-0.23003072, 0.74505627]]
          layer2As = fromLists [[1.78862847, 0.43650985]]
          layer2Zs = fromLists [[0.64667545, -0.35627076]]
          steps = [createInputForwardPropStep inputs, createForwardPropStep layer1As layer1Zs, createForwardPropStep layer2As layer2Zs]

          dW1 = fromLists [ [0.4101000190122487, 7.807203346853248e-2, 0.13798443685274048, 0.10502167417988162]
                          , [0, 0, 0, 0]
                          , [5.2836516249770524e-2, 1.0058654166727896e-2, 1.77776556985907e-2, 1.3530795262454464e-2]]
          db1 = fromLists [[-0.22007063350033446], [0], [-2.835348711039787e-2]]
          dW2 = fromLists [[-0.39202432174003965, -0.13325854898009654, -4.601088848081533e-2]]
          db2 = fromLists [[0.15187860742650333]]
        in
          nnBackward nn steps examples `shouldBe` [(dW1, db1), (dW2, db2)]

    describe "updateNNParams" $ do
      it "calculates correctly for logreg case" $
        let
          nn = buildNNFromList (2, [LayerDefinition ReLU 1]) [1, 2, 3]
          grads = [(fromLists [[10, 100]], fromLists [[1000]])]
          learningRate = 0.1
          expected = buildNNFromList (2, [LayerDefinition ReLU 1]) [0, -8, -97]
        in
          updateNNParams nn grads learningRate `shouldBe` expected

max_stars_repo_name: danielholmes/neural-net

max_issues_repo_name: danielholmes/neural-net

max_forks_repo_name: danielholmes/neural-net


text:
module Eval where

import Control.Monad.Except
import Data.Array
import Data.Complex
import Data.IORef
import Data.Maybe (isNothing)
import Data.Ratio
import DataTypes
       (Arity(..), Env, EnvFrame(..), IOPrimitiveFunc, IOThrowsError,
        LispError(..), LispVal(..), PrimitiveFunc, ThrowsError,
        extractValue, showVal, trapError)
import Parse
import ParserCombinators
import Paths_hascheme (getDataFileName)
import Primitives (eqv, ioPrimitives, primitives)
import System.Console.Haskeline
import System.IO
import Util (liftThrows)
import Vector (outOfBoundsError)

--------------
-- Eval
--------------
eval :: Env -> LispVal -> IOThrowsError LispVal
eval _ val@Void = return val
eval _ val@(String _) = return val
eval _ val@(Character _) = return val
eval _ val@(Integer _) = return val
eval _ val@(Vector _) = return val
eval _ (Rational val) =
  if denominator val == 1
    then return $ Integer $ numerator val
    else return $ Rational val
eval _ val@(Float _) = return val
eval _ (Complex val) =
  if imagPart val == 0
    then return $ Float $ realPart val
    else return $ Complex val
eval _ val@(Bool _) = return val
eval env (Atom id') = getVar env id'
eval _ (List [Atom "quote", val]) = return val
eval env (List [Atom "if", pred', conseq, alt]) = do
  result <- eval env pred'
  case result of
    Bool False -> eval env alt
    _ -> eval env conseq
eval env (List ((Atom "cond"):xs)) = evalCond xs
    -- TODO: [test-expr => proc-expr]
  where
    evalCond :: IOPrimitiveFunc
    evalCond [] = return Void
    evalCond [List (Atom "else":xs')] = evalList env xs'
    evalCond (List (Atom "else":_):_) =
      throwError $ Default "cond: bad syntax (`else` clause must be last)"
    evalCond ((List (pred':conseqs)):xs') = do
      result <- eval env pred'
      case result of
        Bool False -> evalCond xs'
        _ ->
          case conseqs of
            [] -> return result
            _ -> evalList env conseqs
    evalCond a = throwError $ Default (show a)
eval env (List ((Atom "case"):(key:clauses))) = do
  result <- eval env key
  evalClauses result clauses
  where
    evalClauses :: LispVal -> IOPrimitiveFunc
    evalClauses _ [List []] = return Void
    evalClauses key' (List (datum:exprs):rest) = do
      match <- liftThrows $ inList [key', datum]
      case match of
        Bool True -> evalList env exprs
        Bool False -> evalClauses key' rest
        _ -> throwError $ Default "case: bad syntax"
    evalClauses _ _ = throwError $ Default "case: bad syntax"
    inList :: PrimitiveFunc
    inList [_, List []] = return $ Bool False
    inList [key', List (x:xs)] = do
      eq <- eqv [x, key']
      case eq of
        Bool True -> return $ Bool $ True
        _ -> inList [key, List xs]
    inList _ = throwError $ Default "case: bad syntax"
eval _ (List ((Atom "case"):_)) =
  throwError $ Default "case: bad syntax in: (case)"
eval env (List [Atom "set!", Atom var, form]) = eval env form >>= setVar env var
eval env (List [Atom "define", Atom var, form]) =
  eval env form >>= defineVar env var
eval env (List (Atom "define":List (Atom var:params'):body')) =
  makeNormalFunc var env params' body' >>= defineVar env var
eval env (List (Atom "define":DottedList (Atom var:params') varargs:body')) =
  makeVarArgs var varargs env params' body' >>= defineVar env var
eval env (List (Atom "lambda":List params':body')) =
  makeNormalFunc "λ" env params' body'
eval env (List (Atom "lambda":varargs@(Atom _):body')) =
  makeVarArgs "λ" varargs env [] body'
eval env (List (Atom "lambda":DottedList params' varargs:body')) =
  makeVarArgs "λ" varargs env params' body'
eval env (List [Atom "let", List pairs, body']) = do
  List atoms <- getHeads pairs
  _ <- mapM_ ensureAtom atoms
  List vals <- getTails pairs
  vals' <- mapM (eval env) vals
  env' <-
    liftIO $
    bindVars env (Prelude.zipWith (\a b -> (extractVar a, b)) atoms vals')
  eval env' body'
eval env (List [Atom "let*", List pairs, body']) = do
  List atoms <- getHeads pairs
  _ <- mapM_ ensureAtom atoms
  List vals <- getTails pairs
  env' <- buildEnv env atoms vals
  eval env' body'
  where
    buildEnv :: Env -> [LispVal] -> [LispVal] -> IOThrowsError Env
    buildEnv env' [] [] = return env'
    buildEnv env' (atomH:atomT) (valH:valT) = do
      val <- eval env' valH
      env'' <- liftIO $ bindVars env' [(extractVar atomH, val)]
      buildEnv env'' atomT valT
    buildEnv _ _ _ = throwError $ Default "let*: bad syntax"
eval env (List [Atom "letrec", List pairs, body']) = do
  List atoms <- getHeads pairs
  _ <- mapM_ ensureAtom atoms
  List vals <- getTails pairs
  env' <- liftIO $ bindVars env (map (\a -> (extractVar a, Void)) atoms)
  vals' <- mapM (eval env') vals
  _ <-
    mapM_
      (\(n, v) -> setVar env' n v)
      (Prelude.zipWith (\a b -> (extractVar a, b)) atoms vals')
  eval env' body'
eval env (List [Atom "or", expr1, expr2]) = do
  result <- eval env expr1
  case result of
    Bool True -> return result
    _ -> eval env expr2
eval env (List [Atom "and", expr1, expr2]) = do
  result <- eval env expr1
  case result of
    Bool False -> return result
    _ -> eval env expr2
eval env (List [Atom "load", String filename]) = do
  f <- load filename
  _ <- mapM_ (eval env) f
  return Void
eval env (List [Atom "load-print", String filename]) =
  load filename >>= fmap List . mapM (eval env)
eval env (List (Atom "vector-set!":args)) =
  case args of
    [Atom var, Integer n, v] -> do
      Vector vec <- getVar env var
      if n < (fromIntegral $ length vec)
        then setVar env var $ Vector $ vec // [(n, v)]
        else throwError $ outOfBoundsError "vector-ref" n vec
    [a, b, c] ->
      throwError $ TypeMismatch "vector, integer, integer" $ List [a, b, c]
    a -> throwError $ NumArgs (MinMax 3 3) (length args) a
eval env (List (function:args)) = do
  func <- eval env function
  argVals <- mapM (eval env) args
  apply func argVals
eval _ badForm = throwError $ BadSpecialForm "Unrecognized special form" badForm

evalList :: Env -> IOPrimitiveFunc
evalList _ [] = return Void
evalList env [a] = eval env a
evalList env (y:ys) = eval env y >> evalList env ys

getHeads :: IOPrimitiveFunc
getHeads [] = return $ List []
getHeads (List (x:_):ys) = do
  List result <- getHeads ys
  return $ List (x : result)
getHeads _ = throwError $ Default "Unexpected error in let"

getTails :: IOPrimitiveFunc
getTails [] = return $ List []
getTails (List [_, xs]:ys) = do
  List result <- getTails ys
  return $ List (xs : result)
getTails _ = throwError $ Default "Unexpected error in let"

ensureAtom :: LispVal -> IOThrowsError LispVal
ensureAtom n@(Atom _) = return n
ensureAtom _ = throwError $ Default "Type error"

extractVar :: LispVal -> String
extractVar (Atom atom) = atom

makeFunc ::
     Monad m
  => String
  -> Maybe String
  -> Env
  -> [LispVal]
  -> [LispVal]
  -> m LispVal
makeFunc name' varargs env params' body' =
  return $ Func name' (map showVal params') varargs body' env

makeNormalFunc :: String -> Env -> [LispVal] -> IOPrimitiveFunc
makeNormalFunc name' = makeFunc name' Nothing

makeVarArgs :: String -> LispVal -> Env -> [LispVal] -> IOPrimitiveFunc
makeVarArgs name' = (makeFunc name') . Just . showVal

isBoundLocal :: EnvFrame -> String -> IO Bool
isBoundLocal envFrame var =
  case lookup var (bindings envFrame) of
    Nothing -> return False
    Just _ -> return True

isBound :: Env -> String -> IO Bool
isBound envRef var = do
  env <- readIORef envRef
  localBound <- isBoundLocal env var
  if localBound
    then return True
    else case env of
           Global _ -> return False
           Frame {parent = parent'} -> isBound parent' var

getVar :: Env -> String -> IOThrowsError LispVal
getVar envRef var = do
  envFrame <- liftIO $ readIORef envRef
  case lookup var (bindings envFrame) of
    Just a -> liftIO $ readIORef a
    Nothing ->
      case envFrame of
        Global _ -> throwError $ UnboundVar "Getting an unbound variable" var
        Frame {parent = parent'} -> getVar parent' var

setVar :: Env -> String -> LispVal -> IOThrowsError LispVal
setVar envRef var value = do
  envFrame <- liftIO $ readIORef envRef
  case lookup var (bindings envFrame) of
    Just a -> do
      liftIO $ writeIORef a value
      return Void
    Nothing ->
      case envFrame of
        Global _ -> throwError $ UnboundVar "Getting an unbound variable" var
        Frame {parent = parent'} -> setVar parent' var value

defineVar :: Env -> String -> LispVal -> IOThrowsError LispVal
defineVar envRef var value = do
  envFrame <- liftIO $ readIORef envRef
  alreadyLocallyDefined <- liftIO $ isBoundLocal envFrame var
  if alreadyLocallyDefined
    then throwError $
         Default $ "Duplicate definition for identifier in: " ++ var
    else liftIO $ do
           valueRef <- newIORef value
           env <- readIORef envRef
           liftIO $
             writeIORef envRef $
             envFrame {bindings = ((var, valueRef) : (bindings env))}
           return Void

bindVars :: Env -> [(String, LispVal)] -> IO Env
bindVars envRef bindings' = do
  extended <- extendEnv bindings' envRef
  newIORef extended
  where
    extendEnv :: [(String, LispVal)] -> Env -> IO EnvFrame
    extendEnv bindings'' env = do
      newBindings <- liftIO $ mapM addBinding bindings''
      return $ Frame {parent = env, bindings = newBindings}
    addBinding (var, value) = do
      ref <- newIORef value
      return (var, ref)

--------------
-- Apply
--------------
apply :: LispVal -> IOPrimitiveFunc
apply (PrimitiveFunc _ func) args = liftThrows $ func args
apply (Func _ params' varargs body' closure') args =
  if num params' /= num args && isNothing varargs
    then throwError $ NumArgs (Min $ length params') (length args) args
    else (liftIO $ bindVars closure' $ zip params' args) >>= bindVarArgs varargs >>=
         evalBody
  where
    remainingArgs = drop (length params') args
    num = toInteger . length
    evalBody env = fmap last $ mapM (eval env) body'
    bindVarArgs arg env =
      case arg of
        Just argName -> liftIO $ bindVars env [(argName, List $ remainingArgs)]
        Nothing -> return env
apply (IOFunc _ func) args = func args
apply f _ =
  throwError $
  Default $
  "application: not a procedure; " ++
  "expected a procedure that can be applied to arguments; given: " ++ show f

applyProc :: IOPrimitiveFunc
applyProc [func, List args] = apply func args
applyProc (func:args) = apply func args
applyProc _ = throwError $ Default "applyProc error"

readProc :: IOPrimitiveFunc
readProc [] = readProc [Port stdin]
readProc [Port port] = (liftIO $ hGetLine port) >>= liftThrows . readExpr
readProc _ = throwError $ Default "readProc error"

load :: String -> IOThrowsError [LispVal]
load filename = (liftIO $ readFile filename) >>= liftThrows . readExprList

readAll :: IOPrimitiveFunc
readAll [String filename] = fmap List $ load filename
readAll _ = throwError $ Default "readAll error"

--------------
-- Env
-- Bindings
-- Primitives
--------------
nullEnv :: IO Env
nullEnv = newIORef $ Global {bindings = []}

primitiveBindings :: IO Env -- TODO: Move. Here for now just to prevent circular dependency
primitiveBindings =
  nullEnv >>=
  (flip bindVars $
   map
     (makeFunc' IOFunc)
     (ioPrimitives ++
      [ ("apply", applyProc)
      , ("read", readProc)
      , ("read-all", readAll)
      , ("read", readProc)
      , ("read-all", readAll)
      ]) ++
   map (makeFunc' PrimitiveFunc) primitives)
  where
    makeFunc' ::
         (String -> ([LispVal] -> a) -> LispVal)
      -> (String, ([LispVal] -> a))
      -> (String, LispVal)
    makeFunc' constructor (var, func) = (var, constructor var func)

loadStdLib :: Env -> IO Env
loadStdLib env = do
  filename <- liftIO $ getDataFileName "lib/stdlib.scm"
  _ <- evalString env $ "(load \"" ++ filename ++ "\")"
  return env

--------------
-- Run
-- Repl
--------------
runIOThrows :: IOThrowsError String -> IO String
runIOThrows action = runExceptT (trapError action) >>= return . extractValue

readOrThrow :: Parser a -> String -> ThrowsError a
readOrThrow parser input =
  case parse parser input of
    Left (err, _) -> throwError $ DataTypes.Parser err
    Right [(val, _)] -> return val
    _ -> throwError $ Default "Read error"

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

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

evalAndPrint :: Env -> String -> IO ()
evalAndPrint env expr = evalString env expr >>= putStrLn

evalString :: Env -> String -> IO String
evalString env expr =
  runIOThrows $ fmap show $ (liftThrows $ readExpr expr) >>= eval env

runOne :: [String] -> IO ()
runOne args = do
  env <-
    primitiveBindings >>=
    flip bindVars [("args", List $ map String $ drop 1 args)]
  _ <- loadStdLib env
  result <-
    runIOThrows $
    fmap show' $ eval env (List [Atom "load-print", String (args !! 0)])
  putStrLn result
  where
    show' :: LispVal -> String
    show' (List contents) = unwordsList' contents
    show' _ = ""
    unwordsList' :: [LispVal] -> String
    unwordsList' = unlines . map showValNewline . filter printable
    showValNewline :: LispVal -> String
    showValNewline Void = ""
    showValNewline v = showVal v
    printable :: LispVal -> Bool
    printable Void = False
    printable _ = True

runRepl :: InputT IO ()
runRepl = do
  env <- lift $ primitiveBindings >>= loadStdLib
  repl env
  where
    repl :: Env -> InputT IO ()
    repl env = do
      minput <- getInputLine "Lisp>>> "
      case minput of
        Nothing -> return ()
        Just "quit" -> return ()
        Just input -> do
          lift $ evalAndPrint env input
          repl env

max_stars_repo_name: ebenpack/hascheme

max_issues_repo_name: ebenpack/hascheme

max_forks_repo_name: ebenpack/hascheme


text:
{-# LANGUAGE Haskell2010 #-}
import Control.Monad
import Data.Array
import Data.Bits
import Data.Char
import Data.Complex
import Data.Int
import Data.Ix
import Data.List
import Data.Maybe
import Data.Ratio
import Data.Word
import Foreign
import Foreign.C
import Foreign.C.Error
import Foreign.C.String
import Foreign.C.Types
import Foreign.ForeignPtr
import Foreign.Marshal
import Foreign.Marshal.Alloc
import Foreign.Marshal.Array
import Foreign.Marshal.Error
import Foreign.Marshal.Utils
import Foreign.Ptr
import Foreign.StablePtr
import Foreign.Storable
import Numeric
import Prelude
import System.Environment
import System.Exit
import System.IO
import System.IO.Error

main :: IO ()
main = return ()

max_stars_repo_name: jimcrayne/jhc

max_issues_repo_name: spl/jhc

max_forks_repo_name: spl/jhc


text:
module Language.Hakaru.Util.Coda where

import Statistics.Autocorrelation
import qualified Numeric.LinearAlgebra as V
import qualified Data.Vector.Generic as G

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

meanVariance :: Fractional a => [a] -> (a,a)
meanVariance lst = (av,sigma2)
  where
    n   = fromIntegral $ length lst
    av  = sum lst / n
    sigma2 = (foldr (\x acc -> sqr (x - av) + acc) 0 lst) / (n - 1)
    sqr x = x * x

max_stars_repo_name: SamuelSchlesinger/hakaru

max_issues_repo_name: SamuelSchlesinger/hakaru

max_forks_repo_name: SamuelSchlesinger/hakaru


text:
{-# LANGUAGE TemplateHaskell, RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables         #-}
{-# LANGUAGE MultiParamTypeClasses       #-}
{-# LANGUAGE TypeOperators               #-}
{-# LANGUAGE BangPatterns                #-}

module Chapter2Bandit
    ( Bandit(..)
    , Policy(..)
    , mkBandit
    , loopSteps
    ) where

import           Control.Monad
import           Control.Monad.IO.Class
import           Control.Monad.Trans.State

import           Control.Lens (makeLenses, over, element, (+~), (&), (.~), (%~))
import           Data.List(take, repeat)

import           Data.Random
import           Data.Random.Distribution
import           Data.Random.Distribution.Uniform
import           Data.Random.Shuffle.Weighted
import           Data.Random.Source.IO

import           Language.Haskell.TH
import           Numeric.LinearAlgebra (Vector, Matrix)
import qualified Numeric.LinearAlgebra as LA

-- project
import           Utils

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

-- | Bandit Problem
--   Three methods (Policies): epsilon-greedy, UCB, gradient bandit.
--   Policy includes two parts: action selection & state-value estimation update.

-- | For action selections (exploit & explore) of above methods:
--   1. epsilon-greedy: randomly do exploration with probability epsilon;
--   2. UCB using: Qt(a) + c sqrt(log t / Nt(a)) to balance exploit and explore
--   3. Gradient using Qt(a) as Preference value, select action by softmax;

-- | For Qt(a) (Value Estimate) Update:
--   1. 'epsilon-greedy' & 'UCB' methods using form of:
--         newVal = oldVal + stepSize * (curReward - oldVal)
--      here, 'stepSize' parameter in configuration file:
--        if < 0, use 'sample average', then it is stationary;
--        if > 0, use 'exponential, recency-weighted average',
--                then it targets non-stationary problem, weight recent reward more.
--   2. Gradient bandit using gradient approximation, estimate value as preference:
--           newValue = oldValue + stepSize * (curReward - baselineVal) * ( 1 - softmax)
--           for the selected action's state value update
--       and newValue = oldValue - stepSize * (curReward - baselineVal) * softmax
--           for all other non-selected actions'
--       where 'baseline' could be 0, or the average value of all action rewards until the time
--       and 'stepSize' is the same as above.

------------------------------------------------------------------------------------------
data Policy = EGreedy (RVar Bool) -- Bernoulli distribution with p = epsilon
            | UCB Double -- ucb explore parameter 'c'
            | Gradient Bool -- gradient baseline param, if < 0, use average reward value

data Bandit = Bandit {
     _kArms :: Int
    ,_bestValueIdx :: Int -- it is pre-known, for OptimalAction statistics
    ,_qValues :: ![Double] -- estimate value of each action
    ,_nActions :: ![Int] -- count of each atcion has taken
    ,_stepSize :: Double
    ,_totalReward :: !Double
    ,_curStep :: !Int
    ,_bestTakes :: ![Double]
    -- input randoms
    ,_srcRVars :: ![RVar Double]
    -- exploit & explore methods
    ,_policy :: Policy
    }

makeLenses ''Bandit

mkBandit :: Int -> Int -> Double -> Double -> [Double] -> Policy -> Bandit
mkBandit karm totalStep initValue stepSize trueValues policy =
  let (maxValueIdx, _) = maxWithIndex (zip [0..] trueValues)
  in  (Bandit karm maxValueIdx (take karm $ repeat initValue)
              (take karm $ repeat 0) stepSize 0.0 0 []
              (map (flip normal 1.0) trueValues) policy
      )

--------------------------------------------------------------------------------
---- learning

loopSteps :: Int -> StateT Bandit IO [Double]
loopSteps times = replicateM times step

step :: StateT Bandit IO Double
step = selectOneAction >>= takeOneAction

selectOneAction :: StateT Bandit IO Int
selectOneAction = do
  bandit <- get
  actionN <- case _policy bandit of
    EGreedy epsilonRVar -> do
      bExplore <- liftIO $ sample epsilonRVar
      if bExplore then liftIO $ sample (randomElement [0..((_kArms bandit) - 1)])
         else pure . fst . maxWithIndex $ (zip [0..] (_qValues bandit))
    UCB c -> do
      let !ucbEstValues = zipWith (calcUCB c (_curStep bandit)) (_nActions bandit) (_qValues bandit)
      pure $ fst $ maxWithIndex (zip [0..] ucbEstValues)
    Gradient _ -> do
      let gradientValues = map exp (_qValues bandit)
          gradientProbs = map ( / sum gradientValues) gradientValues
          -- Sampling according to distribution (element probability)
          -- The following line will produce good result, even 0.4 without baseline. why?
          -- Is it because we do weighted random extract twice ?
          weightedRVar = weightedShuffle $ zip gradientProbs [0..]
      liftIO $ head <$> sample weightedRVar
      -- pure $ fst $ maxWithIndex (zip [0..] gradientProbs)
  pure actionN
  where
  calcUCB :: Double -> Int -> Int -> Double -> Double
  calcUCB c total n val = val + c * sqrt ((log $ fromIntegral total) / fromIntegral n)

takeOneAction :: Int -> StateT Bandit IO Double
takeOneAction actionN = do
  bandit <- get
  reward <- liftIO $ sample (_srcRVars bandit !! actionN)
  let !bestTake = (_bestValueIdx bandit == actionN) ? (1.0, 0.0)
      !bandit' = bandit & (nActions . element actionN +~ 1)
                        & (totalReward +~ reward)
                        & (curStep +~ 1)
                        & (bestTakes %~ (++ [bestTake]))
      !bandit'' = stateUpdate bandit' actionN reward
  put bandit''
  pure reward
  where
  stateUpdate :: Bandit -> Int -> Double -> Bandit
  stateUpdate bandit actionN reward =
    let !actionTakes = (_nActions bandit) !! actionN
        !ss = (_stepSize bandit < 0) ? (1.0 / fromIntegral actionTakes, _stepSize bandit)
        !oldEstValue = (_qValues bandit) !! actionN
    in  case _policy bandit of
          Gradient bBaseline ->
            let baseline = bBaseline ? (_totalReward bandit / (fromIntegral $ _curStep bandit), 0.0)
                !gradientValues = map exp (_qValues bandit)
                !gradientProbs = map ( / sum gradientValues) gradientValues
                !newEstValues = zipWith3 (updateGradientPreference reward baseline ss)
                                         (_qValues bandit) gradientProbs [0..]
            in  bandit & (qValues .~ newEstValues)
          -- epsilone greedy & ucb use the same updating formular
          _ -> let !newEstValue = oldEstValue + (reward - oldEstValue) * ss
               in  bandit & (qValues . element actionN .~ newEstValue)

  updateGradientPreference :: Double -> Double -> Double -> Double -> Double -> Int -> Double
  updateGradientPreference reward baseline ss oldVal prob actionIdx =
    (actionN == actionIdx) ? (oldVal + ss * (reward - baseline) * (1 - prob),
                              oldVal - ss * (reward - baseline) * prob)
--------------------------------------------------------------------------------

max_stars_repo_name: Xingtao/ReinforceLearningIntro

max_issues_repo_name: Xingtao/ReinforceLearningIntro

max_forks_repo_name: Xingtao/ReinforceLearningIntro


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE LambdaCase #-}
module Taiji.Utils.Clustering
    ( evalClusters
    , computeClusterMetrics
    , computeReproducibility
    , readKNNGraph
    , optimalParam
    , visualizeCluster
    , sampleCells
    , silhouette
    , confusionTable
    , batchCorrect
    , leiden
    , ari
    ) where

import qualified Data.ByteString.Char8 as B
import Data.Binary (encodeFile, decodeFile)
import Bio.Utils.Misc
import qualified Data.Text as T
import qualified Data.HashSet as S
import qualified Data.HashMap.Strict as M
import Data.Hashable (Hashable)
import Control.Arrow (first, second)
import Conduit
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as VM
import qualified Data.Vector.Unboxed as U
import Data.List.Ordered (nubSort)
import Shelly (shelly, run_)
import Data.Conduit.Zlib (multiple, ungzip)
import AI.Clustering.Hierarchical.Types (computeDists, (!))
import AI.Clustering.Hierarchical (euclidean)
import Statistics.Sample (mean)
import System.Random.MWC.Distributions
import Language.Javascript.JMacro
import System.Random.MWC

import IGraph.Random
import IGraph
import qualified IGraph.Algorithms.Community as I

import Taiji.Prelude
import Taiji.Utils.Plot.ECharts
import Taiji.Utils.Plot

-- | Evaluating the clustering results
evalClusters :: FilePath
             -> Optimizer  -- ^ Optimizer
             -> Double     -- ^ Resolution
             -> FilePath   -- ^ knn
             -> IO ([FilePath], FilePath)
evalClusters dir optimizer res knn = do
    gr <- readKNNGraph knn
    gen <- create

    perturbed <- forM [1::Int ..5] $ \i -> do
        let output' = dir <> "/perturbed_clustering_result_" <> show i <> ".bin"
        r <- mutateGraph gr gen >>= leiden res optimizer
        encodeFile output' $ map snd $ sort $ concat $ zipWith (\a b -> zip a $ repeat b) r [0::Int ..]
        return output'

    let output = dir <> "/clustering_result.bin"
    (filter ((>100) . length) <$> leiden res optimizer gr) >>= encodeFile output
    return (perturbed, output)
  where
    mutateGraph gr gen = do
        let n = truncate $ fromIntegral (nEdges gr) * (0.02 :: Double)
        xs <- take n . V.toList <$> uniformShuffle (V.fromList $ edges gr) gen
        return $ delEdges xs gr

computeReproducibility :: [[Int]] -- ^ cluster
                       -> [[Int]]  -- ^ perturbed clusters
                       -> [Double]
computeReproducibility a bs = map (mean . U.fromList) $ transpose $ flip map bs $ \b ->
    let b' = U.fromList b
    in map (reproducibility . map (b' U.!)) a

computeClusterMetrics :: ([FilePath], FilePath)
                      -> FilePath
                      -> IO (Int, Double, Double)
computeClusterMetrics (perturbed, cl) coordinate = do
    stability <- (mean . U.fromList . map ari . comb) <$> mapM decodeFile perturbed
    clusters <- decodeFile cl
    points <- fmap V.fromList $ runResourceT $ runConduit $ sourceFile coordinate .|
        multiple ungzip .| linesUnboundedAsciiC .|
        mapC (U.fromList . map readDouble . B.split '\t') .| sinkList
    gen <- create

    let samplingRatio = 10000 / fromIntegral (V.length points) :: Double
    clusters' <- forM clusters $ \c -> do
        let n = max 2 $ truncate $ fromIntegral (length c) * samplingRatio
        V.toList . V.take n <$> uniformShuffle (V.fromList c) gen

    let sil = silhouette $ (map . map) (points V.!) clusters'
    return (length clusters, sil, stability)
  where
    comb (x:xs) = zip (repeat x) xs ++ comb xs
    comb _ = []

readKNNGraph :: FilePath -> IO (Graph 'U () Double)
readKNNGraph fl = runResourceT $ runConduit $ sourceFile fl .|
    multiple ungzip .| linesUnboundedAsciiC .| sink
  where
    sink = do
        n <- headC >>= \case
            Nothing -> error ""
            Just x -> return $ readInt x
        es <- mapC f .| sinkList
        return $ mkGraph (replicate n ()) $ replaceInf es
    f x = let (a:b:c:_) = B.split '\t' x
          in ((readInt a, readInt b), readDouble' c)
    readDouble' x | x == "inf" = Nothing
                  | otherwise = Just $ readDouble x
    replaceInf es = let m = maximum $ mapMaybe snd es
                    in map (second (fromMaybe m)) es

optimalParam :: FilePath
             -> [(Double, (Int, Double, Double))]
             -> IO Double
optimalParam output input = do
    savePlots output [] plt
    return optimal
  where
    optimal = fst $ maximumBy (comparing (^._2._2)) $ case filter (\x -> x^._2._3 >= 0.9) input of
        [] -> input
        x -> x
    (res, dat) = unzip $ flip map (sortBy (comparing fst) input) $ \(r, (n, sil, stab)) ->
        (r, (fromIntegral n, sil, stab))
    (num, sils, stabs) = unzip3 dat
    plt = map (setDim 400 300 . addAttr toolbox)
        [ addAttr (yAxisLabel "number of clusters") $
                addAttr (xAxisLabel "resolution") $
                lineplot $ zip res num
        , addAttr (yAxisLabel "silhouette width") $
                addAttr (xAxisLabel "resolution") $
                lineplot $ zip res sils
        , addAttr (yAxisLabel "stability") $
                addAttr (xAxisLabel "resolution") $
                lineplot $ zip res stabs ]
    lineplot input = mkEChart [jmacroE| {
        grid: { containLabel: true },
        xAxis: {
            type: "value",
            axisLine: {onZero: false}
        },
        tooltip: {
            trigger: "axis"
        },
        yAxis: {
            type: "value",
            axisLine: {onZero: false}
        },
        series: [{
            data: `dataset`,
            markLine: {symbolSize: 0, precision: 10, data: [{xAxis: `optimal`}]},
            type: "line"
        }]
        } |]
      where
        dataset = map (\(x,y) -> [x,y]) input

visualizeCluster :: [CellCluster] -> [EChart]
visualizeCluster cls =
    [ addAttr toolbox $ scatter' $ flip map cls $ \(CellCluster nm cells _) ->
        (B.unpack nm, map _cell_2d cells) ] ++ if noName then [] else
          [ addAttr toolbox $ scatter' $ map (first head . unzip) $
              groupBy ((==) `on` fst) $ sortBy (comparing fst) $ concatMap
              (map (\x -> (getName $ _cell_barcode x, _cell_2d x)) . _cluster_member) cls
          ]
  where
    noName = null $ getName $ _cell_barcode $ head $
        _cluster_member $ head cls
    getName x = let prefix = fst $ B.breakEnd (=='+') x
                in if B.null prefix then "" else B.unpack $ B.init prefix

-- | Random sample 30,000 cells.
sampleCells :: [CellCluster] -> IO [CellCluster]
sampleCells clusters
    | ratio >= 1 = return clusters
    | otherwise = do
        gen <- create
        forM clusters $ \c -> do
            s <- sampling gen ratio $ V.fromList $ _cluster_member c
            return $ c {_cluster_member = V.toList s}
  where
    n = foldl1' (+) $ map (length . _cluster_member) clusters
    ratio = 1 / (fromIntegral n / 30000) :: Double
    sampling gen frac v = V.take n' <$> uniformShuffle v gen
      where
        n' = max 200 $ truncate $ frac * fromIntegral (V.length v)

silhouette :: [[U.Vector Double]]   -- ^ Groups of objects
           -> Double
silhouette input
    | length input <= 1 = 0
    | otherwise = mean $ U.fromList $ flip concatMap clusters $ \(i, cl) ->
        let cls = map snd $ filter ((/= i) . fst) clusters
        in flip map cl $ \x ->
            let a = getD x cl
                b = minimum $ map (getD x) cls
            in (b - a) / max a b
  where
    getD i cl = mean $ U.fromList $ map (\j -> distMat ! (i, j)) $ filter (/= i) cl
    clusters = zip [0 :: Int ..] $ go 0 input
      where
        go i (x:xs) = let n = length x in take n [i..] : go (i + n) xs
        go _ _ = []
    distMat = computeDists euclidean $ V.fromList $ concat input
{-# INLINE silhouette #-}

confusionTable :: (Eq a, Hashable a) => [[a]] -> [[a]] -> [[Int]]
confusionTable inputA inputB = flip map inputA' $ \a -> flip map inputB' $ \b ->
    S.size $ S.intersection a b
  where
    inputA' = map S.fromList inputA
    inputB' = map S.fromList inputB
{-# INLINE confusionTable #-}

batchCorrect :: [Maybe (B.ByteString, Maybe B.ByteString)]  -- ^ label and group
             -> V.Vector (U.Vector Double)    -- ^ Data
             -> IO (V.Vector (U.Vector Double))
batchCorrect batchInfo input = do
    vec <- V.thaw input
    forM_ batchGroups $ \(idx, label) -> when (length (nubSort label) > 1) $ do
        res <- batchAveraging label (map (input V.!) idx)
        forM_ (zip idx res) $ \(i, r) -> VM.unsafeWrite vec i r
    V.unsafeFreeze vec
  where
    batchGroups = map (unzip . map (\(i, (x, _)) -> (i, x))) $
        groupBy ((==) `on` (snd . snd)) $ sortBy (comparing (snd . snd)) $
        mapMaybe (\(i, b) -> maybe Nothing (\x -> Just (i,x)) b) $ zip [0..] batchInfo
{-# INLINE batchCorrect #-}

batchAveraging :: [B.ByteString]    -- ^ Batch labels
               -> [U.Vector Double]    -- ^ Data
               -> IO [U.Vector Double]
batchAveraging labels input = withTempDir Nothing $ \dir -> do
    let labelFl = dir <> "/label.txt"
        matFl = dir <> "/mat.txt"
        tmp = dir <> "/tmp.txt"
    B.writeFile labelFl $ B.unlines labels
    B.writeFile matFl $ B.unlines $ map (B.intercalate "\t" . map toShortest . U.toList) input
    shelly $ run_ "taiji-utils" ["correct", "--label", T.pack labelFl, T.pack matFl, T.pack tmp]
    map (U.fromList . map readDouble . B.words) . B.lines <$> B.readFile tmp
{-# INLINE batchAveraging #-}

leiden :: Double -> Optimizer -> Graph 'U () Double -> IO [[Int]]
leiden resolution optimizer gr = withSeed 9304 $ fmap (sortBy (flip (comparing length))) .
    I.findCommunity gr nodeWeight edgeWeight I.leiden{I._resolution=res}
  where
    edgeWeight
        | optimizer `elem` [RBConfiguration, CPM] = Nothing
        | otherwise = Just id
    (nodeWeight, res)
        | optimizer `elem` [RBConfiguration, RBConfigurationWeighted] =
            ( Just $ \i _ -> foldl1' (+) $ map (\j -> edgeLab gr (i, j)) $ neighbors gr i
            , resolution / (2 * (foldl1' (+) $ map snd $ labEdges gr)))
        | otherwise = (Nothing, resolution)

-- |
reproducibility :: [Int]   -- membership
                -> Double
reproducibility xs = fromIntegral sameGroup / fromIntegral (c2 $ length xs)
  where
    counts = M.elems $ M.fromListWith (+) $ zip xs $ repeat (1 :: Int)
    sameGroup = foldl1' (+) $ map c2 counts
    c2 x = (x * (x - 1)) `div` 2

-- | Adjusted Rand Index: <http://en.wikipedia.org/wiki/Rand_index>
ari :: ([Int], [Int]) -> Double
ari (a, b) | length a /= length b = error "unequal length"
           | a == b = 1
           | otherwise = ari' $ counts a b
  where
    ari' (Counts cxy cx cy) =  (sum1 - sum2*sum3/choicen2)
                            / (1/2 * (sum2+sum3) - (sum2*sum3) / choicen2)
      where choicen2 = choice (sum . M.elems $ cx) 2
            sum1 = sum [ choice nij 2 | nij <- M.elems cxy ]
            sum2 = sum [ choice ni 2 | ni <- M.elems cx ]
            sum3 = sum [ choice nj 2 | nj <- M.elems cy ]
    -- | Creates count table 'Counts'
    counts xs = foldl' f mempty . zipWith ((,)) xs
        where f cs@(Counts cxy cx cy) p@(x,y) =
                cs { joint       = M.insertWith (+) p 1 cxy
                , marginalFst = M.insertWith (+) x 1 cx
                , marginalSnd = M.insertWith (+) y 1 cy }
    -- | The binomial coefficient: C^n_k = PROD^k_i=1 (n-k-i)\/i
    choice :: Double -> Double -> Double
    choice n k = foldl' (*) 1 [n-k+1 .. n] / foldl' (*) 1 [1 .. k]

-- | Count table
data Counts = Counts
    { joint :: !(M.HashMap (Int, Int) Double) -- ^ Counts of both components
    , marginalFst :: !(M.HashMap Int Double) -- ^ Counts of the first component
    , marginalSnd :: !(M.HashMap Int Double) -- ^ Counts of the second component
    }

instance Monoid Counts where
    mempty = Counts M.empty M.empty M.empty
    c `mappend` k =
        Counts { joint = unionPlus (joint c) (joint k)
               , marginalFst = unionPlus (marginalFst c) (marginalFst k)
               , marginalSnd = unionPlus (marginalSnd c) (marginalSnd k)
               }
      where
        unionPlus m = M.foldlWithKey' (\z k' v -> M.insertWith (+) k' v z) m

instance Semigroup Counts where
    c <> k =
        Counts { joint = unionPlus (joint c) (joint k)
               , marginalFst = unionPlus (marginalFst c) (marginalFst k)
               , marginalSnd = unionPlus (marginalSnd c) (marginalSnd k)
               }
      where
        unionPlus m = M.foldlWithKey' (\z k' v -> M.insertWith (+) k' v z) m

max_stars_repo_name: Taiji-pipeline/Taiji-utils

max_issues_repo_name: Taiji-pipeline/Taiji-utils

max_forks_repo_name: Taiji-pipeline/Taiji-utils


text:
-- Filename: fft.hs
-- Created by: Daniel Winograd-Cort
-- Created on: unknown
-- Last Modified by: Daniel Winograd-Cort
-- Last Modified on: 12/12/2013

-- -- DESCRIPTION --
-- This code was inspired by Euterpea.  It uses UISF to present a GUI that
-- shows the sum of two waves (whose frequencies are specified by the user)
-- as well as the Fast Fourier Transform of that sum.
--
-- This module requires the array and pure-fft packages.

{-# LANGUAGE Arrows #-}
module FRP.UISF.Examples.FFT where
import FRP.UISF
import Numeric.FFT (fft)
import Data.Complex
import Data.Map (Map)
import Data.Maybe (listToMaybe, catMaybes)
import qualified Data.Map as Map

import Data.Array.Unboxed

-- | Alternative for working with Math.FFT instead of Numeric.FFT
--import qualified Math.FFT as FFT
--import Data.Array.IArray
--import Data.Array.CArray
--myFFT n lst = elems $ (FFT.dft) (listArray (0, n-1) lst)

--------------------------------------
-- Sine wave oscillators
--------------------------------------

-- Table definition
type Table = UArray Int Double

-- Sine table generator. Takes an integer representing the number of samples to generate
-- and a list of relative intensities for the overtones of the sine wave.

tableSinesN :: Int -> [Double] -> Table
tableSinesN size amps =
    let wave x   = sum (zipWith (*) [sin (2*pi*x*n) | n <- [1..]] amps)
        delta    = 1 / fromIntegral size
        waveform = take size $ map wave [0,delta..]
        divisor  = (maximum . map abs) waveform
     in listArray (0,size) (map (/divisor) waveform)

-- Two example sine tables.

tab1, tab2 :: Table
tab1 = tableSinesN 4096 [1] -- Basic sine wave
tab2 = tableSinesN 4096 [1.0,0.5,0.33]

-- Table-driven oscillator

osc :: ArrowCircuit a => Table -> Double -> a Double Double
osc table sr = proc freq -> do
    rec
      let delta = 1 / sr * freq
          phase = if next > 1 then frac next else next
      next <- delay 0 -< frac (phase + delta)
    returnA -< ((table!).(`mod` size).round.(*rate)) phase
  where (_,size)    = bounds table
        rate        = fromIntegral size
        frac x = if x > 1 then x - fromIntegral (truncate x) else x

--------------------------------------
-- Fast Fourier Transform
--------------------------------------

-- | Returns n samples of type b from the input stream at a time,
--   updating after k samples.  This function is good for chunking
--   data and is a critical component to fftA
quantize :: ArrowCircuit a => Int -> Int -> a b (SEvent [b])
quantize n k = proc d -> do
    rec (ds,c) <- delay ([],0) -< (take n (d:ds), c+1)
    returnA -< if c >= n && c `mod` k == 0 then Just ds else Nothing

-- | Converts the vector result of a dft into a map from frequency to magnitude.
--   One common use is:
--      fftA >>> arr (fmap $ presentFFT clockRate)
presentFFT :: Double -> [Double] -> Map Double Double
presentFFT clockRate a = Map.fromList $ zipWith (curry mkAssoc) [0..] a where
    mkAssoc (i,c) = (freq, c) where
        samplesPerPeriod = fromIntegral (length a)
        freq = fromIntegral i * (clockRate / samplesPerPeriod)

-- | Given a quantization frequency (the number of samples between each
--   successive FFT calculation) and a fundamental period, this will decompose
--   the input signal into its constituent frequencies.
--   NOTE: The fundamental period must be a power of two!
fftA :: ArrowCircuit a => Int -> Int -> a Double (SEvent [Double])
fftA qf fp = proc d -> do
    carray <- quantize fp qf -< d :+ 0
    returnA -< fmap (map magnitude . take (fp `div` 2) . fft) carray

--------------------------------------
-- UISF Example
--------------------------------------

-- This example shows off the histogram and realtimeGraph widgets by
-- summing two sin waves and displaying them.  Additionally, it makes
-- use of two horizontal sliders.
-- This example also shows off asyncUISFV and how to take a SigFun,
-- of the type used to create sound, and convert it to a UISF.
fftEx :: UISF () ()
fftEx = proc _ -> do
    f1 <- hSlider (1, 2000) 440 -< ()
    _ <- leftRight (label "Freq 1: " >>> display) -< f1
    f2 <- hSlider (1, 2000) 440 -< ()
    _ <- leftRight (label "Freq 2: " >>> display) -< f2
    d <- asyncVT sr 0.1 myAutomaton -< (f1, f2)
    let fft = listToMaybe $ catMaybes $ map (snd . fst) d
        s = map (\((s, _), t) -> (s,t)) d
    _ <- histogram (makeLayout (Stretchy 10) (Fixed 150)) -< fft
    _ <- realtimeGraph (makeLayout (Stretchy 10) (Fixed 150)) 2 Black -< s
    returnA -< ()
  where
    sr = 1000 -- signal rate
    myAutomaton = proc (f1, f2) -> do
        s1 <- osc tab1 sr -< f1
        s2 <- osc tab2 sr -< f2
        let s = (s1 + s2)/2
        fftData <- fftA 100 256 -< s
        returnA -< (s, fftData)

-- This test is run separately from the others.
main :: IO ()
main = runUI (defaultUIParams {uiSize=(500, 600), uiTitle="FFT Example"}) fftEx

max_stars_repo_name: dwincort/UISF

max_issues_repo_name: dwincort/UISF

max_forks_repo_name: dwincort/UISF


text:
module Numeric.FFT.Plan ( plan, planFromFactors ) where

import qualified Control.Monad          as CM
import           Control.Monad.IO.Class
import           Criterion.Measurement              (measure)
import           Criterion.Measurement.Types hiding (measure)
import           Data.Complex
import           Data.Function          (on)
import qualified Data.IntMap.Strict     as IM
import           Data.List              ((\\))
import qualified Data.List              as L
import qualified Data.Set               as S
import qualified Data.Vector            as V
import           Data.Vector.Unboxed
import           Prelude                hiding (any, concatMap, enumFromTo,
                                         filter, length, map, maximum, null,
                                         reverse, scanl, sum, zip, zipWith,
                                         (++))
import qualified Prelude                as P
import           System.Directory
import           System.Environment
import           System.FilePath

import           Numeric.FFT.Execute
import           Numeric.FFT.Special
import           Numeric.FFT.Types
import           Numeric.FFT.Utils

-- | Number of plans to test empirically.
nTestPlans :: Int
nTestPlans = 50

-- | Plan calculation for a given problem size.
plan :: Int -> IO Plan
plan 1 = return $ Plan V.empty Nothing (SpecialBase 1)
plan n = do
  wis <- readWisdom n
  let fixRader p = case plBase p of
        bpl@(RaderBase _ _ _ _ csz _) -> do
          cplan <- liftIO $ plan csz
          return $ p { plBase = bpl { raderConvPlan = cplan } }
        _ -> return p
  pret <- case wis of
    Just (p, _) -> planFromFactors n p
    Nothing -> do
      let ps = testPlans n nTestPlans
          v = generate n (\i -> sin (2 * pi * fromIntegral i / 511) :+ 0)
      tps <- CM.forM ps $ \p -> do
        ptest <- liftIO $ planFromFactors n p >>= fixRader
        let niters = fromIntegral $ 50000 `div` n
        (meas, _) <- measure (nf (execute ptest Forward) v) niters
        let t = measTime meas / fromIntegral (measIters meas)
            _iters = measIters meas
        return (t, p)
      let (rest, resp) = L.minimumBy (compare `on` fst) tps
      liftIO $ writeWisdom n resp rest
      liftIO $ planFromFactors n resp
  fixRader pret

-- | Get execution time for plan for a given problem size.
planTime :: Int -> IO Double
planTime n = do
  wis <- readWisdom n
  case wis of
    Just (_, t) -> return t
    Nothing -> do
      -- Force generation of wisdom if it's not already there.
      _ <- plan n
      Just (_, t) <- readWisdom n
      return t

-- | Plan calculation for a given problem factorisation.
planFromFactors :: Int -> (Int, Vector Int) -> IO Plan
planFromFactors n (lastf, fs) = do
  -- Base transform.
  (base, mextraperm) <- makeBase lastf

  -- Include permutation of base transform if needed.
  let perm = case (digperm, mextraperm) of
        (Just dp, Just ep) -> Just $ dupperm n ep %.% dp
        (Nothing, Just ep) -> Just $ dupperm n ep
        (Just dp, Nothing) -> Just dp
        (Nothing, Nothing) -> Nothing

  return $ Plan dlinfo perm base
  where
    -- Input data "digit reversal" permutation.
    digperm = digrev n fs

    -- Size information for Danielson-Lanczos steps.
    wfacs = map (n `div`) $ scanl (*) 1 fs
    vwfacs = convert wfacs
    vfs = convert fs
    dmatps = V.zipWith (dmat 1) vwfacs vfs
    dmatms = V.zipWith (dmat (-1)) vwfacs vfs
    dlinfo = V.reverse $ V.zip4 vwfacs vfs dmatps dmatms

    -- Calculate diagonal matrix entries used in Danielson-Lanczos steps.
    dmat :: Int -> Int -> Int -> VVVCD
    dmat sign wfac split =
      let ns = wfac `div` split
          w = omega $ sign * wfac
      in V.generate split $
         \r -> V.generate split $
           \c -> map ((w^(ns*r*c) *) . ((w^^) . (c *))) $ enumFromN 0 ns

-- | Read from wisdom for a given problem size.
readWisdom :: Int -> IO (Maybe ((Int, Vector Int), Double))
readWisdom n = do
  home <- getEnv "HOME"
  let wisf = home </> ".fft-plan" </> show n
  ex <- doesFileExist wisf
  if ex then
    (do wist <- readFile wisf
        let ((wisb, wisfs), wistim) = read wist :: ((Int, [Int]), Double)
        return $ Just ((wisb, fromList wisfs), wistim))
    else return Nothing

-- | Write wisdom for a given problem size.
writeWisdom :: Int -> (Int, Vector Int) -> Double -> IO ()
writeWisdom n (b, fs) tim = do
  home <- getEnv "HOME"
  let wisd = home </> ".fft-plan"
      wisf = wisd </> show n
  createDirectoryIfMissing True wisd
  writeFile wisf $ show ((b, toList fs), tim) P.++ "\n"

-- | Make base transform for a given sub-problem size.
makeBase :: Int -> IO (BaseTransform, Maybe VI)
makeBase sz
  | sz `IM.member` specialBases = return (SpecialBase sz, Nothing)
  | isPrime sz                  = makeRaderBase sz
  | otherwise                   = return (makeDFTBase sz, Nothing)

-- | Generate digit reversal permutation using elementary "modulo"
-- permutations: last digit is not permuted to match with using a
-- simple DFT or instance of Rader's algorithm at the "bottom" of the
-- overall algorithm.
digrev :: Int -> VI -> Maybe VI
digrev n fs
  | null fs   = Nothing
  | otherwise = Just $ V.foldl1' (%.%) $ V.map (dupperm n) subperms
  where
    vfs = convert fs

    -- Sizes of the individual permutations that we need, one per
    -- factor.
    sizes = V.scanl div n vfs

    -- Partial sub-permutations, one per factor.
    subperms = V.reverse $ V.zipWith perm sizes vfs

    -- Generate a single "modulo" permutation.
    perm sz fac = concatMap doone $ enumFromN 0 fac
      where doone i = generate (sz `div` fac) (\j -> j * fac + i)

-- | Pre-computation plan for basic DFT transform.
makeDFTBase :: Int -> BaseTransform
makeDFTBase sz = DFTBase sz wsfwd wsinv
  where w = omega sz
        wsfwd = generate sz (w ^)
        wsinv = map (1 /) wsfwd

-- | Pre-compute plan for prime-length Rader FFT transform.
makeRaderBase :: Int -> IO (BaseTransform, Maybe VI)
makeRaderBase sz = do
  -- Forward FFT with embedded plan calculation.
  let fft p xs = execute p Forward xs

  -- Times for unpadded and padded convolution transforms.
  unpadtime <- planTime sz1
  padtime <- planTime pow2sz

  -- Should we use a padded or an unpadded transform for the
  -- convolution?
  let pad = padtime < unpadtime
      csz = if pad then pow2sz else sz1

  -- Plan for convolution transforms.
  cplan <- plan csz

  -- FFT transforms of b sequences for use in convolution, one for
  -- forward transform and one for inverse transform.  Either just the
  -- basic root of unity powers for the convolution (if we're not
  -- padding), or the powers cyclically repeated to make a vector of
  -- the next power-of-two length.
  let convb =
        fft cplan $ if pad
                    then generate csz (\idx -> bs ! (idx `mod` sz1))
                    else bs
      convbinv =
        fft cplan $ if pad
                    then generate csz (\idx -> bsinv ! (idx `mod` sz1))
                    else bsinv

  return (RaderBase sz outperm convb convbinv csz cplan, Just inperm)
  where
    -- Convolution length.
    sz1 = sz - 1

    -- Convolution length padded to next greater power of two.
    pow2sz = if sz1 == 2^ log2 sz1
             then sz1
             else 2 ^ (1 + log2 (2 * sz1 - 3))

    -- Group generator and inverse group generator.
    g = primitiveRoot sz
    ig = invModN sz g

    -- Input value permutation according to group generator indexing.
    inperm = 0 `cons` iterateN sz1 (\n -> (g * n) `mod` sz) 1

    -- Index vector based on inverse group generator ordering.
    outperm = iterateN sz1 (\n -> (ig * n) `mod` sz) 1

    -- Root of unity powers based on inverse group generator indexing,
    -- for forward and inverse transform.
    w = omega sz
    bs = backpermute (map (w ^^) $ enumFromTo 0 sz1) outperm
    bsinv = backpermute (map ((w ^^) . negate) $ enumFromTo 0 sz1) outperm

-- | Base transform type with heuristic ordering.
data BaseType = Special Int | Rader Int deriving (Eq, Show)

-- | Newtype wrapper for custom sorting.
newtype SPlan = SPlan (BaseType, Vector Int) deriving (Eq, Show)

-- | Base transform size.
bSize :: BaseType -> Int
bSize (Special b) = b
bSize (Rader b)   = b

-- | Heuristic ordering for base transform types: special bases come
-- first, then prime bases using Rader's algorithm, ordered according
-- to size compensating for padding needed in the Rader's algorithm
-- convolution.
instance Ord BaseType where
  compare (Special _)  (Rader _)    = LT
  compare (Rader _)    (Special _)  = GT
  compare (Special s1) (Special s2) = compare s2 s1
  compare (Rader r1)   (Rader r2)   = case (isPow2 $ r1 - 1, isPow2 $ r2 - 1) of
    (True, True)   -> compare r1 r2
    (True, False)  -> compare r1 (2 * r2)
    (False, True)  -> compare (2 * r1) r2
    (False, False) -> compare r1 r2

-- | Heuristic ordering for full plans, based first on base type, then
-- on the maximum size of Danielson-Lanczos step.
instance Ord SPlan where
  compare (SPlan (b1, fs1)) (SPlan (b2, fs2)) = case compare b1 b2 of
    LT -> LT
    EQ -> compare (maximum fs2) (maximum fs1)
    GT -> GT

-- | Generate test plans for a given input size, sorted in heuristic
-- order.
testPlans :: Int -> Int -> [(Int, Vector Int)]
testPlans n nplans = L.take nplans $ L.map clean $ L.sort okplans
  where vfs = allFactors n
        bs = usableBases n vfs
        doone b = basePlans n vfs b
        clean (SPlan (b, fs)) = (bSize b, fs)
        allplans = P.concatMap doone bs
        okplans = case L.filter (not . ridiculous) allplans of
          []  -> L.filter (not . reallyRidiculous) allplans
          oks -> oks
        ridiculous (SPlan (_, fs)) = any (> 128) fs
        reallyRidiculous (SPlan (_, fs)) =
          any (> 128) $ filter (not . isPrime) fs

-- | List plans from a single base.
basePlans :: Int -> Vector Int -> BaseType -> [SPlan]
basePlans _n vfs bt = if null lfs
                     then [SPlan (bt, empty)]
                     else P.map (\v -> SPlan (bt, v)) $ leftOvers lfs
  where lfs = fromList $ toList vfs \\ toList (allFactors b)
        b = bSize bt

-- | Produce all distinct permutations and compositions constructable
-- from a given list of factors.
leftOvers :: Vector Int -> [Vector Int]
leftOvers fs =
  if null fs
  then []
  else S.toList $ L.foldl' go S.empty (multisetPerms fs)
  where n = length fs
        go fset perm = foldl' doone fset (enumFromN 0 (2^(n - 1)))
          where doone s i = S.insert (makeComp perm i) s

-- | Usable base transform sizes.
usableBases :: Int -> Vector Int -> [BaseType]
usableBases n fs = P.map Special bs P.++ P.map Rader ps
  where bs = toList $ filter ((== 0) . (n `mod`)) specialBaseSizes
        ps = toList $ filter isPrime $ filter (> maxPrimeSpecialBaseSize) fs

max_stars_repo_name: ian-ross/arb-fft

max_issues_repo_name: ian-ross/arb-fft

max_forks_repo_name: ian-ross/arb-fft


text:
module Feldspar
  ( module Prelude.EDSL
  , module Control.Monad
    -- * Types
    -- ** Syntax
  , Data
  , Syntax
  , Comp
    -- ** Object-language types
  , module Data.Int
  , module Data.Word
  , Complex (..)
  , PrimType'
  , PrimType
  , Type
  , Length
  , Index
  , Ref, DRef
  , Arr, DArr
  , IArr, DIArr
  , Inhabited
  , Syntactic
  , Domain
  , Internal
    -- * Front end
  , eval
  , showAST
  , drawAST
  , module Feldspar.Frontend
  , Bits
  , FiniteBits
  , Integral
  , Ord
  , RealFloat
  , RealFrac
  , Border (..)
  , IxRange
  , AssertionLabel (..)
  ) where

import Prelude.EDSL

import Control.Monad hiding (foldM)

import Data.Bits (Bits, FiniteBits)
import Data.Complex (Complex (..))
import Data.Int
import Data.Word

import Language.Syntactic (Syntactic, Domain, Internal)
import qualified Language.Syntactic as Syntactic

import Language.Embedded.Imperative (Border (..), IxRange)

import Data.Inhabited
import Feldspar.Primitive.Representation
import Feldspar.Representation
import Feldspar.Frontend
import Feldspar.Optimize

-- | Show the syntax tree using Unicode art
--
-- Only user assertions will be included.
showAST :: (Syntactic a, Domain a ~ FeldDomain) => a -> String
showAST = Syntactic.showAST . optimize onlyUserAssertions . Syntactic.desugar

-- | Draw the syntax tree on the terminal using Unicode art
--
-- Only user assertions will be included.
drawAST :: (Syntactic a, Domain a ~ FeldDomain) => a -> IO ()
drawAST = Syntactic.drawAST . optimize onlyUserAssertions . Syntactic.desugar

max_stars_repo_name: Abhiroop/mu-feldspar

max_issues_repo_name: Abhiroop/mu-feldspar

max_forks_repo_name: Abhiroop/mu-feldspar


text:
{-# LANGUAGE FlexibleInstances #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Foreign.BLAS.Level2
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Matrix-Vector operations.
--

module Foreign.BLAS.Level2 (
    BLAS2(..),
    ) where

import Data.Complex
import Foreign( Ptr, Storable, with )

import Foreign.BLAS.Types
import Foreign.BLAS.Level1
import Foreign.BLAS.Double
import Foreign.BLAS.Zomplex

-- | Types with matrix-vector operations.
class (BLAS1 a) => BLAS2 a where
    gbmv :: Trans -> Int -> Int -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    gemv :: Trans -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    gerc :: Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    geru :: Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()

    hbmv :: Uplo -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    hemv :: Uplo -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    her  :: Uplo -> Int -> Double -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    her2 :: Uplo -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    hpmv :: Uplo -> Int -> a -> Ptr a -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    hpr  :: Uplo -> Int -> Double -> Ptr a -> Int -> Ptr a -> IO ()
    hpr2 :: Uplo -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> IO ()

    tbmv :: Uplo -> Trans -> Diag -> Int -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    tbsv :: Uplo -> Trans -> Diag -> Int -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    tpmv :: Uplo -> Trans -> Diag -> Int -> Ptr a -> Ptr a -> Int -> IO ()
    tpsv :: Uplo -> Trans -> Diag -> Int -> Ptr a -> Ptr a -> Int -> IO ()
    trmv :: Uplo -> Trans -> Diag -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    trsv :: Uplo -> Trans -> Diag -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()

withEnum :: (Enum a, Storable a) => Int -> (Ptr a -> IO b) -> IO b
withEnum = with . toEnum
{-# INLINE withEnum #-}

instance BLAS2 Double where
    gemv transa m n alpha pa lda px incx beta py incy =
        withTrans transa $ \ptransa ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            dgemv ptransa pm pn palpha pa plda px pincx pbeta py pincy
    {-# INLINE gemv #-}

    gbmv transa m n kl ku alpha pa lda px incx beta py incy =
        withTrans transa $ \ptransa ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum kl $ \pkl ->
        withEnum ku $ \pku ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            dgbmv ptransa pm pn pkl pku palpha pa plda px pincx pbeta py pincy
    {-# INLINE gbmv #-}

    trmv uplo trans diag n pa lda px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
            dtrmv puplo ptrans pdiag pn pa plda px pincx
    {-# INLINE trmv #-}

    tpmv uplo trans diag n pap px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
            dtpmv puplo ptrans pdiag pn pap px pincx
    {-# INLINE tpmv #-}

    tpsv uplo trans diag n pap px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
            dtpsv puplo ptrans pdiag pn pap px pincx
    {-# INLINE tpsv #-}

    tbmv uplo trans diag n k pa lda px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
            dtbmv puplo ptrans pdiag pn pk pa plda px pincx
    {-# INLINE tbmv #-}

    trsv uplo trans diag n pa lda px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
            dtrsv puplo ptrans pdiag pn pa plda px pincx
    {-# INLINE trsv #-}

    tbsv uplo trans diag n k pa lda px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
            dtbsv puplo ptrans pdiag pn pk pa plda px pincx
    {-# INLINE tbsv #-}

    hemv uplo n alpha pa lda px incx beta py incy =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            dsymv puplo pn palpha pa plda px pincx pbeta py pincy
    {-# INLINE hemv #-}

    hbmv uplo n k alpha pa lda px incx beta py incy =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            dsbmv puplo pn pk palpha pa plda px pincx pbeta py pincy
    {-# INLINE hbmv #-}

    gerc m n alpha px incx py incy pa lda =
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
        withEnum lda $ \plda ->
            dger pm pn palpha px pincx py pincy pa plda
    {-# INLINE gerc #-}

    geru = gerc
    {-# INLINE geru #-}

    her uplo n alpha px incx pa lda =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum lda $ \plda ->
            dsyr puplo pn palpha px pincx pa plda
    {-# INLINE her #-}

    her2 uplo n alpha px incx py incy pa lda =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
        withEnum lda $ \plda ->
            dsyr2 puplo pn palpha px pincx py pincy pa plda
    {-# INLINE her2 #-}

    hpmv uplo n alpha pap px incx beta py incy =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            dspmv puplo pn palpha pap px pincx pbeta py pincy
    {-# INLINE hpmv #-}

    hpr uplo n alpha px incx pap =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
            dspr puplo pn palpha px pincx pap
    {-# INLINE hpr #-}

    hpr2 uplo n alpha px incx py incy pap =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
            dspr2 puplo pn palpha px pincx py pincy pap
    {-# INLINE hpr2 #-}

instance BLAS2 (Complex Double) where
    gemv transa m n alpha pa lda px incx beta py incy =
        withTrans transa $ \ptransa ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            zgemv ptransa pm pn palpha pa plda px pincx pbeta py pincy
    {-# INLINE gemv #-}

    gbmv transa m n kl ku alpha pa lda px incx beta py incy =
        withTrans transa $ \ptransa ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum kl $ \pkl ->
        withEnum ku $ \pku ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            zgbmv ptransa pm pn pkl pku palpha pa plda px pincx pbeta py pincy
    {-# INLINE gbmv #-}

    trmv uplo trans diag n pa lda px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
            ztrmv puplo ptrans pdiag pn pa plda px pincx
    {-# INLINE trmv #-}

    tpmv uplo trans diag n pap px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
            ztpmv puplo ptrans pdiag pn pap px pincx
    {-# INLINE tpmv #-}

    tpsv uplo trans diag n pap px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
            ztpsv puplo ptrans pdiag pn pap px pincx
    {-# INLINE tpsv #-}

    tbmv uplo trans diag n k pa lda px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
            ztbmv puplo ptrans pdiag pn pk pa plda px pincx
    {-# INLINE tbmv #-}

    trsv uplo trans diag n pa lda px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
            ztrsv puplo ptrans pdiag pn pa plda px pincx
    {-# INLINE trsv #-}

    tbsv uplo trans diag n k pa lda px incx =
        withUplo uplo $ \puplo ->
        withTrans trans $ \ptrans ->
        withDiag diag $ \pdiag ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
            ztbsv puplo ptrans pdiag pn pk pa plda px pincx
    {-# INLINE tbsv #-}

    hemv uplo n alpha pa lda px incx beta py incy =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            zhemv puplo pn palpha pa plda px pincx pbeta py pincy
    {-# INLINE hemv #-}

    hbmv uplo n k alpha pa lda px incx beta py incy =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            zhbmv puplo pn pk palpha pa plda px pincx pbeta py pincy
    {-# INLINE hbmv #-}

    gerc m n alpha px incx py incy pa lda =
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
        withEnum lda $ \plda ->
            zgerc pm pn palpha px pincx py pincy pa plda
    {-# INLINE gerc #-}

    geru m n alpha px incx py incy pa lda =
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
        withEnum lda $ \plda ->
            zgeru pm pn palpha px pincx py pincy pa plda
    {-# INLINE geru #-}

    her uplo n alpha px incx pa lda =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum lda $ \plda ->
            zher puplo pn palpha px pincx pa plda
    {-# INLINE her #-}

    her2 uplo n alpha px incx py incy pa lda =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
        withEnum lda $ \plda ->
            zher2 puplo pn palpha px pincx py pincy pa plda
    {-# INLINE her2 #-}

    hpmv uplo n alpha pap px incx beta py incy =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        with beta $ \pbeta ->
        withEnum incy $ \pincy ->
            zhpmv puplo pn palpha pap px pincx pbeta py pincy
    {-# INLINE hpmv #-}

    hpr uplo n alpha px incx pap =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
            zhpr puplo pn palpha px pincx pap
    {-# INLINE hpr #-}

    hpr2 uplo n alpha px incx py incy pap =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
            zhpr2 puplo pn palpha px pincx py pincy pap
    {-# INLINE hpr2 #-}

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas