text:
module STCR2S1 where

import           Control.Monad             as M
import           Control.Monad.Parallel    as MP
import           Data.Array.Repa           as R
import           Data.Binary               (decodeFile)
import           Data.Complex
import           Data.List                 as L
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.MonteCarlo
import           Image.IO
import           STC.CompletionFieldR2S1
import           System.Directory
import           System.Environment
import           System.FilePath
import           Types
import           Utils.Array
import           Utils.Time
import STC.Utils
import Text.Printf
import           FokkerPlanck.FourierSeries
import FokkerPlanck.GreensFunction

main = do
  args@(numPointStr:numOrientationStr:sigmaStr:taoStr:thresholdStr:numTrailStr:maxTrailStr:initDistStr:initOriStr:initSpeedStr:rStr:locationStr:deltaStr:cornerWeightStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      sigma = read sigmaStr :: Double
      tao = read taoStr :: Double
      threshold = read thresholdStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      initDist = read initDistStr :: [R2S1RPPoint]
      initOri = read initOriStr :: Double
      initSpeed = read initSpeedStr :: Double
      r = read rStr :: Double
      location = read locationStr :: (Int, Int)
      delta = read deltaStr :: Double
      cornerWeight = read cornerWeightStr :: Double
      numThread = read numThreadStr :: Int
      sourceDist =
        computeInitialDistribution numPoint numPoint numOrientation delta . L.take 1 $
        initDist
      sinkDist =
        computeInitialDistribution numPoint numPoint numOrientation delta . L.drop 1 $
        initDist
      dist =
        computeInitialDistribution numPoint numPoint numOrientation delta initDist
      folderPath = "output/test/STCR2S1"
  createDirectoryIfMissing True folderPath
  -- Compute the Green's function
  -- arrG <-
  --   solveMonteCarloR2S1
  --     numThread
  --     numTrail
  --     maxTrail
  --     numPoint
  --     numPoint
  --     numOrientation
  --     sigma
  --     tao
  --     r
  --     initSpeed
  --     ""
  let deltaTheta = 2 * pi / fromIntegral numOrientation
  -- arrG <-
  --   sampleR2S1
  --     numPoint
  --     numPoint
  --     delta
  --     deltaTheta
  --     initSpeed
  --     sigma
  --     tao
  --     [0 .. numOrientation - 1]
  arrG <-
    sampleR2S1Corner
      numPoint
      numPoint
      delta
      deltaTheta
      initSpeed
      sigma
      tao
      threshold
      [0 .. numOrientation - 1]
      cornerWeight
  printCurrentTime "Done."
  plan <- makeR2S1Plan emptyPlan arrG
  -- Source Field
  source <- shareWeightST plan sourceDist arrG
  plotImageRepa
    (folderPath </> "Source.png")
    (ImageRepa 8 .
     computeS .
     R.extend (Z :. (1 :: Int) :. All :. All) .
     reduceContrast 100 . R.sumS . rotate3D $
     source)
  plotThetaDimension folderPath "R2S1_" location .
    R.backpermute
      (extent source)
      (\(Z :. k :. i :. j) ->
         (Z :.
          (mod (k + numOrientation - (div numOrientation 2)) numOrientation :: Int) :.
          i :.
          j)) $
    source
  -- MP.mapM_
  --   (\i ->
  --      plotImageRepa
  --        (folderPath </> printf "Source_%d.png" i)
  --        (ImageRepa 8 .
  --         computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.slice source $
  --         (Z :. i :. All :. All)))
  --   [0 .. numOrientation - 1]
  -- Sink Field
  sink <- shareWeightST plan sinkDist arrG -- . rotateST arrG $ div numOrientation 2
  plotImageRepa
    (folderPath </> "Sink.png")
    (ImageRepa 8 .
     computeS . R.extend (Z :. (1 :: Int) :. All :. All) . reduceContrast 100 . R.sumS . rotate3D $
     sink)
  -- MP.mapM_
  --   (\i ->
  --      plotImageRepa
  --        (folderPath </> printf "Sink_%d.png" i)
  --        (ImageRepa 8 .
  --         computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.slice sink $
  --         (Z :. i :. All :. All)))
  --   [0 .. numOrientation - 1]
  -- Completion Field
  let completion = R.zipWith (*) source . timeReversal $ sink
  plotImageRepa
    (folderPath </> "Completion.png")
    (ImageRepa 8 .
     computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.sumS . rotate3D $
     completion)

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 TypeOperators #-}
{-# LANGUAGE IncoherentInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MonoLocalBinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE GADTs #-}
module FFT where

import Prelude ( Double, ($), Integer, Num(..), Int, (^), Bool(..), otherwise )
-- import qualified Prelude

import Data.C ( CVal )
import Data.Complex
import Control.CCat
import Control.CArr
import Control.CArr.CSyn ( SINat, Sing(..), Tree
                         , cdictTree, withCDict, toInteger, withSize )
import Language.SPar.Skel ( (:->), liftAlg )

tsplit' :: forall a n f. (CVal a, CArrLift (:->) f, PAlg f) => SINat n -> Int
        -> (Int -> Int -> a :-> (a, a))
        -> f a (Tree n a)
tsplit' SZ _acc _f = newProc
tsplit' (SS n) acc f =
  withCDict (cdictTree @a n) $ lift (f wl wr) >>>
  (tsplit' n wl f *** tsplit' n wr f)
  where
    wl = acc
    wr = acc + 2 ^ toInteger n

tsplit :: forall a n f. (CArrLift (:->) f, PAlg f, CVal a)
       => SINat n
       -> (Int -> Int -> a :-> (a, a))
       -> f a (Tree n a)
tsplit n = tsplit' n 0

tfold  :: forall f a n. (PAlg f, CVal a)
       => SINat n -> f (a, a) a -> f (Tree n a) a
tfold  SZ _f = id
tfold (SS n) f =
  withCDict (cdictTree @a n) $ (tfold n f *** tfold n f) >>> f

zipTree :: forall a b n f. (CVal a, CVal b, PAlg f, CArrLift (:->) f)
        => SINat n
        -> Prelude.Bool
        -> Int
        -> Int
        -> ((Int, Int), (a, a)) :-> b
        -> f (Tree n a, Tree n a) (Tree n b)
zipTree SZ b l w f
  | b = lift ((lit l &&& lit w) &&& id >>> f)
  | otherwise = (snd &&& fst) >>> lift ((lit l &&& lit w) &&& id >>> f)
zipTree (SS x) b l w f =
  withCDict (cdictTree @a x) $
  withCDict (cdictTree @b x) $
  let swap = ((fst >>> fst) &&& (snd >>> fst)) &&&
             ((fst >>> snd) &&& (snd >>> snd))
  in swap >>> (zipTree x b l w f *** zipTree x b l (w + 2 ^ toInteger x) f)

type D a = a

fmapT :: PAlg f
        => SINat n
        -> f (D [Complex Double]) (D [Complex Double])
        -> f (Tree n (D [Complex Double])) (Tree n (D [Complex Double]))
fmapT SZ f = f
fmapT (SS x) f
  = withCDict (cdictTree @(D [Complex Double]) x) $ fmapT x f *** fmapT x f

fmapTIx :: PAlg f
        => SINat n
        -> f (Int, D [Complex Double]) (D [Complex Double])
        -> Int
        -> f (Tree n (D [Complex Double])) (Tree n (D [Complex Double]))
fmapTIx SZ f k = lit k &&& id >>> f
fmapTIx (SS x) f k
  = withCDict (cdictTree @(D [Complex Double]) x) $
  fmapTIx x f k *** fmapTIx x f (k + (2 ^ (toInteger x :: Integer)))

-- addPadding :: PAlg f => f (D [Complex Double]) (D [Complex Double])
-- addPadding = prim "add_padding"

deinterleave :: Int -> Int -> (D [Complex Double]) :-> (D [Complex Double], D [Complex Double])
deinterleave wl wr = lit wl &&& lit wr >>> prim "deinterleave"

intlit :: (PAlg f, CVal a) => SINat n -> f a Int
intlit i = fromInteger $ toInteger i + 1

fftTree :: (CArrLift (:->) f, PAlg f)
        => SINat n
        -> Int
        -> f (Tree n (D [Complex Double])) (Tree n (D [Complex Double]))
fftTree SZ w
  = liftAlg (intlit SZ &&& (lit w &&& id) >>> prim "baseFFT")
fftTree (SS x) w
  = withCDict (cdictTree @(D [Complex Double]) x) $
    (fftTree x w {-EVENS-} *** fftTree x (w + 2^ toInteger x){-ODDS-}) >>>
    id *** fmapTIx x (lit ps2x &&& id >>> prim "map_exp") 0 {- Multiply left side by exponential -} >>>
    zipTree x True lvl w addc &&& zipTree x False lvl (w + 2^ toInteger x) subc
  where
    lvl :: Int
    lvl = fromInteger (toInteger (SS x) + 1)
    ps2x :: Int
    ps2x = 2 ^ toInteger (SS x)
    addc :: ((Int, Int), (D [Complex Double], D [Complex Double])) :-> (D [Complex Double])
    addc = prim "zip_add"
    subc :: ((Int, Int), (D [Complex Double], D [Complex Double])) :-> (D [Complex Double])
    subc = id *** (snd &&& fst) >>> prim "zip_sub"

fft :: (CArrLift (:->) f, PAlg f)
    => SINat n -> f (D [Complex Double]) (D [Complex Double])
fft n =
  withCDict (cdictTree @(D [Complex Double]) n) $
  -- addPadding >>>
  tsplit n deinterleave >>> fftTree n 0 >>> tfold n (runAt 0 >>> prim "cat")

fft0 :: (CArrLift (:->) f, PAlg f)
     => f (D [Complex Double]) (D [Complex Double])
fft0 = withSize 0 fft

fft1 :: (CArrLift (:->) f, PAlg f)
     => f (D [Complex Double]) (D [Complex Double])
fft1 = withSize 1 fft

fft2 :: (CArrLift (:->) f, PAlg f)
     => f (D [Complex Double]) (D [Complex Double])
fft2 = withSize 2 fft

fft3 :: (CArrLift (:->) f, PAlg f)
     => f (D [Complex Double]) (D [Complex Double])
fft3 = withSize 3 fft

fft4 :: (CArrLift (:->) f, PAlg f)
     => f (D [Complex Double]) (D [Complex Double])
fft4 = withSize 4 fft

fft5 :: (CArrLift (:->) f, PAlg f)
     => f (D [Complex Double]) (D [Complex Double])
fft5 = withSize 5 fft

fft6 :: (CArrLift (:->) f, PAlg f)
     => f (D [Complex Double]) (D [Complex Double])
fft6 = withSize 6 fft

fft7 :: (CArrLift (:->) f, PAlg f)
     => f (D [Complex Double]) (D [Complex Double])
fft7 = withSize 7 fft

fft8 :: (CArrLift (:->) f, PAlg f)
     => f (D [Complex Double]) (D [Complex Double])
fft8 = withSize 8 fft
--
-- fft8 :: PAlg f => f (D [Complex Double]) (D [Complex Double])
-- fft8 = withSize 8 fft
  --where
  --  p2sx :: Integer
  --  p2sx = 2 ^ (fromINat (SS x) :: Integer)

--fft :: PAlg f => SINat n
--    -> f (Tree n [Complex Double]) (Tree n [Complex Double])

-- Below fails for some reason!
--fft :: PAlg f => SINat n
--    -> f (Tree n [Complex Double]) (Tree n [Complex Double])
--fft SZ = cfun $ prim "baseFFT"
--fft (SS n) = cfun $ \x ->
--  vlet (withCDict (cdictTree @[Complex Double] n) $ app (fft n) $ fst x) $ \l ->
--  vlet (app (fft n) $ snd x) $ \r -> _
  -- vlet (par (fft n) $ fst x) $ \l -> _
--fft n = cfun $ \z ->
--  case n of
--    SZ -> prim "baseFFT" z
--    (SS x) ->
--      vlet (par (app $ fft x) $ fst z) $ \l ->
--      vlet (par (app $ fft x) $ snd z) $ \r ->
--      _

max_stars_repo_name: session-arr/session-arr

max_issues_repo_name: session-arr/session-arr

max_forks_repo_name: session-arr/session-arr


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UndecidableInstances #-}

-- |
-- Module      :  Spiral.Array.Repr.Complex
-- Copyright   :  (c) 2016 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Spiral.Array.Repr.Complex (
    RE,
    toRealArray,

    CMPLX,
    toComplexArray,

    Array(..)
  ) where

import Prelude hiding (read)

import Data.Complex (Complex)
import Data.Monoid ((<>))
import Text.PrettyPrint.Mainland
import Text.PrettyPrint.Mainland.Class

import Spiral.Array
import Spiral.Exp

data RE r

toRealArray :: Array r (sh :. Int) (Exp (Complex a)) -> Array (RE r) (sh :. Int) (Exp a)
toRealArray = RE

instance IsArray r (sh :. Int) (Exp (Complex a)) => IsArray (RE r) (sh :. Int) (Exp a) where
    data Array (RE r) (sh :. Int) (Exp a) = RE (Array r (sh :. Int) (Exp (Complex a)))

    extent (RE a) = sh :. 2*n
      where
        sh :. n = extent a

instance Pretty (Array r (sh :. Int) (Exp (Complex a))) => Pretty (Array (RE r) (sh :. Int) (Exp a)) where
    ppr (RE a) = text "Re" <> parens (ppr a)

instance (RealFloatConst a, Shape sh, IArray r (sh :. Int) (Exp (Complex a))) => IArray (RE r) (sh :. Int) (Exp a) where
    index (RE a) (sh :. i)
        | i `rem` 2 == 0 = re
        | otherwise      = im
      where
        (re, im) = unComplexE (index a (sh :. i `quot` 2))

instance (RealFloatConst a, SArray r (sh :. Int) (Exp (Complex a))) => SArray (RE r) (sh :. Int) (Exp a) where
    indexS (RE a) (sh :. i)
        | i `rem` 2 == 0 = re
        | otherwise      = im
      where
        (re, im) = unComplexE (indexS a (sh :. i `quot` 2))

instance (RealFloatConst a, MArray r (sh :. Int) (Exp (Complex a))) => MArray (RE r) (sh :. Int) (Exp a) where
    read (RE a) (sh :. i) = do
        (re, im) <- unComplexE <$> read a (sh :. i `quot` 2)
        if i `rem` 2 == 0
          then return re
          else return im

    write (RE a) (sh :. i) e = do
        (re, im) <- unComplexE <$> read a (sh :. i `quot` 2)
        if i `rem` 2 == 0
          then write a (sh :. i `quot` 2) (ComplexE e im)
          else write a (sh :. i `quot` 2) (ComplexE re e)

instance (RealFloatConst a, SArray r (sh :. Int) (Exp (Complex a))) => Computable (RE r) (sh :. Int) (Exp a) where
    computeP a b =
        forShapeP (extent b) $ \ix ->
            write a ix (indexS b ix)

data CMPLX r

toComplexArray :: Array r (sh :. Int) (Exp a) -> Array (CMPLX r) (sh :. Int) (Exp (Complex a))
toComplexArray = CMPLX

instance IsArray r (sh :. Int) (Exp a) => IsArray (CMPLX r) (sh :. Int) (Exp (Complex a)) where
    data Array (CMPLX r) (sh :. Int) (Exp (Complex a)) = CMPLX (Array r (sh :. Int) (Exp a))

    extent (CMPLX a) = sh :. n `quot` 2
      where
        sh :. n = extent a

instance Pretty (Array r (sh :. Int) (Exp a)) => Pretty (Array (CMPLX r) (sh :. Int) (Exp (Complex a))) where
    ppr (CMPLX a) = text "Cmplx" <> parens (ppr a)

instance (Typed a, Num (Exp a), IArray r (sh :. Int) (Exp a)) => IArray (CMPLX r) (sh :. Int) (Exp (Complex a)) where
    index (CMPLX a) (sh :. i) = ComplexE (index a (sh :. 2*i)) (index a (sh :. 2*i+1))

instance (Typed a, Num (Exp a), SArray r (sh :. Int) (Exp a)) => SArray (CMPLX r) (sh :. Int) (Exp (Complex a)) where
    indexS (CMPLX a) (sh :. i) = ComplexE (indexS a (sh :. 2*i)) (indexS a (sh :. 2*i+1))

instance (RealFloatConst a, MArray r (sh :. Int) (Exp a)) => MArray (CMPLX r) (sh :. Int) (Exp (Complex a)) where
    read (CMPLX a) (sh :. i) = ComplexE <$> read a (sh :. 2*i) <*> read a (sh :. 2*i+1)

    write (CMPLX a) (sh :. i) e = do
        write a (sh :. 2*i) er
        write a (sh :. 2*i+1) ei
      where
        (er, ei) = unComplexE e

instance (RealFloatConst a, Computable r (sh :. Int) (Exp a)) => Computable (CMPLX r) (sh :. Int) (Exp (Complex a)) where
    computeP a (CMPLX b) = computeP (toRealArray a) b

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
module Numeric.Digamma (digamma) where

import Numeric.SpecFunctions (digamma)

max_stars_repo_name: bgamari/digamma

max_issues_repo_name: bgamari/digamma

max_forks_repo_name: bgamari/digamma


text:
{-# LANGUAGE FlexibleInstances #-}

import Debug.Trace
-- import Control.Monad.State.Lazy

import System.Environment (getArgs)
import System.IO (readFile)
import Data.Map.Strict (Map, (!), insert, elems, fromList, toList, findWithDefault, size, empty, member, findMin, findMax, singleton, filter)
import qualified Data.Map.Strict as M
import Data.Set (Set)
import qualified Data.Set as S
--import qualified Data.Array as A
import Data.List (find, intercalate, intersperse, permutations, inits, tails, isPrefixOf)
import Data.List.Split (splitOn)
-- import Data.Complex (Complex((:+)), realPart, imagPart) -- define my own complex

-- import UI.NCurses
-- import Data.Char (chr, ord)
-- -- import Data.Complex (Complex((:+)))
import Data.Maybe (fromJust)
type Instructions = Map Integer Integer

data ComputerState = Ready | Running | Blocked | Done
  deriving (Show, Eq)

data Computer = Computer {
  state :: ComputerState,
  memory :: Map Integer Integer,
  iptr :: Integer, -- instruction pointer
  base :: Integer, -- base offset
  input :: [Integer],
  output :: [Integer],
  time :: Integer}

computer0 = Computer { state = Ready, memory = empty, iptr = 0, base = 0, input = [], output = [], time = 0 }

run1 :: Computer -> Computer -- (Integer, Integer) -> Instructions -> [Integer] -> [Integer]
run1 c
  | state c == Ready = run1 $ c {state = Running}
  | state c == Done = c
  | state c == Blocked = if null (input c) then c else run1 $ c {state = Running}
  | otherwise = -- running
    --traceShow (iptr c, base c, memory c) $
    case instr `mod` 100 of
      1 -> -- add
        addTime $ c {iptr = i+4, memory = insert (addr 3) (arg 1 + arg 2) instructions}
      2 -> -- multiply
        addTime $ c {iptr = i+4, memory = insert (addr 3) (arg 1 * arg 2) instructions}
      3 -> -- read input
        if (null $ input c)
        then c {state = Blocked}
        else addTime $ c {iptr = i+2, memory = insert (addr 1) (head $ input c) instructions, input = tail $ input c}
      4 -> -- output
        --traceShow (arg 1) $
        addTime $ c {iptr = i+2, output = output c ++ [arg 1]}
      5 -> -- jump-if-true
        addTime $ c {iptr = if arg 1 == 0 then i+3 else arg 2}
      6 -> -- jump-if-false
        addTime $ c {iptr = if arg 1 == 0 then arg 2 else i+3}
      7 -> -- less than
        addTime $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 < arg 2 then 1 else 0) instructions}
      8 -> -- equals
        addTime $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 == arg 2 then 1 else 0) instructions}
      9 -> -- set relative base
        addTime $ c {iptr = i+2, base = base c + arg 1}
      99 -> -- halt
        addTime $ c {state = Done}
      _ -> error "unknown opcode"
  where instructions = memory c
        i = iptr c
        instr = instructions!i
        ii x = findWithDefault 0 x instructions
        arg :: Integer -> Integer
        arg n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ ii $ i+n
          1 -> ii $ i+n
          2 -> ii ((ii $ i+n) + base c)
          _ -> error "bad argument mode"
        addr n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ i+n
          1 -> error "address in mode 1"
          2 -> --trace "address in mode 2" $
            ii (i+n) + base c
          _ -> error $ "address in unknown mode"
        addTime c = c{time = time c + 1}

run :: Computer -> Computer
run c | state c == Done = c
      | state c == Blocked = if null (input c) then c else run $ run1 c
      | length (output c) == 3 = c
      | otherwise = run $ run1 c

runNetwork :: Map Int Computer -> Set (Integer, Int) -> Integer
runNetwork computers heap
  | isDone = y
  | otherwise = runNetwork computers'' heap''
  where
    ((_, i), heap') = S.deleteFindMin heap -- find the computer with lowest time
    c = computers!i
    hasOutput = 3 <= length (output c)
    [addr_, x, y] = take 3 $ output c -- laziness means this won't be executed unless needed
    addr = fromIntegral addr_
    isDone = hasOutput && addr == 255
    c' = if hasOutput then c{output = drop 3 $ output c} else c
    c'' = if state c' == Blocked && null (input c') then c'{input = [-1]} else c'
    c''' = run c''
    c2 = computers!addr
    computers' = if hasOutput then M.insert addr (c2{input = input c2 ++ [x,y]}) computers else computers
    computers'' = M.insert i c''' computers'
    heap'' = S.insert (time c''', i) heap'

main = do
  -- [instructionFile] <- getArgs
  instructionStrings <- readFile "23.input.txt" -- instructionFile
  let instructions = fromList . zip [0 ..] $ map read $ splitOn "," instructionStrings

  putStrLn "Part 1"
  let addrs = [0 .. 49]
      computers = M.fromList [(i, computer0{memory=instructions, input=[fromIntegral i]}) | i <- addrs]
      heap = S.fromList [(0,i) | i <- addrs]
  print $ runNetwork computers heap


text:
-- | The method of moments can be used to estimate a number of commonly used distributions.  This module is still under construction as I work out the best way to handle morphisms from the Moments3 type to types of other distributions.  For more information, see the wikipedia entry: <https://en.wikipedia.org/wiki/Method_of_moments_(statistics)>

module HLearn.Models.Distributions.Univariate.Poisson
    ( Poisson
    )
    where

import Control.DeepSeq
import GHC.TypeLits
import qualified Data.Vector.Unboxed as U
import Data.Vector.Unboxed.Deriving

import qualified Statistics.Distribution as S
import Statistics.Distribution.Binomial
import Statistics.Distribution.Poisson

import HLearn.Algebra
import HLearn.Models.Distributions.Common
import HLearn.Models.Distributions.Univariate.Internal.Moments
import HLearn.Models.Distributions.Visualization.Gnuplot

-------------------------------------------------------------------------------
-- Poisson

newtype Poisson sample prob = Poisson {  pmoments :: (Moments3 sample) }
    deriving (Read,Show,Eq,Ord,Monoid,Group)

instance (Num sample) => HomTrainer (Poisson sample prob) where
    type Datapoint (Poisson sample prob) = sample
    train1dp dp = Poisson $ train1dp dp

instance (Num sample) => Probabilistic (Poisson sample prob) where
    type Probability (Poisson sample prob) = prob

-- instance (Integral sample, Floating prob) => PDF (Poisson sample prob) where
--     pdf (Poisson dist) dp
--         | dp < 0    = 0
--         | dp > 100  = pdf (Normal $ Moments3 (fromIntegral $ m0 dist) (fromIntegral $ m1 dist) (fromIntegral $ m2 dist)) $ fromIntegral dp
--         | otherwise = (lambda^^dp)*(exp $ (-1)*lambda)/(fromIntegral $ factorial dp)
--         where
--             lambda = (fromIntegral $ m1 dist) / (fromIntegral $ m0 dist)
-- factorial :: (Integral a) => a -> a
-- factorial 0 = 1
-- factorial n = n*(factorial $ n-1)

instance (Integral sample, Floating prob) => PDF (Poisson sample Double) where
    pdf (Poisson dist) dp = S.probability (poisson lambda) $ fromIntegral dp
        where
            lambda = (fromIntegral $ m1 dist) / (fromIntegral $ m0 dist)

instance
    ( PDF (Poisson sample prob)
--     , PlottableDataPoint sample
    , Show prob
    , Show sample
    , Ord sample
    , Ord prob
    , Fractional prob
    , Integral sample
    ) => PlottableDistribution (Poisson sample prob)
-- instance PlottableDistribution (Poisson Int Double)
        where

    plotType _ = Points

    samplePoints dist = [min..max]
        where
            min = 0
            max = maximum [20,floor $ 3*lambda]
            lambda = (fromIntegral $ m1 $ pmoments dist) / (fromIntegral $ m0 $ pmoments dist)

-- instance (Fractional prob) => Mean (Poisson sample prob) where
--     mean (Poisson dist) = (fromIntegral $ m1 $ pmoments dist) / (fromIntegral $ m0 $ pmoments dist)
--
-- instance (Fractional prob) => Variance (Poisson sample prob) where
--     variance dist = mean dist


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}

import           Control.DeepSeq
import           Control.Exception
import           Control.Monad
import           Control.Monad.IO.Class
import           Control.Monad.Trans.Maybe
import           Control.Monad.Trans.State
import           Data.Bifunctor
import           Data.Bitraversable
import           Data.Default
import           Data.Finite
import           Data.IDX
import           Data.List.Split
import           Data.Time.Clock
import           Data.Traversable
import           Data.Tuple
import           Data.Type.Product
import           Learn.Neural
import           Numeric.BLAS.HMatrix
import           Numeric.LinearAlgebra.Static
import           Text.Printf
import qualified Data.Vector                     as V
import qualified Data.Vector.Generic             as VG
import qualified Data.Vector.Unboxed             as VU
import qualified System.Random.MWC               as MWC
import qualified System.Random.MWC.Distributions as MWC

loadMNIST
    :: FilePath
    -> FilePath
    -> IO (Maybe [(HM '[784], HM '[10])])
loadMNIST fpI fpL = runMaybeT $ do
    i <- MaybeT          $ decodeIDXFile       fpI
    l <- MaybeT          $ decodeIDXLabelsFile fpL
    d <- MaybeT . return $ labeledIntData l i
    r <- MaybeT . return $ for d (bitraverse mkImage mkLabel . swap)
    liftIO . evaluate $ force r
  where
    mkImage :: VU.Vector Int -> Maybe (HM '[784])
    mkImage = fmap HM . create . VG.convert . VG.map (\i -> fromIntegral i / 255)
    mkLabel :: Int -> Maybe (HM '[10])
    mkLabel = fmap (oneHot . only) . packFinite . fromIntegral

main :: IO ()
main = MWC.withSystemRandom $ \g -> do
    Just train <- loadMNIST "data/train-images-idx3-ubyte" "data/train-labels-idx1-ubyte"
    Just test  <- loadMNIST "data/t10k-images-idx3-ubyte"  "data/t10k-labels-idx1-ubyte"
    putStrLn "Loaded data."
    net0 :: Network 'FeedForward HM ( '[784] :~ FullyConnected )
                                   '[ '[300] :~ LogitMap
                                    , '[300] :~ FullyConnected
                                    , '[100] :~ LogitMap
                                    , '[100] :~ FullyConnected
                                    , '[10 ] :~ SoftMax '[10]
                                    ]
                                    '[10] <- initDefNet g
    let dout = alongNet net0 $ Nothing
                           :&% Just 0.2
                           :&% Nothing
                           :&% Just 0.2
                           :&% Nothing
                           :&% DOExt
    flip evalStateT net0 . forM_ [1..] $ \e -> do
      train' <- liftIO . fmap V.toList $ MWC.uniformShuffle (V.fromList train) g
      liftIO $ printf "[Epoch %d]\n" (e :: Int)

      forM_ ([1..] `zip` chunksOf batch train') $ \(b, chnk) -> StateT $ \n0 -> do
        printf "(Batch %d)\n" (b :: Int)

        t0 <- getCurrentTime
        -- n' <- evaluate $ optimizeList_ (bimap only_ only_ <$> chnk) n0
        --                                -- (sgdOptimizer rate netOpPure crossEntropy)
        --                                (adamOptimizer def netOpPure crossEntropy)
        n' <- optimizeListM_ (bimap only_ only_ <$> chnk) n0
                             (adamOptimizerM def (netOpDOPure dout g) crossEntropy)
        t1 <- getCurrentTime
        printf "Trained on %d points in %s.\n" batch (show (t1 `diffUTCTime` t0))

        let trainScore = testNetList maxTest (someNet n') chnk
            testScore  = testNetList maxTest (someNet n') test
        printf "Training error:   %.2f%%\n" ((1 - trainScore) * 100)
        printf "Validation error: %.2f%%\n" ((1 - testScore ) * 100)

        return ((), n')
  where
    -- rate :: Double
    -- rate  = 0.02
    batch :: Int
    batch = 2500

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
module EvalSpec where

import Data.Array
import Test.Hspec
import Interpreter
import Expression
import Operation
import Prelude hiding ((+), (*))
import qualified Prelude as P
import Data.Complex

spec :: Spec
spec = do
    describe "eval basics" $ do
        specify "real scalar 01" $ do
            let x = var "x"
                y = var "y"
                valMap = subs [("x", 1), ("y", 3)] []
            eval valMap x `shouldBe` 1
            eval valMap y `shouldBe` 3
            eval valMap (x + y) `shouldBe` 4
            eval valMap (x * y) `shouldBe` 3
            eval valMap (x .* y) `shouldBe` 3
            eval valMap (x <.> y) `shouldBe` 3
        specify "real scalar 02" $ do
            let x = var "x"
                y = var "y"
                valMap = subs [("x", 100), ("y", 300)] []
            eval valMap x `shouldBe` 100
            eval valMap y `shouldBe` 300
            eval valMap (x + y) `shouldBe` 400
            eval valMap (x * y) `shouldBe` 30000
            eval valMap (x .* y) `shouldBe` 30000
            eval valMap (x <.> y) `shouldBe` 30000
        specify "complex scalar 01" $ do
            let x = var "x"
                y = var "y"
                z1 = x +: y
                z2 = y +: x
                valMap = subs [("x", 1), ("y", 3)] []
            eval valMap z1 `shouldBe` (1 :+ 3)
            eval valMap z2 `shouldBe` (3 :+ 1)
            eval valMap (y * z1) `shouldBe` 3 P.* (1 :+ 3)
            eval valMap (z2 * z1) `shouldBe` (3 :+ 1) P.* (1 :+ 3)
            eval valMap (z2 .* z1) `shouldBe` (3 :+ 1) P.* (1 :+ 3)
            eval valMap (z2 <.> z1) `shouldBe` (3 :+ 1) P.* (1 :+ 3)
        specify "real 1d 01" $ do
            let x = var1d 4 "x"
                y = var1d 4 "y"
                s = var "s"
                valMap = subs [("s", 3)] [("x", listArray (0, 3) [1, 1, 1, 1]), ("y", listArray (0, 3) [1, 1, 1, 1])]
            elems (eval valMap x) `shouldBe` [1, 1, 1, 1]
            elems (eval valMap (x + y)) `shouldBe` [2, 2, 2, 2]
            elems (eval valMap (s * x)) `shouldBe` [3, 3, 3, 3]
            elems (eval valMap (x .* y)) `shouldBe` [1, 1, 1, 1]
            eval valMap (x <.> y) `shouldBe` 4
        specify "complex 1d 01" $ do
            let x = var1d 4 "x"
                y = var1d 4 "y"
                s = var "s"
                valMap = subs [("s", 3)] [("x", listArray (0, 3) [1, 1, 1, 1]), ("y", listArray (0, 3) [2, 2, 2, 2])]
                z = x +: y
                listZ = [1 :+ 2, 1 :+ 2, 1 :+ 2, 1 :+ 2]
            elems (eval valMap z) `shouldBe` listZ
            elems (eval valMap (z + z)) `shouldBe` zipWith (P.+) listZ listZ
            elems (eval valMap (s * z)) `shouldBe` map (P.* 3) listZ
            elems (eval valMap (z .* z)) `shouldBe` zipWith (P.*) listZ listZ
            eval valMap (z <.> z) `shouldBe` P.sum (zipWith (P.*) listZ listZ)

max_stars_repo_name: dandoh/NeatExpression

max_issues_repo_name: dandoh/NeatExpression

max_forks_repo_name: dandoh/NeatExpression


text:
import Test.Tasty
import Test.Tasty.HUnit

import qualified Data.Vector.Storable.Mutable as MV
import qualified Data.Vector.Storable as V
import Data.Vector.Storable ((!))
import Numeric.LinearAlgebra.HArpack
import Control.Monad
import Foreign.C.Types
import System.IO.Unsafe
import Data.Complex

import Data.List
import Data.Ord
import Data.Maybe
import Data.Function (on)

main = defaultMain tests

tests :: TestTree
tests = testGroup "Tests" [unitTests]

unitTests = testGroup "Unit tests"
  [ testCase "Diagonal Matrix LM"   $ diagMatLMCheck   @?= True
  , testCase "Diagonal Matrix SM"   $ diagMatSMCheck   @?= True
  , testCase "Non-Symmetric Sanity" $ complexEigsCheck @?= True ]

diagMat :: SparseMatrix
diagMat = SparseMatrix 1000
            [((n-1,n-1),fromIntegral n) | n <- [1 .. 1000]]

complexEigs :: SparseMatrix
complexEigs = SparseMatrix 7
  [((0,0),3)
  ,((0,1),-2)
  ,((1,0),4)
  ,((1,1),-1)
  ,((2,2),2)
  ,((3,3),1)
  ,((4,4),0.5)
  ,((5,5),0.25)
  ,((6,6),0.125)]

diagMatLMCheck :: Bool
diagMatLMCheck =
  let RReal ePairs' = eigs diagMat (Symmetric SLM) 6
      ePairs = decreasingR ePairs'
      evals = map fst ePairs
      trueEvals = [1000,999..]
      diff = map abs $ zipWith (-) evals trueEvals
  in all (< 1e-10) diff

decreasingC = sortBy (compare `on` (\x -> - realPart (abs (fst x))))

decreasingR = sortBy (compare `on` (\x -> - abs (fst x)))

diagMatSMCheck :: Bool
diagMatSMCheck =
  let RComplex ePairs = eigs diagMat (NonSymmetric NSM) 6
      evals = map fst ePairs
      trueEvals = map (:+ 0) [1,2..]
      diff = map (realPart . abs) $ zipWith (-) evals trueEvals
  in all (< 1e-10) diff

complexEigsCheck :: Bool
complexEigsCheck =
  let RComplex ePairs' = eigs complexEigs (NonSymmetric NLM) 5
      ePairs = decreasingC ePairs'
      evals = map fst ePairs
      trueEvals = [1 :+ 2, 1 :+ (-2), 2, 1, 0.5]
      diff = map (realPart . abs) $ zipWith (-) evals trueEvals
  in all (< 1e-10) diff

max_stars_repo_name: guaraqe/harpack

max_issues_repo_name: guaraqe/harpack

max_forks_repo_name: guaraqe/harpack


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

module School.Unit.Test.RecLin
( recLinTest ) where

import Control.Monad.IO.Class (liftIO)
import Numeric.LinearAlgebra ((><), fromLists, size, toLists)
import School.TestUtils (diffInput, matIndexes, randomMatrix)
import School.Types.FloatEq ((~=))
import School.Unit.RecLin
import School.Unit.Unit (Unit(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitGradient (UnitGradient(..))
import School.Unit.UnitParams (UnitParams(..))
import School.Utils.LinearAlgebra (compareDoubleMatrix, oneMatrix)
import Test.Tasty (TestTree)
import Test.Tasty.QuickCheck hiding ((><), scale)
import Test.Tasty.TH
import Test.QuickCheck.Monadic (assert, monadicIO, pre)

eps :: Double
eps = 1e-5

prec :: Double
prec = 1e-4

prop_apply :: (Positive Int) -> (Positive Int) -> Property
prop_apply (Positive bSize) (Positive fSize) = monadicIO $ do
  m <- liftIO $ randomMatrix bSize fSize
  let input = BatchActivation m
  let params = EmptyParams
  let (BatchActivation result) = apply recLin params input
  let check = fromLists $ map (map (max 0)) $ toLists m
  assert $ result ~= check

prop_gradient_dims :: (Positive Int) -> Positive Int -> Property
prop_gradient_dims (Positive bSize) (Positive fSize) = monadicIO $ do
  input <- liftIO $ fmap BatchActivation $ randomMatrix bSize fSize
  gradIn <- liftIO $ fmap BatchGradient $ randomMatrix bSize fSize
  let params = EmptyParams
  let (BatchGradient gradOut, _) = deriv recLin params gradIn input
  assert $ (size gradOut == (bSize, fSize))

prop_numerical_gradient :: (Positive Int) -> Positive Int -> Property
prop_numerical_gradient (Positive bSize) (Positive fSize) = monadicIO $ do
  pre $ bSize < 25 && fSize < 25
  input <- liftIO $ fmap BatchActivation $ randomMatrix bSize fSize
  let params = EmptyParams
  let idxs = matIndexes bSize fSize
  let num = map (\idx -> diffInput recLin params input eps idx)
                idxs
  let numGrad = (bSize >< fSize) num
  let inGrad = BatchGradient $ oneMatrix bSize fSize
  let (BatchGradient result, _) = deriv recLin params inGrad input
  assert $ compareDoubleMatrix prec result numGrad

recLinTest :: TestTree
recLinTest = $(testGroupGenerator)

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
module Noether.Lemmata.Prelude
  ( module X
  , fromInteger
  , fromString
  ) where

import           Data.Complex as X
import           Data.Monoid  as X ((<>))
import           Prelude      as X hiding (Eq, Monoid, fromInteger, negate,
                                    recip, (&&), (*), (+), (-), (/), (==), (||))

import           Data.Ratio   as X
import           Data.Int     as X
import qualified Data.String  as S
import qualified Prelude      as P

fromInteger :: Num a => Integer -> a
fromInteger = P.fromInteger

fromString :: S.IsString a => String -> a
fromString = S.fromString

max_stars_repo_name: evertedsphere/noether

max_issues_repo_name: evertedsphere/noether

max_forks_repo_name: evertedsphere/noether


text:
module Main (main) where

import Data.Complex (Complex(..))
import Data.Functor.Classes (readPrec1, showsPrec1)
import Text.ParserCombinators.ReadPrec (readPrec_to_S)
import Text.Read (Read(..))

comp :: Complex Int
comp = 1 :+ 1

compareInstances :: Int -> IO ()
compareInstances p = do
  let precBanner = " (at precedence " ++ show p ++ ")"
  putStrLn $ "Read vs. Read1" ++ precBanner
  print (readPrec_to_S readPrec  p "1 :+ 1" :: [(Complex Int, String)])
  print (readPrec_to_S readPrec1 p "1 :+ 1" :: [(Complex Int, String)])
  putStrLn ""
  putStrLn $ "Show vs. Show1" ++ precBanner
  putStrLn $ showsPrec  p comp ""
  putStrLn $ showsPrec1 p comp ""
  putStrLn ""

main :: IO ()
main = do
  compareInstances 6
  compareInstances 7

max_stars_repo_name: JakobBruenker/ghc

max_issues_repo_name: JakobBruenker/ghc

max_forks_repo_name: JakobBruenker/ghc


text:
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE EmptyCase #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveGeneric #-}

-- {-# OPTIONS_GHC -O -ddump-rule-firings #-}

-- | Provides the class `ParamSet` which is used to represent the set of
-- parameters of a particular model.  The goal of SGD is then to find the
-- parameter values which minimize a given objective function.

module Numeric.SGD.ParamSet
  (
  -- * Class
    ParamSet(..)
  -- * Generics
  , GPMap
  , GAdd
  , GSub
  , GDiv
  , GMul
  , GNorm2
  ) where

import           GHC.Generics
import           GHC.TypeNats (KnownNat)

import           Prelude hiding (div)

import qualified Data.Map.Strict as M

import qualified Numeric.LinearAlgebra.Static as LA

-- | Class of types that can be treated as parameter sets.  It provides basic
-- element-wise operations (addition, multiplication, mapping) which are
-- required to perform stochastic gradient descent.  Many of the operations
-- (`add`, `mul`, `sub`, `div`, etc.) have the same interpretation and follow
-- the same laws (e.g. associativity) as the corresponding operations in `Num`
-- and `Fractional`.
--
-- `zero` takes a parameter set as argument and "zero out"'s all its elements
-- (as in the backprop library).  This allows instances for `Maybe`, `M.Map`,
-- etc., where the structure of the parameter set is dynamic.  This leads to
-- the following property:
--
--     @add (zero x) x = x@
--
-- However, `zero` does not have to obey @(add (zero x) y = y)@.
--
-- A `ParamSet` can be also seen as a (structured) vector, hence `pmap` and
-- `norm_2`.  The latter is not strictly necessary to perform SGD, but it is
-- useful to control the training process.
--
-- `pmap` should obey the following law:
--
--     @pmap id x = x@
--
-- If you leave the body of an instance declaration blank, GHC Generics will be
-- used to derive instances if the type has a single constructor and each field
-- is an instance of `ParamSet`.
class ParamSet a where
  -- | Element-wise mapping
  pmap :: (Double -> Double) -> a -> a

  -- | Zero-out all elements
  zero :: a -> a
  zero = pmap (const 0.0)

--   -- | Element-wise negation
--   neg :: a -> a
--   neg = pmap (\x -> -x)

  -- | Element-wise addition
  add :: a -> a -> a
  -- | Elementi-wise substruction
  sub :: a -> a -> a

  -- | Element-wise multiplication
  mul :: a -> a -> a
  -- | Element-wise division
  div :: a -> a -> a

  -- | L2 norm
  norm_2 :: a -> Double

--   default zero :: (Generic a, GZero (Rep a)) => a -> a
--   zero = genericZero
--   {-# INLINE zero #-}

  default pmap
    :: (Generic a, GPMap (Rep a))
    => (Double -> Double) -> a -> a
  pmap = genericPMap
  {-# INLINE pmap #-}

  default add :: (Generic a, GAdd (Rep a)) => a -> a -> a
  add = genericAdd
  {-# INLINE add #-}

  default sub :: (Generic a, GSub (Rep a)) => a -> a -> a
  sub = genericSub
  {-# INLINE sub #-}

  default mul :: (Generic a, GMul (Rep a)) => a -> a -> a
  mul = genericMul
  {-# INLINE mul #-}

  default div :: (Generic a, GDiv (Rep a)) => a -> a -> a
  div = genericDiv
  {-# INLINE div #-}

  default norm_2 :: (Generic a, GNorm2 (Rep a)) => a -> Double
  norm_2 = genericNorm2
  {-# INLINE norm_2 #-}

{-# RULES
"ParamSet pmap/pmap" forall f g p. pmap f (pmap g p) = pmap (f . g) p
  #-}

-- {-# RULES
-- "ParamSet pmap/add/pmap" forall f g p h q.
--   pmap f (add (pmap g p) (pmap h q))
--   = add (pmap (f . g) p) (pmap (f . h) q)
--   #-}

-- -- | 'add' using GHC Generics; works if all fields are instances of
-- -- 'ParamSet', but only for values with single constructors.
-- genericZero :: (Generic a, GZero (Rep a)) => a -> a
-- genericZero x = to $ gzero (from x)
-- {-# INLINE genericZero #-}

-- | 'add' using GHC Generics; works if all fields are instances of
-- 'ParamSet', but only for values with single constructors.
genericAdd :: (Generic a, GAdd (Rep a)) => a -> a -> a
genericAdd x y = to $ gadd (from x) (from y)
{-# INLINE genericAdd #-}

-- | 'sub' using GHC Generics; works if all fields are instances of
-- 'ParamSet', but only for values with single constructors.
genericSub :: (Generic a, GSub (Rep a)) => a -> a -> a
genericSub x y = to $ gsub (from x) (from y)
{-# INLINE genericSub #-}

-- | 'div' using GHC Generics; works if all fields are instances of
-- 'ParamSet', but only for values with single constructors.
genericDiv :: (Generic a, GDiv (Rep a)) => a -> a -> a
genericDiv x y = to $ gdiv (from x) (from y)
{-# INLINE genericDiv #-}

-- | 'mul' using GHC Generics; works if all fields are instances of
-- 'ParamSet', but only for values with single constructors.
genericMul :: (Generic a, GMul (Rep a)) => a -> a -> a
genericMul x y = to $ gmul (from x) (from y)
{-# INLINE genericMul #-}

-- | 'norm_2' using GHC Generics; works if all fields are instances of
-- 'ParamSet', but only for values with single constructors.
genericNorm2 :: (Generic a, GNorm2 (Rep a)) => a -> Double
genericNorm2 x = gnorm_2 (from x)
{-# INLINE genericNorm2 #-}

-- | 'pmap' using GHC Generics; works if all fields are instances of
-- 'ParamSet', but only for values with single constructors.
genericPMap :: (Generic a, GPMap (Rep a)) => (Double -> Double) -> a -> a
genericPMap f x = to $ gpmap f (from x)
{-# INLINE genericPMap #-}

--------------------------------------------------
-- Generics
--
-- Partially borrowed from the backprop library
--------------------------------------------------

-- -- | Helper class for automatically deriving 'add' using GHC Generics.
-- class GZero f where
--     gzero :: f t -> f t
--
-- instance ParamSet p => GZero (K1 i p) where
--     gzero (K1 x) = K1 (zero x)
--     {-# INLINE gzero #-}
--
-- instance (GZero f, GZero g) => GZero (f :*: g) where
--     gzero (x1 :*: y1) = x2 :*: y2
--       where
--         !x2 = gzero x1
--         !y2 = gzero y1
--     {-# INLINE gzero #-}
--
-- instance GZero V1 where
--     gzero = \case {}
--     {-# INLINE gzero #-}
--
-- instance GZero U1 where
--     gzero _ = U1
--     {-# INLINE gzero #-}
--
-- instance GZero f => GZero (M1 i c f) where
--     gzero (M1 x) = M1 (gzero x)
--     {-# INLINE gzero #-}
--
-- -- instance GZero f => GZero (f :.: g) where
-- --     gzero = Comp1 gzero
-- --     {-# INLINE gzero #-}

-- | Helper class for automatically deriving 'add' using GHC Generics.
class GAdd f where
    gadd :: f t -> f t -> f t

instance ParamSet a => GAdd (K1 i a) where
    gadd (K1 x) (K1 y) = K1 (add x y)
    {-# INLINE gadd #-}

instance (GAdd f, GAdd g) => GAdd (f :*: g) where
    gadd (x1 :*: y1) (x2 :*: y2) = x3 :*: y3
      where
        !x3 = gadd x1 x2
        !y3 = gadd y1 y2
    {-# INLINE gadd #-}

instance GAdd V1 where
    gadd = \case {}
    {-# INLINE gadd #-}

instance GAdd U1 where
    gadd _ _ = U1
    {-# INLINE gadd #-}

instance GAdd f => GAdd (M1 i c f) where
    gadd (M1 x) (M1 y) = M1 (gadd x y)
    {-# INLINE gadd #-}

-- instance GAdd f => GAdd (f :.: g) where
--     gadd (Comp1 x) (Comp1 y) = Comp1 (gadd x y)
--     {-# INLINE gadd #-}

-- | Helper class for automatically deriving 'sub' using GHC Generics.
class GSub f where
    gsub :: f t -> f t -> f t

instance ParamSet a => GSub (K1 i a) where
    gsub (K1 x) (K1 y) = K1 (sub x y)
    {-# INLINE gsub #-}

instance (GSub f, GSub g) => GSub (f :*: g) where
    gsub (x1 :*: y1) (x2 :*: y2) = x3 :*: y3
      where
        !x3 = gsub x1 x2
        !y3 = gsub y1 y2
    {-# INLINE gsub #-}

instance GSub V1 where
    gsub = \case {}
    {-# INLINE gsub #-}

instance GSub U1 where
    gsub _ _ = U1
    {-# INLINE gsub #-}

instance GSub f => GSub (M1 i c f) where
    gsub (M1 x) (M1 y) = M1 (gsub x y)
    {-# INLINE gsub #-}

-- instance GSub f => GSub (f :.: g) where
--     gsub (Comp1 x) (Comp1 y) = Comp1 (gsub x y)
--     {-# INLINE gsub #-}

-- | Helper class for automatically deriving 'mul' using GHC Generics.
class GMul f where
    gmul :: f t -> f t -> f t

instance ParamSet a => GMul (K1 i a) where
    gmul (K1 x) (K1 y) = K1 (mul x y)
    {-# INLINE gmul #-}

instance (GMul f, GMul g) => GMul (f :*: g) where
    gmul (x1 :*: y1) (x2 :*: y2) = x3 :*: y3
      where
        !x3 = gmul x1 x2
        !y3 = gmul y1 y2
    {-# INLINE gmul #-}

instance GMul V1 where
    gmul = \case {}
    {-# INLINE gmul #-}

instance GMul U1 where
    gmul _ _ = U1
    {-# INLINE gmul #-}

instance GMul f => GMul (M1 i c f) where
    gmul (M1 x) (M1 y) = M1 (gmul x y)
    {-# INLINE gmul #-}

-- instance GMul f => GMul (f :.: g) where
--     gmul (Comp1 x) (Comp1 y) = Comp1 (gmul x y)
--     {-# INLINE gmul #-}

-- | Helper class for automatically deriving 'div' using GHC Generics.
class GDiv f where
    gdiv :: f t -> f t -> f t

instance ParamSet a => GDiv (K1 i a) where
    gdiv (K1 x) (K1 y) = K1 (div x y)
    {-# INLINE gdiv #-}

instance (GDiv f, GDiv g) => GDiv (f :*: g) where
    gdiv (x1 :*: y1) (x2 :*: y2) = x3 :*: y3
      where
        !x3 = gdiv x1 x2
        !y3 = gdiv y1 y2
    {-# INLINE gdiv #-}

instance GDiv V1 where
    gdiv = \case {}
    {-# INLINE gdiv #-}

instance GDiv U1 where
    gdiv _ _ = U1
    {-# INLINE gdiv #-}

instance GDiv f => GDiv (M1 i c f) where
    gdiv (M1 x) (M1 y) = M1 (gdiv x y)
    {-# INLINE gdiv #-}

-- instance GDiv f => GDiv (f :.: g) where
--     gdiv (Comp1 x) (Comp1 y) = Comp1 (gdiv x y)
--     {-# INLINE gdiv #-}

-- | Helper class for automatically deriving 'norm_2' using GHC Generics.
class GNorm2 f where
    gnorm_2 :: f t -> Double

instance ParamSet a => GNorm2 (K1 i a) where
    gnorm_2 (K1 x) = norm_2 x
    {-# INLINE gnorm_2 #-}

instance (GNorm2 f, GNorm2 g) => GNorm2 (f :*: g) where
    gnorm_2 (x1 :*: y1) =
      sqrt ((x2 ^ (2 :: Int)) + (y2 ^ (2 :: Int)))
      where
        !x2 = gnorm_2 x1
        !y2 = gnorm_2 y1
    {-# INLINE gnorm_2 #-}

instance GNorm2 V1 where
    gnorm_2 = \case {}
    {-# INLINE gnorm_2 #-}

instance GNorm2 U1 where
    gnorm_2 _ = 0
    {-# INLINE gnorm_2 #-}

instance GNorm2 f => GNorm2 (M1 i c f) where
    gnorm_2 (M1 x) = gnorm_2 x
    {-# INLINE gnorm_2 #-}

-- -- TODO: Make sure this makes sense
-- instance GNorm2 f => GNorm2 (f :.: g) where
--     gnorm_2 (Comp1 x) = gnorm_2 x
--     {-# INLINE gnorm_2 #-}

-- | Helper class for automatically deriving 'pmap' using GHC Generics.
class GPMap f where
    gpmap :: (Double -> Double) -> f t -> f t

instance ParamSet a => GPMap (K1 i a) where
    gpmap f (K1 x) = K1 (pmap f x)
    {-# INLINE gpmap #-}

instance (GPMap f, GPMap g) => GPMap (f :*: g) where
    gpmap f (x1 :*: y1) = x2 :*: y2
      where
        !x2 = gpmap f x1
        !y2 = gpmap f y1
    {-# INLINE gpmap #-}

instance GPMap V1 where
    gpmap _ = \case {}
    {-# INLINE gpmap #-}

instance GPMap U1 where
    gpmap _ _ = U1
    {-# INLINE gpmap #-}

instance GPMap f => GPMap (M1 i c f) where
    gpmap f (M1 x) = M1 (gpmap f x)
    {-# INLINE gpmap #-}

-- instance GPMap f => GPMap (f :.: g) where
--     gpmap f (Comp1 x) = Comp1 (gpmap f x)
--     {-# INLINE gpmap #-}

--------------------------------------------------
-- Basic instances
--------------------------------------------------

instance ParamSet Double where
  zero = const 0
  pmap = id
  add = (+)
  sub = (-)
  mul = (*)
  div = (/)
  norm_2 = abs

instance (ParamSet a, ParamSet b) => ParamSet (a, b) where
  pmap f (x, y) = (pmap f x, pmap f y)
  add (x1, y1) (x2, y2) = (x1 `add` x2, y1 `add` y2)
  sub (x1, y1) (x2, y2) = (x1 `sub` x2, y1 `sub` y2)
  mul (x1, y1) (x2, y2) = (x1 `mul` x2, y1 `mul` y2)
  div (x1, y1) (x2, y2) = (x1 `div` x2, y1 `div` y2)
  norm_2 (x, y)
    = sqrt . sum . map ((^(2::Int)))
    $ [norm_2 x, norm_2 y]

instance (KnownNat n) => ParamSet (LA.R n) where
  zero = const 0
  pmap = LA.dvmap
  add = (+)
  sub = (-)
  mul = (*)
  div = (/)
  norm_2 = LA.norm_2

instance (KnownNat n, KnownNat m) => ParamSet (LA.L n m) where
  zero = const 0
  pmap = LA.dmmap
  add = (+)
  sub = (-)
  mul = (*)
  div = (/)
  norm_2 = LA.norm_2

-- | `Nothing` represents a deactivated parameter set component. If `Nothing`
-- is given as an argument to one of the `ParamSet` operations, the result is
-- `Nothing` as well.
--
-- This differs from the corresponding instance in the backprop library, where
-- `Nothing` is equivalent to `Just 0`.  However, the implementation below
-- seems to correspond adequately enough to the notion that a particular
-- component is either active or not in both the parameter set and the
-- gradient, hence it doesn't make sense to combine `Just` with `Nothing`.
instance (ParamSet a) => ParamSet (Maybe a) where
  zero = fmap zero
  pmap = fmap . pmap

  add (Just x) (Just y) = Just (add x y)
  add _ _ = Nothing

  sub (Just x) (Just y) = Just (sub x y)
  sub _ _ = Nothing

  mul (Just x) (Just y) = Just (mul x y)
  mul _ _ = Nothing

  div (Just x) (Just y) = Just (div x y)
  div _ _ = Nothing

  norm_2 = maybe 0 norm_2

-- | A map with different parameter sets (of the same type) assigned to the
-- individual keys.
--
-- When combining two maps with different sets of keys, only their intersection
-- is preserved.
instance (Ord k, ParamSet a) => ParamSet (M.Map k a) where
  zero = fmap zero
  pmap f = fmap (pmap f)
  add = M.intersectionWith add
  sub = M.intersectionWith sub
  mul= M.intersectionWith mul
  div= M.intersectionWith div
  norm_2 = sqrt . sum . map ((^(2::Int)) . norm_2)  . M.elems

max_stars_repo_name: kawu/sgd

max_issues_repo_name: kawu/sgd

max_forks_repo_name: kawu/sgd


text:
import Test.QuickCheck
import Execution
import InputOutput
import Types
import Data.List.Split
import System.Random
import System.Random.Shuffle
import Numeric.LinearAlgebra.HMatrix
import Numeric.LinearAlgebra.Data
import Training
import Execution

-- Exemplo de como usar:
--
-- $ ghci Tests.hs
--  > import Test.QuickCheck
--  > quickCheck prop_reverseReverse

prop_reverseReverse :: [Int] -> Bool
prop_reverseReverse xs = reverse (reverse xs) == xs


text:
{-# LANGUAGE FlexibleContexts #-}
module LearnHaskDA.Chapter8 where
import Numeric.LinearAlgebra as LA
import Data.List as L
import Data.HashMap.Strict as HM
import qualified LearnHaskDA.Chapter4 as Ch4
import qualified LearnHaskDA.Chapter6 as Ch6
import qualified LearnHaskDA.Chapter7 as Ch7 (Counts, MyWord, clean, frequency, wordFrequencyByFst)

import Test.QuickCheck -- This is not imported by the book (and that is a bad sign: the book talk almost nothing about code testing)
import Data.Function (on)

test140 = do
    tweetsEnglish <- Ch4.queryDatabase "tweets.sql" "SELECT user, message FROM tweets WHERE language='en'"
    let tweets = Ch4.readTwoColumns 0 1 tweetsEnglish :: [(String, String)]
    let freqTable = Ch7.frequency tweets
    print $ HM.size freqTable
    let uniqueTweets = keys freqTable
    let cleanedTweets = L.map (\(user, mess) -> (user, Ch7.clean mess)) uniqueTweets
    let wordFrequency = Ch7.frequency $ concatMap snd cleanedTweets
    return (cleanedTweets, wordFrequency)

test141 = do
    (cleanedTweets, wordFrequency) <- test140
    let topWords = take 25 $ sortBy (\(_, c1) (_, c2) -> compare c2 c1) $ HM.toList wordFrequency
    return topWords

test142 = do
    (cleanedTweets, wordFrequency) <- test140
    let notStopWords = L.filter (\(word, _) -> notElem word stopWords) (HM.toList wordFrequency)
    let topWords = take 25 . L.map fst $ sortBy (\(_, c1) (_, c2) -> compare c2 c1) notStopWords
    return topWords

stopWords :: [String]
stopWords = ["", "a", "about", "above", "above", "across", "after",
  "afterwards", "again", "against", "all", "almost", "alone", "along",
  "already", "also", "although", "always", "am", "among", "amongst",
  "amoungst", "amount", "an", "and", "another", "any", "anyhow",
  "anyone", "anything", "anyway", "anywhere", "are", "around", "as",
  "at", "back", "be", "became", "because", "become", "becomes",
  "becoming", "been", "before", "beforehand", "behind", "being",
  "below", "beside", "besides", "between", "beyond", "bill", "both",
  "bottom", "but", "by", "call", "can", "cannot", "cant", "co", "con",
  "could", "couldnt", "cry", "de", "describe", "detail", "do", "done",
  "dont", "down", "due", "during", "each", "eg", "eight", "either",
  "eleven", "else", "elsewhere", "empty", "enough", "etc", "even",
  "ever", "every", "everyone", "everything", "everywhere", "except",
  "few", "fifteen", "fifty", "fill", "find", "fire", "first", "five",
  "for", "former", "formerly", "forty", "found", "four", "from",
  "front", "full", "further", "get", "give", "go", "got", "had", "has",
  "hasnt", "have", "he", "hence", "her", "here", "hereafter", "hereby",
  "herein", "hereupon", "hers", "herself", "him", "himself", "his",
  "how", "however", "https", "hundred", "i", "ie", "if", "im", "in", "inc",
  "indeed", "interest", "into", "is", "it", "its", "itself", "just",
  "keep", "last", "latter", "latterly", "least", "less", "ltd", "made",
  "many", "may", "me", "meanwhile", "might", "mill", "mine", "more",
  "moreover", "most", "mostly", "move", "much", "must", "my", "myself",
  "name", "namely", "neither", "need", "never", "nevertheless",
  "next", "nine", "no", "nobody", "none", "noone", "nor", "not",
  "nothing", "now", "nowhere", "of", "off", "often", "on", "once",
  "one", "only", "onto", "or", "other", "others", "otherwise", "our",
  "ours", "ourselves", "out", "over", "own", "part", "per", "perhaps",
  "please", "put", "rather", "re", "same", "see", "seem", "seemed",
  "seeming", "seems", "serious", "several", "she", "should", "show",
  "side", "since", "sincere", "six", "sixty", "so", "some", "somehow",
  "someone", "something", "sometime", "sometimes", "somewhere", "still",
  "such", "system", "take", "ten", "than", "that", "the", "their",
  "them", "themselves", "then", "thence", "there", "thereafter",
  "thereby", "therefore", "therein", "thereupon", "these", "they",
  "thick", "thin", "third", "this", "those", "though", "three",
  "through", "throughout", "thru", "thus", "to", "together", "too",
  "top", "toward", "towards", "twelve", "twenty", "two", "un", "under",
  "until", "up", "upon", "us", "very", "via", "want", "was", "we",
  "well", "were", "what", "whatever", "when", "whence", "whenever",
  "where", "whereafter", "whereas", "whereby", "wherein", "whereupon",
  "wherever", "whether", "which", "while", "whither", "who", "whoever",
  "whole", "whom", "whose", "why", "will", "with", "within", "without",
  "would", "yet", "you", "your", "youre", "yours", "yourself",
  "yourselves", "rt", "mt", "u"]

bookTopWords = ["like", "amp", "day", "good", "new", "love", "time", "follow", "great",
  "today", "make", "lot", "people", "video", "know", "life", "happy",
  "look", "think", "girl", "win", "photo", "way", "little", "really"]

myTopWords = ["amp","day","new","like","love","good","great","video",
              "make","win","follow","life","pill","today","happy","time",
              "chance","photo","posted","birthday","lot","brain","know",
              "limitless","think"]

topWords = bookTopWords

-- page 143
type UserName = String
type UserWordCount = [[Double]]
test143 :: IO (HashMap UserName (Ch7.Counts Ch7.MyWord), UserWordCount)
test143 = do
    (cleanedTweets, wordFrequency) <- test140
    let wordFrequencyByUser = Ch7.wordFrequencyByFst cleanedTweets
    let getFrequencyOfWordByUser user word = fromIntegral $ (HM.lookupDefault 0 word (wordFrequencyByUser HM.! user)) :: Double
    let allUsers = keys wordFrequencyByUser
    let wordMatrix = L.map (\user -> L.map (getFrequencyOfWordByUser user) topWords) allUsers
    return (wordFrequencyByUser, wordMatrix)

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

test150 = do
    homeRecord <- Ch4.queryDatabase "winloss.sql" "SELECT homeTeam, SUM(homescore > awayscore), SUM(homescore), COUNT(*) FROM winloss GROUP BY homeTeam;"
    awayRecord <- Ch4.queryDatabase "winloss.sql" "SELECT awayTeam, SUM(awayscore > homescore), SUM(awayscore), COUNT(*) FROM winloss GROUP BY awayTeam;"
    let totalWins     = zipWith (+) (Ch4.readDoubleColumn homeRecord 1) (Ch4.readDoubleColumn awayRecord 1)
    let totalRuns     = zipWith (+) (Ch4.readDoubleColumn homeRecord 2) (Ch4.readDoubleColumn awayRecord 2)
    let totalGames    = zipWith (+) (Ch4.readDoubleColumn homeRecord 3) (Ch4.readDoubleColumn awayRecord 3)
    let winPercentage = zipWith (/) totalWins totalGames
    let runsPerGame   = zipWith (/) totalRuns totalGames
    let baseball      = L.map (\(a,b) -> [a,b]) $ zip winPercentage runsPerGame
    let p@(baseballMean, baseballCovMatrix) = meanCov $ fromLists baseball
    return p

-- eigSH :: Field t => Herm t -> (Vector Double, Matrix t)
test151 = do
  (baseballMean, baseballCovMatrix) <- test150
  let p@(baseballEvalues, baseballEvectors) = eigSH (sym baseballCovMatrix)
  return p

-- page 154
-- principalComponentAnalysis :: [[Double]] -> Integer -> Matrix Double
principalComponentAnalysis records top = featureVectors <> topOrderedEvectors
-- Some unnecessary tr in the book:
--    tr $ mul (tr topOrderedEvectors) (tr featureVectors)
  where  featureVectors  = fromLists records
         (_, covMatrix)  = meanCov featureVectors
         (_, evectors)   = eigSH (trustSym covMatrix)
         topOrderedEvectors = takeColumns (fromInteger top) evectors

{-
-- Done: pull request for hmatrix changing the type of meanCov to use Herm
https://github.com/albertoruiz/hmatrix/pull/157
https://github.com/albertoruiz/hmatrix/blob/master/packages/base/src/Internal/Container.hs

meanCovH :: Matrix Double -> (Vector Double, Herm Double)
meanCovH x = (med, cov) where
    r    = rows x
    k    = 1 / fromIntegral r
    med  = konst k r <> x
    meds = konst 1 r `outer` med
    xc   = x - meds
    cov  = scale (recip (fromIntegral (r-1))) (mTm xc)
-}

-- let pcaMatrix = principalComponentAnalysis wordMatrix 5

-- euclideanDistanceSqrdFromRow :: Matrix Double -> Int -> Vector Double
euclideanDistanceSqrdFromRow :: (Num t, Num (Vector t), Container Vector t) => Matrix t -> Int -> Vector t
euclideanDistanceSqrdFromRow records row = sumOfSquaresVector
  where
    d = cols records
    n = rows records
    theRow = subMatrix (row, 0) (1, d) records
    copyMatrix = repmat theRow n 1
    diffMatrix      = records    - copyMatrix  -- Pointwise operation
    diffSqrdMatrix  = diffMatrix * diffMatrix  -- Pointwise operation
    sumOfSquaresVector = sum $ toColumns diffSqrdMatrix

orderBookSimplified :: Eq a => [a] -> [a] -> [Int]
orderBookSimplified unsorted = concatMap (`L.elemIndices` unsorted)

test156' :: Ord a => [a] -> [Int]
test156' xs = orderBookSimplified xs (nub $ sort xs)

-- Much faster (better complexity):
test156 :: Ord a => [a] -> [Int]
test156 = L.map snd . sort . (`zip` [0..])

qc156 xs = test156 xs == test156' xs

findClosestFeaturesToRow :: Matrix Double -> Int -> Int -> [Int]
findClosestFeaturesToRow records row knn = take knn orderOfFeatures
  where  orderOfFeatures = test156 sumOfSquares
         sumOfSquares    = LA.toList $ euclideanDistanceSqrdFromRow records row

-- The users "most active" in the common words: [8225,15747,556, 1044, 1517, 2133]

test157 = do
  (tab, wordMatrix) <- test143
  let pcaMatrix = principalComponentAnalysis wordMatrix 5
  return (tab, findClosestFeaturesToRow pcaMatrix 8225 6)

----------------------------------------------------------------
{-
  Experiment to locate a few users with many words in the tweets database:
  (tab, close) <- test157
  let s = HM.map HM.size tab
  sortBy (compare `on` snd) (HM.toList s)

Top users (at one point):

("WindowMouse",60),("kellyboy1977",62),("cassidyhanft",72),("love_the_rock2",80),("Gvccizay",81),("revraamess",97)]

-}

max_stars_repo_name: patrikja/LearningHaskellDataAnalysis

max_issues_repo_name: patrikja/LearningHaskellDataAnalysis

max_forks_repo_name: patrikja/LearningHaskellDataAnalysis


text:
-- |
-- Module      : Southpaw.Picasso.Render
-- Description : Vector rendering in 2D with Cairo
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created September 1 2015

-- TODO | - Rename (eg. Draw, Shapes) (?)
--        - Text, typography (find good library)
--        - Debug versions (cf. Occlusion.Render)
--        - Use BoundingBox (?)

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC directives
--------------------------------------------------------------------------------------------------------------------------------------------
{-# LANGUAGE TupleSections #-}

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module Southpaw.Picasso.Render where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.Complex
import Control.Monad (forM, forM_, liftM, liftM2, void)

import Control.Monad.IO.Class

import qualified Graphics.Rendering.Cairo                         as Cairo
-- import qualified Graphics.Rendering.Cairo.Internal.Surfaces.Image as Image

import qualified Southpaw.Picasso.Palette     as Palette
import qualified Southpaw.Picasso.Shapes      as Shapes
import           Southpaw.Picasso.RenderUtils
import           Southpaw.Math.Constants
import           Southpaw.Cartesian.Plane.Utilities

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------
-- General ---------------------------------------------------------------------------------------------------------------------------------
-- |
-- TODO: Rename (eg. plot)
-- TODO: Generalise (not just circles)
trail :: Palette.Colour Double -> [Complex Double] -> Cairo.Render ()
trail fill trail = forM_ trail $ \dot -> do
    choose fill
    circle dot 3

-- |
grid :: Int -> Int -> Double -> Cairo.Render ()
grid cols rows size = do
    -- TODO: Figure out how to use fill AND stroke
    Cairo.setLineWidth 4
    gridM_ cols rows $ \ cl rw -> tilePath cl rw >> Cairo.fill   --
    gridM_ cols rows $ \ cl rw -> tilePath cl rw >> Cairo.stroke --
    where
      chooseColour cl rw = if (cl `mod` 2) == (rw `mod` 2) then 0.3 else 0.75 -- TODO: This should be a utility function
      tilePath cl rw     = Cairo.rectangle (fromIntegral cl*size) (fromIntegral rw*size) size size >> Cairo.setSourceRGBA 0.22 0.81 (chooseColour cl rw) 0.32

-- Primitives ------------------------------------------------------------------------------------------------------------------------------

-- |
line :: Complex Double -> Complex Double -> Cairo.Render ()
line (fr:+om) (t:+o) = Cairo.moveTo fr om >> Cairo.lineTo t o

-- | Renders a path of connected lines
-- TODO: Options for colour, width, closed/open, etcairo.
linepath :: [Complex Double] -> Cairo.Render ()
linepath []      = return ()
linepath (e:dge) = void $ vectorise Cairo.moveTo e >> forM dge (vectorise Cairo.lineTo)

-- |
-- TODO: Support asymmetrical crosshairs (?)
crosshairs :: Complex Double -> Complex Double -> Cairo.Render ()
crosshairs centre size = do
  line (centre - (hdx:+0)) (centre + (hdx:+0))
  line (centre - (0:+hdy)) (centre + (0:+hdy))
  where
    (hdx:+hdy) = 0.5*size

-- Shapes ----------------------------------------------------------------------------------------------------------------------------------

-- |
-- TODO: Extract argument conversion logic (centre/size vectors to unpacked left-top/dx/dy)
rectangle :: Complex Double -> Complex Double -> Cairo.Render ()
rectangle (cx:+cy) (dx:+dy) = Cairo.rectangle (cx-dx/2) (cy-dy/2) dx dy

-- | A rectangle with an 'anchor' (in normalised coordinates) that is relative to the centre
-- TODO: Should anchor point be relative to centre or topleft corner, use normalised or absolute coords (?)
-- TODO: Less confusing terminology...
-- TODO: Refactor
-- TODO: Test
anchoredRectangle :: Complex Double -> Complex Double -> Complex Double -> Cairo.Render ()
anchoredRectangle p size anchor = rectangle (p-dotwise (*) (anchor+(0.5:+0.5)) size) size

-- |
-- TODO: Add arguments for colour, stroke, etcairo.
-- TODO: Maybe it'd be better if we stuck to the normal pattern of path-config-action that Cairo follows
-- TODO: Make polymorphic
polygon :: Integral int => int -> Double -> Complex Double -> (Double, Double, Double, Double) -> Bool -> Cairo.Render ()
polygon sides radius origin (r,g,b,a) filled = do
    -- TODO: Refine 'wrap-around logic'
    Cairo.moveTo fx fy
    forM_ rest $ \(x:+y) -> Cairo.lineTo x y

    Cairo.setSourceRGBA r g b a
    Cairo.setLineWidth 12
    if filled
        then Cairo.fill
        else Cairo.stroke
    where
      ((fx:+fy):rest) = Shapes.polygon sides radius origin ++ [fx:+fy]

-- |
-- TODO: Options for fill/stroke, colour, width, etcairo.
circle :: Complex Double -> Double -> Cairo.Render ()
circle (cx:+cy) radius = do
    Cairo.arc cx cy radius 0 τ
    -- Cairo.fill

-- Composite -------------------------------------------------------------------------------------------------------------------------------

-- |
arrow :: Complex Double -> Complex Double -> Double -> Double -> Double -> Cairo.Render ()
arrow from to sl sw hw = do
    let (first:rest) = closePath $ Shapes.arrow from to sl sw hw
    vectorise Cairo.moveTo first
    forM_ rest $ vectorise Cairo.lineTo

-- |
-- Ugh, I hate underscores so much
-- TODO: Make polymorphic
circlearc :: Int -> Complex Double -> Double -> Double -> Double -> Double -> Cairo.Render ()
circlearc
    count    -- Number of small circles
    (ox:+oy) -- Centre of the 'arc' (pixels?)
    spread   -- Radius of the 'arc'
    radius   -- Radius of the small circles
    begin    -- Start angle of the arc
    extent = forM_ [1..count] $ \ n -> do
        let n' = fromIntegral n
        let θ  = begin + n'*extent/fromIntegral count
        Cairo.arc (ox - spread*cos θ) (oy - spread*sin θ) radius 0 τ
        Cairo.setSourceRGBA (0.5 * (1 + sin θ)) (0.1*n') (1/n') 0.95
        Cairo.fill

-- |
bezier :: Complex Double -> Complex Double -> Complex Double -> Cairo.Render ()
bezier (x1:+y1) (x2:+y2) (x3:+y3) = Cairo.curveTo x1 y1 x2 y2 x3 y3

-- |
-- TODO: Generic rounded polygon
roundrect :: Complex Double -> Complex Double -> Double -> Cairo.Render ()
roundrect centre@(cx:+cy) size@(dx:+dy) radius = forM_ (zip [real, imag, real, imag] [(-dx):+(-dy), (dx):+(-dy), (dx):+(dy), (-dx):+(dy)]) $ \(dir, delta@(dx':+dy')) -> do

  -- TODO: Finish refactoring

	--
	-- let dir = (signum (dx*dy))

	-- Line segment
	vectorise Cairo.moveTo (centre + delta + dir radius)
	vectorise Cairo.lineTo (centre - flipx size/2 - real radius)

	-- Curve

	-- -- First line segment
	-- vectorise Cairo.moveTo (centre - size/2       + real radius)
	-- vectorise Cairo.lineTo (centre - flipx size/2 - real radius)

	-- -- Curve
	-- let (cx':+cy') = (centre - flipx size/2 + ((-radius):+radius)) in Cairo.arc cx' cy' radius (3*π/2) (4*π/2)

	-- -- Second line segment
	-- vectorise Cairo.moveTo (centre - flipx size/2  + imag radius)
	-- vectorise Cairo.lineTo (centre + size/2        - imag radius)

	-- -- Curve
	-- let (cx':+cy') = (centre + size/2 - (radius:+radius)) in Cairo.arc cx' cy' radius 0 (π/2)

	-- -- Third line segment
	-- vectorise Cairo.moveTo (centre + size/2       - real radius)
	-- vectorise Cairo.lineTo (centre - flipy size/2 + real radius)

	-- -- Curve
	-- let (cx':+cy') = (centre - flipy size/2 + flipy (radius:+radius)) in Cairo.arc cx' cy' radius (π/2) π

	-- -- Fourth line segment
	-- vectorise Cairo.moveTo (centre - flipy size/2 - imag radius)
	-- vectorise Cairo.lineTo (centre - size/2       + imag radius)

	-- -- Curve
	-- let (cx':+cy') = (centre - size/2 + (radius:+radius)) in Cairo.arc cx' cy' radius π (3*π/2)

-- Images ----------------------------------------------------------------------------------------------------------------------------------

-- |
-- TODO: Factor out clip logic, document properly (cf. other anchored functions in this module)
-- TODO: Wrapper for images (?)
-- TODO: Factor out clip area (?)
anchoredImageWithClip :: (Complex Double -> Complex Double -> Cairo.Render ()) -> Complex Double -> Complex Double -> Cairo.Surface -> Cairo.Render ()
anchoredImageWithClip clip p anchor im = do
  size   <- imageSurfaceSize im
  centre <- return $ p - dotwise (*) anchor size
  clip centre size
  Cairo.clip
  vectorise Cairo.translate $ (centre - size*0.5)
  Cairo.setSourceSurface im 0 0
  -- vectorise Cairo.moveTo (centre - size*0.5)
  Cairo.paint
  vectorise Cairo.translate $ negate (centre - size*0.5)
  Cairo.resetClip

-- |
anchoredImage :: Complex Double -> Complex Double -> Cairo.Surface -> Cairo.Render ()
anchoredImage p anchor im = anchoredImageWithClip rectangle p anchor im

-- |
imageWithClip :: (Complex Double -> Complex Double -> Cairo.Render ()) -> Complex Double -> Cairo.Surface -> Cairo.Render ()
imageWithClip clip centre im = anchoredImageWithClip clip centre (0.0:+0.0) im

-- |
image :: Complex Double -> Cairo.Surface -> Cairo.Render ()
image p im = anchoredImageWithClip rectangle p (0.0:+0.0) im

-- Typography ------------------------------------------------------------------------------------------------------------------------------

-- | Renders the given string (with an arbitrary function 'draw'). The position of the centre of the bounding box is
--   given by the pin point 'p' and an 'anchor' (whose coordinates are normalised with respect to the text bounds).
-- TODO: Extract anchoring logic
-- TODO: Define anchor constants (eg. left:+top, centre:+bottom) (?)
anchoredText :: Complex Double -> Complex Double -> (String -> Cairo.Render a) -> String -> Cairo.Render a
anchoredText p anchor draw text = do
  extents <- textsize text
  vectorise Cairo.moveTo $ p - dotwise (*) (anchor + (0.5:+0.5)) extents
  draw text

-- |
centredText :: Complex Double -> (String -> Cairo.Render a) -> String -> Cairo.Render a
centredText centre draw text = anchoredText centre (0.0:+0.0) draw text

max_stars_repo_name: SwiftsNamesake/Southpaw

max_issues_repo_name: SwiftsNamesake/Southpaw

max_forks_repo_name: SwiftsNamesake/Southpaw


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, DeriveFunctor, MultiParamTypeClasses, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : First attempt at approximate mean
-- | Creator: Xiao Ling
-- | Created: 11/29/2015
-- | see    : https://github.com/snoyberg/conduit
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module Play1129 where

import System.Random

import Data.List
import Data.Random
import Data.Conduit
import qualified Data.Conduit.List as Cl

import Data.Vector (Vector,fromList,cons,snoc,(!?))
import qualified Data.Vector       as V

import Control.Monad.Identity
import Control.Monad.State

import Test.QuickCheck ((==>), (.&&.))
import qualified Test.QuickCheck as T

import Core
import Statistics

{-----------------------------------------------------------------------------
   1. make some examples
   2. make
   3. see if cassava can be put into a conduit
   4. output Row type

   see : http://stackoverflow.com/questions/2376981/haskell-types-frustrating-a-simple-average-function
------------------------------------------------------------------------------}

-- * first figure out how to sample from a list with probability p
-- * it's a dice with 1/p faces

--type MaxSize      = Int
--type Store a      = (Vector a, MaxSize)

--toMedian :: Counter -> Sink a (StateT (Store a) RVar) (Vector a)
--toMedian c = do
--  mi <- await
--  case mi of
--    Nothing -> lift get >>= return . fst
--    Just i  -> do
--      (v,s) <- lift get
--      case V.length v < s of
--        True -> do
--          let v' = i `cons` v
--          put

{-----------------------------------------------------------------------------
    II. Approximate Median try ii
        too much unnnecessary stuff that adds no value
------------------------------------------------------------------------------}

-- * TODO: make this into an array since we'll be doing (!!) which is O(n) in list
-- * A reservoir is a list, its max allowed size, and its current size

-- * this deign is really bad, cannot gaurantee that currentSize is upToDate with length [a]
-- * dependent types would be good, or we can jsut get rid of current size
-- * or hide it in Core
data Reservoir a = R { run :: [a], size :: Int, currentSize :: Int }
  deriving (Show,Functor,Eq)

instance Num a => T.Arbitrary (Reservoir a) where

  -- * by construction, an arbitrary reservoir is always 1 item away from full
  arbitrary = do
    s <- T.arbitrary :: T.Gen (T.Positive Int)
    let s' = T.getPositive s
    let n = s' - 1   :: Int
    return $ R (replicate n 0) s' n

  shrink (R xs s n) = (\s' -> let n = max 0 (s' - 1) in (R (replicate n 0) s' n)) <$> T.shrink s

-- * constructor an empty resoivor of max size `s`
reservoir :: Int -> Reservoir a
reservoir s = R [] s 0

-- * check if reservoir is full
full :: Reservoir a -> Bool
full (R xs s n) = n >= s

-- * insert item, note the order
(>+) :: Reservoir a -> a -> Reservoir a
(>+) t@(R xs s n) x | full t     = t
                    | otherwise  = R (x:xs) s $ succ n

-- * remove the `i`th item from `r`
-- * if i > size then no item removed
(>-) :: Reservoir a -> Int -> Reservoir a
r@(R xs s n) >- i | i > n     = r
                  | otherwise = undefined      -- * really need to use array here

-- * sample an item `x` from the stream with uniform probability place in reservoir
-- * keep a conter `c` of elements seen
sampl' :: Counter -> Sink a (StateT (Reservoir a) RVar) [a]
sampl' i = do
  ma <- await
  case ma of
    Nothing -> lift get >>= \r -> return $ run r
    Just a  -> do
      lift $ store' i a
      sampl' $ succ i

store' :: Counter -> a -> StateT (Reservoir a) RVar ()
store' i a = do
  r <- get
  case full r of
    False -> put $ r >+ a
    _     -> return ()
      where p = (read . show $ size r :: Float)/i

-- * test m/2 - e*m < rank(y) < m/2 + e*m

-- * quickCheck Resoivor

-- * full resoivor is full
prop_fullRes :: Reservoir Int -> T.Property
prop_fullRes r = size r > 0 ==> full (r>+1) == True

---- * inserting into non-full resoivor resutls in additional element
prop_insertNotFull :: Reservoir Int -> T.Property
prop_insertNotFull r = size r > 0 ==> currentSize (r>+1) == currentSize r + 1

-- * inserting into full resoivor leaves it unchanged
prop_insertFull :: Reservoir Int -> T.Property
prop_insertFull r = size r > 0 ==> currentSize (r'>+1) == currentSize r'
  where r' = r >+ 1

props :: T.Property
props =  prop_fullRes .&&. prop_insertFull .&&. prop_insertFull

testReservoir :: IO ()
testReservoir = T.quickCheck props

{-----------------------------------------------------------------------------
    II. Approximate Median Naive
------------------------------------------------------------------------------}

-- * setting problem: can you make this streaming? assuming you know m?
-- * can you do running uniform distribution?

-- * Find approximate eps e-median of list `xs`
-- * Setting: list given before hand for some parameter t
-- * Use : runRVar (naive xs e d) StdRandom
medianNaive :: (Floating a, Ord a) => [a] -> Eps -> Delta -> RVar a
medianNaive xs e d | m <= t    = return $ median xs
             | otherwise = fmap median $ sampl'' xs m t
              where
                m = length xs
                t = round $ 7/(e^2) * log (2/d)

-- * use vector here due to (!!)
-- * sample `t` items from the list `[0..m-1]` *with replacement*
-- * Use : runRVar (sample m t) StdRandom
sampl'' :: [a] -> Int -> Int -> RVar [a]
sampl'' xs m t | m < 0 || t < 0   = return []
          | otherwise        = (\is -> [ xs !! i | i <- is]) <$> mis
            where mis = fmap sort . replicateM t $ uniform 0 (m-1)


text:
{-# LANGUAGE ScopedTypeVariables #-}

module Math.Probably.Unscented where

import Math.Probably.Sampler
import Math.Probably.FoldingStats
import Math.Probably.MCMC (empiricalMean, empiricalCovariance)
import Debug.Trace

import Numeric.LinearAlgebra hiding (find)

-- smells because it is inefficient. For testing.
smellyTransform :: Sampler a -> Int -> (a -> Vector Double) -> Sampler (Vector Double, Matrix Double)
smellyTransform dist n f = do
   ys  <- sequence $ replicate n $ fmap f dist
   let n = dim $ head ys
   return (empiricalMean ys, empiricalCovariance ys)

alphaUT = 1e-1
kappaUT = 0
betaUT = 2

unscentedTransform :: (Vector Double, Matrix Double) -- ^ random variable
                   -> (Vector Double -> Vector Double) -- ^ nonlinear transformation
                   -> (Vector Double, Matrix Double, [Vector Double])
unscentedTransform (xmean, xcov) f = (ymean, ycov, xs) where
   l = dim xmean
   lr = (realToFrac l) :: Double
   lambda = traceit "lambda " $alphaUT * alphaUT * (lr + kappaUT) - lr
   matSqrt = traceit "SQRT " $ toRows $ chol $ scale (lr+lambda) xcov
   xs = concat [[xmean],
                [xmean + (matSqrt!!i)| i <- [0..(l-1)]],
                [xmean - (matSqrt!!i)| i <- [0..(l-1)]]]
   wm = traceit "wm " $ (lambda / (lr+lambda)) : replicate (2*l) (1/(2*(lr+lambda)))
   wc = traceit "wc " $ ((lambda / (lr+lambda)) + 1 - alphaUT*alphaUT +betaUT) : replicate (2*l) (1/(2*(lr+lambda)))
   ys = map f xs
   ymean = sum $ zipWith (scale) wm ys
   ycov = sum [scale wi $ (yi - ymean) `outer` (yi - ymean) | (wi,yi) <- zip wc ys]

unscentedKalmanFilterAdditive :: (Vector Double -> Vector Double) -- ^ processes nonlinearity
                      -> (Vector Double -> Vector Double) -- ^ observation nonlinearity
                      -> Matrix Double -- ^ process noise
                      -> Matrix Double -- ^ observation noise
                      -> (Vector Double, Matrix Double) -- ^ initial state estimate
                      -> [Vector Double] -- ^ observations
                      -> [(Vector Double, Matrix Double)]
unscentedKalmanFilterAdditive procF obsF procCov obsCov (xmn0, xcov0) [] = []
unscentedKalmanFilterAdditive procF obsF procCov obsCov (xmn0, xcov0) (yobs:yobss) =
           (xmn1, xcov1) : unscentedKalmanFilterAdditive procF obsF procCov obsCov (xmn1, xcov1) (yobss) where
  nx = dim xmn0
  nrx = realToFrac $ nx
  kx = 3-nrx
  l = dim xmn0
  lr = (realToFrac l) :: Double
  lambda = alphaUT * alphaUT * (lr + kappaUT) - lr

  --xmn1zz = ymn1pred
  --xcov1zz = ycov1pred

  (xmn1pred, xcov1pred', _) = unscentedTransform (xmn0, xcov0) procF
  xcov1pred = xcov1pred' + procCov

  (ymn1pred, ycov1pred', chis) = unscentedTransform (xmn1pred, xcov1pred) obsF
  ycov1pred = ycov1pred' + obsCov

  ny = dim ymn1pred
  nry = realToFrac $ ny
  ky = 3-nry

  wc = ((lambda / (lr+lambda)) + 1 - alphaUT*alphaUT +betaUT) : replicate (2*l) (1/(2*(lr+lambda)))

  ws = (ky / (nry+ky)) : replicate (2*ny) (1/(2*(nry+ky)))
-- $ trace (show $ map dim [chii,xmn1pred, obsF chii, ymn1pred])
  crossCov = traceit "cross " $ sum [scale wi
                                      $ (chii - xmn1pred) `outer`
                                        (obsF chii - ymn1pred)
                      | (wi,chii) <- zip ws chis]

  gain = crossCov `multiply` inv ycov1pred

  xmn1 = xmn1pred + gain `mXv` (yobs - ymn1pred)
  xcov1 = traceit "newxcov " $ xcov1pred - gain `multiply` ycov1pred `multiply` trans gain

traceit s x = trace (s++" : "++show x) x

simDynamical :: Int
             -> (Vector Double -> Vector Double) -- ^ processes nonlinearity
             -> (Vector Double -> Vector Double) -- ^ observation nonlinearity
             -> Matrix Double -- ^ process noise
             -> Matrix Double -- ^ observation noise
             -> Vector Double -- ^ initial state
             -> Sampler [(Vector Double, Vector Double)]

simDynamical 0 procF obsF procCov obsCov x0 = return []
simDynamical n procF obsF procCov obsCov x0 = do
        xnoise <- multiNormal (constant 0$ dim x0) procCov
        let x1 = procF $ join [ x0, xnoise]
        ynoise <- multiNormal (constant 0$ rows obsCov) obsCov
        let y1 = obsF $ join [ x1, ynoise]
        rest <- simDynamical (n-1) procF obsF procCov obsCov x1
        return $ (x1, y1) : rest

unscentedKalmanFilter :: (Vector Double -> Vector Double) -- ^ processes nonlinearity
                      -> (Vector Double -> Vector Double) -- ^ observation nonlinearity
                      -> Matrix Double -- ^ process noise
                      -> Matrix Double -- ^ observation noise
                      -> (Vector Double, Matrix Double) -- ^ initial state estimate
                      -> [Vector Double] -- ^ observations
                      -> [(Vector Double, Matrix Double)]
unscentedKalmanFilter procF obsF procCov obsCov (xmn0, xcov0) [] = []
unscentedKalmanFilter procF obsF procCov obsCov (xmn0, xcov0) (yobs:yobss) =
           (xmn1, xcov1) : unscentedKalmanFilter procF obsF procCov obsCov (xmn1, xcov1) (yobss) where
  nx = dim xmn0
  nrx = realToFrac $ nx
  kx = 3-nrx

  xmn1zz = ymn1pred
  xcov1zz = ycov1pred

  xmn0aug = traceit "xm0aug" $ join [xmn0, constant 0 nx]
  xcov0aug = fromBlocks [[xcov0, 0],
                         [0,     procCov]]
  (xmn1pred, xcov1pred, _) = unscentedTransform (xmn0aug, xcov0aug) procF

  xmn1aug = traceit "xm1aug"$ join [xmn1pred, constant 0 nx]
  xcov1aug = fromBlocks [[xcov1pred, 0],
                         [0,         obsCov]]

  (ymn1pred, ycov1pred, chis) = unscentedTransform (xmn1aug, xcov1aug) obsF

  ny = dim ymn1pred
  nry = realToFrac $ ny
  ky = 3-nry

  ws = (ky / (nry+ky)) : replicate (2*ny) (1/(2*(nry+ky)))
-- $ trace (show $ map dim [chii,xmn1pred, obsF chii, ymn1pred])
  crossCov = traceit "cross " $ sum [scale wi
                                      $ (subVector 0 nx chii - xmn1pred) `outer`
                                        (obsF chii - ymn1pred)
                      | (wi,chii) <- zip ws chis]

  gain = crossCov `multiply` inv ycov1pred

  xmn1 = xmn1pred + gain `mXv` (yobs - ymn1pred)
  xcov1 = traceit "newxcov " $ xcov1pred - gain `multiply` ycov1pred `multiply` trans gain

max_stars_repo_name: glutamate/probably

max_issues_repo_name: glutamate/probably

max_forks_repo_name: glutamate/probably


text:
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE LambdaCase #-}

module LinearAlgebra (linearAlgebraTest) where

import Test.QuickCheck
import Test.Tasty.QuickCheck

import Test.Tasty

import System.Exit

import Numeric.LinearAlgebra (rank)
import qualified Numeric.LinearAlgebra.Data as Matrix

import Math.Tensor.Internal.LinearAlgebra (independentColumns, independentColumnsMat)

data SmallInt = S0 | S1 deriving (Show, Ord, Eq, Enum, Bounded)

toSmall :: Int -> SmallInt
toSmall 0 = S0
toSmall 1 = S1
toSmall i = error $ "cannot convert " ++ show i ++ " to SmallInt"

fromSmall :: Num a => SmallInt -> a
fromSmall S0 = 0
fromSmall S1 = 1

instance Arbitrary SmallInt where
    arbitrary = arbitraryBoundedEnum

data MatrixData a = MatrixData (Positive Int) (Positive Int) [a] deriving Show

instance Arbitrary a => Arbitrary (MatrixData a) where
    arbitrary = do
      m@(Positive m') <- arbitrary
      n@(Positive n') <- arbitrary
      xs <- vector (m'*n')
      return $ MatrixData m n xs

prop_smallValues :: MatrixData SmallInt -> Bool
prop_smallValues (MatrixData (Positive rows) (Positive cols) xs) =
    rank mat' == rank mat
  where
    mat  = (rows Matrix.>< cols) $ map fromSmall xs
    mat' = independentColumnsMat mat

prop_ints :: MatrixData Int -> Bool
prop_ints (MatrixData (Positive rows) (Positive cols) xs) =
    rank mat' == rank mat
  where
    mat  = (rows Matrix.>< cols) $ map fromIntegral xs
    mat' = independentColumnsMat mat

prop_doubles :: MatrixData Double -> Bool
prop_doubles (MatrixData (Positive rows) (Positive cols) xs) =
    rank mat' == rank mat
  where
    mat  = (rows Matrix.>< cols) xs
    mat' = independentColumnsMat mat

prop_consec :: Positive Int -> Int -> Bool
prop_consec (Positive dim') start =
    independentColumns mat == [0,1]
  where
    dim = dim' + 100
    mat = (dim Matrix.>< dim) $ map fromIntegral [start..]

testCase1 = testProperty "prop_smallValues" prop_smallValues
testCase2 = testProperty "prop_ints" prop_ints
testCase3 = testProperty "prop_doubles" prop_doubles
testCase4 = testProperty "prop_consec" prop_consec

linearAlgebraTest = testGroup "LinearAlgebraTest" [testCase1, testCase2, testCase3, testCase4]


text:
{-# LANGUAGE DeriveGeneric             #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE InstanceSigs              #-}
{-# LANGUAGE LambdaCase                #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE PartialTypeSignatures     #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE TypeApplications          #-}
{-# LANGUAGE TypeFamilies              #-}
{-# LANGUAGE TypeInType                #-}
{-# LANGUAGE UndecidableInstances      #-}

module Learn.Neural.Layer.Recurrent.LSTM (
  LSTM
  ) where

import           Data.Kind
import           Data.Singletons.Prelude
import           GHC.Generics                   (Generic)
import           GHC.Generics.Numeric
import           GHC.TypeLits
import           Learn.Neural.Layer
import           Numeric.BLAS
import           Numeric.Backprop
import           Numeric.Backprop.Iso
import           Statistics.Distribution
import           Statistics.Distribution.Normal
import qualified Generics.SOP                   as SOP

data LSTM :: Type

instance ( BLAS b
         , KnownNat i
         , KnownNat o
         , Floating (b '[o])
         , Floating (b '[o,i])
         , Floating (b '[o,o])
         )
        => Component LSTM b '[i] '[o] where

    data CParam LSTM b '[i] '[o] =
        LP { _lpForgetInp     :: !(b '[o,i])
           , _lpForgetState   :: !(b '[o,o])
           , _lpForgetBias    :: !(b '[o])
           , _lpRememberInp   :: !(b '[o,i])
           , _lpRememberState :: !(b '[o,o])
           , _lpRememberBias  :: !(b '[o])
           , _lpCommitInp     :: !(b '[o,i])
           , _lpCommitState   :: !(b '[o,o])
           , _lpCommitBias    :: !(b '[o])
           , _lpOutInp        :: !(b '[o,i])
           , _lpOutState      :: !(b '[o,o])
           , _lpOutBias       :: !(b '[o])
           }
      deriving (Generic)
    data CState LSTM b '[i] '[o] =
        LS { _lsCellState   :: !(b '[o])
           , _lsHiddenState :: !(b '[o])
           }
      deriving (Generic)
    type CConstr LSTM b '[i] '[o] = (Num (b '[o,i]), Num (b '[o,o]), Num (b '[o]))
    data CConf   LSTM b '[i] '[o] = forall d. ContGen d => LC d

    componentOp :: forall s. (Num (b '[i]), Num (b '[o]), CConstr LSTM b '[i] '[o])
        => OpB s '[b '[i], CParam LSTM b '[i] '[o], CState LSTM b '[i] '[o]] '[b '[o], CState LSTM b '[i] '[o]]
    componentOp = bpOp . withInps $ \(x :< p :< s :< Ø) -> do
        fI :< fS :< fB :< rI :< rS :< rB :< cI :< cS :< cB :< oI :< oS :< oB :< Ø <- gTuple #<~ p
        sC :< sH :< Ø <- gTuple #<~ s
        let forget   = sum [ matVecOp .$ (fI :< x  :< Ø)
                           , matVecOp .$ (fS :< sH :< Ø)
                           , fB
                           ]
            remember = sum [ matVecOp .$ (rI :< x  :< Ø)
                           , matVecOp .$ (rS :< sH :< Ø)
                           , rB
                           ]
            commit   = sum [ matVecOp .$ (cI :< x  :< Ø)
                           , matVecOp .$ (cS :< sH :< Ø)
                           , cB
                           ]
            out      = sum [ matVecOp .$ (oI :< x  :< Ø)
                           , matVecOp .$ (oS :< sH :< Ø)
                           , oB
                           ]
        forget'   <- tmapOp logistic ~$ (forget   :< Ø)
        remember' <- tmapOp logistic ~$ (remember :< Ø)
        commit'   <- tmapOp tanh     ~$ (commit   :< Ø)
        out'      <- tmapOp logistic ~$ (out      :< Ø)
        sC'       <- tmapOp tanh ~$ ((forget' * sC + remember' * commit') :< Ø)
        finalOut  <- bindVar $ out' * sC'
        s' :< Ø   <- isoVar (from (gTuple @(CState LSTM b '[i] '[o])) . tup1)
                            (sC' :< finalOut :< Ø)
        return $ finalOut :< s' :< Ø
      where
        logistic :: Floating a => a -> a
        logistic x = 1 / (1 + exp (-x))

    defConf = LC (normalDistr 0 0.01)
    initParam = \case
      i `SCons` SNil -> \case
        so@(o `SCons` SNil) -> \(LC d) g -> do
          LP <$> genA (o `SCons` (i `SCons` SNil)) (\_ -> realToFrac <$> genContVar d g)
             <*> genA (o `SCons` (o `SCons` SNil)) (\_ -> realToFrac <$> genContVar d g)
             <*> genA so                           (\_ -> realToFrac <$> genContVar d g)
             <*> genA (o `SCons` (i `SCons` SNil)) (\_ -> realToFrac <$> genContVar d g)
             <*> genA (o `SCons` (o `SCons` SNil)) (\_ -> realToFrac <$> genContVar d g)
             <*> genA so                           (\_ -> realToFrac <$> genContVar d g)
             <*> genA (o `SCons` (i `SCons` SNil)) (\_ -> realToFrac <$> genContVar d g)
             <*> genA (o `SCons` (o `SCons` SNil)) (\_ -> realToFrac <$> genContVar d g)
             <*> genA so                           (\_ -> realToFrac <$> genContVar d g)
             <*> genA (o `SCons` (i `SCons` SNil)) (\_ -> realToFrac <$> genContVar d g)
             <*> genA (o `SCons` (o `SCons` SNil)) (\_ -> realToFrac <$> genContVar d g)
             <*> genA so                           (\_ -> realToFrac <$> genContVar d g)
        _ -> error "inaccessible"
      _ -> error "inaccessible"
    initState _ so (LC d) g =
        LS <$> genA so (\_ -> realToFrac <$> genContVar d g)
           <*> genA so (\_ -> realToFrac <$> genContVar d g)

instance ( BLAS b
         , KnownNat i
         , KnownNat o
         , Floating (b '[o])
         , Floating (b '[o,i])
         , Floating (b '[o,o])
         )
      => ComponentLayer 'Recurrent LSTM b '[i] '[o] where
    componentRunMode = RMNotFF

instance SOP.Generic (CState LSTM b '[i] '[o])
instance SOP.Generic (CParam LSTM b '[i] '[o])

instance (Num (b '[o,i]), Num (b '[o,o]), Num (b '[o])) => Num (CParam LSTM b '[i] '[o]) where
    (+)         = genericPlus
    (-)         = genericMinus
    (*)         = genericTimes
    negate      = genericNegate
    abs         = genericAbs
    signum      = genericSignum
    fromInteger = genericFromInteger

instance (Fractional (b '[o,i]), Fractional (b '[o,o]), Fractional (b '[o])) => Fractional (CParam LSTM b '[i] '[o]) where
    (/)          = genericDivide
    recip        = genericRecip
    fromRational = genericFromRational

instance (Floating (b '[o,i]), Floating (b '[o,o]), Floating (b '[o])) => Floating (CParam LSTM b '[i] '[o]) where
    pi      = genericPi
    exp     = genericExp
    (**)    = genericPower
    log     = genericLog
    logBase = genericLogBase
    sin     = genericSin
    cos     = genericCos
    tan     = genericTan
    asin    = genericAsin
    acos    = genericAcos
    atan    = genericAtan
    sinh    = genericSinh
    cosh    = genericCosh
    tanh    = genericTanh
    asinh   = genericAsinh
    acosh   = genericAcosh
    atanh   = genericAtanh

instance (Num (b '[o,i]), Num (b '[o,o]), Num (b '[o])) => Num (CState LSTM b '[i] '[o]) where
    (+)         = genericPlus
    (-)         = genericMinus
    (*)         = genericTimes
    negate      = genericNegate
    abs         = genericAbs
    signum      = genericSignum
    fromInteger = genericFromInteger

instance (Fractional (b '[o,i]), Fractional (b '[o,o]), Fractional (b '[o])) => Fractional (CState LSTM b '[i] '[o]) where
    (/)          = genericDivide
    recip        = genericRecip
    fromRational = genericFromRational

instance (Floating (b '[o,i]), Floating (b '[o,o]), Floating (b '[o])) => Floating (CState LSTM b '[i] '[o]) where
    pi      = genericPi
    exp     = genericExp
    (**)    = genericPower
    log     = genericLog
    logBase = genericLogBase
    sin     = genericSin
    cos     = genericCos
    tan     = genericTan
    asin    = genericAsin
    acos    = genericAcos
    atan    = genericAtan
    sinh    = genericSinh
    cosh    = genericCosh
    tanh    = genericTanh
    asinh   = genericAsinh
    acosh   = genericAcosh
    atanh   = genericAtanh