text:
module Backprop
  (BackpropNet, buildBackpropNet, logisticSigmoidAS, tanhAS, identityAS)
    where

import MatrixPlus as P
import NeuralNet

import Numeric.LinearAlgebra as N
import Test.QuickCheck
import Control.Exception
import Debug.Trace

-- {{{ Layers

-- | An individual layer in a neural network, prior to propagation
data Layer = Layer
    {
      -- The weights for this layer
      lW :: Matrix Double,
      -- The activation specification for this layer
      lAS :: ActivationSpec
    }

instance Show Layer where
    show layer = "w=" ++ show (lW layer) ++ ", activation spec=" ++ show (lAS layer)

-- | Generate a layer of the specified "size", with arbitrary values.
-- | To see sample values (of "size" 4), in GHCi type: sample (sizedArbLayer 4)
sizedArbLayer :: Int -> Gen Layer
sizedArbLayer n = do
    -- QuickCheck passes a value n >= 0, but we can't create a vector of length 0, so use n+1
    r <- choose(1, n+1)
    c <- choose(1, n+1)
    w <- sizedArbWeightMatrix r c
    s <- (arbitrary :: Gen ActivationSpec)
    return Layer{ lW=w, lAS=s }

-- | Generate a layer of the specified "size" and the specified input width, with arbitrary values.
-- | To see sample values (of input width 3, "size" 4), in GHCi type: sample (sizedArbLayer' 3 4)
sizedArbLayer' :: Int -> Int -> Gen Layer
sizedArbLayer' c n = do
    -- QuickCheck passes a value n >= 0, but we can't create a vector of length 0, so use n+1
    r <- choose(1, n+1)
    w <- sizedArbWeightMatrix r c
    s <- (arbitrary :: Gen ActivationSpec)
    return Layer{ lW=w, lAS=s }

-- | Generate a layer of arbitrary "size", with arbitrary values.
-- | To see some sample values, in GHCi type: sample' arbitrary :: IO [Layer]
instance Arbitrary Layer where
    arbitrary = sized sizedArbLayer

inputWidth :: Layer -> Int
inputWidth = cols . lW

outputWidth :: Layer -> Int
outputWidth = rows . lW

-- }}}

-- {{{ Propagation

-- | An individual layer in a neural network, after propagation but prior to backpropagation
data PropagatedLayer
    = PropagatedLayer
        {
          -- The input to this layer
          pIn :: ColumnVector Double,
          -- The output from this layer
          pOut :: ColumnVector Double,
          -- The value of the first derivative of the activation function for this layer
          pF'a :: ColumnVector Double,
          -- The weights for this layer
          pW :: Matrix Double,
          -- The activation specification for this layer
          pAS :: ActivationSpec
        }
    | PropagatedSensorLayer
        {
          -- The output from this layer
          pOut :: ColumnVector Double
        }

instance Show PropagatedLayer where
    show (PropagatedLayer x y f'a w s) =
        "in=" ++ show x
        ++ ", out=" ++ show y
        ++ ", f'(a)=" ++ show f'a
        ++ ", w=" ++ show w
        ++ ", " ++ show s
    show (PropagatedSensorLayer x) = "out=" ++ show x

    -- {{{ Testing

{-
-- | Generate a propagated hidden or output layer of the specified "size", with arbitrary values.
-- | To see sample values (of "size" 4), in GHCi type: sample (sizedArbPropagatedLayer 4)
sizedArbPropagatedLayer :: Int -> Gen PropagatedLayer
sizedArbPropagatedLayer n = do
    -- QuickCheck passes a value n >= 0, but we can't create a vector of length 0, so use n+1
    r <- choose(1, n+1)
    c <- choose(1, n+1)
    w <- sizedArbWeightMatrix r c
    y <- sizedArbColumnVector r
    return PropagatedLayer{ pOut=y
                          , pF'a=(zeroColumnVector r)
                          , pW=w
                          , pF=id
                          , pF'=id'
                          }

-- | Generate a propagated sensor layer of the specified "size", with arbitrary values.
-- | To see sample values (of "size" 4), in GHCi type: sample (sizedArbPropagatedLayer2 4)
sizedArbPropagatedLayer0 :: Int -> Gen PropagatedLayer
sizedArbPropagatedLayer0 n = do
    -- QuickCheck passes a value n >= 0, but we can't create a vector of length 0, so use n+1
    r <- choose(1, n+1)
    y <- sizedArbColumnVector r
    return PropagatedSensorLayer{ pOut=y }

-- | Generate a propagated layer of arbitrary "size", with arbitrary values.
-- | To see some sample values, in GHCi type: sample' arbitrary :: IO [PropagatedLayer]
instance Arbitrary PropagatedLayer where
    arbitrary = frequency
        [ (4, sized sizedArbPropagatedLayer)
        , (1, sized sizedArbPropagatedLayer0) ]
-}
    -- }}}

-- | Propagate the inputs through this layer to produce an output.
propagate :: PropagatedLayer -> Layer -> PropagatedLayer
propagate layerJ layerK = PropagatedLayer
        {
          pIn = x,
          pOut = y,
          pF'a = f'a,
          pW = w,
          pAS = lAS layerK
        }
  where x = pOut layerJ
        w = lW layerK
        a = w <> x
        f = asF ( lAS layerK )
        y = P.mapMatrix f a
        f' = asF' ( lAS layerK )
        f'a = P.mapMatrix f' a

-- }}}

-- {{{ Backpropagation

-- | An individual layer in a neural network, after backpropagation
data BackpropagatedLayer = BackpropagatedLayer
    {
      -- Del-sub-z-sub-l of E
      bpDazzle :: ColumnVector Double,
      -- The error due to this layer
      bpErrGrad :: ColumnVector Double,
      -- The value of the first derivative of the activation
      --   function for this layer
      bpF'a :: ColumnVector Double,
      -- The input to this layer
      bpIn :: ColumnVector Double,
      -- The output from this layer
      bpOut :: ColumnVector Double,
      -- The weights for this layer
      bpW :: Matrix Double,
      -- The activation specification for this layer
      bpAS :: ActivationSpec
    }

instance Show BackpropagatedLayer where
    show layer =
        "dazzle=" ++ show (bpDazzle layer)
        ++ ", grad=" ++ show (bpErrGrad layer)
        ++ ", in=" ++ show (bpIn layer)
        ++ ", out=" ++ show (bpOut layer)
        ++ ", w=" ++ show (bpW layer)
        ++ ", activationFunction=?, activationFunction'=?"

backpropagateFinalLayer ::
    PropagatedLayer -> ColumnVector Double -> BackpropagatedLayer
backpropagateFinalLayer l t = BackpropagatedLayer
    {
      bpDazzle = dazzle,
      bpErrGrad = errorGrad dazzle f'a (pIn l),
      bpF'a = pF'a l,
      bpIn = pIn l,
      bpOut = pOut l,
      bpW = pW l,
      bpAS = pAS l
    }
    where dazzle =  pOut l - t
          f'a = pF'a l

errorGrad :: ColumnVector Double -> ColumnVector Double -> ColumnVector Double
    -> ColumnVector Double
errorGrad dazzle f'a input = (dazzle * f'a) <> trans input

-- | Propagate the inputs backward through this layer to produce an output.
backpropagate :: PropagatedLayer -> BackpropagatedLayer -> BackpropagatedLayer
backpropagate layerJ layerK = BackpropagatedLayer
    {
      bpDazzle = dazzleJ,
      bpErrGrad = errorGrad dazzleJ f'aJ (pIn layerJ),
      bpF'a = pF'a layerJ,
      bpIn = pIn layerJ,
      bpOut = pOut layerJ,
      bpW = pW layerJ,
      bpAS = pAS layerJ
    }
    where dazzleJ = wKT <> (dazzleK * f'aK)
          dazzleK = bpDazzle layerK
          wKT = trans ( bpW layerK )
          f'aK = bpF'a layerK
          f'aJ = pF'a layerJ

-- }}}

-- {{{ Adjusting weights after backpropagation

update :: Double -> BackpropagatedLayer -> Layer
update rate layer = Layer
        {
          lW = wNew,
          lAS = bpAS layer
        }
    where wOld = bpW layer
          delW = rate `scale` bpErrGrad layer
          wNew = wOld - delW

-- }}}

-- {{{ Building a network

data BackpropNet = BackpropNet
    {
      layers :: [Layer],
      learningRate :: Double
    } deriving Show

buildBackpropNet ::
  -- The learning rate
  Double ->
  -- The weights for each layer
  [Matrix Double] ->
  -- The activation specification (used for all layers)
  ActivationSpec ->
  -- The network
  BackpropNet
buildBackpropNet lr ws s = BackpropNet { layers=ls, learningRate=lr }
  where checkedWeights = scanl1 checkDimensions ws
        ls = map buildLayer checkedWeights
        buildLayer w = Layer { lW=w, lAS=s }

checkDimensions :: Matrix Double -> Matrix Double -> Matrix Double
checkDimensions w1 w2 =
  if rows w1 == cols w2
       then w2
       else error "Inconsistent dimensions in weight matrix"

propagateNet :: ColumnVector Double -> BackpropNet -> [PropagatedLayer]
propagateNet input net = tail calcs
  where calcs = scanl propagate layer0 (layers net)
        layer0 = PropagatedSensorLayer{ pOut=validatedInputs }
        validatedInputs = validateInput net input

validateInput :: BackpropNet -> ColumnVector Double -> ColumnVector Double
validateInput net = validateInputValues . validateInputDimensions net

validateInputDimensions ::
    BackpropNet ->
    ColumnVector Double ->
    ColumnVector Double
validateInputDimensions net input =
  if got == expected
       then input
       else error ("Input pattern has " ++ show got ++ " bits, but " ++ show expected ++ " were expected")
           where got = rows input
                 expected = inputWidth (head (layers net))

validateInputValues :: ColumnVector Double -> ColumnVector Double
validateInputValues input =
  if (min >= 0) && (max <= 1)
       then input
       else error "Input bits outside of range [0,1]"
       where min = minimum ns
             max = maximum ns
             ns = toList ( flatten input )

backpropagateNet ::
  ColumnVector Double -> [PropagatedLayer] -> [BackpropagatedLayer]
backpropagateNet target layers = scanr backpropagate layerL hiddenLayers
  where hiddenLayers = init layers
        layerL = backpropagateFinalLayer (last layers) target

-- }}}

-- {{{ Define BackpropNet to be an instance of Neural Net

instance NeuralNet BackpropNet where
  evaluate = evaluateBPN
  train = trainBPN

evaluateBPN :: BackpropNet -> [Double] -> [Double]
evaluateBPN net input = columnVectorToList( pOut ( last calcs ))
  where calcs = propagateNet x net
        x = listToColumnVector (1:input)

trainBPN :: BackpropNet -> [Double] -> [Double] -> BackpropNet
trainBPN net input target = BackpropNet { layers=newLayers, learningRate=rate }
  where newLayers = map (update rate) backpropagatedLayers
        rate = learningRate net
        backpropagatedLayers = backpropagateNet (listToColumnVector target) propagatedLayers
        propagatedLayers = propagateNet x net
        x = listToColumnVector (1:input)

-- }}}

-- {{{ General Testing

-- | A layer with suitable input and target vectors, suitable for testing.
data LayerTestData = LTD (ColumnVector Double) Layer (ColumnVector Double)
  deriving Show

-- | Generate a layer with suitable input and target vectors, of the specified "size",
-- | with arbitrary values.
-- | To see sample values (of "size" 4), in GHCi type: sample (sizedLayerTestData 4)
sizedLayerTestData :: Int -> Gen LayerTestData
sizedLayerTestData n = do
    l <- sizedArbLayer n
    x <- sizedArbColumnVector (inputWidth l)
    t <- sizedArbColumnVector (outputWidth l)
    return (LTD x l t)

instance Arbitrary LayerTestData where
  -- | To see sample values, in GHCi type: sample arbLayerTestData
  arbitrary = sized sizedLayerTestData

-- | Training reduces error in the final (output) layer
prop_trainingReducesFinalLayerError :: LayerTestData -> Property
prop_trainingReducesFinalLayerError (LTD x l t) =
    -- (collect l) . -- uncomment to view test data
    (classifyRange "len x " n 0 25) .
    (classifyRange "len x " n 26 50) .
    (classifyRange "len x " n 51 75) .
    (classifyRange "len x " n 76 100) $
    errorAfter < errorBefore || errorAfter < 0.01
        where n = inputWidth l
              pl0 = PropagatedSensorLayer{ pOut=x }
              pl = propagate pl0 l
              bpl = backpropagateFinalLayer pl t
              errorBefore = P.magnitude (t - pOut pl)
              lNew = update 0.0000000001 bpl -- make sure we don't overshoot the mark
              plNew = propagate pl0 lNew
              errorAfter =  P.magnitude (t - pOut plNew)

iterateTraining n pl0 r ltd = pOut plFinal
    where iterations = iterate (trainOneLayer pl0 r) ltd
          (LTD _ lFinal _) = last (take n iterations)
          plFinal = propagate pl0 lFinal

trainOneLayer :: PropagatedLayer -> Double -> LayerTestData -> LayerTestData
trainOneLayer pl0 r (LTD x l t) = LTD x lNew t
    where pl = propagate pl0 l
          bpl = backpropagateFinalLayer pl t
          lNew = update r bpl

-- | Testable property:
-- | Training a single layer with the same input repeatedly will eventually yield the target
prop_trainingOneLayerWithOneInputYieldsPerfection :: LayerTestData -> Property
prop_trainingOneLayerWithOneInputYieldsPerfection (LTD x l t) =
    -- (collect l) . -- uncomment to view test data
    (classifyRange "len x " n 0 25) .
    (classifyRange "len x " n 26 50) .
    (classifyRange "len x " n 51 75) .
    (classifyRange "len x " n 76 100) $
    e < 0.1
        where n = inputWidth l
              r = 1
              pl0 = PropagatedSensorLayer{ pOut=x }
              y = iterateTraining 100 pl0 r (LTD x l t)
              e = P.magnitude (t - y)

-- | A layer with suitable input and target vectors, suitable for testing.
data TwoLayerTestData =
  TLTD (ColumnVector Double) Layer Layer (ColumnVector Double)
    deriving Show

-- | Generate a layer with suitable input and target vectors, of the specified "size",
-- | with arbitrary values.
-- | To see sample values (of "size" 4), in GHCi type: sample (sizedTwoLayerTestData 4)
sizedTwoLayerTestData :: Int -> Gen TwoLayerTestData
sizedTwoLayerTestData n = do
    l1 <- sizedArbLayer n
    l2 <- sizedArbLayer' (outputWidth l1) n
    x <- sizedArbColumnVector (inputWidth l1)
    t <- sizedArbColumnVector (outputWidth l2)
    return (TLTD x l1 l2 t)

instance Arbitrary TwoLayerTestData where
  -- | To see sample values, in GHCi type: sample arbTwoLayerTestData
  arbitrary = sized sizedTwoLayerTestData

-- | Training reduces error in a hidden layer
prop_trainingReducesHiddenLayerError :: TwoLayerTestData -> Property
prop_trainingReducesHiddenLayerError (TLTD x l1 l2 t)=
    -- (collect l) . -- uncomment to view test data
    (classifyRange "len x " n 0 25) .
    (classifyRange "len x " n 26 50) .
    (classifyRange "len x " n 51 75) .
    (classifyRange "len x " n 76 100) $
    errorAfter < errorBefore || errorAfter < 0.01
        where n = inputWidth l1
              pl0 = PropagatedSensorLayer{ pOut=x }
              pl1 = propagate pl0 l1
              pl2 = propagate pl1 l2
              bpl2 = backpropagateFinalLayer pl2 t
              bpl1 = backpropagate pl1 bpl2
              errorBefore = P.magnitude (t - pOut pl2)
              l1New = update 0.00000001 bpl1 -- make sure we don't overshoot the mark
              -- leave layer 2 alone, we're only interested in layer 1
              pl1New = propagate pl0 l1New
              pl2New = propagate pl1New l2
              errorAfter = P.magnitude (t - pOut pl2New)

classifyRange :: Testable a => String -> Int -> Int -> Int -> a -> Property
classifyRange s n n0 n1 =
    classify (n >= n0 && n <= n1) (s ++ show n0 ++ ".." ++ show n1)

-- }}}

-- {{{ Common activation functions

data ActivationSpec = ActivationSpec
    {
      asF :: Double -> Double,
      asF' :: Double -> Double,
      desc :: String
    }

instance Show ActivationSpec where
    show = desc

identityAS = ActivationSpec
    {
      asF = id,
      asF' = const 1,
      desc = "identity"
    }

logisticSigmoidAS :: Double -> ActivationSpec
logisticSigmoidAS c = ActivationSpec
    {
        asF = logisticSigmoid c,
        asF' = logisticSigmoid' c,
        desc = "logistic sigmoid, c=" ++ show c
    }

arbitraryLogisticAS :: Gen ActivationSpec
arbitraryLogisticAS = do
    c <- choose(0,1) -- TODO Can it be > 1?
    return (logisticSigmoidAS c)

instance Arbitrary ActivationSpec where
--    arbitrary = oneof [ return identityAS, arbitraryLogisticAS ]
    arbitrary = return identityAS

logisticSigmoid :: (Field a, Floating a) => a -> a -> a
logisticSigmoid c a = 1 / (1 + exp((-c) * a))

logisticSigmoid' :: (Field a, Floating a) => a -> a -> a
logisticSigmoid' c a = (c * f a) * (1 - f a)
  where f = logisticSigmoid c

tanhAS :: ActivationSpec
tanhAS = ActivationSpec
    {
      asF = tanh,
      asF' = tanh',
      desc = "tanh"
    }

tanh' x = 1 - (tanh x)^2

-- }}}

max_stars_repo_name: mhwombat/backprop-example

max_issues_repo_name: mhwombat/backprop-example

max_forks_repo_name: mhwombat/backprop-example


text:
{-# LANGUAGE BangPatterns #-}

import Gauge.Main (defaultMain, bench, whnf)

import Data.Int (Int64)
import Data.List (unfoldr)
import Data.Maybe (fromJust)
import Data.Primitive.Contiguous (PrimArray, fromListN)
import Statistics.Array.Types (AscList)
import System.Random (random, mkStdGen)

import qualified Data.Primitive as PM
import qualified Statistics.Array as Stats
import qualified Statistics.Array.Types as Asc

main :: IO ()
main = defaultMain
  [ bench "sort-16" $ whnf fromAscInt input16
  , bench "sort-128" $ whnf fromAscInt input128
  , bench "sort-1024" $ whnf fromAscInt input1024
  , bench "sort-65536" $ whnf fromAscInt input65536
  , bench "list-derivative" $ whnf performListDerivativeInt asc1024
  , bench "median-of-absolute-deviations" $ whnf Stats.mad asc1024
  ]
  where
  input16, input128, input1024, input65536 :: PrimArray Int64
  !input16 = fromListN 16 $ take 16 . drop 0 $ inputInf
  !input128 = fromListN 128 $ take 128 . drop 16 $ inputInf
  !input1024 = fromListN 1024 $ take 1024 . drop (16+128) $ inputInf
  !input65536 = fromListN 65536 $ take 65536 . drop (16+128+1024) $ inputInf
  asc1024 :: AscList Int64
  asc1024 = fromJust $ Asc.fromArray input1024

fromAscInt :: PrimArray Int64 -> Int64
{-# noinline fromAscInt #-}
fromAscInt !x = case Asc.fromArray x of
  Nothing -> 0
  Just (Asc.AscList y) -> fromIntegral (PM.sizeofPrimArray y)

-- This is here to make it easy to inspect Core to confirm that specialization
-- works correctly.
performListDerivativeInt :: AscList Int64 -> PrimArray Int64
{-# noinline performListDerivativeInt #-}
performListDerivativeInt !x = Stats.listDerivative x

inputInf :: [Int64]
inputInf = unfoldr (Just . random) seed
  where
  seed = mkStdGen 1234567

ext: hs

max_stars_repo_name: Zankoku-Okuno/array-statistics

max_issues_repo_name: Zankoku-Okuno/array-statistics

max_forks_repo_name: Zankoku-Okuno/array-statistics


text:
module Hw5
    (Signal,
    Vec(..),
    rd,
    absolute,
    range,
    dft,
    idft,
    realV,
    low_pass',
    low_pass,
    fft,
    ifft,
    makeNoise)
    where

import Data.Complex
import Debug.Trace
import GHC.Base (liftA2)
import System.Random
-- import Control.Monad.Reader

-- {-# LANGUAGE  KindSignatures  #-}
-- {-# LANGUAGE  MultiParamTypeClasses  #-}

newtype Vec a = Vec {runVec :: [a]}
newtype  Reader r a = Reader { runReader  :: r -> a }
type Signal = Vec (Complex Double)

imagV :: Num a => Vec a -> Vec (Complex a)
realV :: Num a => Vec a -> Vec (Complex a)
range :: Double -> Double -> Double -> [Double]
absolute :: Vec (Complex Double) -> Vec Double
rd :: Int -> Vec Double -> Vec Double
twiddle :: Double -> Double -> Signal
dft :: Signal -> Signal
idft :: Signal -> Signal
mask :: Int -> Int -> Signal
fft :: Signal -> Signal
ifft :: Signal -> Signal
makeNoise :: Random a => Int -> Int -> a -> a -> Vec a

makeNoise seed n low high =
    let x = mkStdGen seed
        genNum l h = randomR (l, h)
        genList 1 l h g = [c] where (c,y) = genNum l h g
        genList n l h g = c : genList (n-1) l h y where (c,y) = genNum l h g
        in
            Vec $ genList n low high x

length' (Vec x) = length x

rd n = fmap (\x -> fromIntegral (round $ c * x) / c)
    where c = 10^n

absolute = fmap (\(r:+i) -> sqrt(r*r + i*i))

range from to count = map (\x -> from + x * step) [0..count-1]
    where step = (to - from)/count

imagV (Vec a) = Vec $ map (0:+) a
realV (Vec a) = Vec $ map (:+0) a

instance Show a => Show (Vec a) where
    show (Vec lst) = "[" ++ drop 1 lst' ++ "]"
        where
            lst' = mconcat $ map (\x -> " " ++ show x) lst

instance Functor Vec where
    fmap f (Vec x) = Vec $ map f x

instance Applicative Vec where
    pure = Vec . repeat
    (Vec f) <*> (Vec x) = Vec $ zipWith ($) f x
    liftA2 f (Vec x) (Vec y) = Vec $ zipWith f x y

instance Foldable Vec where
    foldr f c (Vec a) = Prelude.foldr f c a

instance Semigroup (Vec a) where
    Vec a <> Vec b = Vec $ a ++ b

instance Monoid (Vec a) where
    mempty = Vec []

instance Num a => Num (Vec a) where
    (+) = liftA2 (+)
    (-) = liftA2 (-)
    (*) = liftA2 (*)
    negate = fmap negate
    abs = fmap abs
    signum = fmap signum
    fromInteger x = pure $ fromInteger x

instance (Floating a) => Fractional (Vec a) where
    (/) = liftA2 (/)
    fromRational x = pure $ fromRational x

instance (Floating a) => Floating (Vec a) where
    pi = pure pi
    exp = fmap exp
    log = fmap log
    sin = fmap sin
    cos = fmap cos
    asin = fmap asin
    acos = fmap acos
    atan = fmap atan
    sinh = fmap sinh
    cosh = fmap cosh
    asinh = fmap asinh
    acosh = fmap acosh
    atanh = fmap atanh

mask freq n = realV $ Vec $ one ++ zero ++ one where
    one = (take freq $ repeat 1)
    zero = (take (n-freq*2) $ repeat 0)

low_pass' :: Int -> Signal -> Signal

low_pass' v sig =
    let
        m = mask v $ length sig
    in idft $ m * dft sig

twiddle n k =
    let a = [0.. (n - 1)]
        f x = exp (0:+ x * (-2*pi*k / n))
    in Vec $ fmap f a

dft x =
    let bigN = fromIntegral $ length x -- length of x
        n = Vec [0 .. bigN-1]
        f k = sum $ x * twiddle bigN k
    in fmap f n

idft x =
    let bigN = fromIntegral $ length x -- length of x
        cX = fmap conjugate x
        n = Vec [0 .. bigN-1]
        diC = (1.0 / bigN) :+ 0
        f k = sum $ cX * twiddle bigN k
    in conj diC $ fmap f n

conj :: Complex Double -> Signal -> Signal
conj b (Vec a) = fmap (\y -> b * conjugate y) (Vec a)

twidd n k = exp( imagV $ k * pure (-2 * pi / fromIntegral n) )
twiddI n k = exp( imagV $ k * pure (2 * pi / fromIntegral n) )

trace_fft m x = trace (m ++ ": " ++ (show $ rd 2 $ absolute x)) x -- debugging helper to print fft results
fft x
    | n <= 16 = dft x
    | n <= 1 = x
    | otherwise =
        let (even, odd) = split $ runVec x
            (e, o) = (fft $ Vec even, fft $ Vec odd)
            k = Vec [0..]
            p = twidd n k
        in
            (e + p * o) <> (e - p * o)

    where
        n = length x
        split [] = ([], [])
        split [a] = ([a], [])
        split (a:b:c) = (a:x, b:y)
            where (x,y) = split c

ifft x
    | n <= 16 = idft x
    | n <= 1 = x
    | otherwise =
        let (even, odd) = split $ runVec x
            (e, o) = (fft $ Vec even, fft $ Vec odd)
            -- k = Vec [0..]
            bigN = fromIntegral n
            cE = fmap conjugate e
            cO = fmap conjugate o
            diC = (-1.0 / bigN) :+ 0
            p = twiddI n $ Vec [0..]
        in
            conj diC (e + p * o) <> conj diC (e - p * o)

    where
        n = length x
        split [] = ([], [])
        split [a] = ([a], [])
        split (a:b:c) = (a:x, b:y)
            where (x,y) = split c

low_pass :: Int -> Signal -> Signal

low_pass v sig =
    let
        m = mask v $ length sig
    in ifft $ m * fft sig

max_stars_repo_name: NikolaPeevski/Haskell-stuff

max_issues_repo_name: NikolaPeevski/Haskell-stuff

max_forks_repo_name: NikolaPeevski/Haskell-stuff


text:
module MachineLearning.LeastSquaresModelTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus

import MachineLearning.DataSets (dataset1)

import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning as ML
import MachineLearning.Optimization (checkGradient)
import MachineLearning.Model
import MachineLearning.LeastSquaresModel
import MachineLearning.Regularization (Regularization(..))

(x, y) = ML.splitToXY dataset1

x1 = ML.addBiasDimension x
initialTheta :: LA.Vector LA.R
initialTheta = LA.konst 1000 (LA.cols x1)
zeroTheta :: LA.Vector LA.R
zeroTheta = LA.konst 0 (LA.cols x1)

tests = [ testGroup "model" [
            testCase "cost, no reg"      $ assertApproxEqual "" 1e-5 1.6190245331702874e12 (cost LeastSquares RegNone x1 y initialTheta)
            , testCase "cost, lambda = 0"      $ assertApproxEqual "" 1e-5 1.6190245331702874e12 (cost LeastSquares (L2 0) x1 y initialTheta)
            , testCase "cost, lambda = 1"    $ assertApproxEqual "" 1e-5 1.619024554446883e12 (cost LeastSquares (L2 1) x1 y initialTheta)
            , testCase "cost, lambda = 1000" $ assertApproxEqual "" 1e-5 1.619045809766032e12 (cost LeastSquares (L2 1000) x1 y initialTheta)
            , testCase "gradient, no reg" $ assertVector "" 1e-5 gradient_l0 (gradient LeastSquares RegNone x1 y initialTheta)
            , testCase "gradient, lambda = 0" $ assertVector "" 1e-5 gradient_l0 (gradient LeastSquares (L2 0) x1 y initialTheta)
            , testCase "gradient, lambda = 1" $ assertVector "" 1e-5 gradient_l1 (gradient LeastSquares (L2 1) x1 y initialTheta)
            , testCase "gradient, lambda = 1000" $ assertVector "" 1e-5 gradient_l1000 (gradient LeastSquares (L2 1000) x1 y initialTheta)
            ]
          , testGroup "gradient checking" [
              testCase "non-zero theta, no reg" $ assertApproxEqual "" 1 0 (checkGradient LeastSquares RegNone x1 y initialTheta 1e-4)
              , testCase "non-zero theta, non-zero lambda" $ assertApproxEqual "" 1 0 (checkGradient LeastSquares (L2 2) x1 y initialTheta 1e-4)
              , testCase "zero theta, non-zero lambda" $ assertApproxEqual "" 1 0 (checkGradient LeastSquares (L2 2) x1 y zeroTheta 1e-4)
              , testCase "non-zero theta, zero lambda" $ assertApproxEqual "" 1 0 (checkGradient LeastSquares (L2 0) x1 y initialTheta 1e-4)
              , testCase "zero theta, zero lambda" $ assertApproxEqual "" 1 0(checkGradient LeastSquares (L2 0) x1 y zeroTheta 1e-4)
              ]
        ]

gradient_l0    = LA.vector [1664438.4042553192,3.865303999468085e9,5567440.808510638]
gradient_l1    = LA.vector [1664438.4042553192,3.865304020744681e9,5567462.085106383]
gradient_l1000 = LA.vector [1664438.4042553192,3.86532527606383e9, 5588717.404255319]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
import Data.Complex

main :: IO ()
main = do let res = dft (map (\x -> x :+ 0) [1..10])
       	  putStrLn $ show res

dft :: RealFloat a => [Complex a] -> [Complex a]
dft x = map (\(_, k) -> singleDFT k idxList) idxList
  where idxList = zip x [0..]

singleDFT :: RealFloat a => Int -> [(Complex a, Int)] -> Complex a
singleDFT k x = sum $ map (\(e, n) -> (e *) $ fCoeff k n $ length x) $ x

fCoeff :: RealFloat a => Int -> Int -> Int -> Complex a
fCoeff k n n_tot = exp (0.0 :+ ((-2.0) * pi * (fromIntegral $ k * n) / (fromIntegral n_tot)))

max_stars_repo_name: fredmorcos/attic

max_issues_repo_name: fredmorcos/attic

max_forks_repo_name: fredmorcos/attic


text:
{-|
Module: MachineLearning.NeuralNetwork.MultiSvmLoss
Description: Multi SVM Loss.
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Multi SVM Loss.
-}

module MachineLearning.NeuralNetwork.MultiSvmLoss
(
  scores
  , gradient
  , loss
)

where

import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA
import MachineLearning.Types (R, Matrix)
import MachineLearning.Utils (sumByRows, reduceByRows)

-- | SVM Delta
svmD :: R
svmD = 1.0

scores = id

gradient :: Matrix -> Matrix -> Matrix
gradient scores y =
    let correct_scores = sumByRows $ scores*(LA.step y)
        margins = scores - (correct_scores - (LA.scalar svmD))
        margins' = (1-y)*(LA.step margins)
        k = sumByRows margins'
    in margins' - (y * k)

loss :: Matrix -> Matrix -> R
loss scores y =
  let nSamples = fromIntegral $ LA.rows scores
      correct_scores = sumByRows $ scores*(LA.step y)
      margins = scores - (correct_scores - (LA.scalar svmD))
      margins' = margins * (LA.step margins)
      loss = LA.sumElements margins' - nSamples * svmD
  in loss / nSamples


text:
--
-- Benchmarks for LUSolve.
--

module LUSolveBenchmark (
    benchmarks
    ) where

import           Criterion                     (Benchmark, bench, bgroup, env,
                                                nf)

import           Numeric.LinearAlgebra.LUSolve (luFactor)

import qualified Data.Matrix.Generic           as M
import qualified Data.Vector.Unboxed           as V
import           System.Random

bundle :: Int -> [ a ] -> [[ a ]]
bundle _ [] = []
bundle n xs = take n xs : bundle n (drop n xs)

type Mat = M.Matrix V.Vector Double

runLUFactor :: Mat -> Mat
runLUFactor = (\(x, _, _) -> x) . luFactor

setupEnv :: IO (Mat, Mat, Mat)
setupEnv = do
  let mVals = randoms (mkStdGen 1) -- not a top level CAF, so will be GC'd promptly
      randomSquareMatrices :: Int -> [ Mat ]
      randomSquareMatrices n = Prelude.map (\vs -> M.fromLists (bundle n vs)) (bundle (n * n) mVals)
      m100:_ = randomSquareMatrices 100
      m500:_ = randomSquareMatrices 500
      m1000:_= randomSquareMatrices 1000
  return (m100, m500, m1000)

benchmarks :: Benchmark
benchmarks = env setupEnv $ \ ~(m100, m500, m1000) ->
  bgroup "luFactor"
    [ bench "100 x 100"   $ nf runLUFactor m100
    , bench "500 x 500"   $ nf runLUFactor m500
    , bench "1000 x 1000" $ nf runLUFactor m1000
    ]


text:
{- | Test of learning

In this example, two networks are used : 'simple' which is the reference and 'wrong' which is a wrong start.
The goal is to use test patterns to learn the right 'simple' network from 'wrong'. Only the values are learnt.
The topology of both networks is the same.

@
simple :: (['TDV' Bool],'SBN' 'CPT')
simple = 'runBN' $ do
    a <- 'variable' \"a\" ('t' :: Bool)
    b <- 'variable' \"b\" ('t' :: Bool)
--
    'proba' a '~~' [0.4,0.6]
    'cpt' b [a] '~~' [0.8,0.2,0.2,0.8]
--
    return [a,b]
@

and 'wrong' where the probability for a is wrong.

@
wrong :: (['TDV' Bool],'SBN' 'CPT')
wrong = 'runBN' $ do
    a <- 'variable' \"a\" ('t' :: Bool)
    b <- 'variable' \"b\" ('t' :: Bool)
--
    'proba' a '~~' [0.2,0.8]
    'cpt' b [a] '~~' [0.8,0.2,0.2,0.8]
--
    return [a,b]
@

So, the first thing to do is generate test patterns. We are using the 'discreteAncestralSampler' for this. This function is
generating a sequence of graphs. We are just interested in the values. So, we get the values with 'allVertexValues'.

@
generatePatterns :: IO [[DVI]]
generatePatterns = do
    let (vars\@[a,b],exampleG) = simple
    r <- 'runSampling' 5000 0 ('discreteAncestralSampler' exampleG)
    return (map 'allVertexValues' r)
@

Once we have the data, we can try to learn the network:

@
emTest = do
  samples <- generatePatterns
  let (_,simpleG) = simple
      (_,wrongG) = wrong
  print simpleG
  'printGraphValues' simpleG
  'printGraphValues' wrongG
--
  'printGraphValues' ('learnEM' samples wrongG)
@

First, we display the topology of the graph and the values for the reference graph and the wrong one.
Then, we use the 'learnEM' function to learn a new network from the samples. And, we print the new values to check.

-}
module Bayes.Examples.EMTest(
    -- * Test function
    emTest
  ) where

import Bayes.Sampling
import System.Random.MWC.CondensedTable
import qualified Data.Vector as V
import Bayes.Factor
import Bayes
import Bayes.FactorElimination
import Data.Maybe(fromJust)
import Bayes.BayesianNetwork
import Bayes.Factor.CPT
import Statistics.Sample.Histogram
import Bayes.EM

simple :: ([TDV Bool],SBN CPT)
simple = runBN $ do
    a <- variable "a" (t :: Bool)
    b <- variable "b" (t :: Bool)

    proba a ~~ [0.4,0.6]
    cpt b [a] ~~ [0.8,0.2,0.2,0.8]

    return [a,b]

wrong :: ([TDV Bool],SBN CPT)
wrong = runBN $ do
    a <- variable "a" (t :: Bool)
    b <- variable "b" (t :: Bool)

    proba a ~~ [0.2,0.8]
    cpt b [a] ~~ [0.8,0.2,0.2,0.8]

    return [a,b]

setJust vertex v = Just v

generatePatterns :: IO [[DVI]]
generatePatterns = do
    let (vars@[a,b],exampleG) = simple
    r <- runSampling 5000 0 (discreteAncestralSampler exampleG)
    return (map allVertexValues r)

emTest = do
  samples <- generatePatterns
  let (_,simpleG) = simple
      (_,wrongG) = wrong
  print simpleG
  printGraphValues simpleG
  printGraphValues wrongG

  printGraphValues (learnEM samples wrongG)

max_stars_repo_name: sid-kap/hbayes

max_issues_repo_name: sid-kap/hbayes

max_forks_repo_name: sid-kap/hbayes


text:
module Syntax where

import qualified Data.Complex as ComplexData
import qualified Data.Ord as OrdData

type Identifier = String
type Package = Identifier
type SliceStart = Integer
type SliceEnd = Integer
type SliceStep = Integer
type Ellipsis = ()
type ArrayLength = Integer
type ReadChannel = ()
type WriteChannel = ()
type Receiver = [FunctionParameterDecl]
type ImportDecl = [ImportSpec]

data QualifiedIdentifier = Qualified Package Identifier
  deriving (Eq, Ord, Show)

data Declaration
  = ConstDecl [Identifier] (Maybe Type) [Expression]
  | TypeDecl Identifier Type
  | VarDecl [Identifier] (Maybe Type) [Expression]
  deriving (Eq, Ord, Show)

data TopLevelDecl
  = Declaration Declaration
  | FunctionDecl Identifier FunctionSignature (Maybe [Statement])
  | ExternFunctionDecl Identifier FunctionSignature    -- this is not part of go RFC
  | MethodDecl Receiver Identifier FunctionSignature (Maybe [Statement])
  deriving (Eq, Ord, Show)

data Expression
  = Unary UnaryExpr
  | Binary BinaryOp Expression Expression
  deriving (Eq, Ord, Show)

data UnaryExpr
  = PrimaryExpr PrimaryExpr
  -- FIXME: had to change this from 'UnaryExpr' to 'Expression'
  -- don't know if it will cause any problems
  | UnaryOpExpr UnaryOp Expression
  deriving (Eq, Ord, Show)

data PrimaryExpr
  = Operand Operand
  | Conversion Conversion
  | Selector PrimaryExpr Identifier
  | Index PrimaryExpr Expression
  | Slice PrimaryExpr (Maybe SliceStart) (Maybe SliceEnd) (Maybe SliceStep)
  | TypeAssertion PrimaryExpr Type
  | Arguments PrimaryExpr [Expression] (Maybe Type) (Maybe Ellipsis)
  deriving (Eq, Ord, Show)

data BinaryOp
  = Or        -- ||
  | And       -- &&
  | RelBinOp RelOp
  | AddBinOp AddOp
  | MulBinOp MulOp
  deriving (Eq, Ord, Show)

data RelOp
  = Equals    -- ==
  | NotEquals -- !=
  | Less      -- <
  | LessEq    -- <=
  | Greater   -- >
  | GreaterEq -- >=
  deriving (Eq, Ord, Show)

data AddOp
  = Plus      -- +
  | Minus     -- -
  | BitOr     -- |
  | BitXor    -- ^
  deriving (Eq, Ord, Show)

data MulOp
  = Mult      -- *
  | Div       -- /
  | Mod       -- %
  | LShift    -- <<
  | RShift    -- >>
  | BitAnd    -- &
  | AndXor    -- &^ - FIXME: no idea what is that
  deriving (Eq, Ord, Show)

data UnaryOp
  = UnPlus    -- +
  | UnMinus   -- -
  | Not       -- !
  | UnCaret   -- ^
  | Deref     -- *
  | Pointerof -- &
  | ChanRead  -- <-
  deriving (Eq, Ord, Show)

data TypeName
  = Identifier Identifier
  | QualifiedIdentifier QualifiedIdentifier
  deriving (Eq, Ord, Show)

data TypeWithEllipsis = TypeWithEllipsis Type (Maybe Ellipsis)

data Operand
  = Literal Literal
  | OperandName OperandName
  | MethodExpr Type Identifier
  | ExpressionInParens Expression
  deriving (Eq, Ord, Show)

data OperandName
  = OperandIdentifier Identifier
  | OperandQualified QualifiedIdentifier
  deriving (Eq, Ord, Show)

data Conversion = TypeConversion Type Expression
  deriving (Eq, Ord, Show)

data Type
  = TypeName TypeName
  | TypeLiteral TypeLiteral
  deriving (Eq, Ord, Show)

data TypeLiteral
  = ArrayType ArrayLength Type
  | StructType [FieldDecl]
  | PointerType Type
  | FunctionType FunctionSignature
  | InterfaceType [MethodSpec]
  | SliceType Type
  | MapType Type Type
  | ChannelType Type (Maybe ReadChannel) (Maybe WriteChannel)
  deriving (Eq, Ord, Show)

data FieldDecl
  = Field [Identifier] Type
  | AnonymousField TypeName
  deriving (Eq, Ord, Show)

data MethodSpec
  = Method Identifier FunctionSignature
  | Interface TypeName
  deriving (Eq, Ord, Show)

data FunctionSignature = Signature [FunctionParameterDecl] (Maybe FunctionResult)
  deriving (Eq, Ord, Show)

data FunctionParameterDecl = ParameterDecl [Identifier] Type (Maybe Ellipsis)
  deriving (Eq, Ord, Show)

data FunctionResult
  = FunctionTupleResult [FunctionParameterDecl]
  | FunctionSingleResult Type
  deriving (Eq, Ord, Show)

data Statement
  = DeclarationStmt Declaration
  | LabeledStmt Identifier Statement
  | SimpleStmt SimpleStmt
  | GoStmt Expression
  | ReturnStmt [Expression]
  | BreakStmt (Maybe Identifier)
  | ContinueStmt (Maybe Identifier)
  | GotoStmt Identifier
  | FallthroughStmt
  | Block [Statement]
  | IfStmt IfStmt
  | SwitchStmt SwitchStmt
  | SelectStmt [CommClauseStmt]
  | ForStmt ForStmt
  | DeferStmt Expression
  deriving (Eq, Ord, Show)

data SimpleStmt
  = Empty
  | ExpressionStmt Expression
  | SendStmt SendStmt
  | IncStmt Expression
  | DecStmt Expression
  | AssignmentStmt [Expression] AssignOp [Expression]
  | ShortVarDeclStmt [Identifier] [Expression]
  deriving (Eq, Ord, Show)

data IfStmt = If (Maybe SimpleStmt) Expression [Statement] (Maybe ElseStmt)
  deriving (Eq, Ord, Show)

data ElseStmt
  = ElseIf IfStmt
  | Else [Statement]
  deriving (Eq, Ord, Show)

data SwitchStmt
  = ExprSwitchStmt (Maybe SimpleStmt) (Maybe Expression) [ExprCaseClauseStmt]
  | TypeSwitchStmt (Maybe SimpleStmt) TypeSwitchGuardStmt [TypeCaseClauseStmt]
  deriving (Eq, Ord, Show)

data ExprCaseClauseStmt
  = ExprCaseClause [Expression] [Statement]
  | DefaultCaseClause [Statement]
  deriving (Eq, Ord, Show)

data TypeSwitchGuardStmt = TypeSwitchGuard (Maybe Identifier) PrimaryExpr
  deriving (Eq, Ord, Show)

data TypeCaseClauseStmt
  = TypeCaseClause [Type] [Statement]
  | DefaultTypeCaseClause [Statement]
  deriving (Eq, Ord, Show)

data CommClauseStmt = CommClause CommCaseStmt [Statement]
  deriving (Eq, Ord, Show)

data CommCaseStmt
  = SendCommClause SendStmt
  | RecvCommClause RecvStmt
  | DefaultCommClause
  deriving (Eq, Ord, Show)

data SendStmt = Send Expression Expression
  deriving (Eq, Ord, Show)

data RecvStmt = Recv RecvTarget Expression
  deriving (Eq, Ord, Show)

data RecvTarget
  = RecvAssignment [Expression]
  | RecvShortVarDecl [Identifier]
  deriving (Eq, Ord, Show)

data ForStmt = For ForDefinition [Statement]
  deriving (Eq, Ord, Show)

data ForDefinition
  = Condition Expression
  | ForClause (Maybe SimpleStmt) (Maybe Expression) (Maybe SimpleStmt)
  | RangeClause RangeTarget Expression
  deriving (Eq, Ord, Show)

data RangeTarget
  = RangeAssignment [Expression]
  | RangeShortVarDecl [Identifier]
  deriving (Eq, Ord, Show)

data AssignOp
  = AssignAddOp AddOp
  | AssignMulOp MulOp
  | Assign
  deriving (Eq, Ord, Show)

data Literal
  = IntLiteral Integer
  | FloatLiteral Double
  | ImaginaryLiteral ComplexValue
  | RuneLiteral Char
  | StringLiteral String
  | CompositeLiteral LiteralType LiteralValue
  | FunctionLiteral FunctionSignature [Statement]
  deriving (Eq, Ord, Show)

data LiteralType
  = StructTypeLiteral [FieldDecl]
  | ArrayTypeLiteral ArrayLength Type
  | AutoSizeArrayTypeLiteral Type
  | SliceTypeLiteral Type
  | MapTypeLiteral Type Type
  | TypeNameLiteral TypeName
  deriving (Eq, Ord, Show)

data LiteralValue = LiteralValue [LiteralElement]
  deriving (Eq, Ord, Show)

data LiteralElement = LiteralElement (Maybe LiteralKey) Value
  deriving (Eq, Ord, Show)

data LiteralKey
  = FieldName Identifier
  | KeyExpr Expression
  | KeyLiteralValue LiteralValue
  deriving (Eq, Ord, Show)

data Value
  = ValueExpr Expression
  | ValueLiteralValue LiteralValue
  deriving (Eq, Ord, Show)

data SourceFile = SourceFile Identifier [ImportDecl] [TopLevelDecl]
  deriving (Eq, Ord, Show)

data ImportSpec = ImportSpec Identifier String
  deriving (Eq, Ord, Show)

type ComplexValue = ComplexData.Complex Double

instance (Eq t, Ord t, RealFloat t) => Ord (ComplexData.Complex t) where
  compare a b =
    case byRealPart of
      OrdData.EQ -> byImagPart
      _ -> byRealPart
    where
      byRealPart = compare (ComplexData.realPart a) (ComplexData.realPart b)
      byImagPart = compare (ComplexData.imagPart a) (ComplexData.imagPart b)

max_stars_repo_name: waterlink/hgp

max_issues_repo_name: waterlink/hgp

max_forks_repo_name: waterlink/hgp


text:
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Numeric.Loss.MeanSquaredError where

import Numeric.Loss
import Dependent.Size
import qualified Numeric.Vector.Sized as Sized
import Numeric.Vector.Sized(SizedVector',SizedVector)
import Data.Proxy
import GHC.TypeLits(KnownNat,Nat,natVal)

newtype MeanSquaredError (n :: Nat) = MeanSquaredError ()

instance forall n. KnownNat n => Loss (ZZ ::. n) (MeanSquaredError n) where
  -- forward (MeanSquaredError targets) x
  --   = (err, x)
  --     where
  --       err = (normalize . Sized.foldAll (+) 0 . Sized.map square) (Sized.zipWith (-) x (Sized.use targets))
  --       square x = x ^2
  --       n = fromInteger $ natVal (Proxy :: Proxy (n :: Nat))
  --       normalize total = Sized.map (\t -> t / n) total

  lossDerivative _ x targets = dx
    where
      dx = normalize (Sized.zipWith (-) x targets)
      n = fromInteger $ natVal (Proxy :: Proxy (n :: Nat))
      normalize zs = Sized.map (\z -> 2/n * z) zs

max_stars_repo_name: mixed-signals/mixed-signals

max_issues_repo_name: mixed-signals/mixed-signals

max_forks_repo_name: mixed-signals/mixed-signals


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
-- | Goodness of fit runs tests.
module Statistics.Test.Runs (
    runsTest
  , module Statistics.Test.Types
  ) where

import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import Numeric.SpecFunctions (choose)
import Prelude hiding (sum)

import Statistics.Sample.Internal (sum)
import Statistics.Test.Types
import Statistics.Types

-- | Goodness of fit test for binned data. It uses only sign of
--   deviations of observations from their expectations. Null
--   hypothesis is that all possible patterns of sign occurrences are
--   equiprobable.
--
--   It's asymptotically independent from chi-square test. So their
--   results could be directly combined
runsTest :: (G.Vector v Bool) => v Bool -> Test ()
{-# INLINE runsTest #-}
runsTest v
  = Test { testSignificance = mkPValue $ cumulativeProb n m r
         , testStatistics   = fromIntegral r
         , testDistribution = ()
         }
  where
    (n,m,r) = computeRuns v

-- Compute number of positive elements, negative elements and runs
computeRuns :: (G.Vector v Bool) => v Bool -> (Int,Int,Int)
{-# INLINE computeRuns #-}
computeRuns v
  = fini $ G.foldl' step (0,0,0,Nothing) v
  where
    step (!nP,!nM,!nR,!old) f =
      ( if f then nP+1 else nP
      , if f then nM   else nM+1
      , if old == Just f then nR else nR+1
      , Just f
      )
    fini (nP,nM,nR,_) = (nP,nM,nR)

-- Compute denormalized probability of getting R runs given N positive
-- and M positive elements
denormProbability :: Int -> Int -> Int -> Double
denormProbability n m r
  | even r    = 2 * ((m-1) `choose` (s-1)) * ((n-1) `choose` (s-1))
  | otherwise = ((m-1) `choose` (s-1)) * ((n-1) `choose` (s-2))
              + ((m-1) `choose` (s-2)) * ((n-1) `choose` (s-1))
  where
    s = r `quot` 2

-- Probability of getting R<=R[observed]
cumulativeProb :: Int -> Int -> Int -> Double
cumulativeProb n m r
  = min 1
  $ sum (U.map (denormProbability n m) $ U.enumFromTo 1 r)
  / ((n+m) `choose` m)


text:
{-| This module defines the 'Consume' typeclass, used for incrementally
    destructing inputs to random non-strict functions.

    Calling 'consume' on some value lazily returns an abstract type of 'Input',
    which contains all the entropy present in the original value. Paired with
    'Test.StrictCheck.Produce', these @Input@ values can be used to generate
    random non-strict functions, whose strictness behavior is dependent on the
    values given to them.
-}
module Test.StrictCheck.Consume
  ( -- * Incrementally consuming input
    Input
  , Inputs
  , Consume(..)
  -- * Manually writing 'Consume' instances
  , constructor
  , normalize
  , consumeTrivial
  , consumePrimitive
  -- * Generically deriving 'Consume' instances
  , GConsume
  , gConsume
  ) where

import Test.QuickCheck
import Generics.SOP
import Generics.SOP.NS

import Test.StrictCheck.Internal.Inputs

import Data.Complex

import Data.Foldable as Fold
import Data.List.NonEmpty (NonEmpty(..))
import Data.Tree     as Tree
import Data.Set      as Set
import Data.Map      as Map
import Data.Sequence as Seq
import Data.IntMap   as IntMap
import Data.IntSet   as IntSet

-- | Lazily monomorphize some input value, by converting it into an @Input@.
-- This is an incremental version of QuickCheck's @CoArbitrary@ typeclass.
-- It can also be seen as a generalization of the @NFData@ class.
--
-- Instances of @Consume@ can be derived automatically for any type implementing
-- the @Generic@ class from "GHC.Generics". Using the @DeriveAnyClass@
-- extension, we can say:
--
-- > import GHC.Generics as GHC
-- > import Generics.SOP as SOP
-- >
-- > data D x y
-- >   = A
-- >   | B (x, y)
-- >   deriving (GHC.Generic, SOP.Generic, Consume)
--
-- This automatic derivation follows these rules, which you can follow too if
-- you're manually writing an instance for some type which is not @Generic@:
--
-- For each distinct constructor, make a single call to 'constructor' with
-- a distinct @Int@, and a list of @Input@s, each created by recursively calling
-- 'consume' on every field in that constructor. For abstract types (e.g. sets),
-- the same procedure can be used upon an extracted list representation of the
-- contents.
class Consume a where
  -- | Convert an @a@ into an @Input@ by recursively destructing it using calls
  -- to @consume@
  consume :: a -> Input
  default consume :: GConsume a => a -> Input
  consume = gConsume

-- | Reassemble pieces of input into a larger Input: this is to be called on the
-- result of @consume@-ing subparts of input
constructor :: Int -> [Input] -> Input
constructor n !is =
  Input (Variant (variant n)) is

-- | Use the CoArbitrary instance for a type to consume it
--
-- This should only be used for "flat" types, i.e. those which contain no
-- interesting consumable substructure, as it's fully strict (non-incremental)
consumePrimitive :: CoArbitrary a => a -> Input
consumePrimitive !a =
  Input (Variant (coarbitrary a)) []

-- | Consume a type which has no observable structure whatsoever
--
-- This should only be used for types for which there is only one inhabitant, or
-- for which inhabitants cannot be distinguished at all.
consumeTrivial :: a -> Input
consumeTrivial !_ =
  Input mempty []

-- | Fully normalize something which can be consumed
normalize :: Consume a => a -> ()
normalize (consume -> input) = go input
  where
    go (Input _ is) = Fold.foldr seq () (fmap go is)

--------------------------------------------
-- Deriving Consume instances generically --
--------------------------------------------

-- | The constraints necessary to generically @consume@ something
type GConsume a = (Generic a, All2 Consume (Code a))

-- | Generic 'consume'
gConsume :: GConsume a => a -> Input
gConsume !(from -> sop) =
  constructor (index_SOP sop)
  . hcollapse
  . hcliftA (Proxy @Consume) (K . consume . unI)
  $ sop

---------------
-- Instances --
---------------

instance Consume (a -> b)  where consume = consumeTrivial
instance Consume (Proxy p) where consume = consumeTrivial

instance Consume Char     where consume = consumePrimitive
instance Consume Word     where consume = consumePrimitive
instance Consume Int      where consume = consumePrimitive
instance Consume Double   where consume = consumePrimitive
instance Consume Float    where consume = consumePrimitive
instance Consume Rational where consume = consumePrimitive
instance Consume Integer  where consume = consumePrimitive
instance (CoArbitrary a, RealFloat a) => Consume (Complex a) where
  consume = consumePrimitive

instance Consume ()
instance Consume Bool
instance Consume Ordering
instance Consume a => Consume (Maybe a)
instance (Consume a, Consume b) => Consume (Either a b)
instance Consume a => Consume [a]

instance Consume a => Consume (NonEmpty a) where
  consume (a :| as) = constructor 0 [consume a, consume as]

instance Consume a => Consume (Tree a) where
  consume (Node a as) = constructor 0 [consume a, consume as]

instance Consume v => Consume (Map k v) where
  consume = constructor 0 . fmap (consume . snd) . Map.toList

consumeContainer :: (Consume a, Foldable t) => t a -> Input
consumeContainer = constructor 0 . fmap consume . Fold.toList

instance Consume v => Consume (Seq v)    where consume = consumeContainer
instance Consume v => Consume (Set v)    where consume = consumeContainer
instance Consume v => Consume (IntMap v) where consume = consumeContainer
instance Consume IntSet where
  consume = consumeContainer . IntSet.toList

-- TODO: instances for the rest of Containers

instance (Consume a, Consume b) => Consume (a, b)
instance (Consume a, Consume b, Consume c) => Consume (a, b, c)
instance (Consume a, Consume b, Consume c, Consume d) => Consume (a, b, c, d)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e
         ) => Consume
  (a, b, c, d, e)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         ) => Consume
  (a, b, c, d, e, f)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g
         ) => Consume
  (a, b, c, d, e, f, g)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h
         ) => Consume
  (a, b, c, d, e, f, g, h)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i
         ) => Consume
  (a, b, c, d, e, f, g, h, i)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v, Consume w
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v, Consume w, Consume x
          ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v, Consume w, Consume x
         , Consume y
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v, Consume w, Consume x
         , Consume y, Consume z
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)

max_stars_repo_name: janpath/StrictCheck

max_issues_repo_name: janpath/StrictCheck

max_forks_repo_name: janpath/StrictCheck


text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-|
Module      : Grenade.Layers.Concat
Description : Concatenation layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental

This module provides the concatenation layer, which runs two chilld layers in parallel and combines their outputs.
-}
module Grenade.Layers.Concat
  ( Concat(..)
  -- , SpecConcat(..)
  ) where

import           Data.Serialize

import           Control.DeepSeq              (NFData (..))
import           Data.Singletons
import           GHC.Generics                 (Generic)
import           GHC.TypeLits
import           Numeric.LinearAlgebra.Static (R, row, split, splitRows, unrow, ( # ),
                                               (===))
#if MIN_VERSION_base(4,9,0)
import           Data.Kind                    (Type)
#endif

import           Grenade.Core
import           Grenade.Onnx

-- | A Concatenating Layer.
--
-- This layer shares it's input state between two sublayers, and concatenates their output.
--
-- With Networks able to be Layers, this allows for very expressive composition of complex Networks.
--
-- The Concat layer has a few instances, which allow one to flexibly "bash" together the outputs.
--
-- Two 1D vectors, can go to a 2D shape with 2 rows if their lengths are identical.
-- Any 2 1D vectors can also become a longer 1D Vector.
--
-- 3D images become 3D images with more channels. The sizes must be the same, one can use Pad
-- and Crop layers to ensure this is the case.
data Concat :: Shape -> Type -> Shape -> Type -> Type where
  Concat :: !x -> !y -> Concat m x n y
  deriving (Generic, NFData)

instance (Show x, Show y) => Show (Concat m x n y) where
  show (Concat x y) = "Concat\n" ++ show x ++ "\n" ++ show y

-- | Run two layers in parallel, combining their outputs.
instance (UpdateLayer x, UpdateLayer y) => UpdateLayer (Concat m x n y) where
  type Gradient (Concat m x n y) = (Gradient x, Gradient y)
  runUpdate opt (Concat x y) (x', y') = Concat (runUpdate opt x x') (runUpdate opt y y')
  reduceGradient grads = (reduceGradient @x xs, reduceGradient @y ys)
    where
      (xs, ys) = unzip grads

instance (RandomLayer x, RandomLayer y) => RandomLayer (Concat m x n y) where
  createRandomWith m gen = Concat <$> createRandomWith m gen <*> createRandomWith m gen

instance ( SingI i
         , Layer x i ('D1 o)
         , Layer y i ('D1 o)
         ) => Layer (Concat ('D1 o) x ('D1 o) y) i ('D2 2 o) where
  type Tape (Concat ('D1 o) x ('D1 o) y) i ('D2 2 o) = (Tape x i ('D1 o), Tape y i ('D1 o))

  runForwards (Concat x y) input =
    let (xT, xOut :: S ('D1 o)) = runForwards x input
        (yT, yOut :: S ('D1 o)) = runForwards y input
    in case (xOut, yOut) of
        (S1D xOut', S1D yOut') ->
            ((xT, yT), S2D (row xOut' === row yOut'))

  runBackwards (Concat x y) (xTape, yTape) (S2D o) =
    let (ox, oy) = splitRows o
        (x', xB :: S i) = runBackwards x xTape (S1D $ unrow ox)
        (y', yB :: S i) = runBackwards y yTape (S1D $ unrow oy)
    in  ((x', y'), xB + yB)

instance ( SingI i
         , Layer x i ('D1 m)
         , Layer y i ('D1 n)
         , KnownNat o
         , KnownNat m
         , KnownNat n
         , o ~ (m + n)
         , n ~ (o - m)
         , (m <=? o) ~ 'True
         ) => Layer (Concat ('D1 m) x ('D1 n) y) i ('D1 o) where
  type Tape (Concat ('D1 m) x ('D1 n) y) i ('D1 o) = (Tape x i ('D1 m), Tape y i ('D1 n))

  runForwards (Concat x y) input =
    let (xT, xOut :: S ('D1 m)) = runForwards x input
        (yT, yOut :: S ('D1 n)) = runForwards y input
    in case (xOut, yOut) of
        (S1D xOut', S1D yOut') ->
            ((xT, yT), S1D (xOut' # yOut'))

  runBackwards (Concat x y) (xTape, yTape) (S1D o) =
    let (ox :: R m , oy :: R n) = split o
        (x', xB :: S i) = runBackwards x xTape (S1D ox)
        (y', yB :: S i) = runBackwards y yTape (S1D oy)
    in  ((x', y'), xB + yB)

-- | Concat 3D shapes, increasing the number of channels.
instance ( SingI i
         , Layer x i ('D3 rows cols m)
         , Layer y i ('D3 rows cols n)
         , KnownNat (rows * n)
         , KnownNat (rows * m)
         , KnownNat (rows * o)
         , KnownNat o
         , KnownNat m
         , KnownNat n
         , ((rows * m) + (rows * n)) ~ (rows * o)
         , ((rows * o) - (rows * m)) ~ (rows * n)
         , ((rows * m) <=? (rows * o)) ~ 'True
         ) => Layer (Concat ('D3 rows cols m) x ('D3 rows cols n) y) i ('D3 rows cols o) where
  type Tape (Concat ('D3 rows cols m) x ('D3 rows cols n) y) i ('D3 rows cols o) = (Tape x i ('D3 rows cols m), Tape y i ('D3 rows cols n))

  runForwards (Concat x y) input =
    let (xT, xOut :: S ('D3 rows cols m)) = runForwards x input
        (yT, yOut :: S ('D3 rows cols n)) = runForwards y input
    in case (xOut, yOut) of
        (S3D xOut', S3D yOut') ->
            ((xT, yT), S3D (xOut' === yOut'))

  runBackwards (Concat x y) (xTape, yTape) (S3D o) =
    let (ox, oy) = splitRows o
        (x', xB :: S i) = runBackwards x xTape (S3D ox :: S ('D3 rows cols m))
        (y', yB :: S i) = runBackwards y yTape (S3D oy :: S ('D3 rows cols n))
    in  ((x', y'), xB + yB)

instance (Serialize a, Serialize b) => Serialize (Concat sa a sb b) where
  put (Concat a b) = put a *> put b
  get = Concat <$> get <*> get

-------------------- OnnxLoadableParallel instance --------------------
instance OnnxLoadableParallel (Concat m x n y) x y where
  mkParallelLayer = Concat

instance OnnxOperator (Concat m x n y) where
  onnxOpTypeNames _ = ["Concat"]

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE ImplicitParams      #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications    #-}

module Main ( main ) where

import           Data.Int
import           Data.List (sort)
import           Data.List.NonEmpty (NonEmpty(..))
import           Data.Massiv.Array (getComp, index')
import qualified Data.Set as S
import qualified Data.Vector as V
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Unboxed as U
import           Data.Word
import           Geography.MapAlgebra
import qualified Numeric.LinearAlgebra as LA
import           Prelude as P
import           Test.HUnit.Approx
import qualified Test.QuickCheck.Arbitrary as QC
import           Test.Tasty
import           Test.Tasty.HUnit
import           Test.Tasty.QuickCheck

---

main :: IO ()
main = do
  img <- fromRight <$> fromGray "data/gray512.tif"
  defaultMain $ suite img

suite :: Raster S p 512 512 Word8 -> TestTree
suite r = testGroup "Unit Tests"
  [ testGroup "Raster Creation"
    [ testCase "constant (256x256)"     $ length (lazy small) @?= 65536
    , testCase "constant (2^16 x 2^16)" $ length lazybig @?= 4294967296
    , testCase "Image Reading (RGBA)"   $ do
        i <- fileRGBA
        fmap (getComp . _array . _red) i @?= Right Par
    ]
  , testGroup "Typeclass Ops"
    [ testCase "(==)" $ assertBool "(==) doesn't work" (small == small)
    , testCase "(+)"  $ strict P (lazy one + lazy one) @?= two
    ]
  , testGroup "Folds"
    [ testCase "sum (small)" $ P.sum (lazy small) @?= 327680
    -- , testCase "sum (large)" $ P.sum lazybig @?= 21474836480
    ]
  , testGroup "Local Ops"
    [ testCase "(+)"       $ P.sum (lazy small + lazy small) @?= (327680 * 2)
    , testCase "lmin"      $ strict P (lmin one two) @?= one
    , testCase "lvariety"  $ (strict P . lvariety . fmap lazy $ one :| [two]) @?= two
    , testCase "lmajority" $ (strict P . lmajority . fmap lazy $ one :| [one, two]) @?= one
    , testCase "lminority" $ (strict P . lminority . fmap lazy $ one :| [one, two]) @?= two
    -- , testCase "(+) big"   $ strict P (lazy big + lazy big) @?= bog
    ]
  , testGroup "Focal Ops"
    [ testCase "fvariety" $ strict P (fvariety one) @?= one
    , testCase "fmax"     $ strict P (fmax one) @?= one
    , testCase "fmin"     $ strict P (fmin one) @?= one
    , testGroup "flinkage"
      [ testCase "single point" singlePoint
      , testCase "2x2 same" twoByTwoSame
      , testCase "2x2 diff" twoByTwoDiff
      , testCase "3x3" threeByThree
      ]
    , testCase "flength" flengthTest
    , testCase "fpartition" fpartitionTest
    , testCase "fshape" fshapeTest
    , testCase "ffrontage" ffrontageTest
    , testGroup "farea"
      [ testCase "3x3 Open" fareaOpen
      , testCase "3x3 Centre" fareaCentre
      , testCase "4x4 Complex" fareaComplex
      ]
    , testGroup "fvolume"
      [ testCase "3x3 Flat" fvolumeFlat
      , testCase "3x3 Hill" fvolumeHill
      ]
    , testProperty "Least Squares" leastSquares
    , testGroup "fgradient"
      [ testCase "3x3 Flat" fgradientFlat
      , testCase "3x3 (tau/8)" fgradient45
      ]
    , testGroup "faspect"
      [ testCase "3x3 Flat" faspectFlat
      , testCase "3x3 East" faspectEast
      , testCase "3x3 South" faspect45
      ]
    , testGroup "fdownstream"
      [ testCase "3x3 Spikey" fdownstream4
      , testCase "3x3 Flat" fdownstreamFlat
      , testCase "3x3 Peak" fdownstreamPeak
      , testCase "3x3 Pit"  fdownstreamPit
      ]
    , testGroup "fupstream"
      [ testCase "3x3 Peak" fupstreamPeak
      , testCase "3x3 Flat" fupstreamFlat
      ]
    ]
  , testGroup "Histograms"
    [ testCase "Total Sum"     $ VS.sum (_histogram $ histogram r) @?= 262144
    , testCase "10 Breaks"     $ length (breaks $ histogram r) @?= 10
    , testCase "Sorted Breaks" $ do
        let bs = breaks $ histogram r
        sort bs @?= bs
    ]
  ]

fromRight :: Either a b -> b
fromRight (Right b) = b
fromRight _         = error "Was Left"

one :: Raster P p 7 7 Word
one = constant P Seq 1

two :: Raster P p 7 7 Word
two = constant P Seq 2

small :: Raster P p 256 256 Int
small = constant P Seq 5

lazybig :: Raster D p 65536 65536 Int
lazybig = constant D Par 5

-- big :: Raster P p 65536 65536 Word8
-- big = constant P Par 5

-- bog :: Raster P p 65536 65536 Word8
-- bog = constant P Par 10

-- indices :: Raster S p 512 512 Word8
-- indices = fromFunction S Par (\(r :. c) -> fromIntegral $ r + c)

fileRGBA :: IO (Either String (RGBARaster p 512 512 Word8))
fileRGBA = fromRGBA "data/512x512.tif"

singlePoint :: Assertion
singlePoint = actual @?= expected
  where expected :: Raster B p 1 1 Line
        expected = constant B Seq (Line 0)
        actual :: Raster B p 1 1 Line
        actual = strict B . flinkage $ constant P Seq (1 :: Int)

twoByTwoSame :: Assertion
twoByTwoSame = actual @?= expected
  where expected :: Raster S p 2 2 Line
        expected = fromRight . fromVector Seq . VS.fromList
          $ P.map (Line . _drain . drainage . S.fromList) [ [ East, South ]
                                                          , [ West, South ]
                                                          , [ North,East ]
                                                          , [ West, North ] ]
        actual :: Raster S p 2 2 Line
        actual = fromRight . fmap (strict S . flinkage) . fromVector Seq $ U.fromList ([1,1,1,1] :: [Int])

twoByTwoDiff :: Assertion
twoByTwoDiff = actual @?= expected
  where expected :: Raster S p 2 2 Line
        expected = fromRight . fromVector Seq . VS.fromList
          $ P.map (Line . _drain . drainage . S.fromList) [ [ SouthEast ]
                                                          , [ SouthWest ]
                                                          , [ NorthEast ]
                                                          , [ NorthWest ] ]
        actual :: Raster S p 2 2 Line
        actual = fromRight . fmap (strict S . flinkage) . fromVector Seq $ U.fromList ([1,2,2,1] :: [Int])

threeByThree :: Assertion
threeByThree = actual @?= expected
  where expected :: Raster S p 3 3 Line
        expected = fromRight . fromVector Seq . VS.fromList
          $ P.map (Line . _drain . drainage . S.fromList) [ [ ]
                                                          , [ South ]
                                                          , [ ]
                                                          , [ East ]
                                                          , [ North, West, South, East ]
                                                          , [ West ]
                                                          , [ ]
                                                          , [ North ]
                                                          , [ ] ]
        actual :: Raster S p 3 3 Line
        actual = fromRight . fmap (strict S . flinkage) . fromVector Seq $ U.fromList ([1,2,1,2,2,2,1,2,1] :: [Int])

flengthTest :: Assertion
flengthTest = actual @?= expected
  where actual :: Raster U p 3 3 Double
        actual = strict U . flength . flinkage . fromRight . fromVector Seq $ VS.fromList ([1,2,1,2,2,2,1,2,1] :: [Int])
        expected :: Raster U p 3 3 Double
        expected = fromRight . fromVector Seq $ U.fromList [ 0, 0.5, 0, 0.5, 2, 0.5, 0, 0.5, 0 ]

fpartitionTest :: Assertion
fpartitionTest = actual @?= expected
  where expected :: Raster B p 2 2 Corners
        expected = fromRight . fromVector Seq $ V.fromList [ Corners Open Open Open Open
                                                           , Corners Open Open Open Open
                                                           , Corners OneSide Open OneSide Complete
                                                           , Corners Open Open Open Open ]
        actual :: Raster B p 2 2 Corners
        actual = strict B . fpartition . fromRight . fromVector Seq $ U.fromList ([1,1,2,1] :: [Int])

fshapeTest :: Assertion
fshapeTest = actual @?= expected
 where expected :: Raster B p 3 3 Corners
       expected = fromRight . fromVector Seq $ V.fromList [ Corners Open Open OutFlow Open
                                                          , Corners Open Open Open Open
                                                          , Corners Open OutFlow Open Open
                                                          , Corners Open Open Open Open
                                                          , Corners Complete Complete Complete Complete
                                                          , Corners Open Open Open Open
                                                          , Corners Open Open Open OutFlow
                                                          , Corners Open Open Open Open
                                                          , Corners OutFlow Open Open Open ]
       actual :: Raster B p 3 3 Corners
       actual = strict B . fshape . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,0,1,1,1,1] :: [Int])

ffrontageTest :: Assertion
ffrontageTest = let ?epsilon = 0.001 in actual @?~ expected
  where expected :: Double
        expected = 1 + (1 / sqrt 2)
        actual :: Double
        actual = flip index' (1 :. 1) . _array . strict S $ ffrontage rast
        rast :: Raster DW p 4 4 Corners
        rast = fshape . fromRight . fromVector Seq $ U.fromList ( [1,1,1,0
                                                                  ,1,0,0,0
                                                                  ,1,0,0,1
                                                                  ,1,0,1,1] :: [Int] )

fareaOpen :: Assertion
fareaOpen = actual @?= expected
  where expected :: Raster U p 3 3 Double
        expected = fromRight . fromVector Seq $ U.fromList [1,1,1,1,1,1,1,1,1]
        actual :: Raster U p 3 3 Double
        actual = strict U . farea . fshape . fromRight . fromVector Seq $ U.fromList ([0,0,0,0,0,0,0,0,0] :: [Int])

fareaCentre :: Assertion
fareaCentre = actual @?= expected
  where expected :: Raster U p 3 3 Double
        expected = fromRight . fromVector Seq $ U.fromList [ 1 + 1/8, 1, 1 + 1/8
                                                           , 1, 1/2, 1
                                                           , 1 + 1/8, 1, 1 + 1/8 ]
        actual :: Raster U p 3 3 Double
        actual = strict U . farea . fshape . fromRight . fromVector Seq $ U.fromList ([0,0,0,0,1,0,0,0,0] :: [Int])

fareaComplex :: Assertion
fareaComplex = let ?epsilon = 0.001 in actual @?~ (7 / 8)
  where actual :: Double
        actual = flip index' (1 :. 1) . _array . strict P $ farea rast
        rast :: Raster DW p 4 4 Corners
        rast = fshape . fromRight . fromVector Seq $ U.fromList ( [1,1,1,0
                                                                  ,1,0,0,0
                                                                  ,1,0,0,1
                                                                  ,1,0,1,1] :: [Int] )

fvolumeFlat :: Assertion
fvolumeFlat = strict U (fvolume expected) @?= expected
  where expected :: Raster U p 3 3 Double
        expected = fromRight . fromVector Seq $ U.fromList [8,8,8,8,8,8,8,8,8]

fvolumeHill :: Assertion
fvolumeHill = index' (_array actual) (1 :. 1) @?= expected
  where expected :: Double
        expected = P.sum [20,20,16,20,16,16,16,16,12,16,12,12] / 12
        actual :: Raster U p 3 3 Double
        actual = strict U . fvolume @Double . fromRight . fromVector Seq $ U.fromList [24,24,24
                                                                                ,16,16,16
                                                                                ,8,8,8]

newtype Vec = Vec [Double] deriving (Show)

instance Arbitrary Vec where
  arbitrary = Vec <$> QC.vector 9

-- | A QuickCheck property to test whether my custom Least Squares is as
-- accurate as the one provided by HMatrix.
leastSquares :: Vec -> Bool
leastSquares (Vec vs) = f 0 && f 1 && f 2
  where m = head . LA.toColumns $ LA.linearSolveLS zing (LA.col vs)
        v = leftPseudo LA.#> LA.vector vs
        f i = (m LA.! i) =~ (v LA.! i)

-- | Approximate Equality.
(=~) :: Double -> Double -> Bool
a =~ b = abs (a - b) < 0.0001

zing :: LA.Matrix Double
zing = LA.matrix 3 [ -0.5, -0.5, 1
                   , -0.5, 0, 1
                   , -0.5, 0.5, 1
                   , 0, -0.5, 1
                   , 0, 0, 1
                   , 0, 0.5, 1
                   , 0.5, -0.5, 1
                   , 0.5, 0, 1
                   , 0.5, 0.5, 1 ]

fgradientFlat :: Assertion
fgradientFlat = actual @?= expected
  where expected :: Raster U p 3 3 Double
        expected = fromRight . fromVector Seq $ U.fromList [0,0,0,0,0,0,0,0,0]
        actual :: Raster U p 3 3 Double
        actual = strict U . fgradient . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,1,1,1,1,1] :: [Double])

fgradient45 :: Assertion
fgradient45 = let ?epsilon = 0.0001 in index' (_array actual) (1 :. 1) @?~ (tau / 8)
  where actual :: Raster U p 3 3 Double
        actual = strict U . fgradient . fromRight . fromVector Seq $ U.fromList ([3,3,3,2,2,2,1,1,1] :: [Double])

faspectFlat :: Assertion
faspectFlat = index' (_array actual) (1 :. 1) @?= Nothing
  where actual :: Raster B p 3 3 (Maybe Double)
        actual = strict B . faspect . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,1,1,1,1,1] :: [Double])

faspect45 :: Assertion
faspect45 = index' (_array actual) (1 :. 1) @?= Just (tau / 2)
  where actual :: Raster B p 3 3 (Maybe Double)
        actual = strict B . faspect . fromRight . fromVector Seq $ U.fromList ([3,3,3,2,2,2,1,1,1] :: [Double])

faspectEast :: Assertion
faspectEast = let ?epsilon = 0.0001 in index' (_array actual) (1 :. 1) @?~ (tau / 4)
  where actual :: Raster B p 3 3 Double
        actual = strict B . faspect' . fromRight . fromVector Seq $ U.fromList ([3,2,1,3,2,1,3,2,1] :: [Double])

fdownstream4 :: Assertion
fdownstream4 = index' (_array actual) (1 :. 1) @?= drainage (S.fromList [North,South,East,West])
  where actual :: Raster S p 3 3 Drain
        actual = strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([3,1,3,1,2,1,3,1,3] :: [Double])

fdownstreamFlat :: Assertion
fdownstreamFlat = index' (_array actual) (1 :. 1) @?= drainage (S.fromList [East ..])
  where actual :: Raster S p 3 3 Drain
        actual = strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,1,1,1,1,1] :: [Double])

fdownstreamPeak :: Assertion
fdownstreamPeak = index' (_array actual) (1 :. 1) @?= drainage (S.fromList [NorthEast, NorthWest, SouthWest, SouthEast])
  where actual :: Raster S p 3 3 Drain
        actual = strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,3,1,1,1,1] :: [Double])

fdownstreamPit :: Assertion
fdownstreamPit = index' (_array actual) (1 :. 1) @?= Drain 0
  where actual :: Raster S p 3 3 Drain
        actual = strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([2,2,2,2,1,2,2,2,2] :: [Double])

fupstreamFlat :: Assertion
fupstreamFlat = index' (_array actual) (1 :. 1) @?= drainage (S.fromList [East ..])
  where actual :: Raster S p 3 3 Drain
        actual = strict S . fupstream . strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,1,1,1,1,1] :: [Double])

fupstreamPeak :: Assertion
fupstreamPeak = index' (_array actual) (1 :. 1) @?= Drain 0
  where actual :: Raster S p 3 3 Drain
        actual = strict S . fupstream . strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,3,1,1,1,1] :: [Double])

{-
hists :: Raster S p 512 512 Word8 -> Assertion
hists r = do
  u  <- histU r
  um <- histMut r
  histU' r @?= u
  histMut' r @?= um
  histMut'' r @?= um
  histMut'' r @?= u
  um @?= u
-}

max_stars_repo_name: fosskers/mapalgebra

max_issues_repo_name: fosskers/mapalgebra

max_forks_repo_name: fosskers/mapalgebra


text:
{-# LANGUAGE Strict, OverloadedStrings #-}
module GradeParser (
    outOf
  , outOfSafe
  , eachOutOf
  , eachOutOfSafe
  , stats
  , histOutput
  , printStats
  , printStatsHist
  ) where

import qualified Data.ByteString as B
import Grade
import qualified Data.Vector as T
import Data.Csv hiding ((.=))
import ShowByteString
import qualified Data.HashMap.Strict as M
import Control.Monad.Reader
import Control.Applicative
import Histogram
import Data.List
import Statistics.Sample
import CSVmonad
import System.Process

outOf :: Field -> Double -> CSVmonad Grade
outOf s o = ReaderT (\v -> ((/. o) <$> v .: s) <|> pure (absent o))

eachOutOf s o = ReaderT (\v -> sequence (map (parseGrade o) $ M.elems $ M.filterWithKey (\k _ -> B.isPrefixOf s k) v))

outOfSafe :: Field -> Double -> CSVmonad Grade
outOfSafe s o = ReaderT (\v -> (chkGrade o <$> v .: s) <|> pure (absent o))

eachOutOfSafe s o = ReaderT (\v -> sequence (map (parseGradeSafe o) $ M.elems $ M.filterWithKey (\k _ -> B.isPrefixOf s k) v))

chkGrade :: Double -> Double -> Grade
chkGrade o n = if n <= o && n >= 0 then n /. o else  outoferror

chkGrade' :: Double -> Double -> Grade
chkGrade' _ _ = outoferror

parseGrade o g = ((/.o) <$> parseField g) <|> pure (absent o)

parseGradeSafe o g = (chkGrade o <$> parseField g) <|> pure (absent o)

getRight (Right x) = x
getRight (Left e) = error $ show e

stats :: (T.Vector Double -> b) -> CSVmonad Grade -> T.Vector NamedRecord -> Either GradeErrors b
stats f l v = fmap f $ sequence $ T.filter isNotAbsent $ fmap (marksEither . getRight . runParser .  runReaderT l ) v
  where isNotAbsent (Left (Absent _)) = False
        isNotAbsent _ = True

histOutput :: Double -> FilePath -> CSVmonad Grade -> T.Vector NamedRecord -> IO ()
histOutput m s l v = case stats (histogram (assignBin m) . T.toList) l v of
                        Right a -> histToChart s a >> writeHist "distribution.txt" a >> callCommand "perl -lane 'print $F[0], \"\t\", \"=\" x ($F[1] / 1)' distribution.txt" >> return ()
                        --Right a -> histToChart s a >> writeHist "distribution.txt" a >> createProcess (shell "perl -lane 'print $F[0], \"\t\", \"=\" x ($F[1] / 1)' distribution.txt") >> return ()
                        Left e -> print e

eitherPrintDouble (Right x) = print $ roundTo 1 x
eitherPrintDouble (Left x) = print x

printStat n s l y = putStrLn n >> eitherPrintDouble (stats s l y) >> putStrLn ""

--printStats l y = do printStat "The Mean" mean l y
--                    printStat "The Standard deviation" stdDev l y
--                    printStat "The Median" median l y
--                    histOutput 10 "histogram.png" l y
--
--printStats l y = do let m = getRight $ stats mean l y
--                    let s = getRight $ stats stdDev l y
--                    let me = getRight $ stats median l y
--                    printStat "The Mean" mean l y
--                    printStat "The Standard deviation" stdDev l y
--                    printStat "The Median" median l y
--                    histOutput 5 "histogram.png" l y
--                    let sugg= ("      O: " ++ show (roundTo 1 $ m + 1.65*s) ++"\n      "++ "A: " ++ show (roundTo 1 $ m + 0.85*s) ++"\n      "++ "B: " ++ show (roundTo 1 $ m) ++"\n      "++ "C: " ++ show (roundTo 1 $ m - s) ++"\n      "++ "D: " ++ show (roundTo 1 $ me / 2))
--                    putStrLn sugg
--                    writeFile "suggested" sugg

printStats = printStatsHist 10
printStatsHist  h l y = do let m = getRight $ stats mean l y
                           let s = getRight $ stats stdDev l y
                           let me = getRight $ stats median l y
                           printStat "The Mean" mean l y
                           printStat "The Standard deviation" stdDev l y
                           printStat "The Median" median l y
                           histOutput h "histogram.png" l y
                           let sugg= ("      O: " ++ show (roundTo 1 $ m + 1.65*s) ++"\n      "++ "A: " ++ show (roundTo 1 $ m + 0.85*s) ++"\n      "++ "B: " ++ show (roundTo 1 $ m) ++"\n      "++ "C: " ++ show (roundTo 1 $ m - s) ++"\n      "++ "D: " ++ show (roundTo 1 $ me / 2))
                           putStrLn "  "
                           putStrLn "Suggested grades:"
                           putStrLn sugg
                           writeFile "suggested" sugg

round' x | f < 0.5 = i
         | otherwise = i + 1
         where i = floor x
               f = x - fromIntegral i

roundTo n f = fromInteger (round' $ f * (10^n)) / (10.0^^n)

median :: T.Vector Double -> Double
median xs = sort (T.toList xs) !! n
  where n = T.length xs `div` 2

max_stars_repo_name: sejdm/Grades

max_issues_repo_name: sejdm/Grades

max_forks_repo_name: sejdm/Grades