text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE OverloadedStrings #-}

module Main where

import qualified Statistics.Distribution as Stat
import qualified Statistics.Distribution.Normal as Stat
import qualified Statistics.Distribution.Uniform as Stat
import Control.Monad.Primitive
import qualified System.Random.MWC as Rand
import Control.Monad.State.Strict
import Control.Monad.Reader
import Lens.Micro
import Types
import qualified Data.List as List
import Options.Generic

-- Contraints that allow for random number generation
type BanditM m = (MonadIO m, MonadReader (Rand.Gen (PrimState IO)) m)

-- |Generate a list of random bandits using dist to generate their parameters.
randomBandits ::
  (BanditM m, Stat.ContGen d1, Stat.ContGen d2)
  => Int
  -> d1
  -> d2
  -> m [Bandit BanditDist]
randomBandits count avgDist stdDevDist =
  replicateM count randomBandit
  where
    randomBandit = do
      gen <- ask
      mean <- liftIO $ Stat.genContVar avgDist gen
      stdDev <- liftIO $ Stat.genContVar stdDevDist gen
      return $ Bandit (Stat.normalDistr mean (abs stdDev))

-- |Pulls the handle and returns the reward
pullArm ::
  (BanditM m, Stat.ContGen d)
  => Bandit d
  -> m BanditStats
pullArm (Bandit d) = do
  gen <- ask
  reward <- liftIO $ Stat.genContVar d gen
  return $ BanditStats reward 1

simulate ::
  BanditM m
  => SimState -- ^ Starting state
  -> Int      -- ^ Number of rounds
  -> (SimState -> m Int) -- ^ Choosing function
  -> m SimState
simulate simState 0 _ = return simState
simulate simState rounds f = do
  chosenIndex <- f simState
  let Just (chosenStat, chosenBandit) = simState ^? bandits. ix chosenIndex -- Pick a bandit
  newStat <- pullArm chosenBandit
  let newState = simState
                 & (bandits . ix chosenIndex) . _1 .~ (chosenStat `mappend` newStat) -- Update stats
                 & roundsCount %~ (+1) -- Update rounds count
                 & regretHist %~ (\t -> (regret simState) : t)
  simulate newState (rounds - 1) f

-- | Initialize the state but creating the specified number of bandits and pulling each
-- | one once.
initState ::
  BanditM m
  => Int -- Number of bandits
  -> m SimState
initState count = do
  initBandits <- randomBandits count (Stat.uniformDistr 0 50) (Stat.normalDistr 0 5)
  let initStats = replicate count mempty
  return $ SimState (zip initStats initBandits) 0 []

regret :: SimState -> Double
regret simState =
  let allTimesPulled = fromIntegral <$> (simState ^.. bandits . each . _1 . timesPulled)
      allMeans = simState ^.. bandits . each . _2 . distLens . meanLens
      actualExpectedReturn = sum (zipWith (*) allTimesPulled allMeans)
      maxMean = maximum (simState ^.. bandits . each . _2 . distLens . meanLens)
      rounds = fromIntegral $ simState ^. roundsCount
      maxExpectedReturn = maxMean * rounds
  in
    maxExpectedReturn - actualExpectedReturn

type Concrete a = (ReaderT (Rand.Gen (PrimState IO)) IO) a

run ::
  Int
  -> Int
  -> (SimState -> Concrete Int)
  -> Concrete SimState
run rounds count f = do
  startingState <- initState count
  simulate startingState rounds f

-- | Return the index of the largest item as determined by the function
argMaxIndex :: (Ord b, Num b) => [a] -> (a -> b) -> Maybe Int
argMaxIndex as f =
  let magnitudes = f <$> as
      indexAndMagnitude = zip [0..] magnitudes
      biggestToSmallest = (List.sortOn (negative . snd) indexAndMagnitude)
  in biggestToSmallest ^? _head . _1
  where
    negative x = -x -- Swap the sign, so it's sorted biggest to smallest.

-- | Choose the bandit with the highest Upper Confidence Bound
choose :: SimState -> Int
choose simState =
  let banditStats = simState ^.. bandits . each . _1
      curRound = simState ^. roundsCount
      Just maxIndex = argMaxIndex banditStats (ucb curRound)
  in maxIndex

ucb :: Int -> BanditStats -> Double
ucb curRound banditStats =
  let totalTimesPulled = fromIntegral (banditStats ^. timesPulled)
      empiricalMean = (banditStats ^. totalReward) / totalTimesPulled
      t = fromIntegral curRound
      confidence = sqrt ((2.0 * log (1 + (t * ((log t) ^ 2)))) / totalTimesPulled)
      in
    if totalTimesPulled /= 0 then empiricalMean + confidence else maxDouble
  where
    maxDouble = 100000000 -- TODO: What's the actual max double?

-- | Choose a random bandit.
chooseRandom :: BanditM m => SimState -> m Int
chooseRandom simState = do
    let count = fromIntegral $ length $ simState ^. bandits
    let d = Stat.uniformDistr 0 (count - 1)
    gen <- ask
    rand <- liftIO $ Stat.genContVar d gen
    return $ truncate rand

-- | Pick the bandit with the best empirical mean.
chooseBestMean :: SimState -> Int
chooseBestMean simState =
  let banditStats = simState ^.. bandits . each . _1
      Just maxIndex = argMaxIndex banditStats avgReturn
  in
    maxIndex

-- | Pick the bandit with the highest mean reward.
godMode :: SimState -> Int
godMode simState =
  let bandits_ = simState ^.. bandits . each . _2
      Just maxIndex = argMaxIndex bandits_ (\(Bandit d) -> Stat.mean d)
  in
    maxIndex

avgReturn :: BanditStats -> Double
avgReturn (BanditStats reward pulled) = reward / (fromIntegral pulled)

main :: IO ()
main = do
  (opts :: BanditOpts Unwrapped) <-unwrapRecord "Multi-armed bandit simulator."
  rand <- Rand.createSystemRandom
  finalState <- runReaderT (run (optRounds opts) (optBandits opts) (return . godMode)) rand
  _ <- traverse (putStrLn . show) (reverse (finalState ^. regretHist))
  return ()

data BanditOpts w =
  BanditOpts { optRounds :: w ::: Int <?> "Number of rounds to simulate."
             , optBandits :: w ::: Int <?> "Number of bandits."
             -- , optMeanMean :: w ::: Double <?> "Mean of the mean."
             -- , optStdDevMean :: w ::: Double <?> "StdDev of the mean."
             -- , optMeanStdDev :: w ::: Double <?> "Mean of the StdDev."
             -- , optStdDevStdDev :: w ::: Double <?> "StdDev of the StdDev."
             } deriving (Generic)

instance ParseRecord (BanditOpts Wrapped)

max_stars_repo_name: beala/bandits

max_issues_repo_name: beala/bandits

max_forks_repo_name: beala/bandits


text:
{-#LANGUAGE OverloadedStrings#-}

{-
Project name: Pipeline for analyze DropSeq data
Min Zhang
Date: Feb 5, 2016
Version: v0.1.0
README: Process raw fastq file, QC and mapping
-}

module DropSeq
where

import qualified Data.Text.Lazy as T
import qualified Data.Text.Lazy.IO as TextIO
import qualified Data.Char as C
import Control.Applicative
import qualified Data.List as L
import qualified Data.Set as Set
import Control.Monad (fmap)
import Data.Ord (comparing)
import Data.Function (on)
import qualified Safe as S
import qualified Data.HashMap.Lazy as M
import qualified Data.Maybe as Maybe
import qualified Data.Foldable as F (all)
import Data.Traversable (sequenceA)

import qualified Data.Vector as V
import qualified Statistics.Sample.Histogram as Stats

import MyText
import MyTable
import Util
import IO
import DataTypes

{- Part I
1. Check length of reads
2. Check number of unique cell barcode
3. Check number of unique UMI
-}

testR1Path = "/Users/minzhang/Documents/data/P56_dropseq/MH-L1-32920888/Data/Intensities/BaseCalls/MH-L1_S1_L001_R1_001.fastq"
testR2Path = "/Users/minzhang/Documents/data/P56_dropseq/MH-L1-32920888/Data/Intensities/BaseCalls/MH-L1_S1_L001_R2_001.fastq"

dropSeqQc = do
  inputFastq1 <- take 10000 <$> importFastq testR1Path
  let readLen = (V.fromList $ map (fromIntegral . T.length . fqseq) inputFastq1)
  let results = (Stats.histogram 5 readLen) :: (V.Vector Double, V.Vector Int)
  print results

-- high demand on memory
ifFastqNameMatch fq1 fq2 = and $ zipWith (==) (L.sort $ map fqname fq1) (L.sort $ map fqname fq2)

-- import fq1 only seq, to check the diversity of the barcodes.
testR1SeqPath = "/Users/minzhang/Documents/data/P56_dropseq/MH-L1-32920888/Data/Intensities/BaseCalls/MH-L1_S1_L001_R1_001.test.seq.fastq"
r1SeqPath = "/Users/minzhang/Documents/data/P56_dropseq/MH-L1-32920888/Data/Intensities/BaseCalls/MH-L1_S1_L001_R1_001.seq.fastq"

barcodeQc inputPath = do
  allBarcodes <- map (\x->(T.take 7 $ T.drop 5 x, T.take 8 $ T.drop 12 x)) . T.lines <$> TextIO.readFile inputPath
  let results = (Set.size $ Set.fromList $ map fst allBarcodes, Set.size $ Set.fromList $ map snd allBarcodes)
  print results

sortUMI inputPath = do
  results <- map (T.splitAt 12 . T.init) . T.lines <$> TextIO.readFile inputPath
  let output = T.unlines $ map (T.pack . show . length) $ L.reverse $ L.sortBy (comparing length) $ L.groupBy ((==) `on` snd) $ L.sortBy (comparing snd) results
  TextIO.writeFile "outputUMIsize.txt" output

ext: hs

max_stars_repo_name: Min-/DropSeq

max_issues_repo_name: Min-/DropSeq

max_forks_repo_name: Min-/DropSeq


text:
module STCR2Z2T0S0EndPoint where

import           Control.Monad             as M
import           Data.Array.Repa           as R
import           Data.Binary               (decodeFile)
import           Data.Complex
import           Data.List                 as L
import           DFT.Plan
import           FokkerPlanck.DomainChange
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           Image.Transform           (normalizeValueRange)
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           System.Random
import           Text.Printf
import           Types
import           Utils.Array
import           Utils.Parallel

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:histFilePath:numIterationStr:writeSourceFlagStr:cutoffRadiusEndPointStr:cutoffRadiusStr:reversalFactorStr:cStr:patchNormFlagStr:patchNormSizeStr:approximatedEigenValueStr:shape2DStr:segmentsFilePath:segIdxStr:useFFTWWisdomFlagStr:fftwWisdomFileName:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      cutoffRadiusEndPoint = read cutoffRadiusEndPointStr :: Int
      cutoffRadius = read cutoffRadiusStr :: Int
      reversalFactor = read reversalFactorStr :: Double
      patchNormFlag = read patchNormFlagStr :: Bool
      patchNormSize = read patchNormSizeStr :: Int
      approximatedEigenValue = read approximatedEigenValueStr :: Double
      shape2D@(Points _ minDist _) = read shape2DStr :: Points Shape2D
      segIdx = read segIdxStr :: Int
      useFFTWWisdomFlag = read useFFTWWisdomFlagStr :: Bool
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2Z2T0S0EndPoint"
      a = 20 :: Int
      b = 5 :: Int
      c = read cStr :: Int
      endPointFilePath =
        folderPath </>
        (printf
           "EndPoint_%d_%d_%d_%d_%d_%d_%.2f_%.2f_%d_%d_%d_%f.dat"
           numPoint
           (round thetaFreq :: Int)
           (round scaleFreq :: Int)
           (round maxScale :: Int)
           (round tao :: Int)
           cutoffRadiusEndPoint
           thetaSigma
           scaleSigma
           a
           b
           c
           reversalFactor)
      fftwWisdomFilePath = folderPath </> fftwWisdomFileName
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      (PinwheelHollow0 10)
      (cutoff cutoffRadius radialArr)
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  plan <-
    makeR2Z2T0S0Plan emptyPlan useFFTWWisdomFlag fftwWisdomFilePath arrR2Z2T0S0
  -- (plan, pathNormMethod) <-
  --   makePatchNormFilter plan' numPoint numPoint patchNormFlag patchNormSize
      -- minDist = 8
      -- kanizsaTriangle1 =
      --   makeShape2D $ Points (-30, -30) minDist (Corner 30 60 80) --(PacMan 30 60 50)  --(Ehrenstein 8 15 40)  -- (IncompleteCircle 0 60 50 ) -- (TJunction 45 50) -- (PacMan 0 60 100 ) -- (Corner 0 60 100) --(KanizsaTriangle1 0 480 160 80)
  -- ys <-
  --   (\aa -> aa L.!! segIdx) <$> decodeFile segmentsFilePath :: IO [(Int, Int)]
  let pointSet = makeShape2D shape2D
      shapeArr =
        getShape2DRepaArray
          numPoint
          numPoint
          (L.map
             (\(x, y) ->
                (x + fromIntegral numPoint / 2, y + fromIntegral numPoint / 2))
             pointSet)
      xs =
        L.map (\(x, y) -> R2S1RPPoint (x, y, 0, 1)) . getShape2DIndexList $
        pointSet
      -- (xAvg, yAvg) =
      --   (\(as, bs) ->
      --      ( round $
      --        (fromIntegral . L.sum $ as) / (fromIntegral . L.length $ as)
      --      , round $
      --        (fromIntegral . L.sum $ bs) / (fromIntegral . L.length $ bs))) .
      --   L.unzip $
      --   ys
      -- centeredYs = L.map (\(x, y) -> (x - xAvg, y - yAvg)) ys
      -- xs = L.map (\(x, y) -> R2S1RPPoint (x, y, 0, 1)) centeredYs
      -- shapeArr =
      --   L.head . cluster2Array numPoint numPoint $
      --   [L.map (\(x, y) -> (x + div numPoint 2, y + div numPoint 2)) centeredYs]
      -- pointSet = L.map (\(x, y) -> (fromIntegral x, fromIntegral y)) centeredYs
  let bias = computeBiasR2T0S0 numPoint numPoint theta0Freqs scale0Freqs xs
      eigenVec =
        computeInitialEigenVectorR2T0S0
          numPoint
          numPoint
          theta0Freqs
          scale0Freqs
          thetaFreqs
          scaleFreqs
          xs
  plotImageRepa (folderPath </> "Shape.png") . ImageRepa 8 $ shapeArr
  endPointFlag <- doesFileExist endPointFilePath
  completionFieldR2Z2'' <-
    if endPointFlag
      then readRepaArray endPointFilePath
      else (do putStrLn "Couldn't find the endpoint data. Start computing..."
               arrR2Z2T0S0EndPoint <-
                 computeUnboxedP $
                 computeR2Z2T0S0ArrayRadial
                   -- pinwheel
                   -- (pinwheelHollowNonzeronCenter 16)
                   (PinwheelHollow0 4)
                   (cutoff cutoffRadiusEndPoint radialArr)
                   numPoint
                   numPoint
                   1
                   maxScale
                   thetaFreqs
                   scaleFreqs
                   theta0Freqs
                   scale0Freqs
               pathNormMethod <-
                 if patchNormFlag
                   then do
                     let points =
                           createIndex2D .
                           L.map
                             (\(i, j) ->
                                (i + div numPoint 2, j + div numPoint 2)) .
                           getShape2DIndexList $
                           pointSet
                         ys =
                           pointCluster
                             (connectionMatrixP
                                (ParallelParams numThread 1)
                                (minDist + 1)
                                points) $
                           points
                     M.zipWithM_
                       (\i ->
                          plotImageRepa
                            (folderPath </> (printf "Cluster%03d.png" i)) .
                          ImageRepa 8)
                       [1 :: Int ..] .
                       cluster2Array numPoint numPoint $
                       ys
                     return . PowerMethodConnection $ ys
                   else return PowerMethodGlobal
               (R.foldAllP max 0 . R.map magnitude $ arrR2Z2T0S0EndPoint) >>=
                 print
               completionFieldR2Z2' <-
                 powerMethodR2Z2T0S0Reversal
                   plan
                   folderPath
                   numPoint
                   numPoint
                   numOrientation
                   thetaFreqs
                   theta0Freqs
                   numScale
                   scaleFreqs
                   scale0Freqs
                   maxScale
                   arrR2Z2T0S0EndPoint
                   pathNormMethod
                   numIteration
                   writeSourceFlag
                   (printf
                      "_%d_%d_%d_%d_%d_%d_%.2f_%.2f_%f_EndPoint"
                      numPoint
                      (round thetaFreq :: Int)
                      (round scaleFreq :: Int)
                      (round maxScale :: Int)
                      (round tao :: Int)
                      cutoffRadiusEndPoint
                      thetaSigma
                      scaleSigma
                      reversalFactor)
                   0.5
                   reversalFactor
                   bias
                   eigenVec
               -- writeRepaArray endPointFilePath completionFieldR2Z2'
               return completionFieldR2Z2')
  let completionFieldR2Z2 = R.zipWith (*) completionFieldR2Z2'' bias
      endpointBias = rotateBiasR2Z2T0S0 180 theta0Freqs completionFieldR2Z2
      -- endpointBias =
      --   R.zipWith
      --     (+)
      --     (rotateBiasR2Z2T0S0 90 theta0Freqs completionFieldR2Z2)
      --     (rotateBiasR2Z2T0S0 (-90) theta0Freqs completionFieldR2Z2)
                    -- rotateBiasR2Z2T0S0 0 theta0Freqs . R.traverse completionFieldR2Z2 id $ \f idx@(Z :. _ :. _ :. i :. j) ->
                    --   if (sqrt . fromIntegral $
                    --       (i - div numPoint 2) ^ 2 + (j - div numPoint 2) ^ 2) >
                    --      35
                    --     then 0
                    --     else f idx
      biasMag =
        R.sumS .
        R.sumS .
        rotate4D .
        rotate4D .
        R.map magnitude .
        r2z2Tor2s1rp numOrientation thetaFreqs numScale scaleFreqs $
        endpointBias
  plotImageRepa (folderPath </> "EndPointBias.png") .
    ImageRepa 8 . computeS . R.extend (Z :. (1 :: Int) :. All :. All) $
    biasMag
  printf "%f %f\n" reversalFactor (R.sumAllS biasMag)
  -- powerMethodR2Z2T0S0BiasReversal
  --   plan
  --   folderPath
  --   numPoint
  --   numPoint
  --   numOrientation
  --   thetaFreqs
  --   theta0Freqs
  --   numScale
  --   scaleFreqs
  --   scale0Freqs
  --   arrR2Z2T0S0
  --   -- numIteration
  --   10
  --   writeSourceFlag
  --   (printf
  --      "_%d_%d_%d_%d_%d_%d_%.2f_%.2f"
  --      numPoint
  --      (round thetaFreq :: Int)
  --      (round scaleFreq :: Int)
  --      (round maxScale :: Int)
  --      (round tao :: Int)
  --      cutoffRadius
  --      thetaSigma
  --      scaleSigma)
  --   0.1
  --   (computeS endpointBias)
  --   (R.fromFunction
  --      (Z :. (L.length thetaFreqs) :. (L.length scaleFreqs) :.
  --       (L.length theta0Freqs) :.
  --       (L.length scale0Freqs) :.
  --       numPoint :.
  --       numPoint) $ \(Z :. k :. l :. _ :. _ :. i :. j) ->
  --      if k == div (L.length thetaFreqs) 2 && l == div (L.length scaleFreqs) 2
  --        then 1 / (fromIntegral $ (L.length theta0Freqs * L.length scale0Freqs)) :+
  --             0
  --        else 0)

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE BangPatterns #-}

module Pact.GasModel.GasModel where

import Control.Exception (bracket)
import Control.Monad (void, replicateM)
import Data.List (foldl', sortOn)
import GHC.Conc (numCapabilities)
import Statistics.Types (Estimate(..))

import qualified Criterion.Main as C
import qualified Criterion as C
import qualified Criterion.Types as C
import qualified Data.ByteString.Lazy.Char8 as BSL8
import qualified Data.Csv as Csv
import qualified Data.HashMap.Strict as HM
import qualified Data.Text as T
import qualified Options.Applicative as O

import Pact.Eval (eval)
import Pact.GasModel.GasTests
import Pact.GasModel.Types
import Pact.GasModel.Utils
import Pact.Types.Lang
import Pact.Types.Runtime hiding (GasPrice)

data Option = Option
  { _oFilter :: Maybe T.Text
  , _oBenchOnly :: Bool
  }
  deriving (Eq,Show)

options :: O.ParserInfo Option
options = O.info (O.helper <*> parser)
          (O.fullDesc <> O.header "Gas Model Benchmarks")
  where
    parser =
      Option
      <$> ((Just <$> O.strOption
            (O.short 'f' <> O.long "filter" <> O.metavar "NATIVE" <>
             O.help "Pipe-delimited list of natives to run")) O.<|>
            pure Nothing)
      <*> (O.flag False True
           (O.short 'b' <> O.long "bench" <> O.help "Just bench and display gas price"))

optionToGasTests
  :: Option
  -> [(NativeDefName, GasUnitTests)]
optionToGasTests opt = HM.toList $
  case _oFilter opt of
    Nothing -> unitTests
    Just ts -> HM.filterWithKey matchingFuns unitTests
      where tlist = NativeDefName <$> T.split (== '|') ts
            matchingFuns k _ = k `elem` tlist

type NanoSeconds = Double

secToNs :: Double -> NanoSeconds
secToNs dInSeconds =
  dInSeconds * 1000000000

newtype Average = Average Double
instance Csv.ToField Average where
  toField (Average avg) = Csv.toField avg

average :: (Foldable f) => f Double -> Average
average nums = Average $! totalSum / (fromIntegral size)
  where (totalSum,size) = foldl' f (0.0, 0 :: Int) nums
        f (!sum', !size') x = (sum' + x, size' + 1)

newtype GasPrice = GasPrice Integer
  deriving (Show)
instance Csv.ToField GasPrice where
  toField (GasPrice p) = Csv.toField p

gasPriceDivisor :: Double
gasPriceDivisor = 2500.0

gasPrice :: [NanoSeconds] -> GasPrice
gasPrice measurements = GasPrice price
  where
    (Average avg) = average measurements
    price = ceiling (avg / gasPriceDivisor)

exec
  :: EvalState
  -> EvalEnv e
  -> [Term Name]
  -> IO (Either PactError [Term Name], EvalState)
exec s e terms = runEval' s e $ mapM eval terms

mockRun
  :: PactExpression
  -> (EvalEnv e, EvalState)
  -> IO (Either PactError [Term Name], EvalState)
mockRun expr (env, state) = do
  terms <- compileCode (_pactExpressionFull expr)
  exec state env terms

setupEnv :: GasSetup e -> IO (EvalEnv e, EvalState)
setupEnv (GasSetup e s _ _) = do
  env <- e
  state <- s
  return (env, state)

mockRuns
  :: GasUnitTests
  -> IO ()
mockRuns tests = do
  void $ mapOverGasUnitTests tests run run
  where
    run expr dbSetup = do
      putStrLn $ "Dry run for " ++ T.unpack (getDescription expr dbSetup)

      (res,_) <- bracket (setupEnv dbSetup)
                         (gasSetupCleanup dbSetup)
                         (mockRun expr)
      eitherDie (getDescription expr dbSetup) res

-- | For debugging purposes
mockRuns'
  :: GasUnitTests
  -> IO [GasTestResult (Either PactError [Term Name], EvalState)]
mockRuns' tests = do
  mapOverGasUnitTests tests run run
  where
    run expr dbSetup = do
      bracket (setupEnv dbSetup) (gasSetupCleanup dbSetup) (mockRun expr)

bench
  :: PactExpression
  -> GasSetup e
  -> IO NanoSeconds
bench expr dbSetup = do
  terms <- compileCode (_pactExpressionFull expr)
  putStrLn $ T.unpack (getDescription expr dbSetup)
  report <- bracket setup teardown $ C.benchmark' . (run terms)
  return $ secToNs $
           estPoint $
           C.anMean $
           C.reportAnalysis report
  where
    setup = do
      s <- setupEnv dbSetup
      return $ NoopNFData s
    teardown (NoopNFData env) = do
      (gasSetupCleanup dbSetup) env
    run terms ~(NoopNFData (env, state)) =
          C.nfIO (exec state env terms)

benchesOnce
  :: GasUnitTests
  -> IO [GasTestResult Average]
benchesOnce tests =
  mapOverGasUnitTests tests run run
  where
    run expr setup = do
      means <- replicateM 1 $ bench expr setup
      return $! average means

-- | Benchmarks each native function's tests three times
benchesMultiple
  :: (NativeDefName, GasUnitTests)
  -> IO (NativeDefName, [GasTestResult Average])
benchesMultiple (funName, tests) = do
  results <- mapOverGasUnitTests tests run run
  return (funName, results)
  where
    batchSize = 3
    run expr setup = do
      means <- replicateM batchSize $ bench expr setup
      return $! average means

meansToCSVEncoding
  :: [GasTestResult Average]
  -> (Csv.Header, [Csv.NamedRecord])
meansToCSVEncoding results =
  (Csv.header headersField, map toNamedRecord results)
  where
    headersField =
      map Csv.toField (["function", "description",
                        "sqliteDbAverage (ns)", "mockDbAverage (ns)"] :: [T.Text])
    toNamedRecord (GasTestResult funName desc sqlAvg mockAvg) =
      Csv.namedRecord $
      zip headersField [Csv.toField (asString funName),
                        Csv.toField desc,
                        Csv.toField sqlAvg,
                        Csv.toField mockAvg
                       ]

writeRawCSV
  :: [GasTestResult Average]
  -> IO ()
writeRawCSV results = do
  let (headers, records) = meansToCSVEncoding results
      content = Csv.encodeByName headers records
  BSL8.writeFile "gas-model-raw-data.csv" content

sqliteGasPrice :: [GasTestResult Average] -> GasPrice
sqliteGasPrice rawResults = gasPrice dataPoints
  where
    dataPoints =
      map (liftAverage . _gasTestResultSqliteDb) onlySimpleTestResults
      where liftAverage (Average d) = d

    onlySimpleTestResults =
      filter simpleTestCheck rawResults
      where
        simpleTestCheck (GasTestResult _ desc _ _) =
          (not $ T.isInfixOf "med" desc) &&
          (not $ T.isInfixOf "long" desc)

encodeGasPrice
  :: [(NativeDefName, GasPrice)]
  -> (Csv.Header, [Csv.NamedRecord])
encodeGasPrice allPrices =
  (headers,
   map encodeFunctionPrice allPrices
  )
  where
    priceHeader =
      "gas price (" <> T.pack (show gasPriceDivisor) <> " ns)"
    functionNameHeader = "function"
    headers = Csv.header $
      map Csv.toField [functionNameHeader, priceHeader]

    encodeFunctionPrice (funName, price) =
      Csv.namedRecord $
      [(Csv.toField functionNameHeader) Csv..= asString funName] <>
      [(Csv.toField priceHeader) Csv..= price]

gasPriceExplanation :: [T.Text]
gasPriceExplanation = [purpose, numCores, os, hardwareNotes,
                       implementation, numOfIterations, backend, calcPrice, ""]
  where
    purpose = "Purpose: Calculate the gas price of Pact native functions using a data-driven model"
    numCores = "Number of cores available for benchmarks: " <> T.pack (show numCapabilities)
    os = "Operating System: <INSERT OS>"
    hardwareNotes = "More information on hardware used: <INSERT MORE INFORMATION ON HARDWARE USED>"
    implementation = "Implementation: For every native function, executes and benchmarks "
                     <> "(using Criterion) simple examples of said function."
    numOfIterations = "Number of iterations: Each function's tests are run three times against each backend type and averaged."
    backend = "Benchmark backend(s): [Sqlite db with `fastNoJournalPragmas`]"
    calcPrice = "From benchmark to a native function's price: For every backend type, all of the means of the function's"
                <>" benchmark examples are converted into nanoseconds, averaged together, divided by "
                <> T.pack (show gasPriceDivisor)
                <> ", and rounded up to the nearest integer."
                <> "Only the benchmark averages of the simple function tests are used when calculating the function's price."

writeGasPriceCSV
  :: [(NativeDefName, [GasTestResult Average])]
  -> IO ()
writeGasPriceCSV allResults = do
  let gasPrices = map (\(n,res) -> (n, sqliteGasPrice res)) allResults
      (headers, records) = encodeGasPrice gasPrices
      content = Csv.encodeByName headers records
      explanation = Csv.encode $
        map (\t -> Csv.record [Csv.toField t]) gasPriceExplanation
  BSL8.writeFile "gas-prices.csv" (explanation <> content)

coverageReport :: IO ()
coverageReport = do
  print $ "Missing benchmark tests for "
           ++ show (length untestedNatives)
           ++ " out of "
           ++ show (length allNatives)
           ++ " native functions."
  mapM_ (print . show) untestedNatives

main :: IO ()
main = do
  opt <- O.execParser options
  let tests = optionToGasTests opt

  -- Enforces that unit tests succeed
  putStrLn "Doing dry run of benchmark tests"
  mapM_ (mockRuns . snd) tests

  putStrLn "Running benchmark(s)"

  if _oBenchOnly opt then
    let displayGasPrice (funName, t) = do
          res <- benchesOnce t
          putStrLn $ (T.unpack (asString funName)) ++ ": " ++ show (sqliteGasPrice res)
          putStrLn ""
    in mapM_ displayGasPrice tests
  else do
    let testsSorted = sortOn fst tests
    allBenches <- mapM benchesMultiple testsSorted

    putStrLn "Exporting raw benchmarks data"
    writeRawCSV (concatMap snd allBenches)

    putStrLn "Exporting data-driven gas prices"
    writeGasPriceCSV allBenches

    putStrLn "Reporting coverage"
    coverageReport

-- | For debugging
runGasTestByName :: T.Text -> IO [GasTestResult (Either PactError [Term Name], EvalState)]
runGasTestByName nname = do
  case (unitTestFromDef (NativeDefName nname)) of
    Nothing -> eitherDie nname (Left ("gas tests not found" :: T.Text))
    Just g -> mockRuns' g

runAllGasTests :: IO ()
runAllGasTests = do
  mapM_ (mockRuns . snd) (HM.toList unitTests)

max_stars_repo_name: edmundnoble/pact

max_issues_repo_name: edmundnoble/pact

max_forks_repo_name: edmundnoble/pact


text:
{-# LANGUAGE DeriveDataTypeable, GeneralizedNewtypeDeriving, Rank2Types #-}
{-# OPTIONS_HADDOCK hide #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Matrix.Base
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--

module Numeric.LinearAlgebra.Matrix.Base (
    Matrix(..),
    dim,

    fromList,
    fromRow,
    fromCol,
    zero,
    constant,

    at,
    unsafeAt,

    indices,
    elems,
    assocs,

    col,
    unsafeCol,
    cols,

    update,
    unsafeUpdate,
    accum,
    unsafeAccum,

    map,
    zipWith,
    unsafeZipWith,

    slice,
    unsafeSlice,

    splitRowsAt,
    dropRows,
    takeRows,

    splitColsAt,
    dropCols,
    takeCols,

    fromVector,
    toVector,

    isContig,
    unsafeFromForeignPtr,
    unsafeToForeignPtr,
    unsafeWith,

    ) where

import Prelude hiding ( read, map, zipWith )
import qualified Prelude as P

import Data.AEq( AEq(..) )
import Data.Typeable( Typeable )
import Foreign( ForeignPtr, Ptr, Storable )
import Text.Printf( printf )

import Numeric.LinearAlgebra.Vector( Vector )
import qualified Numeric.LinearAlgebra.Vector as V

-- | Immutable dense matrices.
data Matrix e =
      Matrix {-# UNPACK #-} !(Vector e) -- matrix data
             {-# UNPACK #-} !Int        -- row dimension
             {-# UNPACK #-} !Int        -- column dimension
             {-# UNPACK #-} !Int        -- leading dimension
  deriving (Typeable)

-- | Get the matrix dimensions (number of rows and number of columns).
dim :: (Storable e) => Matrix e -> (Int,Int)
dim (Matrix _ m n _) = (m,n)
{-# INLINE dim #-}

-- | Indicates if the elements of the matrix are stored contigously
isContig :: (Storable e) => Matrix e -> Bool
isContig (Matrix _ m _ lda) = lda == m || m == 0
{-# INLINE isContig #-}

-- | Create a matrix of the given dimension with elements initialized
-- to the values from the list, in column major order.
fromList :: (Storable e) => (Int,Int) -> [e] -> Matrix e
fromList (m,n) es
    | m < 0 || n < 0 = error $
        printf "fromList (%d,%d): negative dimension" m n
    | otherwise = let
        v = V.fromList (m*n) es
        lda = max 1 m
        in Matrix v m n lda
{-# INLINE fromList #-}

-- | Create a matrix of the given dimension with all elements initialized
-- to the given value
constant :: (Storable e) => (Int,Int) -> e -> Matrix e
constant (m,n) e
    | m < 0 || n < 0 = error $
        printf "constant (%d,%d): negative dimension" m n
    | otherwise = let
        v = V.constant (m*n) e
        lda = max 1 m
        in Matrix v m n lda
{-# INLINE constant #-}

-- | Create a zero of the given dimension with all elements initialized
-- to zero.
zero :: (Storable e, Num e) => (Int,Int) -> Matrix e
zero (m,n)
    | m < 0 || n < 0 = error $
        printf "zero (%d,%d): negative dimension" m n
    | otherwise = let
        v = V.zero (m*n)
        lda = max 1 m
        in Matrix v m n lda
{-# INLINE zero #-}

-- | Returns the element of a matrix at the specified index.
at :: (Storable e) => Matrix e -> (Int,Int) -> e
at a ij@(i,j)
    | i < 0 || i >= m || j < 0 || j >= n = error $
        printf ("at <matrix with dim (%d,%d)> (%d,%d):"
                ++ " invalid index") m n i j
    | otherwise =
        unsafeAt a ij
  where
      (m,n) = dim a
{-# INLINE at #-}

unsafeAt :: (Storable e) => Matrix e -> (Int,Int) -> e
unsafeAt (Matrix v _ _ lda) (i,j) =
    V.unsafeAt v (i + j * lda)
{-# INLINE unsafeAt #-}

-- | Get the indices of the elements in the matrix, in column-major order.
indices :: (Storable e) => Matrix e -> [(Int,Int)]
indices a = [ (i,j) | j <- [ 0..n-1 ], i <- [ 0..m-1 ] ]
  where (m,n) = dim a

-- | Returns a list of the elements of a matrix, in the same order as their
-- indices.
elems :: (Storable e) => Matrix e -> [e]
elems (Matrix v m _ lda)
    | lda == m  = V.elems v
    | otherwise = let
        breakCols [] = []
        breakCols es = let (c,es') = splitAt lda es in c:(breakCols es')

        dropJunk c = take m c

        in concatMap dropJunk $ breakCols (V.elems v)
{-# INLINE elems #-}

-- | Returns the contents of a matrix as a list of associations.
assocs :: (Storable e) => Matrix e -> [((Int,Int),e)]
assocs x = zip (indices x) (elems x)
{-# INLINE assocs #-}

-- | Version of 'update' that doesn't range-check indices.
unsafeUpdate :: (Storable e) => Matrix e -> [((Int,Int),e)] -> Matrix e
unsafeUpdate (Matrix v m n lda) ijes = let
    ies = [ (i + j * lda, e) | ((i,j),e) <- ijes ]
    v' = V.unsafeUpdate v ies
    in Matrix v' m n lda

-- | Create a new matrix by replacing the values at the specified indices.
update :: (Storable e) => Matrix e -> [((Int,Int),e)] -> Matrix e
update (Matrix v m n lda) ijes = let
    ies = [ if i < 0 || i >= m || j < 0 || j >= n
                then error $ printf
                         ("update"
                         ++ " <matrix with dim (%d,%d)>"
                         ++ " [ ..((%d,%d),_).. ]"
                         ++ ": invalid index")
                         m n i j
                else (i + j * lda, e)
          | ((i,j),e) <- ijes
          ]
    v' = V.unsafeUpdate v ies
    in Matrix v' m n lda

-- | Same as 'accum' but does not range-check indices.
unsafeAccum :: (Storable e)
            => (e -> e' -> e)
            -> Matrix e
            -> [((Int,Int), e')]
            -> Matrix e
unsafeAccum f (Matrix v m n lda) ijes = let
    ies = [ (i + j * lda, e) | ((i,j),e) <- ijes ]
    v' = V.unsafeAccum f v ies
    in Matrix v' m n lda

-- | @accum f@ takes a matrix and an association list and accumulates
-- pairs from the list into the matrix with the accumulating function @f@.
accum :: (Storable e)
      => (e -> e' -> e)
      -> Matrix e
      -> [((Int,Int), e')]
      -> Matrix e
accum f (Matrix v m n lda) ijes = let
    ies = [ if i < 0 || i >= m || j < 0 || j >= n
                then error $ printf
                         ("accum"
                         ++ " <matrix with dim (%d,%d)>"
                         ++ " [ ..((%d,%d),_).. ]"
                         ++ ": invalid index")
                         m n i j
                else (i + j * lda, e)
          | ((i,j),e) <- ijes
          ]
    v' = V.unsafeAccum f v ies
    in Matrix v' m n lda

-- | Construct a new matrix by applying a function to every element of
-- a matrix.
map :: (Storable e, Storable e')
    => (e -> e')
    -> Matrix e
    -> Matrix e'
map f a = fromList (dim a) $ P.map f (elems a)
{-# INLINE map #-}

-- | Construct a new matrix by applying a function to every pair of elements
-- of two matrices.  The two matrices must have identical dimensions.
zipWith :: (Storable e, Storable e', Storable f)
        => (e -> e' -> f)
        -> Matrix e
        -> Matrix e'
        -> Matrix f
zipWith f a a'
    | mn /= mn' = error $
        printf ("zipWith"
                ++ " <matrix with dim %s> "
                ++ " <matrix with dim %s>"
                ++ ": dimension mismatch"
                ) (show mn) (show mn')
    | otherwise =
        unsafeZipWith f a a'
  where
    mn  = dim a
    mn' = dim a'
{-# INLINE zipWith #-}

-- | Version of 'zipWith' that does not check if the input matrices
-- have the same dimensions.
unsafeZipWith :: (Storable e, Storable e', Storable f)
              => (e -> e' -> f)
              -> Matrix e
              -> Matrix e'
              -> Matrix f
unsafeZipWith f a a' =
    fromList (dim a') $ P.zipWith f (elems a) (elems a')
{-# INLINE unsafeZipWith #-}

-- | Get the given column of the matrix.
col :: (Storable e) => Matrix e -> Int -> Vector e
col a j
    | j < 0 || j >= n = error $
        printf ("col <matrix with dim (%d,%d)> %d:"
                ++ " index out of range") m n j
    | otherwise =
        unsafeCol a j
  where
    (m,n) = dim a
{-# INLINE col #-}

-- | Version of 'col' that doesn't range-check indices.
unsafeCol :: (Storable e) => Matrix e -> Int -> Vector e
unsafeCol (Matrix v m _ lda) j =
    V.unsafeSlice (j*lda) m v
{-# INLINE unsafeCol #-}

-- | Get a list of the columns of the matrix.
cols :: (Storable e) => Matrix e -> [Vector e]
cols a = P.map (unsafeCol a) [ 0..n-1 ]
  where
    (_,n) = dim a
{-# INLINE cols #-}

-- | @slice (i,j) (m,n) a@ creates a submatrix view of @a@ starting at
-- element @(i,j)@ and having dimensions @(m,n)@.
slice :: (Storable e)
      => (Int,Int)
      -> (Int,Int)
      -> Matrix e
      -> Matrix e
slice (i,j) (m',n') a
    | (i < 0 || m' < 0 || i + m' > m
       || j < 0 || n' < 0 || j + n' > n) = error $
        printf ( "slice"
               ++ " (%d,%d)"
               ++ " (%d,%d)"
               ++ " <matrix with dim (%d,%d)>"
               ++ ": index out of range"
               ) i j m' n' m n
    | otherwise =
        unsafeSlice (i,j) (m',n') a
  where
    (m,n) = dim a
{-# INLINE slice #-}

-- | Version of 'slice' that doesn't range-check indices.
unsafeSlice :: (Storable e)
            => (Int,Int) -> (Int,Int) -> Matrix e -> Matrix e
unsafeSlice (i,j) (m',n') (Matrix v _ _ lda) = let
    o = i + j*lda
    l = if m' == 0
            then 0
            else lda * n'
    v' = V.unsafeSlice o l v
    in Matrix v' m' n' lda
{-# INLINE unsafeSlice #-}

-- | Create a view of a matrix by taking the initial rows.
takeRows :: (Storable e) => Int -> Matrix e -> Matrix e
takeRows i a = slice (0,0) (i,n) a
  where
    (_,n) = dim a

-- | Create a view of a matrix by dropping the initial rows.
dropRows :: (Storable e) => Int -> Matrix e -> Matrix e
dropRows i a = slice (i,0) (m-i,n) a
  where
    (m,n) = dim a

-- | Split a matrix into two blocks and returns views into the blocks.  If
-- @(a1, a2) = splitRowsAt i a@, then
-- @a1 = slice (0,0) (i,n) a@ and
-- @a2 = slice (i,0) (m-i,n) a@, where @(m,n)@ is the dimension of @a@.
splitRowsAt :: (Storable e) => Int -> Matrix e -> (Matrix e, Matrix e)
splitRowsAt i a
    | i < 0 || i > m = error $
        printf ("splitRowsAt %d <matrix with dim (%d,%d)>:"
                ++ " invalid index") i m n
    | otherwise = let
        a1 = unsafeSlice (0,0) (i,n)   a
        a2 = unsafeSlice (i,0) (m-i,n) a
    in (a1,a2)
  where
    (m,n) = dim a
{-# INLINE splitRowsAt #-}

-- | Create a view of a matrix by taking the initial columns.
takeCols :: (Storable e) => Int -> Matrix e -> Matrix e
takeCols j a = slice (0,0) (m,j) a
  where
    (m,_) = dim a

-- | Create a view of a matrix by dropping the initial columns.
dropCols :: (Storable e) => Int -> Matrix e -> Matrix e
dropCols j a = slice (0,j) (m,n-j) a
  where
    (m,n) = dim a

-- | Split a matrix into two blocks and returns views into the blocks.  If
-- @(a1, a2) = splitColsAt j a@, then
-- @a1 = slice (0,0) (m,j) a@ and
-- @a2 = slice (0,j) (m,n-j) a@, where @(m,n)@ is the dimension of @a@.
splitColsAt :: (Storable e) => Int -> Matrix e -> (Matrix e, Matrix e)
splitColsAt j a
    | j < 0 || j > n = error $
        printf ("splitColsAt %d <matrix with dim (%d,%d)>:"
                ++ " invalid index") j m n
    | otherwise = let
        a1 = unsafeSlice (0,0) (m,j)   a
        a2 = unsafeSlice (0,j) (m,n-j) a
    in (a1,a2)
  where
    (m,n) = dim a
{-# INLINE splitColsAt #-}

-- | Convert a matrix to a vector by stacking its columns.
toVector :: (Storable e)
         => Matrix e
         -> Vector e
toVector a@(Matrix v m n _lda)
    | isContig a = v
    | otherwise  = V.fromList (m*n) $ elems a

-- | Cast a vector to a matrix of the given shape.
fromVector :: (Storable e)
           => (Int,Int)
           -> Vector e
           -> Matrix e
fromVector (m,n) v
    | nv /= m * n = error $
        printf ("fromVector"
               ++ " (%d,%d)"
               ++ " <vector with dim %d>"
               ++ ": dimension mismatch"
               ) m n nv
    | otherwise =
        Matrix v m n (max 1 m)
  where
    nv = V.dim v
{-# INLINE fromVector #-}

-- | Cast a vector to a matrix with one column.
fromCol :: (Storable e)
        => Vector e
        -> Matrix e
fromCol v = Matrix v m 1 (max 1 m)
  where
    m = V.dim v

-- | Cast a vector to a matrix with one row.
fromRow :: (Storable e)
        => Vector e
        -> Matrix e
fromRow v = Matrix v 1 n 1
  where
    n = V.dim v

instance (Storable e, Show e) => Show (Matrix e) where
    show x = "fromList " ++ show (dim x) ++ " " ++ show (elems x)
    {-# INLINE show #-}

instance (Storable e, Eq e) => Eq (Matrix e) where
    (==) = compareWith (==)
    {-# INLINE (==) #-}

instance (Storable e, AEq e) => AEq (Matrix e) where
    (===) = compareWith (===)
    {-# INLINE (===) #-}
    (~==) = compareWith (~==)
    {-# INLINE (~==) #-}

compareWith :: (Storable e, Storable e')
            => (e -> e' -> Bool)
            -> Matrix e
            -> Matrix e'
            -> Bool
compareWith cmp a a' =
    dim a == dim a'
    && and (P.zipWith cmp (elems a) (elems a'))
{-# INLINE compareWith #-}

-- | Create a matrix from a 'ForeignPtr' with offset, dimensions, and lda. The
-- data may not be modified through the ForeignPtr afterwards.
unsafeFromForeignPtr :: (Storable e)
                     => ForeignPtr e -- ^ pointer
                     -> Int	         -- ^ offset
                     -> (Int,Int)    -- ^ dimensions
                     -> Int          -- ^ leading dimension (lda)
                     -> Matrix e
unsafeFromForeignPtr p o (m,n) lda = let
    nv = if m == 0 then 0 else n * lda
    v = V.unsafeFromForeignPtr p o nv
    in Matrix v m n lda
{-# INLINE unsafeFromForeignPtr #-}

-- | Yield the underlying 'ForeignPtr' together with the offset to the data
-- the matrix dimensions, and the lda. The data may not be modified through
-- the 'ForeignPtr'.
unsafeToForeignPtr :: (Storable e)
                   => Matrix e
                   -> (ForeignPtr e, Int, (Int,Int), Int)
unsafeToForeignPtr (Matrix v m n lda) = let
    (f,o,_) = V.unsafeToForeignPtr v
    in (f, o, (m,n), lda)
{-# INLINE unsafeToForeignPtr #-}

-- | Execute an 'IO' action with a pointer to the first element in the
-- matrix and the leading dimension (lda).
unsafeWith :: (Storable e) => Matrix e -> (Ptr e -> Int -> IO a) -> IO a
unsafeWith (Matrix v _ _ lda) f =
    V.unsafeWith v $ \p -> f p lda
{-# INLINE unsafeWith #-}

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
module Types where

import qualified Data.Text as T
import Numeric ()
import Data.Complex
import Data.Ratio
import Data.Array
import Data.IORef
import System.IO
import Control.Monad.Except

import Text.ParserCombinators.Parsec (ParseError)

type Env = IORef [(T.Text, IORef LispVal)]

data LispVal
  = Atom T.Text
  | List [LispVal]
  | DottedList [LispVal] LispVal
  | Number Integer
  | String T.Text
  | Bool Bool
  | Char Char
  | Float Double
  | Ratio Rational
  | Complex (Complex Double)
  | Vector (Array Int LispVal)
  | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
  | IOFunc ([LispVal] -> IOThrowsError LispVal)
  | Port Handle
  | Func
      { params :: [T.Text]
      , varargs :: Maybe T.Text
      , body :: [LispVal]
      , closure :: Env
      }

instance Show LispVal where
  show = showVal

showVal :: LispVal -> String
showVal (Atom name) = T.unpack name
showVal (String contents) = "\"" ++ T.unpack contents ++ "\""
showVal (Number contents) = show contents
showVal (Float contents) = show contents
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (Char a) = [a]
showVal (List contents) = "(" ++ unwordsList contents ++ ")"
showVal (DottedList a as) = "(" ++ unwordsList a ++ " . " ++ showVal as ++ ")"
showVal (Ratio rat) = show (numerator rat) ++ "/" ++ show (denominator rat)
showVal (Complex (a :+ b)) = show a ++ "+" ++ show b ++ "i"
showVal (Vector arr) = "#(" ++ unwordsList (elems arr) ++ ")"
showVal (PrimitiveFunc _) = "<primitive>"
showVal (Port _)   = "<IO port>"
showVal (IOFunc _) = "<IO primitive>"
showVal Func {..} =
   "(lambda (" ++ unwords (map show params) ++
      (case varargs of
         Nothing -> ""
         Just arg -> " . " ++ T.unpack arg) ++ ") ...)"

unwordsList :: [LispVal] -> String
unwordsList = unwords . map showVal

instance Eq LispVal where
  (Atom a) == (Atom b) = a == b
  (List a) == (List b) = a == b
  (DottedList a c) == (DottedList b d) = a == b && c == d
  (Number a) == (Number b) = a == b
  (String a) == (String b) = a == b
  (Bool a) == (Bool b) = a == b
  (Char a) == (Char b) = a == b
  (Float a) == (Float b) = a == b
  (Ratio a) == (Ratio b) = a == b
  (Complex a) == (Complex b) = a == b
  (Vector a) == (Vector b) = a == b
  (Port a) == (Port b) = a == b
  -- TODO Add for other types
  _ == _ = False

data LispError
  = NumArgs Integer [LispVal]
  | VariableNumArgs [Integer] [LispVal]
  -- TODO This can be improved with something for min max num of params
  | TypeMismatch T.Text LispVal
  | Parser ParseError
  | BadSpecialForm T.Text LispVal
  | NotFunction T.Text T.Text
  | UnboundVar T.Text T.Text
  | OutOfBounds
  | DivideByZero
  | Default T.Text
  deriving (Eq)

showError :: LispError -> String
showError (UnboundVar message varname)  = T.unpack message ++ ": " ++ T.unpack varname
showError (BadSpecialForm message form) = T.unpack message ++ ": " ++ show form
showError (NotFunction message func)    = T.unpack message ++ ": " ++ show func
showError (NumArgs expected found)      = "Expected " ++ show expected
                                       ++ " args; found values " ++ unwordsList found
showError (VariableNumArgs expe found)  = "Expected " ++ show expe
                                       ++ " args; found values " ++ unwordsList found
showError (TypeMismatch expected found) = "Invalid type: expected " ++ T.unpack expected
                                       ++ ", found " ++ show found
showError (Parser parseErr)             = "Parse error at " ++ show parseErr
showError OutOfBounds                   = "Out of bounds error"
showError DivideByZero                  = "Division by 0 error"
showError (Default string)              = T.unpack string

instance Show LispError where
  show = showError

type ThrowsError = Either LispError
type IOThrowsError = ExceptT LispError IO

 -------------------------
-- Unpacking
-------------------------

unpackNum :: LispVal -> ThrowsError Integer
unpackNum (Number n) = return n
unpackNum notNum     = throwError $ TypeMismatch "number" notNum

unpackFloat :: LispVal -> ThrowsError Double
unpackFloat (Float n) = return n
unpackFloat notNum     = throwError $ TypeMismatch "float" notNum

unpackStr :: LispVal -> ThrowsError T.Text
unpackStr (String s) = return s
unpackStr notString  = throwError $ TypeMismatch "string" notString

unpackBool :: LispVal -> ThrowsError Bool
unpackBool (Bool b) = return b
unpackBool notBool  = throwError $ TypeMismatch "boolean" notBool

max_stars_repo_name: guimcaballero/scheme-hs

max_issues_repo_name: guimcaballero/scheme-hs

max_forks_repo_name: guimcaballero/scheme-hs


text:
module NeuralNet.ActivationSpec (activationSpec) where

import Test.Hspec
import NeuralNet.Activation
import Numeric.LinearAlgebra
import NeuralNet.Matrix

activationSpec :: SpecWith ()
activationSpec =
  describe "NeuralNet.Activation" $ do
    describe "forward" $ do
      it "should correctly work with ReLU" $
        let
          source = fromLists [ [1, 0, 0, 0]
                             , [0, 1, -1, 0]
                             , [0, 0, 1, 0]
                             , [0, 0, 0, 1]]
          expected = fromLists [ [1, 0, 0, 0]
                             , [0, 1, 0, 0]
                             , [0, 0, 1, 0]
                             , [0, 0, 0, 1]]
        in forward ReLU source `shouldBe` expected

      it "should correctly work with Sigmoid" $
        let
          source = zero 4 4
          expected = (4><4) (repeat 0.5)
        in forward Sigmoid source `shouldBe` expected

      it "should correctly work with Tanh" $
        let
          source = zero 4 4
          expected = (4><4) (repeat 0)
        in forward Tanh source `shouldBe` expected

    describe "backward" $ do
      it "should correctly work with ReLU eg 1" $
        let
          dA = fromLists [[-0.41675785, -0.05626683]]
          z = fromLists [[0.04153939, -1.11792545]]
          expected = fromLists [[-0.41675785, 0]]
        in backward ReLU dA z `shouldBe` expected

      it "should correctly work with ReLU eg 2" $
        let
          dA = fromLists [[0.12913162, -0.44014127], [-0.14175655, 0.48317296], [0.01663708, -0.05670698]]
          z = fromLists [[-0.7129932, 0.62524497], [-0.16051336, -0.76883635], [-0.23003072, 0.74505627]]
          expected = fromLists [[0, -0.44014127], [0, 0], [0, -0.05670698]]
        in backward ReLU dA z `shouldBe` expected

--      it "should correctly work with Tanh" $
--        let
--          source = diagonalList 3 0 [1..]
--          expected = fromLists [[1, 0, 0], [0, 4, 0], [0, 0, 9]]
--        in backward Tanh source `shouldBe` expected

      it "should correctly work with Sigmoid eg 1" $
        let
          dA = fromLists [[-0.41675785, -0.05626683]]
          z = fromLists [[0.04153939, -1.11792545]]
          expected = fromLists [[-0.1041445301481718, -1.0447915348382566e-2]]
        in backward Sigmoid dA z `shouldBe` expected

      it "should correctly work with Sigmoid eg 2" $
        let
          dA = fromLists [[-0.5590876, 1.77465392]]
          z = fromLists [[0.64667545, -0.35627076]]
          expected = fromLists [[-0.1261203972684239, 0.4298776138385045]]
        in backward Sigmoid dA z `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 Cas.Internal.Def.Interface
    ( module Cas.Internal.Def.Interface
    , module ReEx
    ) where

import           Data.List
import           Cas.Misc
import           PreludeCustom
import           Cas.Internal.Def as ReEx
import           Data.Complex.Generic

cPi, cE :: T
cPi = TC Pi
cE  = TC E

f :: F -> T
f  = TF

i = TF (imag 1) :: T

zeroF, oneF, aInvOneF :: F
(zeroF, oneF, aInvOneF ) = (0, 1, (-1))

zero, one, aInvOne :: T
(zero, one, aInvOne) = (f 0, f 1, f (-1))

x :: Int -> T
x = TX . X

_x :: Int -> X
_x = X

add, mul :: [T] -> T
add ts' = case length ts of
         0 -> zero
         1 -> head ts
         _ -> TAdd $ Add ts
    where
    ts = filter (not . isZero) ts'

mul ts'' = case length ts of
         0 -> one
         1 -> head ts
         _ -> TMul $ Mul ts
    where
    ts' = filter (not . isOne) ts''
    ts  = if elem zero ts'' then [zero] else ts'

addBin, mulBin,sub, pow :: T -> T -> T
addBin x y = add [x,y]
mulBin x y = mul [x,y]
sub x y = addBin x (aInv y)
pow = TPow

addRaw = TAdd . Add
mulRaw = TMul . Mul
ln :: T -> T
ln = TLn

aInv :: T -> T
aInv =  mul . (:[aInvOne])

mInv :: T -> T
mInv = flip pow aInvOne

tAddTs, tMulTs :: T -> [T]

tXI    = xI    . tX
tAddTs = addTs . tAdd
tMulTs = mulTs . tMul

isZeroF, isOneF :: F -> Bool
isZero, isOne, isC, isF, isX, isAdd, isMul, isLeaf, isFlat :: T   -> Bool
isThisX :: X -> T -> Bool

isZero    = (==zero)
isZeroF   = (==zeroF)
isOne     = (==one)
isOneF    = (==oneF)
isC       = \case (TC _)   -> True; _ -> False
isF       = \case (TF _)   -> True; _ -> False
isX       = \case (TX _)   -> True; _ -> False
isThisX x = \case (TX y)   -> x==y; _ -> False
isAdd     = \case (TAdd _) -> True; _ -> False
isMul     = \case (TMul _) -> True; _ -> False
isPow     = \case (TPow _ _) -> True; _ -> False
isLn      = \case (TLn _)  -> True; _ -> False

isLeaf     = or . ([isF, isC, isX,isLn] <*>) . return
isNonFLeaf = or . ([isX, isC, isLn] <*>) . return
isOper     = or . ([isAdd, isMul, isPow] <*>) . return

isFlat     = isAdd ?>>> all (not . isAdd) . tAddTs
         ||> isMul ?>>- all (not . isMul) . tMulTs
                    ||> True

simpleShow :: T -> String
simpleShow TC {..}      = show $  tC
simpleShow t@TF {..}    = wrap $ showF  tF
simpleShow t@TX {..}    = "x" <> show ( tXI t)
simpleShow t@TAdd {..}  = wrap $ intercalate "+" $ fmap simpleShow (tAddTs t)
simpleShow t@TMul {..}  = wrap $ intercalate "*" $ fmap simpleShow (tMulTs t)
simpleShow t@TPow {..}  = wrap $ simpleShow tPowB <> "^" <> simpleShow tPowE
simpleShow t@TLn {..}   = wrap $ "ln " <> simpleShow tLnT

(+:+) z w = TF $ z :+ w

eqType :: T -> T -> Bool
eqType (TF _) (TF _) = True
eqType (TX _) (TX _) = True
eqType (TC _) (TC _) = True
eqType (TAdd _) (TAdd _) = True
eqType (TMul _) (TMul _) = True
eqType (TPow _ _) (TPow _ _) = True
eqType (TLn _) (TLn _) = True
eqType t r = False

max_stars_repo_name: vcanadi/cas

max_issues_repo_name: vcanadi/cas

max_forks_repo_name: vcanadi/cas


text:
module Vectorspace where

import Data.Complex
import Data.Maybe

newtype Vector a = Vector [a] deriving (Eq)
newtype Matrix a = Matrix [[a]] deriving (Eq)

instance Show a => Show (Vector a) where
        show (Vector v)
            | length v == 0 = "()"
            | length v == 1 = "(" ++ show (v!!0) ++ ")"
            | length v == 2 = start v ++ end v
            | length v >= 3 = start v ++ center v ++ end v where
                start v = vectortop $ (show.head) v
                center v = concatMap (vectorcenter.show) ( (init.tail) v)
                end   v = vectorbottom $ (show.last) v
                vectortop x = "/" ++ x ++ "\\\n"
                vectorcenter x = "|" ++ x ++ "|\n"
                vectorbottom x = "\\" ++ x ++ "/"

instance Show a => Show (Matrix a) where
        show (Matrix m) = start m ++ "\n" ++ center m ++ "\n" ++ end m where
            start m  = (matrixtop.concat) [ cell v | v <- head m ]
            center m = (matrixcenter.concat) [ cell v | v <- (head.init.tail) m ]
            end   m  = (matrixbottom.concat) [ cell v | v <- last m ]
            matrixtop v    = "┌" ++ v ++ "┐"
            matrixcenter v = "│" ++ v ++ "│"
            matrixbottom v = "└" ++ v ++ "┘"
            cell v = [ if i< length (show v) then (show v)!!i else ' ' | i <- [0..cellWidth] ]
            cellWidth =  maximum $ map (length.show) $ concat m

instance Floating a => Num (Vector a) where
        -- componentwise addition
        (Vector v) + (Vector w)
                            | length v == length w = Vector $ zipWith (+) v w
                            | length v >  length w = Vector v + Vector (take (length v) (w++repeat 0))
                            | length v <  length w = Vector w + Vector v
        (Vector v) * (Vector w) = Vector $ zipWith (*) v w
        negate (Vector v) = Vector $ map (*(-1)) v
        abs (Vector v) = Vector [sqrt $ sum (map (**2) v)]
        signum (Vector v) = Vector $ map signum v
        fromInteger n = Vector [fromInteger n]

-- multiplies a vector by a scalar
vecmul :: Integer -> [Integer] -> [Integer]
vecmul a = map (*a)

-- multiplies a matrix by a scalar
matscalmult :: Num a => [[a]] -> a -> [[a]]
matscalmult m s = map (\y -> map (\x -> (m!!y!!x)*s) [0..length (m!!y)-1]) [0..length m-1]

-- multiplies / composes two matrices
matmult :: Num a => [[a]] -> [[a]] -> [[a]]
matmult m n
    | length m > length n         = matmult m (scale n (length m))
    | length m /= length (head m) = matmult (sqrify m) n
    | length n /= length (head n) = matmult m (sqrify n)
    | otherwise                   = [[
                                          sum (zipWith (*) (getColMat m y) (getRowMat n x))
                                      | x<-[0..(length (m!!y)-1)] ]
                                  | y<-[0..(length m-1)]]

-- multiplies a matrix by a vector / applies a matrix to a vector
matvecmult :: Num a => [[a]] -> [a] -> [a]
matvecmult m v = map head $ matmult m (map (\e-> [e]) v)

applyMatrix :: Num a => [[a]] -> [a] -> [a]
applyMatrix m v = map (sum . zipWith (*) v) m

-- adds two matrices together
matadd :: Num a => [[a]] -> [[a]] -> [[a]]
matadd m n
    | length m > length n       = matadd m (scale n (length m))
    | length m /= length (m!!0) = matadd (sqrify m) n
    | length n /= length (n!!0) = matadd m (sqrify n)
    | otherwise                 = [[
                                         (m!!y!!x) + (n!!y!!x)
                                    | x <- [0..((length(m!!0))-1)]]
                                | y <- [0..((length m)-1)]]

-- subtracts one matrix from another
matsub :: Num a => [[a]] -> [[a]] -> [[a]]
matsub m n = matadd m (matscalmult n (-1))

-- creates an identity matrix of dimension n
matid :: Num a => Int -> [[a]]
matid n = matgen n (\ i j -> if i==j then 1 else 0 )

-- maps a matrices entries based on a function that takes the indices as input
matmap :: Num a => [[a]] -> ( Int -> Int -> a -> a ) -> [[a]]
matmap m f = map (\y -> map (\x -> f y x (m!!y!!x) ) [0..(length (m!!y)-1)]) [0..( length m -1)]

-- complex conjugat of a matrix
matcon :: RealFloat a => [[Complex a]] -> [[Complex a]]
matcon m = matmap m (\y x a -> conjugate a)

-- generates a matrix of size t and then maps it with values
matgen :: (Num a, Num t, Enum t) => t -> (t -> t -> a) -> [[a]]
matgen s f = [[
             f x y
        | x<- [0..(s-1)]]
    | y<- [0..(s-1)]]

-- fills in entries with 0's, so that a matrix is square
sqrify :: Num a => [[a]] -> [[a]]
sqrify m = scale m (max (length m) (length (head m)))

-- scales a matrix m to a size s (also makes it a square matrix
scale :: Num a => [[a]] -> Int -> [[a]]
scale m s = [[ if (y<(length m) && x<(length (m!!0))) then(m!!y!!x) else(0) | x <- [0..(s-1)] ] | y <- [0..(s-1)] ]

-- pretty prints a matrix
printMat :: (Foldable t, Show a) => t a -> IO ()
printMat m = mapM_ print m

-- all non zero elements of a matrix are replaces by '.', otherwise they are replaced by ' '
prettifyMat m = [[ if(m!!y!!x == 0)then(' ')else('.') | x<-[0..(length (m!!y)-1)] ] | y<-[0..((length m)-1)]]

indexMatrix :: (Num a, Num b,Enum b)  => [[a]] -> [([(a,b)], b)]
indexMatrix m = (map (\i -> zip i [0..]) m) `zip` [0..]

-- gets the element of the y th row and x th collumn of a matrix
matxy :: Num a => [[a]] -> Int -> Int -> a
matxy m x y = m!!x!!y

-- transposes a matrix
transpose :: Num a => [[a]] -> [[a]]
transpose m = do
    let im = indexMatrix m
    map (\i -> map (\j -> matxy m (snd j) (snd i)) (fst i)) im

-- transposes a vector
vectrans :: Num a => [a] -> [[a]]
vectrans v = [ [x] | x<-v ]

-- adjugate matrix of m
adju :: Num a => [[a]] -> [[a]]
adju m = transpose (cof m)

-- cut out the matrix from x1-x2 and y1-y2 of a matrix m
submatrix :: Num a => [[a]] -> Int -> Int -> Int -> Int -> [[a]]
submatrix m x1 y1 x2 y2 = do
    let ycut = take (y2-y1) (drop y1 m)
    map (\i -> take (x2-x1) (drop x1 i)) ycut

-- gets the elements of a matrix in row r
getRowMat :: Num a => [[a]] -> Int -> [a]
getRowMat m r = [ m!!r!!x | x <- [0..((length (m!!0))-1)]]

-- gets the elements of a matrix in collumn c
getColMat :: Num a => [[a]] -> Int -> [a]
getColMat m c = [ m!!x!!c | x <- [0..((length m)-1)]]

-- return the matrix with collumn c removed
elimRowMat :: Num a => [[a]] -> Int -> [[a]]
elimRowMat m r = map (\i -> take r i ++ drop (r+1) i) m

-- returns the matrix with row r removed
elimColMat :: Num a => [[a]] -> Int -> [[a]]
elimColMat m c = take c m ++ drop (c+1) m

-- eliminates row y and collumn x
elimRowColMat :: (Num a) => [[a]] -> Int -> Int -> [[a]]
elimRowColMat m x y = elimRowMat (elimColMat m y) x

-- takes a matrix and replaces the n th collumn with a vector
replCol :: Num a => [[a]] -> Int -> [a] -> [[a]]
replCol m i v = map (\y -> map (\x -> if(x==i && y<length v)then(v!!y)else(m!!y!!x)) [0..((length (m!!0))-1)]) [0..((length m)-1)]

-- takes a matrix and replaces the n th row with a vector
replRow :: Num a => [[a]] -> Int -> [a] -> [[a]]
replRow m i v = map (\y -> map (\x -> if(y==i && x<length v)then(v!!x)else(m!!y!!x)) [0..((length (m!!0))-1)]) [0..((length m)-1)]

-- calculates the determinant of a matrix
det :: Num a => [[a]] -> a
det m = case (length m, length (m!!0)) of
    (1,1) -> m!!0!!0
    (p,q) -> do
        let grid = [ (m!!0!!x)*(-1)^x | x<-[0..p-1] ]
        let underdet = [ (grid!!x) * det(elimRowColMat m x 0) | x<-[0..p-1]]
        sum underdet

-- returns the minor of a matrix from row x and col y
minor :: Num a => [[a]] -> Int -> Int -> a
minor m x y = det (elimRowColMat m x y)

-- cofactor matrix / i.e. replace each entry with their minor and apply the checkerboard pattern of (1) and (-1)
cof :: Num a => [[a]] -> [[a]]
cof m = [
        [
            ((-1)^(x+y)) * (minor m x y)
        | x<-[0..(length (m!!y))-1]]
    | y<-[0..(length m)-1]]

-- returns if a matrix is invertible or nor
invertible :: (Eq a, Num a) => [[a]]-> Bool
invertible m = (det m) /= 0

-- returns the inverse of a matrix, if it exists
inv :: (Eq a, Fractional a) => [[a]] -> Maybe [[a]]
inv m
    | not (invertible m) = Nothing
    | otherwise = Just (matscalmult (adju m) (1/det(m)))

-- removes duplicates inside a list
remDupl :: (Eq a) => [a] -> [a]
remDupl a
    | a == [] = []
    | not (elem (a!!0) (drop 1 a)) = (a!!0):(remDupl(drop 1 a))
    | otherwise = remDupl (drop 1 a)

-- true if list a and list b are linearily dependent
linDep :: (Fractional a, Eq a, Enum a) => [a] -> [a] -> Bool
linDep a b
    | length a < length b   = linDep b a
    | length a == 0         = True
    | a!!0==b!!0 && a!!0==0  = linDep (tail a) (tail b)
    | a!!0/=0 && b!!0/=0    = all (\e -> (fst e)/b!!0 == (snd e)/a!!0) $ zip a (b++[0,0..])
    | otherwise             = False

-- adds a collumn of 0 at the beginning
addCol :: Num a => [[a]] -> [[a]]
addCol m = map (\e -> 0:e) m

-- reduced row-echelon form (TODO)
rref :: (Fractional a, Eq a) => [[a]] -> [[a]]
rref m = rref_step 0 m

rref_step :: (Fractional a, Eq a) => Int -> [[a]] -> [[a]]
rref_step s m
    | coef s == 0                               = rref_step (s+1) m
    | length (newPivot (normalize s) s) == 0    = normalize s
    | s+1 < length m                            = rref_step (s+1) (elim s)
    | otherwise                                 = normalize s
    where
        coef s          = length ( filter (\i-> i/=0) (getRowMat m s) )
        ones s          = length ( filter (\i-> i==1) (getRowMat m s) )
        emptyUntil s r  = 0 == length ( filter (\i-> i/=0) (take s (getColMat m r)))
        normalize s     = map (\y-> if(m!!y!!s /= 0 && emptyUntil s y)then(map (/m!!y!!s) (m!!y))else(m!!y) ) [0..((length m)-1)]
        newPivot n s    = (take 1 (filter (\y-> (emptyUntil s y) && (n!!y!!s==1)) [0..((length n)-1)]))
        elim s          = do
                            let n = normalize s
                            let np = (newPivot n s)!!0
                            map (\y-> if(n!!y!!s /=0 && y/=np)then( map (\x-> (n!!y!!x)- ((n!!y!!s)*(n!!np!!x))) [0..((length (n!!0))-1)] )else( n!!y )) [0..((length n)-1)]

dimKer :: (Fractional a, Eq a) => [[a]] -> Int
dimKer m = length $ filter (\e-> 0 == (length ( filter (\i-> i/=0) e)) ) (rref m)

--rank :: (Fractional a, Eq a) => [[a]] -> Int
rank m = length m - dimKer m

ker :: (Fractional a, Eq a) => [[a]] -> Maybe [[a]]
ker m
    | det m /= 0 = Nothing
    | otherwise = do
        let r = rref m
        let f = filter (\i-> any (\j-> j/=0) i) r
        Just f

solvMat :: Fractional a => [[a]] -> [a] -> [a]
solvMat m v = map (\x-> det (replCol m x v) / det m) [0..((length (m!!0))-1)]

--jord :: Num a => [[a]] -> [[a]]
--jord m =

--charpoly :: Num a => [[a]] -> a -> a
--charpoly m x = det (matadd m (matscalmult (matid (length m)) x ))
--
--eigenVal :: (RealFrac a, Enum a, Integral b) => [[a]] -> [b]
--eigenVal m = filter (\i-> i/=0) $ remDupl $ map (round) $ filter (\x -> round (charpoly m x)>(-10) && round (charpoly m x)<10) [-20,-19.9 .. 20]
--
--geoMult :: (Fractional a, Eq a) => [[a]] -> a -> Int
--geoMult m e = dimKer $ matadd (matscalmult (matid (length (sqrify m))) e) (matscalmult m (-1))

tr :: Num a => [[a]] -> a
tr m = sum $ map (\i -> (sqrify m)!!i!!i) [0..((length (m!!0))-1)]

polyadd :: Num a => [a] -> [a] -> [a]
polyadd p q
    | length p < length q = polyadd q p
    | otherwise = [ p!!x + if(x<length q)then(q!!x)else(0) | x<- [0..((length p)-1)]]

polymonomult :: Num a => [a] -> (a,Int) -> [a]
polymonomult p m = [0| x<-[1..(snd m)]] ++ map (* (fst m)) p

--polymult :: Num a => [a] -> [a] -> [a]
polymult p q = finish $ map (\i -> polymonomult q (p!!i,i) ) [0..((length p)-1)]
    where
        finish x
            | length x == 1 = x
            | otherwise = finish $ polyadd (x!!0) (x!!1) : tail (tail x)

--jnf :: Num a => [[a,Int,Int]] => [[a]]

--algMult ::
--algMult m es =
--    where s = length (sqrify m)
--
--partitions :: Int -> [Int]
--partitions x =

dotprod :: Num a => [a] -> [a] -> a
dotprod v w = sum $ matvecmult (vectrans v) w

orthoVec :: (Num a, Eq a) => [a] -> [a] -> Bool
orthoVec v w = (dotprod v w) == 0

adjo :: RealFloat a => [[Complex a]] -> [[Complex a]]
adjo m = transpose $ matcon m

ortho :: (Num a, Eq a) => [[a]] -> Bool
ortho m = matid (length (sqrify m)) == matmult m (transpose m)

hermit :: RealFloat a => [[Complex a]] -> Bool
hermit m = adjo m == m

unitar :: RealFloat a => [[Complex a]] -> Bool
unitar m = adjo m == fromMaybe [[0]] (inv m)

kroneck :: Num a => [[a]] -> [[a]] -> [[a]]
kroneck m n = [
            [ (m!!( yPosM y )!!( xPosM x )) * (n!!( yPosN y )!!( xPosN x ))
        | x<-[0..((length (m!!0))*(length (n!!0))-1)]]
    | y<-[0..((length m)*(length n)-1)]]
    where
        yPosM y = y `mod` (length m)
        xPosM x = x `mod` (length (m!!0))
        yPosN y = if (y==0)then(0)else( y `div` (length m) )
        xPosN x = if (x==0)then(0)else( x `div` (length (m!!0)) )

max_stars_repo_name: Quoteme/vectorspace

max_issues_repo_name: Quoteme/vectorspace

max_forks_repo_name: Quoteme/vectorspace


text:
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}

module Grenade.Layers.Internal.BLAS
    ( -- easy interface
      matXVec
    , outerV
    , checkVectors
    , unsafeMemCopyVectorFromTo
    , memCopyVectorFromTo
    , unsafeMemZero
    , memZero
      -- more complicated, but direct, function calls
    , BlasTranspose (..)
    , swapTranspose
    , dgemmUnsafe
    , dgemvUnsafe
    , dgerUnsafe
    ) where

import           Control.Monad
import           Data.IORef
import           Data.Proxy
import qualified Data.Vector.Storable         as V
import           Foreign                      (withForeignPtr)
import           Foreign.C.Types
import           Foreign.Ptr
import           Foreign.Storable             (sizeOf)
import           GHC.IO.Handle.Text           (memcpy)
import           GHC.TypeLits
import qualified Numeric.LinearAlgebra        as LA
import qualified Numeric.LinearAlgebra.Static as LAS
import           System.IO.Unsafe             (unsafePerformIO)

import           Grenade.Layers.Internal.CUDA
import           Grenade.Types
import           Grenade.Utils.Vector

import           Debug.Trace

#define USE_DGEMM_ONLY 0

-- | Computes vec2 <- mat * vec1 + beta * vec2.
matXVec :: BlasTranspose -> V.Vector RealNum -> V.Vector RealNum -> RealNum -> V.Vector RealNum -> IO (V.Vector RealNum)
matXVec trMat mat vec1 beta vec2 =
#if USE_DGEMM_ONLY
  dgemmUnsafe trMat BlasNoTranspose (m, k) (ay, 1) 1.0 mat vec1 beta vec2
# else
  dgemvUnsafe trMat (m, k) 1.0 mat vec1 beta vec2
#endif
  where
    ay = V.length vec1
    ax = V.length vec2
    (m, k) = swapTranspose trMat (ax, ay)

-- | Computes the outer product of two vectors: mat <- vec1 `outer` vec2
outerV :: V.Vector RealNum -> V.Vector RealNum -> IO (V.Vector RealNum)
outerV vec1 vec2 =
#if USE_DGEMM_ONLY
  dgemmUnsafe BlasNoTranspose BlasNoTranspose (o, 1) (1, i) 1.0 vec1 vec2 0 (createVectorUnsafe (i * o)) -- beta = 0 initialises the matrix
#else
  createVector (i * o) >>= memZero >>= dgerUnsafe (o, i) 1.0 vec1 vec2
#endif
  where o = V.length vec1
        i = V.length vec2

-- | Check two vectors if they are equal. For testing purposes.
checkVectors :: V.Vector RealNum -> V.Vector RealNum -> Bool
checkVectors v1 v2 = V.length v1 == V.length v2 && and (zipWith (==) (toStr v1) (toStr v2))
  where
    toStr :: V.Vector RealNum -> [String]
    toStr v = map (show . round . (*10^5)) $ V.toList v
{-# INLINE checkVectors #-}

-- | Newtype holding CINT for Transpose values.
newtype CBLAS_TRANSPOSET =
  CBLAS_TransposeT CInt
  deriving (Eq, Show)

-- | Transpose values
data BlasTranspose
  = BlasNoTranspose
  | BlasTranspose
  | BlasConjTranspose
  | BlasConjNoTranspose
  deriving (Eq, Show)

encodeTransposeIntBool :: BlasTranspose -> Int
encodeTransposeIntBool BlasNoTranspose     = 0
encodeTransposeIntBool BlasTranspose       = 1
encodeTransposeIntBool BlasConjTranspose   = 1
encodeTransposeIntBool BlasConjNoTranspose = 0
{-# INLINE encodeTransposeIntBool #-}

swapTranspose :: BlasTranspose -> (Int, Int) -> (Int, Int)
swapTranspose BlasNoTranspose x        = x
swapTranspose BlasTranspose (a, b)     = (b, a)
swapTranspose BlasConjNoTranspose x    = x
swapTranspose BlasConjTranspose (a, b) = (b, a)
{-# INLINE swapTranspose #-}

-- | Error text
mkDimText :: (Show a1, Show a2, Show a3, Show a4, Show a5, Show a6) => (a1, a2) -> (a3, a4) -> (a5, a6) -> String
mkDimText (ax, ay) (bx, by) (cx, cy) = "resulting dimensions: [" ++ show ax ++ "x" ++ show ay ++ "]*[" ++ show bx ++ "x" ++ show by ++ "]=[" ++ show cx ++ "x" ++ show cy ++ "]"

-- | Computes: C <- alpha*op( A )*op( B ) + beta*C, where op(X) may transpose the matrix X
--
-- dgemm, see http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gaeda3cbd99c8fb834a60a6412878226e1.html for the documentation.
--
-- void cblas_dgemm (
--              const CBLAS_LAYOUT      layout,
--              const CBLAS_TRANSPOSE   TransA,
--              const CBLAS_TRANSPOSE   TransB,
--              const int       M,
--              const int       N,
--              const int       K,
--              const double    alpha,
--              const double *          A,
--              const int       lda,
--              const double *          B,
--              const int       ldb,
--              const double    beta,
--              double *        C,
--              const int       ldc
-- 	)
{-# NOINLINE dgemmUnsafe #-}
dgemmUnsafe :: BlasTranspose    -- ^ Transpose Matrix A
            -> BlasTranspose    -- ^ Transpose Matrix B
            -> (Int, Int)       -- ^ Rows and cols of A on entry (not transposed)
            -> (Int, Int)       -- ^ Rows and Cols of B on entry (not transposed)
            -> RealNum           -- ^ Alpha
            -> V.Vector RealNum  -- ^ A
            -> V.Vector RealNum  -- ^ B
            -> RealNum           -- ^ Beta
            -> V.Vector RealNum  -- ^ C
            -> IO (V.Vector RealNum)  -- ^ Return new C
dgemmUnsafe trA trB (axIn, ayIn) (bxIn, byIn) alpha matrixA matrixB beta matrixC
  | isBadGemm =
    error $!
    "bad dimension args to dgemmUnsafe: ax ay bx by cx cy: " ++
    show [ax, ay, bx, by, ax, by] ++ " matrix C length: " ++ show (V.length matrixC) ++ "\n\t" ++ mkDimText (ax, ay) (bx, by) (ax, by)
  | otherwise = do
      V.unsafeWith matrixA $ \aPtr' ->
        V.unsafeWith matrixB $ \bPtr' ->
          V.unsafeWith matrixC $ \cPtr' ->  do
#ifdef USE_FLOAT
            sgemm_direct
#else
            dgemm_direct
#endif
              (encodeTransposeIntBool trA) -- transpose A
              (encodeTransposeIntBool trB) -- transpose B
              (fromIntegral ax)     -- rows of C = rows of A transposed
              (fromIntegral by)     -- cols of C = cols of B transposed
              (fromIntegral ay)     -- k = cols of A transposed = rows of B transposed
              alpha
              aPtr'
              (fromIntegral axIn) -- LDA
              bPtr'
              (fromIntegral bxIn) -- LDB
              beta
              cPtr'
              (fromIntegral ax)   -- LDC
            return matrixC
  where
    (ax, ay) = swapTranspose trA (axIn, ayIn)
    (bx, by) = swapTranspose trB (bxIn, byIn)
    isBadGemm = minimum [ax, ay, bx, by] <= 0 || not (ax * by == V.length matrixC && ay == bx)

-- | Computes: Y <- alpha*op( A )*X + beta*Y, where op(A) may transpose the matrix A
--
-- dgemv, see http://www.netlib.org/lapack/explore-html/d7/d15/group__double__blas__level2_gadd421a107a488d524859b4a64c1901a9.html#gadd421a107a488d524859b4a64c1901a9
--
-- void cblas_dgemv 	(
--              const CBLAS_LAYOUT  	layout,
-- 		const CBLAS_TRANSPOSE  	TransA,
-- 		const int  	M,
-- 		const int  	N,
-- 		const double  	alpha,
-- 		const double *  	A,
-- 		const int  	lda,
-- 		const double *  	X,
-- 		const int  	incX,
-- 		const double  	beta,
-- 		double *  	Y,
-- 		const int  	incY
-- 	)
{-# NOINLINE dgemvUnsafe #-}
dgemvUnsafe :: BlasTranspose    -- ^ Transpose Matrix
            -> (Int, Int)       -- ^ rows and cols of A on entry (not transposed)
            -> RealNum           -- ^ Alpha
            -> V.Vector RealNum  -- ^ A
            -> V.Vector RealNum  -- ^ X
            -> RealNum           -- ^ Beta
            -> V.Vector RealNum  -- ^ C
            -> IO (V.Vector RealNum)  -- ^ Return new C
dgemvUnsafe trA (m, k) alpha matrixA vecX beta vecY
  | ax /= V.length vecY || ay /= V.length vecX =
    error $!
    "bad dimension args to dgemvUnsafe: ax ay (length vecX) (length vecY): " ++
    show [ax, ay, V.length vecX, V.length vecY] ++ " \n\t" ++ mkDimText (ax, ay) (V.length vecX, 1) (m, 1)
  | otherwise = do
      V.unsafeWith matrixA $ \aPtr' ->
        V.unsafeWith vecX $ \xPtr' ->
          V.unsafeWith vecY $ \yPtr' -> do
#ifdef USE_FLOAT
            sgemv_direct
#else
            dgemv_direct
#endif
              (encodeTransposeIntBool trA)      -- transpose A
              (fromIntegral m)
              (fromIntegral k)
              alpha
              aPtr'
              (fromIntegral m)
              xPtr'
              1
              beta
              yPtr'
              1
            return vecY
  where
    (ax, ay) = swapTranspose trA (m, k)

-- | Computes: A <- alpha*X*Y^T + A
--
-- dger, see http://www.netlib.org/lapack/explore-html/d7/d15/group__double__blas__level2_ga458222e01b4d348e9b52b9343d52f828.html#ga458222e01b4d348e9b52b9343d52f828
--
-- void cblas_dger 	(
--              const CBLAS_LAYOUT  	layout,
-- 		const int  	M,
-- 		const int  	N,
-- 		const double  	alpha,
-- 		const double *  	X,
-- 		const int  	incX,
-- 		const double *  	Y,
-- 		const int  	incY,
-- 		double *  	A,
-- 		const int  	lda
-- 	)
{-# NOINLINE dgerUnsafe #-}
dgerUnsafe :: (Int, Int)           -- ^ Dimensions of matrix A
           -> RealNum               -- ^ Alpha
           -> V.Vector RealNum      -- ^ X
           -> V.Vector RealNum      -- ^ C
           -> V.Vector RealNum      -- ^ A
           -> IO (V.Vector RealNum)  -- ^ Return new C
dgerUnsafe (ax, ay) alpha vecX vecY matrixA
  | ax /= len || ay /= V.length vecY =
    error $! "bad dimension args to dgerUnsafe: X Y ax ay: " ++ show [len, V.length vecY, ax, ay] ++ " \n\t" ++ mkDimText (ax, ay) (V.length vecX, 1) (V.length vecY, 1)
  | otherwise = do
      V.unsafeWith matrixA $ \aPtr' ->
        V.unsafeWith vecX $ \xPtr' ->
          V.unsafeWith vecY $ \yPtr' -> do
#ifdef USE_FLOAT
            sger_direct
#else
            dger_direct
#endif
              (fromIntegral ax)
              (fromIntegral ay)
              alpha
              xPtr'
              1
              yPtr'
              1
              aPtr'
              (fromIntegral ax)
            return matrixA
  where
    len = V.length vecX

-- |  Matrix mult for general dense matrices
type BLASGemmFunFFI scale el
  =  Int -- transpose A: 1, not transpose A: 0
  -> Int -- transpose B: 1, not transpose B: 0
  -> CInt -- m
  -> CInt -- n
  -> CInt -- k
  -> {- scal A * B -} scale  -- alpha
  -> {- Matrix A-} Ptr el    -- A
  -> CInt                    -- LDA
  -> {- B -} Ptr el
  -> CInt
  -> scale                   -- beta
  -> {- C -}  Ptr el
  -> CInt
  -> IO ()

foreign import ccall unsafe "dgemm_direct" dgemm_direct :: BLASGemmFunFFI Double Double
foreign import ccall unsafe "sgemm_direct" sgemm_direct :: BLASGemmFunFFI Float Float

-- |  Matrix mult for general dense matrices
type BLASGemvFunFFI scale el
  =  Int    -- transpose A: 1, not transpose A: 0
  -> CInt   -- m
  -> CInt   -- n
  -> scale  -- alpha
  -> Ptr el -- Matrix A
  -> CInt   -- LDA
  -> Ptr el
  -> CInt
  -> scale -- beta
  -> Ptr el
  -> CInt
  -> IO ()

foreign import ccall unsafe "dgemv_direct" dgemv_direct :: BLASGemvFunFFI Double Double
foreign import ccall unsafe "sgemv_direct" sgemv_direct :: BLASGemvFunFFI Float Float

type BlasGerxFunFFI scale el
  =  CInt
  -> CInt
  -> scale
  -> Ptr el
  -> CInt
  -> Ptr el
  -> CInt
  -> Ptr el
  -> CInt
  -> IO ()

foreign import ccall unsafe "dger_direct" dger_direct :: BlasGerxFunFFI Double Double
foreign import ccall unsafe "sger_direct" sger_direct :: BlasGerxFunFFI Float Float

-- toRows :: Int -> V.Vector Double -> [V.Vector Double]
-- toRows m vec = LA.toRows . reshapeF m . LA.vector . V.toList $ vec
--   where reshapeF r = LA.tr' . LA.reshape r

-- vec1 :: LAS.R 10
-- vec1 = LAS.vector [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.34868882831497655, 0.0, 1.4026932193043212e-2]

-- -- dEdy:
-- -- mmCheck : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-4.8711473950841355e-2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-1.9595481318788644e-3]
-- -- mm' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-0.24602759633121374,-1.2277801012906878e-2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-9.89711206989971e-3,-4.939070836308897e-4,0.0,0.0]

-- vec3 :: V.Vector Double
-- vec3 = (
--   V.fromList [0.0, 0.0, 0.0, 0.0, (-0.13969898085418284)])

-- vec2 :: LAS.R 5
-- vec2 = LAS.vector [0.0, 0.0, 0.0, 0.0, -0.13969898085418284]

-- res = vec1 `LAS.outer` vec2

-- test =
--   toRows 10 $
--   outerV (LAS.extract vec1) (LAS.extract vec2) (V.replicate (LAS.size vec1 * LAS.size vec2) 10)

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
module DFTSpec (
  dftSpec
)
where
import Test.Hspec
import Test.QuickCheck
import Control.Exception (evaluate)
import DFT (dft, idft)
import Data.Complex
import Data.List (zip, repeat)

dftSpec :: Spec
dftSpec = do
  describe "dft and inverse dft tests..." $ do
    it "checks whether composition of inverse dft and dft returns the original sequence" $
      property $ prop_dft_correct

prop_dft_correct :: [Float] -> Bool
prop_dft_correct xs = are_numbers_nearly_equal ((idft . dft . map (\x -> x :+ 0) ) ys) ys
  where ys = preprocess xs

preprocess :: [Float] -> [Float]
preprocess [] = []
preprocess xs = take targetLength (xs ++ (repeat 0))
  where targetLength = 2 ^ ceiling (log ((fromIntegral . length) xs) / log 2)

are_numbers_nearly_equal :: [Complex Float] -> [Float] -> Bool
are_numbers_nearly_equal [] [] = True
are_numbers_nearly_equal (x:xs) (y:ys)
  | realPartsEq = are_numbers_nearly_equal xs ys
  | otherwise   = False
  where realPartsEq = nearly_equal (realPart x) y

nearly_equal :: Float -> Float -> Bool
nearly_equal x y = abs (x - y) <= 1e-4

max_stars_repo_name: gboduljak/dft

max_issues_repo_name: gboduljak/dft

max_forks_repo_name: gboduljak/dft


text:
{-# LANGUAGE CPP                     #-}
{-# LANGUAGE FlexibleContexts        #-}
{-# LANGUAGE TypeApplications        #-}
{-# LANGUAGE ScopedTypeVariables     #-}
{-# LANGUAGE ConstraintKinds         #-}
{-# LANGUAGE TypeFamilies            #-}
{-# LANGUAGE TypeOperators           #-}
{-# LANGUAGE DataKinds               #-}

-- For OkLC as a class
{-# LANGUAGE UndecidableInstances    #-}
{-# LANGUAGE FlexibleInstances       #-}
{-# LANGUAGE MultiParamTypeClasses   #-}

{-# OPTIONS_GHC -Wall #-}

{-# OPTIONS -Wno-type-defaults #-}

{-# OPTIONS_GHC -Wno-missing-signatures #-}
{-# OPTIONS_GHC -Wno-unused-imports #-}

{-# OPTIONS_GHC -fsimpl-tick-factor=500 #-}
{-# OPTIONS_GHC -dsuppress-idinfo #-}
{-# OPTIONS_GHC -fdicts-strict #-}
{-# OPTIONS_GHC -Wno-orphans #-}

module Main where

import Prelude hiding (unzip,zip,zipWith) -- (id,(.),curry,uncurry)
import qualified Prelude as P

import Data.Monoid (Sum(..))
import Data.Foldable (fold)
import Control.Applicative (liftA2)
import Control.Arrow (second)
import Control.Monad ((<=<))
import Data.List (unfoldr)  -- TEMP
import Data.Complex (Complex)
import GHC.Float (int2Double)

import qualified ConCat.AltCat as A
import ConCat.AltCat
 (toCcc,toCcc',unCcc,unCcc',conceal,(:**:)(..),Ok,Ok2,U2,equal)

import ConCat.Rebox

import ConCat.Circuit (GenBuses,(:>))
import ConCat.Syntactic (Syn,render)
import ConCat.RunCircuit (run)

type EC = Syn :**: (:>)

runU2 :: U2 a b -> IO ()
runU2 = print
{-# INLINE runU2 #-}

type GO a b = (GenBuses a, Ok2 (:>) a b)

runSyn :: Syn a b -> IO ()
runSyn syn = putStrLn ('\n' : render syn)
{-# INLINE runSyn #-}

runSynCirc :: GO a b => String -> EC a b -> IO ()
runSynCirc nm (syn :**: circ) = runSyn syn >> runCirc nm circ
{-# INLINE runSynCirc #-}

runCirc :: GO a b => String -> (a :> b) -> IO ()
runCirc nm circ = run nm [] circ
{-# INLINE runCirc #-}

add5 :: Double -> Double
add5 x = x + 5

decide :: Int -> Bool
decide 4 = True
decide _ = False

fix :: (a -> a) -> a
fix f = let {x = f x} in x

fac :: Int -> Int
fac 1 = 1
fac n = n*fac(n-1)

main :: IO ()
--main = print "hello world!"

main = sequence_ [
   putChar '\n' -- return ()

--  -- Circuit graphs
  , runSynCirc "add"         $ toCcc $ (+) @Double
  , runSynCirc "add5"        $ toCcc add5
  , runSynCirc "add-uncurry" $ toCcc $ uncurry ((+) @Double)
  , runSynCirc "decide"      $ toCcc $ decide
--  , runSynCirc "fac"         $ toCcc (fix (\rec n -> if n == 0 then 1 else n * rec (n-1)) :: Int -> Int)

--  , runSynCirc "fac"         $ toCcc $ fac
--  , runSynCirc "dup"         $ toCcc $ A.dup @(->) @Int
--  , runSynCirc "fst"         $ toCcc $ fst @R @R
--  , runSynCirc "twice"       $ toCcc $ twice @R
--  , runSynCirc "sqr"         $ toCcc $ sqr @R
--  , runSynCirc "complex-mul" $ toCcc $ uncurry ((*) @C)
--  , runSynCirc "magSqr"      $ toCcc $ magSqr @R
--  , runSynCirc "cosSinProd"  $ toCcc $ cosSinProd @R
--  , runSynCirc "xp3y"        $ toCcc $ \ (x,y) -> x + 3 * y :: R
--  , runSynCirc "horner"      $ toCcc $ horner @R [1,3,5]
--  , runSynCirc "cos-2xx"     $ toCcc $ \ x -> cos (2 * x * x) :: R
--
  ]

max_stars_repo_name: thma/ccc

max_issues_repo_name: thma/ccc

max_forks_repo_name: thma/ccc


text:
module Number where
import Data.Complex --Complex
import Data.Ratio --Rational

data LispNumber = Integer Integer -- Stores a Haskell Integer
                | Real Double | Rational Rational
                | Complex (Complex Double) deriving(Eq)

instance Show LispNumber where show = showNumber

instance Num LispNumber where
  (Integer a)  + (Integer b) = Integer $ a + b
  (Integer a)  + (Real b) = Real $ fromInteger a + b
  (Integer a)  + (Rational b) = Rational $ fromInteger  a + b
  (Integer a)  + (Complex b) = Complex $ fromInteger a + b

  (Real a)  + (Real b) = Real $ a + b
  (Real a)  + (Integer b) = Real $ a + fromInteger b
  (Real a)  + (Rational b) = Real $ a + fromRational b
  --(Real a)  + (Complex b) = Complex $

  (Rational a)  + (Rational b) = Rational $ a + b
  (Complex a)  + (Complex b) = Complex $ a + b

  (Integer a)  - (Integer b) = Integer $ a - b
  (Integer a)  - (Real b) = Real $ fromInteger a - b
  (Integer a)  - (Rational b) = Rational $ fromInteger  a - b
  (Integer a)  - (Complex b) = Complex $ fromInteger a - b

  (Real a)  - (Real b) = Real $ a - b
  (Rational a)  - (Rational b) = Rational $ a - b
  (Complex a)  - (Complex b) = Complex $ a - b

  (Integer a)  * (Integer b) = Integer $ a * b
  (Real a)  * (Real b) = Real $ a * b
  (Rational a)  * (Rational b) = Rational $ a * b
  (Complex a)  * (Complex b) = Complex $ a * b

  abs (Integer a)  = Integer $ abs a
  abs (Real a)     = Real $ abs a
  abs (Rational a) = Rational $ abs a
  abs (Complex a)   = Complex $ abs a

  fromInteger  = Integer

  signum (Integer a)  = Integer $ signum a
  signum (Real a)     = Real $ signum a
  signum (Rational a) = Rational $ signum a
  signum (Complex a)   = Complex $ signum a

instance Fractional LispNumber where
  (Integer a) / (Integer b)   = Rational $ a % b
  (Real a) / (Real b)         = Real $ a / b
  (Complex a) / (Complex b)   = Complex $ a / b
  (Rational a) / (Rational b) = Rational $ a / b

  fromRational = Rational

instance Ord LispNumber where
  compare (Integer a) (Integer b) = compare a b
  compare (Integer a) (Rational b) = compare (fromIntegral a) b
  compare (Integer a) (Real b) = compare (fromIntegral a) b
  compare (Integer a) (Complex b) = error "Cant order complex numbers"

  compare (Rational b) (Integer a)= compare b (fromIntegral a)
  compare (Real b) (Integer a)= compare b (fromIntegral a)
  compare (Complex b) (Integer a) = error "Cant order complex numbers"

  compare (Real a) (Real b) = compare a b
  compare (Rational a) (Rational b) = compare a b
  compare (Complex _) (Complex _) = error "Cant order complex numbers"

showNumber :: LispNumber -> String
showNumber (Integer a)  = show a
showNumber (Complex a)  = show a
showNumber (Rational a) = show (numerator a) ++ "/" ++ show (denominator a)
showNumber (Real a)     = show a

max_stars_repo_name: Unaimend/schemehs

max_issues_repo_name: Unaimend/schemehs

max_forks_repo_name: Unaimend/schemehs


text:
{-# LANGUAGE DataKinds                  #-}
{-# LANGUAGE FlexibleContexts           #-}
{-# LANGUAGE FlexibleInstances          #-}
{-# LANGUAGE GADTs                      #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE PolyKinds                  #-}
{-# LANGUAGE RankNTypes                 #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE StandaloneDeriving         #-}
{-# LANGUAGE TypeApplications           #-}
{-# LANGUAGE TypeFamilies               #-}
{-# LANGUAGE TypeOperators              #-}
{-# LANGUAGE TypeSynonymInstances       #-}
{-# LANGUAGE UndecidableInstances       #-}

{-
A lot of error handling needs to be added to the code. For example
0 < hopsize < windowsize < fftsize needs to be reinforced.
-}

module Lib.STFT where

import qualified Data.Complex                  as C
import           Data.Proxy
import           Data.Type.Equality
import qualified Data.Vector                   as V
import           GHC.TypeNats
import           Numeric.FFT.Vector.Invertible

-- A small number
epsilon :: Double
epsilon = 2.2204460492503131e-16

-- |The magnitude spectrum of a signal
newtype MagnitudeSpectrum = MkMag (V.Vector Double) deriving Show

-- |The phase spectrum of a signal
newtype PhaseSpectrum = MkPhase (V.Vector Double) deriving Show

-- |Tuple of magnitude and phase Spectrum
newtype Spectrums = MkSpect (MagnitudeSpectrum, PhaseSpectrum) deriving Show

-- |A complex signal
newtype Signal = MkSignal (V.Vector (C.Complex Double))

-- |A Window
newtype WindowFn n = Window { runWindow :: Winsz n -> V.Vector (C.Complex Double) }

type family IsPo2 (m :: Nat) :: Bool where
  IsPo2 n = 2 ^ Log2 n == n

data FFTsz n where
  MkFFTsz :: (KnownNat n, IsPo2 n ~ 'True) => Proxy n -> FFTsz n

data Winsz n where
  MkWinsz :: KnownNat n => Proxy n -> Winsz n

data Hopsz n where
  MkHopsz :: KnownNat n => Proxy n -> Hopsz n

data Config h w f where
  MkConfig :: forall h w f. (KnownNat f, KnownNat w, KnownNat h, h <= w, w <= f) => Hopsz h -> Winsz w -> FFTsz f -> Config h w f

config :: Config 256 2048 2048
config = MkConfig hop win fft where
  hop = (MkHopsz (Proxy @256))
  win = (MkWinsz (Proxy @2048))
  fft = (MkFFTsz (Proxy @2048))

hamming :: forall n. KnownNat n => WindowFn n
hamming = Window $ \(MkWinsz _) -> V.generate size goHamm
  where
    size = fromIntegral $ natVal (Proxy @n) :: Int
    goHamm :: Int -> C.Complex Double
    goHamm i = (C.:+ 0) $ 0.54 - 0.46 * cos (2 * pi * fromIntegral i / (fromIntegral size - 1))

-- |Phase unwrapping algorithm. Converting to and from list for pattern matching is a little cumbersome.
unwrapPhase :: PhaseSpectrum -> PhaseSpectrum
unwrapPhase (MkPhase xs) = MkPhase $ V.fromList $ (diff (V.toList xs) 0)
  where
    diff (x:y:ys) accm
      | (y+accm) - x > pi    = x: diff (y + accm - (2 * pi):ys) (accm-(2*pi))
      | (y+accm) - x < (-pi) = x: diff (y+accm+(2*pi):ys) (accm+(2*pi))
      | otherwise             = x: diff (y+accm:ys) accm
    diff [x] _     = [x]
    diff []  _     = []

halfInt, halfIntPlus :: Int -> Int
halfIntPlus win = (win + 1) `div` 2
halfInt win = win `div` 2

-- |Takes a complex signal and fftsz performs the zero phase windowing
-- https://ccrma.stanford.edu/~jos/sasp/Zero_Phase_Zero_Padding.html
zeroPhaseZeroPad :: forall n. KnownNat n => FFTsz n -> Signal -> Signal
zeroPhaseZeroPad _ (MkSignal sig) = MkSignal $ positive V.++ zeros V.++ negative
  where
    fftsz = fromIntegral $ natVal (Proxy @n)
    win = V.length sig
    halfSig = halfInt win
    halfSigPlus = halfIntPlus win
    positive = V.slice halfSig halfSigPlus sig
    negative = V.slice 0 halfSig sig
    zeros = V.replicate (fftsz - halfSig - halfSigPlus) (0 C.:+ 0)

nearZero :: Double -> Bool
nearZero a = abs a <= 1e-12

-- |Calculates the magnitude spectrum in Db
magSpect :: Signal -> MagnitudeSpectrum
magSpect (MkSignal sig) = MkMag $ V.map (dB. near0 . C.magnitude) sig
  where
    dB x = 20 * logBase 10 x
    near0 x = if nearZero x then epsilon else x

-- |Calculates the phase spectrum
phaseSpect :: Signal -> PhaseSpectrum
phaseSpect (MkSignal vec) = unwrapPhase $ MkPhase $ V.map (C.phase . to0) vec
  where
    to0 x = go (C.realPart x) (C.imagPart x)
    go x y | nearZero x && nearZero y = 0 C.:+ 0
           | nearZero x               = 0 C.:+ y
           | nearZero y               = x C.:+ 0
           | otherwise                = x C.:+ y

-- | Takes a windowed signal and FFT size returns a the
-- magnitude and phase spectrums. Magnitude in Db, phase unwrapped,
-- both positive halfs of the spectrums
dftAnal :: forall n. KnownNat n => FFTsz n -> Signal -> Spectrums
dftAnal fftsz sig = MkSpect (magSpect trsf, phaseSpect trsf)
  where
      (MkSignal inputVec) = zeroPhaseZeroPad fftsz sig
      -- ^Positive spectrum includes sample 0
      positiveSpect = ((fromIntegral $ natVal (Proxy @n)) `div` 2) + 1
      trsf = MkSignal $ V.slice 0 positiveSpect (run dft inputVec)

-- | Splits the vector into chunks of size n that are
-- m distance apart. Static check that chunksize is larger than hop sz.
divvy :: forall n m a. (KnownNat n, KnownNat m) => Winsz n
      -> Hopsz m -> V.Vector a -> [V.Vector a]
divvy _ _ = go n' m'
  where
    n' = fromIntegral $ natVal (Proxy :: Proxy n)
    m' = fromIntegral $ natVal (Proxy :: Proxy m)

    go :: Int -> Int -> V.Vector a -> [V.Vector a]
    go x y vec
      | V.length vec <= x = []
      | otherwise = V.take x vec : go x y (V.drop y vec)

normalize :: Fractional a => V.Vector a -> V.Vector a
normalize win = fmap (/ V.foldl1' (+) win) win

-- | Takes an input signal, a window, an fftsize and a hop size and
-- computes the short time Fourier transform, returning a list of
-- magnitude and phase spectrums.
stft :: forall h w f. (KnownNat h, KnownNat w, KnownNat f) => Config h w f -> Signal
     -> [Spectrums]
stft (MkConfig hopsz winsz fftsz) (MkSignal sig) = fmap winDFT splitSig
  where
    win = runWindow hamming winsz
    halfWindow = floor (fromIntegral (natVal (Proxy @w)) / 2 :: Double)
    zs = V.replicate halfWindow (0 C.:+ 0)
    sigzs = zs V.++ sig V.++ zs
    splitSig = divvy winsz hopsz sigzs

    winDFT :: V.Vector (C.Complex Double) -> Spectrums
    winDFT vec = dftAnal fftsz (MkSignal $ V.zipWith (*) vec (normalize win))

-- Takes a tuple of vectors (Magnitude, Phase) and a window size
-- and returns the original signal
dftSynth :: forall n. KnownNat n => Winsz n
         -> Spectrums -> V.Vector (C.Complex Double)
dftSynth _ (MkSpect (MkMag magVec, MkPhase phaseVec)) =
  let win = fromIntegral $ natVal (Proxy @n)
      hwin = halfInt win
      hwinp = halfIntPlus win
      vec = V.zipWith (\x y -> (x, y)) magVec phaseVec
      posFreqs = fmap f vec where
        f (mag, phase) = (10**(mag/20) C.:+ 0) *
                         exp ((phase C.:+ 0) *
                         (0 C.:+ 1))
      negFreqs = (fmap f . V.reverse) (g vec) where
        g = (V.init . V.tail)
        f (mag, phase) = (10**(mag/20) C.:+ 0) *
                         exp ((phase C.:+ 0) *
                         (0 C.:+ (-1)))
      resultidft = run idft (posFreqs V.++ negFreqs)
  in V.slice (V.length resultidft - hwin)
                    hwin resultidft V.++ V.take hwinp resultidft

stftSynth :: forall h w. (KnownNat h, KnownNat w) => [Spectrums]
          -> Hopsz h -> Winsz w -> Signal
stftSynth magphase _ winsz =
  let hwinp = halfIntPlus $ fromIntegral $ natVal (Proxy @w)
      signalFrames = fmap (fmap (fromIntegral (natVal (Proxy @h)) * ) . dftSynth winsz ) magphase
      signalTuples = fmap (V.splitAt hwinp) signalFrames
      overlapAdd (x1, x2) (y1, y2) = (x1 V.++ V.zipWith (+) x2 y1, y2)
   in MkSignal $ V.drop hwinp $ fst (foldl overlapAdd (Prelude.head signalTuples) (Prelude.tail signalTuples))

max_stars_repo_name: davlum/haskell-twm

max_issues_repo_name: davlum/haskell-twm

max_forks_repo_name: davlum/haskell-twm