text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.StudentT
-- Copyright : (c) 2011 Aleksey Khudyakov
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Student-T distribution
module Statistics.Distribution.StudentT (
    StudentT
    -- * Constructors
  , studentT
  , studentTE
  , studentTUnstandardized
    -- * Accessors
  , studentTndf
  ) where

import Control.Applicative
import Data.Aeson          (FromJSON(..), ToJSON, Value(..), (.:))
import Data.Binary         (Binary(..))
import Data.Data           (Data, Typeable)
import GHC.Generics        (Generic)
import Numeric.SpecFunctions (
  logBeta, incompleteBeta, invIncompleteBeta, digamma)

import qualified Statistics.Distribution as D
import Statistics.Distribution.Transform (LinearTransform (..))
import Statistics.Internal

-- | Student-T distribution
newtype StudentT = StudentT { studentTndf :: Double }
                   deriving (Eq, Typeable, Data, Generic)

instance Show StudentT where
  showsPrec i (StudentT ndf) = defaultShow1 "studentT" ndf i
instance Read StudentT where
  readPrec = defaultReadPrecM1 "studentT" studentTE

instance ToJSON StudentT
instance FromJSON StudentT where
  parseJSON (Object v) = do
    ndf <- v .: "studentTndf"
    maybe (fail $ errMsg ndf) return $ studentTE ndf
  parseJSON _ = empty

instance Binary StudentT where
  put = put . studentTndf
  get = do
    ndf <- get
    maybe (fail $ errMsg ndf) return $ studentTE ndf

-- | Create Student-T distribution. Number of parameters must be positive.
studentT :: Double -> StudentT
studentT ndf = maybe (error $ errMsg ndf) id $ studentTE ndf

-- | Create Student-T distribution. Number of parameters must be positive.
studentTE :: Double -> Maybe StudentT
studentTE ndf
  | ndf > 0   = Just (StudentT ndf)
  | otherwise = Nothing

errMsg :: Double -> String
errMsg _ = modErr "studentT" "non-positive number of degrees of freedom"

instance D.Distribution StudentT where
  cumulative      = cumulative
  complCumulative = complCumulative

instance D.ContDistr StudentT where
  density    d@(StudentT ndf) x = exp (logDensityUnscaled d x) / sqrt ndf
  logDensity d@(StudentT ndf) x = logDensityUnscaled d x - log (sqrt ndf)
  quantile = quantile

cumulative :: StudentT -> Double -> Double
cumulative (StudentT ndf) x
  | x > 0     = 1 - 0.5 * ibeta
  | otherwise = 0.5 * ibeta
  where
    ibeta = incompleteBeta (0.5 * ndf) 0.5 (ndf / (ndf + x*x))

complCumulative :: StudentT -> Double -> Double
complCumulative (StudentT ndf) x
  | x > 0     = 0.5 * ibeta
  | otherwise = 1 - 0.5 * ibeta
  where
    ibeta = incompleteBeta (0.5 * ndf) 0.5 (ndf / (ndf + x*x))

logDensityUnscaled :: StudentT -> Double -> Double
logDensityUnscaled (StudentT ndf) x =
    log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf)

quantile :: StudentT -> Double -> Double
quantile (StudentT ndf) p
  | p >= 0 && p <= 1 =
    let x = invIncompleteBeta (0.5 * ndf) 0.5 (2 * min p (1 - p))
    in case sqrt $ ndf * (1 - x) / x of
         r | p < 0.5   -> -r
           | otherwise -> r
  | otherwise = modErr "quantile" $ "p must be in [0,1] range. Got: "++show p

instance D.MaybeMean StudentT where
  maybeMean (StudentT ndf) | ndf > 1   = Just 0
                           | otherwise = Nothing

instance D.MaybeVariance StudentT where
  maybeVariance (StudentT ndf) | ndf > 2   = Just $! ndf / (ndf - 2)
                               | otherwise = Nothing

instance D.Entropy StudentT where
  entropy (StudentT ndf) =
    0.5 * (ndf+1) * (digamma ((1+ndf)/2) - digamma(ndf/2))
    + log (sqrt ndf)
    + logBeta (ndf/2) 0.5

instance D.MaybeEntropy StudentT where
  maybeEntropy = Just . D.entropy

instance D.ContGen StudentT where
  genContVar = D.genContinuous

-- | Create an unstandardized Student-t distribution.
studentTUnstandardized :: Double -- ^ Number of degrees of freedom
                       -> Double -- ^ Central value (0 for standard Student T distribution)
                       -> Double -- ^ Scale parameter
                       -> LinearTransform StudentT
studentTUnstandardized ndf mu sigma
  | sigma > 0 = LinearTransform mu sigma $ studentT ndf
  | otherwise = modErr "studentTUnstandardized" $ "sigma must be > 0. Got: " ++ show sigma

modErr :: String -> String -> a
modErr fun msg = error $ "Statistics.Distribution.StudentT." ++ fun ++ ": " ++ msg


text:
-- {-# LANGUAGE DuplicateRecordFields #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
-- {-# LANGUAGE TypeSynonymInstances #-} -- needed fo GenIO
-- {-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE AllowAmbiguousTypes #-} -- MWC.Gen create
{-# LANGUAGE UndecidableInstances #-} -- for UniformDist
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE PolyKinds #-}

{-# OPTIONS_HADDOCK show-extensions #-}

-- |
-- Module      : PSO.Random
-- Description : Uniform interface for "Control.Random.Mersenne" and
--               "Control.Random.MWC"
-- Description : Particle Swarm Optimisation
-- Copyright   : (c) Tom Westerhout, 2017
-- License     : BSD3
-- Maintainer  : t.westerhout@student.ru.nl
-- Stability   : experimental

module PSO.Random
    ( Generator(..)
    , Randomisable(..)
    , UniformDist(..)
    , RandomScalable(..)
    , mkMTGen
    , mkMWCGen
    , mkMWCGenST
    , uniformList
    , uniformVector
    , randomSpin
    ) where

import System.IO.Unsafe
import Data.Word(Word8,Word32)
import Data.Proxy
import Data.Complex(Complex(..))
import Control.Monad.Primitive
import Control.Monad.ST
import Control.Monad.Reader
import Data.Vector(singleton)
import GHC.Float
import qualified Data.Vector.Storable as V
import qualified Data.Vector.Generic as GV
import Foreign.C.Types(CDouble, CFloat)
import qualified Foreign
import qualified System.Random.Mersenne as Mersenne
import qualified System.Random.MWC as MWC

class (PrimMonad m) => Generator m g where
    create :: Maybe Word32 -> m g

instance Generator IO Mersenne.MTGen where
    create = Mersenne.newMTGen

instance Generator IO (MWC.Gen RealWorld) where
    create x = case x of
        (Just n) -> MWC.initialize . singleton $ n
        Nothing  -> MWC.createSystemRandom

instance Generator (ST s) (MWC.Gen s) where
    create x = case x of
        (Just n) -> MWC.initialize . singleton $ n
        Nothing  -> MWC.create

class (Monad m) => Randomisable m a where
    random :: m a

instance (Mersenne.MTRandom a)
  => Randomisable (ReaderT Mersenne.MTGen IO) a where
    random = ask >>= lift . Mersenne.random

instance {-# OVERLAPS #-} Randomisable (ReaderT Mersenne.MTGen IO) Float where
  random = ask >>= lift . liftM double2Float . Mersenne.random

instance (MWC.Variate a)
  => Randomisable (ReaderT (MWC.Gen RealWorld) IO) a where
    random = ask >>= lift . MWC.uniform

instance (MWC.Variate a)
  => Randomisable (ReaderT (MWC.Gen s) (ST s)) a where
    random = ask >>= lift . MWC.uniform

instance {-# OVERLAPS #-} (MWC.Variate a)
  => Randomisable (ReaderT (MWC.Gen RealWorld) IO) (Complex a) where
    random = liftM2 (:+) random random

mkMTGen :: Maybe Word32 -> IO Mersenne.MTGen
mkMTGen = create

mkMWCGen :: Maybe Word32 -> IO MWC.GenIO
mkMWCGen = create

mkMWCGenST :: Maybe Word32 -> ST s (MWC.GenST s)
mkMWCGenST = create

class (Monad m) => UniformDist m a where
  uniform :: (a, a) -> m a

class (Monad m) => RandomScalable m a where
  randScale :: a -> m a

uniformFloating :: (RealFloat a, Randomisable m a)
  => (a, a) -> m a
uniformFloating (low, high) = return (\x -> low + (high - low) * x) `ap` random

uniformIntegralMWC ::
     (Integral a, MWC.Variate a, PrimMonad m)
  => (a, a) -> ReaderT (MWC.Gen (PrimState m)) m a
uniformIntegralMWC bounds = ask >>= lift . MWC.uniformR bounds

uniformIntegralMT ::
     (Integral a, Bounded a, Mersenne.MTRandom a)
  => (a, a) -> ReaderT Mersenne.MTGen IO a
uniformIntegralMT (low, high) = ask >>= lift . Mersenne.random >>= \x ->
  return $ low + x `mod` (high - low + 1)

instance (Monad m, Randomisable m Float) => UniformDist m Float where
    uniform = uniformFloating

instance (Monad m, Randomisable m Double) => UniformDist m Double where
    uniform = uniformFloating

instance (Monad m, Randomisable m CFloat) => UniformDist m CFloat where
    uniform = uniformFloating

instance (Monad m, Randomisable m CDouble) => UniformDist m CDouble where
    uniform = uniformFloating

instance (RealFloat a, UniformDist m a)
  => UniformDist m (Complex a) where
    uniform ((rlow :+ ilow), (rhigh :+ ihigh)) =
      liftM2 (:+) (uniform (rlow, rhigh)) (uniform (ilow, ihigh))

instance (UniformDist m a) => UniformDist m [a] where
  uniform (low, high) = mapM uniform $ zip low high

instance (Foreign.Storable a, UniformDist m a)
  => UniformDist m (V.Vector a) where
    uniform (low, high) = V.zipWithM (\l h -> uniform (l, h)) low high

instance UniformDist (ReaderT (MWC.Gen RealWorld) IO) Word8 where
    uniform = uniformIntegralMWC

instance UniformDist (ReaderT (MWC.Gen RealWorld) IO) Int where
    uniform = uniformIntegralMWC

instance UniformDist (ReaderT Mersenne.MTGen IO) Word8 where
    uniform = uniformIntegralMT

instance UniformDist (ReaderT Mersenne.MTGen IO) Int where
    uniform = uniformIntegralMT

randScaleFloating :: (RealFloat a, Randomisable m a) => a -> m a
randScaleFloating x = liftM (* x) random

instance (Monad m, Randomisable m Float) => RandomScalable m Float where
    randScale = randScaleFloating

instance (Monad m, Randomisable m Double) => RandomScalable m Double where
    randScale = randScaleFloating

instance (Monad m, Randomisable m CFloat) => RandomScalable m CFloat where
    randScale = randScaleFloating

instance (Monad m, Randomisable m CDouble) => RandomScalable m CDouble where
    randScale = randScaleFloating

instance {-# OVERLAPS #-} (RandomScalable m a)
  => RandomScalable m (Complex a) where
    randScale (x :+ y) = return (:+) `ap` (randScale x) `ap` (randScale y)

instance (Monad m, Foreign.Storable a, RandomScalable m a)
  => RandomScalable m (V.Vector a) where
    randScale = V.mapM randScale

-- instance (RandomScalable m a, GV.Vector v a) => RandomScalable m (v a) where
--     randScale = GV.mapM randScale

uniformList :: (Monad m, UniformDist m a) => Int -> (a, a) -> m [a]
uniformList n bounds = replicateM n (uniform bounds)

uniformVector ::
     (Monad m, Foreign.Storable a, UniformDist m a)
  => Int -> (a, a) -> m (V.Vector a)
uniformVector n = liftM V.fromList . uniformList n

randomSpin ::
     (Monad m, Foreign.Storable a, Num a, Randomisable m Bool)
  => Int -> m (V.Vector a)
randomSpin n =
  let fromBool True  = 1
      fromBool False = (-1)
   in V.replicateM n (fromBool <$> random)

ext: hs

max_stars_repo_name: twesterhout/tcm-swarm

max_issues_repo_name: twesterhout/tcm-swarm

max_forks_repo_name: twesterhout/tcm-swarm


text:
-- | A module that shows some simple examples demonstrating how to use the package.
module System.Console.Ansigraph.Examples (

    demo
  , legend
  , showColors

-- * Horizontal graphs
  , waveDemoComplex
  , waveDemoReal
  , waveDemoPositive
  , wave

-- * Matrix graphs

  , matDemoReal
  , matDemoComplex
  , unitary

) where

import System.Console.Ansigraph

import System.Console.ANSI
import Control.Monad      (forM_)
import Data.Complex       (Complex (..), cis, realPart)
-- for GHC <= 7.8
import Control.Applicative

---- Wave Demo ----

-- | A complex wave vector, part of the graph of /z(x,t) = exp(ix - it)/ in some units.
wave :: [Complex Double]
wave = cis . (*) (pi/20) <$> [0..79]

deltas :: [Double]
deltas = (*) (-pi/10) <$> [0..80]

waves :: [[Complex Double]]
waves = zipWith (\z -> map (* z)) (cis <$> deltas) $ repeat wave

rwaves :: [[Double]]
rwaves = map (map realPart) waves

pwaves :: [PosGraph]
pwaves = PosGraph . map (+1) <$> rwaves

-- | Display an animation of the positive real function /p(x,t) = cos(x-t) + 1/ in some units.
waveDemoPositive :: IO ()
waveDemoPositive = animate pwaves

-- | Display an animation of the real function /r(x,t) = cos(x-t)/ in the standard style, i.e. with both
--   positive and negative regions.
waveDemoReal :: IO ()
waveDemoReal = animate rwaves

-- | Display an animation of the complex wave /z(x,t) = exp(ix - it)/ in some units.
waveDemoComplex :: IO ()
waveDemoComplex = animate waves

---- Matrix Demos ----

vscale :: Num a => a -> [a] -> [a]
vscale x = map (* x)

mscale :: Num a => a -> [[a]] -> [[a]]
mscale x = map $ map (* x)

-- The following is a quick implementation of the matrix tensor product ('mox').
-- The details are not relevant to the library, but are only used for this example.

fromRealVs :: [Double] -> [Double] -> [Complex Double]
fromRealVs = zipWith (:+)

fromRealMs :: [[Double]] -> [[Double]] -> [[Complex Double]]
fromRealMs = zipWith fromRealVs

vox :: Num a => [a] -> [a] -> [a]
vox v w = concatMap (`vscale` w) v

-- matrix tensor product
mox :: Num a => [[a]] -> [[a]] -> [[a]]
mox m1 m2 = zipWith vox (stepOne m1 m2) (stepTwo m1 m2)
  where stepOne, stepTwo :: [[a]] -> [[a]] -> [[a]]
        stepOne x y = concatMap (replicate (length y)) x
        stepTwo x y = concat $ replicate (length x) y

---- Complex matrix example ----

sx, sz, sI :: [[Double]]
sz  = [[1,0],[0,-1]]
sx  = [[0,1],[1,0]]
sI  = [[1,0],[0,1]]

-- Time-exponentials of pauli matrices
-- exp(itσ) = cos(t)I + i sin(t)σ
sinSX, sinSZ :: Double -> [[Complex Double]]
sinSX t = fromRealMs (mscale (cos t) sI) (mscale (sin t) sx)
sinSZ t = fromRealMs (mscale (cos t) sI) (mscale (sin t) sz)

-- | An example of a time-dependent complex matrix.
unitary :: Double -> [[Complex Double]]
unitary t = sinSZ t `mox` sinSX (2*t)

slowDeltas :: [Double]
slowDeltas = (*) (pi/50) <$> [0..100]

-- | Shows an animation of an example time-dependent matrix formed from Pauli matrices, called
--   'unitary'. Specifically, it is the tensor product of σz and σx exponentiated with different
--   frequencies.
matDemoComplex :: IO ()
matDemoComplex = animate $ unitary <$> slowDeltas

---- Real matrix example ----

ry :: Double -> [[Double]]
ry t = [[cos t, 0, sin t]
       ,[0, 1, 0]
       ,[-(sin t),0,cos t]]

-- | An example real matrix animation.
matDemoReal :: IO ()
matDemoReal = animate $ (\t -> ry t `mox` ry (2*t)) <$> slowDeltas

---- Show ANSI Colors ----

colors = [Black,Red,Green,Yellow,Blue,Magenta,Cyan,White]
intensities = [Dull,Vivid]

ansicolors :: [AnsiColor]
ansicolors = [ AnsiColor i c | c <- colors, i <- intensities ]

-- | Show all of the available 'AnsiColor's with corresponding 'ColorIntensity', 'Color' pairs.
showColors = do
  boldStrLn noColoring "Available colors"
  newline
  forM_ ansicolors $ \c -> do
    let clr = Coloring Nothing (Just c)
    colorStr clr $ replicate 20 ' '
    putStrLn $ "  " ++ show (intensity c) ++ " " ++ show (color c)
  setSGR [Reset]

cb, bc :: Coloring
cb = mkColoring (AnsiColor Dull Black) (AnsiColor Vivid Cyan)
bc = invert cb

newline = putStrLn ""

verticalPad io = do
  newline
  io
  newline
  newline

-- | Displays a legend showing color conventions for supported graph types.
legend = do
  boldStrLn cb "       Legend       "
  newline

  boldStrLn noColoring "Horizontal Graphs"
  newline
  colorStr (fromBG blue) "  "
  putStrLn " Real component (positive and negative)"
  newline
  colorStr (fromBG pink) "  "
  putStrLn " Imag component (positive and negative)"
  newline

  boldStrLn noColoring "Matrix Graphs"
  newline
  putStr "  "
  colorStr (mkColoring white pink) "+i"
  putStrLn "  "
  colorStr (mkColoring white red) "-r"
  putStr "  "
  colorStrLn (mkColoring white blue) "+r"
  putStr "  "
  colorStr (mkColoring white green) "-i"
  putStrLn "  "

-- | Run all of the demos in sequence.
demo = do

  verticalPad $ boldStrLn cb "       Ansigraph demo       "

  putStr "Positive function graph  "
  colorStrLn bc " cos (x - t) + 1 "
  verticalPad waveDemoPositive

  putStr "Real function graph  "
  colorStrLn bc "  cos (x - t)  "
  verticalPad waveDemoReal

  putStr "Complex function graph  "
  colorStrLn bc "  exp (ix - it)  "
  verticalPad waveDemoComplex

  putStr "Real matrix graph  "
  colorStrLn bc "  rotate_Y(t) ⊗ rotate_Y(2t)  "
  verticalPad matDemoReal

  putStr "Complex matrix graph  "
  colorStrLn bc "  exp (it σz) ⊗ exp (2it σx)  "
  verticalPad matDemoComplex

  verticalPad showColors

  verticalPad legend

max_stars_repo_name: fieldstrength/ansigraph

max_issues_repo_name: fieldstrength/ansigraph

max_forks_repo_name: fieldstrength/ansigraph


text:
import Data.Complex

direction :: Char -> Complex Float
direction 'U' = 0 :+ 1
direction 'D' = 0 :+ (-1)
direction 'L' = (-1) :+ 0
direction 'R' = 1 :+ 0

clamp :: Float -> Float
clamp = max (-1) . min 1

advance :: Complex Float -> Complex Float -> Complex Float
advance pos dir =
  let (x :+ y) = pos + dir
  in clamp x :+ clamp y

complexToButton :: Complex Float -> Char
complexToButton ((-1) :+   1)  = '1'
complexToButton (  0  :+   1)  = '2'
complexToButton (  1  :+   1)  = '3'
complexToButton ((-1) :+   0)  = '4'
complexToButton (  0  :+   0)  = '5'
complexToButton (  1  :+   0)  = '6'
complexToButton ((-1) :+ (-1)) = '7'
complexToButton (  0  :+ (-1)) = '8'
complexToButton (  1  :+ (-1)) = '9'

part1 :: String -> String
part1 input = do
  line <- lines input
  let button = complexToButton $ foldl advance 0 $ map direction line
  return button

complexToButton2 :: Complex Float -> Char
complexToButton2 (  0  :+   2)  = '1'
complexToButton2 ((-1) :+   1)  = '2'
complexToButton2 (  0  :+   1)  = '3'
complexToButton2 (  1  :+   1)  = '4'
complexToButton2 ((-2) :+   0)  = '5'
complexToButton2 ((-1) :+   0)  = '6'
complexToButton2 (  0  :+   0)  = '7'
complexToButton2 (  1  :+   0)  = '8'
complexToButton2 (  2  :+   0)  = '9'
complexToButton2 ((-1) :+ (-1)) = 'A'
complexToButton2 (  0  :+ (-1)) = 'B'
complexToButton2 (  1  :+ (-1)) = 'C'
complexToButton2 (  0  :+ (-2)) = 'D'
complexToButton2 _              = '0' -- default

advance2 :: Complex Float -> Complex Float -> Complex Float
advance2 pos dir =
  let pos' = pos + dir
  in if complexToButton2 pos' == '0'
  then pos
  else pos'

part2 :: String -> String
part2 input = do
  line <- lines input
  let button = complexToButton2 $ foldl advance2 0 $ map direction line
  return button


text:
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE EmptyDataDeriving #-}
{-# LANGUAGE EmptyCase #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE ExplicitNamespaces #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE RankNTypes #-}
module Models.Integrals.Types (module Models.Integrals.Types, module TLC.Terms) where

-- import Data.Ratio
import Algebra.Classes
import qualified Algebra.Morphism.Affine as A
import qualified Algebra.Morphism.LinComb as LC
import Prelude hiding (Num(..), Fractional(..), (^), product, sum, pi, sqrt
                      , exp, (**))
import Data.Complex
import TLC.Terms (type (∈)(..), Type(..), type(×), type(⟶))
import qualified Algebra.Expression as E
import Data.Function (on)
import Data.Foldable
import Data.List
--------------------------------------------------------------------------------
-- | Types

type C = Complex Double

data Zero deriving (Show,Eq,Ord)
newtype Number = Number (E.Expr Zero) deriving
  (Additive,Group,AbelianAdditive,Multiplicative,Division,Field,Ring,Roots,Transcendental)
instance Scalable Number Number where (*^) = (*)
type Rat = Number
fromNumber :: Number -> E.Expr Zero
fromNumber (Number x) = x
evalNumber :: Number -> Double
evalNumber (Number e) = E.eval (\case) e
instance Show Number where show = show . fromNumber
instance Eq Rat where
  (==) = (==) `on` evalNumber
instance Ord Rat where
  compare = compare `on` evalNumber
instance DecidableZero Rat where
  isZero  = (== 0) . evalNumber
type RatLike α = (Ring α, Ord α, DecidableZero α)

type Var γ = 'R ∈ γ
data Dir = Min | Max deriving (Eq,Ord,Show)
data Elem γ = Vari (Var γ)
            | Supremum Dir [Ret γ] -- either minimum or maximum of the arguments
  deriving (Eq, Ord, Show)
type Available α γ = α ∈ γ
type Expr γ = A.Affine (Var γ) Rat
type Ret γ = E.Expr (Elem γ)
type Cond γ = Cond' (Expr γ)
data Cond' e = IsNegative { condExpr :: e }
              -- Meaning of this constructor: expression ≤ 0
              | IsZero { condExpr :: e }
              -- Meaning of this constructor: expression = 0
   deriving (Eq, Show, Functor, Foldable, Traversable, Ord)

data Domain γ = Domain { domainLoBounds, domainHiBounds :: [Expr γ] }
  deriving (Show, Eq, Ord)

-- | needs to be 1st a order representation for optimisations to be
-- implementable
data P (γ :: Type) where
  Done :: Ret γ -> P γ
  Cond :: Cond γ -> P γ -> P γ
  Integrate :: Domain γ -> P (γ × 'R) -> P γ
  Add :: P γ -> P γ -> P γ
  Power :: P γ -> Rat -> P γ
  Mul :: [P γ] -> P γ
  -- Can this replaced by "Scale"? No, because we do integration in
  -- normalisation factors as well.
  -- Scale :: Ret γ -> P γ -> P γ
  deriving (Ord, Eq)

lift' :: Applicative f
  => (forall v. Available v γ -> f (Available v δ))
  -> (forall v. Available v (γ × α) -> f (Available v (δ × α)))
lift' _ (Get) = pure Get
lift' f (Weaken x) = Weaken <$> (f x)

class VarTraversable t where
  varTraverse :: (Applicative f)
    => (forall x. Available x γ -> f (Available x δ))
    -> t γ -> f (t δ)

instance VarTraversable Domain where
 varTraverse f (Domain los his)
  = Domain <$> traverse (A.traverseVars f) los <*>
               traverse (A.traverseVars f) his

instance VarTraversable P where
  varTraverse f = \case
    Power e k -> Power <$> varTraverse f e <*> pure k
    Done x -> Done <$> traverse (varTraverse f) x
    Mul xs -> Mul <$> traverse (varTraverse f) xs
    Integrate d e ->
      Integrate <$> (varTraverse f d) <*> (varTraverse (lift' f) e)
    Cond e x -> Cond <$> traverse (A.traverseVars f) e <*> varTraverse f x
    Add x y  -> Add <$> varTraverse f x <*> varTraverse f y

instance VarTraversable Elem where
  varTraverse f = \case
    Vari x -> Vari <$> f x
    Supremum d es -> Supremum d <$> traverse (traverse (varTraverse f)) es

deriving instance Show (P γ)

----------------------------------
-- | Smart constructors

conds_ :: [Cond γ] -> P γ -> P γ
conds_ cs e = foldr Cond e cs

isPositive,isNegative :: Expr γ -> Cond γ
isPositive e = isNegative (negate e)
isNegative e = IsNegative e

-- Domain without restriction
full :: Domain γ
full = Domain [] []

lessThan, greaterThan :: Expr γ -> Expr γ -> Cond γ
t `lessThan` u = isNegative (t - u)
t `greaterThan` u = u `lessThan` t

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

instance Multiplicative (P γ) where
  one = Done one
  x * y = Mul [x,y]

instance Division (P γ) where
  recip x = Power x (negate one)

instance AbelianAdditive (P γ)
instance Group (P γ) where
  negate = (negate (one :: Ret γ) *^)
instance Scalable (E.Expr (Elem γ)) (P γ) where
  k *^ e = Done k * e

pattern PZero :: P γ
pattern PZero <- Done E.Zero

instance Additive (P γ) where
  zero =  Done E.Zero
  PZero + x = x
  x + PZero = x
  x + y = Add x y

setHere :: f γ -> (Var γ -> f γ) -> Var (γ × α) -> f γ
setHere a f = \case
  Get -> a
  Weaken x -> f x

----------------------------------------------------------------
-- Normalising substitutions of variables to affine expressions

type SubstE γ δ = Var γ -> Expr δ

wkSubst :: SubstE γ δ -> SubstE (γ × α) (δ × α)
wkSubst f = \case
  Get -> A.var Get
  Weaken x -> A.mapVars Weaken (f x)

substExpr :: SubstE γ δ ->  Expr γ -> Expr δ
substExpr = A.subst

substCond :: SubstE γ δ -> Cond γ -> Cond δ
substCond f = fmap (substExpr f)

substDomain :: SubstE γ δ -> Domain γ -> Domain δ
substDomain f (Domain lo hi) = Domain
                                 (nub (substExpr f <$> lo))
                                 (nub (substExpr f <$> hi))

wkP :: P γ -> P (γ × α)
wkP = substP $ \i -> A.var (Weaken i)

substElem :: forall γ ζ. SubstE γ ζ -> Elem γ -> Ret ζ
substElem v = \case
  Supremum dir es -> supremum dir (substRet v <$> es)
  Vari x -> exprToPoly (v x)

substRet  :: forall γ ζ. SubstE γ ζ -> Ret γ -> Ret ζ
substRet v = E.eval (substElem v)

substP :: SubstE γ δ -> P γ -> P δ
substP f p0 = case p0 of
  Done e -> Done (E.eval (substElem f) e)
  Add p1 p2 -> substP f p1 + substP f p2
  Power p k -> Power (substP f p) k
  Mul ps -> Mul (substP f <$> ps)
  Cond c p -> Cond (substCond f c) (substP f p)
  Integrate d p -> Integrate (substDomain f d) (substP (wkSubst f) p) -- integrations are never simplified by substitution

swap2P :: P (γ × α × β) -> P (γ × β × α)
swap2P = substP $ \case
  Get -> A.var (Weaken Get)
  Weaken Get -> A.var Get
  Weaken (Weaken x) -> A.var (Weaken (Weaken x))

wkExpr :: Expr γ -> Expr (γ × β)
wkExpr = substExpr (A.var . Weaken)

condVars :: Cond γ -> [Var γ]
condVars c = case condExpr c of
   (A.Affine _ e) -> map fst (LC.toList e)

retVars :: Ret γ -> [Var γ]
retVars x = concatMap elemVars (toList x)

elemVars :: Elem γ -> [Var γ]
elemVars = \case
   Vari x -> [x]
   Supremum _ es -> concatMap retVars es

supremum :: Dir -> [Ret γ] -> Ret γ
supremum _ [e] = e
supremum dir es =
  case traverse (fmap retToNumber . traverse (varTraverse (const Nothing))) es of
    Just cs | not (null es) ->
       constPoly ((case dir of
                     Max -> maximum
                     Min -> minimum)
                   cs)
    _ -> E.Var (Supremum dir es)

constPoly :: Number -> Ret γ
constPoly (Number n) = (\case) <$> n

varPoly :: 'R ∈ γ -> Ret γ
varPoly = E.Var . Vari

retToNumber :: Ret 'Unit -> Number
retToNumber = E.eval $ \case
     Vari x -> case x of
     Supremum d xs -> (case d of Max -> maximum; Min -> minimum) (fmap retToNumber xs)

numberToRet :: Number -> Ret γ
numberToRet (Number x) = fmap (\case) x
exprToPoly :: Expr γ -> Ret γ
exprToPoly = A.eval (fmap (\case) .fromNumber) (E.Var .  Vari)

mkSuprema :: Domain γ -> (Ret γ, Ret γ)
mkSuprema (Domain los his) = (supremum Max $ map exprToPoly los,
                              supremum Min $ map exprToPoly his)


text:
module Levity.Class.Functor
  ( FunctorL(..)
  ) where

import Control.Applicative (ZipList)
import Control.Exception (Handler)
import Control.Monad.ST (ST)
import Data.Complex (Complex)
import Data.Functor.Const (Const)
import Data.Functor.Identity (Identity)
import Data.List.NonEmpty (NonEmpty)
import Data.Monoid qualified as Monoid
import Data.Proxy (Proxy)
import Data.Semigroup qualified as Semigroup
import GHC.Base (Type)
import GHC.Conc (STM)
import GHC.Exts
import GHC.Generics qualified as Generic
import Prelude
import System.Console.GetOpt (ArgDescr, OptDescr, ArgOrder)
import Text.ParserCombinators.ReadP (ReadP)
import Text.ParserCombinators.ReadPrec (ReadPrec)

import Levity.Types.Base

type FunctorL :: (Type -> TYPE r) -> Constraint
class FunctorL f where
  {-# minimal fmapL #-}
  fmapL :: (a -> b) -> (f a -> f b)

instance Prelude.Functor f => FunctorL (Base1 f) where
  fmapL = Prelude.fmap

-- Manually-defined instances are because of: https://gitlab.haskell.org/ghc/ghc/-/issues/20680
deriving via (Base1 []) instance FunctorL []
deriving via (Base1 Prelude.Maybe) instance FunctorL Prelude.Maybe
deriving via (Base1 Prelude.IO) instance FunctorL Prelude.IO
deriving via (Base1 Generic.Par1) instance FunctorL Generic.Par1
deriving via (Base1 NonEmpty) instance FunctorL NonEmpty
deriving via (Base1 ReadP) instance FunctorL ReadP
deriving via (Base1 ReadPrec) instance FunctorL ReadPrec
deriving via (Base1 Down) instance FunctorL Down
deriving via (Base1 (Semigroup.Arg a)) instance FunctorL (Semigroup.Arg a)
deriving via (Base1 (ST s)) instance FunctorL (ST s)
deriving via (Base1 ((,) a)) instance FunctorL ((,) a)
deriving via (Base1 (Prelude.Either a)) instance FunctorL (Prelude.Either a)
deriving via (Base1 Complex) instance FunctorL Complex
deriving via (Base1 Semigroup.Min) instance FunctorL Semigroup.Min
deriving via (Base1 Semigroup.Max) instance FunctorL Semigroup.Max
deriving via (Base1 Semigroup.First) instance FunctorL Semigroup.First
deriving via (Base1 Semigroup.Last) instance FunctorL Semigroup.Last
deriving via (Base1 ArgOrder) instance FunctorL ArgOrder
deriving via (Base1 OptDescr) instance FunctorL OptDescr
deriving via (Base1 ArgDescr) instance FunctorL ArgDescr
deriving via (Base1 Identity) instance FunctorL Identity
deriving via (Base1 Handler) instance FunctorL Handler
deriving via (Base1 STM) instance FunctorL STM
deriving via (Base1 Monoid.First) instance FunctorL Monoid.First
deriving via (Base1 Monoid.Last) instance FunctorL Monoid.Last
deriving via (Base1 Monoid.Dual) instance FunctorL Monoid.Dual
deriving via (Base1 Monoid.Sum) instance FunctorL Monoid.Sum
deriving via (Base1 Monoid.Product) instance FunctorL Monoid.Product
deriving via (Base1 ZipList) instance FunctorL ZipList

deriving via (Base1 Generic.U1) instance FunctorL Generic.U1
deriving via (Base1 Generic.V1) instance FunctorL Generic.V1
instance (Prelude.Functor f) => FunctorL (Generic.Rec1 f) where { fmapL = Prelude.fmap }
deriving via (Base1 (Generic.URec (Ptr ()))) instance FunctorL (Generic.URec (Ptr ()))
deriving via (Base1 (Generic.URec Char)) instance FunctorL (Generic.URec Char)
deriving via (Base1 (Generic.URec Double)) instance FunctorL (Generic.URec Double)
deriving via (Base1 (Generic.URec Float)) instance FunctorL (Generic.URec Float)
deriving via (Base1 (Generic.URec Int)) instance FunctorL (Generic.URec Int)
deriving via (Base1 (Generic.K1 i c)) instance FunctorL (Generic.K1 i c)
deriving via (Base1 Proxy) instance FunctorL Proxy
deriving via (Base1 (Const m)) instance FunctorL (Const m)
instance (Functor f) => FunctorL (Monoid.Alt f) where { fmapL = fmap }
instance (Functor f) => FunctorL (Monoid.Ap f) where { fmapL = fmap }
deriving via (Base1 ((->) r)) instance FunctorL ((->) r)
instance (Functor f, Functor g) => FunctorL (f Generic.:+: g) where { fmapL = fmap }
instance (Functor f, Functor g) => FunctorL (f Generic.:*: g) where { fmapL = fmap }
instance FunctorL ((# , #) a) where { fmapL f (# a, b #) = (# a, f b #) }
instance FunctorL ((# ,, #) a b) where { fmapL f (# a, b, c #) = (# a, b, f c #) }

max_stars_repo_name: chessai/levity

max_issues_repo_name: chessai/levity

max_forks_repo_name: chessai/levity


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TupleSections         #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# OPTIONS_GHC -fno-warn-orphans  #-}

module Nudge (
    nudgeProp
  , nudgeProp2
  ) where

import           Control.Monad
import           Data.Bifunctor
import           Data.Finite
import           Data.Kind
import           Data.Maybe
import           Data.Proxy
import           GHC.TypeLits
import           Hedgehog
import           Lens.Micro
import           Lens.Micro.Platform                   ()
import           Numeric.Backprop
import qualified Data.Ix                               as Ix
import qualified Data.Vector.Sized                     as SV
import qualified Hedgehog.Gen                          as Gen
import qualified Hedgehog.Range                        as Range
import qualified Numeric.LinearAlgebra                 as HU
import qualified Numeric.LinearAlgebra.Static          as H
import qualified Numeric.LinearAlgebra.Static.Backprop as B

nudge :: Double
nudge = 1e-6

eps :: Double
eps = 1e-11

class (Backprop c, Show c, Show (TIx c)) => Testing c where
    type TIx c :: Type
    allIx  :: c -> [TIx c]
    ixLens :: TIx c -> Lens' c Double
    scalarize :: Reifies s W => BVar s c -> BVar s Double
    genTest :: Gen c

sized
    :: forall s t d. H.Sized t s d
    => Lens' s (d t)
sized f = fmap (fromJust . H.create) . f . H.extract

ixContainer
    :: forall t d. HU.Container d t
    => HU.IndexOf d
    -> Lens' (d t) t
ixContainer i = lens (`HU.atIndex` i)
                     (\xs x -> HU.accum xs (\_ _ -> x) [(i, x)])

instance Testing Double where
    type TIx Double = ()
    allIx _ = [()]
    ixLens _ = id
    scalarize = abs
    genTest = Gen.filter ((> eps) . (**2)) $
         Gen.double (Range.linearFracFrom 0 (-5) 5)

instance KnownNat n => Testing (H.R n) where
    type TIx (H.R n) = Int
    allIx v = [0 .. H.size v - 1]
    ixLens i = sized . ixContainer i
    scalarize = B.norm_2V
    genTest = H.vector <$> replicateM n genTest
      where
        n = fromInteger $ natVal (Proxy @n)

instance (KnownNat n, KnownNat m) => Testing (H.L n m) where
    type TIx (H.L n m) = (Int, Int)
    allIx m = Ix.range ((0,0), bimap pred pred (H.size m))
    ixLens i = sized . ixContainer i
    scalarize = sqrt . B.sumElements . (**2)
    genTest = H.matrix <$> replicateM nm genTest
      where
        nm = fromInteger $ natVal (Proxy @n) * natVal (Proxy @m)

instance Testing (HU.Vector Double) where
    type TIx (HU.Vector Double) = Int
    allIx v = [0 .. HU.size v - 1]
    ixLens = ixContainer
    scalarize = liftOp1 . op1 $ \xs -> (HU.sumElements xs, (`HU.konst` HU.size xs))
    genTest = HU.fromList <$> replicateM 3 genTest

instance Testing (HU.Matrix Double) where
    type TIx (HU.Matrix Double) = (Int, Int)
    allIx m = Ix.range ((0,0), bimap pred pred (HU.size m))
    ixLens = ixContainer
    scalarize = liftOp1 . op1 $ \xs -> (HU.sumElements xs, (`HU.konst` HU.size xs))
    genTest = HU.fromLists <$> (replicateM 3 . replicateM 2) genTest

instance (KnownNat n, Testing a, Num a) => Testing (SV.Vector n a) where
    type TIx (SV.Vector n a) = (Finite n, TIx a)
    allIx = fst . SV.imapM (\i x -> ((fromIntegral i,) <$> allIx x , x))
    ixLens (i,j) = SV.ix i . ixLens j
    scalarize = scalarize . liftOp1 o . (^ (2 :: Int))
      where
        o :: Op '[SV.Vector n a] a
        o = op1 $ \xs -> (SV.sum xs, SV.replicate)
    genTest = SV.replicateM genTest

instance (Testing a, Testing b) => Testing (a, b) where
    type TIx (a, b) = Either (TIx a) (TIx b)
    allIx (x, y) = (Left  <$> allIx x)
                ++ (Right <$> allIx y)
    ixLens (Left  i) = _1 . ixLens i
    ixLens (Right j) = _2 . ixLens j
    scalarize t = B.norm_2V (B.vec2 (scalarize (t ^^. _1))
                                    (scalarize (t ^^. _2))
                            )
    genTest = (,) <$> genTest <*> genTest

instance (Testing a, Testing b, Testing c) => Testing (a, b, c) where
    type TIx (a, b, c) = Either (TIx a) (Either (TIx b) (TIx c))
    allIx (x, y, z) = (Left          <$> allIx x)
                   ++ (Right . Left  <$> allIx y)
                   ++ (Right . Right <$> allIx z)
    ixLens (Left         i ) = _1 . ixLens i
    ixLens (Right (Left  j)) = _2 . ixLens j
    ixLens (Right (Right k)) = _3 . ixLens k
    scalarize t = B.norm_2V (B.vec3 (scalarize (t ^^. _1))
                                    (scalarize (t ^^. _2))
                                    (scalarize (t ^^. _3))
                            )
    genTest = (,,) <$> genTest <*> genTest <*> genTest

validGrad
    :: Monad m
    => Lens' c Double
    -> c
    -> c
    -> (c -> Double)
    -> PropertyT m (Double, Double)
validGrad l x0 g f = forAll $ Gen.double (Range.constantFrom 0 (-nudge) nudge) <&> \d ->
    let x   = x0 & l %~ (+d)
        old = f x0 + (g ^. l) * d
        new = f x
    in  (old, new)

nudgeProp
    :: forall c d. (Testing c, Testing d)
    => (forall s. Reifies s W => BVar s c -> BVar s d)
    -> Property
nudgeProp f = property $ do
    (inp, i) <- forAll $ do
      inp <- genTest
      i   <- Gen.element (allIx inp)
      return (inp, i)
    let (r,gr) = backprop (scalarize . f) inp
    when (r**2 < eps) discard
    (old, new) <- validGrad (ixLens i) inp gr (evalBP (scalarize . f))
    footnoteShow (r, gr, old, new, (old - new)**2, ((old - new)/old)**2)
    assert $ ((old - new)/old)**2 < eps

nudgeProp2
    :: forall c d e. (Testing c, Testing d, Testing e)
    => (forall s. Reifies s W => BVar s c -> BVar s d -> BVar s e)
    -> Property
nudgeProp2 f = property $ do
    (inpC, inpD, i) <- forAll $ do
      inpC <- genTest
      inpD <- genTest
      i    <- Gen.element (allIx (inpC, inpD))
      return (inpC, inpD, i)
    let (r, gr) = backprop2 (\x -> scalarize . f x) inpC inpD
    when (r**2 < eps) discard
    (old, new) <- validGrad (ixLens i) (inpC, inpD) gr
          (evalBP (\t -> scalarize $ f (t ^^. _1) (t ^^. _2)))
    footnoteShow (r, gr, old, new, (old - new)**2, ((old - new)/old)**2)
    assert $ ((old - new)/old)**2 < eps

instance (HU.Container HU.Vector a, Num a) => Backprop (HU.Matrix a) where
    -- TODO: make more efficient?
    zero = HU.cmap (const 0)
    add  = HU.add
    one  = HU.cmap (const 1)

instance (KnownNat n, Num a) => Backprop (SV.Vector n a) where
    zero = (0 <$)
    add  = (+)
    one  = (1 <$)

max_stars_repo_name: mstksg/hmatrix-backprop

max_issues_repo_name: mstksg/hmatrix-backprop

max_forks_repo_name: mstksg/hmatrix-backprop


text:
{-# LANGUAGE NamedFieldPuns #-}

module Zambini where

import Data.Complex
import Data.Fixed (mod')

import Graphics.Gloss.Interface.Pure.Simulate

import ViewUtil

type R = Double
type C = Complex Double

data Zambini = Zambini { m :: R -- physical attributes
                       , i :: R
                       , r :: R
                       , fMax :: R
                       -- input
                       , f1 :: R
                       , f2 :: R
                       -- state
                       , x :: C
                       , v :: C
                       , theta :: R
                       , omega :: R
                       }

step :: R -> Zambini -> Zambini
step dt c@Zambini{m, i, r, fMax, f1, f2, x, v, theta, omega} = c{x=x', v=v', theta=theta', omega=omega'}
  where
    f1' = min 1 f1
    f2' = min 1 f2
    a = mkPolar ((f1' + f2') / m) theta
    v' = v + a * realToFrac dt
    x' = x + v' * realToFrac dt
    torque = r * (f1' - f2')
    omega' = omega + torque / i * dt
    theta' = theta + omega' * dt `mod'` 2*pi

render :: Zambini -> Picture
render Zambini{f1, f2, x, theta} = move x $ rotateRad (realToFrac theta) $ pictures [zambiniPic, control1 f1, control2 f2]
  where
    zambiniPic = color (greyN 0.5) $ pictures [ polygon [(-5,0), (-10, 25), (10, 25), (5, 0), (10, -25), (-10, -25)]
                                              , polygon [(0, 10), (20, 0), (0, -10)]
                                              ]

    control = rotate 180 . arrow (5, 20) (15, 100) (low, high)
      where
        low  = makeColorI 0 64 192 255
        high = makeColorI 0 160 192 255

    control1 = translate (-20) (-20) . control
    control2 = translate (-20) (20) . control

    move (x :+ y) = translate (realToFrac x) (realToFrac y)

max_stars_repo_name: tyehle/control-systems

max_issues_repo_name: tyehle/control-systems

max_forks_repo_name: tyehle/control-systems


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards, PatternGuards, BangPatterns #-}
-----------------------------------------------------------------------------
--
-- Module      :  AI.Network.RNN.RNN
-- Copyright   :  (c) JP Moresmau
-- License     :  BSD3
--
-- Maintainer  :  JP Moresmau <jp@moresmau.fr>
-- Stability   :  experimental
-- Portability :
--
-- | Recurrent neural network
--
-----------------------------------------------------------------------------

module AI.Network.RNN.RNN  where

import Numeric.LinearAlgebra.HMatrix hiding ((|>))
import Control.DeepSeq
import Control.Monad.Random hiding (fromList)

import AI.Network.RNN.Types
import AI.Network.RNN.Util

data RNNDimensions = RNNDimensions
  { rnndInput    :: Int
  , rnndInternal :: Int
  , rnndOutput   :: Int
  , rnndWithBack :: Bool
  }  deriving (Show,Read,Eq,Ord)

instance NFData RNNDimensions where
    rnf RNNDimensions{..} = rnf (rnndInput,rnndInternal,rnndOutput,rnndWithBack)

data RNNetwork = RNNetwork
  { rnnDimensions :: RNNDimensions
  , rnnMIn        :: Matrix Double
  , rnnM          :: Matrix Double
  , rnnMOut       :: Matrix Double
  , rnnMBack      :: Maybe (Matrix Double)
  , rnnState      :: Vector Double
  , rnnOutput     :: Vector Double
  } deriving (Show,Read,Eq)

networkDimensions :: RNNetwork -> RNNDimensions
networkDimensions = rnnDimensions

instance NFData RNNetwork where
    rnf RNNetwork{..} = rnf (rnnDimensions,rnnMIn,rnnM,rnnMOut,rnnMBack,rnnState,rnnOutput)

createNetwork
    :: RNNDimensions
    -> Matrix Double
    -> Matrix Double
    -> Matrix Double
    -> Maybe (Matrix Double)
    -> Vector Double
    -> Vector Double
    -> Either [String] RNNetwork
createNetwork dim mIn m mOut mmback st out =
    case checkNetworkDimensions dim mIn m mOut mmback st of
        []  -> Right $ RNNetwork dim mIn m mOut mmback st out
        err -> Left err

createNetworkFromArray
    :: RNNDimensions
    -> [Double]
    -> Either String RNNetwork
createNetworkFromArray dim@RNNDimensions{..} v =
    let
        enough = length v >= totalDataLength dim
    in if enough
        then Right $ fromVector dim $ fromList v
        else Left $ "Not enough data in array (needs at least"++show (totalDataLength dim)++")"

networkToArray :: RNNetwork -> [Double]
networkToArray RNNetwork{..} =
       concat (toLists rnnMIn)
    ++ concat (toLists rnnM)
    ++ concat (toLists rnnMOut)
    ++ (case rnnMBack of
            Just b -> concat $ toLists b
            _      -> [])
    ++ toList rnnState
    ++ toList rnnOutput

createNetworkFromVector
    :: RNNDimensions
    -> Vector Double
    -> Either String RNNetwork
createNetworkFromVector dim@RNNDimensions{..} v =
    let
        enough = size v >= totalDataLength dim
    in if enough
        then {-# SCC "createNetworkFromVector" #-} Right $ fromVector dim v
        else Left $ "Not enough data in vector (needs at least"++ show (totalDataLength dim) ++")"

inputMatrixLength :: RNNDimensions -> Int
inputMatrixLength RNNDimensions{..} = rnndInput * rnndInternal

internalMatrixLength :: RNNDimensions -> Int
internalMatrixLength RNNDimensions{..} = rnndInternal * rnndInternal

outputMatrixLength :: RNNDimensions -> Int
outputMatrixLength RNNDimensions{..} = rnndOutput * (rnndInput + rnndInternal)

backMatrixLength :: RNNDimensions -> Int
backMatrixLength RNNDimensions{..} = if rnndWithBack then rnndInternal * rnndOutput else 0

totalDataLength :: FullSize RNNetwork
totalDataLength dim@RNNDimensions{..} =
    inputMatrixLength dim + internalMatrixLength dim + outputMatrixLength dim
    + backMatrixLength dim
    + rnndInternal + rnndOutput

randNetwork :: (Monad m,RandomGen g) =>  RNNDimensions -> RandT g m RNNetwork
randNetwork dim@RNNDimensions{..} = do
    s <- getRandom
    -- vs <- sequence (replicate (totalDataLength dim) getRandom)
    return $ fromVector dim $ randomVector s Gaussian (totalDataLength dim)

collectErrors :: [(Bool,a)] ->  [a]
collectErrors = map snd . filter fst

checkDimensions
    :: RNNDimensions
    -> [String]
checkDimensions RNNDimensions{..} = collectErrors
    [(rnndInput    < 0, "input must be >=0")
    ,(rnndInternal < 1, "internal must be >0")
    ,(rnndOutput   < 1, "output must be >0")]

checkNetworkDimensions
    :: RNNDimensions
    -> Matrix Double
    -> Matrix Double
    -> Matrix Double
    -> Maybe (Matrix Double)
    -> Vector Double
    -> [String]
checkNetworkDimensions dim@RNNDimensions{..} mIn m mOut mmback st = checkDimensions dim ++ collectErrors (
    [(cols mIn /= rnndInput, "input matrix column count doesn't match input neurons count")
    ,(rows mIn /= rnndInternal, "input matrix row count doesn't match internal neurons count")
    ,(cols m /= rnndInternal, "internal matrix column count doesn't match internal neurons count")
    ,(rows m /= rnndInternal, "internal matrix row count doesn't match internal neurons count")
    ,(cols mOut /= (rnndInput + rnndInternal), "internal matrix column count doesn't match input + internal neurons count")
    ,(rows mOut /= rnndOutput, "output matrix row count doesn't match output neurons count")
    ,(size st /= rnndInternal,"internal state length doesn't match internal neurons count")]
    ++(case mmback of
        Nothing -> []
        Just mback ->
            [(rows mback /= rnndInternal,"back matrix row count doesn't match internal neurons count")
            ,(cols m /= rnndOutput,"back matrix column count doesn't match output neurons count")
            ]))

instance RNNEval RNNetwork where
    type Size RNNetwork = RNNDimensions
    evalStep rnn@RNNetwork{..} iv =
        let
            sum1 = (rnnMIn #> iv) + (rnnM #> rnnState)
            sum2 = case rnnMBack of
                Just mback -> sum1 + (mback #> rnnOutput)
                _          -> sum1
            s2 = cmap tanh sum2
            out = cmap sigmoid (rnnMOut #> vjoin [iv,s2])
        in (rnn{rnnState=s2,rnnOutput=out},out)
    toVector RNNetwork{..} = vjoin
        [ flatten rnnMIn
        , flatten rnnM
        , flatten rnnMOut
        , maybe (fromList []) flatten rnnMBack
        , rnnState
        , rnnOutput
        ]
    fromVector dim@RNNDimensions{..} v =
        let
            lenIn   = inputMatrixLength dim
            len     = internalMatrixLength dim
            lenOut  = outputMatrixLength dim
            lenBack = backMatrixLength dim
            [v1,v2,v3,v4,v5,v6] = takesV [lenIn,len,lenOut,lenBack,rnndInternal,rnndOutput] v
            mIn     = reshape rnndInput v1
            m       = reshape rnndInternal v2
            mOut    = reshape (rnndInput + rnndInternal) v3
            mback = if rnndWithBack
                            then Just $ reshape rnndOutput v4
                            else Nothing
        in RNNetwork dim mIn m mOut mback v5 v6
    rnnsize = rnnDimensions
    fullSize = totalDataLength . rnnDimensions

max_stars_repo_name: JPMoresmau/rnn

max_issues_repo_name: JPMoresmau/rnn

max_forks_repo_name: JPMoresmau/rnn


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ConstraintKinds #-}
module Noise(
      histogram
    , Structure(..)
    , quantizationNoise
    ) where

import qualified Statistics.Sample.Histogram as H
import qualified Data.Vector.Unboxed as U
import Plot
import Graphics.PDF
import Common
import Generators
import Signal
import System.Random
import Text.Printf

nbBins = 20

drawHist :: Int
         -> U.Vector Double
         -> U.Vector Double
         -> Int
         -> Int
         -> CoordinateMapping Double Double
         -> Draw ()
drawHist nb v b wi hi (toPoint,_) = do
    let binWidth :: Double
        binWidth = fromIntegral wi / fromIntegral nb
        d = binWidth / 2.0
        drawBin p = do
           fillColor (Rgb 0.8 0.8 1.0)
           let (x :+ y) = toPoint p
               r = Rectangle ((x - d) :+ 0) ((x + d) :+ y)
           fill r
           stroke r
           return ()
    mapM_ drawBin $ (U.toList $ U.zip v b)

histogram :: (U.Unbox a, HasDoubleRepresentation a)
          => [a]
          -> StyledSignal Double Double
histogram l = do
    let (s,v1) = H.histogram nbBins (U.fromList . map toDouble $ l)
        v = U.map (/ (U.sum v1)) v1
        nb = U.length s
        mas = U.maximum s
        mis = U.minimum s
        d = (mas - mis) / fromIntegral nb / 2.0
        detailed x = printf "%.2g" x
        style = defaultPlotStyle { title = Just "Histogram"
                                 , horizontalLabel = Nothing
                                 , verticalLabel = Nothing
                                 , horizontalBounds = Just (mis-d,mas+d)
                                 , verticalBounds = Just (U.minimum v, U.maximum v)
                                 , epilog = drawHist nb s v
                                 , axis = False
                                 , horizontalTickRepresentation = detailed
                                 }
    discreteSignalsWithStyle (U.length s) style []

class Structure m where
    doubleVersion :: (Sample i, Sample o) => m f i o -> m Double Double Double
    transferFunction :: (Sample i, Sample o) => m f i o -> Signal i -> Signal o

quantizationNoise :: (Structure m , Random i, Sample o, Sample i)
                  => Signal i
                  -> m p i o
                  -> IO (Signal Double)
quantizationNoise r structure = do
    let ds = transferFunction (doubleVersion structure) (mapS toDouble r)
        fs = mapS toDouble . transferFunction structure $ r
    return $ zipWithS (-) ds fs

max_stars_repo_name: alpheccar/Signal

max_issues_repo_name: cpehle/Signal

max_forks_repo_name: cpehle/Signal


text:
import Math.Combinat.Partitions (_partitions)
import Numeric.LinearAlgebra
import System.IO (putStrLn)
import Text.Printf (printf)

getX :: Matrix Double -> Matrix Double
getX = takeColumns 2

getY :: Matrix Double -> Vector Double
getY = head . toColumns . dropColumns 2

nonLinear :: Matrix Double -> Matrix Double
nonLinear = fromRows . map (fromList . nonLinearRow . toList) . toRows
  where
    nonLinearRow [x1, x2] = [1, x1, x2, x1^2, x2^2, x1 * x2, abs (x1 - x2), abs (x1 + x2)]

classificationError :: Vector Double -> Vector Double -> Double
classificationError v y = let sv = mapVector signum v
                              diffs = zipVectorWith (((fromIntegral . fromEnum).) . (/=)) sv y
                          in  foldVector (+) 0 diffs / fromIntegral (dim y)

regularizedRegression :: Matrix Double -> Vector Double -> Double -> Vector Double
regularizedRegression x y lambda = let xt = trans x
                                       xtx = multiply xt x
                                       n = cols x
                                       lambdaI = diag (constant lambda n)
                                       z = multiply (pinv (xtx + lambdaI)) xt
                                   in  z `mXv` y

experiment :: Matrix Double -> Vector Double -> Matrix Double -> Vector Double -> Double -> (Double, Double)
experiment x y x' y' lambda = let z = regularizedRegression x y lambda
                                  xz = x `mXv` z
                                  ein = classificationError xz y
                                  x'z = x' `mXv` z
                                  eout = classificationError x'z y'
                              in  (ein, eout)

main = do
         inputMatrix <- fromFile "in.dta" (35, 3)
         testMatrix <- fromFile "out.dta" (250, 3)
         let x = nonLinear (getX inputMatrix)
             y = getY inputMatrix
             x' = nonLinear (getX testMatrix)
             y' = getY testMatrix
         let run = experiment x y x' y'

         let (ein, eout) = run 0
         putStrLn $ printf "2)\tE_{in} = %.5f\n\tE_{out} = %.5f" ein eout

         let (ein, eout) = run (10**(-3))
         putStrLn $ printf "3)\tE_{in} = %.5f\n\tE_{out} = %.5f" ein eout

         let (ein, eout) = run (10**3)
         putStrLn $ printf "4)\tE_{in} = %.5f\n\tE_{out} = %.5f" ein eout

         let lambdas = [2, 1, 0, -1, -2] :: [Double]
             eouts = map (snd . run . (10**)) lambdas
             best = minimum (zip eouts lambdas)
         putStrLn $ printf "5)\t%d" (floor $ snd best :: Int)
         putStrLn $ printf "6)\t%f" (fst best)

         let layerings = filter (all (>= 2)) $ _partitions 36
             cost p = cost' $ 10:p
               where
                cost' [x] = x
                cost' (x:y:xs) = x * (y - 1) + cost' (y:xs)
             costs = map cost layerings
         putStrLn $ printf "9)\t%d" (minimum costs)
         putStrLn $ printf "10)\t%d" (maximum costs)


text:
{-# Language RankNTypes #-}

module Parsing where

import Prelude hiding (negate)

-- Inspiration from
-- https://github.com/justinethier/husk-scheme/blob/master/hs-src/Language/Scheme/Parser.hs

import Text.ParserCombinators.Parsec hiding (spaces)
import Text.Parsec.Language
import qualified Text.Parsec.Token as P
import Data.Functor.Identity (Identity)
import Text.Parsec.Prim (ParsecT)

import qualified Data.Char as DC
import Numeric
import Data.Complex
import Data.Ratio
import Data.Array
import Control.Monad.Except

import qualified Data.Text as T

import Types

style :: LanguageDef ()
style =
  emptyDef
    { P.commentStart = "#|"
    , P.commentEnd = "|#"
    , P.commentLine = ";"
    , P.nestedComments = True
    , P.identStart = letter <|> symbol
    , P.identLetter = letter <|> digit <|> symbol
    , P.reservedNames = []
    , P.caseSensitive = True
    }

lexer :: P.GenTokenParser String () Data.Functor.Identity.Identity
lexer = P.makeTokenParser style

dot :: ParsecT String () Identity String
dot = P.dot lexer

parens :: ParsecT String () Identity a -> ParsecT String () Identity a
parens = P.parens lexer

braces :: ParsecT String () Identity a -> ParsecT String () Identity a
braces = P.braces lexer

brackets :: ParsecT String () Identity a -> ParsecT String () Identity a
brackets = P.brackets lexer

bars :: ParsecT String () Identity a -> ParsecT String () Identity a
bars = try . between (symbol "|") (symbol "|")
  where
    symbol name = lexeme (string name)
    lexeme p = do
      x <- p
      whiteSpace
      return x

identifier :: ParsecT String () Identity String
identifier = P.identifier lexer

whiteSpace :: ParsecT String () Identity ()
whiteSpace = P.whiteSpace lexer

lexeme :: ParsecT String () Identity a -> ParsecT String () Identity a
lexeme = P.lexeme lexer

symbol :: Parser Char
symbol = oneOf "!$%&*+-/:<=>?@^_~."

-- |Parse an atom (scheme symbol)
parseAtom :: Parser LispVal
parseAtom = do
  atom <- identifier
  if atom == "."
     then pzero -- Do not match this form
     else return $ Atom $ T.pack atom

parseBool :: Parser LispVal
parseBool = do
  _ <- char '#'
  (char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False))

parseEscapedChar :: forall st .
                    GenParser Char st Char
parseEscapedChar = do
  _ <- char '\\'
  c <- anyChar
  case c of
    'a' -> return '\a'
    'b' -> return '\b'
    'n' -> return '\n'
    't' -> return '\t'
    'r' -> return '\r'
    'x' -> do
        num <- many $ letter <|> digit
        _ <- char ';'
        parseHexScalar num
    _ -> return c

parseString :: Parser LispVal
parseString = do
  _ <- char '"'
  x <- many (parseEscapedChar <|> noneOf "\"")
  _ <- char '"'
  return $ String $ T.pack x

parseChar :: Parser LispVal
parseChar = do
  _ <- try (string "#\\")
  c <- anyChar
  r <- many (letter <|> digit)
  let pchr = c : r
  case pchr of
    "space"     -> return $ Char ' '
    "newline"   -> return $ Char '\n'
    "alarm"     -> return $ Char '\a'
    "backspace" -> return $ Char '\b'
    "delete"    -> return $ Char '\DEL'
    "escape"    -> return $ Char '\ESC'
    "null"      -> return $ Char '\0'
    "return"    -> return $ Char '\n'
    "tab"       -> return $ Char '\t'
    _ -> case (c : r) of
        [ch] -> return $ Char ch
        ('x' : hexs) -> do
            rv <- parseHexScalar hexs
            return $ Char rv
        _ -> pzero

parseHexScalar :: String -> GenParser Char st Char
parseHexScalar num = do
    let ns = Numeric.readHex num
    case ns of
        [] -> fail $ "Unable to parse hex value " ++ show num
        _ -> return $ DC.chr $ fst $ head ns

parseNumber :: Parser LispVal
parseNumber = parseDecimal1 <|> parseDecimal2 <|> parseHex <|> parseOct <|> parseBin

negate :: Num a => String -> a -> a
negate "-" num = -num
negate "+" num = num
negate "" num  = num
negate _ _ = undefined

parseDecimal1 :: Parser LispVal
parseDecimal1 = do
  sign <- many (oneOf "-+")
  x <- many1 digit
  (return . Number . negate sign . read) x

parseDecimal2 :: Parser LispVal
parseDecimal2 = do
  _ <- try $ string "#d"
  sign <- many (oneOf "-+")
  x <- many1 digit
  (return . Number . negate sign . read) x

parseHex :: Parser LispVal
parseHex = do
  _ <- try $ string "#x"
  sign <- many (oneOf "-+")
  x <- many1 hexDigit
  return $ Number $ negate sign (hex2dig x)

parseOct :: Parser LispVal
parseOct = do
  _ <- try $ string "#o"
  sign <- many (oneOf "-+")
  x <- many1 octDigit
  return $ Number $ negate sign (oct2dig x)

parseBin :: Parser LispVal
parseBin = do
  _ <- try $ string "#b"
  sign <- many (oneOf "-+")
  x <- many1 (oneOf "10")
  return $ Number $ negate sign (bin2dig x)

oct2dig :: (Eq a, Num a) => String -> a
oct2dig x = fst $ head (readOct x)
hex2dig :: (Eq a, Num a) => String -> a
hex2dig x = fst $ head (readHex x)
bin2dig :: String -> Integer
bin2dig  = bin2dig' 0
bin2dig' :: Num t => t -> String -> t
bin2dig' digint "" = digint
bin2dig' digint (x:xs) = let old = 2 * digint + (if x == '0' then 0 else 1) in
                         bin2dig' old xs

parseFloat :: Parser LispVal
parseFloat = do
  sign <- many (oneOf "-+")
  x <- many1 digit
  _ <- char '.'
  y <- many1 digit
  return $ Float $ negate sign  (fst . head $ readFloat (x ++ "." ++ y))

parseRatio :: Parser LispVal
parseRatio = do
  sign <- many (oneOf "-+")
  x <- many1 digit
  _ <- char '/'
  y <- many1 digit
  return $ Ratio (negate sign (read x) % read y)

toDouble :: LispVal -> Double
toDouble (Float f) = realToFrac f
toDouble (Number n) = fromIntegral n
toDouble _ = undefined

parseComplex :: Parser LispVal
parseComplex = do
  x <- try parseFloat <|> parseNumber
  _ <- char '+'
  _ <- many $ char '('
  y <- try parseFloat <|> parseNumber
  _ <- many $ char ')'
  _ <- char 'i'
  return $ Complex (toDouble x :+ toDouble y)

parseList :: Parser LispVal
parseList = List <$> sepBy parseExpr whiteSpace

parseDottedList :: Parser LispVal
parseDottedList = do
  phead <- endBy parseExpr whiteSpace
  case phead of
    [] -> pzero -- car is required; no match
    _ -> do
      ptail <- dot >> parseExpr
      case ptail of
        DottedList ls l -> return $ DottedList (phead ++ ls) l
        -- Issue #41
        -- Improper lists are tricky because if an improper list ends in a
        -- proper list, then it becomes proper as well. The following cases
        -- handle that, as well as preserving necessary functionality when
        -- appropriate, such as for unquoting.
        --
        -- FUTURE: I am not sure if this is complete, in fact the "unquote"
        -- seems like it could either be incorrect or one special case among
        -- others. Anyway, for the 3.3 release this is good enough to pass all
        -- test cases. It will be revisited later if necessary.
        --
        List (Atom "unquote" : _) -> return $ DottedList phead ptail
        List ls -> return $ List $ phead ++ ls
        {- Regarding above, see
           http://community.schemewiki.org/?scheme-faq-language#dottedapp

           Note, however, that most Schemes expand literal lists occurring in
           function applications, e.g. (foo bar . (1 2 3)) is expanded into
           (foo bar 1 2 3) by the reader. It is not entirely clear whether this
           is a consequence of the standard - the notation is not part of the
           R5RS grammar but there is strong evidence to suggest a Scheme
           implementation cannot comply with all of R5RS without performing this
           transformation. -}
        _ -> return $ DottedList phead ptail

parseQuoted :: Parser LispVal
parseQuoted = do
    _ <- lexeme $ char '\''
    x <- parseExpr
    return $ List [Atom "quote", x]

parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
  _ <- lexeme $ char '`'
  x <- parseExpr
  return $ List [Atom "quasiquote", x]

parseUnquoted :: Parser LispVal
parseUnquoted = do
  _ <- try (lexeme $ char ',')
  x <- parseExpr
  return $ List [Atom "unquote", x]

parseUnquoteSpliced :: Parser LispVal
parseUnquoteSpliced = do
  _ <- try (lexeme $ string ",@")
  x <- parseExpr
  return $ List [Atom "unquote-splicing", x]

-- |Parse a vector
parseVector :: Parser LispVal
parseVector = do
  vals <- sepBy parseExpr whiteSpace
  return $ Vector (listArray (0, length vals - 1) vals)

parseExpr :: Parser LispVal
parseExpr =
      try (lexeme parseComplex)
  <|> try (lexeme parseRatio)
  <|> try (lexeme parseFloat)
  <|> try (lexeme parseNumber)
  <|> lexeme parseChar
  <|> parseUnquoteSpliced
  <|> do _ <- try (lexeme $ string "#(")
         x <- parseVector
         _ <- lexeme $ char ')'
         return x
  -- <|> do _ <- try (lexeme $ string "#u8(")
  --        x <- parseByteVector
  --        _ <- lexeme $ char ')'
  --        return x
--  <|> do _ <- try (lexeme $ string "#hash(")
--         x <- parseHashTable
--         _ <- lexeme $ char ')'
--         return x
  <|> try parseNil
  <|> try parseAtom
  <|> lexeme parseString
  <|> lexeme parseBool
  <|> parseQuoted
  <|> parseQuasiQuoted
  <|> parseUnquoted
  <|> try (parens parseList)
  <|> parens parseDottedList
  <|> try (brackets parseList)
  <|> brackets parseDottedList
  <|> try (braces parseList)
  <|> braces parseDottedList
  <|> try (bars parseList)
  <|> bars parseDottedList
  <?> "Expression"

parseNil :: Parser LispVal
parseNil = do
  string "nil"
  return $ List []

-- |Initial parser used by the high-level parse functions
mainParser :: Parser LispVal
mainParser = do
    _ <- whiteSpace
    parseExpr

-- |Use a parser to parse the given text, throwing an error
--  if there is a problem parsing the text.
readOrThrow :: Parser a -> String -> ThrowsError a
readOrThrow parser input = case parse parser "lisp" input of
  Left err -> throwError $ Parser err
  Right val -> return val

-- |Parse an expression from a string of text
readExpr :: String -> ThrowsError LispVal
readExpr = readOrThrow mainParser

-- |Parse many expressions from a string of text
readExprList :: String -> ThrowsError [LispVal]
readExprList = readOrThrow (endBy mainParser whiteSpace)

max_stars_repo_name: guimcaballero/scheme-hs

max_issues_repo_name: guimcaballero/scheme-hs

max_forks_repo_name: guimcaballero/scheme-hs


text:
module ActivationFunction
    ( step'
    , sigmoid
    , relu
    , softMax'
    ) where

import Numeric.LinearAlgebra

step' :: (Element b, Container c a, Ord a, Num a, Num b) => c a -> c b
step' = cmap (\x -> if x>0 then 1 else 0)

sigmoid :: (Container c b, Floating b) => c b -> c b
sigmoid = cmap (\x -> 1 / (1 + exp(-x)))

relu :: (Container c b, Ord b, Num b) => c b -> c b
relu = cmap (max 0)

-- This function is deprecated, causes overflow.
softMax :: (Container c b, Floating b) => c b -> c b
softMax xs = cmap (/sumCexp) $ cexp xs
    where sumCexp = (sumElements . cexp) xs

softMax' :: (Container c b, Floating b) => c b -> c b
softMax' xs = cmap (/sumCexp) $ cexp' xs
    where sumCexp = (sumElements . cexp') xs

cexp :: (Container c b, Floating b) => c b -> c b
cexp = cmap exp

cexp' :: (Container c b, Floating b) => c b -> c b
cexp' xs = cexp $ cmap (subtract $ maxElement xs) xs


text:
module STCR2Z2T0S0KoffkaCross where

import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           DFT.Plan
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           STC.PowerMethod
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoDecayStr:taoReversalStr:taoCornerStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:histFilePath:numIterationStr:writeSourceFlagStr:dStr:wStr:sigmaStr:thetaStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      taoDecay = read taoDecayStr :: Double
      taoReversal = read taoReversalStr :: Double
      taoCorner = read taoCornerStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      d = read dStr :: Int
      w = read wStr :: Int
      sigma = read sigmaStr :: Double
      theta = read thetaStr :: Double
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2Z2T0S0KoffkaCross"
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0ReversalCornerRadial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          taoDecay
          taoReversal
          taoCorner
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      radialArr
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  plan <- makeR2Z2T0S0Plan emptyPlan arrR2Z2T0S0
  let xs =
        [ R2S1RPPoint (d, w, theta, 1)
        , R2S1RPPoint (d, -w, theta, 1)
        , R2S1RPPoint (w, d, theta + 90, 1)
        , R2S1RPPoint (-w, d, theta + 90, 1)
        , R2S1RPPoint (-d, w, theta + 180, 1)
        , R2S1RPPoint (-d, -w, theta + 180, 1)
        , R2S1RPPoint (w, -d, theta + 270, 1)
        , R2S1RPPoint (-w, -d, theta + 270, 1)
        ]
      bias =
        -- computeBiasR2T0S0Gaussian
        --   numPoint
        --   numPoint
        --   theta0Freqs
        --   scale0Freqs
        --   90
        --   sigma
        --   xs
        computeS $
        R.zipWith
          (+)
          (computeBiasR2T0S0Gaussian
             numPoint
             numPoint
             theta0Freqs
             scale0Freqs
             90
             sigma
             xs)
          (computeBiasR2T0S0Gaussian
             numPoint
             numPoint
             theta0Freqs
             scale0Freqs
             (-90)
             sigma
             xs)
      eigenVec =
        computeInitialEigenVectorR2T0S0
          numPoint
          numPoint
          theta0Freqs
          scale0Freqs
          thetaFreqs
          scaleFreqs
          xs
  powerMethodR2Z2T0S0Bias
    plan
    folderPath
    numPoint
    numPoint
    numOrientation
    thetaFreqs
    theta0Freqs
    numScale
    scaleFreqs
    scale0Freqs
    arrR2Z2T0S0
    numIteration
    writeSourceFlag
    -- (printf "_%d" (round r :: Int))
    (printf
       "_%d_%d_%d_%d_%d_%.2f_%.2f_%d_%d"
       numPoint
       (round maxScale :: Int)
       (round taoDecay :: Int)
       (round taoReversal :: Int)
       (round taoCorner :: Int)
       thetaSigma
       scaleSigma
       d
       w)
    0.5
    bias
    eigenVec

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:
-- |
module Main where

import Control.Monad
import Control.Monad.ST (runST)

import qualified Data.Vector.Unboxed as U

import Criterion.Main
import System.Random.MWC
--
import Numeric.Sum
import Statistics.Monoid
import Statistics.Monoid.Extra

sampleD10,sampleD100,sampleD1000 :: U.Vector Double
sampleD10   = runST $ U.replicateM 10   . uniform =<< create
sampleD100  = runST $ U.replicateM 100  . uniform =<< create
sampleD1000 = runST $ U.replicateM 1000 . uniform =<< create

sampleI10,sampleI100,sampleI1000 :: U.Vector Int
sampleI10   = runST $ U.replicateM 10   . uniform =<< create
sampleI100  = runST $ U.replicateM 100  . uniform =<< create
sampleI1000 = runST $ U.replicateM 1000 . uniform =<< create

main :: IO ()
main = defaultMain
  [ bgroup "Count"
      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> Count) sampleD10
      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> Count) sampleD100
      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> Count) sampleD1000
      ]
  , bgroup "KBN"
      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> KBNSum) sampleD10
      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> KBNSum) sampleD100
      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> KBNSum) sampleD1000
      ]
  , bgroup "Kahan"
      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> KahanSum) sampleD10
      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> KahanSum) sampleD100
      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> KahanSum) sampleD1000
      ]
  , bgroup "MeanKBN"
      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> MeanKBN) sampleD10
      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> MeanKBN) sampleD100
      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> MeanKBN) sampleD1000
      , bench "10   I" $ whnf (reduceSampleVec :: U.Vector Int    -> MeanKBN) sampleI10
      , bench "100  I" $ whnf (reduceSampleVec :: U.Vector Int    -> MeanKBN) sampleI100
      , bench "1000 I" $ whnf (reduceSampleVec :: U.Vector Int    -> MeanKBN) sampleI1000
      ]
  , bgroup "WelfordMean"
      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> WelfordMean) sampleD10
      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> WelfordMean) sampleD100
      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> WelfordMean) sampleD1000
      , bench "10   I" $ whnf (reduceSampleVec :: U.Vector Int    -> WelfordMean) sampleI10
      , bench "100  I" $ whnf (reduceSampleVec :: U.Vector Int    -> WelfordMean) sampleI100
      , bench "1000 I" $ whnf (reduceSampleVec :: U.Vector Int    -> WelfordMean) sampleI1000
      ]
  , bgroup "Variance"
      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> Variance) sampleD10
      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> Variance) sampleD100
      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> Variance) sampleD1000
      , bench "10   I" $ whnf (reduceSampleVec :: U.Vector Int    -> Variance) sampleI10
      , bench "100  I" $ whnf (reduceSampleVec :: U.Vector Int    -> Variance) sampleI100
      , bench "1000 I" $ whnf (reduceSampleVec :: U.Vector Int    -> Variance) sampleI1000
      ]
  ]