text:
{-# LANGUAGE RecordWildCards, FlexibleContexts, GADTs, ExistentialQuantification #-}

{-| FIR filtering, decimation and resampling.

    FIR filters (and decimators, resamplers) work by taking successive dot products between the filter coefficients and the input data at increasing offsets. Sometimes the dot product fits entirely within one input buffer and other times it spans two input buffers (but never more because we assume that the filter length is less than the buffer size).

    We divide the filtering code by these two cases. Each filter (or decimator, resampler) is described by a data structure such as `Filter` with two functions, one for filtering within a single buffer and one that crosses buffers.

    The user must first create one of these data structures using the helper functions and pass this data structure to one of `firFilter`, `firDecimator`, or `firResampler` to create the `Pipe` that does the filtering. For example:

    > decimatorStruct   <- fastDecimatorC cpuInfo decimation coeffs
    > let decimatorPipe :: Pipe (Vector (Complex Float)) (Vector (Complex Float)) IO ()
    >     decimatorPipe =  firDecimator decimatorStruct outputSize

    There are polymorphic Haskell only implementations of filtering, decimation and resampling, for example, `haskellFilter`. In addition, there are optimised C implementations that use SIMD instructions on x86 machines, such as `fastFilterR`. These are always specialized to either real or complex numbers. There are also even faster implementations specialized for the case where the filter coefficients are symmetric as in a linear phase filter such as `fastFilterSymR`.

    The Haskell implementations are reasonably fast due to the Vector library and GHC's LLVM backend, however, if speed is important you are much better off with the C implementations.

    In the future we may avoid the cross buffer filtering function by mapping the buffers consecutively in memory as (I believe) GNU Radio does.

    An extensive benchmark suite exists in the /benchmarks subdirectory of this package.
-}
module SDR.Filter (
    -- * Types
    Filter(..),
    Decimator(..),
    Resampler(..),

    -- * Helper Functions
    -- ** Filters
    haskellFilter,

    -- *** Real Data
    fastFilterCR,
    fastFilterSSER,
    fastFilterAVXR,
    fastFilterR,

    -- *** Complex Data
    fastFilterCC,
    fastFilterSSEC,
    fastFilterAVXC,
    fastFilterC,

    -- *** Linear Phase Real Data
    fastFilterSymSSER,
    fastFilterSymAVXR,
    fastFilterSymR,

    -- ** Decimators
    haskellDecimator,

    -- *** Real Data
    fastDecimatorCR,
    fastDecimatorSSER,
    fastDecimatorAVXR,
    fastDecimatorR,

    -- *** Complex Data
    fastDecimatorCC,
    fastDecimatorSSEC,
    fastDecimatorAVXC,
    fastDecimatorC,

    -- *** Linear Phase Real Data
    fastDecimatorSymSSER,
    fastDecimatorSymAVXR,
    fastDecimatorSymR,

    -- ** Resamplers
    haskellResampler,

    -- *** Real Data
    fastResamplerCR,
    fastResamplerSSER,
    fastResamplerAVXR,
    fastResamplerR,

    -- *** Complex Data
    fastResamplerCC,
    fastResamplerSSEC,
    fastResamplerAVXC,
    fastResamplerC,

    -- * Filter
    firFilter,

    -- * Decimate
    firDecimator,

    -- * Resample
    firResampler,

    -- * DC Blocking Filter
    dcBlockingFilter
    ) where

import           Control.Applicative
import           Data.Complex
import           Control.Exception           hiding (assert)
import qualified Data.Vector.Generic         as VG
import qualified Data.Vector.Generic.Mutable as VGM
import qualified Data.Vector.Storable        as VS
import           Control.Monad.Primitive

import           Pipes

import           SDR.Util
import           SDR.FilterInternal          hiding (mkResampler, mkResamplerC)
import           SDR.CPUID

{- | A `Filter` contains all of the information needed by the `filterr`
     function to perform filtering. i.e. it contains the filter coefficients
     and pointers to the functions to do the actual filtering.
-}
data Filter m v vm a = Filter {
    numCoeffsF    :: Int,
    filterOne     :: Int -> v a -> vm (PrimState m) a -> m (),
    filterCross   :: Int -> v a -> v a -> vm (PrimState m) a -> m ()
}

{- | A `Decimator` contains all of the information needed by the `decimate`
     function to perform decimation i.e. it contains the filter coefficients
     and pointers to the functions to do the actual decimation.
-}
data Decimator m v vm a = Decimator {
    numCoeffsD    :: Int,
    decimationD   :: Int,
    decimateOne   :: Int -> v a -> vm (PrimState m) a -> m (),
    decimateCross :: Int -> v a -> v a -> vm (PrimState m) a -> m ()
}

{- | A `Resampler` contains all of the information needed by the `resample`
     function to perform resampling i.e. it contains the filter coefficients
     and pointers to the functions to do the actual resampling.
-}
data Resampler m v vm a = forall dat. Resampler {
    numCoeffsR     :: Int,
    decimationR    :: Int,
    interpolationR :: Int,
    startDat       :: dat,
    resampleOne    :: dat -> Int -> v a -> vm (PrimState m) a -> m (dat, Int),
    resampleCross  :: dat -> Int -> v a -> v a -> vm (PrimState m) a -> m (dat, Int)
}

duplicate :: [a] -> [a]
duplicate = concatMap func
    where func x = [x, x]

{-# INLINE haskellFilter #-}
-- | Returns a slow Filter data structure entirely implemented in Haskell
haskellFilter :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a)
              => [b]                  -- ^ The filter coefficients
              -> IO (Filter m v vm a) -- ^ The `Filter` data structure
haskellFilter coeffs = do
    let vCoeffs     = VG.fromList coeffs
    evaluate vCoeffs
    let filterOne   = filterHighLevel      vCoeffs
        filterCross = filterCrossHighLevel vCoeffs
        numCoeffsF  = length coeffs
    return Filter {..}

mkFilter :: Int
         -> FilterRR
         -> [Float]
         -> IO (Filter IO VS.Vector VS.MVector Float)
mkFilter sizeMultiple filterFunc coeffs = do
    let l          = length coeffs
        numCoeffsF = roundUp l sizeMultiple
        diff       = numCoeffsF - l
        vCoeffs    = VG.fromList $ coeffs ++ replicate diff 0
    evaluate vCoeffs
    let filterOne   = filterFunc           vCoeffs
        filterCross = filterCrossHighLevel vCoeffs
    return Filter {..}

-- | Returns a fast Filter data structure implemented in C. For filtering real data with real coefficients.
fastFilterCR :: [Float]                                   -- ^ The filter coefficients
             -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure
fastFilterCR = mkFilter 1 filterCRR

-- | Returns a fast Filter data structure implemented in C using SSE instructions. For filtering real data with real coefficients.
fastFilterSSER :: [Float]                                -- ^ The filter coefficients
            -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure
fastFilterSSER = mkFilter 4 filterCSSERR

-- | Returns a fast Filter data structure implemented in C using AVX instructions. For filtering real data with real coefficients.
fastFilterAVXR :: [Float]                                -- ^ The filter coefficients
            -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure
fastFilterAVXR = mkFilter 8 filterCAVXRR

-- | Returns a fast Filter data structure implemented in C using the fastest SIMD instruction set your processor supports. For filtering real data with real coefficients.
fastFilterR :: CPUInfo                                   -- ^ The CPU's capabilities
            -> [Float]                                   -- ^ The filter coefficients
            -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure
fastFilterR info = featureSelect info fastFilterCR [(hasAVX, fastFilterAVXR), (hasSSE42, fastFilterSSER)]

mkFilterC :: Int
          -> FilterRC
          -> [Float]
          -> IO (Filter IO VS.Vector VS.MVector (Complex Float))
mkFilterC sizeMultiple filterFunc coeffs = do
    let l           = length coeffs
        numCoeffsF  = roundUp sizeMultiple l
        diff        = numCoeffsF - l
        vCoeffs     = VG.fromList $ duplicate $ coeffs ++ replicate diff 0
        vCoeffs2    = VG.fromList $ coeffs ++ replicate diff 0
    evaluate vCoeffs
    let filterOne   = filterFunc           vCoeffs
        filterCross = filterCrossHighLevel vCoeffs2
    return Filter {..}

-- | Returns a fast Filter data structure implemented in C For filtering complex data with real coefficients.
fastFilterCC :: [Float]                                             -- ^ The filter coefficients
             -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Filter` data structure
fastFilterCC = mkFilterC 1 filterCRC

-- | Returns a fast Filter data structure implemented in C using SSE instructions. For filtering complex data with real coefficients.
fastFilterSSEC :: [Float]                                          -- ^ The filter coefficients
               -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Filter` data structure
fastFilterSSEC = mkFilterC 2 filterCSSERC

-- | Returns a fast Filter data structure implemented in C using AVX instructions. For filtering complex data with real coefficients.
fastFilterAVXC :: [Float]                                          -- ^ The filter coefficients
               -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Filter` data structure
fastFilterAVXC = mkFilterC 4 filterCAVXRC

-- | Returns a fast Filter data structure implemented in C using the fastest SIMD instruction set your processor supports. For filtering complex data with real coefficients.
fastFilterC :: CPUInfo                                             -- ^ The CPU's capabilities
            -> [Float]                                             -- ^ The filter coefficients
            -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Filter` data structure
fastFilterC info = featureSelect info fastFilterCC [(hasAVX, fastFilterAVXC), (hasSSE42, fastFilterSSEC)]

mkFilterSymR :: FilterRR
             -> [Float]
             -> IO (Filter IO VS.Vector VS.MVector Float)
mkFilterSymR filterFunc coeffs = do
    let vCoeffs     = VG.fromList coeffs
    let vCoeffs2    = VG.fromList $ coeffs ++ reverse coeffs
    evaluate vCoeffs
    evaluate vCoeffs2
    let filterOne   = filterFunc          vCoeffs
        filterCross = filterCrossHighLevel vCoeffs2
        numCoeffsF  = length coeffs * 2
    return Filter {..}

-- | Returns a fast Filter data structure implemented in C using SSE instructions. For filtering real data with real coefficients. For filters with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4.
fastFilterSymSSER :: [Float]                                   -- ^ The first half of the filter coefficients
                  -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure
fastFilterSymSSER = mkFilterSymR filterCSSESymmetricRR

-- | Returns a fast Filter data structure implemented in C using AVX instructions. For filtering real data with real coefficients. For filters with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4.
fastFilterSymAVXR :: [Float]                                   -- ^ The first half of the filter coefficients
                  -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure
fastFilterSymAVXR = mkFilterSymR filterCAVXSymmetricRR

-- | Returns a fast Filter data structure implemented in C using the fastest SIMD instruction set your processor supports. For filtering complex data with real coefficients. For filters with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4.
fastFilterSymR :: CPUInfo                                   -- ^ The CPU's capabilities
               -> [Float]                                   -- ^ The filter coefficients
               -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure
fastFilterSymR info = featureSelect info (error "At least SSE4.2 required") [(hasAVX, fastFilterSymAVXR), (hasSSE42, fastFilterSymSSER)]

{-# INLINE haskellDecimator #-}
-- | Returns a slow Decimator data structure entirely implemented in Haskell
haskellDecimator :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a)
                 => Int                     -- ^ The decimation factor
                 -> [b]                     -- ^ The filter coefficients
                 -> IO (Decimator m v vm a) -- ^ The `Decimator` data structure
haskellDecimator decimationD coeffs = do
    let vCoeffs     = VG.fromList coeffs
    evaluate vCoeffs
    let decimateOne   = decimateHighLevel      decimationD vCoeffs
        decimateCross = decimateCrossHighLevel decimationD vCoeffs
        numCoeffsD    = length coeffs
    return $ Decimator {..}

mkDecimator :: Int
            -> DecimateRR
            -> Int
            -> [Float]
            -> IO (Decimator IO VS.Vector VS.MVector Float)
mkDecimator sizeMultiple filterFunc decimationD coeffs = do
    let l          = length coeffs
        numCoeffsD = roundUp l sizeMultiple
        diff       = numCoeffsD - l
        vCoeffs    = VG.fromList $ coeffs ++ replicate diff 0
    evaluate vCoeffs
    let decimateOne   = filterFunc             decimationD vCoeffs
        decimateCross = decimateCrossHighLevel decimationD vCoeffs
    return Decimator {..}

-- | Returns a fast Decimator data structure implemented in C. For decimating real data with real coefficients.
fastDecimatorCR :: Int                                          -- ^ The decimation factor
                -> [Float]                                      -- ^ The filter coefficients
                -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure
fastDecimatorCR = mkDecimator 1 decimateCRR

-- | Returns a fast Decimator data structure implemented in C using SSE instructions. For decimating real data with real coefficients.
fastDecimatorSSER :: Int                                          -- ^ The decimation factor
                  -> [Float]                                      -- ^ The filter coefficients
                  -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure
fastDecimatorSSER = mkDecimator 4 decimateCSSERR

-- | Returns a fast Decimator data structure implemented in C using AVX instructions. For decimating real data with real coefficients.
fastDecimatorAVXR :: Int                                          -- ^ The decimation factor
                  -> [Float]                                      -- ^ The filter coefficients
                  -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure
fastDecimatorAVXR = mkDecimator 8 decimateCAVXRR

-- | Returns a fast Decimator data structure implemented in C using the fastest SIMD instruction set your processor supports. For decimating real data with real coefficients.
fastDecimatorR :: CPUInfo                                      -- ^ The CPU's capabilities
               -> Int                                          -- ^ The decimation factor
               -> [Float]                                      -- ^ The filter coefficients
               -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure
fastDecimatorR info = featureSelect info fastDecimatorCR [(hasAVX, fastDecimatorAVXR), (hasSSE42, fastDecimatorSSER)]

mkDecimatorC :: Int
             -> DecimateRC
             -> Int
             -> [Float]
             -> IO (Decimator IO VS.Vector VS.MVector (Complex Float))
mkDecimatorC sizeMultiple filterFunc decimationD coeffs = do
    let l          = length coeffs
        numCoeffsD = roundUp l sizeMultiple
        diff       = numCoeffsD - l
        vCoeffs    = VG.fromList $ duplicate $ coeffs ++ replicate diff 0
        vCoeffs2   = VG.fromList $ coeffs ++ replicate diff 0
    evaluate vCoeffs
    let decimateOne   = filterFunc             decimationD vCoeffs
        decimateCross = decimateCrossHighLevel decimationD vCoeffs2
    return Decimator {..}

-- | Returns a fast Decimator data structure implemented in C. For decimating complex data with real coefficients.
fastDecimatorCC :: Int                                                    -- ^ The decimation factor
                -> [Float]                                                -- ^ The filter coefficients
                -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Decimator` data structure
fastDecimatorCC = mkDecimatorC 1 decimateCRC

-- | Returns a fast Decimator data structure implemented in C using SSE instructions. For decimating complex data with real coefficients.
fastDecimatorSSEC :: Int                                                    -- ^ The decimation factor
                  -> [Float]                                                -- ^ The filter coefficients
                  -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Decimator` data structure
fastDecimatorSSEC = mkDecimatorC 2 decimateCSSERC

-- | Returns a fast Decimator data structure implemented in C using AVX instructions. For decimating complex data with real coefficients.
fastDecimatorAVXC :: Int                                                 -- ^ The decimation factor
               -> [Float]                                                -- ^ The filter coefficients
               -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Decimator` data structure
fastDecimatorAVXC = mkDecimatorC 4 decimateCAVXRC

-- | Returns a fast Decimator data structure implemented in C using the fastest SIMD instruction set your processor supports. For decimating complex data with real coefficients.
fastDecimatorC :: CPUInfo                                                -- ^ The CPU's capabilities
               -> Int                                                    -- ^ The decimation factor
               -> [Float]                                                -- ^ The filter coefficients
               -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Decimator` data structure
fastDecimatorC info = featureSelect info fastDecimatorCC [(hasAVX, fastDecimatorAVXC), (hasSSE42, fastDecimatorSSEC)]

mkDecimatorSymR :: DecimateRR
                -> Int
                -> [Float]
                -> IO (Decimator IO VS.Vector VS.MVector Float)
mkDecimatorSymR filterFunc decimationD coeffs = do
    let vCoeffs    = VG.fromList coeffs
    let vCoeffs2   = VG.fromList $ coeffs ++ reverse coeffs
    evaluate vCoeffs
    evaluate vCoeffs2
    let decimateOne   = filterFunc             decimationD vCoeffs
        decimateCross = decimateCrossHighLevel decimationD vCoeffs2
        numCoeffsD    = length coeffs * 2

    return Decimator {..}
-- | Returns a fast Decimator data structure implemented in C using SSE instructions. For decimating real data with real coefficients. For decimators with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4.
fastDecimatorSymSSER :: Int                                          -- ^ The decimation factor
                     -> [Float]                                      -- ^ The first half of the filter coefficients
                     -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure
fastDecimatorSymSSER = mkDecimatorSymR decimateCSSESymmetricRR

-- | Returns a fast Decimator data structure implemented in C using AVX instructions. For decimating real data with real coefficients. For decimators with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4.
fastDecimatorSymAVXR :: Int                                          -- ^ The decimation factor
                     -> [Float]                                      -- ^ The first half of the filter coefficients
                     -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure
fastDecimatorSymAVXR = mkDecimatorSymR decimateCAVXSymmetricRR

-- | Returns a fast Decimator data structure implemented in C using the fastest SIMD instruction set your processor supports. For decimating real data with real coefficients. For decimators with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4.
fastDecimatorSymR :: CPUInfo                                      -- ^ The CPU's capabilities
                  -> Int                                          -- ^ The decimation factor
                  -> [Float]                                      -- ^ The filter coefficients
                  -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure
fastDecimatorSymR info = featureSelect info (error "at least AVX required") [(hasAVX, fastDecimatorSymAVXR), (hasSSE42, fastDecimatorSymSSER)]

{-# INLINE haskellResampler #-}
-- | Returns a slow Resampler data structure entirely implemented in Haskell
haskellResampler :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a)
                  => Int                     -- ^ The interpolation factor
                  -> Int                     -- ^ The decimation factor
                  -> [b]                     -- ^ The filter coefficients
                  -> IO (Resampler m v vm a) -- ^ The `Resampler` data structure
haskellResampler interpolationR decimationR coeffs = do
    let vCoeffs     = VG.fromList coeffs
    evaluate vCoeffs
    let resampleOne   v w x y   = func <$> resampleHighLevel      interpolationR decimationR vCoeffs v w x y
        resampleCross v w x y z = func <$> resampleCrossHighLevel interpolationR decimationR vCoeffs v w x y z
        numCoeffsR            = length coeffs
        func          x       = (x, x)
        startDat              = 0
    return Resampler {..}

mkResampler :: Int
            -> ResampleRR
            -> Int
            -> Int
            -> [Float]
            -> IO (Resampler IO VS.Vector VS.MVector Float)
mkResampler sizeMultiple filterFunc interpolationR decimationR coeffs = do
    let vCoeffs     = VG.fromList coeffs
    evaluate vCoeffs
    resamp <- filterFunc interpolationR decimationR coeffs
    let resampleOne   v w x y   = func1 <$> resamp (fst v) w x y
        resampleCross (group, offset) count x y z = do
            offset' <- resampleCrossHighLevel interpolationR decimationR vCoeffs offset count x y z
            return (((group + count) `mod` interpolationR, offset'), offset')
        numCoeffsR              = roundUp (length coeffs) (interpolationR * sizeMultiple)
        func1 group             = let offset = interpolationR - 1 - ((interpolationR + group * decimationR - 1) `mod` interpolationR) in ((group, offset), offset)
        startDat                = (0, 0)
    return Resampler {..}

mkResamplerC :: Int
             -> ResampleRC
             -> Int
             -> Int
             -> [Float]
             -> IO (Resampler IO VS.Vector VS.MVector (Complex Float))
mkResamplerC sizeMultiple filterFunc interpolationR decimationR coeffs = do
    let vCoeffs     = VG.fromList coeffs
    evaluate vCoeffs
    resamp <- filterFunc interpolationR decimationR coeffs
    let resampleOne   v w x y   = func1 <$> resamp (fst v) w x y
        resampleCross (group, offset) count x y z = do
            offset' <- resampleCrossHighLevel interpolationR decimationR vCoeffs offset count x y z
            return (((group + count) `mod` interpolationR, offset'), offset')
        numCoeffsR              = roundUp (length coeffs) (interpolationR * sizeMultiple)
        func1 group             = let offset = interpolationR - 1 - ((interpolationR + group * decimationR - 1) `mod` interpolationR) in ((group, offset), offset)
        startDat                = (0, 0)
    return Resampler {..}

-- | Returns a fast Resampler data structure implemented in C. For filtering real data with real coefficients.
fastResamplerCR :: Int                                          -- ^ The interpolation factor
                -> Int                                          -- ^ The decimation factor
                -> [Float]                                      -- ^ The filter coefficients
                -> IO (Resampler IO VS.Vector VS.MVector Float) -- ^ The `Resampler` data structure
fastResamplerCR = mkResampler 1 resampleCRR2

-- | Returns a fast Resampler data structure implemented in C using SSE instructions. For filtering real data with real coefficients.
fastResamplerSSER :: Int                                          -- ^ The interpolation factor
                  -> Int                                          -- ^ The decimation factor
                  -> [Float]                                      -- ^ The filter coefficients
                  -> IO (Resampler IO VS.Vector VS.MVector Float) -- ^ The `Resampler` data structure
fastResamplerSSER = mkResampler 4 resampleCSSERR

-- | Returns a fast Resampler data structure implemented in C using AVX instructions. For filtering real data with real coefficients.
fastResamplerAVXR :: Int                                          -- ^ The interpolation factor
                  -> Int                                          -- ^ The decimation factor
                  -> [Float]                                      -- ^ The filter coefficients
                  -> IO (Resampler IO VS.Vector VS.MVector Float) -- ^ The `Resampler` data structure
fastResamplerAVXR = mkResampler 8 resampleCAVXRR

-- | Returns a fast Resampler data structure implemented in C using the fastest SIMD instruction set your processor supports. For resampling real data with real coefficients.
fastResamplerR :: CPUInfo                                      -- ^ The CPU's capabilities
               -> Int                                          -- ^ The interpolation factor
               -> Int                                          -- ^ The decimation factor
               -> [Float]                                      -- ^ The filter coefficients
               -> IO (Resampler IO VS.Vector VS.MVector Float) -- ^ The `Resampler` data structure
fastResamplerR info = featureSelect info fastResamplerCR [(hasAVX, fastResamplerAVXR), (hasSSE42, fastResamplerSSER)]

-- | Returns a fast Resampler data structure implemented in C. For filtering complex data with real coefficients.
fastResamplerCC :: Int                                          -- ^ The interpolation factor
                -> Int                                          -- ^ The decimation factor
                -> [Float]                                      -- ^ The filter coefficients
                -> IO (Resampler IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Resampler` data structure
fastResamplerCC = mkResamplerC 1 resampleCRC

-- | Returns a fast Resampler data structure implemented in C using SSE instructions. For filtering complex data with real coefficients.
fastResamplerSSEC :: Int                                          -- ^ The interpolation factor
                  -> Int                                          -- ^ The decimation factor
                  -> [Float]                                      -- ^ The filter coefficients
                  -> IO (Resampler IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Resampler` data structure
fastResamplerSSEC = mkResamplerC 4 resampleCSSERC

-- | Returns a fast Resampler data structure implemented in C using AVX instructions. For filtering complex data with real coefficients.
fastResamplerAVXC :: Int                                          -- ^ The interpolation factor
                  -> Int                                          -- ^ The decimation factor
                  -> [Float]                                      -- ^ The filter coefficients
                  -> IO (Resampler IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Resampler` data structure
fastResamplerAVXC = mkResamplerC 8 resampleCAVXRC

-- | Returns a fast Resampler data structure implemented in C using the fastest SIMD instruction set your processor supports. For resampling complex data with real coefficients.
fastResamplerC :: CPUInfo                                      -- ^ The CPU's capabilities
               -> Int                                          -- ^ The interpolation factor
               -> Int                                          -- ^ The decimation factor
               -> [Float]                                      -- ^ The filter coefficients
               -> IO (Resampler IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Resampler` data structure
fastResamplerC info = featureSelect info fastResamplerCC [(hasAVX, fastResamplerAVXC), (hasSSE42, fastResamplerSSEC)]

data Buffer v a = Buffer {
    buffer :: v a,
    offset :: Int
}

space Buffer{..} = VGM.length buffer - offset

newBuffer :: (PrimMonad m, VGM.MVector vm a) => Int -> m (Buffer (vm (PrimState m)) a)
newBuffer size = do
    buf <- VGM.new size
    return $ Buffer buf 0

advanceOutBuf :: (PrimMonad m, VG.Vector v a) => Int -> Buffer (VG.Mutable v (PrimState m)) a -> Int -> Pipe b (v a) m (Buffer (VG.Mutable v (PrimState m)) a)
advanceOutBuf blockSizeOut buf@(Buffer bufOut offsetOut) count =
    if count == space buf then do
        bufOutF <- lift $ VG.unsafeFreeze bufOut
        yield bufOutF
        lift $ newBuffer blockSizeOut
    else
        return $ Buffer bufOut (offsetOut + count)

-- | My own assert implementation since the GHC one doesnt seem to work even with optimisations disabled and using -fno-ignore-asserts
assert loc False = error loc
assert loc True  = return ()

--Filtering
{-# INLINE firFilter #-}
{-| Create a pipe that performs filtering -}
firFilter :: (PrimMonad m, Functor m, VG.Vector v a, Num a)
        => Filter m v (VG.Mutable v) a -- ^ The `Filter` data structure
        -> Int                         -- ^ The output block size
        -> Pipe (v a) (v a) m ()       -- ^ The `Pipe` that does the filtering
firFilter Filter{..} blockSizeOut = do
    inBuf  <- await
    outBuf <- lift $ newBuffer blockSizeOut
    simple inBuf outBuf

    where

    simple bufIn bufferOut@(Buffer bufOut offsetOut) = do
        assert "filter 1" (VG.length bufIn >= numCoeffsF)

        let count = min (VG.length bufIn - numCoeffsF + 1) (space bufferOut)
        lift $ filterOne count bufIn (VGM.unsafeDrop offsetOut bufOut)

        bufferOut' <- advanceOutBuf blockSizeOut bufferOut count
        let bufIn' =  VG.drop count bufIn

        case VG.length bufIn' < numCoeffsF of
            False -> simple bufIn' bufferOut'
            True  -> do
                next <- await
                crossover bufIn' next bufferOut'

    crossover bufLast bufNext bufferOut@(Buffer bufOut offsetOut) = do
        assert "filter 2" (VG.length bufLast < numCoeffsF)
        assert "filter 3" (VG.length bufLast > 0)

        let count = min (VG.length bufLast) (space bufferOut)
        lift $ filterCross count bufLast bufNext (VGM.unsafeDrop offsetOut bufOut)

        bufferOut' <- advanceOutBuf blockSizeOut bufferOut count

        case VG.length bufLast == count of
            True  -> simple bufNext bufferOut'
            False -> crossover (VG.drop count bufLast) bufNext bufferOut'

--Decimation
{-# INLINE firDecimator #-}
{-| Create a pipe that performs decimation -}
firDecimator :: (PrimMonad m, Functor m, VG.Vector v a, Num a)
         => Decimator m v (VG.Mutable v) a -- ^ The `Decimator` data structure
         -> Int                            -- ^ The output block size
         -> Pipe (v a) (v a) m ()          -- ^ The `Pipe` that does the decimation
firDecimator Decimator{..} blockSizeOut = do
    inBuf  <- await
    outBuf <- lift $ newBuffer blockSizeOut
    simple inBuf outBuf

    where

    simple bufIn bufferOut@(Buffer bufOut offsetOut) = do
        assert "decimate 1" (VG.length bufIn >= numCoeffsD)

        let count = min (((VG.length bufIn - numCoeffsD) `quot` decimationD) + 1) (space bufferOut)
        lift $ decimateOne count bufIn (VGM.unsafeDrop offsetOut bufOut)

        bufferOut' <- advanceOutBuf blockSizeOut bufferOut count
        let bufIn' = VG.drop (count * decimationD) bufIn

        case VG.length bufIn' < numCoeffsD of
            False -> simple bufIn' bufferOut'
            True  -> do
                next <- await
                crossover bufIn' next bufferOut'

    crossover bufLast bufNext bufferOut@(Buffer bufOut offsetOut) = do
        assert "decimate 2" (VG.length bufLast < numCoeffsD)
        assert "decimate 3" (VG.length bufLast > 0)

        let count = min (VG.length bufLast `quotUp` decimationD) (space bufferOut)
        lift $ decimateCross count bufLast bufNext (VGM.unsafeDrop offsetOut bufOut)

        bufferOut' <- advanceOutBuf blockSizeOut bufferOut count

        case VG.length bufLast <= count * decimationD of
            True  -> simple (VG.drop (count * decimationD - VG.length bufLast) bufNext) bufferOut'
            False -> crossover (VG.drop (count * decimationD) bufLast) bufNext bufferOut'

{-
Rational Downsampling:

Input upsampled by 3:    |**|**|**|**|**|**|**|**|**|**|**|
Output downsampled by 7: |******|******|******|******|*****

                  Consider here ^
                  Next output is here  ^

Filter offset is 2

k is number of used inputs

filterOffset + k*interpolation = decimation + filterOffset'
where
    k > 0
    0 <= filterOffset, filterOffset' < interpolation

k*interpolation - filterOffset' = decimation - filterOffset
k*interpolation - filterOffset' - 1 = decimation - filterOffset - 1

(k-1) * interpolation + (interpolation - filterOffset' - 1) = decimation - filterOffset - 1

k             = (decimation - filterOffset - 1) / interpolation + 1
filterOffset' = interpolation - 1 - (decimation - filterOffset - 1) % interpolation

Only works if decimation > interpolation

-}

{-
Rational Upsampling:

Input upsampled by 7:    |******|******|******|******|*****
Output downsampled by 3: |**|**|**|**|**|**|**|**|**|**|**|

              Consider Here ^
              Next sample is   ^

Filter offset is 4

filterOffset + k * interpolation = decimation + filterOffset'
where
    k = {0, 1}
    0 <= filterOffset, filterOffset' < interpolation

k * interpolation + (interpolation - filterOffset' - 1) = decimation - filterOffset + interpolation - 1

k = (decimation - filterOffset + interpolation - 1) / interpolation

============================

Or, equivalently,

k = 0 | filterOffset >= decimation
    1 | otherwise

o = o - decimation + k * interpolation

-}

--Rational resampling
quotUp q d = (q + (d - 1)) `quot` d

{-# INLINE firResampler #-}
{-| Create a pipe that performs resampling -}
firResampler :: (PrimMonad m, VG.Vector v a, Num a)
         => Resampler m v (VG.Mutable v) a -- ^ The `Resampler` data structure
         -> Int                            -- ^ The output block size
         -> Pipe (v a) (v a) m ()          -- ^ The `Pipe` that does the resampling
firResampler Resampler{..} blockSizeOut = do
    inBuf  <- await
    outBuf <- lift $ newBuffer blockSizeOut
    simple inBuf outBuf startDat 0

    where

    simple bufIn bufferOut@(Buffer bufOut offsetOut) dat filterOffset = do
        assert "resample 1" (VG.length bufIn * interpolationR >= numCoeffsR - filterOffset)
        --available number of samples == interpolation * num_input
        --required number of samples  == decimation * (num_output - 1) + filter_length - filter_offset
        let count = min (((VG.length bufIn * interpolationR - numCoeffsR + filterOffset) `quot` decimationR) + 1) (space bufferOut)
        (dat, endOffset) <- lift $ resampleOne dat count bufIn (VGM.unsafeDrop offsetOut bufOut)
        assert "resample 2" ((count * decimationR + endOffset - filterOffset) `rem` interpolationR == 0)
        bufferOut' <- advanceOutBuf blockSizeOut bufferOut count
        --samples no longer needed starting from filterOffset == count * decimation - filterOffset
        --inputs lying in this region                         == (count * decimation - filterOffset) / interpolation (rounding up)
        let usedInput = (count * decimationR - filterOffset) `quotUp` interpolationR
            bufIn'    = VG.drop usedInput bufIn

        case VG.length bufIn' * interpolationR < numCoeffsR - endOffset of
            False -> simple bufIn' bufferOut' dat endOffset
            True  -> do
                next <- await
                --TODO: why is this not needed in filter and decimator
                case VG.length bufIn' == 0 of
                    True ->  simple    next bufferOut' dat endOffset
                    False -> crossover bufIn' next bufferOut' dat endOffset

    crossover bufLast bufNext bufferOut@(Buffer bufOut offsetOut) dat filterOffset = do
        assert "resample 3" (VG.length bufLast * interpolationR < numCoeffsR - filterOffset)
        --outputsComputable is the number of outputs that need to be computed for the last buffer to no longer be needed
        --outputsComputable * decimation == numInput * interpolation + filterOffset + k
        let outputsComputable = (VG.length bufLast * interpolationR + filterOffset) `quotUp` decimationR
            count = min outputsComputable (space bufferOut)
        assert "resample 4" (count /= 0)
        (dat, endOffset) <- lift $ resampleCross dat count bufLast bufNext (VGM.unsafeDrop offsetOut bufOut)
        assert "resample 5" ((count * decimationR + endOffset - filterOffset) `rem` interpolationR == 0)
        bufferOut' <- advanceOutBuf blockSizeOut bufferOut count

        let inputUsed = (count * decimationR - filterOffset) `quotUp` interpolationR

        case inputUsed >= VG.length bufLast of
            True  -> simple (VG.drop (inputUsed - VG.length bufLast) bufNext) bufferOut' dat endOffset
            False -> crossover (VG.drop inputUsed bufLast) bufNext bufferOut' dat endOffset

-- | A DC blocking filter
dcBlockingFilter :: Pipe (VS.Vector Float) (VS.Vector Float) IO ()
dcBlockingFilter = func 0 0
    where
    func lastSample lastOutput = do
        dat <- await
        out <- lift $ VGM.new (VG.length dat)
        (lastSample, lastOutput) <- lift $ dcBlocker (VG.length dat) lastSample lastOutput dat out
        outF <- lift $ VG.unsafeFreeze out
        yield outF
        func lastSample lastOutput

max_stars_repo_name: adamwalker/sdr

max_issues_repo_name: peixian/sdr

max_forks_repo_name: adamwalker/sdr


text:
{-
  this file is taken from the hommage library
-}
{-# LANGUAGE FlexibleContexts #-}
module WavFile
 (
 -- * Binary Files
   writeDataFile
 , readDataFile
 , openDataFile

 -- * WAV-Files
 , writeWavFile
 , writeWavFileMono
 , writeWavFileStereo
 , readWavFile
 , openWavFile

 , writeWavFiles

 -- * Cast between Int16 and Double representation of WAV-data
 , wavInt16ToDouble
 , wavDoubleToInt16

 -- * Low-Level implemetation
 -- ** Arrays and Files
 , readArrayFromFile
 , writeArrayToFileWithHeader
 , writeArrayToFile

 -- ** Single Stream
 , openSingleInputFile
 , openSingleOutputFile

 , openSingleInputWavFile
 , openSingleOutputWavFileMono
 , openSingleOutputWavFileStereo

 -- ** Buffered Stream
 , openInputFile
 , openOutputFile
 , openOutputFileWithHeader
 , openInputWavFile
 , openOutputWavFileMono
 , openOutputWavFileStereo

 -- ** Header Stuff and others
 , HeaderFun
 , HeaderSize
 , noHeader

 , wavHeaderFunMono
 , wavHeaderFunStereo
 , wavHeaderSize

 , initWriteWavHeaderMono
 , initWriteWavHeaderStereo
 , initReadWavHeader
 , closeWriteWavHeader
 , encode
 , decode
 , encodeWavLengt
 , initWavHeaderMono
 , initWavHeaderStereo
 , sizeOfArrayElements
 , inferSizeOfArrayElements
 , inferSizeOfArrayElements'
 )
 where

import System.IO
import Control.Monad
import Control.Monad.Fix
import Data.Char
import Data.Int
import GHC.IO.Handle
import GHC.IO
import GHC.Base hiding (mapM)
import Data.Complex

import Foreign.Storable
import Data.Array.Storable
import Data.Array.IArray
import Data.IORef

import System.IO.Unsafe
import GHC.Weak

---------------------------------------------------------------------------------------------------
wavInt16ToDouble :: Int16 -> Double
wavInt16ToDouble i = (fromIntegral i) / 32768.0

wavDoubleToInt16 :: Double -> Int16
wavDoubleToInt16 d = round (d * 32767.0)
---------------------------------------------------------------------------------------------------
readDataFile :: Storable a => FilePath -> IO [a]
readDataFile fp = init
 where
  init = do (c,r,_) <- openSingleInputFile 10000 0 fp
            mkWeak r () (Just c)
            loop r
  loop r = r >>= maybe (return []) (\a -> return (a : unsafePerformIO (loop r)))

openDataFile :: Storable a => FilePath -> [a]
openDataFile fp = unsafePerformIO $ readDataFile fp

writeDataFile :: Storable a => FilePath -> [a] -> IO ()
writeDataFile fp xs = do
 (c,w) <- openSingleOutputFile 4000 fp
 let loop (x:r) = w x >> loop r
     loop []    = c
 loop xs

---------------------------------------------------------------------------------------------------
readWavFile :: FilePath -> IO (Either [Int16] [(Int16,Int16)])
readWavFile fp = {- putStrLn ("opening " ++ fp) >> -} init
 where
  init = do (c,r,_) <- openSingleInputWavFile 100000 fp
            either (\rm -> do mkWeak rm () (Just c) --(putStrLn ("close " ++ fp) >> c))
                              fmap Left $ loop rm)
                   (\rs -> do mkWeak rs () (Just c) --(putStrLn ("close " ++ fp) >> c))
                              fmap Right $ loop rs)
                   r
  loop r = r >>= maybe (return []) (\a -> return (a : unsafePerformIO (loop r)))

openWavFile :: FilePath -> Either [Int16] [(Int16,Int16)]
openWavFile fp = unsafePerformIO $ readWavFile fp

writeWavFile :: FilePath -> Either [Int16] [(Int16, Int16)] -> IO ()
writeWavFile fp (Left w)  = writeWavFileMono fp w
writeWavFile fp (Right w) = writeWavFileStereo fp w

writeWavFileMono :: FilePath -> [Int16] -> IO ()
writeWavFileMono fp xs = do
 (c,w) <- openSingleOutputWavFileMono 100000 fp
 let loop (x:r) = w x >> loop r
     loop []    = c
 loop xs

writeWavFileStereo :: FilePath -> [(Int16,Int16)] -> IO ()
writeWavFileStereo fp xs = do
 (c,w) <- openSingleOutputWavFileStereo 100000 fp
 let loop (x:r) = w x >> loop r
     loop []    = c
 loop xs

writeWavFiles :: FilePath -> String  -> [Either [Int16] [(Int16, Int16)]] -> IO ()
writeWavFiles fp fn sns = do
  let init count (Left x : xs) = do (c,w) <- openSingleOutputWavFileMono 100000 (fp ++ fn ++ "_" ++ show count ++ ".wav")
                                    r <- init (count+1) xs
                                    return (Left (c,w,x) : r)
      init count (Right x : xs) = do (c,w) <- openSingleOutputWavFileStereo 100000 (fp ++ fn ++ "_" ++ show count ++ ".wav")
                                     r <- init (count+1) xs
                                     return (Right (c,w,x) : r)
      init _ _ = return []
  ss <- init 0 sns
  let step (Left (c,w,x:xs) : r) = do w x
                                      r <- step r
                                      return (Left (c,w,xs) : r)
      step (Left (c,w,[]) : r)   = c >> step r
      step (Right (c,w,x:xs) : r) = do w x
                                       r <- step r
                                       return (Right (c,w,xs) : r)
      step (Right (c,w,[]) : r)   = c >> step r
      step []                     = return []

      loop [] = return ()
      loop xs = step xs >>= loop
  loop ss
---------------------------------------------------------------------------------------------------
-- Header

-- | The first Action is applied after opening the file. Then the data
--   bytes are written. The second Action is then called with the number of
--   bytes of the data. It must close the handle.
type HeaderFun = (Handle -> IO (), Handle -> Int -> IO ())
type HeaderSize = Int

noHeader :: HeaderFun
noHeader = (const $ return (), \h _ -> return ())

-------------------------------------------------
-- WavHeader

-- | The first Action writes the Header of the wav-file.
--   Then the wav-data is written. The second action moves
--   the Handle to the positions where the lenght of the wav-file
--   are encoded (in the header) and writes the right number, which
--   is unknown before all data is written.
--   Afterwards it closes the file.
wavHeaderFunMono :: HeaderFun
wavHeaderFunMono = (initWriteWavHeaderMono, closeWriteWavHeader)

wavHeaderFunStereo :: HeaderFun
wavHeaderFunStereo = (initWriteWavHeaderStereo, closeWriteWavHeader)

wavHeaderSize :: HeaderSize
wavHeaderSize = 44
-------------------------------------------------
initWriteWavHeaderMono :: Handle -> IO ()
initWriteWavHeaderMono h = mapM_ (hPutChar h) initWavHeaderMono

initWriteWavHeaderStereo :: Handle -> IO ()
initWriteWavHeaderStereo h = mapM_ (hPutChar h) initWavHeaderStereo

initReadWavHeader :: Handle -> IO (Maybe (Bool, Int, Int, Int))
initReadWavHeader handle =
 let check [] = return True
     check (x:xs) = hGetChar handle >>= \y -> if x /= ord y then return False else check xs
     m >>? n = m >>= \b -> if b then n else return Nothing
 in hSeek handle AbsoluteSeek 0 >>
    check [82,73,70,70] >>?
    (replicateM 4 (hGetChar handle) >>= \l1 ->
    check [87, 65, 86, 69, 102,109,116,32, 16, 0,  0,  0, 1,  0] >>?
    (hGetChar handle >>= \s1 ->
    check [0, 68, 172, 0, 0] >>?
    (replicateM 4 (hGetChar handle) >>= \x1 ->
    (hGetChar handle >>= \s2 ->
    check [0,  16, 0, 100, 97,116,97] >>?
    (replicateM 4 (hGetChar handle) >>= \l2 ->
    let b | ord s1 == 1 && ord s2 == 2 = False
          | ord s1 == 2 && ord s2 == 4 = True
          | otherwise          = error ("invalid wav file " ++ show s1 ++ ", " ++ show s2)
    in return $ Just (b, decode l1, decode l2, decode x1) ) ) ) ) )

closeWriteWavHeader :: Handle -> Int -> IO ()
closeWriteWavHeader h len = do
 let (sa,sb) = encodeWavLengt len
 hSeek h AbsoluteSeek 4
 mapM (hPutChar h) sa
 hSeek h AbsoluteSeek 40
 mapM (hPutChar h) sb
 hClose h

encode :: Int -> String
encode a = map chr
 [ mod a 256
 , mod (div a 256) 256
 , mod (div a 65536) 256
 , mod (div a 16777216) 256
 ]

decode :: String -> Int
decode [a0,a1,a2,a3] = sum
 [ord a0, 256 * ord a1, 65536 * ord a2, 16777216 * ord a3]

encodeWavLengt :: Int -> (String, String)
encodeWavLengt len = (encode (len + 8), encode len)

initWavHeaderMono :: [Char]
initWavHeaderMono =
 map chr
 [82,73,70,70,
  0,0,0,0,
  87, 65, 86, 69,
  102,109,116,32,
  16, 0,  0,  0,
  1,  0,  1,  0,
  68, 172,0,  0,
  136,88, 1,  0,
  2,  0,  16, 0,
  100, 97,116,97,
  0,0,0,0]

initWavHeaderStereo :: [Char]
initWavHeaderStereo =
 map chr
 [82,73,70,70,
  0,0,0,0,
  87, 65, 86, 69,
  102,109,116,32,
  16, 0,  0,  0,
  1,  0,  2,  0,
  68, 172,0,  0,
  16,177, 2,  0,
  4,  0,  16, 0,
  100, 97,116,97,
  0,0,0,0]

---------------------------------------------------------------------------------------------------
sizeOfArrayElements :: Storable a => StorableArray Int a -> Int
sizeOfArrayElements x = sizeOf (foo x)
 where
  foo :: StorableArray Int a -> a
  foo x = undefined

inferSizeOfArrayElements :: Storable a => (StorableArray Int a -> IO Int) -> StorableArray Int a -> Int
inferSizeOfArrayElements _ arr = sizeOfArrayElements arr

inferSizeOfArrayElements' :: Storable a => (StorableArray Int a -> Int -> IO ()) -> StorableArray Int a -> Int
inferSizeOfArrayElements' _ arr = sizeOfArrayElements arr

inferSizeOfArrayElements'' :: Storable a => (StorableArray Int a -> Int -> IO Int) -> StorableArray Int a -> Int
inferSizeOfArrayElements'' _ arr = sizeOfArrayElements arr

---------------------------------------------------------------------------------------------------
-- Array / File
---------------------------------------------------------------------------------------------------
readArrayFromFile :: Storable a => HeaderSize -> FilePath -> IO (StorableArray Int a, Int)
readArrayFromFile offset filepath = do
 handle <- openBinaryFile filepath ReadMode
 hSeek handle AbsoluteSeek $ fromIntegral offset
 bytes <- hFileSize handle
 let bytesize = fromIntegral bytes - offset
 arr <- mfix $ \arr -> do let bufsize = div bytesize (sizeOfArrayElements arr)
                          newArray_ (0, bufsize-1)
 withStorableArray arr $ \ptr -> hGetBuf handle ptr bytesize
 hClose handle
 return (arr, div bytesize (sizeOfArrayElements arr))

writeArrayToFileWithHeader :: (MArray StorableArray a IO, Storable a) => HeaderFun -> FilePath -> StorableArray Int a -> IO ()
writeArrayToFileWithHeader (initH, closeH) filepath arr = do
 handle <- openBinaryFile filepath WriteMode
 initH handle
 (0, n) <- getBounds arr
 let elemsize = sizeOfArrayElements arr
     bytesize = (1+n) * elemsize
 withStorableArray arr (\ptr -> hPutBuf handle ptr bytesize)
 closeH handle bytesize

writeArrayToFile :: (MArray StorableArray a IO, Storable a) => FilePath -> StorableArray Int a -> IO ()
writeArrayToFile filepath arr = do
 handle <- openBinaryFile filepath WriteMode
 (0, n) <- getBounds arr
 let elemsize = sizeOfArrayElements arr
     bytesize = (1+n) * elemsize
 withStorableArray arr (\ptr -> hPutBuf handle ptr bytesize)
 hClose handle

{-
writeArrayToFileWithHeader :: Storable a => HeaderFun -> FilePath -> StorableArray Int a -> IO ()
writeArrayToFileWithHeader (initH, closeH) filepath arr = do
 handle <- openBinaryFile filepath WriteMode
 initH handle
 let (0,n) = bounds arr
     elemsize = sizeOfArrayElements arr
     bytesize = (1+n) * elemsize
 withStorableArray arr (\ptr -> hPutBuf handle ptr bytesize)
 closeH handle bytesize

writeArrayToFile :: Storable a => FilePath -> StorableArray Int a -> IO ()
writeArrayToFile filepath arr = do
 handle <- openBinaryFile filepath WriteMode
 let (0,n) = bounds arr
     elemsize = sizeOfArrayElements arr
     bytesize = (1+n) * elemsize
 withStorableArray arr (\ptr -> hPutBuf handle ptr bytesize)
 hClose handle
-}
---------------------------------------------------------------------------------------------------
-- Single Wav
---------------------------------------------------------------------------------------------------
openSingleInputWavFile :: Int -> FilePath -> IO (IO (), Either (IO (Maybe Int16)) (IO (Maybe (Int16,Int16))), Int)
openSingleInputWavFile n filepath = do
 handle <- openBinaryFile filepath ReadMode
 initReadWavHeader handle >>= maybe (error "unknown wav format") (\(isstereo, len1, len2, x) -> do
  hSeek handle AbsoluteSeek 44
  let stereo = do let bufsize = n*2
                  rpos   <- newIORef bufsize
                  rmax   <- newIORef bufsize
                  buf    <- newArray_ (0, bufsize-1)
                  let elemsize = sizeOfArrayElements buf
                      bytesize = elemsize * bufsize
                      next = readIORef rpos >>= \pos ->
                             readIORef rmax >>= \max ->
                             if pos >= max
                              then if max < bufsize
                                    then return Nothing
                                    else withStorableArray buf (\ptr -> hGetBuf handle ptr bytesize) >>= \bts ->
                                         writeIORef rmax (div bts elemsize) >>
                                         writeIORef rpos 0 >>
                                         next
                              else readArray buf pos >>= \v1 ->
                                   readArray buf (pos+1) >>= \v2 ->
                                   writeIORef rpos (pos + 2) >>
                                   return (Just (v1,v2))
                      close = hClose handle
                  return (close, Right next,div len2 4)
      mono   = do let bufsize = n
                  rpos   <- newIORef bufsize
                  rmax   <- newIORef bufsize
                  buf    <- newArray_ (0, bufsize-1)
                  let elemsize = sizeOfArrayElements buf
                      bytesize = elemsize * bufsize
                      next = readIORef rpos >>= \pos ->
                             readIORef rmax >>= \max ->
                             if pos >= max
                              then if max < bufsize
                                    then return Nothing
                                    else withStorableArray buf (\ptr -> hGetBuf handle ptr bytesize) >>= \bts ->
                                         writeIORef rmax (div bts elemsize) >>
                                         writeIORef rpos 0 >>
                                         next
                              else readArray buf pos >>= \v ->
                                   writeIORef rpos (pos + 1) >>
                                   return (Just v)
                      close = hClose handle
                  return (close, Left next,div len2 2)
  if isstereo then stereo else mono )
---------------------------------------------------------------------------------------------------
openSingleOutputWavFileMono :: Int -> FilePath -> IO (IO (), Int16 -> IO ())
openSingleOutputWavFileMono bufsize filepath = do
 handle <- openBinaryFile filepath WriteMode
 rpos   <- newIORef 0
 buf    <- newArray_ (0, bufsize-1)
 rcnt   <- newIORef 0
 let elemsize = sizeOfArrayElements buf
     bytesize = elemsize * bufsize
     write v = readIORef rpos >>= \pos ->
               if pos >= bufsize
                then withStorableArray buf (\ptr -> hPutBuf handle ptr bytesize) >>
                     modifyIORef rcnt (+bytesize) >>
                     writeIORef rpos 1 >> writeArray buf 0 v
                else writeArray buf pos v >> writeIORef rpos (pos+1)
     close = readIORef rpos >>= \pos ->
             withStorableArray buf (\ptr -> hPutBuf handle ptr (pos*elemsize)) >>
             readIORef rcnt >>= \cnt -> closeWriteWavHeader handle (elemsize*pos+cnt)
 initWriteWavHeaderMono handle
 return (close, write)
---------------------------------------------------------------------------------------------------
openSingleOutputWavFileStereo :: Int -> FilePath -> IO (IO (), (Int16, Int16) -> IO ())
openSingleOutputWavFileStereo bufsizearg filepath = do
 let bufsize | bufsizearg < 2  = 4
             | even bufsizearg = bufsizearg
             | otherwise       = bufsizearg + 1
 handle <- openBinaryFile filepath WriteMode
 rpos   <- newIORef 0
 buf    <- newArray_ (0, bufsize-1)
 rcnt   <- newIORef 0
 let elemsize = sizeOfArrayElements buf
     bytesize = elemsize * bufsize
     write (v1,v2) = readIORef rpos >>= \pos ->
                     if pos >= bufsize
                      then withStorableArray buf (\ptr -> hPutBuf handle ptr bytesize) >>
                           modifyIORef rcnt (+bytesize) >>
                           writeIORef rpos 2 >> writeArray buf 0 v1 >> writeArray buf 1 v2
                      else writeArray buf pos v1 >> writeArray buf (pos+1) v2 >> writeIORef rpos (pos+2)
     close = readIORef rpos >>= \pos ->
             withStorableArray buf (\ptr -> hPutBuf handle ptr (pos*elemsize)) >>
             readIORef rcnt >>= \cnt -> closeWriteWavHeader handle (elemsize*pos+cnt)
 initWriteWavHeaderStereo handle
 return (close, write)
---------------------------------------------------------------------------------------------------
-- Single Storable
---------------------------------------------------------------------------------------------------
openSingleInputFile :: Storable a => Int -> HeaderSize -> FilePath -> IO (IO (), IO (Maybe a), Int)
openSingleInputFile bufsize offset filepath = do
 handle <- openBinaryFile filepath ReadMode
 bytes <- hFileSize handle
 hSeek handle AbsoluteSeek $ fromIntegral offset
 rpos   <- newIORef bufsize
 rmax   <- newIORef bufsize
 buf    <- newArray_ (0, bufsize-1)
 let elemsize = sizeOfArrayElements buf
     bytesize = elemsize * bufsize
     next = readIORef rpos >>= \pos ->
            readIORef rmax >>= \max ->
            if pos >= max
             then if max < bufsize
                   then return Nothing
                   else withStorableArray buf (\ptr -> hGetBuf handle ptr bytesize) >>= \bts ->
                        writeIORef rmax (div bts elemsize) >>
                        writeIORef rpos 0 >>
                        next
             else readArray buf pos >>= \v ->
                  writeIORef rpos (pos + 1) >>
                  return (Just v)
     close = hClose handle
 return (close, next,div (fromIntegral bytes - offset) elemsize)
---------------------------------------------------------------------------------------------------
openSingleOutputFile :: Storable a => Int -> FilePath -> IO (IO (), a -> IO ())
openSingleOutputFile bufsize filepath = do
 handle <- openBinaryFile filepath WriteMode
 rpos   <- newIORef 0
 buf    <- newArray_ (0, bufsize-1)
 let elemsize = sizeOfArrayElements buf
     bytesize = elemsize * bufsize
     write v = readIORef rpos >>= \pos ->
               if pos >= bufsize
                then withStorableArray buf (\ptr -> hPutBuf handle ptr bytesize) >>
                     writeIORef rpos 1 >> writeArray buf 0 v
                else writeArray buf pos v >> writeIORef rpos (pos+1)
     close = readIORef rpos >>= \pos ->
             withStorableArray buf (\ptr -> hPutBuf handle ptr (pos*elemsize)) >>
             hClose handle
 return (close, write)
---------------------------------------------------------------------------------------------------
-- Buffered Storable
---------------------------------------------------------------------------------------------------
openInputFile :: Storable a => HeaderSize -> FilePath -> Int -> IO (IO (), StorableArray Int a -> Int -> IO Int, Int)
openInputFile offset filepath arrsize = do
 handle <- openBinaryFile filepath ReadMode
 hSeek handle AbsoluteSeek $ fromIntegral offset
 bytes <- hFileSize handle
 let nr_bytes  = fromIntegral bytes - offset
     elemsize  = inferSizeOfArrayElements'' read undefined
     read arr k = withStorableArray arr $ \ptr ->
                  hGetBuf handle ptr (k * elemsize) >>= return . flip div elemsize
     close     = hClose handle
 return $ (close, read, div nr_bytes elemsize)

openOutputFile :: Storable a => FilePath -> IO (IO (), StorableArray Int a -> Int -> IO ())
openOutputFile filepath = do
 handle <- openBinaryFile filepath WriteMode
 let elemsize = inferSizeOfArrayElements' write undefined
     write arr k = withStorableArray arr (\ptr -> hPutBuf handle ptr (k*elemsize))
     close = hClose handle
 return $ (close, write)

openOutputFileWithHeader :: Storable a => HeaderFun -> FilePath -> IO (IO (), StorableArray Int a -> Int -> IO ())
openOutputFileWithHeader (initH, closeH) filepath = do
 handle <- openBinaryFile filepath WriteMode
 rlen   <- newIORef 0
 initH handle
 let elemsize = inferSizeOfArrayElements' write undefined
     write arr k = modifyIORef rlen (+k) >> withStorableArray arr (\ptr -> hPutBuf handle ptr (k*elemsize))
     close =  readIORef rlen >>= \l -> closeH handle (elemsize*l)
 return $ (close, write)
---------------------------------------------------------------------------------------------------
-- Buffered Wav
---------------------------------------------------------------------------------------------------

openInputWavFile :: FilePath -> IO (IO (), StorableArray Int Int16 -> Int -> IO (), Int, Bool)
openInputWavFile filepath = do
 handle <- openBinaryFile filepath ReadMode
 initReadWavHeader handle >>= maybe (error "unknown wav format") (\(isstereo, len1, len2, x) -> do
  hSeek handle AbsoluteSeek 44
  let elemsize  = inferSizeOfArrayElements' read undefined
      read arr k = withStorableArray arr $ \ptr ->
                   hGetBuf handle ptr (k * elemsize) >> return ()
      close     = hClose handle
  return $ (close, read, div len2 2, isstereo) )

openOutputWavFileMono :: FilePath -> IO (IO (), StorableArray Int Int16 -> Int -> IO ())
openOutputWavFileMono filepath = do
 handle <- openBinaryFile filepath WriteMode
 rlen   <- newIORef 0
 let elemsize = inferSizeOfArrayElements' write undefined
     write arr k = withStorableArray arr (\ptr -> hPutBuf handle ptr (k*elemsize)) >> modifyIORef rlen (+k)
     close = readIORef rlen >>= \l -> closeWriteWavHeader handle (elemsize*l)
 initWriteWavHeaderMono handle
 return $ (close, write)

openOutputWavFileStereo :: FilePath -> IO (IO (), StorableArray Int Int16 -> Int -> IO ())
openOutputWavFileStereo filepath = do
 handle <- openBinaryFile filepath WriteMode
 rlen   <- newIORef 0
 let elemsize = inferSizeOfArrayElements' write undefined
     write arr k = withStorableArray arr (\ptr -> hPutBuf handle ptr (k*elemsize)) >> modifyIORef rlen (+k)
     close = readIORef rlen >>= \l -> closeWriteWavHeader handle (elemsize*l)
 initWriteWavHeaderStereo handle
 return $ (close, write)
---------------------------------------------------------------------------------------------------

ext: hs

max_stars_repo_name: aische/type-of-sound

max_issues_repo_name: aische/type-of-sound

max_forks_repo_name: aische/type-of-sound


text:
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE InstanceSigs         #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE UndecidableInstances #-}

module Network.Network
( Network(..) ) where

import           Network.Layer
import           Numeric.LinearAlgebra
import           System.Random

-- | A network is
class Network a where
  type Parameters :: *

  predict :: Vector Double -> a -> Vector Double
  createNetwork :: (RandomGen g) => RandomTransform -> g -> Parameters -> a

max_stars_repo_name: mckeankylej/hwRecog

max_issues_repo_name: mckeankylej/hwRecog

max_forks_repo_name: mckeankylej/hwRecog


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

module Foreign.BLAS.Level1 (
    BLAS1(..),
    ) where

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

import Foreign.VMath.VFractional( VFractional )
import Foreign.BLAS.Double
import Foreign.BLAS.Zomplex

-- | Types with vector-vector operations.
class (Storable a, VFractional a) => BLAS1 a where
    copy  :: Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    swap  :: Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    dotc  :: Int -> Ptr a -> Int -> Ptr a -> Int -> IO a
    dotu  :: Int -> Ptr a -> Int -> Ptr a -> Int -> IO a
    nrm2  :: Int -> Ptr a -> Int -> IO Double
    asum  :: Int -> Ptr a -> Int -> IO Double
    iamax :: Int -> Ptr a -> Int -> IO Int

    scal  :: Int -> a -> Ptr a -> Int -> IO ()

    axpy  :: Int -> a -> Ptr a -> Int -> Ptr a -> Int -> IO ()

    rotg  :: Ptr a -> Ptr a -> Ptr a -> Ptr a -> IO ()
    rot   :: Int -> Ptr a -> Int -> Ptr a -> Int -> Double -> Double -> IO ()

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

instance BLAS1 Double where
    copy n px incx py incy =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
            dcopy pn px pincx py pincy
    {-# INLINE copy #-}

    swap n px incx py incy =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
            dswap pn px pincx py pincy
    {-# INLINE swap #-}

    dotc n px incx py incy =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
            ddot pn px pincx py pincy
    {-# INLINE dotc #-}

    dotu = dotc
    {-# INLINE dotu #-}

    nrm2 n px incx =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
            dnrm2 pn px pincx
    {-# INLINE nrm2 #-}

    asum n px incx =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
            dasum pn px pincx
    {-# INLINE asum #-}

    iamax n px incx =
        withEnum n $ \pn ->
        withEnum incx $ \pincx -> do
            i <- idamax pn px pincx
            return $! fromEnum (i - 1)
    {-# INLINE iamax #-}

    axpy n alpha px incx py incy =
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
            daxpy pn palpha px pincx py pincy
    {-# INLINE axpy #-}

    scal n alpha px incx =
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
            dscal pn palpha px pincx
    {-# INLINE scal #-}

    rotg = drotg
    {-# INLINE rotg #-}

    rot n px incx py incy c s =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
        with c $ \pc ->
        with s $ \ps ->
            drot pn px pincx py pincy pc ps
    {-# INLINE rot #-}

instance BLAS1 (Complex Double) where
    copy n px incx py incy =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
            zcopy pn px pincx py pincy
    {-# INLINE copy #-}

    swap n px incx py incy =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
            zswap pn px pincx py pincy
    {-# INLINE swap #-}

    dotc n px incx py incy =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
        with 0 $ \pdotc -> do
            zdotc pdotc pn px pincx py pincy
            peek pdotc
    {-# INLINE dotc #-}

    dotu n px incx py incy =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
        with 0 $ \pdotu -> do
            zdotu pdotu pn px pincx py pincy
            peek pdotu
    {-# INLINE dotu #-}

    nrm2 n px incx =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
            znrm2 pn px pincx
    {-# INLINE nrm2 #-}

    asum n px incx =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
            zasum pn px pincx
    {-# INLINE asum #-}

    iamax n px incx =
        withEnum n $ \pn ->
        withEnum incx $ \pincx -> do
            i <- izamax pn px pincx
            return $! fromEnum (i - 1)
    {-# INLINE iamax #-}

    axpy n alpha px incx py incy =
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
            zaxpy pn palpha px pincx py pincy
    {-# INLINE axpy #-}

    scal n alpha px incx =
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum incx $ \pincx ->
            zscal pn palpha px pincx
    {-# INLINE scal #-}

    rotg = zrotg
    {-# INLINE rotg #-}

    rot n px incx py incy c s =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        withEnum incy $ \pincy ->
        with c $ \pc ->
        with s $ \ps ->
            zdrot pn px pincx py pincy pc ps
    {-# INLINE rot #-}

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE TupleSections #-}
module Geo.Computations.Trail
       ( -- * Types
         AvgMethod(..)
       , Selected(..)
       , PointGrouping
       , TransformGrouping
         -- * Utility Functions
       , isSelected
       , isNotSelected
       , onSelected
       , selLength
         -- * Trail Functions
         -- ** Queries
       , totalDistance
       , totalTime
       , avgSpeeds
       , slidingAverageSpeed
       , closestDistance
       , convexHull
         -- ** Transformations
       , bezierCurveAt
       , bezierCurve
       , linearTime
       , filterPoints
         -- ** Grouping Methods
       , betweenSpeeds
       , restLocations
       , spansTime
       , everyNPoints
         -- ** Group Transformations
       , intersectionOf
       , invertSelection
       , firstGrouping
       , lastGrouping
       , unionOf
       , refineGrouping
       , (/\), (\/)
         -- ** Composite Operations (Higher Level)
       , smoothRests
       , smoothTrail
        -- * Misc
       , bezierPoint
         ) where

import Text.Show.Functions ()
import Geo.Computations.Basic

import Control.Arrow (first)
import Control.Monad
import Data.Fixed (mod')
import Data.Function (on)
import Data.List as L
import Data.Maybe
import Data.Ord
import Data.Time

import Statistics.Function as F
import Statistics.Sample
import qualified Data.Vector.Unboxed as V

takeWhileEnd :: (a -> Bool) -> [a] -> (Maybe a, [a],[a])
takeWhileEnd p xs = go xs Nothing
  where
  go [] e = (e, [], [])
  go (a:as) e
    | p a = let (e',xs2,zs) = go as (Just a) in (e',a:xs2,zs)
    | otherwise = (e,[], a:as)

data AvgMethod c
  = AvgMean              -- ^ Obtain the 'mean' of the considered points
  | AvgHarmonicMean      -- ^ Obtain the 'harmonicMean'
  | AvgGeometricMean     -- ^ Obtain the 'geometricMean'
  | AvgMedian            -- ^ Obtain the median of the considered points
  | AvgEndPoints         -- ^ Compute the speed considering only the given endpoints
  | AvgMinOf [AvgMethod c] -- ^ Take the minimum of the speeds from the given methods
  | AvgWith ([c] -> Speed)

-- | @avgSpeeds n points@
-- Average speed using a window of up to @n@ seconds and averaging by taking the
-- Median ('AvgMedian').
avgSpeeds :: NominalDiffTime -> Trail Point -> [(UTCTime, Speed)]
avgSpeeds = slidingAverageSpeed AvgHarmonicMean

-- | @slidingAverageSpeed m n@ Average speed using a moving window of up to @n@ seconds
-- and an 'AvgMethod' of @m@.
slidingAverageSpeed :: AvgMethod Point -> NominalDiffTime -> Trail Point -> [(UTCTime, Speed)]
slidingAverageSpeed _ _ [] = []
slidingAverageSpeed m minTime xs =
  let pts   = map unSelect (spansTime minTime xs)
      spds  = map (getAvg m) pts
      times = map getAvgTimes pts
  in concatMap maybeToList $ zipWith (\t s -> fmap (,s) t) times spds
  where
  getTimeDiff a b = on (liftM2 diffUTCTime) pntTime a b

  --  getAvg :: [] -> AvgMethod -> Speed
  getAvg _ [] = 0
  getAvg _ [_] = 0
  getAvg m2 cs =
    let ss = getSpeedsV cs
    in case m2 of
        AvgMean -> mean ss
        AvgHarmonicMean  -> harmonicMean ss
        AvgGeometricMean -> geometricMean ss
        AvgMedian ->
          let ss' = F.sort $ getSpeedsV cs
              len = V.length ss'
              mid = len `div` 2
          in if V.length ss' < 3
             then mean ss'
             else if odd len then ss' V.! mid else mean (V.slice mid 2 ss')
        AvgEndPoints -> fromMaybe 0 $ speed (head cs) (last cs)
        AvgMinOf as -> minimum $ map (flip getAvg cs) as
        AvgWith f -> f cs
  getAvgTimes [] = Nothing
  getAvgTimes [x] = pntTime x
  getAvgTimes ps = getAvgTime (head ps) (last ps)
  getAvgTime a b = liftM2 addUTCTime (getTimeDiff b a) (pntTime a)
  getSpeedsV = V.fromList . getSpeeds
  getSpeeds zs = concatMap maybeToList $ zipWith speed zs (drop 1 zs)

-- | A PointGrouping is a function that selects segments of a trail.
--
-- Grouping point _does not_ result in deleted points. It is always true that:
--
--     forall g :: PointGrouping c -->
--     concatMap unSelect (g ts) == ts
--
-- The purpose of grouping is usually for later processing.  Any desire to drop
-- points that didn't meet a particular grouping criterion can be filled with
-- a composition with 'filter' (or directly via 'filterPoints').
type PointGrouping c = Trail c -> [Selected (Trail c)]

-- | Given a selection of coordinates, transform the selected
-- coordinates in some way (while leaving the non-selected
-- coordinates unaffected).
type TransformGrouping c = [Selected (Trail c)] -> [Selected (Trail c)]

-- | When grouping points, lists of points are either marked as 'Select' or 'NotSelect'.
data Selected a = Select {unSelect :: a} | NotSelect {unSelect :: a}
  deriving (Eq, Ord, Show)

isSelected :: Selected a -> Bool
isSelected (Select _) = True
isSelected _ = False

isNotSelected :: Selected a -> Bool
isNotSelected = not . isSelected

selLength :: Selected [a] -> Int
selLength = length . unSelect

onSelected :: (a -> b) -> (a -> b) -> Selected a -> b
onSelected f _ (Select a) = f a
onSelected _ g (NotSelect a) = g a

instance Functor Selected where
  fmap f (Select x) = Select $ f x
  fmap f (NotSelect x) = NotSelect $ f x

dropExact :: Int -> [Selected [a]] -> [Selected [a]]
dropExact _ [] = []
dropExact i (x:xs) =
  case compare (selLength x) i of
    EQ -> xs
    LT -> dropExact (i - selLength x) xs
    GT -> fmap (drop i) x : xs

-- | Groups trail segments into contiguous points within the speed
-- and all others outside of the speed.  The "speed" from point p(i)
-- to p(i+1) is associated with p(i) (except for the first speed
-- value, which is associated with both the first and second point)
betweenSpeeds :: Double -> Double -> PointGrouping Point
betweenSpeeds low hi ps =
  let spds = concatMap maybeToList $ zipWith speed ps (drop 1 ps)
      psSpds = [(p,s) | p <- ps, s <- maybeToList (listToMaybe spds) ++ spds]
      inRange x = x >= low && x <= hi
      chunk [] = []
      chunk xs@(x:_) =
        let op p = if inRange (snd x) then first Select . span p else first NotSelect . break p
            (r,rest) = op (inRange . snd) xs
        in r : chunk rest
  in map (fmap (map fst)) $ chunk psSpds

-- | A "rest point" means the coordinates remain within a given distance
-- for at least a particular amount of time.
restLocations :: Distance -> NominalDiffTime -> PointGrouping Point
restLocations d s ps =
  let consToFirst x [] = [NotSelect [x]]
      consToFirst x (a:as) = (fmap (x:) a) : as
      go [] [] = []
      go [] nonRests = [NotSelect $ reverse nonRests]
      go (a:as) nonRests =
        case takeWhileEnd ((<=) d . distance a) as of
          (Just l, close, far) ->
            case (pntTime a, pntTime l) of
              (Just t1, Just t2) ->
                let diff = diffUTCTime t2 t1
                in if diff >= s then NotSelect (reverse nonRests) : Select (a:close) : go far [] else go as (a:nonRests)
              _ -> consToFirst a $ go as nonRests
          _ -> consToFirst a $ go as nonRests
  in go ps []

-- | Chunks points into groups spanning at most the given time
-- interval.
spansTime :: NominalDiffTime -> PointGrouping Point
spansTime n ps =
  let times  = mkTimePair ps
      chunk [] = []
      chunk xs@(x:_) =
        let (good,rest) = span ((<= addUTCTime n (snd x)) . snd) xs
        in if null good then [xs] else good : chunk rest
  in map (Select . map fst) $ chunk times

-- | Intersects the given groupings
intersectionOf :: [PointGrouping Point] -> PointGrouping Point
intersectionOf gs ps =
  let groupings = map ($ ps) gs
      -- chunk :: [[Selected [pnts]]] -> pnts -> [pnts]
      chunk _ [] = []
      chunk ggs xs =
        let minLen = max 1 . minimum . concatMap (take 1) $ map (map selLength) ggs   -- FIXME this is all manner of broken
            sel = if all isSelected (concatMap (take 1) ggs) then Select else NotSelect
            (c,rest) = splitAt minLen xs
        in sel c : chunk (filter (not . null) $ map (dropExact minLen) ggs) rest
  in chunk groupings ps

-- | Union all the groupings
unionOf :: [PointGrouping Point] -> PointGrouping Point
unionOf gs ps =
  let groupings = map ($ ps) gs
      chunk _ [] = []
      chunk ggs xs =
        let getSegs = concatMap (take 1)
            segs = getSegs ggs
            len =
              if any isSelected segs
                 then max 1 . maximum . getSegs . map (map selLength) . map (filter isSelected) $ ggs
                 else max 1 . minimum . getSegs . map (map selLength) $ ggs
            sel = if any isSelected segs then Select else NotSelect
            (c,rest) = splitAt len xs
        in sel c : chunk (filter (not . null) $ map (dropExact len) ggs) rest
  in chunk groupings ps

-- | Intersection binary operator
(/\) :: [Selected (Trail a)] -> TransformGrouping a
(/\) _ [] = []
(/\) [] _ = []
(/\) xsL@(Select x:_) ysL@(Select y:_) =
  let z = if length x < length y then x else y
      xs' = selListDrop (length z) xsL
      ys' = selListDrop (length z) ysL
  in Select z : (xs' /\ ys')
(/\) xs (NotSelect y:ys) = NotSelect y : (selListDrop (length y) xs /\ ys)
(/\) (NotSelect x:xs) ys = NotSelect x : (xs /\ selListDrop (length x) ys)

-- | Union binary operator
(\/) :: [Selected (Trail a)] -> TransformGrouping a
(\/) xs [] = xs
(\/) [] ys = ys
(\/) (Select x:xs) (Select y : ys) =
  let xLen = length x
      yLen = length y
  in if xLen < yLen
       then (Select y :) (selListDrop (yLen - xLen) xs \/ ys)
       else (Select x :) (xs \/ selListDrop (xLen - yLen) ys)
(\/) (Select x:_) ys = Select x : selListDrop (length x) ys
(\/) xs (Select y:_) = Select y : selListDrop (length y) xs
(\/) xsL@(NotSelect x:xs) ysL@(NotSelect y:ys) =
  let xLen = length x
      yLen = length y
  in if xLen < yLen
        then (NotSelect x:) (xs \/ selListDrop xLen ysL)
        else (NotSelect y:) (selListDrop yLen xsL \/ ys)

selListDrop :: Int -> [Selected [a]] -> [Selected [a]]
selListDrop 0 xs = xs
selListDrop _ [] = []
selListDrop n (x:xs) =
  let x' = drop n (unSelect x)
  in fmap (const x') x : selListDrop (n - (selLength x - length x')) xs

-- | Inverts the selected/nonselected segments
invertSelection :: TransformGrouping a
invertSelection = map (onSelected NotSelect Select)

-- | @firstGrouping f ps@ only the first segment remains 'Select'ed, and only
-- if it was already selected by @f@.
firstGrouping ::  TransformGrouping a
firstGrouping ps = take 1 ps ++ map (NotSelect . unSelect) (drop 1 ps)

-- | Only the last segment, if any, is selected (note: the current
-- implementation is inefficient, using 'reverse')
lastGrouping ::  TransformGrouping a
lastGrouping ps  = let ps' = reverse ps in reverse $ take 1 ps' ++ map (NotSelect . unSelect) (drop 1 ps')

-- | Chunks the trail into groups of N points
everyNPoints ::  Int -> PointGrouping a
everyNPoints n ps
  | n <= 0 = [NotSelect ps]
  | otherwise = go ps
    where
      go [] = []
      go xs = let (h,t) = splitAt n xs in Select h : go t

-- | For every selected group, refine the selection using the second
-- grouping method.  This differs from 'IntersectionOf' by restarting
-- the second grouping algorithm at the beginning each group selected
-- by the first algorithm.
refineGrouping ::  PointGrouping a -> TransformGrouping a
refineGrouping b = concatMap (onSelected b (\x -> [NotSelect x]))

-- | Remove all points that remain 'NotSelect'ed by the given grouping algorithm.
filterPoints :: PointGrouping a -> Trail a -> Trail a
filterPoints g = concatMap unSelect . filter isSelected . g

-- Extract the time from each coordinate.  If no time is available then
-- the coordinate is dropped!
mkTimePair :: Trail Point -> [(Point,UTCTime)]
mkTimePair xs =
  let timesM = map (\x-> fmap (x,) $ pntTime x) xs
  in concatMap maybeToList timesM

-- | Construct a bezier curve using the provided trail.  Construct a
-- new trail by sampling the given bezier curve at the given times.
-- The current implementation assumes the times of the input
-- coordinates are available and all equal (Ex: all points are 5
-- seconds apart), the results will be poor if this is not the case!
bezierCurveAt :: [UTCTime] -> Trail Point -> Trail Point
bezierCurveAt _ [] = []
bezierCurveAt selectedTimes xs =
  let timesDef = mkTimePair xs
      end = last timesDef
      top = head timesDef
      tTime  = diffUTCTime (snd end) (snd top)
      times = if null selectedTimes then map snd timesDef else selectedTimes
      diffTimes = [diffUTCTime t (snd top) / tTime | t <- times]
      queryTimes = map realToFrac diffTimes
  in if tTime <= 0 || any (\x -> x < 0 || x > 1) queryTimes
        then xs -- error "bezierCurveAt has a out-of-bound time!"
        else
         if null timesDef || any (\x -> x < 0 || x > 1) queryTimes
         then xs
         else let curvePoints = (map (bezierPoint xs) queryTimes)
                  newTimes = [addUTCTime t (snd top) | t <- diffTimes]
              in zipWith (\t p -> p { pntTime = Just t}) newTimes curvePoints

bezierPoint :: [Point] -> Double -> Point
bezierPoint pnts t   = go pnts
  where
  go [] = error "GPS Package: Can not create a bezier point from an empty list"
  go [p] = p
  go ps = interpolate (go (init ps)) (go (tail ps)) t

-- | Interpolate selected points onto a bezier curve.  Note this gets
-- exponentially more expensive with the length of the segment being
-- transformed - it is not advisable to perform this operation on
-- trail segements with more than ten points!
bezierCurve ::  [Selected (Trail Point)] -> Trail Point
bezierCurve = concatMap (onSelected (bezierCurveAt []) Prelude.id)

-- | Filter out any points that go backward in time (thus must not be
-- valid if this is a trail)
linearTime :: [Point] -> [Point]
linearTime [] = []
linearTime (p:ps) = go (pntTime p) ps
  where
  go _ [] = []
  go t (x:xs) = if pntTime x < t then go t xs else x : go (pntTime x) xs

-- | Return the closest distance between two trails (or Nothing if a
-- trail is empty).  Inefficient implementation:
-- O( (n * m) * log (n * m) )
closestDistance :: Trail Point -> Trail Point -> Maybe Distance
closestDistance as bs = listToMaybe $ L.sort [distance a b | a <- as, b <- bs]

-- | Find the total distance traveled
totalDistance :: [Point] -> Distance
totalDistance as = sum $ zipWith distance as (drop 1 as)

totalTime :: Trail Point -> NominalDiffTime
totalTime [] = 0
totalTime xs@(x:_) = fromMaybe 0 $ liftM2 diffUTCTime (pntTime $ last xs) (pntTime x)

-- | Uses Grahams scan to compute the convex hull of the given points.
-- This operation requires sorting of the points, so don't try it unless
-- you have notably more memory than the list of points will consume.
convexHull :: [Point] -> [Point]
convexHull lst =
        let frst = southMost lst
        in case frst of
                Nothing -> []
                Just f  ->
                     let sorted = L.sortBy (comparing (eastZeroHeading f)) (filter (/= f) lst)
                     in case sorted of
                        (a:b:cs) -> grahamScan (b:a:f:[]) cs
                        cs       -> f : cs
  where
  grahamScan [] _ = []
  grahamScan ps [] = ps
  grahamScan (x:[]) _ = [x]
  grahamScan (p2:p1:ps) (x:xs) =
        case turn p1 p2 x of
                LeftTurn  -> grahamScan (x:p2:p1:ps) xs
                Straight  -> grahamScan (x:p2:p1:ps) xs
                _         -> grahamScan (p1:ps) (x:xs)

eastZeroHeading :: Point -> Point -> Heading
eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s

data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum)

turn :: Point -> Point -> Point -> Turn
turn a b c =
        let h1 = eastZeroHeading a b
            h2 = eastZeroHeading b c
            d  = h2 - h1
        in if d >= 0 && d < pi then LeftTurn else RightTurn

-- | Find the southmost point
southMost :: [Point] -> Maybe Point
southMost []  = Nothing
southMost cs = Just . minimumBy (comparing pntLat) $ cs

---------- COMPOSIT OPERATIONS ---------------
-- These operations are simply implemented using the previously
-- defined functions. They can serve either for concise use for novice
-- users or as instructional examples.
------------------------------------------

-- | Smooth points with rest areas using a bezierCurve.
--
-- Parameters: rest for 1 minute within 30 meters get smoothed
-- in a bezier curve over every 8 points.
smoothRests :: Trail Point -> Trail Point
smoothRests = bezierCurve . refineGrouping (everyNPoints 8) . restLocations 30 60

-- | Smooth every 7 points using a bezier curve
smoothTrail :: Trail Point -> Trail Point
smoothTrail = gSmoothSome 7

-- | Smooth every n points using a bezier curve
gSmoothSome :: Int -> Trail Point -> Trail Point
gSmoothSome n = bezierCurve . everyNPoints n

max_stars_repo_name: beezee/gps

max_issues_repo_name: beezee/gps

max_forks_repo_name: beezee/gps


text:
module Fractals
    ( julia
    , mandlebrot
    , drawMandlebrot
    , Coordinate
    , drawJuliaAtPoint
    , drawJuliaDefault
    )
where

import           Data.Complex
import           Graphics.GD

type Coordinate = (Double, Double)

iterateSet
    :: Int
    -> Complex Double -- c
    -> Complex Double -- z
    -> Int
    -> Int
iterateSet maxIter c z iter = if iter > maxIter
    then 0
    else
        let z' = z ^ 2 + c
        in  if magnitude z' > 2
                then iter
                else iterateSet maxIter c z' (iter + 1)

julia :: Complex Double -> Complex Double -> Int
julia c z = iterateSet 64 c z 0

mandlebrot :: Complex Double -> Complex Double -> Int
mandlebrot c z = iterateSet 127 c z 0

drawMandlebrot :: Coordinate -> Color
drawMandlebrot (x, y) = pointColour 127 $ mandlebrot (x :+ y) (0 :+ 0)

drawJuliaAtPoint :: Complex Double -> Coordinate -> Color
drawJuliaAtPoint c (x, y) = pointColour 64 $ julia c (x :+ y)

drawJuliaDefault :: Coordinate -> Color
drawJuliaDefault (x, y) = drawJuliaAtPoint ((-0.423) :+ 0.745) (x, y)

pointColour :: Int -> Int -> Color
pointColour _ 0 = rgb 255 255 255
pointColour maxIter x =
    let colour = iterToColour maxIter x in rgb colour colour colour

iterToColour :: Int -> Int -> Int
iterToColour maxIter iter =
    floor (255.0 - 255.0 * fromIntegral iter / fromIntegral maxIter)

max_stars_repo_name: EwanGilligan/FractalGenerator

max_issues_repo_name: EwanGilligan/FractalGenerator

max_forks_repo_name: EwanGilligan/FractalGenerator


text:
module Sound.Hommage.Filter
 (
 -- * Filters
   lowpass
 , highpass
 , bandpass
 , stretchpass
 , morphpass
-- , vocoder

-- , lowpassrect

 -- * FilterSpec
-- , specfilter
 , FilterSpec
 , shiftFilterSpec

 -- * FilterSpec and Lists
 , filterSpecToList
 , filterSpecsToLists
 , interpolFilterLists

 -- * Other
 , ffreq
 , ffreqi
 , average512

 -- * Filterfunctions
 , lowpassFS
 , highpassFS
 , bandpassFS
 , lowp0FS
 , lowp1FS
 , bandp0FS
 , bandp1FS
 , highp0FS
 , highp1FS

 )
 where

--import Sound.Hommage.DFTFilter
import Data.Complex

--import Sound.Hommage.FFT
---------------------------------------------------------------------------------------------------
lowpass :: [Double] -> [Double] -> [Double]
lowpass s c = concat $ map filterSpecToList $
 zipWith lowpassFS (average512 s) (average512 c)

highpass :: [Double] -> [Double] -> [Double]
highpass s c = concat $ map filterSpecToList $
 zipWith highpassFS (average512 s) (average512 c)

bandpass :: [Double] -> [Double] -> [Double] -> [Double]
bandpass r s c = concat $ map filterSpecToList $
 zipWith3 bandpassFS (average512 r) (average512 s) (average512 c)

stretchpass:: FilterSpec -> [Double] -> [Double]
stretchpass fs c = concat $ map filterSpecToList $
 map (flip shiftFilterSpec fs) $ map (2.0**) $ average512 c

morphpass:: FilterSpec -> FilterSpec -> [Double] -> [Double]
morphpass f1 f2 c =
 let l1 = filterSpecToList f1
     l2 = filterSpecToList f2
 in concat $ map (interpolFilterLists l1 l2) $ average512 c
---------------------------------------------------------------------------------------------------
{-

-- | first argument is width of filter, second is cutoff
lowpassfilter :: [Double] -> [Double] -> [Double] -> [Double]
lowpassfilter s c =

lowpassrect :: [Double] -> [Double] -> [Double] -> [Double]
lowpassrect s c =
 let i = ([]:) $ map filterSpecToList $ zipWith lowpassFS (average512 s) (average512 c)
 in rect_filter' 9 i

{-
specfilterFFT :: [FilterSpec] -> [Double] -> [Double]
specfilterFFT fs = fftt 8 (concat . zipWith fun fs )

fun :: FilterSpec -> [Complex Double] -> [Complex Double]
fun fs ((c0 :+ _) : cs) = zipWith f ((c0 :+ 0) : cs) (filterSpecToList fs)
 where
  f (x:+y) v = (x * v) :+ (y * v)
-}

-- | first argument is width of filter, second is cutoff
highpassfilter :: [Double] -> [Double] -> [Double] -> [Double]
highpassfilter s c =
 let i = concat $ map filterSpecToList $ zipWith highpassFS (average512 s) (average512 c)
 in dftfilter i

-- | first argument is width of filter, second is slope, third is cutoff
bandpassfilter :: [Double] -> [Double] -> [Double] -> [Double] -> [Double]
bandpassfilter r s c =
 let i = concat $ map filterSpecToList $ zipWith3 bandpassFS (average512 r) (average512 s) (average512 c)
 in dftfilter i

stretchfilter :: FilterSpec -> [Double] -> [Double] -> [Double]
stretchfilter fs c =
 let i = concat $ map filterSpecToList $ map (flip shiftFilterSpec fs) $ map (2.0**) $ average512 c
 in dftfilter i

morphfilter :: FilterSpec -> FilterSpec -> [Double] -> [Double] -> [Double]
morphfilter f1 f2 c =
 let l1 = filterSpecToList f1
     l2 = filterSpecToList f2
     i = concat $ map (interpolFilterLists l1 l2) $ average512 c
 in dftfilter i

---------------------------------------------------------------------------------------------------
specfilter :: [FilterSpec] -> [Double] -> [Double]
specfilter fs =
 let i = concat $ map filterSpecToList fs
 in dftfilter i
-}
---------------------------------------------------------------------------------------------------
-- | for filter coefficient number: (0.0 .. 1.0) -> (1.0 .. 512.0)
ffreq :: Double -> Double
ffreq k = r
 where
  k' = abs k
  r | k' > 1.0  = 512.0
    | otherwise = 2 ** (9.0 * k')

-- | floor of 'ffreq'
ffreqi :: Double -> Int
ffreqi = floor . ffreq

-- | the resulting list is 1 \/ 512 as long as the input list.
--   512 elements are read and their average is the next output value.
average512 :: [Double] -> [Double]
average512 = loop 512 0.0
 where
  loop 0 s xs     = (s / 512.0) : loop 512 0.0 xs
  loop k s (x:xs) = loop (k-1) (s+x) xs
  loop _ _ _      = []
--mkFilterFun :: (Double -> FilterSpec) -> (Double, Double) -> FilterFun Double
--mkFilterFun filt (cutoff, fmod) = withFilterEnv . mkFilterEnv . filt . (+cutoff) . (*fmod)
---------------------------------------------------------------------------------------------------
-- | Filtering starts at coeff nr 0 and value 0.0 (constant coeff is always zero).
--   The elements in FilterSpec define the next value and how many
--   coeffs it takes to reach this value.
type FilterSpec   = [(Int, Double)]

-- | 0..1: lower, 1..: higher
shiftFilterSpec :: Double -> FilterSpec -> FilterSpec
shiftFilterSpec d fs = map f fs
 where
  f (n,v) = (round (d * fromIntegral n), v)

-- | Converting a FilterSpec to a list of 512 Doubles
filterSpecToList :: FilterSpec -> [Double]
filterSpecToList xs = take 512 $ loop 0.0 xs
 where
  loop c []          = repeat c
  loop c ((n,d): ls) = let f = (d-c) / fromIntegral n in take n (iterate (+f) c) ++ loop d ls

-- | Converting and interpolating a sequence of FilterSpecs to a list of lists. The first elemet of the
-- tuple, an Int, describes how many lists it takes to reach the given FilterSpec.
filterSpecsToLists :: [(Int, FilterSpec)] -> [[Double]]
filterSpecsToLists = loop (filterSpecToList  [])
 where
  loop l []         = []
  loop l ((n,f):xs) = let f' = filterSpecToList f
                      in fun n l f' ++ loop f' xs
  fun n f1 f2       = let d = 1.0 / fromIntegral n
                      in take n $ map (interpolFilterLists f1 f2) $ iterate (+d) 0.0

interpolFilterLists :: [Double] -> [Double] -> Double -> [Double]
interpolFilterLists f1 f2 n | n >= 1.0  = f2
                            | n <= 0.0  = f1
                            | otherwise = let n' = 1.0 - n in
                                          zipWith (\a b -> a * n' + b * n) f1 f2

---------------------------------------------------------------------------------------------------
lowpassFS :: Double -> Double -> FilterSpec
lowpassFS d c = [(1,1.0), (round fd, 1.0), (round (fd * abs d), 0.0)]
 where
  fd = ffreq c

highpassFS :: Double -> Double -> FilterSpec
highpassFS d c = [(round fd, 0.0), (round (fd * abs d), 1.0)]
 where
  fd = ffreq c

-- | range, curve, cutoff
bandpassFS :: Double -- ^ width of freq-window
           -> Double --
           -> Double -- ^ cutoff
           -> FilterSpec
bandpassFS r s f = [(l1, 0.0), (l2, 1.0), (l3, 1.0), (l4, 0.0)]
 where

  m1 = abs r + 1
  m2 = m1 + abs s
  a = ffreqi (f / m2)
  b = ffreqi (f / m1)
  c = ffreqi (f * m1)
  d = ffreqi (f * m2)

  l1 = a
  l2 = b - a
  l3 = c - b
  l4 = d - c

bandpassFS' :: Double -- ^ width of freq-window
           -> Double --
           -> Double -- ^ cutoff
           -> FilterSpec
bandpassFS' r s f = [(1 + round a, 0.0), (1 + round b, 1.0), (1 + round c, 1.0), (1 + round d, 0.0)]
 where
  r' = 1.0 + abs r
  s' = 1.0 + abs s
  fd = ffreq f
  a  = fd
  b  = fd * s' - a
  c  = b * r'  - b
  d  = c * s'  - c
---------------------------------------------------------------------------------------------------
lowp0FS :: Double -> FilterSpec
lowp0FS f = [(0, 1.0), (fr, 1.0), (1, 0.0)]
 where
  fr = ffreqi f

lowp1FS :: Double -> FilterSpec
lowp1FS f = [(0, 1.0), (fr, 1.0), (fr, 0.0)]
 where
  fr = ffreqi f

bandp0FS :: Double -> FilterSpec
bandp0FS f = [(f0-1, 0.0), (f0, 1.0), (fr, 1.0), (1, 0.0)]
 where
  fr = ffreqi f
  f0 = fr - (div fr 2) + 1

bandp1FS :: Double -> FilterSpec
bandp1FS f = [(fh, 0.0), (fh, 1.0), (fh, 0.0)]
 where
  fr = ffreqi f
  fh = div fr 2

highp0FS :: Double -> FilterSpec
highp0FS f = [(fr, 0.0), (1, 1.0)]
 where
  fr = ffreqi f

highp1FS :: Double -> FilterSpec
highp1FS f = [(fh, 0.0), (fh, 1.0)]
 where
  fr = ffreqi f
  fh = div fr 2
---------------------------------------------------------------------------------------------------
{-
-- | first argument is control signal
vocoder :: [Double] -> [Double] -> [Double]
vocoder xs ys = dftsynthese $ zipWith f (dftanalyse xs) (dftanalyse ys)
 where
  f x (a :+ b) = let z = magnitude x in (a*z :+ b*z)
-}
---------------------------------------------------------------------------------------------------

max_stars_repo_name: aische/hommage

max_issues_repo_name: aische/hommage

max_forks_repo_name: aische/hommage


text:
-- https://drive.google.com/open?id=0BxJP-hCBgo5OcU9td2Fhc2xHek0

module Numeric.Combinations
  ( combIndex
  , combIndexSorted
  )
  where

import Control.Monad.ST (runST)

import Data.Vector.Unboxed ((!))
import qualified Data.Vector.Unboxed as UV
import Data.Vector.Algorithms.Merge (sort)

import Numeric.SpecFunctions (choose)

import Debug.Trace

-- | Calculates the combination index. Assumes that neither alpha nor x
--   contain any duplicate elements.
combIndex :: (UV.Unbox a, Ord a)
          => UV.Vector a -- Alphabet of length N
          -> UV.Vector a -- element of length K
          -> Maybe Int -- Index in C_K^A
combIndex alpha x | UV.length x > UV.length alpha = Nothing -- Require that the alphabet is bigger than the element!
combIndex alpha x = combIndexSorted alpha_sorted x_sorted where
  alpha_sorted = runST $ do
    m_alpha <- UV.thaw alpha
    sort m_alpha
    UV.unsafeFreeze m_alpha
  x_sorted = runST $ do
    m_x <- UV.thaw x
    sort m_x
    UV.unsafeFreeze m_x

-- | Calculates the combination index assuming both alpha and x are pre-sorted.
--   This will give incorrect answers if the inputs are not sorted!
combIndexSorted :: (UV.Unbox a, Ord a)
                => UV.Vector a -- Alphabet of length N
                -> UV.Vector a -- element of length K
                -> Maybe Int -- Index in C_K^A
combIndexSorted alpha x | UV.length x > UV.length alpha = Nothing
combIndexSorted alpha x = do
    x_idx <- m_x_idx
    return $ UV.sum
           . UV.map (flip inner_sum x_idx)
           $ UV.enumFromN 0 cap_k
  where
    m_x_idx = UV.mapM (flip UV.elemIndex alpha) x
    cap_k = UV.length x
    cap_n = UV.length alpha
    i_k (-1) _ = -1
    i_k k x_idx = x_idx ! k
    inner_sum k x_idx =
          UV.sum
        . UV.map (\q -> floor
                      $ choose (cap_n - (i1 + 1) - q) (cap_k - k - 1))
        $ UV.enumFromN 1 (i0 - (i1 + 1))
      where
        i0 = i_k k x_idx
        i1 = i_k (k-1) x_idx