text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE DefaultSignatures #-}
-----------------------------------------------------------------------------
-- |
-- Copyright   :  (C) 2012-2015 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  provisional
-- Portability :  portable
--
-- Operations on free vector spaces.
-----------------------------------------------------------------------------
module Linear.Vector
  ( Additive(..)
  , E(..)
  , negated
  , (^*)
  , (*^)
  , (^/)
  , sumV
  , basis
  , basisFor
  , scaled
  , outer
  , unit
  ) where

import Control.Applicative
import Control.Lens
import Data.Complex
import Data.Foldable as Foldable (forM_, foldl')
import Data.Functor.Compose
import Data.Functor.Product
import Data.HashMap.Lazy as HashMap
import Data.Hashable
import Data.IntMap as IntMap
import Data.Map as Map
import qualified Data.Vector as Vector
import Data.Vector (Vector)
import qualified Data.Vector.Mutable as Mutable
import GHC.Generics
import Linear.Instances ()

-- $setup
-- >>> import Linear.V2

-- | Basis element
newtype E t = E { el :: forall x. Lens' (t x) x }

infixl 6 ^+^, ^-^
infixl 7 ^*, *^, ^/

class GAdditive f where
  gzero :: Num a => f a
  gliftU2 :: (a -> a -> a) -> f a -> f a -> f a
  gliftI2 :: (a -> b -> c) -> f a -> f b -> f c

instance GAdditive U1 where
  gzero = U1
  {-# INLINE gzero #-}
  gliftU2 _ U1 U1 = U1
  {-# INLINE gliftU2 #-}
  gliftI2 _ U1 U1 = U1
  {-# INLINE gliftI2 #-}

instance (GAdditive f, GAdditive g) => GAdditive (f :*: g) where
  gzero = gzero :*: gzero
  {-# INLINE gzero #-}
  gliftU2 f (a :*: b) (c :*: d) = gliftU2 f a c :*: gliftU2 f b d
  {-# INLINE gliftU2 #-}
  gliftI2 f (a :*: b) (c :*: d) = gliftI2 f a c :*: gliftI2 f b d
  {-# INLINE gliftI2 #-}

instance (Additive f, GAdditive g) => GAdditive (f :.: g) where
  gzero = Comp1 $ gzero <$ (zero :: f Int)
  {-# INLINE gzero #-}
  gliftU2 f (Comp1 a) (Comp1 b) = Comp1 $ liftU2 (gliftU2 f) a b
  {-# INLINE gliftU2 #-}
  gliftI2 f (Comp1 a) (Comp1 b) = Comp1 $ liftI2 (gliftI2 f) a b
  {-# INLINE gliftI2 #-}

instance Additive f => GAdditive (Rec1 f) where
  gzero = Rec1 zero
  {-# INLINE gzero #-}
  gliftU2 f (Rec1 g) (Rec1 h) = Rec1 (liftU2 f g h)
  {-# INLINE gliftU2 #-}
  gliftI2 f (Rec1 g) (Rec1 h) = Rec1 (liftI2 f g h)
  {-# INLINE gliftI2 #-}

instance GAdditive f => GAdditive (M1 i c f) where
  gzero = M1 gzero
  {-# INLINE gzero #-}
  gliftU2 f (M1 g) (M1 h) = M1 (gliftU2 f g h)
  {-# INLINE gliftU2 #-}
  gliftI2 f (M1 g) (M1 h) = M1 (gliftI2 f g h)
  {-# INLINE gliftI2 #-}

instance GAdditive Par1 where
  gzero = Par1 0
  gliftU2 f (Par1 a) (Par1 b) = Par1 (f a b)
  {-# INLINE gliftU2 #-}
  gliftI2 f (Par1 a) (Par1 b) = Par1 (f a b)
  {-# INLINE gliftI2 #-}

-- | A vector is an additive group with additional structure.
class Functor f => Additive f where
  -- | The zero vector
  zero :: Num a => f a
#ifndef HLINT
  default zero :: (GAdditive (Rep1 f), Generic1 f, Num a) => f a
  zero = to1 gzero
#endif

  -- | Compute the sum of two vectors
  --
  -- >>> V2 1 2 ^+^ V2 3 4
  -- V2 4 6
  (^+^) :: Num a => f a -> f a -> f a
  (^+^) = liftU2 (+)
  {-# INLINE (^+^) #-}

  -- | Compute the difference between two vectors
  --
  -- >>> V2 4 5 ^-^ V2 3 1
  -- V2 1 4
  (^-^) :: Num a => f a -> f a -> f a
  x ^-^ y = x ^+^ negated y

  -- | Linearly interpolate between two vectors.
  lerp :: Num a => a -> f a -> f a -> f a
  lerp alpha u v = alpha *^ u ^+^ (1 - alpha) *^ v
  {-# INLINE lerp #-}

  -- | Apply a function to merge the 'non-zero' components of two vectors, unioning the rest of the values.
  --
  -- * For a dense vector this is equivalent to 'liftA2'.
  --
  -- * For a sparse vector this is equivalent to 'unionWith'.
  liftU2 :: (a -> a -> a) -> f a -> f a -> f a
#ifndef HLINT
  default liftU2 :: Applicative f => (a -> a -> a) -> f a -> f a -> f a
  liftU2 = liftA2
  {-# INLINE liftU2 #-}
#endif

  -- | Apply a function to the components of two vectors.
  --
  -- * For a dense vector this is equivalent to 'liftA2'.
  --
  -- * For a sparse vector this is equivalent to 'intersectionWith'.
  liftI2 :: (a -> b -> c) -> f a -> f b -> f c
#ifndef HLINT
  default liftI2 :: Applicative f => (a -> b -> c) -> f a -> f b -> f c
  liftI2 = liftA2
  {-# INLINE liftI2 #-}
#endif

instance (Additive f, Additive g) => Additive (Product f g) where
  zero = Pair zero zero
  liftU2 f (Pair a b) (Pair c d) = Pair (liftU2 f a c) (liftU2 f b d)
  liftI2 f (Pair a b) (Pair c d) = Pair (liftI2 f a c) (liftI2 f b d)
  Pair a b ^+^ Pair c d = Pair (a ^+^ c) (b ^+^ d)
  Pair a b ^-^ Pair c d = Pair (a ^-^ c) (b ^-^ d)
  lerp alpha (Pair a b) (Pair c d) = Pair (lerp alpha a c) (lerp alpha b d)

instance (Additive f, Additive g) => Additive (Compose f g) where
  zero = Compose $ zero <$ (zero :: f Int)
  {-# INLINE zero #-}
  Compose a ^+^ Compose b = Compose $ liftU2 (^+^) a b
  {-# INLINE (^+^) #-}
  Compose a ^-^ Compose b = Compose $ liftU2 (^-^) a b
  {-# INLINE (^-^) #-}
  liftU2 f (Compose a) (Compose b) = Compose $ liftU2 (liftU2 f) a b
  {-# INLINE liftU2 #-}
  liftI2 f (Compose a) (Compose b) = Compose $ liftI2 (liftI2 f) a b
  {-# INLINE liftI2 #-}

instance Additive ZipList where
  zero = ZipList []
  {-# INLINE zero #-}
  liftU2 f (ZipList xs) (ZipList ys) = ZipList (liftU2 f xs ys)
  {-# INLINE liftU2 #-}
  liftI2 = liftA2
  {-# INLINE liftI2 #-}

instance Additive Vector where
  zero = mempty
  {-# INLINE zero #-}
  liftU2 f u v = case compare lu lv of
    LT | lu == 0   -> v
       | otherwise -> Vector.modify (\ w -> Foldable.forM_ [0..lu-1] $ \i -> Mutable.unsafeWrite w i $ f (Vector.unsafeIndex u i) (Vector.unsafeIndex v i)) v
    EQ -> Vector.zipWith f u v
    GT | lv == 0   -> u
       | otherwise -> Vector.modify (\ w -> Foldable.forM_ [0..lv-1] $ \i -> Mutable.unsafeWrite w i $ f (Vector.unsafeIndex u i) (Vector.unsafeIndex v i)) u
    where
      lu = Vector.length u
      lv = Vector.length v
  {-# INLINE liftU2 #-}
  liftI2 = Vector.zipWith
  {-# INLINE liftI2 #-}

instance Additive Maybe where
  zero = Nothing
  {-# INLINE zero #-}
  liftU2 f (Just a) (Just b) = Just (f a b)
  liftU2 _ Nothing ys = ys
  liftU2 _ xs Nothing = xs
  {-# INLINE liftU2 #-}
  liftI2 = liftA2
  {-# INLINE liftI2 #-}

instance Additive [] where
  zero = []
  {-# INLINE zero #-}
  liftU2 f = go where
    go (x:xs) (y:ys) = f x y : go xs ys
    go [] ys = ys
    go xs [] = xs
  {-# INLINE liftU2 #-}
  liftI2 = Prelude.zipWith
  {-# INLINE liftI2 #-}

instance Additive IntMap where
  zero = IntMap.empty
  {-# INLINE zero #-}
  liftU2 = IntMap.unionWith
  {-# INLINE liftU2 #-}
  liftI2 = IntMap.intersectionWith
  {-# INLINE liftI2 #-}

instance Ord k => Additive (Map k) where
  zero = Map.empty
  {-# INLINE zero #-}
  liftU2 = Map.unionWith
  {-# INLINE liftU2 #-}
  liftI2 = Map.intersectionWith
  {-# INLINE liftI2 #-}

instance (Eq k, Hashable k) => Additive (HashMap k) where
  zero = HashMap.empty
  {-# INLINE zero #-}
  liftU2 = HashMap.unionWith
  {-# INLINE liftU2 #-}
  liftI2 = HashMap.intersectionWith
  {-# INLINE liftI2 #-}

instance Additive ((->) b) where
  zero   = const 0
  {-# INLINE zero #-}
  liftU2 = liftA2
  {-# INLINE liftU2 #-}
  liftI2 = liftA2
  {-# INLINE liftI2 #-}

instance Additive Complex where
  zero = 0 :+ 0
  {-# INLINE zero #-}
  liftU2 f (a :+ b) (c :+ d) = f a c :+ f b d
  {-# INLINE liftU2 #-}
  liftI2 f (a :+ b) (c :+ d) = f a c :+ f b d
  {-# INLINE liftI2 #-}

instance Additive Identity where
  zero = Identity 0
  {-# INLINE zero #-}
  liftU2 = liftA2
  {-# INLINE liftU2 #-}
  liftI2 = liftA2
  {-# INLINE liftI2 #-}

-- | Compute the negation of a vector
--
-- >>> negated (V2 2 4)
-- V2 (-2) (-4)
negated :: (Functor f, Num a) => f a -> f a
negated = fmap negate
{-# INLINE negated #-}

-- | Sum over multiple vectors
--
-- >>> sumV [V2 1 1, V2 3 4]
-- V2 4 5
sumV :: (Foldable f, Additive v, Num a) => f (v a) -> v a
sumV = Foldable.foldl' (^+^) zero
{-# INLINE sumV #-}

-- | Compute the left scalar product
--
-- >>> 2 *^ V2 3 4
-- V2 6 8
(*^) :: (Functor f, Num a) => a -> f a -> f a
(*^) a = fmap (a*)
{-# INLINE (*^) #-}

-- | Compute the right scalar product
--
-- >>> V2 3 4 ^* 2
-- V2 6 8
(^*) :: (Functor f, Num a) => f a -> a -> f a
f ^* a = fmap (*a) f
{-# INLINE (^*) #-}

-- | Compute division by a scalar on the right.
(^/) :: (Functor f, Fractional a) => f a -> a -> f a
f ^/ a = fmap (/a) f
{-# INLINE (^/) #-}

-- | Produce a default basis for a vector space. If the dimensionality
-- of the vector space is not statically known, see 'basisFor'.
basis :: (Additive t, Traversable t, Num a) => [t a]
basis = basisFor (zero :: Additive v => v Int)

-- | Produce a default basis for a vector space from which the
-- argument is drawn.
basisFor :: (Traversable t, Num a) => t b -> [t a]
basisFor = \t ->
   ifoldMapOf traversed ?? t $ \i _ ->
     return                  $
       iover  traversed ?? t $ \j _ ->
         if i == j then 1 else 0
{-# INLINABLE basisFor #-}

-- | Produce a diagonal (scale) matrix from a vector.
--
-- >>> scaled (V2 2 3)
-- V2 (V2 2 0) (V2 0 3)
scaled :: (Traversable t, Num a) => t a -> t (t a)
scaled = \t -> iter t (\i x -> iter t (\j _ -> if i == j then x else 0))
  where
  iter :: Traversable t => t a -> (Int -> a -> b) -> t b
  iter x f = iover traversed f x
{-# INLINE scaled #-}

-- | Create a unit vector.
--
-- >>> unit _x :: V2 Int
-- V2 1 0
unit :: (Additive t, Num a) => ASetter' (t a) a -> t a
unit l = set' l 1 zero

-- | Outer (tensor) product of two vectors
outer :: (Functor f, Functor g, Num a) => f a -> g a -> f (g a)
outer a b = fmap (\x->fmap (*x) b) a

max_stars_repo_name: ekmett/linear

max_issues_repo_name: ekmett/linear

max_forks_repo_name: ekmett/linear


text:
{-# LANGUAGE BangPatterns     #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE Strict           #-}
module STC.InitialDistribution where

import           Array.UnboxedArray       as AU
import           Data.Array.Repa          as R
import           Data.Complex
import           Data.List                as L
import           Data.Vector.Generic      as VG
import           Data.Vector.Storable     as VS
import           DFT.Plan
import           FourierPinwheel
import           Pinwheel.FourierSeries2D
import           STC.Convolution
import           STC.DFTArray
import           STC.Point
import           STC.Utils
import           Utils.Distribution
import           Utils.List
import           Utils.Parallel
import Filter.Utils

{-# INLINE computeInitialDistribution #-}
computeInitialDistribution ::
     Int -> Int -> [Double] -> [Double] -> Double -> [Point] -> DFTArray
computeInitialDistribution rows cols phiFreqs rhoFreqs halfLogPeriod points =
  let (!minR, !maxR) = computeRange rows
      (!minC, !maxC) = computeRange cols
  in if L.any (\(Point _ _ _ s) -> s <= 0) points
       then error $ "computeInitialDistribution: initial scale <= 0."
       else DFTArray rows cols phiFreqs rhoFreqs .
            parMap
              rdeepseq
              (\(!rFreq, !thetaFreq) ->
                 VG.convert .
                 toUnboxedVector .
                 AU.accum (+) 0 ((minC, minR), (maxC, maxR)) .
                 L.map
                   (\(Point x y theta scale) ->
                      ( (round x, round y)
                      , cis $
                        (log scale) * (-rFreq) - (thetaFreq * theta * pi / 180))) $
                 points) $
            (,) <$> rhoFreqs <*> phiFreqs

{-# INLINE computeInitialDistribution' #-}
computeInitialDistribution' ::
     Int -> Int -> [Double] -> [Double] -> Double -> [Point] -> DFTArray
computeInitialDistribution' rows cols phiFreqs rhoFreqs  halfLogPeriod points =
  let (!minR, !maxR) = computeRange rows
      (!minC, !maxC) = computeRange cols
  in if L.any (\(Point _ _ _ s) -> s <= 0) points
       then error $ "computeInitialDistribution: initial scale <= 0."
       else DFTArray rows cols phiFreqs rhoFreqs .
            parMap
              rdeepseq
              (\(!thetaFreq) ->
                 VG.convert .
                 toUnboxedVector .
                 AU.accum (+) 0 ((minC, minR), (maxC, maxR)) .
                 L.map
                   (\(Point x y theta _) ->
                      ((round x, round y), cis $ -(thetaFreq * theta * pi / 180))) $
                 points) $
            phiFreqs

{-# INLINE computeInitialDistributionPowerMethod #-}
computeInitialDistributionPowerMethod ::
     Int -> Int -> [Double] -> [Double] -> Double -> [Point] -> DFTArray
computeInitialDistributionPowerMethod rows cols thetaFreqs rFreqs halfLogPeriod points =
  let (!minR, !maxR) = computeRange rows
      (!minC, !maxC) = computeRange cols
      !zeroVec = VG.replicate (rows * cols) 0
  in if L.any (\(Point _ _ _ s) -> s <= 0) points
       then error $ "computeInitialDistributionPowerMethod: initial scale <= 0."
       else DFTArray rows cols thetaFreqs rFreqs .
            L.map
              (\(!rFreq, !thetaFreq) ->
                 if rFreq == 0 && thetaFreq == 0
                   then VG.convert .
                        toUnboxedVector .
                        AU.accum (+) 0 ((minC, minR), (maxC, maxR)) .
                        L.map (\(Point x y theta scale) -> ((round x, round y), 1)) $
                        points
                   else zeroVec) $
            (,) <$> rFreqs <*> thetaFreqs

{-# INLINE computeInitialDistributionPowerMethod' #-}
computeInitialDistributionPowerMethod' ::
     Int -> Int -> [Double] -> [Double] -> [Point] -> DFTArray
computeInitialDistributionPowerMethod' rows cols phiFreqs rhoFreqs points =
  let (!minR, !maxR) = computeRange rows
      (!minC, !maxC) = computeRange cols
      !zeroVec = VG.replicate (rows * cols) 0
  in if L.any (\(Point _ _ _ s) -> s <= 0) points
       then error $ "computeInitialDistributionPowerMethod: initial scale <= 0."
       else DFTArray rows cols phiFreqs rhoFreqs .
            L.map
              (\(!thetaFreq) ->
                 if thetaFreq == 0
                   then VG.convert .
                        toUnboxedVector .
                        AU.accum (+) 0 ((minC, minR), (maxC, maxR)) .
                        L.map (\(Point x y theta scale) -> ((round x, round y), 1)) $
                        points
                   else zeroVec) $
            phiFreqs

{-# INLINE computeInitialDistributionPowerMethodSparse #-}
computeInitialDistributionPowerMethodSparse ::
     [Double] -> [Double] -> [Point] -> [R.Array U DIM2 (Complex Double)]
computeInitialDistributionPowerMethodSparse thetaFreqs rFreqs points =
  let !numThetaFreq = L.length thetaFreqs
      !numRFreq = L.length rFreqs
      !thetaCenter = div numThetaFreq 2
      !rCenter = div numRFreq 2
      !initArr =
        computeUnboxedS . fromFunction (Z :. numRFreq :. numThetaFreq) $ \(Z :. i :. j) ->
          if i == rCenter && j == thetaCenter
            then 1
            else 0
  in L.replicate (L.length points) initArr

{-# INLINE computeInitialDistributionPowerMethodSparse' #-}
computeInitialDistributionPowerMethodSparse' ::
     [Double] -> [Double] -> [Point] -> [R.Array U DIM1 (Complex Double)]
computeInitialDistributionPowerMethodSparse' thetaFreqs rFreqs points =
  let !numThetaFreq = L.length thetaFreqs
      !thetaCenter = div numThetaFreq 2
      !initArr =
        computeUnboxedS . fromFunction (Z :. numThetaFreq) $ \(Z :. j) ->
          if j == thetaCenter
            then 1
            else 0
      -- !initArr =
      --   computeUnboxedS .
      --   R.traverse (fromListUnboxed (Z :. numThetaFreq) thetaFreqs) id $ \f (Z :. thetaFreq) ->
      --     cis $ -(fromIntegral thetaFreq * 36 * pi / 180)
  in L.replicate (L.length points) initArr

computeInitialDistributionFull' ::
     Int -> Double -> Int -> Int -> [Point] -> DFTArray
computeInitialDistributionFull' numR2Freq period phiFreq rhoFreq points =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
      phiFreqs = L.map fromIntegral [-phiFreq .. phiFreq]
      rhoFreqs = L.map fromIntegral [-rhoFreq .. rhoFreq]
  in DFTArray
       numR2Freq
       numR2Freq
       phiFreqs
       rhoFreqs
       [ VG.fromList
         [ L.foldl'
           (\b (Point x y theta _) ->
              b +
              cis
                (-(angularFreq * theta * pi / 180 +
                   (freqX * x + freqY * y) * 2 * pi / period)))
           0
           points
         | freqY <- r2Freqs
         , freqX <- r2Freqs
         ]
       | angularFreq <- phiFreqs
       ]

computeInitialDistributionPowerMethodPinwheelBasis' ::
     Int -> Double -> Double -> Int -> Int -> [Point] -> DFTArray
computeInitialDistributionPowerMethodPinwheelBasis' numR2Freq sigma period phiFreq rhoFreq points =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
      phiFreqs = L.map fromIntegral [-phiFreq .. phiFreq]
      rhoFreqs = L.map fromIntegral [-rhoFreq .. rhoFreq]
      zeroVec = VG.replicate (numR2Freq ^ 2) 0
      envelope =
        VG.convert .
        toUnboxed . computeS . fromFunction (Z :. numR2Freq :. numR2Freq) $ \(Z :. i :. j) ->
          (exp $
           (fromIntegral $ (i - div numR2Freq 2) ^ 2 + (j - div numR2Freq 2) ^ 2) *
           (sigma ^ 2) /
           (-2)) *
          (sigma ^ 2) /
          (2 * pi) :+
          0
  in DFTArray
       numR2Freq
       numR2Freq
       phiFreqs
       rhoFreqs
       [ if angularFreq == 0
         then VG.zipWith (*) envelope . VG.fromList $
              [ L.foldl'
                (\b (Point x y _ _) ->
                   b + cis (-(freqX * x + freqY * y) * 2 * pi / period))
                0
                points
              | freqY <- r2Freqs
              , freqX <- r2Freqs
              ]
         else zeroVec
       | angularFreq <- phiFreqs
       ]

computeInitialDistributionPowerMethodPinwheelBasis ::
     Int -> Double -> Double -> Int -> Int -> [Point] -> DFTArray
computeInitialDistributionPowerMethodPinwheelBasis numR2Freq sigma period phiFreq rhoFreq points =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
      phiFreqs = L.map fromIntegral [-phiFreq .. phiFreq]
      rhoFreqs = L.map fromIntegral [-rhoFreq .. rhoFreq]
      zeroVec = VG.replicate (numR2Freq ^ 2) 0
      envelope =
        VG.convert .
        toUnboxed . computeS . fromFunction (Z :. numR2Freq :. numR2Freq) $ \(Z :. i :. j) ->
          (exp $
           (fromIntegral $ (i - div numR2Freq 2) ^ 2 + (j - div numR2Freq 2) ^ 2) *
           (sigma ^ 2) /
           (-2)) *
          (sigma ^ 2) /
          (2 * pi) :+
          0
      -- envelope =
      --   VG.convert . toUnboxed . computeUnboxedS $
      --   analyticalFourierCoefficients2 numR2Freq 1 0 0 sigma period (period * sqrt 2)
  in DFTArray
       numR2Freq
       numR2Freq
       phiFreqs
       rhoFreqs
       [ if angularFreq == 0 && radialFreq == 0
         then VG.zipWith (*) envelope .
              VG.fromList $
              [ L.foldl'
                (\b (Point x y _ _) ->
                   b + cis (-(freqX * x + freqY * y) * 2 * pi / period))
                0
                points
              | freqY <- r2Freqs
              , freqX <- r2Freqs
              ]
         else zeroVec
       | radialFreq <- rhoFreqs
       , angularFreq <- phiFreqs
       ]

computeInitialDistributionFourierPinwheel ::
     (VG.Vector vector (Complex Double), NFData (vector (Complex Double)))
  => Int
  -> Double
  -> Double
  -> Int
  -> Int
  -> Int
  -> Int
  -> [Point]
  -> FPArray (vector (Complex Double))
computeInitialDistributionFourierPinwheel numR2Freqs period periodEnv phiFreq rhoFreq thetaFreq rFreq points =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freqs
      thetaFreqs = L.map fromIntegral [-thetaFreq .. thetaFreq]
      rFreqs = L.map fromIntegral [-rFreq .. rFreq]
   in FPArray
        numR2Freqs
        numR2Freqs
        (2 * rFreq + 1)
        (2 * thetaFreq + 1)
        (2 * rhoFreq + 1)
        (2 * phiFreq + 1) .
      L.map
        (\radialFreq ->
           VG.concat .
           parMap
             rdeepseq
             (\angularFreq ->
                VG.convert . toUnboxed . computeS . makeFilter2D . fromListUnboxed (Z :. numR2Freqs :. numR2Freqs) $
                   [ L.foldl'
                     (\b (Point x y theta scale) ->
                        b +
                        cis
                          (-(angularFreq * theta * pi / 180 +
                             (freqX * x + freqY * y) * 2 * pi / period +
                             2 * pi / log periodEnv * log scale * radialFreq)))
                     0
                     points
                   | freqY <- r2Freqs
                   , freqX <- r2Freqs
                   ]) $
           thetaFreqs) $
      rFreqs

computeInitialDistributionPowerMethodFourierPinwheel ::
     Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> VS.Vector (Complex Double)
  -> FPArray (VS.Vector (Complex Double))
computeInitialDistributionPowerMethodFourierPinwheel numR2Freqs phiFreq rhoFreq thetaFreq rFreq bias =
  let thetaFreqs = L.map fromIntegral [-thetaFreq .. thetaFreq]
      rFreqs = L.map fromIntegral [-rFreq .. rFreq]
      numThetaFreq = 2 * thetaFreq + 1
      zeroVec = VG.replicate (numThetaFreq * numR2Freqs ^ 2) 0
   in FPArray
        numR2Freqs
        numR2Freqs
        (2 * rFreq + 1)
        numThetaFreq
        (2 * rhoFreq + 1)
        (2 * phiFreq + 1) .
      L.map
        (\radialFreq ->
           if radialFreq == 0
             then bias
             else zeroVec) $
      rFreqs

computeInitialDistributionPowerMethodFourierPinwheelFull ::
     Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> VS.Vector (Complex Double)
  -> FPArray (VS.Vector (Complex Double))
computeInitialDistributionPowerMethodFourierPinwheelFull numR2Freqs phiFreq rhoFreq thetaFreq rFreq bias =
  let numThetaFreq = 2 * thetaFreq + 1
      numRFreq = 2 * rFreq + 1
   in FPArray
        numR2Freqs
        numR2Freqs
        (2 * rFreq + 1)
        numThetaFreq
        (2 * rhoFreq + 1)
        (2 * phiFreq + 1) .
      parMap
        rdeepseq
        (\radialFreq ->
           VG.convert .
           toUnboxed .
           computeS .
           R.slice
             (fromUnboxed
                (Z :. numRFreq :. numThetaFreq :. numR2Freqs :. numR2Freqs) .
              VG.convert $
              bias) $
           (Z :. radialFreq :. All :. All :. All)) $
      [0 .. 2 * rFreq]

ext: hs

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE AllowAmbiguousTypes                      #-}
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE MultiParamTypeClasses                    #-}
{-# LANGUAGE PatternSynonyms                          #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module Backprop.Learn.Model.Neural (
    -- * Feed-forward
    FCp, fc, fca
  , fcWeights, fcBias
    -- * Recurrent
  , fcr, fcra
  , FCRp, fcrBias, fcrInputWeights, fcrStateWeights
  ) where

import           Backprop.Learn.Model.Combinator
import           Backprop.Learn.Model.Regression
import           Backprop.Learn.Model.State
import           Backprop.Learn.Model.Types
import           Data.Tuple
import           GHC.TypeNats
import           Lens.Micro
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import qualified Numeric.LinearAlgebra.Static          as H

-- | Parameters for fully connected feed-forward layer with bias.
type FCp = LRp

fcWeights :: Lens (FCp i o) (FCp i' o) (L o i) (L o i')
fcWeights = lrBeta

fcBias :: forall i o. Lens' (FCp i o) (R o)
fcBias = lrAlpha

-- | Fully connected feed-forward layer with bias.  Parameterized by its
-- initialization distribution.
--
-- Note that this has no activation function; to use as a model with
-- activation function, chain it with an activation function using 'RMap',
-- ':.~', etc.; see 'FCA' for a convenient type synonym and constructor.
--
-- Without any activation function, this is essentially a multivariate
-- linear regression.
--
-- With the logistic function as an activation function, this is
-- essentially multivariate logistic regression. (See 'logReg')
fc  :: (KnownNat i, KnownNat o)
    => Model ('Just (FCp i o)) 'Nothing (R i) (R o)
fc = linReg

-- | Convenient synonym for an 'fC' post-composed with a simple
-- parameterless activation function.
fca :: (KnownNat i, KnownNat o)
    => (forall z. Reifies z W => BVar z (R o) -> BVar z (R o))
    -> Model ('Just (FCp i o)) 'Nothing (R i) (R o)
fca f = funcD f <~ linReg

-- | Fully connected recurrent layer with bias.
fcr :: (KnownNat i, KnownNat o, KnownNat s)
    => (forall z. Reifies z W => BVar z (R o) -> BVar z (R s))      -- ^ store
    -> Model ('Just (FCRp s i o)) ('Just (R s)) (R i) (R o)
fcr s = recurrent H.split (H.#) s fc

-- | Convenient synonym for an 'fcr' post-composed with a simple
-- parameterless activation function.
fcra
    :: (KnownNat i, KnownNat o, KnownNat s)
    => (forall z. Reifies z W => BVar z (R o) -> BVar z (R o))
    -> (forall z. Reifies z W => BVar z (R o) -> BVar z (R s))      -- ^ store
    -> Model ('Just (FCRp s i o)) ('Just (R s)) (R i) (R o)
fcra f s = funcD f <~ recurrent H.split (H.#) s fc

-- | Parameter for fully connected recurrent layer.
type FCRp s i o = FCp (i + s) o

lensIso :: (s -> (a, x)) -> ((b, x) -> t) -> Lens s t a b
lensIso f g h x = g <$> _1 h (f x)

fcrInputWeights
    :: (KnownNat s, KnownNat i, KnownNat i', KnownNat o)
    => Lens (FCRp s i o) (FCRp s i' o) (L o i) (L o i')
fcrInputWeights = fcWeights
                . lensIso H.splitCols (uncurry (H.|||))

fcrStateWeights
    :: (KnownNat s, KnownNat s', KnownNat i, KnownNat o)
    => Lens (FCRp s i o) (FCRp s' i o) (L o s) (L o s')
fcrStateWeights = fcWeights
                . lensIso (swap . H.splitCols) (uncurry (H.|||) . swap)

fcrBias :: forall s i o. Lens' (FCRp s i o) (R o)
fcrBias = fcBias @(i + s) @o

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
module Edges.Prelude
(
  module Exports,
  UnboxedVector,
  modifyTVar',
  forMInAscendingRange_,
  forMInDescendingRange_,
  showText,
  -- * Optics
  Lens,
  Lens',
  Prism,
  Prism',
  lens,
  prism,
  _Left,
  _Just,
)
where

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

-- contravariant
-------------------------
import Data.Functor.Contravariant as Exports
import Data.Functor.Contravariant.Divisible as Exports

-- pointed
-------------------------
import Data.Pointed as Exports
import Data.Copointed as Exports

-- profunctors
-------------------------
import Data.Profunctor.Unsafe as Exports
import Data.Profunctor.Choice as Exports
import Data.Profunctor.Strong as Exports

-- semigroupoids
-------------------------
import Data.Semigroupoid as Exports

-- cereal
-------------------------
import Data.Serialize as Exports (Serialize)

-- deferred-folds
-------------------------
import DeferredFolds.Unfoldl as Exports (Unfoldl(..))
import DeferredFolds.UnfoldlM as Exports (UnfoldlM(..))

-- foldl
-------------------------
import Control.Foldl as Exports (Fold(..), FoldM(..))

-- hashable
-------------------------
import Data.Hashable as Exports (Hashable)

-- primitive
-------------------------
import Data.Primitive as Exports
import Control.Monad.Primitive as Exports

-- primitive-extras
-------------------------
import PrimitiveExtras.PrimMultiArray as Exports (PrimMultiArray)
import PrimitiveExtras.TVarArray as Exports (TVarArray)

-- unordered-containers
-------------------------
import Data.HashMap.Strict as Exports (HashMap)

-- vector
-------------------------
import Data.Vector as Exports (Vector)

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

-- deepseq
-------------------------
import Control.DeepSeq as Exports

-- data-dword
-------------------------
import Data.DoubleWord as Exports

-- cereal-vector
-------------------------
import Data.Vector.Serialize ()

-- cereal-data-dword
-------------------------
import Cereal.DataDWord.Serialize ()

import qualified Data.Vector.Unboxed as UnboxedVector
import qualified Data.Vector.Unboxed.Deriving as UnboxedVectorDeriving

type UnboxedVector = UnboxedVector.Vector

{-# INLINE modifyTVar' #-}
modifyTVar' :: TVar a -> (a -> a) -> STM ()
modifyTVar' var f = do
    x <- readTVar var
    writeTVar var $! f x

{-# INLINE forMInAscendingRange_ #-}
forMInAscendingRange_ :: Applicative m => Int -> Int -> (Int -> m a) -> m ()
forMInAscendingRange_ !startN !endN f =
  ($ startN) $ fix $ \loop !n -> if n < endN then f n *> loop (succ n) else pure ()

{-# INLINE forMInDescendingRange_ #-}
forMInDescendingRange_ :: Applicative m => Int -> Int -> (Int -> m a) -> m ()
forMInDescendingRange_ !startN !endN f =
  ($ pred startN) $ fix $ \loop !n -> if n >= endN then f n *> loop (pred n) else pure ()

showText = fromString . show

-- * Optics
-------------------------

type Lens s t a b = forall f. Functor f => (a -> f b) -> s -> f t

type Lens' s a = Lens s s a a

type Prism s t a b = forall p f. (Choice p, Applicative f) => p a (f b) -> p s (f t)

type Prism' s a = Prism s s a a

{-# INLINE lens #-}
lens :: (s -> a) -> (s -> b -> t) -> Lens s t a b
lens sa sbt afb s = sbt s <$> afb (sa s)

{-# INLINE prism #-}
prism :: (b -> t) -> (s -> Either t a) -> Prism s t a b
prism bt seta = dimap seta (either pure (fmap bt)) . right'

{-# INLINE _Left #-}
_Left :: Prism' (Either a b) a
_Left = prism Left (either Right (Left . Right))

{-# INLINE _Just #-}
_Just :: Prism' (Maybe a) a
_Just = prism Just (maybe (Left Nothing) Right)

-- * Instances for unboxed vectors
-------------------------

UnboxedVectorDeriving.derivingUnbox "Word96"
  [t| Word96 -> (Word32, Word64) |]
  [| \ (Word96 a b) -> (a, b) |]
  [| \ (a, b) -> Word96 a b |]

UnboxedVectorDeriving.derivingUnbox "Word128"
  [t| Word128 -> (Word64, Word64) |]
  [| \ (Word128 a b) -> (a, b) |]
  [| \ (a, b) -> Word128 a b |]

UnboxedVectorDeriving.derivingUnbox "Word256"
  [t| Word256 -> (Word128, Word128) |]
  [| \ (Word256 a b) -> (a, b) |]
  [| \ (a, b) -> Word256 a b |]

max_stars_repo_name: metrix-ai/edges

max_issues_repo_name: metrix-ai/edges

max_forks_repo_name: metrix-ai/edges


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

-- |
-- Module      :  Test.Factorization
-- Copyright   :  (c) 2016-2020 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Test.Codegen (
    codegenTests,
    ditCodegenTests,
    difCodegenTests,
    splitRadixCodegenTests,
    improvedSplitRadixCodegenTests,
    searchCodegenTests,
    modCodegenTests
  ) where

import Control.Monad (mzero,
                      replicateM)
import Data.Complex
import Data.Modular
import Data.Proxy (Proxy(..))
import Data.Typeable (Typeable)
import qualified Data.Vector.Storable as VS
import Data.Word (Word32)
import GHC.TypeLits (KnownNat,
                     natVal)
import Test.Hspec
import Test.Hspec.Core.Spec (Result)
import Test.QuickCheck ((===),
                        Arbitrary(..),
                        Gen,
                        Property,
                        counterexample,
                        forAll)
import Test.QuickCheck.Modifiers (Positive(..))

import Spiral (Config(..))
import Spiral.Driver
import Spiral.Exp
import Spiral.FFT.CooleyTukey
import Spiral.Monad
import Spiral.RootOfUnity
import Spiral.SPL
import Spiral.Search
import Spiral.Search.FFTBreakdowns
import Spiral.Search.OpCount

import qualified Test.FFTW as FFTW
import Test.Gen
import Test.Instances ()

codegenTests :: Config
             -> [Int]
             -> Spec
codegenTests conf sizes = do
    ditCodegenTests conf sizes
    difCodegenTests conf sizes
    splitRadixCodegenTests conf sizes
    improvedSplitRadixCodegenTests conf sizes
    modCodegenTests conf

ditCodegenTests :: Config -> [Int] -> Spec
ditCodegenTests conf sizes =
    describe "Generated DIT" $
    mkCodegenTests conf "DIT" (return . Re . dit) sizes

difCodegenTests :: Config -> [Int] -> Spec
difCodegenTests conf sizes =
    describe "Generated DIF" $
    mkCodegenTests conf "DIF" (return . Re . dif) sizes

splitRadixCodegenTests :: Config -> [Int] -> Spec
splitRadixCodegenTests conf sizes =
    describe "Generated splitRadix DFT" $
    mkCodegenTests conf "SplitRadix" (runSearch () f . Re . DFT) sizes
  where
    f :: (Typeable a, Typed a, MonadSpiral m)
      => SPL (Exp a)
      -> S s m (SPL (Exp a))
    f (F n w) = splitRadixBreakdown n w
    f _       = mzero

improvedSplitRadixCodegenTests :: Config -> [Int] -> Spec
improvedSplitRadixCodegenTests conf sizes =
    describe "Generated improved split radix DFT" $
    mkCodegenTests conf "ImprovedSplitRadix" (runSearch () f . Re . DFT) sizes
  where
    f :: (Typeable a, Typed a, MonadSpiral m)
      => SPL (Exp a)
      -> S s m (SPL (Exp a))
    f (F n w) = improvedSplitRadixBreakdown n w
    f _       = mzero

searchCodegenTests :: Config -> [Int] -> Spec
searchCodegenTests conf sizes =
    describe "Generated opcount-optimized DFT" $
    mkCodegenTests conf "OpcountSearch" (searchOpCount . Re . DFT) sizes

mkCodegenTests :: Config
               -> String
               -> (Int -> Spiral (SPL (Exp Double)))
               -> [Int]
               -> Spec
mkCodegenTests conf desc f = mapM_ dftTest
  where
    dftTest :: Int -> Spec
    dftTest n = it (desc ++ "(" ++ show n ++ ")") $ \(k :: Property -> IO Result) -> do
        Re e <- runSpiralWith mempty $ f n
        withComplexTransform conf ("dft" ++ show n) e $ \dft ->
          k $ forAll (uniformVectorsOfSize n) $ \v -> epsDiff (dft v) (FFTW.fft n v)

-- | Generate vectors of a given size with uniformly distributed elements.
uniformVectorsOfSize :: (Arbitrary (Uniform01 a), VS.Storable a)
                     => Int
                     -> Gen (VS.Vector a)
uniformVectorsOfSize n = do
    xs <- replicateM n arbitrary
    return $ VS.fromList (map unUniform01 xs)

-- | A uniform value i the range [0,1]
newtype Uniform01 a = Uniform01 { unUniform01 :: a }
  deriving (Eq, Ord, Show)

instance Arbitrary (Uniform01 a) => Arbitrary (Uniform01 (Complex a)) where
    arbitrary = do
        Uniform01 r <- arbitrary
        Uniform01 i <- arbitrary
        return $ Uniform01 $ r :+ i

instance Arbitrary (Uniform01 Double) where
    -- This code take from the dsp Haskell library, authored by Matt Donadio.
    -- 53 bits in [0,1], i.e., 64-bit IEEE 754 in [0,1]
    -- 67108864 = 2^26
    -- 9007199254740991 = 2^53 - 1
    arbitrary = do
        u1 :: Word32 <- arbitrary
        u2 :: Word32 <- arbitrary
        -- 27 bits
        let a = fromIntegral u1 / 32.0
        -- 26 bits
        let b = fromIntegral u2 / 64.0
        return $ Uniform01 $ (a * 67108864.0 + b) / 9007199254740991.0

-- | Return 'True' if the maximum pairwise difference between two vectors is
-- less than epsilon.
epsDiff :: forall a . (Show a, RealFloat a, VS.Storable a)
        => VS.Vector (Complex a)
        -> VS.Vector (Complex a)
        -> Property
epsDiff v1 v2 =
    counterexample ("Max delta: " ++ show maxDelta) $
    maxDelta < eps
  where
    maxDelta :: a
    maxDelta = VS.maximum $ VS.zipWith (\x y -> magnitude (x - y)) v1 v2

    eps :: a
    eps = 1e-12

modCodegenTests :: Config
                -> Spec
modCodegenTests conf =
    describe "Generated modular DFT" $ do
    describe "ℤ/17" $
        mkModCodegenTests (Proxy :: Proxy 17) conf "DIT" (return . dit) [2^n | n <- [1..4::Int]]
    describe "ℤ/2013265921" $
        mkModCodegenTests (Proxy :: Proxy 2013265921) conf "DIT" (return . dit) [2^n | n <- [1..12::Int]]

mkModCodegenTests :: forall p . KnownNat p
                  => Proxy p
                  -> Config
                  -> String
                  -> (Int -> Spiral (SPL (Exp (ℤ/p))))
                  -> [Int]
                  -> Spec
mkModCodegenTests proxy conf desc f = mapM_ dftTest
  where
    p :: Integer
    p = natVal proxy

    dftTest :: Int -> Spec
    dftTest n = it (desc ++ "(" ++ show n ++ ")") $ \(k :: Property -> IO Result) -> do
        e <- runSpiralWith mempty $ f n
        withModularTransform conf ("moddft_" ++ show p ++ "_" ++ show n) e $ \moddft ->
          k $ forAll basis $ \i -> moddft (e_i i) === res i
      where
        basis :: Gen Int
        basis = do
            Positive x <- arbitrary
            return $ x `mod` n

        e_i :: Int -> VS.Vector (ℤ/p)
        e_i i = VS.generate n $ \j -> if i == j then 1 else 0

        res :: Int -> VS.Vector (ℤ/p)
        res i = VS.generate n $ \j -> a^j
          where
            a :: ℤ/p
            a = w^i

        w :: ℤ/p
        w = omega n

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
-- nonlinear least-squares fitting

import Spinell.Fitting
import Spinell.Graphing
import Numeric.LinearAlgebra

model [a,lambda,b] t = a * exp (-lambda * t)+b --the function we will be fitting

jac [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]

xs = [0..100]

sigma = replicate (length xs) 0.1 -- will create an array of repeating '0.1': [0.1,0.1...]

noise = toList $ scalar 0.1 * (randomVector 0 Gaussian (length xs)) --a list of random small values, that will play the role of the noise in the data

ys = map (model [5,0.1,1]) xs --we plainly apply the model to the data

ysNoisy = [ y + err | (y,err) <- zip ys noise] -- and then add the noise, just to make it look more realistic

res = fit xs ys sigma model [1,1,1] defFitOpt{jacob = ManualJacob jac}

main :: IO()
main = plotFit xs ysNoisy sigma model (getParams res) defGraphOpt

max_stars_repo_name: Magalame/Spinell

max_issues_repo_name: Magalame/Spinell

max_forks_repo_name: Magalame/Spinell


text:
{-# LANGUAGE FlexibleContexts #-}
module CW (parseOpts, run) where

import           Control.Monad
import           Control.Monad.Error.Class
import           Data.Char
import           Data.Complex
import           Data.List
import qualified Filter as F
import           IQ
import           Options.Applicative
import           Pipes
import qualified Pipes.Prelude as P
import           Render
import           System.IO
import           ZeroStuff

data CWOptions = CWOptions
  { sampleRate :: Int
  , wpm :: Int
  , text :: String
  , outputFile :: String }

parseOpts :: Parser CWOptions
parseOpts = CWOptions
         <$> option auto
             ( long "sample-rate"
            <> short 's'
            <> metavar "SAMPLE_RATE"
            <> help "Sample rate of the produced IQ file" )
         <*> option auto
             ( long "wpm"
            <> short 'w'
            <> value 25
            <> showDefault
            <> metavar "WPM_RATE"
            <> help "Rate of keying in Words per Minute" )
         <*> strOption
             ( long "in-text"
            <> short 't'
            <> metavar "INPUT_TEXT"
            <> help "Input text message" )
         <*> strOption
             ( long "out-file"
            <> short 'o'
            <> metavar "OUTPUT_FILE"
            <> help "Output IQ file" )

data CW = Dah
        | Dit

cwMap :: MonadError String m => Char -> m [CW]
cwMap c
  | c == 'a' = pure [ Dah, Dit ]
  | c == 'b' = pure [ Dah, Dit, Dit, Dit ]
  | c == 'c' = pure [ Dah, Dit, Dah, Dit ]
  | c == 'd' = pure [ Dah, Dit, Dit ]
  | c == 'e' = pure [ Dit ]
  | c == 'f' = pure [ Dit, Dit, Dah, Dit ]
  | c == 'g' = pure [ Dah, Dah, Dit ]
  | c == 'h' = pure [ Dit, Dit, Dit, Dit ]
  | c == 'i' = pure [ Dit, Dit ]
  | c == 'j' = pure [ Dit, Dah, Dah, Dah ]
  | c == 'k' = pure [ Dah, Dit, Dah ]
  | c == 'l' = pure [ Dit, Dah, Dit, Dit ]
  | c == 'm' = pure [ Dah, Dah ]
  | c == 'n' = pure [ Dah, Dit ]
  | c == 'o' = pure [ Dah, Dah, Dah ]
  | c == 'p' = pure [ Dit, Dah, Dah, Dit ]
  | c == 'q' = pure [ Dah, Dah, Dit, Dah ]
  | c == 'r' = pure [ Dit, Dah, Dit ]
  | c == 's' = pure [ Dit, Dit, Dit ]
  | c == 't' = pure [ Dah ]
  | c == 'u' = pure [ Dit, Dit, Dah ]
  | c == 'v' = pure [ Dit, Dit, Dit, Dah ]
  | c == 'w' = pure [ Dit, Dah, Dah ]
  | c == 'x' = pure [ Dah, Dit, Dit, Dah ]
  | c == 'y' = pure [ Dah, Dit, Dah, Dah ]
  | c == 'z' = pure [ Dah, Dah, Dit, Dit ]
  | c == '1' = pure [ Dit, Dah, Dah, Dah, Dah ]
  | c == '2' = pure [ Dit, Dit, Dah, Dah, Dah ]
  | c == '3' = pure [ Dit, Dit, Dit, Dah, Dah ]
  | c == '4' = pure [ Dit, Dit, Dit, Dit, Dah ]
  | c == '5' = pure [ Dit, Dit, Dit, Dit, Dit ]
  | c == '6' = pure [ Dah, Dit, Dit, Dit, Dit ]
  | c == '7' = pure [ Dah, Dah, Dit, Dit, Dit ]
  | c == '8' = pure [ Dah, Dah, Dah, Dit, Dit ]
  | c == '9' = pure [ Dah, Dah, Dah, Dah, Dit ]
  | c == '0' = pure [ Dah, Dah, Dah, Dah, Dah ]
  | otherwise = throwError $ "Unknown CW character: '" ++ [c] ++ "'"

data KeyState = On | Off
data Keying = Keying
  { units :: Int
  , keyState :: KeyState }

cwKeyingMap :: CW -> Keying
cwKeyingMap Dah = Keying 3 On
cwKeyingMap Dit = Keying 1 On

keyLetter :: MonadError String m => Char -> m [Keying]
keyLetter l = intersperse (Keying 1 Off) <$> (map cwKeyingMap <$> cwMap l)

keyWord :: MonadError String m => String -> m [Keying]
keyWord word = intercalate  [Keying 3 Off] <$> sequence (keyLetter <$> word)

keyText :: String -> Either String [Keying]
keyText text = intercalate [Keying 7 Off] <$> sequence (keyWord <$> words text)

synthesizeKeying :: MonadIO m => Int -> Int -> Pipe Keying IQ m ()
synthesizeKeying sr wpm = forever $ do
  keying <- await

  let iq = case keying of
        Keying _ On -> 1.0 :+ 0.0
        Keying _ Off  -> 0.0 :+ 0.0
      n = round $ fromIntegral sr * (60.0 / (50.0 * fromIntegral wpm)) * fromIntegral (units keying)
  replicateM_ n $ yield iq

run :: CWOptions -> IO ()
run opts = case keyText $ toLower <$> text opts of
  Left err -> putStr err
  Right keyings -> do
    putStr "Synthesizing samples... "
    hFlush stdout

    let kernel = F.lowPass (sampleRate opts) $ F.FilterOptions 100 1000

    withBinaryFile (outputFile opts) WriteMode $ \f ->
      runEffect $ each keyings
           >-> synthesizeKeying (sampleRate opts) (wpm opts)
           >-> F.convolve kernel
           >-> cfileSink f

    putStrLn "Done."

max_stars_repo_name: hexagonal-sun/ayeQ

max_issues_repo_name: hexagonal-sun/ayeQ

max_forks_repo_name: hexagonal-sun/ayeQ


text:
-- ch03 example 3.8

import Statistics.Distribution
import Statistics.Distribution.Binomial

b :: BinomialDistribution
b = binomial 10 0.3

main = do
  putStrLn $ "If, at least nine of ten will recover:  "
          ++ show (1-cumulative b 8.00)
  print $ complCumulative b 8.00
  print $ probability b 9 + probability b 10

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
module Titanic where

import Numeric.LinearAlgebra                         (col, disp, Matrix, R, rand, scale, toLists, tr, Vector, (><))
import Numeric.LinearAlgebra.HMatrix                 (mul, matrix, sumElements)
import Text.Read                                     (readMaybe)

type Gender = Double
type Age = Double
type SibSp = Double
type Embarked = Double
type PClass = Double
type Survived = Double

data Row = Row Gender Age PClass Embarked

processRow :: [[Char]] -> [Double]
processRow row = [gender, age] ++ pclass ++ embarked where
    gender = getGender row
    age = getAge row
    pclass = getPClass row
    embarked = getEmbarked row

processFeatures :: [[[Char]]] -> Matrix Double
processFeatures inputs = matrix 8 $ concat [processRow x | x <- inputs]

processGender :: [[[Char]]] -> Matrix Gender
processGender inputs = do
    col [getGender x | x <- inputs]

processTargets :: [[[Char]]] -> Matrix Survived
processTargets inputs = do
    col [getTarget x | x <- inputs]

getAge :: [[Char]] -> Age
getAge row =  case (readMaybe (row !! 6)) of
    -- 80 is the maximum age in the dataset. We divide the age by it to
    -- keep it between 0.0 and 1.0
    Just x -> x / 80
    Nothing -> 0.0

getPClass :: [[Char]] -> [PClass]
getPClass row =  case (row !! 2) of
    "1" -> [1.0, 0.0, 0.0]
    "2" -> [0.0, 1.0, 0.0]
    "3" -> [0.0, 0.0, 1.0]
    "" -> [0.0, 0.0, 0.0]

getEmbarked :: [[Char]] -> [Embarked]
getEmbarked row = case (readMaybe (row !! 12)) of
    Just x -> case x of
        "C" -> [1.0, 0.0, 0.0]
        "Q" -> [0.0, 1.0, 0.0]
        "S" -> [0.0, 0.0, 1.0]
    Nothing -> [0.0, 0.0, 0.0]

getGender :: [[Char]] -> Gender
getGender row = case (row !! 5) of
    "male"   -> -1.0
    "female" -> 1.0

getTarget :: [[Char]] -> Survived
getTarget row = fromIntegral $ read $ row !! 1

max_stars_repo_name: staturecrane/haskell-ml

max_issues_repo_name: staturecrane/haskell-ml

max_forks_repo_name: staturecrane/haskell-ml


text:
-- | Dynamic Programming module.
module DP where

import Data.Map (Map)
import qualified Data.Map.Strict as M
import qualified Data.Set as S
import Numeric.LinearAlgebra.Data
       (Matrix, Vector, (!), (¿), (><), maxIndex)
import Numeric.LinearAlgebra.HMatrix as H

import MDP
import Planning

import Debug.Trace

data MDPM s a = MDPM
    { statesm :: Map s Int
    , actionsm :: Map a Int
    , transitionm :: Map a (Matrix Double)
    , rewardm :: Matrix Double
    , gammam :: Double             -- discount
    , sizeStates :: Int           -- S
    , sizeActions :: Int          -- A
    , actionsr :: Map Int a
    }

prepareMDP :: (Ord s, Ord a) => MDP s a -> MDPM s a
prepareMDP mdp =
    MDPM
    { statesm = statesM
    , actionsm = actionsM
    , transitionm = transitionM
    , rewardm = rewardM
    , gammam = gamma mdp
    , sizeStates = S.size (states mdp)
    , sizeActions = S.size (actions mdp)
    , actionsr = M.fromList (zip [0..] (S.toList (actions mdp)))
    }
  where
    statesM = M.fromList (zip (S.toList (states mdp)) [0 ..])
    actionsM = M.fromList (zip (S.toList (actions mdp)) [0 ..])
    transitionM = transitionMap
    transitionMap =
        M.fromList
            (map
                 (\a ->
                       (a, createMatrix a))
                 (S.toList (actions mdp)))
    createMatrix a =
        H.fromLists
            (map
                 (\s ->
                       map (transition mdp a s) (S.toList (states mdp)))
                 (S.toList (states mdp)))
    rewardM = createRewardMatrix
    createRewardMatrix =
        H.fromLists
            (map
                 (\s ->
                       map
                           (\a ->
                                 reward mdp a s)
                           (S.toList (actions mdp)))
                 (S.toList (states mdp)))

valueIterationM :: (Ord a) => MDPM s a -> Int -> Matrix Double -> Matrix Double
valueIterationM mdpm k v
  | k == 0 = v
  | otherwise = valueIterationM mdpm (k - 1) (calc v)
  where
    calc v' =
        rewardm mdpm +
        gammam mdpm `H.scale`
        fromColumns
            (map
                 (\a ->
                       flatten
                           (transitionm mdpm M.! a H.<>
                            (v' ¿ [actionsm mdpm M.! a])))
                 (M.keys (actionsm mdpm)))

valueIteration :: (Ord a, Ord s) => MDP s a -> Int -> s -> a
valueIteration mdp k s = actionsr mdpm M.! maxIndex (vk ! (statesm mdpm M.! s))
  where
    mdpm = prepareMDP mdp
    v0 = (sizeStates mdpm >< sizeActions mdpm) [0,0 ..]
    vk = valueIterationM mdpm k v0

max_stars_repo_name: DbIHbKA/hmdp

max_issues_repo_name: DbIHbKA/hmdp

max_forks_repo_name: DbIHbKA/hmdp


text:
{-# LANGUAGE FlexibleInstances #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Foreign.LAPACK
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Types with a Linear Algebra PACKage.
--
module Foreign.LAPACK (
    -- * LAPACK typeclass
    LAPACK(..),
    -- * Enums
    module Foreign.LAPACK.Types,
    ) where

import Control.Monad( forM_, when )
import Control.Exception( assert )
import Data.Complex( Complex )
import Foreign
import Foreign.BLAS
import Foreign.LAPACK.Types
import Foreign.LAPACK.Double
import Foreign.LAPACK.Zomplex

-- | Types with LAPACK operations.
class (BLAS3 e) => LAPACK e where
    geqrf :: Int -> Int -> Ptr e -> Int -> Ptr e -> IO ()
    gelqf :: Int -> Int -> Ptr e -> Int -> Ptr e -> IO ()

    heevr :: EigJob -> EigRange -> Uplo -> Int -> Ptr e -> Int -> Double
          -> Ptr Double -> Ptr e -> Int -> Ptr Int -> IO Int

    larfg :: Int -> Ptr e -> Ptr e -> Int -> IO e
    potrf :: Uplo -> Int -> Ptr e -> Int -> IO Int
    potrs :: Uplo -> Int -> Int -> Ptr e -> Int -> Ptr e -> Int -> IO ()
    pptrf :: Uplo -> Int -> Ptr e -> IO Int
    pptrs :: Uplo -> Int -> Int -> Ptr e -> Ptr e -> Int -> IO ()

    unmqr :: Side -> Trans -> Int -> Int -> Int -> Ptr e -> Int -> Ptr e -> Ptr e -> Int -> IO ()
    unmlq :: Side -> Trans -> Int -> Int -> Int -> Ptr e -> Int -> Ptr e -> Ptr e -> Int -> IO ()

callWithWork :: (Storable e) => (Ptr e -> Ptr LAInt -> Ptr LAInt -> IO ()) -> IO LAInt
callWithWork call =
    alloca $ \pquery ->
    with (-1) $ \pnone ->
    alloca $ \pinfo -> do
        call pquery pnone pinfo
        lwork <- (ceiling . max 1) `fmap` (peek (castPtr pquery) :: IO Double)
        allocaArray lwork $ \pwork ->
            withEnum lwork $ \plwork -> do
                call pwork plwork pinfo
                peek pinfo

callWithWorkIWork :: (Storable e)
    => (Ptr e -> Ptr LAInt -> Ptr LAInt -> Ptr LAInt -> Ptr LAInt -> IO ())
    -> IO LAInt
callWithWorkIWork call =
    alloca $ \pquery ->
    alloca $ \piquery ->
    with (-1) $ \pnone ->
    alloca $ \pinfo -> do
        call pquery pnone piquery pnone pinfo
        lwork <- (ceiling . max 1) `fmap` (peek (castPtr pquery) :: IO Double)
        liwork <- max 1 `fmap` peek piquery
        allocaArray lwork $ \pwork ->
            withEnum lwork $ \plwork ->
                allocaArray (fromEnum liwork) $ \piwork ->
                    with liwork $ \pliwork -> do
                        call pwork plwork piwork pliwork pinfo
                        peek pinfo

callWithWorkRWorkIWork :: (Storable e)
    => (Ptr e -> Ptr LAInt -> Ptr Double -> Ptr LAInt -> Ptr LAInt -> Ptr LAInt -> Ptr LAInt -> IO ())
    -> IO LAInt
callWithWorkRWorkIWork call =
    alloca $ \pquery ->
    alloca $ \prquery ->
    alloca $ \piquery ->
    with (-1) $ \pnone ->
    alloca $ \pinfo -> do
        call pquery pnone prquery pnone piquery pnone pinfo
        lwork <- (ceiling . max 1) `fmap` (peek (castPtr pquery) :: IO Double)
        lrwork <- (ceiling . max 1) `fmap` peek prquery
        liwork <- max 1 `fmap` peek piquery
        allocaArray lwork $ \pwork ->
            withEnum lwork $ \plwork ->
                allocaArray lrwork $ \prwork ->
                    withEnum lrwork $ \plrwork ->
                        allocaArray (fromEnum liwork) $ \piwork ->
                            with liwork $ \pliwork -> do
                                call pwork plwork prwork plrwork piwork pliwork pinfo
                                peek pinfo

checkInfo :: (Num info) => info -> IO ()
checkInfo info = assert (info == 0) $ return ()

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

instance LAPACK Double where
    geqrf m n pa lda ptau =
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
            checkInfo =<< callWithWork (dgeqrf pm pn pa plda ptau)

    gelqf m n pa lda ptau =
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
            checkInfo =<< callWithWork (dgelqf pm pn pa plda ptau)

    heevr jobz range uplo n pa lda abstol pw pz ldz psuppz =
        withEigJob jobz $ \pjobz ->
        withEigRange range $ \prange pvl pvu pil piu ->
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
        with abstol $ \pabstol ->
        alloca $ \pm ->
        withEnum ldz $ \pldz ->
        allocaArray nsuppz $ \psuppz' -> do
            checkInfo =<<
                (callWithWorkIWork $
                    dsyevr pjobz prange puplo pn pa plda pvl pvu pil piu
                           pabstol pm pw pz pldz psuppz')
            m <- fromEnum `fmap` peek pm
            when (psuppz /= nullPtr) $
                forM_ [ 0..2*m - 1 ] $ \i -> do
                    s <- peekElemOff psuppz' i
                    pokeElemOff psuppz i (fromEnum s)
            return m
      where
        nsuppz = if psuppz == nullPtr then 0 else 2 * max 1 n

    unmqr side trans m n k pa lda ptau pc ldc =
        withSide side $ \pside ->
        withTrans trans $ \ptrans ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        withEnum lda $ \plda ->
        withEnum ldc $ \pldc ->
            checkInfo =<< callWithWork (dormqr pside ptrans pm pn pk pa plda ptau pc pldc)

    unmlq side trans m n k pa lda ptau pc ldc =
        withSide side $ \pside ->
        withTrans trans $ \ptrans ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        withEnum lda $ \plda ->
        withEnum ldc $ \pldc ->
            checkInfo =<< callWithWork (dormlq pside ptrans pm pn pk pa plda ptau pc pldc)

    larfg n palpha px incx =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        alloca $ \ptau -> do
            dlarfg pn palpha px pincx ptau
            peek ptau

    potrf uplo n pa lda =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
        alloca $ \pinfo -> do
            dpotrf puplo pn pa plda pinfo
            info <- fromEnum `fmap` peek pinfo
            assert (info >= 0) $ return info

    potrs uplo n nrhs pa lda pb ldb =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum nrhs $ \pnrhs ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        alloca $ \pinfo -> do
            dpotrs puplo pn pnrhs pa plda pb pldb pinfo
            checkInfo =<< peek pinfo

    pptrf uplo n pa =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        alloca $ \pinfo -> do
            dpptrf puplo pn pa pinfo
            info <- fromEnum `fmap` peek pinfo
            assert (info >= 0) $ return info

    pptrs uplo n nrhs pa pb ldb =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum nrhs $ \pnrhs ->
        withEnum ldb $ \pldb ->
        alloca $ \pinfo -> do
            dpptrs puplo pn pnrhs pa pb pldb pinfo
            checkInfo =<< peek pinfo

instance LAPACK (Complex Double) where
    geqrf m n pa lda ptau =
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
            checkInfo =<< callWithWork (zgeqrf pm pn pa plda ptau)

    gelqf m n pa lda ptau =
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
            checkInfo =<< callWithWork (zgelqf pm pn pa plda ptau)

    heevr jobz range uplo n pa lda abstol pw pz ldz psuppz =
        withEigJob jobz $ \pjobz ->
        withEigRange range $ \prange pvl pvu pil piu ->
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
        with abstol $ \pabstol ->
        alloca $ \pm ->
        withEnum ldz $ \pldz ->
        allocaArray nsuppz $ \psuppz' -> do
            checkInfo =<<
                (callWithWorkRWorkIWork $
                    zheevr pjobz prange puplo pn pa plda pvl pvu pil piu
                           pabstol pm pw pz pldz psuppz')
            m <- fromEnum `fmap` peek pm
            when (psuppz /= nullPtr) $
                forM_ [ 0..2*m - 1 ] $ \i -> do
                    s <- peekElemOff psuppz' i
                    pokeElemOff psuppz i (fromEnum s)
            return m
      where
        nsuppz = if psuppz == nullPtr then 0 else 2 * max 1 n

    unmqr side trans m n k pa lda ptau pc ldc =
        withSide side $ \pside ->
        withTrans trans $ \ptrans ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        withEnum lda $ \plda ->
        withEnum ldc $ \pldc ->
            checkInfo =<< callWithWork (zunmqr pside ptrans pm pn pk pa plda ptau pc pldc)

    unmlq side trans m n k pa lda ptau pc ldc =
        withSide side $ \pside ->
        withTrans trans $ \ptrans ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        withEnum lda $ \plda ->
        withEnum ldc $ \pldc ->
            checkInfo =<< callWithWork (zunmlq pside ptrans pm pn pk pa plda ptau pc pldc)

    larfg n palpha px incx =
        withEnum n $ \pn ->
        withEnum incx $ \pincx ->
        alloca $ \ptau -> do
            zlarfg pn palpha px pincx ptau
            peek ptau

    potrf uplo n pa lda =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum lda $ \plda ->
        alloca $ \pinfo -> do
            zpotrf puplo pn pa plda pinfo
            info <- fromEnum `fmap` peek pinfo
            assert (info >= 0) $ return info

    potrs uplo n nrhs pa lda pb ldb =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum nrhs $ \pnrhs ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        alloca $ \pinfo -> do
            zpotrs puplo pn pnrhs pa plda pb pldb pinfo
            checkInfo =<< peek pinfo

    pptrf uplo n pa =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        alloca $ \pinfo -> do
            zpptrf puplo pn pa pinfo
            info <- fromEnum `fmap` peek pinfo
            assert (info >= 0) $ return info

    pptrs uplo n nrhs pa pb ldb =
        withUplo uplo $ \puplo ->
        withEnum n $ \pn ->
        withEnum nrhs $ \pnrhs ->
        withEnum ldb $ \pldb ->
        alloca $ \pinfo -> do
            zpptrs puplo pn pnrhs pa pb pldb pinfo
            checkInfo =<< peek pinfo

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: patperry/hs-linear-algebra

max_forks_repo_name: patperry/hs-linear-algebra


text:
{-# LANGUAGE ViewPatterns #-}
-- | Calculation of confidence intervals
module Statistics.ConfidenceInt (
    poissonCI
  , poissonNormalCI
  , binomialCI
  , naiveBinomialCI
    -- * References
    -- $references
  ) where

import Statistics.Distribution
import Statistics.Distribution.ChiSquared
import Statistics.Distribution.Beta
import Statistics.Types

-- | Calculate confidence intervals for Poisson-distributed value
-- using normal approximation
poissonNormalCI :: Int -> Estimate NormalErr Double
poissonNormalCI n
  | n < 0     = error "Statistics.ConfidenceInt.poissonNormalCI negative number of trials"
  | otherwise = estimateNormErr n' (sqrt n')
  where
    n' = fromIntegral n

-- | Calculate confidence intervals for Poisson-distributed value for
--   single measurement. These are exact confidence intervals
poissonCI :: CL Double -> Int -> Estimate ConfInt Double
poissonCI cl@(significanceLevel -> p) n
  | n <  0    = error "Statistics.ConfidenceInt.poissonCI: negative number of trials"
  | n == 0    = estimateFromInterval m (m1,m2) cl
  | otherwise = estimateFromInterval m (m1,m2) cl
  where
    m  = fromIntegral n
    m1 = 0.5 * quantile      (chiSquared (2*n  )) (p/2)
    m2 = 0.5 * complQuantile (chiSquared (2*n+2)) (p/2)

-- | Calculate confidence interval using normal approximation. Note
--   that this approximation breaks down when /p/ is either close to 0
--   or to 1. In particular if @np < 5@ or @1 - np < 5@ this
--   approximation shouldn't be used.
naiveBinomialCI :: Int         -- ^ Number of trials
                -> Int         -- ^ Number of successes
                -> Estimate NormalErr Double
naiveBinomialCI n k
  | n <= 0 || k < 0 = error "Statistics.ConfidenceInt.naiveBinomialCI: negative number of events"
  | k > n           = error "Statistics.ConfidenceInt.naiveBinomialCI: more successes than trials"
  | otherwise       = estimateNormErr p σ
  where
    p = fromIntegral k / fromIntegral n
    σ = sqrt $ p * (1 - p) / fromIntegral n

-- | Clopper-Pearson confidence interval also known as exact
--   confidence intervals.
binomialCI :: CL Double
           -> Int               -- ^ Number of trials
           -> Int               -- ^ Number of successes
           -> Estimate ConfInt Double
binomialCI cl@(significanceLevel -> p) ni ki
  | ni <= 0 || ki < 0 = error "Statistics.ConfidenceInt.binomialCI: negative number of events"
  | ki > ni           = error "Statistics.ConfidenceInt.binomialCI: more successes than trials"
  | ki == 0           = estimateFromInterval eff (0, ub) cl
  | ni == ki          = estimateFromInterval eff (lb,0 ) cl
  | otherwise         = estimateFromInterval eff (lb,ub) cl
  where
    k   = fromIntegral ki
    n   = fromIntegral ni
    eff = k / n
    lb  = quantile      (betaDistr  k      (n - k + 1)) (p/2)
    ub  = complQuantile (betaDistr (k + 1) (n - k)    ) (p/2)

-- $references
--
--  * Clopper, C.; Pearson, E. S. (1934). "The use of confidence or
--    fiducial limits illustrated in the case of the
--    binomial". Biometrika 26: 404–413. doi:10.1093/biomet/26.4.404
--
--  * Brown, Lawrence D.; Cai, T. Tony; DasGupta, Anirban
--    (2001). "Interval Estimation for a Binomial Proportion". Statistical
--    Science 16 (2): 101–133. doi:10.1214/ss/1009213286. MR 1861069.
--    Zbl 02068924.

max_stars_repo_name: intricate/statistics

max_issues_repo_name: intricate/statistics

max_forks_repo_name: intricate/statistics


text:
module RandomUtil.RandomMatrixTest
  (
    testsRandomMatrix
  ) where

import Test.HUnit

import Numeric.LinearAlgebra.HMatrix hiding (corr)
import System.Random

import RandomUtil.RandomMatrix
import Matrix.State

mat :: Matrix Double
mat = reshape 3 $ fromList $ take 9 [1..]

data Two = Two

instance RandomGen Two where
  next g  = (2, g)
  split g = (g, g)

testsRandomMatrix :: Test
testsRandomMatrix = TestList [ TestLabel "selectRandomMatrixElement" testsSelectRandomMatrixElement,
                               TestLabel "randomMatrixIndex" testsRandomMatrixIndex
                             ]

testsSelectRandomMatrixElement :: Test
testsSelectRandomMatrixElement = TestList [TestLabel "test1" testSelects]
testsRandomMatrixIndex :: Test
testsRandomMatrixIndex = TestList [TestLabel "test1" testInBounds]

-- testsSelectRandomMatrixElement
testSelects :: Test
testSelects = TestCase (assertEqual "selectRandomMatrixElement - Select"
                         5
                         (fst $ selectRandomMatrixElement Two mat))

-- testsRandomMatrixIndex
testInBounds :: Test
testInBounds = TestCase (do gen <- getStdGen
                            let listI  = getListRandomIndex gen 100 mat
                            let result = checkAllTrue listI (isValid mat)
                            assertEqual "randomMatrixIndex - Indices in bounds" True result)

-- testsRandomMatrixIndex helper function
isValid :: Matrix a -> MatrixIndex -> Bool
isValid m (i,j) = i >= 0 && j >= 0 && i < rows m && j < cols m

getListRandomIndex :: (Element a) => StdGen -> Int -> Matrix a -> [MatrixIndex]
getListRandomIndex _ 0 _ = []
getListRandomIndex g n m = rIndex : getListRandomIndex g' (n-1) m
  where
    (rIndex,g') = randomMatrixIndex g m

checkAllTrue :: [a] -> (a -> Bool) -> Bool
checkAllTrue xs f = foldr (&&) True $ map f xs

max_stars_repo_name: eklinkhammer/neural-algorithms

max_issues_repo_name: eklinkhammer/neural-algorithms

max_forks_repo_name: eklinkhammer/neural-algorithms


text:
{-# LANGUAGE Strict, StrictData #-}
{-# LANGUAGE FlexibleInstances, UndecidableInstances #-}
{-# LANGUAGE DeriveGeneric, DeriveAnyClass #-}

------------------------------------------------
-- |
-- Module    :  Text.TeXout
-- Copyright :  (c) Jun Yoshida 2019
-- License   :  BSD3
--
-- Export TeX codes.
--
------------------------------------------------

module Text.TeXout where

import GHC.Generics (Generic)
import Control.DeepSeq (NFData)

import Control.Monad
import Control.Monad.ST
import Data.STRef

import Data.Text (Text)
import qualified Data.Text as T

import Data.Map.Strict (Map)
import qualified Data.Map.Strict as Map

import Numeric (showFFloat)

import Numeric.LinearAlgebra (Z)
import Numeric.F2 (F2)
import Numeric.Algebra.FreeModule
import Numeric.Algebra.Frobenius

-- | Class of types which are exhibited in LaTeX display maths.
class TeXMathShow a where
  texMathShow :: a -> Text

-- * Instances
instance (TeXMathShow a) => TeXMathShow (Maybe a) where
  texMathShow = maybe T.empty texMathShow

instance TeXMathShow String where
  texMathShow = T.pack

instance TeXMathShow Int where
  texMathShow = T.pack . show

instance TeXMathShow Double where
  texMathShow x = T.pack $ showFFloat (Just 4) x ""

instance TeXMathShow Z where
  texMathShow = T.pack . show

instance TeXMathShow F2 where
  texMathShow = T.pack . show

instance TeXMathShow SL2B where
  texMathShow SLI = T.singleton '1'
  texMathShow SLX = T.singleton 'X'

instance (TeXMathShow a, TeXMathShow b) => TeXMathShow (a,b) where
  texMathShow (x,y)
    = T.pack "\\left(" <>
      texMathShow x <> T.pack "," <> texMathShow y <>
      T.pack "\\right)"

instance (TeXMathShow a, Ord b, TeXMathShow b) => TeXMathShow (FreeMod a b) where
  texMathShow fm
    = case Map.toList (termMap fm) of
        [] -> texMathShow (0 :: Int)
        ts -> T.intercalate (T.singleton '+') (fmap showProd ts)
    where
      showProd (term,coeff) = texMathShow coeff <> texMathShow term

instance (TeXMathShow a) => TeXMathShow (Map (Int,Int) a) where
  texMathShow mp
    = case foldl rangeFinder Nothing (Map.keys mp) of
        Nothing -> T.empty
        (Just (imin,imax,jmin,jmax))
          -> runST $ do
          stTeX <- newSTRef T.empty
          writeSTRef stTeX $ beginEnvLn "array" [FixArg ("r|" ++ replicate (jmax-jmin+1) 'c')]
          modifySTRef' stTeX (<> mkArrayLn "i\\backslash j" [jmin..jmax])
          modifySTRef' stTeX (<> T.pack "\\\\\\hline\n")
          forM_ [imin..imax] $ \i -> do
            let arrayLn = mkArrayLn i (fmap (\j-> mp Map.!?(i,j)) [jmin..jmax])
            modifySTRef stTeX (<> arrayLn)
            if i < imax
              then (modifySTRef' stTeX (<> T.pack "\\\\\n"))
              else (modifySTRef' stTeX (<> T.singleton '\n'))
          modifySTRef' stTeX (<> endEnv "array")
          readSTRef stTeX
    where
      rangeFinder Nothing (i,j) = Just (i,i,j,j)
      rangeFinder (Just (imin,imax,jmin,jmax)) (i,j) =
        let imin' = if i < imin then i else imin
            imax' = if i > imax then i else imax
            jmin' = if j < jmin then j else jmin
            jmax' = if j > jmax then j else jmax
        in Just (imin',imax',jmin',jmax')

------------------------------
-- * Macros argment control
------------------------------
embrace :: String -> Text
embrace str = '{' `T.cons` T.pack str `T.snoc` '}'

data TeXArg = OptArg String | FixArg String
  deriving (Eq, Ord, Generic, NFData)

encloseArg :: TeXArg -> Text
encloseArg (OptArg arg)
  = case arg of
      []        -> T.empty
      otherwise -> '[' `T.cons` T.pack arg `T.snoc` ']'
encloseArg (FixArg arg) = embrace arg

encloseArgs :: (Foldable t) => t TeXArg -> Text
encloseArgs = foldl (\enc arg -> enc <> encloseArg arg) T.empty

----------------------------
-- * Headers on TeX source
----------------------------
argEnclose :: String -> Text
argEnclose arg = '{' `T.cons` T.pack arg `T.snoc` '}'

optArgEnclose :: String -> Text
optArgEnclose []  = T.empty
optArgEnclose arg = '[' `T.cons` T.pack arg `T.snoc` ']'

documentClass :: String -> String -> Text
documentClass opts cls
  = T.pack "\\documentclass" <> optArgEnclose opts <> argEnclose cls

usePackage :: String -> String -> Text
usePackage opts pkg
  = T.pack "\\usepackage" <> optArgEnclose opts <> argEnclose pkg

------------------------------
-- * Macros and Environments
------------------------------
macro :: (Foldable t) => String -> t TeXArg -> Text
macro csname args
  = '\\' `T.cons` T.pack csname <> encloseArgs args

beginEnv :: (Foldable t) => String -> t TeXArg -> Text
beginEnv env args
  = T.pack "\\begin" <> encloseArg (FixArg env) <> encloseArgs args

beginEnvLn :: (Foldable t) => String -> t TeXArg -> Text
beginEnvLn env args = beginEnv env args `T.snoc` '\n'

endEnv :: String -> Text
endEnv env
  = T.pack "\\end" <> encloseArg (FixArg env)

endEnvLn :: String -> Text
endEnvLn env
  = endEnv env `T.snoc` '\n'

---------------------
-- * Miscellaneous
---------------------
horizontalLine :: Text
horizontalLine
  = T.pack "\\leavevmode\n\\\n\\noindent\\makebox[\\linewidth]{\\rule{\\paperwidth}{0.4pt}}\n"

mkArrayLn :: (TeXMathShow a, TeXMathShow b) => a -> [b] -> Text
mkArrayLn x ys
  = T.intercalate (T.pack " & ") (texMathShow x: fmap texMathShow ys)

mathBB :: String -> Text
mathBB c = macro "mathbb" [FixArg c]

cyclicGrp :: (Integral a) => a -> Text
cyclicGrp n =
  mathBB "Z" <> T.singleton '/' <> texMathShow (fromIntegral n :: Int)

ofRank :: (Integral a, TeXMathShow a) => Text -> a -> Text
ofRank txt 1 = txt
ofRank txt n = txt <> T.pack "^{\\oplus " <> texMathShow n `T.snoc` '}'

flatZip :: Eq a => [a] -> [(Int,a)]
flatZip [] = []
flatZip xs@(_:_) = uncurry (:) $ foldr bin ((1,last xs),[]) (init xs)
  where
    bin y ((n,z),zs)
      | y==z      = ((n+1,z),zs)
      | otherwise = ((1,y),(n,z):zs)

finAbGroup :: (Integral a, Integral a') => a -> [a'] -> Text
finAbGroup freeRk torsions =
  let dsum = if freeRk > 0 && not (null torsions)
             then T.pack "\\oplus "
             else T.empty
      freepart = if freeRk > 0
                 then mathBB "Z" `ofRank` (fromIntegral freeRk :: Int)
                 else T.empty
      torGrps = fmap (\tor -> cyclicGrp (snd tor) `ofRank` fst tor) (flatZip torsions)
      torpart = T.intercalate (T.pack "\\oplus ") torGrps
  in freepart <> dsum <> torpart

max_stars_repo_name: Junology/linkedgram

max_issues_repo_name: Junology/linkedgram

max_forks_repo_name: Junology/linkedgram


text:
{-# LANGUAGE AllowAmbiguousTypes     #-}
{-# LANGUAGE ConstrainedClassMethods #-}
{-# LANGUAGE ConstraintKinds         #-}
{-# LANGUAGE DataKinds               #-}
{-# LANGUAGE FlexibleContexts        #-}
{-# LANGUAGE FlexibleInstances       #-}
{-# LANGUAGE GADTs                   #-}
{-# LANGUAGE InstanceSigs            #-}
{-# LANGUAGE KindSignatures          #-}
{-# LANGUAGE LambdaCase              #-}
{-# LANGUAGE MultiParamTypeClasses   #-}
{-# LANGUAGE PolyKinds               #-}
{-# LANGUAGE RankNTypes              #-}
{-# LANGUAGE ScopedTypeVariables     #-}
{-# LANGUAGE StandaloneDeriving      #-}
{-# LANGUAGE TypeApplications        #-}
{-# LANGUAGE TypeFamilies            #-}
{-# LANGUAGE TypeFamilyDependencies  #-}
{-# LANGUAGE TypeInType              #-}
{-# LANGUAGE TypeOperators           #-}
{-# LANGUAGE UndecidableInstances    #-}
{-# LANGUAGE ViewPatterns            #-}

module TensorOps.Types where

import           Control.Category
import           Control.Monad.Primitive
import           Data.Finite
import           Data.Kind
import           Data.Singletons
import           Data.Type.Combinator
import           Data.Type.Length                  as TCL
import           Data.Type.Product                 as TCP
import           Data.Type.Product.Util            as TCP
import           Data.Type.Sing
import           Data.Type.Uniform
import           Data.Type.Vector
import           Prelude hiding                    ((.), id)
import           Statistics.Distribution
import           System.Random.MWC
import           TensorOps.NatKind
import           Type.Class.Higher
import           Type.Class.Known
import           Type.Class.Witness
import           Type.Family.List

{-# RULES
"realToFrac/Double->Double" realToFrac = id :: Double -> Double
"realToFrac/Float->Float" realToFrac = id :: Float -> Float
    #-}

class NatKind k => Tensor (t :: [k] -> Type) where
    type ElemT t  :: Type

    -- TODO: can we detach Vec from liftT ?
    liftT   :: SingI o
            => (Vec n (ElemT t) -> ElemT t)
            -> Vec n (t o)
            -> t o
    gmul    :: (SingI (Reverse os ++ ns), SingI (ms ++ ns))
            => Length ms
            -> Length os
            -> Length ns
            -> t (ms         ++ os)
            -> t (Reverse os ++ ns)
            -> t (ms         ++ ns)
    -- is a list really the best data structure here?
    -- maybe Foldable f?
    sumT    :: SingI o => [t o] -> t o
    scaleT  :: SingI o => ElemT t -> t o -> t o
    transp  :: (SingI ns, SingI (Reverse ns))
            => t ns
            -> t (Reverse ns)
    -- sumRow  :: Remove ns n ms
    --         -> t ns
    --         -> t ms
    mapRows :: SingI (ns ++ ms)
            => Length ns
            -> (t ms -> t ms)
            -> t (ns ++ ms)
            -> t (ns ++ ms)
    sumRows :: (SingI (n ': ns), SingI ns)
            => t (n ': ns)
            -> t ns
    diag    :: SingI (n ': ns)
            => Uniform n ns
            -> t '[n]
            -> t (n ': ns)
    getDiag :: SingI n
            => Uniform n ns
            -> t (n ': n ': ns)
            -> t '[n]
    genRand :: (ContGen d, PrimMonad m, SingI ns)
            => d
            -> Gen (PrimState m)
            -> m (t ns)
    generateA :: (Applicative f, SingI ns)
              => (Prod (IndexN k) ns -> f (ElemT t))
              -> f (t ns)
    ixRows
        :: (Applicative f, SingI (ms ++ os))
        => Length ms
        -> Length os
        -> (Prod (IndexN k) ms -> t ns -> f (t os))
        -> t (ms ++ ns)
        -> f (t (ms ++ os))
    (!) :: t ns
        -> Prod (IndexN k) ns
        -> ElemT t

-- type TensorOp = OpPipe TOp

-- | Function and gradient
data VFunc n
    = VF { vfFunc :: !(forall a. RealFloat a => Vec n a -> a      )
         , vfGrad :: !(forall a. RealFloat a => Vec n a -> Vec n a)
         }

-- -- | A kludge to get around lack of impredicative types in Haskell
-- newtype VFunc n = VF { getVF :: forall a. RealFloat a => Vec n a -> a }

data TOp :: [[k]] -> [[k]] -> Type where
    TOp :: { runTOp   :: !(forall t. (Tensor t, RealFloat (ElemT t)) => Prod t ns -> Prod t ms)
           , gradTOp' :: !(forall t. (Tensor t, RealFloat (ElemT t)) => Prod t ns -> Prod t ms -> Prod t ns)
           } -> TOp ns ms

gradTOp
    :: (Tensor t, RealFloat (ElemT t))
    => TOp ns '[ '[] ]
    -> Prod t ns
    -> Prod t ns
gradTOp o xs = gradTOp' o xs (only (getI $ generateA (\_ -> I 1)))

instance Category TOp where
    id = TOp id
             (flip const)
    {-# INLINE id #-}

    (.) :: forall as bs cs. ()
        => TOp bs cs
        -> TOp as bs
        -> TOp as cs
    TOp f2 g2 . TOp f1 g1 = TOp f3 g3
      where
        f3  :: forall t. (Tensor t, RealFloat (ElemT t))
            => Prod t as
            -> Prod t cs
        f3 = f2 . f1
        {-# INLINE f3 #-}
        g3  :: forall t. (Tensor t, RealFloat (ElemT t))
            => Prod t as
            -> Prod t cs
            -> Prod t as
        g3 xs ds = g1 xs (g2 (f1 xs) ds)
        {-# INLINE g3 #-}
    {-# INLINE (.) #-}

idOp
    :: forall ns. ()
    => TOp ns ns
idOp = id
{-# INLINE idOp #-}

firstOp
    :: forall os ns ms. (Known Length ns, Known Length ms)
    => TOp ns ms
    -> TOp (ns ++ os) (ms ++ os)
firstOp (TOp f g) = TOp f' g'
  where
    f'  :: forall t. (Tensor t, RealFloat (ElemT t))
        => Prod t (ns ++ os)
        -> Prod t (ms ++ os)
    f' = overProdInit @os known f
    {-# INLINE f' #-}
    g'  :: forall t. (Tensor t, RealFloat (ElemT t))
        => Prod t (ns ++ os)
        -> Prod t (ms ++ os)
        -> Prod t (ns ++ os)
    g' (takeProd @os known -> xs) = overProdInit @os known (g xs)
    {-# INLINE g' #-}
{-# INLINE firstOp #-}

secondOp
    :: forall os ns ms. Known Length os
    => TOp ns ms
    -> TOp (os ++ ns) (os ++ ms)
secondOp (TOp f g) = TOp f' g'
  where
    f'  :: forall t. (Tensor t, RealFloat (ElemT t))
        => Prod t (os ++ ns)
        -> Prod t (os ++ ms)
    f' = overProdTail @os known f
    {-# INLINE f' #-}
    g'  :: forall t. (Tensor t, RealFloat (ElemT t))
        => Prod t (os ++ ns)
        -> Prod t (os ++ ms)
        -> Prod t (os ++ ns)
    g' (dropProd @os known -> xs) = overProdTail @os known (g xs)
    {-# INLINE g' #-}
{-# INLINE secondOp #-}

(*>>)
    :: forall as bs cs ds. (Known Length as, Known Length bs)
    => TOp as bs
    -> TOp (bs ++ cs) ds
    -> TOp (as ++ cs) ds
t1 *>> t2 = firstOp @cs t1 >>> t2
infixr 0 *>>
{-# INLINE (*>>) #-}

(<<*)
    :: forall as bs cs ds. (Known Length as, Known Length bs)
    => TOp (bs ++ cs) ds
    -> TOp as bs
    -> TOp (as ++ cs) ds
(<<*) = flip ((*>>) @as @bs @cs @ds)
infixr 2 <<*
{-# INLINE (<<*) #-}

(***)
    :: forall as bs cs ds. (Known Length as, Known Length cs)
    => TOp as cs
    -> TOp bs ds
    -> TOp (as ++ bs) (cs ++ ds)
TOp f1 g1 *** TOp f2 g2 = TOp f3 g3
  where
    f3  :: forall t. (Tensor t, RealFloat (ElemT t))
        => Prod t (as ++ bs)
        -> Prod t (cs ++ ds)
    f3 = overProdSplit known f1 f2
    {-# INLINE f3 #-}
    g3  :: forall t. (Tensor t, RealFloat (ElemT t))
        => Prod t (as ++ bs)
        -> Prod t (cs ++ ds)
        -> Prod t (as ++ bs)
    g3 (splitProd known->(xs, ys)) = overProdSplit known (g1 xs) (g2 ys)
    {-# INLINE g3 #-}
{-# INLINE (***) #-}

(&&&)
    :: forall as bs cs. (Known Length bs, SingI as)
    => TOp as bs
    -> TOp as cs
    -> TOp as (bs ++ cs)
TOp f1 g1 &&& TOp f2 g2 = TOp f3 g3
  where
    f3  :: forall t. (Tensor t, RealFloat (ElemT t))
        => Prod t as
        -> Prod t (bs ++ cs)
    f3 = TCP.append' <$> f1 <*> f2
    {-# INLINE f3 #-}
    g3  :: forall t. (Tensor t, RealFloat (ElemT t))
        => Prod t as
        -> Prod t (bs ++ cs)
        -> Prod t as
    g3 xs (splitProd known->(dtdys,dtdzs)) =
        zipProdWith3 (\s gxy gxz -> sumT [gxy,gxz] \\ s)
                     (singProd sing)
                     (g1 xs dtdys)
                     (g2 xs dtdzs)
    {-# INLINE g3 #-}
{-# INLINE (&&&) #-}

-- -- | TODO: replace with `syntactic`?
-- data OpPipe :: ([k] -> [k] -> Type) -> [k] -> [k] -> Type where
--     OPØ   :: OpPipe f a a
--     Pop   :: !(Sing a)
--           -> !(Sing b)
--           -> !(Sing d)
--           -> !(f a b)
--           -> !(OpPipe f (b ++ d) c)
--           -> OpPipe f (a ++ d) c

-- pappend
--     :: forall a b c d f. ()
--     => Sing a
--     -> Sing b
--     -> Sing d
--     -> OpPipe f a b
--     -> OpPipe f (b ++ d) c
--     -> OpPipe f (a ++ d) c
-- pappend _ sB sD = \case
--     OPØ -> id
--     Pop (sA' :: Sing a')
--         (sB' :: Sing b')
--         (sD' :: Sing d')
--         (x   :: f a' b'  )
--         (xs  :: OpPipe f (b' ++ d') b)
--           -> \ys -> let lD' :: Length d'
--                         lD' = singLength sD'
--                     in  Pop sA' sB' (sD' %:++ sD) x (pappend (sB' %:++ sD') sB sD xs ys)
--                           \\ appendAssoc (singLength sA') lD' lD
--                           \\ appendAssoc (singLength sB') lD' lD
--   where
--     lD :: Length d
--     lD = singLength sD

-- pipe
--     :: forall t a b. (SingI a, SingI b)
--     => t a b
--     -> OpPipe t a b
-- pipe o = Pop sing sing SNil o OPØ
--            \\ appendNil (singLength (sing :: Sing a))
--            \\ appendNil (singLength (sing :: Sing b))

-- pop :: forall a b c d f. (SingI a, SingI b, SingI d)
--     => Length d
--     -> f a b
--     -> OpPipe f (b ++ d) c
--     -> OpPipe f (a ++ d) c
-- pop _ = Pop (sing :: Sing a) (sing :: Sing b) (sing :: Sing d)

-- infixr 4 ~.
-- (~.)
--     :: (SingI a, SingI b, SingI d)
--     => (Length a, Length d, f a b)
--     -> OpPipe f (b ++ d) c
--     -> OpPipe f (a ++ d) c
-- (_, lD, x) ~. y = pop lD x y

instance Eq1 Finite

max_stars_repo_name: mstksg/tensor-ops

max_issues_repo_name: mstksg/tensor-ops

max_forks_repo_name: mstksg/tensor-ops


text:
module FFT.FFTControl(
    fftbflySPipeline,
    fftdc,
    fftdcPipeline
) where
import Data.Complex
import FFT.Samples
import System.Environment
import Control.Parallel
import Strategies
import Criterion.Main
import FFT.Orig

-- twiddle factors
tw :: Int -> Int -> Complex Float
tw n k = cis (-2 * pi * fromIntegral k / fromIntegral n)

fftdc as = dc split threshold combine worker as
	where
                combine [xs] = xs
		combine [xs,ys] = interleave xs ys
		threshold = (\x -> length x <= 1) --(floor $ sqrt $ fromIntegral as))
		worker a = [a]
		split xs = [cs, ds]
			where (cs, ds) = bflyS xs

fftdcPipeline as = dc split threshold combine worker as
	where
                combine [xs] = xs
		combine [xs,ys] = interleave xs ys
		threshold = (\x -> length x <= 1) --(floor $ sqrt $ fromIntegral as))
		worker a = [a]
		split xs = [cs, ds]
			where (cs, ds) = bflySPipeline xs

fftbflySPipeline as = interleave ls rs
    where
        (cs,ds) = bflySPipeline as
        ls = fft cs
        rs = fft ds

interleave :: [Complex Float] -> [Complex Float] -> [Complex Float]
interleave [] bs = bs
interleave (a:as) bs = a : interleave bs as

bflyS :: [Complex Float] -> ([Complex Float], [Complex Float])
bflyS as = (los,rts)
  where
    (ls,rs) = halve as
    los = zipWith (+) ls rs
    ros = zipWith (-) ls rs
    rts = zipWith (*) ros [tw (length as) i | i <- [0..(length ros) - 1]]

bflySPipeline :: [Complex Float] -> ([Complex Float], [Complex Float])
bflySPipeline as = (los,rts)
    where
        (ls, rs) = halve as
        los = zipWith (+) ls rs
        rts = parpipeline2 (zipWith (-) ls) (zipWith (*) [tw (length as) i | i <- [0..(length rs) - 1]]) rs

-- split the input into two halves
halve as = splitAt n' as
  where
    n' = div (length as + 1) 2

max_stars_repo_name: BNJHope/parallel-fourier-transforms

max_issues_repo_name: BNJHope/parallel-fourier-transforms

max_forks_repo_name: BNJHope/parallel-fourier-transforms