text:
--
-- This file was automatically generated by the husk scheme compiler (huskc)
--
--  http://justinethier.github.io/husk-scheme
--  (c) 2010 Justin Ethier
--  Version 3.16
--
module Main where
import Language.Scheme.Core
import Language.Scheme.Numerical
import Language.Scheme.Primitives
import Language.Scheme.Types     -- Scheme data types
import Language.Scheme.Variables -- Scheme variable operations
import Control.Monad.Error
import Data.Array
import  qualified Data.ByteString as BS
import Data.Complex
import  qualified Data.Map
import Data.Ratio
import Data.Word
import System.IO
import Debug.Trace

-- |Get variable at runtime
getRTVar env var = do
  v <- getVar env var
  return $ case v of
    List _ -> Pointer var env
    DottedList _ _ -> Pointer var env
    String _ -> Pointer var env
    Vector _ -> Pointer var env
    ByteVector _ -> Pointer var env
    HashTable _ -> Pointer var env
    _ -> v

applyWrapper env cont (Nil _) (Just (a:as))  = do
  apply cont a as

applyWrapper env cont value (Just (a:as))  = do
  apply cont a $ as ++ [value]

getDataFileName' :: FilePath -> IO FilePath
getDataFileName' name = return $ "/home/justin/.cabal/share/husk-scheme-3.16/" ++ name

exec55_3 env cont _ _ = do
  liftIO $ registerExtensions env getDataFileName'
  continueEval env (makeCPSWArgs env cont exec []) (Nil "")

main :: IO ()
main = do
  env <- r5rsEnv
  result <- (runIOThrows $ liftM show $ hsInit env (makeNullContinuation env) (Nil "") Nothing)
  case result of
    Just errMsg -> putStrLn errMsg
    _ -> return ()

hsInit env cont _ _ = do
  _ <- defineVar env " modules " $ HashTable $ Data.Map.fromList []
  run env cont (Nil "") (Just [])

exec _ _ _ _ = return $ Nil ""

run :: Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal
run env cont _ _  = do
  result <- makeNormalHFunc env (["obj","alist"]) defineFuncEntryPt480
  _ <- defineVar env "assv" result
  (trace ("run") continueEval) env (makeCPSWArgs env cont f479 []) result
defineFuncEntryPt480 env cont value (Just args) = do
  value <- getRTVar env "foldl"
  (trace ("defineFuncEntryPt480 " ++ (show value) ++ " " ++ (show args)) continueEval) env (makeCPSWArgs env (makeCPSWArgs env cont applyNextArg482 $ args ++ [value] ++ []) applyFirstArg481 []) (Nil "")
applyFirstArg481 env cont value (Just args) = do
  value <- getRTVar env "mem-helper"
  (trace ("applyFirstArg481" ++ (show value) ++ " " ++ (show args)) continueEval) env (makeCPSWArgs env (makeCPSWArgs env cont applyNextArg484 $ args ++ [value] ++ []) applyFirstArg483 []) (Nil "")

applyFirstArg483 :: Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal
applyFirstArg483 env cont _ _  = do
  (trace ("applyFirstArg483") continueEval) env (makeCPSWArgs env (makeCPSWArgs env cont applyNextF486 []) applyStubF485 []) $ Nil""

applyStubF485 :: Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal
applyStubF485 env cont _ _  = do
  val <- getRTVar env "curry"
  (trace ("applyStubF485") continueEval) env cont val

applyNextF486 :: Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal
applyNextF486 env cont value _  = do
  v0 <- getRTVar env "eqv?"
  v1 <- getRTVar env "obj"
  (trace ("applyNextF486 " ++ (show value)) apply) cont value [v0,v1]
applyNextArg484 env cont value (Just args) = do

  (trace ("applyNextArg484 " ++ (show value) ++ " " ++ (show args)) continueEval) env (makeCPSWArgs env (makeCPSWArgs env cont applyWrapper $ args ++ [value] ++ []) applyFirstArg487 []) (Nil "")

applyFirstArg487 :: Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal
applyFirstArg487 env cont value _  = do
  val <- getRTVar env "car"
  (trace ("applyFirstArg487 " ++ (show value) ++ " []") continueEval) env cont val
applyNextArg482 env cont value (Just args) = do

  (trace ("applyNextArg482 " ++ (show value) ++ " " ++ (show args)) continueEval) env (makeCPSWArgs env (makeCPSWArgs env cont applyNextArg489 $ args ++ [value] ++ []) applyFirstArg488 []) (Nil "")

applyFirstArg488 :: Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal
applyFirstArg488 env cont value _  = do
  x1 <-   return $ Bool False
  (trace ("applyFirstArg488 " ++ (show value) ++ " " ++ (show "[]")) continueEval) env cont x1
applyNextArg489 env cont value (Just args) = do

  (trace ("applyNextArg489 " ++ (show value) ++ " " ++ (show args)) continueEval) env (makeCPSWArgs env (makeCPSWArgs env cont applyWrapper $ args ++ [value] ++ []) applyFirstArg490 []) (Nil "")

applyFirstArg490 :: Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal
applyFirstArg490 env cont value _  = do
  val <- getRTVar env "alist"
  (trace ("" ++ (show value) ++ " " ++ (show "[]")) continueEval) env cont val
f479 env cont value (Just args) = do
  value <- getRTVar env "assv"
  (trace ("f479 " ++ (show value) ++ " " ++ (show args)) continueEval) env (makeCPSWArgs env (makeCPSWArgs env cont applyNextArg492 $ args ++ [value] ++ []) applyFirstArg491 []) (Nil "")

applyFirstArg491 :: Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal
applyFirstArg491 env cont value _  = do
  x1 <-  return $ Atom "b"
  (trace ("applyFirstArg491 " ++ (show value) ++ " " ++ (show "[]")) continueEval) env cont x1
applyNextArg492 env cont value (Just args) = do

  (trace ("applyNextArg492 " ++ (show value) ++ " " ++ (show args)) continueEval) env (makeCPSWArgs env (makeCPSWArgs env cont applyWrapper $ args ++ [value] ++ []) applyFirstArg493 []) (Nil "")

applyFirstArg493 :: Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal
applyFirstArg493 env cont value _  = do
  x1 <-  return $ List [List [Atom "a",Number (1)]]
  (trace ("applyFirstArg493 " ++ (show value) ++ " " ++ (show "[]")) continueEval) env cont x1

max_stars_repo_name: svenkeidel/husk-scheme

max_issues_repo_name: egisatoshi/husk-scheme

max_forks_repo_name: egisatoshi/husk-scheme


text:
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeOperators #-}
module Filter(
      FIR(..)
    , bode
    , phasePlot
    ) where

import Signal
import Math.Polynomial
import Common
import Noise
import Fixed
import GHC.TypeLits
import SpecialInt
import Data.Bits
import Plot
import Data.Complex
import Text.Printf

import Debug.Trace

debug a = trace (show a) a

data FIR p i o where
 FIR :: (SingI (nc + ni), SingI nc, SingI ni, SingI no, SingI s, SingI r, HasDoubleRepresentation (Fixed a nc s r), Num (Fixed ag (nc + ni) s r), Num (Fixed a nc s r), Eq a, Conversion ag ao, ConvertConstraint ag ao, AMulConstraint a ag)
     => Poly (Fixed a nc s r)
     -> FIR ((Fixed ag (nc + ni) s r),(Fixed a nc s r)) (Fixed a ni s r) (Fixed ao no s r)
 FIRD :: Poly Double
      -> FIR Double Double Double

polyRepresentation :: FIR p i o -> Poly Double
polyRepresentation (FIRD p) = p
polyRepresentation (FIR f) = fmap toDouble f

instance Show p => Show (FIR p i o) where
  show (FIR c) = show c
  show (FIRD d) = show d

delay :: s -> Signal s -> Signal s
delay a s = consS a s

guardType :: FIR ((Fixed ag ng s r),(Fixed a nc s r)) i (Fixed ao no s r)
          -> Fixed ag n s r
          -> Fixed ag n s r
guardType _ = id

filterDWith :: [Double]
            -> Signal Double
            -> Signal Double
filterDWith (a:l) s | null l = mapS (a*) s
                    | otherwise = zipWithS (+) (mapS (a*) s) (filterDWith l (delay 0 s))
filterDWith [] s = s

_filterWith :: (AMulConstraint a ag, SingI (na + ni), SingI s, SingI r, SingI na, SingI ni, Num (Fixed a ni s r), Num (Fixed ag (na + ni) s r))
            => FIR ((Fixed ag ng s r),(Fixed a na s r)) (Fixed a ni s r) (Fixed ao no s r)
            -> [Fixed a na s r]
            -> Signal (Fixed a ni s r)
            -> Signal (Fixed ag (na + ni) s r)
_filterWith f (a:l) s | null l = mapS (amul a) s
                      | otherwise = zipWithS (+) (mapS (amul a) s) (_filterWith f l (delay 0 s))
_filterWith _ [] s = error "Empty polynomial can't be used for filtering"

filterWith :: (Conversion ag ao, ConvertConstraint ag ao, Num (Fixed a ni s r), Num (Fixed ag (na + ni) s r), SingI r, SingI s, SingI (na + ni), SingI ni, SingI na, SingI no, AMulConstraint a ag)
           => FIR ((Fixed ag ng s r),(Fixed a na s r)) (Fixed a ni s r) (Fixed ao no s r)
           -> [Fixed a na s r]
           -> Signal (Fixed a ni s r)
           -> Signal (Fixed ao no s r)
filterWith f l = mapS convert . _filterWith f l

instance Structure FIR where
    doubleVersion (FIR p) = FIRD $ fmap toDouble p
    transferFunction f@(FIR p) s =
        let c = polyCoeffs LE $ p
        in
        filterWith f c s
    transferFunction (FIRD p) s =
        let c = polyCoeffs LE p
        in
        filterDWith c s

bode :: (Sample o,Sample i)
     => FIR p i o -> StyledSignal Double Double
bode f =
    let fd = polyRepresentation f
        log10 x = log x / log 10
        polyF = evalPoly (fmap (:+ 0) fd)
        transferFunction = \t -> 20 * log10 (magnitude . polyF . cis $ -t)
        dw = 0.01
        values = map transferFunction [0,dw..pi]
        nb = length values
        ticksWithPhase n = printf "%.2f pi" (n / fromIntegral nb :: Double)
    in
    (plotSpectrum nb [ AS $ fromListS 0 values])
                 `withNewStyle` (\s -> s { verticalLabel = Just "Filter Gain (dB)"
                                         , title = Just "Bode Plot"
                                         , horizontalTickRepresentation = ticksWithPhase
                                         })

phasePlot :: (Sample o,Sample i)
     => FIR p i o -> StyledSignal Double Double
phasePlot f =
    let fd = polyRepresentation f
        log10 x = log x / log 10
        polyF = evalPoly (fmap (:+ 0) fd)
        transferFunction = \t -> phase . polyF . cis $ -t
        dw = 0.01
        values = map transferFunction [0,dw..pi]
        nb = length values
        ticksWithPhase n = printf "%.2f pi" (n / fromIntegral nb :: Double)
    in
    (plotSpectrum nb [ AS $ fromListS 0 values])
                 `withNewStyle` (\s -> s { verticalLabel = Just "Filter Phase"
                                         , title = Just "Phase Plot"
                                         , horizontalTickRepresentation = ticksWithPhase
                                         })

ext: hs

max_stars_repo_name: alpheccar/Signal

max_issues_repo_name: cpehle/Signal

max_forks_repo_name: cpehle/Signal


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-|
Module      : Grenade.Layers.LeakyRelu
Description : Rectifying linear unit layer
Copyright   : (c) Manuel Schneckenreither, 2020
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.LeakyRelu (
    LeakyRelu (..)
  , SpecLeakyRelu (..)
  , specLeakyRelu1D
  , specLeakyRelu2D
  , specLeakyRelu3D
  , leakyRelu
  ) where

import           Control.DeepSeq                (NFData (..), force)
import           Data.Constraint                (Dict (..))
import           Data.Reflection                (reifyNat)
import           Data.Serialize
import           Data.Singletons
import           GHC.Generics                   (Generic)
import           GHC.TypeLits
import qualified Numeric.LinearAlgebra.Static   as LAS
import           Unsafe.Coerce                  (unsafeCoerce)

import           Grenade.Core
import           Grenade.Dynamic
import           Grenade.Dynamic.Internal.Build
import           Grenade.Types
import           Grenade.Utils.Conversion       (toLayerShape)
import           Grenade.Utils.Vector

-- | A rectifying linear unit.
--   A layer which can act between any shape of the same dimension, acting as a
--   diode on every neuron individually.
data LeakyRelu = LeakyRelu
  deriving (Generic, NFData, Show)

instance UpdateLayer LeakyRelu where
  type Gradient LeakyRelu = ()
  runUpdate _ _ _ = LeakyRelu

instance RandomLayer LeakyRelu where
  createRandomWith _ _ = return LeakyRelu

instance Serialize LeakyRelu where
  put _ = return ()
  get = return LeakyRelu

instance (KnownNat i) => Layer LeakyRelu ('D1 i) ('D1 i) where
  type Tape LeakyRelu ('D1 i) ('D1 i) = S ('D1 i)

  runForwards _ (S1D y) = (S1D y, S1D (relu y))
    where
      relu = LAS.dvmap (\a -> if a < 0 then alpha * a else a)
  runForwards _ (S1DV y) = (S1DV y, S1DV (mapVectorInPlace relu y))
    where
      relu = (\a -> if a < 0 then alpha * a else a) -- we do not change the sign, so we can use inPlace here
  runBackwards _ (S1D y) (S1D dEdy) = ((), S1D (relu' y * dEdy))
    where
      relu' = LAS.dvmap (\a -> if a < 0 then alpha else 1)
  runBackwards _ (S1DV y) (S1DV dEdy) = ((), zipWithVectorInPlaceSnd reluDif y dEdy `seq` S1DV dEdy)
  runBackwards x y dEdy = runBackwards x y (toLayerShape y dEdy)

instance (KnownNat i, KnownNat j) => Layer LeakyRelu ('D2 i j) ('D2 i j) where
  type Tape LeakyRelu ('D2 i j) ('D2 i j) = S ('D2 i j)

  runForwards _ (S2D y) = (S2D y, S2D (relu y))
    where
      relu = LAS.dmmap (\a -> if a < 0 then alpha * a else a)
  runForwards _ (S2DV y) = (S2DV y, S2DV (mapVectorInPlace relu y))
    where
      relu = (\a -> if a < 0 then alpha * a else a)
  runBackwards _ (S2D y) (S2D dEdy) = ((), S2D (relu' y * dEdy))
    where
      relu' = LAS.dmmap (\a -> if a < 0 then alpha else 1)
  runBackwards _ (S2DV y) (S2DV dEdy) = ((), zipWithVectorInPlaceSnd reluDif y dEdy `seq` S2DV dEdy)
  runBackwards x y dEdy = runBackwards x y (toLayerShape y dEdy)

instance (KnownNat i, KnownNat j, KnownNat k) => Layer LeakyRelu ('D3 i j k) ('D3 i j k) where

  type Tape LeakyRelu ('D3 i j k) ('D3 i j k) = S ('D3 i j k)

  runForwards _ (S3D y) = (S3D y, S3D (relu y))
    where
      relu = LAS.dmmap (\a -> if a < 0 then alpha * a else a)
  runBackwards _ (S3D y) (S3D dEdy) = ((), S3D (relu' y * dEdy))
    where
      relu' = LAS.dmmap (\a -> if a < 0 then alpha else 1)

alpha :: RealNum
alpha = 0.02

reluDif :: RealNum -> RealNum -> RealNum
reluDif a g
  | a < 0 = alpha * g
  | otherwise = g

-------------------- DynamicNetwork instance --------------------

instance FromDynamicLayer LeakyRelu where
  fromDynamicLayer inp _ _ = SpecNetLayer $ SpecLeakyRelu (tripleFromSomeShape inp)

instance ToDynamicLayer SpecLeakyRelu where
  toDynamicLayer _ _ (SpecLeakyRelu (rows, cols, depth)) =
     reifyNat rows $ \(_ :: (KnownNat rows) => Proxy rows) ->
     reifyNat cols $ \(_ :: (KnownNat cols) => Proxy cols) ->
     reifyNat depth $ \(_ :: (KnownNat depth) => Proxy depth) ->
     case (rows, cols, depth) of
         (_, 1, 1)    -> return $ SpecLayer LeakyRelu (sing :: Sing ('D1 rows)) (sing :: Sing ('D1 rows))
         (_, _, 1) -> return $ SpecLayer LeakyRelu (sing :: Sing ('D2 rows cols)) (sing :: Sing ('D2 rows cols))
         _    -> case (unsafeCoerce (Dict :: Dict()) :: Dict (KnownNat (rows GHC.TypeLits.* depth))) of
           Dict -> return $ SpecLayer LeakyRelu (sing :: Sing ('D3 rows cols depth)) (sing :: Sing ('D3 rows cols depth))

-- | Create a specification for a elu layer.
specLeakyRelu1D :: Integer -> SpecNet
specLeakyRelu1D i = specLeakyRelu3D (i, 1, 1)

-- | Create a specification for a elu layer.
specLeakyRelu2D :: (Integer, Integer) -> SpecNet
specLeakyRelu2D (i, j) = specLeakyRelu3D (i, j, 1)

-- | Create a specification for a elu layer.
specLeakyRelu3D :: (Integer, Integer, Integer) -> SpecNet
specLeakyRelu3D = SpecNetLayer . SpecLeakyRelu

-- | Add a LeakyRelu layer to your build.
leakyRelu :: BuildM ()
leakyRelu = buildGetLastLayerOut >>= buildAddSpec . SpecNetLayer . SpecLeakyRelu

-------------------- GNum instances --------------------

instance GNum LeakyRelu where
  _ |* LeakyRelu = LeakyRelu
  _ |+ LeakyRelu = LeakyRelu

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


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

-- | Functions to represent a 'Vector' on disk in efficient, if
-- unportable, ways.
--
-- This module uses memory-mapping, a feature of all modern
-- operating-systems, to mirror the disk contents in memory. There are
-- quite a few advantages to memory-mapping files instead of reading
-- the files traditionally:
--
--  * Speed: memory-mapping is often much faster than traditional
--    reading.
--
--  * Memory efficiency: Memory-mapped files are loaded into RAM
--    on-demand, and easily swapped out. The upside is that the
--    program can work with data-sets larger than the available RAM,
--    as long as they are accessed carefully.
--
-- The caveat to using memory-mapping is that it makes the files
-- specific to the current architecture because of the endianness of
-- the data. For more information, see the description in
-- "System.IO.MMap"
--
-- If you wish to write the contents in a portable fashion, either use
-- the ASCII load and save functions in "Numeric.Container", or use
-- the binary serialization in "Data.Binary".

module Data.Packed.Vector.MMap (
  -- * Memory-mapping 'Vector' from disk
  unsafeMMapVector,
  unsafeLazyMMapVectors,

  -- * Writing 'Vector' to disk

  -- | These functions write the 'Vector' in a way suitable for
  -- reading back with 'unsafeMMapVector'.
  hPutVector,
  writeVector
) where

import Control.Monad (when)

import System.IO
import System.IO.MMap
import System.IO.Unsafe

import Foreign.ForeignPtr
import Foreign.Ptr
import Foreign.Storable

import qualified Numeric.LinearAlgebra.Devel as I
import qualified Numeric.LinearAlgebra.HMatrix as I
-- import qualified Data.Packed.Vector as I
import Data.Int

---------------------------
-- Memory-Mapping 'Vector' from disk

-- | Map a file into memory (read-only) as a 'Vector'.
--
-- It is considered unsafe because changes to the underlying file may
-- (or may not) be reflected in the 'Vector', which breaks referential
-- transparency.
unsafeMMapVector :: forall a. Storable a => FilePath -- ^ Path of the file to map
                                         -> Maybe (Int64, Int) -- ^ 'Nothing' to map entire file into memory, otherwise 'Just (fileOffset, elementCount)'
                                         -> IO (I.Vector a)
unsafeMMapVector path range =
  do (foreignPtr, offset, size) <- mmapFileForeignPtr path ReadOnly $
        case range of
          Nothing -> Nothing
          Just (start, length) -> Just (start, length * sizeOf (undefined :: a))
     return $ I.unsafeFromForeignPtr foreignPtr offset (size `div` sizeOf (undefined :: a))

-- | Map a file into memory as a lazy-list of equal-sized 'Vector',
-- even if they can't all fit in the address space at the same time.
--
-- > (numVectors,vectors) <- unsafeLazyMMapVectors filename Nothing vectorSize
--
-- Commonly, a data file will contain multiple vectors of equal length
-- (matrix). This function is convenient for those uses, but it plays
-- a more important role: supporting data-sets that cannot fit in the
-- address space of the current machine.
--
-- On 32-bit machines the address space is only 4GB, and it is
-- actually pretty easy to find data-sets that are too large to be
-- represented, even in virtual memory.
--
-- This function loads the data in chunks, and as long as you drop
-- your reference to the vectors as you consume the data, the old
-- chunks will be unmapped before mapping the next chunk.
--
-- The number of vectors in the list is returned because it's often
-- needed, yet calculating it using 'length' would demand the whole
-- list.
unsafeLazyMMapVectors :: forall a. (Storable a, I.Container I.Vector a) => FilePath -- ^ Path of the file to map
                      -> Maybe (Int64, Int64)
                      -- ^ 'Nothing' to map entire file into memory,
                      -- otherwise @'Just' (fileOffset, totalElementCount)@
                      -> Int -- ^ The number of elements in each 'Vector'
                      -> IO (Int64,[I.Vector a]) -- ^ Return @(numberOfVectors,vectors)@
unsafeLazyMMapVectors path range vsize = do
  when (vecSize > maxChunkSize) vecTooBigError
  filesize <- withFile path ReadMode hFileSize
  let filesize' :: Int64
      filesize' = fI filesize
  imgs <- unsafeInterleaveIO $ unsafeLazyMMapVectors' filesize' path range vsize
  return (nimages range filesize', imgs)
      where
        nimages :: Maybe (Int64, Int64) -> Int64 -> Int64
        nimages Nothing fsz = fsz `div` imageSize
        nimages (Just (_,sz)) _ = sz `div` imageSize
        imageSize = fI vsize * eltSize
        eltSize = fI (sizeOf (undefined :: a))
        vecSize = fI vsize * eltSize
        vecTooBigError = fail "The requested vector size can't be mapped into memory"

unsafeLazyMMapVectors' :: forall a. (Storable a, I.Container I.Vector a) => Int64
                       -> FilePath
                       -> Maybe (Int64, Int64)
                       -> Int
                       -> IO [I.Vector a]
unsafeLazyMMapVectors' fileSize
                       fileName
                       fileRange
                       numEltsPerVec
                           | mapSize < maxChunkSize = mmapAll
                           | otherwise = mmapChunks 0
    where
      mapSize, eltSize, vecSize, chunkSize, baseOffset :: Int64
      eltSize = fI $ sizeOf (undefined :: a)
      (baseOffset,mapSize) = case fileRange of
                               Just (off,nelts) -> (off,nelts*eltSize)
                               _ -> (0,fileSize)
      vecSize = fI numEltsPerVec * eltSize
      chunkSize = (maxChunkSize `div` vecSize) * vecSize

      fileRange' = do
        (offset, nelts) <- fileRange
        return (offset, fI nelts)

      splitVecs :: I.Vector a -> [I.Vector a]
      splitVecs bigVec = let nvecs = I.size bigVec `div` numEltsPerVec
                         in I.takesV (replicate nvecs numEltsPerVec) bigVec

      mmapAll :: IO [I.Vector a]
      mmapAll = do
        allVecs <- unsafeMMapVector fileName fileRange'
        return $ splitVecs allVecs

      mmapChunks :: Int64 -> IO [I.Vector a]
      mmapChunks offs | remaining <= 0 = return []
                      | otherwise = do
        chunk <- unsafeMMapVector fileName mmapRange
        rest <- unsafeInterleaveIO $ mmapChunks (offs+chunkSize')
        return $ splitVecs chunk ++ rest
          where
            mmapRange = Just (baseOffset+offs,fI (chunkSize' `div` eltSize))
            remaining = mapSize-offs
            chunkSize' = min chunkSize remaining

-- Maximum size for chunks
maxChunkSize :: Int64
maxChunkSize = fI (maxBound `div` 256 :: Int)

-- Handy alias for 'fromIntegral'
fI :: (Integral a, Num b) => a -> b
fI = fromIntegral

---------------------------
-- Writing 'Vector' to disk

-- | Write out a vector verbatim into an open file handle.
hPutVector :: forall a. Storable a => Handle -> I.Vector a -> IO ()
hPutVector h v = withForeignPtr fp $ \p -> hPutBuf h (p `plusPtr` offset) sz
      where
        (fp, offset, n) = I.unsafeToForeignPtr v
        eltsize = sizeOf (undefined :: a)
        sz = n * eltsize

-- | Write the vector verbatim to a file.
writeVector :: forall a. Storable a => FilePath -> I.Vector a -> IO ()
writeVector fp v = withFile fp WriteMode $ \h -> hPutVector h v

max_stars_repo_name: haedosa/hmatrix-mmap

max_issues_repo_name: haedosa/hmatrix-mmap

max_forks_repo_name: haedosa/hmatrix-mmap


text:
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE TupleSections #-}

module HackageBenchmark (
    hackageBenchmarkMain

-- Exposed for testing:
  , CabalResult(..)
  , isSignificantTimeDifference
  , combineTrialResults
  , isSignificantResult
  , shouldContinueAfterFirstTrial
  ) where

import Control.Monad (forM_, replicateM, unless, when)
import qualified Data.ByteString as B
import Data.List (nub, unzip4)
import Data.Maybe (isJust)
import Data.Monoid ((<>))
import Data.String (fromString)
import Data.Time (NominalDiffTime, diffUTCTime, getCurrentTime)
import qualified Data.Vector.Unboxed as V
import Options.Applicative
import Statistics.Sample (mean, stdDev)
import Statistics.Test.MannWhitneyU ( PositionTest(..), TestResult(..)
                                    , mannWhitneyUCriticalValue
                                    , mannWhitneyUtest)
import Statistics.Types (PValue, mkPValue)
import System.Directory (getTemporaryDirectory)
import System.Exit (ExitCode(..), exitFailure)
import System.FilePath ((</>))
import System.IO ( BufferMode(LineBuffering), hPutStrLn, hSetBuffering, stderr
                 , stdout)
import System.Process ( StdStream(CreatePipe), CreateProcess(..), callProcess
                      , createProcess, readProcess, shell, waitForProcess )
import Text.Printf (printf)

import Distribution.Package (PackageName, mkPackageName, unPackageName)

data Args = Args {
    argCabal1                      :: FilePath
  , argCabal2                      :: FilePath
  , argCabal1Flags                 :: [String]
  , argCabal2Flags                 :: [String]
  , argPackages                    :: [PackageName]
  , argMinRunTimeDifferenceToRerun :: Double
  , argPValue                      :: PValue Double
  , argTrials                      :: Int
  , argPrintTrials                 :: Bool
  , argPrintSkippedPackages        :: Bool
  , argTimeoutSeconds              :: Int
  }

data CabalTrial = CabalTrial NominalDiffTime CabalResult

data CabalResult
  = Solution
  | NoInstallPlan
  | BackjumpLimit
  | Unbuildable
  | UnbuildableDep
  | ComponentCycle
  | ModReexpIssue
  | PkgNotFound
  | Timeout
  | Unknown
  deriving (Eq, Show)

hackageBenchmarkMain :: IO ()
hackageBenchmarkMain = do
  hSetBuffering stdout LineBuffering
  args@Args {..} <- execParser parserInfo
  checkArgs args
  printConfig args
  pkgs <- getPackages args
  putStrLn ""

  let -- The maximum length of the heading and package names.
      nameColumnWidth :: Int
      nameColumnWidth =
          maximum $ map length $ "package" : map unPackageName pkgs
      runCabal1 = runCabal argTimeoutSeconds argCabal1 argCabal1Flags
      runCabal2 = runCabal argTimeoutSeconds argCabal2 argCabal2Flags

  -- When the output contains both trails and summaries, label each row as
  -- "trial" or "summary".
  when argPrintTrials $ putStr $ printf "%-16s " "trial/summary"
  putStrLn $
      printf "%-*s %-14s %-14s %11s %11s %11s %11s %11s"
             nameColumnWidth "package" "result1" "result2"
             "mean1" "mean2" "stddev1" "stddev2" "speedup"

  forM_ pkgs $ \pkg -> do
    let printTrial msgType result1 result2 time1 time2 =
            putStrLn $
            printf "%-16s %-*s %-14s %-14s %10.3fs %10.3fs"
                   msgType nameColumnWidth (unPackageName pkg)
                   (show result1) (show result2)
                   (diffTimeToDouble time1) (diffTimeToDouble time2)

    CabalTrial t1 r1 <- runCabal1 pkg
    CabalTrial t2 r2 <- runCabal2 pkg
    if not $
       shouldContinueAfterFirstTrial argMinRunTimeDifferenceToRerun t1 t2 r1 r2
    then when argPrintSkippedPackages $
         if argPrintTrials
         then printTrial "trial (skipping)" r1 r2 t1 t2
         else putStrLn $ printf "%-*s (first run times were too similar)"
                                nameColumnWidth (unPackageName pkg)
    else do
      when argPrintTrials $ printTrial "trial" r1 r2 t1 t2
      (ts1, ts2, rs1, rs2) <- (unzip4 . ((t1, t2, r1, r2) :) <$>)
                            . replicateM (argTrials - 1) $ do
        CabalTrial t1' r1' <- runCabal1 pkg
        CabalTrial t2' r2' <- runCabal2 pkg
        when argPrintTrials $ printTrial "trial" r1' r2' t1' t2'
        return (t1', t2', r1', r2')

      let result1 = combineTrialResults rs1
          result2 = combineTrialResults rs2
          times1 = V.fromList (map diffTimeToDouble ts1)
          times2 = V.fromList (map diffTimeToDouble ts2)
          mean1 = mean times1
          mean2 = mean times2
          stddev1 = stdDev times1
          stddev2 = stdDev times2
          speedup = mean1 / mean2

      when argPrintTrials $ putStr $ printf "%-16s " "summary"
      if isSignificantResult result1 result2
          || isSignificantTimeDifference argPValue ts1 ts2
      then putStrLn $
           printf "%-*s %-14s %-14s %10.3fs %10.3fs %10.3fs %10.3fs %10.3f"
                  nameColumnWidth (unPackageName pkg)
                  (show result1) (show result2) mean1 mean2 stddev1 stddev2 speedup
      else when (argPrintTrials || argPrintSkippedPackages) $
           putStrLn $
           printf "%-*s (not significant)" nameColumnWidth (unPackageName pkg)
  where
    checkArgs :: Args -> IO ()
    checkArgs Args {..} = do
      let die msg = hPutStrLn stderr msg >> exitFailure
      unless (argTrials > 0) $ die "--trials must be greater than 0."
      unless (argMinRunTimeDifferenceToRerun >= 0) $
          die "--min-run-time-percentage-difference-to-rerun must be non-negative."
      unless (isSampleLargeEnough argPValue argTrials) $
          die "p-value is too small for the number of trials."

    printConfig :: Args -> IO ()
    printConfig Args {..} = do
      putStrLn "Comparing:"
      putStrLn $ "1: " ++ argCabal1 ++ " " ++ unwords argCabal1Flags
      callProcess argCabal1 ["--version"]
      putStrLn $ "2: " ++ argCabal2 ++ " " ++ unwords argCabal2Flags
      callProcess argCabal2 ["--version"]
      -- TODO: Print index state.
      putStrLn "Base package database:"
      callProcess "ghc-pkg" ["list"]

    getPackages :: Args -> IO [PackageName]
    getPackages Args {..} = do
      pkgs <-
          if null argPackages
          then do
            putStrLn $ "Obtaining the package list (using " ++ argCabal1 ++ ") ..."
            list <- readProcess argCabal1 ["list", "--simple-output"] ""
            return $ nub [mkPackageName $ head (words line) | line <- lines list]
          else do
            putStrLn "Using given package list ..."
            return argPackages
      putStrLn $ "Done, got " ++ show (length pkgs) ++ " packages."
      return pkgs

runCabal :: Int -> FilePath -> [String] -> PackageName -> IO CabalTrial
runCabal timeoutSeconds cabal flags pkg = do
  ((exitCode, err), time) <- timeEvent $ do
    tmpDir <- getTemporaryDirectory

    let timeout = "timeout --foreground -sINT " ++ show timeoutSeconds
        cabalCmd = unwords $
            [ cabal

              -- A non-existent store directory prevents cabal from reading the
              -- store, which would cause the size of the store to affect run
              -- time.
            , "--store-dir=" ++ (tmpDir </> "non-existent-store-dir")

            , "v2-install"

              -- These flags prevent a Cabal project or package environment from
              -- affecting the install plan.
            , "--ignore-project"
            , "--package-env=non-existent-package-env"

              -- --lib allows solving for packages with libraries or
              -- executables.
            , "--lib"

            , unPackageName pkg

            , "--dry-run"

              -- The test doesn't currently handle stdout, so we suppress it
              -- with silent. nowrap simplifies parsing the errors messages.
            , "-vsilent+nowrap"]

             ++ flags
        cmd = (shell (timeout ++ " " ++ cabalCmd)) { std_err = CreatePipe }

    -- TODO: Read stdout and compare the install plans.
    (_, _, Just errh, ph) <- createProcess cmd
    err <- B.hGetContents errh
    (, err) <$> waitForProcess ph
  let exhaustiveMsg =
          "After searching the rest of the dependency tree exhaustively"
      result
        | exitCode == ExitSuccess                                                          = Solution
        | exitCode == ExitFailure 124                                                      = Timeout
        | fromString exhaustiveMsg `B.isInfixOf` err                                       = NoInstallPlan
        | fromString "Backjump limit reached" `B.isInfixOf` err                            = BackjumpLimit
        | fromString "none of the components are available to build" `B.isInfixOf` err     = Unbuildable
        | fromString "Dependency on unbuildable" `B.isInfixOf` err                         = UnbuildableDep
        | fromString "Dependency cycle between the following components" `B.isInfixOf` err = ComponentCycle
        | fromString "Problem with module re-exports" `B.isInfixOf` err                    = ModReexpIssue
        | fromString "There is no package named" `B.isInfixOf` err                         = PkgNotFound
        | otherwise                                                                        = Unknown
  return (CabalTrial time result)

isSampleLargeEnough :: PValue Double -> Int -> Bool
isSampleLargeEnough pvalue trials =
    -- mannWhitneyUCriticalValue, which can fail with too few samples, is only
    -- used when both sample sizes are less than or equal to 20.
    trials > 20 || isJust (mannWhitneyUCriticalValue (trials, trials) pvalue)

isSignificantTimeDifference :: PValue Double -> [NominalDiffTime] -> [NominalDiffTime] -> Bool
isSignificantTimeDifference pvalue xs ys =
  let toVector = V.fromList . map diffTimeToDouble
  in case mannWhitneyUtest SamplesDiffer pvalue (toVector xs) (toVector ys) of
       Nothing             -> error "not enough data for mannWhitneyUtest"
       Just Significant    -> True
       Just NotSignificant -> False

-- Should we stop after the first trial of this package to save time? This
-- function skips the package if the results are uninteresting and the times are
-- within --min-run-time-percentage-difference-to-rerun.
shouldContinueAfterFirstTrial :: Double
                              -> NominalDiffTime
                              -> NominalDiffTime
                              -> CabalResult
                              -> CabalResult
                              -> Bool
shouldContinueAfterFirstTrial 0                            _  _  _       _       = True
shouldContinueAfterFirstTrial _                            _  _  Timeout Timeout = False
shouldContinueAfterFirstTrial maxRunTimeDifferenceToIgnore t1 t2 r1      r2      =
    isSignificantResult r1 r2
 || abs (t1 - t2) / min t1 t2 >= realToFrac (maxRunTimeDifferenceToIgnore / 100)

isSignificantResult :: CabalResult -> CabalResult -> Bool
isSignificantResult r1 r2 = r1 /= r2 || not (isExpectedResult r1)

-- Is this result expected in a benchmark run on all of Hackage?
isExpectedResult :: CabalResult -> Bool
isExpectedResult Solution       = True
isExpectedResult NoInstallPlan  = True
isExpectedResult BackjumpLimit  = True
isExpectedResult Timeout        = True
isExpectedResult Unbuildable    = True
isExpectedResult UnbuildableDep = True
isExpectedResult ComponentCycle = True
isExpectedResult ModReexpIssue  = True
isExpectedResult PkgNotFound    = False
isExpectedResult Unknown        = False

-- Combine CabalResults from multiple trials. Ignoring timeouts, all results
-- should be the same. If they aren't the same, we returns Unknown.
combineTrialResults :: [CabalResult] -> CabalResult
combineTrialResults rs
  | allEqual rs                          = head rs
  | allEqual [r | r <- rs, r /= Timeout] = Timeout
  | otherwise                            = Unknown
  where
    allEqual :: Eq a => [a] -> Bool
    allEqual xs = length (nub xs) == 1

timeEvent :: IO a -> IO (a, NominalDiffTime)
timeEvent task = do
  start <- getCurrentTime
  r <- task
  end <- getCurrentTime
  return (r, diffUTCTime end start)

diffTimeToDouble :: NominalDiffTime -> Double
diffTimeToDouble = fromRational . toRational

parserInfo :: ParserInfo Args
parserInfo = info (argParser <**> helper)
     ( fullDesc
    <> progDesc ("Find differences between two cabal commands when solving"
                   ++ " for all packages on Hackage.")
    <> header "hackage-benchmark" )

argParser :: Parser Args
argParser = Args
    <$> strOption
         ( long "cabal1"
        <> metavar "PATH"
        <> help "First cabal executable")
    <*> strOption
         ( long "cabal2"
        <> metavar "PATH"
        <> help "Second cabal executable")
    <*> option (words <$> str)
         ( long "cabal1-flags"
        <> value []
        <> metavar "FLAGS"
        <> help "Extra flags for the first cabal executable")
    <*> option (words <$> str)
         ( long "cabal2-flags"
        <> value []
        <> metavar "FLAGS"
        <> help "Extra flags for the second cabal executable")
    <*> option (map mkPackageName . words <$> str)
         ( long "packages"
        <> value []
        <> metavar "PACKAGES"
        <> help ("Space separated list of packages to test, or all of Hackage"
                   ++ " if unspecified"))
    <*> option auto
         ( long "min-run-time-percentage-difference-to-rerun"
        <> showDefault
        <> value 0.0
        <> metavar "PERCENTAGE"
        <> help ("Stop testing a package when the difference in run times in"
                   ++ " the first trial are within this percentage, in order to"
                   ++ " save time"))
    <*> option (mkPValue <$> auto)
         ( long "pvalue"
        <> showDefault
        <> value (mkPValue 0.05)
        <> metavar "DOUBLE"
        <> help ("p-value used to determine whether to print the results for"
                   ++ " each package"))
    <*> option auto
         ( long "trials"
        <> showDefault
        <> value 10
        <> metavar "N"
        <> help "Number of trials for each package")
    <*> switch
         ( long "print-trials"
        <> help "Whether to include the results from individual trials in the output")
    <*> switch
         ( long "print-skipped-packages"
        <> help "Whether to include skipped packages in the output")
    <*> option auto
         ( long "timeout"
        <> showDefault
        <> value 90
        <> metavar "SECONDS"
        <> help "Maximum time to run a cabal command, in seconds")

max_stars_repo_name: strake/cabal

max_issues_repo_name: strake/cabal

max_forks_repo_name: strake/cabal


text:
module Numeric.Tensile.Operations.Linear.Internal where

import Data.Vector.Sized (Vector)
import Numeric.Tensile.Dimensions
import Numeric.Tensile.Tensor.Internal
import Unsafe.Coerce (unsafeCoerce)

import qualified Data.Finite as F
import qualified Data.Vector.Sized as N
import qualified Data.Vector.Storable as S
import qualified Data.Vector.Storable.Mutable as M

-- import qualified Numeric.LinearAlgebra.HMatrix as Ma

-- f :: (Dims d' -> Perm n -> Perm n) -> Perm n -> Tensor d e -> Tensor d e'
-- f dim2Idx perm t = Tensor $ reifySomeDims (permuteDims perm (dims @d)) $ \p ->
--   modifyIdx (reflect p) (modify (permuteIdxs (dim2Idx (reflect p) perm) _)) (reflect p)) t -- basically make user derive the Idxs d' -> Idxs d'

-- dim2Idx :: Rank d ~ n => Dims d -> Perm n -> Perm n
-- takes a perm on dimensions and derives a perm in indices, eg
-- dim2Idx d p = lowerPerm' ...
-- --
-- otherwise consider using the raw index fold and lowerPerm???
-- could also create :  Perm d d'

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

transpose
  :: Elt e
  => Permutable d d'
  => Dims d -> Perm (Rank d) -> Tensor d e -> Tensor d' e
transpose d p (Tensor v) = Tensor v'
  where v' = modifyIdxs d v $ \i m ->
               remapIdxs p d i $ \d' i' ->
                 M.modify m (const $ v S.! fromIdxs d' (unsafePermute p i)) (fromIdxs d' i')

{-
gen#
  :: forall s t d. PrimBytes t
  => Int# -- ^ number of elements, not checked!
             --   Avoid using this argument if possible.
  -> (s -> (# s, t #))
  -> s -> (# s, ArrayBase t d #)
gen# n f z0 = go (byteSize @e undefined *# n)
  where
    go bsize = case runRW#
     ( \s0 -> case newByteArray# bsize s0 of
         (# s1, mba #) -> case loop0 mba 0# z0 s1 of
           (# s2, z1 #) -> case unsafeFreezeByteArray# mba s2 of
             (# s3, ba #) -> (# s3, (# z1, ba #) #)
     ) of (# _, (# z1, ba #) #) -> (# z1, ArrayBase (# | (# 0# , n , ba #) #) #)
    {-# NOINLINE go #-}
    loop0 mba i z s
      | isTrue# (i ==# n) = (# s, z #)
      | otherwise = case f z of
          (# z', x #) -> loop0 mba (i +# 1#) z' (writeArray mba i x s)

3-D tensor `a` w shape=[2, 2, 3]
[[[ 1,  2,  3],
  [ 4,  5,  6]],
 [[ 7,  8,  9],
  [10, 11, 12]]]

3-D tensor `b` w shape=[2, 3, 2]
[[[13, 14],
  [15, 16],
  [17, 18]],
 [[19, 20],
  [21, 22],
  [23, 24]]]
                )

mul `a` `b` has shape=[2,2,2]
[[[ 94, 100],
  [229, 244]],
 [[508, 532],
  [697, 730]]]
-}

{-

-- <#
-- same as tf.matmul
productN
  :: forall a b c x. All KnownDim '[a, b, c]
  => KnownDims x
  => T (x +: a +: b)
  -> T (x +: b +: c)
  -> T (x +: a +: c)
productN = undefined
  where
    d = dims @x
    --s = unsafeCoerce a

product
  :: forall m x y. KnownDim m
  => KnownDims x
  => KnownDims y
  => T (x +: m) -> T (m :+ y) -> T (x ++ y)
product t u
    | I# m <- fromIntegral $ dimVal @m
    , I# n <- fromIntegral $ totalDim' @x
    , I# k <- fromIntegral $ totalDim' @y
    , nk <- n *# k
    = let loop1 i j l r | isTrue# (l ==# m) = r
                        | otherwise = loop1 i j (l +# 1#)
                          (r + ix# (i +# n *# l) t * ix# (l +# m *# j) u)

          loop2 (T# i j) | isTrue# (j ==# k) = (# T# i j, 0 #)
                         | isTrue# (i ==# n) = loop2 (T# 0# (j +# 1#))
                         | otherwise = (# T# (i +# 1#) j, loop1 i j 0# 0 #)
      in case gen# nk loop2 (T# 0# 0#) of
          (# _, r #) -> r

data T# = T# Int# Int#

-- mode-i tensor-matrix product
-- see http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=4A0663C7848627DADDBA6A243BC43E78?doi=10.1.1.130.782&rep=rep1&type=pdf
product
  :: forall m x y. KnownDim m
  => KnownDims x
  => KnownDims y
  => T (x ++ [m] ++ y)
  -> T '[m, n]
  -> T (x ++ n ++ y)

-- #>
productR
  :: All KnownDim '[a, b, c]
  => KnownDims x
  => T (a :+ b :+ x)
  -> T (b :+ c :+ x)
  -> T (a :+ c :+ x)
productR = undefined

-- <#
-- same as tf.matmul
productN
  :: All KnownDim '[a, b, c]
  => KnownDims x
  => T (x +: a +: b)
  -> T (x +: b +: c)
  -> T (x +: a +: c)
productN = undefined

transpose
  :: forall m n x. All KnownDim '[m, n]
  => T '[n, m] --(n :+ m :+ x)
  -> T '[m, n] --(m :+ n :+ x)
transpose t = case elemSize0 t of
  0# -> broadcast (ix# 0# t)
  nm | I# m <- fromIntegral $ dimVal @m
     , I# n <- fromIntegral $ dimVal @n
     -> let f ( I# i,  I# j )
              | isTrue# (i ==# m) = f ( 0 , I# (j +# 1#) ) -- skip to next col
              | otherwise         = (# ( I# (i +# 1#), I# j ), ix# (i *# n +# j) t #) --col-major indexing
        in case gen# nm f (0,0) of
          (# _, r #) -> r

-}

max_stars_repo_name: cmk/tensile

max_issues_repo_name: cmk/tensile

max_forks_repo_name: cmk/tensile


text:
module Main where

import qualified Statistics.Matrix as M
import qualified Statistics.Matrix.Fast as F
import qualified Statistics.Matrix.Fast.Algorithms as FA
import           Statistics.Matrix (Matrix (..))
import qualified Statistics.Matrix.Algorithms as A

import qualified Data.Vector.Unboxed         as U
import           Data.Vector.Unboxed         (Vector)

import qualified System.Random.MWC as Mwc

import qualified Chronos.Bench as C

n :: Int
n = 100

testVector :: IO (Vector Double)
testVector = do
    gen <- Mwc.create
    Mwc.uniformVector gen (n*n)

testMatrix :: IO Matrix
testMatrix = do
    vec <- testVector
    return $ Matrix n n vec

runtimelight :: Vector Double -> Matrix -> IO ()
runtimelight v a = do
    let
      v2 = U.take n v

    C.defaultMainWith (C.defaultConfig {C.timeout = Just 3}) [
                   C.bench "norm" M.norm v2,
                   C.bench "Fast.norm" F.norm v2,

                   C.bench "multiplyV" (M.multiplyV a) (v2),
                   C.bench "Fast.multiplyV" (F.multiplyV a) (v2),

                   C.bench "transpose" M.transpose a,
                   C.bench "Fast.transpose" F.transpose a ,
                   C.bench "ident" M.ident n,
                   C.bench "diag" M.diag v2
                   ]

runtimeheavy :: Matrix -> Matrix -> IO ()
runtimeheavy a b = do

    C.defaultMainWith (C.defaultConfig {C.timeout = Just 1}) [
                   C.bench "multiply" (M.multiply a) b,
                   C.bench "Fast.multiply" (F.multiply a) b,
                   C.bench "qr" A.qr a,
                   C.bench "Fast.qr" FA.qr a
                   ]

main :: IO ()
main = do
    v <- testVector
    a <- testMatrix
    b <- testMatrix

    --
    putStrLn "---Benchmarking light operations---"
    -- we split heavy and light, we lose some precision in the bar plots from chronos
    runtimelight v a
    putStrLn "---Benchmarking heavy operations---"
    runtimeheavy a b

max_stars_repo_name: arvindd/dh-core

max_issues_repo_name: arvindd/dh-core

max_forks_repo_name: arvindd/dh-core


text:
-- | Recursive-descent r7rs Scheme parser, built on top of Alex
-- and an incremental stream architecture.
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TupleSections #-}
module Parser.Parser (readDatum, parseDatum, parseDatumSeq) where

import Control.Monad.Trans (lift)
import Control.Monad.ST.Trans
import Data.Complex (Complex(..))
import Data.Ratio (numerator, denominator, (%))
import Text.Parsec as P hiding (runParser, token, tokenPrim, eof)
import Text.Parsec.Error

import Val
import EvaluationMonad (panic)
import Parser.Monad hiding (token)
import Parser.Lexer

-- Unfortunately, we have to define the parts of the monad that go on top
-- of Alex here, and not in Monad.hs. We could get /most/ of the way there
-- by reorganizing tokens and lexemes, but we can't move 'alexMonadScan'
-- there because it requires things defined in Lexer.hs.
-- I guess in an ideal world, we'd define all this stuff in Monad.hs and
-- get it to compile with a Lexer.hs-boot file.

-- | A lazily-constructed, cached, incremental 'Lexeme' stream type.
-- Haskell's own laziness cannot (obviously, perhaps there is a way)
-- be used to construct a stream of 'Lexeme's because we do not know when
-- the stream will end until the parser decides that it is done parsing.
newtype S s = S (STRef s (Maybe (Maybe (Lexeme, S s))))

unconsS :: S s -> M s (Maybe (Lexeme, S s))
unconsS (S ref) = do
  cache <- readSTRef ref
  case cache of
    Just res -> return res
    Nothing -> do
      lx <- lift alexMonadScan
      restS <- S <$> newSTRef Nothing
      -- this will put a literal TokEOF in the stream at the end.
      let res = Just (lx, restS)
      writeSTRef ref $ Just res
      return res

instance Stream (S s) (M s) Lexeme where
  uncons = unconsS

-- | The monad underlying ParsecT.
type M s = STT s Alex

-- | The parsec parsing monad.
type Parser s = ParsecT (S s) () (M s)

-- | Run a parser, given a label and an input port. The output is:
-- Left (message, first-token-was-eof?)
-- Right parsed-value
runParser :: (forall s. Parser s a) -> String -> Port -> Either (String, Bool) a
runParser p label inputPort
  = case runAlex' alex label inputPort of
    Right (Left (parseError, b)) -> Left (show parseError, b)
    Right (Right v) -> Right v
    Left lexError -> Left (lexError, False)
  where
    alex = runSTT $ do
      initialS <- S <$> newSTRef Nothing
      let parse = runParserT p () label initialS
      firstUncons <- unconsS initialS
      case firstUncons of
        Just (L _ TokEOF, _) -> mapLeft (,True)  <$> parse
        _other               -> mapLeft (,False) <$> parse

runParser_ :: (forall s. Parser s a) -> String -> Port -> Either String a
runParser_ p lbl port = mapLeft fst $ runParser p lbl port

mapLeft :: (e1 -> e2) -> Either e1 b -> Either e2 b
mapLeft f (Left  x) = Left (f x)
mapLeft _ (Right x) = Right x

liftAlex :: Alex a -> Parser s a
liftAlex = lift . lift

-- ----------------------------------------------------------------------------
-- Entrypoints

-- | Read the next datum parseable from the given port.
readDatum :: String -> Port -> Either (String, Bool) Val
readDatum = runParser datum

-- | Parse an input port that contains exactly one datum.
parseDatum :: String -> Port -> Either String Val
parseDatum = runParser_ $ datum <* eof

parseDatumSeq :: String -> Port -> Either String [Val]
parseDatumSeq = runParser_ $ datumSeq <* eof
  -- enforce that we reach the end of the stream, since we only
  -- use this entrypoint to parse whole source files.

-- ----------------------------------------------------------------------------
-- Token parsers

-- | this is like 'Text.Parsec.Prim.tokenPrim', except that 'posFromTok' should
-- give the position of the _matched_ token, not the next one. We need
-- something like this because we absolutely cannot allow a lookahead token.
-- This has the side effect of leaving Parsec's internally stored position
-- on the previously-lexed token, instead of the start of the next one like
-- parsec is expecting. This is /more reasonable/ behavior, at least in our
-- case, and results in better error messages in 'validateByte'.
--
-- The previous 'SourcePos' is available to the position-extraction function
-- so that it can build a 'SourcePos' with parsec's 'setXXX' functions, since
-- the constructor is not exposed. It's unlikely that you'll actually need to
-- use the previous 'SourcePos' in a meaningful way.
tokenPrim :: Stream s m t
          => (t -> String) -> (SourcePos -> t -> SourcePos) -> (t -> Maybe a)
          -> ParsecT s u m a
tokenPrim showTok posFromTok test = mkPT $
  \(State input pos usr) -> do
    r <- uncons input
    pure $ case r of
      Nothing -> Empty $ pure $ P.Error $ unexpectError "" pos
      Just (c,rest) -> let newpos = posFromTok pos c
                       in case test c of
        Nothing -> Empty $ pure $ P.Error $ unexpectError (showTok c) newpos
        Just x -> let newState = State rest newpos usr
                  in newpos `seq` newState `seq`
                     Consumed $ pure $ Ok x newState (newErrorUnknown newpos)
  where
    unexpectError msg pos = newErrorMessage (SysUnExpect msg) pos

token :: (Token -> Maybe a) -> Parser s a
token match = tokenPrim displayLexeme movePos match'
  where
    movePos oldPos (L (AlexPn _addr line col) _)
      = setSourceColumn (setSourceLine oldPos line) col

    match' (L _ tok) = match tok

tokenWithPos :: (Token -> Maybe a) -> Parser s (a, AlexPosn)
tokenWithPos match = tokenPrim displayLexeme movePos match'
  where
    movePos oldPos (L (AlexPn _addr line col) _)
      = setSourceColumn (setSourceLine oldPos line) col

    match' (L pos tok) = (,pos) <$> match tok

displayLexeme :: Lexeme -> String
displayLexeme (L _ t) = "'" ++ show t ++ "'"

eof :: Parser s ()
eof = token $ \case
  TokEOF -> Just ()
  _other -> Nothing

simpleDatum :: Parser s Val
simpleDatum = token $ \case
  TokSymbol sym -> Just $ Symbol $ pack sym
  TokNumber num -> Just $ mkNumber num
  TokBool b     -> Just $ Bool b
  TokChar c     -> Just $ Char c
  TokString str -> Just $ IString $ pack str
  _other        -> Nothing

bracketed :: Parser s a -> Parser s a
bracketed inside = do
  rightBracket <- leftBracket
  r <- inside
  rightBracket
  return r
  where
    leftBracket = token $ \case
      TokLParen -> mkRight $ \case
        TokRParen -> Just ()
        _other    -> Nothing
      TokLBracket -> mkRight $ \case
        TokRBracket -> Just ()
        _other      -> Nothing
      TokLBrace -> mkRight $ \case
        TokRBrace -> Just ()
        _other    -> Nothing
      _other -> Nothing
    mkRight = return . token

sugarToken :: Parser s Token
sugarToken = token $ \case
  TokQuote -> Just TokQuote
  TokComma -> Just TokComma
  TokCommaAt -> Just TokCommaAt
  TokBackquote -> Just TokBackquote
  TokHashQuote -> Just TokHashQuote
  TokHashBackquote -> Just TokHashBackquote
  TokHashComma -> Just TokHashComma
  TokHashCommaAt -> Just TokHashCommaAt
  _other -> Nothing

dot :: Parser s ()
dot = token $ \case
  TokDot -> Just ()
  _other -> Nothing

datumCommentStart :: Parser s ()
datumCommentStart = do
  token $ \case
    TokDatumComment -> Just ()
    _other -> Nothing
  liftAlex enterDatumComment

datumCommentEnd :: Parser s ()
datumCommentEnd = liftAlex exitDatumComment

startVector :: Parser s ()
startVector = token $ \case
  TokLVector -> Just ()
  _other -> Nothing

startByteVector :: Parser s ()
startByteVector = token $ \case
  TokLByteVector -> Just ()
  _other -> Nothing

startLabelDef :: Parser s Int
startLabelDef = do
  (lbl, pos) <- tokenWithPos $ \case
    TokLabelDef lbl -> Just lbl
    _other -> Nothing
  liftAlex $ startDef pos lbl
  return lbl

labelReference :: Parser s Val
labelReference = do
  lbl <- token $ \case
    TokLabelRef lbl -> Just lbl
    _other -> Nothing
  liftAlex $ lookupLabel lbl

-- | Parses a label reference, but additionally ensures
-- that the parser's internal 'defining_label' field does
-- not match the parsed label. If it does, raises an error
-- (at the Alex level).
safeLabelReference :: Parser s Val
safeLabelReference = do
  lbl <- token $ \case
    TokLabelRef lbl -> Just lbl
    _other -> Nothing
  liftAlex $ testRef lbl >> lookupLabel lbl

-- ----------------------------------------------------------------------------
-- Core parser

-- trying to be a little bit more aggressive about having a simple grammar
-- with parsec, since we don't have to satisfy shift/reduce conflicts.
-- We still have to stratify label definitions, to check well-definedness.

datum :: Parser s Val
datum = simpleDatum
        <|> compoundDatum
        <|> definedDatum
        <|> labelReference
        <|> (datumComment >> datum)
        <?> "datum"

-- This parser is not (fully) incremental. It will give up on EOF,
-- of course, but if the input stream contains only some data
-- followed by non-data text, it will try to read one "token"
-- out of the non-data text before returning. This is unavoidable,
-- otherwise it can't know when to quit.
datumSeq :: Parser s [Val]
datumSeq = many datum

datumComment :: Parser s ()
datumComment = do
  datumCommentStart
  _ <- datum
  datumCommentEnd

compoundDatum :: Parser s Val
compoundDatum = (prefixSugar <$> sugarToken <*> datum)
                <|> list
                <|> vector

list :: Parser s Val
list = labeled "list" $ bracketed listContents

listContents :: Parser s Val
listContents = do
  prefix <- many datum
  maybeDot <- optionMaybe $ dot >> datum
  case maybeDot of
    Nothing -> return $ makeImmutableList prefix
    Just d  -> return $ makeImproperImmutableList prefix d

vector :: Parser s Val
vector = normalVector <|> byteVector

normalVector :: Parser s Val
normalVector = labeled "vector" $ do
  startVector
  contents <- many datum
  closeParen
  return $ makeVector contents

byteVector :: Parser s Val
byteVector = labeled "bytevector" $ do
  startByteVector
  bytes <- many byte
  closeParen
  return $ makeByteVector bytes

closeParen :: Parser s ()
closeParen = token $ \case
  TokRParen -> Just ()
  _other -> Nothing

byte :: Parser s Byte
byte = labeled "byte" $ do
  (num, pos) <- tokenWithPos $ \case
    TokNumber n -> Just n
    _other -> Nothing
  validateByte pos num

definedDatum :: Parser s Val
definedDatum = do
  lbl <- startLabelDef
  _ <- fail $ "At the moment, cyclic references will hang, so Thcheme"
              ++ " refuses all label definitions."
  val <- defDatum
  liftAlex $ do
    endDef
    unlessDatumComment $ defineLabel lbl val
    return val

defDatum :: Parser s Val
defDatum = simpleDatum
           <|> safeLabelReference
           <|> (prefixSugar <$> sugarToken <*> datum)
           <|> vector
           <|> defList

-- | like 'list', but if the list contains only a dot, and that
-- dot is a reference to a label that is defined as the list itself,
-- then raises an error as the label is ill-defined.
defList :: Parser s Val
-- Note 'listContents -> . d' is a production of the grammar, but it
-- doesn't matter as long as we test for that case first.
defList = labeled "list" $ bracketed $ (dot >> defDatum) <|> listContents

-- ----------------------------------------------------------------------------
-- Utilities for constructing/validating Vals

labeled :: String -> Parser s a -> Parser s a
labeled = flip label

mkNumber :: Numeric -> Val
mkNumber (LitInteger i)  = Number $ Real $ Bignum i
mkNumber (LitRational r) = Number $ Real $ Ratnum r
mkNumber (LitFloating f) = Number $ Real $ Flonum f
mkNumber (LitComplex l r) = Number $ Complex $
  numericToNumber l :+ numericToNumber r

numericToNumber :: Numeric -> RealNumber
numericToNumber (LitInteger i)  = Bignum i
numericToNumber (LitRational r) = Ratnum r
numericToNumber (LitFloating f) = Flonum f
numericToNumber LitComplex{} = panic "numericToNumber: nested complex"

prefixSugar :: Token -> Val -> Val
prefixSugar tok v = makeImmutableList [Symbol (desugar tok), v]
  where
    desugar TokQuote         = "quote"
    desugar TokBackquote     = "quasiquote"
    desugar TokComma         = "unquote"
    desugar TokCommaAt       = "unquote-splicing"
    desugar TokHashQuote     = "syntax"
    desugar TokHashBackquote = "quasisyntax"
    desugar TokHashComma     = "unsyntax"
    desugar TokHashCommaAt   = "unsyntax-splicing"
    desugar _ = panic "desugar: not a sugar token!"

validateByte :: AlexPosn -> Numeric -> Parser s Byte
validateByte p n = case n of
  LitInteger i -> checkSize i
  LitRational r
    | denominator r == 1 -> checkSize $ numerator r
    | otherwise -> boom
  LitComplex real imag -> do checkExactZero imag
                             validateByte p real
  _other -> boom
  where
    checkSize :: Integer -> Parser s Byte
    checkSize i
      | i >= 0 && i <= 255 = return $ fromInteger i
      | otherwise = boom

    checkExactZero :: Numeric -> Parser s ()
    checkExactZero (LitInteger 0) = pure ()
    checkExactZero (LitRational r)
      | r == 0 % 1 = pure ()
      | otherwise = boom
    checkExactZero _ = boom

    boom :: Parser s a
    boom = fail $ show n ++ " is not an exact byte value"

max_stars_repo_name: JKTKops/Thcheme

max_issues_repo_name: JKTKops/Thcheme

max_forks_repo_name: JKTKops/Thcheme


text:
{-# LANGUAGE TemplateHaskell, RankNTypes #-}

module Geometry.LRS (
    -- * Types
    Dictionary,
    -- * Functions
    getDictionary,
    selectPivot,
    reverseRS,
    lrs,
    polytope2LiftedCone,
    liftedCone2polyotpe,
    lrsFacetEnumeration,

    colFromList
    )
 where

import Geometry.Vertex
import Geometry.Matrix

import Data.Matrix as M hiding (trace)
import qualified Data.Vector as V
import Data.Maybe
import Data.List
import Util
import Control.Lens
import Debug.Trace
import Numeric.LinearProgramming hiding (simplex)
import Numeric.LinearAlgebra.HMatrix as HM
import Data.List

type Row = M.Matrix Rational
type Col = M.Matrix Rational

data Dictionary = Dict {
                            __B :: [Int],
                            __N :: [Int],
                            _dict :: M.Matrix Rational
                        }
    deriving (Show, Eq)

makeLenses ''Dictionary

numRows :: Dictionary -> Int
numRows dictionary = nrows $ dictionary^.dict

numCols :: Dictionary -> Int
numCols dictionary = ncols $ dictionary^.dict

(|*|) = multStrassen

colFromList :: [a] -> M.Matrix a
colFromList = colVector . V.fromList

rowFromList :: [a] -> M.Matrix a
rowFromList = rowVector . V.fromList

-- | Wrapper of submatrix function. This function works for indices starting from 0 rather than 1 in the submatrix function
submatrix' ::
    (Int, Int)      -- ^ Rows indices
    -> (Int, Int)   -- ^ Cols indices
    -> M.Matrix a
    -> M.Matrix a
submatrix' (ro,rk) (co,ck) m = submatrix (ro+1) (rk+1) (co+1) (ck+1) m

-- | Lens for columns
colAt :: Int -> Lens' (M.Matrix Rational) Col
colAt j = lens (getCol' j) (\m c -> setCol' j c m)

getCol' :: Int -> M.Matrix a -> M.Matrix a
getCol' idx = colVector . (getCol (idx+1))

setCol' :: Int -> Col -> M.Matrix Rational -> M.Matrix Rational
setCol' idx col mat
    | idx == 0 = col <|> right
    | idx == (ncols mat) - 1 = left <|> col
    | otherwise = left <|> col <|> right
    where
        left = submatrix' (0,(nrows mat)-1) (0,idx-1) mat
        right = submatrix' (0,(nrows mat)-1) (idx+1, (ncols mat)-1) mat

-- | Lens for rows
rowAt :: Int -> Lens' (M.Matrix Rational) Row
rowAt i = lens (getRow' i) (\m r -> setRow' i r m)

getRow' :: Int -> M.Matrix a -> M.Matrix a
getRow' idx = rowVector . (getRow (idx+1))

setRow' :: Int -> Col -> M.Matrix Rational -> M.Matrix Rational
setRow' idx row mat
    | idx == 0 = row <-> down
    | idx == (nrows mat) - 1 = up <-> row
    | otherwise = up <-> row <-> down

    where
        up = submatrix' (0,idx-1) (0,(ncols mat)-1)  mat
        down = submatrix' (idx+1, (nrows mat)-1) (0,(ncols mat)-1) mat

-- | Lens for accessing elements

elemAt :: (Int, Int) -> Lens' (M.Matrix a) a
elemAt (i, j) = lens (getElem' (i,j)) (\m x -> setElem' x (i,j) m)

getElem' :: (Int, Int) -> M.Matrix a -> a
getElem' (row, col) = getElem (row+1) (col+1)

setElem' :: a -> (Int, Int) -> M.Matrix a -> M.Matrix a
setElem' elem (row, col) mat
    | row < 0 || col < 0 || row >= nrows mat || col >= ncols mat = error ("setElem': Trying to set at (" ++ show row ++ "," ++ show col ++ ") in a " ++ show (nrows mat) ++ "x" ++ show (ncols mat) ++ " matrix.")
    | otherwise = setElem elem (row+1, col+1) mat

mapCol' :: (Int -> a -> a) -> Int -> M.Matrix a -> M.Matrix a
mapCol' f col m
    | col < 0 || col > (ncols m)-1 = error "mapCol': index out of bounds"
    | otherwise = mapCol f (col+1) m

mapRow' :: (Int -> a -> a) -> Int -> M.Matrix a -> M.Matrix a
mapRow' f row m
    | row < 0 || row > (nrows m)-1 = error "mapCol': index out of bounds"
    | otherwise = mapRow f (row+1) m

getOptimumVertex :: M.Matrix Rational -> Col -> Maybe Vertex
getOptimumVertex mat col = fmap (map toRational) $ feasNOpt result
    where
        problem = Maximize $ replicate (ncols mat) (-1)
        constraints = Dense $ safeZipWith (:<=:) (map (map fromRational) $ M.toLists mat) (map fromRational $ concat $ M.toLists col)
        result = exact problem constraints (map Free [1..ncols mat])
        feasNOpt (Optimal (_, vertex)) = Just vertex
        feasNOpt _ = Nothing

{-
    Dictionary form

                p21     p31
    invA_b
                p22     p32

 -}

toHMatrix :: M.Matrix Rational -> HM.Matrix Double
toHMatrix matrix = HM.matrix cols (map (fromRational) (concat $ M.toLists matrix))
    where
        cols = ncols matrix

findBasis :: (Maybe (M.Matrix Rational), Maybe (M.Matrix Rational)) -> M.Matrix Rational -> Int -> (Maybe (M.Matrix Rational), Maybe (M.Matrix Rational))
findBasis (basic, cobasic) original colNum
    | currRank == rows && colNum < cols = findBasis (basic, Just newCobasic) original (colNum + 1)-- Just ((fromJust cobasic) <|> (submatrix' (0,rows-1) (colNum, (ncols original)-1) original)) )
    | currRank == rows = (basic, cobasic)-- Just ((fromJust cobasic) <|> (submatrix' (0,rows-1) (colNum, (ncols original)-1) original)) )
    | currRank == candRank = findBasis (basic, Just newCobasic) original (colNum + 1)
    | currRank < candRank = findBasis (Just newBasic, cobasic) original (colNum + 1)
    where
        rows = nrows original
        cols = ncols original
        currRank = if isNothing basic then 0 else rank (toHMatrix $ fromJust basic)
        candVec = original ^. colAt colNum
        candMat = if isNothing basic then candVec else (fromJust basic) <|> candVec
        candRank = rank $ toHMatrix candMat

        newBasic = if isNothing basic then candVec else (fromJust basic) <|> candVec
        newCobasic = if isNothing cobasic then candVec else (fromJust cobasic) <|> candVec

sortSystem :: M.Matrix Rational -> (M.Matrix Rational, M.Matrix Rational)
sortSystem matrix = trace ("Matrix : "++ show matrix ++ "\n\nSorted " ++ show (findBasis (Nothing, Nothing) matrix 0) ) (((,) <$> (fromJust . fst) <*> (fromJust . snd)) (findBasis (Nothing, Nothing) matrix 0))
-- sortSystem :: M.Matrix Rational -> Col -> Vertex -> (M.Matrix Rational, Col)

-- sortSystem mat col vertex
--     | cols /=  (rows `div` 2) = (,) mat col
--     | otherwise = (trace $ show $ M.fromLists (upper ++ lower)) (,) (M.fromLists (upper ++ lower)) col
--     where
--         matLists = M.toLists mat

--         rows = nrows mat
--         cols = ncols mat
--         slackMatr = submatrix' (0,rows-1) (0, rows - cols -2) $ identity rows
--         lower  = getIndependent matLists [] 0 rows slackMatr
--         upper = matLists \\ lower

-- getDictionary :: Matrix Rational -> Col -> Vertex -> Dictionary
-- getDictionary _A b vertex = trace ("In dictionary" ++ show (rows,cols, nrows dictionary, ncols dictionary) ++ show newDict) Dict [0..rows] [rows+1..rows+cols] newDict
--     -- trace ("In dictionary" ++ show _A_B)
--     where
--         (newA, newb) = sortSystem _A b vertex
--         rows = nrows newA
--         cols = ncols newA
--         slack = identity rows
--        -- dictionary = newA <|> slack
--         topRow = rowFromList $ 1 : replicate rows 0 ++  replicate cols (1)
--         leftCol = colFromList $ replicate rows 0
--         dictionary = topRow <-> (leftCol <|> slack <|> newA )

--         -- c_B = submatrix' (0,0) (1,rows) topRow
--         -- c_N = submatrix' (0,0) (rows+1, rows+cols) topRow
--         -- _A_B = submatrix' (0,rows-1) (0,rows-1) dictionary
--         -- _A_N = submatrix' (0,rows-1) (rows, rows+cols-1) dictionary
--         _A_B = submatrix' (0,rows) (0,rows) dictionary
--         _A_N = submatrix' (0,rows) (rows+1, rows+cols) dictionary
--         Right invA_B = inverse _A_B

--         p2 = - invA_B |*| _A_N
--         p3 = invA_B |*| newb
--         -- p21 = (c_N - c_B |*| invA_B |*| _A_N)
--         -- p22 = invA_B |*| _A_N
--         -- p31 = -c_B |*| invA_B |*| newb
--         -- p32 = invA_B |*| newb
--         newDict = ((identity (rows+1)) <|> p2 <|> p3)

-- | Constructs initial dictionary given a set of inequalities and an optimum initial vertex
-- getDictionary :: M.Matrix Rational -> Col -> Vertex -> Dictionary
-- getDictionary _A b vertex = Dict [0..rows] [rows+1..rows+cols] newDict
--     where
--         -- (newA, newb) = sortSystem _A b vertex
--         initialDict = _A <|> slack
--         (basic, cobasic) = sortSystem initialDict
--         (newA, newb) = (,) _A b
--         rows = nrows newA
--         cols = ncols newA
--         slack = identity rows
--         dictionary = newA <|> slack
--         topRow = rowFromList $ 1 :  replicate rows 0 ++  replicate cols 1
--         c_B = submatrix' (0,0) (1,rows) topRow
--         c_N = submatrix' (0,0) (rows+1, rows+cols) topRow
--         _A_B = submatrix' (0,rows-1) (0,rows-1) dictionary
--         _A_N = submatrix' (0,rows-1) (rows, rows+cols-1) dictionary
--         Right invA_B = inverse _A_B
--         p21 = (c_N - c_B |*| invA_B |*| _A_N)
--         p22 = invA_B |*| _A_N
--         p31 = -c_B |*| invA_B |*| newb
--         p32 = invA_B |*| newb

--         newDict = ((identity (rows+1)) <|> (p21 <-> p22) <|> (p31 <-> p32))

getDictionary :: M.Matrix Rational -> Col -> Vertex -> Dictionary
getDictionary _A b vertex = Dict [0..rows] [rows+1..rows+cols] newDict
    where

        -- (newA, newb) = sortSystem _A b vertex
        initialDict = _A <|> slack
        (basic, cobasic) = sortSystem initialDict
        (newA, newb) = (,) _A b
        rows = nrows newA
        cols = ncols newA
        slack = identity rows
        -- dictionary = newA <|> slack
        topRow = rowFromList $ 1 :  replicate rows 0 ++  replicate cols 1
        c_B = submatrix' (0,0) (1,rows) topRow
        c_N = submatrix' (0,0) (rows+1, rows+cols) topRow
        _A_B = basic --submatrix' (0,rows-1) (0,rows-1) dictionary
        _A_N = cobasic --submatrix' (0,rows-1) (rows, rows+cols-1) dictionary
        Right invA_B = inverse _A_B
        p21 = (c_N - c_B |*| invA_B |*| _A_N)
        p22 = invA_B |*| _A_N
        p31 = -c_B |*| invA_B |*| newb
        p32 = invA_B |*| newb

        newDict = ((identity (rows+1)) <|> (p21 <-> p22) <|> (p31 <-> p32))

enteringVariable :: Dictionary -> Maybe Int
enteringVariable dictionary
    | null negs = Nothing
    | otherwise = trace ("ENTERINGVAR\nBASE: " ++ show (dictionary^._B) ++ "\nCOBASE: " ++ show (dictionary^._N))  (Just $ (fst.head.sort) negs)
    where
        rows = nrows $ dictionary^.dict
        cobasic_0 = zip (dictionary^._N) (map (\j -> dictionary^.dict.elemAt (0,j)) [rows..])
        negs = filter (\(_,value) -> value < 0) cobasic_0

lexMinRatio :: Dictionary -> Int -> Int
lexMinRatio dictionary s
    | dim+1 == rows = 0
    | null indexed_s = 0
    | otherwise = trace ("LEXMINRATIO\nBASE: " ++ show (dictionary^._B) ++ "\nCOBASE: " ++ show (dictionary^._N))  (dictionary ^. _B) !! (fst $ indexed_s !! (fromJust $ elemIndex (minimum ratios) ratios))
    where
        rows = numRows dictionary
        cols = numCols dictionary
        dim = cols - rows - 1
        dictMatrix = dictionary^.dict
        col_s = dictMatrix ^. colAt s
        _D = (dictMatrix ^. colAt (cols-1)) <|> submatrix' (0,rows-1) (0, rows-1) dictMatrix
        slice_s = (concat.M.toLists) $ submatrix' (dim+1, rows-1) (0,0) col_s
        indexed_s = filter (snd.(fmap (>0))) $ safeZipWith (,) [dim+1..rows-1] slice_s
        sub_D = map (head . M.toLists . (\i -> _D ^. rowAt i) . fst) indexed_s
        ratios = safeZipWith (map $) (map ((flip (/)) . snd) indexed_s) sub_D

-- | Selects entering and leaving variables for the next pivot in the dictionary
selectPivot :: Dictionary -> Maybe (Int, Int)
selectPivot dictionary
    | s == Nothing = Nothing
    | otherwise = Just (r, fromJust s)
    where
        s = enteringVariable dictionary
        r = lexMinRatio dictionary (fromJust s)

pivot :: Int ->Int -> Dictionary -> Dictionary
pivot r s dictionary = Dict new_B new_N (_E |*| dictMatrix)
    where
        dictMatrix = dictionary ^. dict
        col_s = dictMatrix ^. colAt (s)
        t = fromJust $ elemIndex r (dictionary ^. _B) -- idxR
        idxS = fromJust $ elemIndex s (dictionary ^. _N)
        a_t = col_s ^. elemAt (t,0)
        eta = (mapCol' (\_ x -> -x/a_t) 0 col_s) & elemAt (t,0) .~ (1/a_t)
        _E = (identity (nrows dictMatrix)) & colAt t .~ eta
        new_B = dictionary^._B & element t .~ s
        new_N = dictionary^._N & element idxS .~ r

simplex :: Dictionary -> Dictionary
simplex dictionary
    | idxsToPivot == Nothing = dictionary
    | otherwise = simplex $ uncurry pivot (fromJust idxsToPivot) dictionary
    where
        idxsToPivot = selectPivot dictionary

-- | Reverse search.
reverseRS ::
        Dictionary
    ->  Int          -- ^ element in N
    ->  Maybe Int
reverseRS dictionary v
    | conditions == False = trace ("REVERSE RS: " ++ show [firstCondition, sndCondition, lastCondition]) Nothing
    | conditions == True = trace ("REVERSE RS: " ++ show v) Just u
    where
        dictMatrix = dictionary^.dict
        rows = nrows dictMatrix
        v_col = dictMatrix ^. colAt (v)
        w_row_0 = dictMatrix ^. rowAt 0
        u = lexMinRatio dictionary v
        i = fromJust $ elemIndex u (dictionary ^. _B)
        w_row_i = mapRow' (\_ x -> (v_col ^. elemAt (0,0))/(v_col ^. elemAt (i,0)) * x) 0 $ dictMatrix ^. rowAt (i)
        diff_ws = (head.M.toLists) $ w_row_0 - w_row_i
        firstCondition = (w_row_0 ^. elemAt (0,v)) > 0
        sndCondition = u /= 0
        lastCondition = all (>=0) [(diff_ws!!(snd idx))|idx <- (zip (dictionary^._N) [rows..]), (fst idx) < u ]
        conditions = firstCondition && sndCondition && lastCondition -- Proposition 6.1

-- reverseRS :: -- reverse with pivot and selectPivot
--     Dictionary ->
--     Int ->
--     Maybe Int
-- reverseRS dictionary v
--     | newPivots == Nothing = Nothing
--     | condition == False = Nothing
--     | condition == True = Just u
--     where
--         u = lexMinRatio dictionary v
--         prev_B = pivot u v dictionary
--         newPivots = selectPivot prev_B
--         condition = fst (fromJust newPivots) == v && (snd (fromJust newPivots)) == u

getVertex :: Dictionary -> Vertex
getVertex dictionary = concat $ M.toLists $ submatrix' (1,dim) (cols-1, cols-1) (dictionary^.dict)
    where
        rows = numRows dictionary
        cols = numCols dictionary
        dim = cols - rows - 1

revSearch :: Dictionary -> [Vertex]
revSearch dictionary@(Dict _B _N dictMatrix) -- = getVertex dictionary : concatMap revSearch pivoted
    | (not.null) possibleRay = possibleRay ++ (concatMap revSearch pivoted)
    | otherwise = getVertex dictionary : concatMap revSearch pivoted
    where
        rows = numRows dictionary
        cols = numCols dictionary
        valid_N = [i | i <- _N, (reverseRS dictionary i) /= Nothing]
        valid_B = map (lexMinRatio dictionary) valid_N
        pivoted = map (\(r, s) ->  pivot r s dictionary) $ zip valid_B valid_N
        possibleRay = hasRay dictionary

-- lexMin :: Dictionary -> Maybe Extremal
-- lexMin dictionary

hasRay :: Dictionary -> [Vertex]
hasRay dictionary = trace ("RAYS: " ++ show (map snd rays) ++ "\n\n" ++ show (map fst rays) ++ "\n" ++ show (dictionary ^. dict) ++ "\n\n") (map snd rays)
    -- idxsPivot == Nothing = Nothing
    -- r == 0 = Just ray
    -- otherwise = Nothing
    where
        dictMatrix = dictionary ^. dict
        rows = numRows dictionary
        dim = cols - rows - 1
        cols = numCols dictionary
        nonPositive column = all (<=0) ((concat.M.toLists) column)
        colsWithRays = if dim+1 == rows then
                            [(j, dictMatrix^.colAt j) | j <- [rows..cols-2]]
                        else [(j, dictMatrix^.colAt j) | j <- [rows..cols-2], nonPositive (submatrix' (dim+1,rows-1) (j,j) dictMatrix )]
        rays =  zip (map fst colsWithRays) (map (concat . M.toLists . (submatrix' (1,dim) (0,0)) . snd ) colsWithRays)

-- | Implementation of Lexicographic Reverse Search
lrs :: M.Matrix Rational -> Col -> Vertex-> [Vertex]
lrs matrix b vertex = (sort.nub) $ revSearch dictionary
    where
        dictionary = simplex $ getDictionary matrix b vertex

-- | Transforms a V polytope into a H pointed cone
polytope2LiftedCone :: [Vertex] -> (M.Matrix Rational, Col)
polytope2LiftedCone points = (colOnes <|> hMatrix , colZeros)
    where
        nrows = length points
        colOnes = colFromList $ replicate nrows 1
        hMatrix = M.fromLists points
        colZeros = colFromList $ replicate nrows 0

-- | Transforms a V pointed cone into a H polyotpe
liftedCone2polyotpe :: [Vertex] -> (M.Matrix Rational, Col)
liftedCone2polyotpe points = (hMatrix, b)
    where
        dim = length (points !! 0)
        rays  = points \\ [(replicate dim 0)] -- remove vertex at the origin
        b = colFromList $ map (negate.head) rays
        hMatrix = M.fromLists $ map tail rays

-- | Implementation of Lexicographic Reverse Search for Facet Enumeration
lrsFacetEnumeration :: [IVertex] -> (M.Matrix Rational, Col)
lrsFacetEnumeration ipoints = hPolytope
    where
        points =  map (map toRational) ipoints
        cols = length (points!!0)
        initialVertex = replicate (cols + 1) 0
        hLiftedCone@(matrix, b) = polytope2LiftedCone points
        vLiftedCone = lrs matrix b initialVertex
        hPolytope = liftedCone2polyotpe vLiftedCone

max_stars_repo_name: tropical-group/facet-enumeration

max_issues_repo_name: tropical-group/facet-enumeration

max_forks_repo_name: tropical-group/facet-enumeration


text:
{-# LANGUAGE Arrows              #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies        #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS -fplugin=Overloaded -fplugin-opt=Overloaded:Categories #-}
module Main where

import Control.Monad      (when)
import Data.Word          (Word64)
import Numeric            (showFFloat)
import System.Environment (getArgs)
import Data.List (intercalate)

import qualified Control.Category
import qualified Numeric.LinearAlgebra  as LA
import qualified System.Random.SplitMix as SM

import Overloaded.Categories
import VectorSpace

-- | A Function which computes value and derivative at the point.
newtype AD a b = AD (a -> (b, L a b))

instance Category AD where
    id = AD (\x -> (x, L id))

    AD g . AD f = AD $ \a ->
        let (b, L f') = f a
            (c, L g') = g b
        in (c, L (g' . f'))

instance CategoryWith1 AD where
    type Terminal AD = ()

    terminal = AD (const ((), terminal))

instance CartesianCategory AD where
    type Product AD = (,)

    proj1 = AD (\x -> (fst x, proj1))
    proj2 = AD (\x -> (snd x, proj2))

    fanout (AD f) (AD g) = AD $ \a ->
        let (b, f') = f a
            (c, g') = g a
        in ((b, c), fanout f' g')

instance GeneralizedElement AD where
    type Object AD a = a

    konst x = AD (\_ -> (x, L $ \_ -> LZ))

ladd :: LinMap r (a, a) -> LinMap r a
ladd (LH f g) = LA f g
ladd (LV f g) = LV (ladd f) (ladd g)
ladd (LA a b) = LA (ladd a) (ladd b)
ladd LZ       = LZ
ladd (LD k)   = LV (LD k) (LD k)

lmult :: Double -> Double -> LinMap r (a, a) -> LinMap r a
lmult x y (LH f g) = LA (lmul y f) (lmul x g)
lmult x y (LV f g) = LV (lmult x y f) (lmult x y g)
lmult x y (LA f g) = LA (lmult x y f) (lmult x y g)
lmult _ _ LZ       = LZ
lmult x y (LD k)   = LV (LD (k * y)) (LD (k * x))

plus :: AD (Double, Double) Double
plus = AD $ \(x,y) -> (x + y, L ladd)

minus :: AD (Double, Double) Double
minus = AD $ \(x,y) -> (x - y, L $ lmult (-1) 1)

mult :: AD (Double, Double) Double
mult = AD $ \(x,y) -> (x * y, L $ lmult x y)

scale :: Double -> AD Double Double
scale k = AD $ \x -> (k * x, linear k)

evaluateAD :: (HasDim a, HasDim b) => AD a b -> a -> (b, LA.Matrix Double)
evaluateAD (AD f) x = let (y, f') = f x in (y, evalL f')

-------------------------------------------------------------------------------
-- Simple examples
-------------------------------------------------------------------------------

ex1 :: AD Double Double
ex1 = plus %% fanout identity identity

ex2 :: AD Double Double
ex2 = mult %% fanout identity identity

-------------------------------------------------------------------------------
-- Quadratic function
-------------------------------------------------------------------------------

quad :: AD (Double, Double) Double
quad = proc (x, y) -> do
    x2  <- mult    -< (x, x)
    y2  <- mult    -< (y, y)
    tmp <- plus    -< (x2, y2)
    z   <- konst 5 -< ()
    plus -< (tmp, z)

-------------------------------------------------------------------------------
-- Newton
-------------------------------------------------------------------------------

findZero :: AD Double Double -> Double -> [Double]
findZero f x0 = take 10 results
  where
    results = iterate go x0

    go :: Double -> Double
    go x =
        let (y, m) = evaluateAD f x
            [[y']] = LA.toLists m
        in x - gamma * (y / y')

    gamma = 0.1

-------------------------------------------------------------------------------
-- Gradient descent
-------------------------------------------------------------------------------

gradDesc :: forall a. VectorSpace a => AD a Double -> a -> [a]
gradDesc f = iterate go where
    go :: a -> a
    go x =
        let (_, m) = evaluateAD f x
            [grad] = LA.toLists $ LA.tr $ LA.scale gamma m

        in fromVector $ zipWith (-) (toVector x) grad

    gamma = 0.1

-------------------------------------------------------------------------------
-- Random
-------------------------------------------------------------------------------

randomDoubles :: Word64 -> [Double]
randomDoubles seed = go (SM.mkSMGen seed) where
    go g = let (d, g') = SM.nextDouble g in d : go g'

-------------------------------------------------------------------------------
-- Dot
-------------------------------------------------------------------------------

class VectorSpace' a where
    sumN :: AD a Double
    multN :: AD (a, a) a

instance (VectorSpace' a, VectorSpace' b) => VectorSpace' (a, b) where
    sumN = proc (x, y) -> do
        x' <- sumN -< x
        y' <- sumN -< y
        plus -< (x', y')

    multN = proc ((x1, x2), (y1, y2)) -> do
        z1 <- multN -< (x1, y1)
        z2 <- multN -< (x2, y2)
        identity -< (z1, z2)

instance VectorSpace' Double where
    sumN  = identity
    multN = mult

dot :: VectorSpace' a => AD (a, a) Double
dot = sumN %% multN

-------------------------------------------------------------------------------
-- ML stuff
-------------------------------------------------------------------------------

tanhAD :: AD Double Double
tanhAD = AD $ \x ->
    let y = tanh x
    in (y, linear (1 - y * y))

sigmoidAD :: AD Double Double
sigmoidAD = AD $ \x ->
    let y = 1 / (1 + exp (- x))
    in (x, linear (y * (1 - y)))

-- | weights for 2x1 connection. Two weights and bias.
type Weights' = ((Double, Double), Double)

-- | Two internal neurons, and final output
type Weights = ((Weights', Weights'), Weights')

startWeights :: Weights
startWeights = fromVector $ randomDoubles 1337

neuron :: AD (Weights', (Double, Double)) Double
neuron = proc ((ws, bias), i) -> do
    o <- dot -< (ws, i)
    tanhAD %% plus -< (o, bias)

network :: AD (Weights, (Double, Double)) Double
network = proc (((w1, w2), w3), xy) -> do
    u <- neuron  -< (w1, xy)
    v <- neuron  -< (w2, xy)
    neuron -< (w3, (u, v))

networkError :: AD Weights Double
networkError = proc ws -> do
    -- xor!
    s1 <- ex 1 1 0 -< ws
    s2 <- ex 0 0 0 -< ws
    s3 <- ex 1 0 1 -< ws
    s4 <- ex 0 1 1 -< ws

    sumN -< ((s1,s2), (s3, s4))
  where
    ex :: Double -> Double -> Double -> AD Weights Double
    ex x y z = proc ws -> do
         x1 <- konst x -< ()
         y1 <- konst y -< ()
         e1 <- konst z -< ()
         a1 <- network -< (ws, (x1, y1))
         r1 <- minus   -< (e1, a1)
         mult -< (r1, r1)

train :: Weights
train = gradDesc networkError startWeights !! 500

-------------------------------------------------------------------------------
-- Main
-------------------------------------------------------------------------------

main :: IO ()
main = do
    putStrLn $ "quad (2,3) = " ++ show (evaluateAD quad (2,3))
    putStrLn $ "gradDesc quad (2,3) = " ++ show (gradDesc quad (2,3) !! 30)

    print $ evaluateAD tanhAD 1
    print $ evaluateAD sigmoidAD 1

    putStrLn "Training the net (for xor)"
    let ws = train
    putStrLn $ "Parameters = " ++ show (toVector ws)
    putStrLn $ "Error = " ++ show (fst $ evaluateAD networkError ws)
    let example xy =
          putStrLn $ "eval " ++ show xy ++ " = " ++ showFFloat (Just 2) (fst $ evaluateAD network (ws, xy)) ""

    example (0, 0)
    example (0, 1)
    example (1, 0)
    example (1, 1)

    args <- getArgs
    when ("plot" `elem` args) $ do
        putStrLn "Outputting plot data: datafile.dat"

        let n = 20 :: Int
        let points = [ fromIntegral x / fromIntegral n | x <- [0..n] ] :: [Double]

        let output :: String
            output = unlines
                [ intercalate "\t"
                    [ show x
                    , show y
                    , show (fst (evaluateAD network (ws, (x, y))))
                    ]
                | x <- points
                , y <- points
                ]

        writeFile "datafile.dat" output

max_stars_repo_name: Mikolaj/overloaded

max_issues_repo_name: Mikolaj/overloaded

max_forks_repo_name: Mikolaj/overloaded


text:
-- Take a Ten Minute Walk
-- http://www.codewars.com/kata/54da539698b8a2ad76000228/

module Codewars.Kata.TenMinuteWalk where

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

isValidWalk :: String -> Bool
isValidWalk walk = case splitAt 10 walk of (ds, []) -> (length walk == 10) &&
                                                       ((==0) . magnitude . sum . map (\c -> case c of
                                                           'n' -> 0 :+ 1
                                                           's' -> 0 :+ (-1)
                                                           'w' -> (-1) :+ 0
                                                           'e' -> 1 :+ 0) $ walk)
                                           (_, rest) -> False

max_stars_repo_path: src/6 kyu/TenMinuteWalk.hs

max_stars_repo_name: gafiatulin/codewars

max_issues_repo_path: src/6 kyu/TenMinuteWalk.hs

max_issues_repo_name: gafiatulin/codewars

max_forks_repo_path: src/6 kyu/TenMinuteWalk.hs

max_forks_repo_name: gafiatulin/codewars


text:
import System.Environment
import Data.Complex
import Paraiso

main = do
  args <- getArgs
  let arch = if "--cuda" `elem` args then
               CUDA 128 128
             else
               X86
  putStrLn $ compile arch code
    where
      code = do
         parallel 16384 $ do
           c <- allocate
           z <- allocate
           c =$ (Rand (-2.0) 2.0) :+ (Rand (-2.0) 2.0)
           z =$ (0 :+ 0 :: Complex (Expr Double))
           cuda $ do
             sequential (65536) $ do
               z =$ z * z + c
           output [realPart c,imagPart c, realPart z, imagPart z]