text:
{-# LANGUAGE ConstraintKinds       #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE IncoherentInstances   #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeSynonymInstances  #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# OPTIONS -fno-warn-orphans #-}

module Data.Array.Accelerate.Data.Complex (

  -- * Rectangular from
  Complex(..),
  real,
  imag,

  -- * Polar form
  mkPolar,
  cis,
  polar,
  magnitude,
  phase,

  -- * Conjugate
  conjugate,

) where

import Prelude
import Data.Complex                             ( Complex(..) )
import Data.Array.Accelerate                    as A
import Data.Array.Accelerate.Smart
import Data.Array.Accelerate.Product
import Data.Array.Accelerate.Array.Sugar

type instance EltRepr (Complex a) = EltRepr (a, a)

instance Elt a => Elt (Complex a) where
  eltType _             = eltType (undefined :: (a,a))
  toElt p               = let (a, b) = toElt p in a :+ b
  fromElt (a :+ b)      = fromElt (a, b)

instance cst a => IsProduct cst (Complex a) where
  type ProdRepr (Complex a) = ProdRepr (a, a)
  fromProd cst (x :+ y) = fromProd cst (x, y)
  toProd cst p          = let (x, y) = toProd cst p in (x :+ y)
  prod cst _            = prod cst (undefined :: (a, a))

instance (Lift Exp a, Elt (Plain a)) => Lift Exp (Complex a) where
  type Plain (Complex a) = Complex (Plain a)
  lift (x1 :+ x2)       = Exp $ Tuple (NilTup `SnocTup` lift x1 `SnocTup` lift x2)

instance Elt a => Unlift Exp (Complex (Exp a)) where
  unlift e
    = let x     = Exp $ SuccTupIdx ZeroTupIdx `Prj` e
          y     = Exp $ ZeroTupIdx `Prj` e
      in
      x :+ y

instance (Elt a, IsFloating a) => Num (Exp (Complex a)) where
  (+)           = lift2 ((+) :: Complex (Exp a) -> Complex (Exp a) -> Complex (Exp a))
  (-)           = lift2 ((-) :: Complex (Exp a) -> Complex (Exp a) -> Complex (Exp a))
  (*)           = lift2 ((*) :: Complex (Exp a) -> Complex (Exp a) -> Complex (Exp a))
  negate        = lift1 (negate :: Complex (Exp a) -> Complex (Exp a))
  signum        = lift1 (signum :: Complex (Exp a) -> Complex (Exp a))
  abs           = lift1 (abs :: Complex (Exp a) -> Complex (Exp a))
  fromInteger n = lift (constant (fromInteger n) :+ 0)

instance (Elt a, IsFloating a) => Fractional (Exp (Complex a)) where
  c / c'
    = let x  :+ y       = unlift c
          x' :+ y'      = unlift c'     :: Complex (Exp a)
          den           = x'^(2 :: Int) + y'^(2 :: Int)
          re            = (x * x' + y * y') / den
          im            = (y * x' - x * y') / den
      in
      lift (re :+ im)

  fromRational x
    = lift (constant (fromRational x) :+ constant 0)

instance (Elt a, IsFloating a, RealFloat a) => Floating (Exp (Complex a)) where
  sqrt z
    = let
          x :+ y        = unlift z
          v'            = abs y / (u'*2)
          u'            = sqrt ((magnitude z + abs x) / 2)
          (u, v)        = unlift ( x A.<* 0 ? ( lift (v',u'), lift (u',v') ) )
      in
      x ==* 0 &&* y ==* 0 ?
        {- then -} ( 0
        {- else -} , lift (u :+ (y A.<* 0 ? (-v,v))) )

  pi            = lift (pi :+ constant 0)
  log z         = lift (log (magnitude z) :+ phase z)
  exp           = lift1 (exp :: Complex (Exp a) -> Complex (Exp a))
  sin           = lift1 (sin :: Complex (Exp a) -> Complex (Exp a))
  cos           = lift1 (cos :: Complex (Exp a) -> Complex (Exp a))
  tan           = lift1 (tan :: Complex (Exp a) -> Complex (Exp a))
  sinh          = lift1 (sinh :: Complex (Exp a) -> Complex (Exp a))
  cosh          = lift1 (cosh :: Complex (Exp a) -> Complex (Exp a))
  tanh          = lift1 (tanh :: Complex (Exp a) -> Complex (Exp a))
  asin          = lift1 (asin :: Complex (Exp a) -> Complex (Exp a))
  acos          = lift1 (acos :: Complex (Exp a) -> Complex (Exp a))
  atan          = lift1 (atan :: Complex (Exp a) -> Complex (Exp a))
  asinh         = lift1 (asinh :: Complex (Exp a) -> Complex (Exp a))
  acosh         = lift1 (acosh :: Complex (Exp a) -> Complex (Exp a))
  atanh         = lift1 (atanh :: Complex (Exp a) -> Complex (Exp a))

-- | The non-negative magnitude of a complex number
--
magnitude :: (Elt a, IsFloating a) => Exp (Complex a) -> Exp a
magnitude c =
  let r :+ i    = unlift c
  in sqrt (r*r + i*i)

-- | The phase of a complex number, in the range @(-'pi', 'pi']@. If the
-- magnitude is zero, then so is the phase.
--
phase :: (Elt a, IsFloating a) => Exp (Complex a) -> Exp a
phase c =
  let x :+ y    = unlift c
  in atan2 y x

-- | The function 'polar' takes a complex number and returns a (magnitude,
-- phase) pair in canonical form: the magnitude is non-negative, and the phase
-- in the range @(-'pi', 'pi']@; if the magnitude is zero, then so is the phase.
--
polar :: (Elt a, IsFloating a) => Exp (Complex a) -> Exp (a,a)
polar z =  lift (magnitude z, phase z)

-- | Form a complex number from polar components of magnitude and phase.
--
mkPolar :: (Elt a, IsFloating a) => Exp a -> Exp a -> Exp (Complex a)
mkPolar r theta  =  lift $ r * cos theta :+ r * sin theta

-- | @'cis' t@ is a complex value with magnitude @1@ and phase @t@ (modulo
-- @2*'pi'@).
--
cis :: (Elt a, IsFloating a) => Exp a -> Exp (Complex a)
cis theta = lift $ cos theta :+ sin theta

-- | Return the real part of a complex number
--
real :: Elt a => Exp (Complex a) -> Exp a
real c =
  let r :+ _    = unlift c
  in  r

-- | Return the imaginary part of a complex number
--
imag :: Elt a => Exp (Complex a) -> Exp a
imag c =
  let _ :+ i    = unlift c
  in  i

-- | Return the complex conjugate of a complex number, defined as
--
-- > conjugate(Z) = X - iY
--
conjugate :: (Elt a, IsNum a) => Exp (Complex a) -> Exp (Complex a)
conjugate z = lift $ real z :+ (- imag z)

max_stars_repo_name: fmma/accelerate

max_issues_repo_name: fmma/accelerate

max_forks_repo_name: fmma/accelerate


text:
{-# LANGUAGE FlexibleInstances, FlexibleContexts #-}
module Main where
import Test.Hspec
import Test.QuickCheck
import Numeric.LinearAlgebra
import qualified Data.Vector.Storable as V
import Test.Gen
import Test.Utils

main = hspec $ do
  describe "Corr Single" $ do
    it "implements corr2 correct.0" $ do
      forAll (pair (squared_real_matrices 3) (squared_real_matrices 4)) $
        \(m1, m2) -> ioProperty $ do
          r <- test_corr2 2 m1 m2
          return $ good_corr2 2 m1 m2 `eqShowWhenFail` r
    it "implements corr2 correct.1" $ do
      forAll (pair (squared_real_matrices 9) (squared_real_matrices 9)) $
        \(m1, m2) -> ioProperty $ do
          r <- test_corr2 4 m1 m2
          return $ good_corr2 4 m1 m2 `eq` r
    it "implements corr2 correct.2" $ do
      forAll (pair (squared_real_matrices 5) (squared_real_matrices 28)) $
        \(m1, m2) -> ioProperty $ do
          r <- test_corr2 2 m1 m2
          return $ eqShowWhenFail (good_corr2 2 m1 m2) r
    it "implements corr2 correct.3" $ do
      forAll (pair (choose (2,30)) (pair small_matrices small_matrices)) $
        \(p, (m1, m2)) -> ioProperty $ do
          r <- test_corr2 p m1 m2
          return $ good_corr2 p m1 m2 `eqShowWhenFail` r
  describe "Corr Many" $ do
    it "with 2 kernels" $ do
      forAll (pair (sequence $ replicate 2 $ squared_real_matrices 3) (squared_real_matrices 7)) $
        \(m1s, m2) -> ioProperty $ do
          rs <- test_corr2_arr 2 m1s m2
          return $ conjoin $ zipWith (\m r -> good_corr2 2 m m2 `eqShowWhenFail` r) m1s rs
    it "with 5 kernels" $ do
      forAll (pair (sequence $ replicate 4 $ squared_real_matrices 15) (squared_real_matrices 88)) $
        \(m1s, m2) -> ioProperty $ do
          rs <- test_corr2_arr 2 m1s m2
          ss <- return $ map (\m -> good_corr2 2 m m2) m1s
          return $ conjoin $ zipWith eq rs ss

eqShowWhenFail m1 m2 =
    whenFail (do let va = flatten m1
                 let vb = flatten m2
                 let err x 0 = x
                     err x y = abs ((x - y) / y)
                 let ev = (V.zipWith err va vb)
                     ei = V.maxIndex ev
                 putStrLn $ "Max error ration: " ++ show (ev V.! ei, va V.! ei, vb V.! ei)
                 putStrLn $ show m1
                 putStrLn $ show m2)
        (m1 `eq` m2)

ext: hs

max_stars_repo_name: pierric/neural-network

max_issues_repo_name: pierric/neural-network

max_forks_repo_name: pierric/neural-network


text:
{-# LANGUAGE BangPatterns #-}
module IllusoryContourShape where

import           Control.Monad               as M
import           Data.Array.Repa             as R
import           Data.Binary
import           Data.Complex
import           Data.List                   as L
import           Data.Vector.Generic         as VG
import           FokkerPlanck
import           Image.IO
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Array
import           Utils.Time
import qualified Data.Array.IArray as IA
import           FokkerPlanck.GreensFunction

main = do
  args@(gpuIDStr:numPointStr:deltaXStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:cutoffStr:deltaLogStr:taoStr:numTrailStr:maxTrailStr:phiFreqsStr:rhoFreqsStr:thetaFreqsStr:scaleFreqsStr:initScaleStr:shape2DStr:histFilePath:writeFlagStr:numIterationStr:suffix:stdStr:stdGStr:numThreadStr:_) <-
    getArgs
  let gpuID = read gpuIDStr :: [Int]
      numPoint = read numPointStr :: Int
      deltaX = read deltaXStr :: Double
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      phiFreq = read phiFreqsStr :: Double
      phiFreqs = [-phiFreq .. phiFreq]
      rhoFreq = read rhoFreqsStr :: Double
      rhoFreqs = [-rhoFreq .. rhoFreq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Double
      scaleFreqs = [-scaleFreq .. scaleFreq]
      initScale = read initScaleStr :: Double
      shape2D@(Points _ minDist _) = read shape2DStr :: Points Shape2D
      numThread = read numThreadStr :: Int
      folderPath =
        "output/test/IllusoryContourShape" </> (takeBaseName histFilePath) L.++
        "_" L.++
        cutoffStr L.++
        "_" L.++
        deltaXStr
      maxScale = read maxScaleStr :: Double
      cutoff = read cutoffStr :: Double
      halfLogPeriod = log maxScale
      deltaLog = read deltaLogStr :: Double
      writeFlag = read writeFlagStr :: Bool
      std = read stdStr :: Double
      stdG = read stdGStr :: Double
      numIteration = read numIterationStr :: Int
  print thetaFreqs
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  hist <-
    if flag
      then do
        printCurrentTime $
          "read Fourier coefficients data from " L.++ histFilePath
        decodeFile histFilePath
      -- else runMonteCarloFourierCoefficientsGPU
      --        gpuID
      --        numThread
      --        numTrail
      --        maxTrail
      --        thetaSigma
      --        scaleSigma
      --        maxScale
      --        tao
      --        phiFreqs
      --        rhoFreqs
      --        thetaFreqs
      --        scaleFreqs
      --        deltaLog
      --        initScale
      --        histFilePath
      --        (emptyHistogram
      --           [ L.length phiFreqs
      --           , L.length rhoFreqs
      --           , L.length thetaFreqs
      --           , L.length scaleFreqs
      --           ]
      --           0)
      else sampleCartesian
             folderPath
             histFilePath
             gpuID
             numPoint
             numPoint
             180
             deltaX
             deltaLog
             initScale
             thetaSigma
             tao
             phiFreqs
             rhoFreqs
             thetaFreqs
             scaleFreqs
  let !points =
        L.map (\(!x, !y) -> Point x y 0 1) . getShape2DIndexList . makeShape2D $
        shape2D
      -- !xs = L.map (\(a, b) -> (a * deltaX, b * deltaX)) . makeShape2D $ shape2D
      !xs = getShape2DIndexList . makeShape2D $ shape2D
      !coefficients =
        normalizeFreqArr std phiFreqs rhoFreqs . getNormalizedHistogramArr $
        hist
       -- = getNormalizedHistogramArr $ hist :: R.Array U DIM4 (Complex Double)
      !thetaRHarmonics =
        computeThetaRHarmonics
          numOrientation
          numScale
          thetaFreqs
          scaleFreqs
          halfLogPeriod
  print xs
  plan <-
    makePlan
      folderPath
      emptyPlan
      numPoint
      numPoint
      (L.length thetaFreqs)
      (L.length scaleFreqs)
  let !initSourceSparse =
        computeInitialDistributionPowerMethodSparse phiFreqs rhoFreqs points
      !harmonicsArraySparse =
        computeHarmonicsArraySparse
          numPoint
          deltaX
          numPoint
          deltaX
          phiFreqs
          rhoFreqs
          thetaFreqs
          scaleFreqs
          halfLogPeriod
          cutoff
  -- let !initSource =
  --       computeInitialDistributionPowerMethod'
  --         numPoint
  --         numPoint
  --         phiFreqs
  --         rhoFreqs
  --         -- thetaFreqs
  --         -- scaleFreqs
  --         points
  --     !bias = computeBias numPoint numPoint points
  gaussian <- gaussianFilter2D plan numPoint stdG
  -- harmonicsArray <-
  --   dftHarmonicsArrayG
  --     plan
  --     numPoint
  --     deltaX
  --     numPoint
  --     deltaX
  --     phiFreqs
  --     rhoFreqs
  --     thetaFreqs
  --     scaleFreqs
  --     halfLogPeriod
  --     cutoff
  --     gaussian
  -- print . IA.indices $ harmonicsArraySparse
  -- completion <-
  computeContourSparse
    plan
    folderPath
    coefficients
    harmonicsArraySparse
    -- thetaRHarmonics
    (fromListUnboxed (Z :. L.length phiFreqs) phiFreqs)
    (fromListUnboxed (Z :. L.length rhoFreqs) rhoFreqs)
    cutoff
    -- gaussian
                      -- std
    xs
    numIteration
    suffix
    initSourceSparse
  -- completion <-
  --   computeContour'
  --    plan
  --    folderPath
  --    writeFlag
  --    coefficients
  --    harmonicsArray
  --    bias
  --    -- gaussian
  --    numIteration
  --    suffix
  --    thetaRHarmonics
  --    initSource
  -- plotDFTArrayThetaR
  --   (folderPath </> (printf "Completion_%s.png" suffix))
  --   numPoint
  --   numPoint
  --   thetaRHarmonics
  --   completion
  -- plotDFTArrayThetaRMag
  --   (folderPath </> (printf "CompletionMax_%s.png" suffix))
  --   numPoint
  --   numPoint
  --   thetaRHarmonics
  --   completion
  -- computeContourSparse'''
  --   plan
  --   folderPath
  --   coefficients
  --   harmonicsArraySparse
  --   thetaRHarmonics
  --   (fromListUnboxed (Z :. L.length phiFreqs) phiFreqs)
  --   (fromListUnboxed (Z :. L.length rhoFreqs) rhoFreqs)
  --   cutoff
  --   xs
  --   numIteration
  --   (suffix L.++ "Test")
  --   initSourceSparse

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 LinearRegressionExercise where

import qualified Graphics.Gnuplot.Advanced as GP

import qualified Graphics.Gnuplot.Graph as Graph

import qualified Graphics.Gnuplot.Plot.ThreeDimensional as Plot3D
import qualified Graphics.Gnuplot.Graph.ThreeDimensional as Graph3D

import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D
import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D
import qualified Graphics.Gnuplot.Terminal.Default as DefaultTerm
import Graphics.Gnuplot.Plot.TwoDimensional (linearScale, )

import LinearRegression
import Numeric.LinearAlgebra
import Prelude hiding ((<>))
import System.Directory
import Util

-- replace this with the path to your data.
exerciseOneDataOne :: IO String
exerciseOneDataOne = do
  home <- getHomeDirectory
  return (home ++ "/code/haskell/CourseraMachineLearning/app/ex1data1.txt")

exerciseOneDataTwo :: IO String
exerciseOneDataTwo = do
  home <- getHomeDirectory
  return (home ++ "/code/haskell/CourseraMachineLearning/app/ex1data2.txt")

exerciseOne = do
  -- Part 1: Basic function
  putStrLn "Running warmUpExercise ..."
  putStrLn "5x5 Identity Matrix: "
  let a = ident 5 :: Matrix Double
  print a

  -- Part 2: Plotting
  putStrLn "Plotting Data ..."
  dataOne <- exerciseOneDataOne
  [x', y'] <- fmap toColumns (loadMatrix dataOne)
  GP.plotSync DefaultTerm.cons (plotPoints (toList x') (toList y'))

  -- Part 3: Cost and Gradient Descent
  let m = size y'
  let x = ones m 1 ||| reshape 1 x'
  let y = asColumn y'
  let xPoints = toList x'
  let yPoints = toList y'
  let theta = zeros 2 1
  let iterations = 1500
  let alpha = 0.01
  putStrLn "Testing the cost function ..."
  let j = computeCost x y theta
  putStrLn ("With theta = [0, 0]\n Cost computed = " ++ show j)
  putStrLn "Expected cost value (approx) 32.07"
  let k = computeCost x y (asColumn $ 2 |> [-1, 2])
  putStrLn ("With theta = [-1, 2]\n Cost computed = " ++ show k)
  putStrLn "Expected cost value (approx) 54.24"

  putStrLn "Running Gradient Descent ..."
  let (costHist, thetaOpt) = gradientDescent  x y theta alpha iterations
  putStrLn ("Theta found by gradient descent: " ++ show thetaOpt)
  putStrLn "Expected theta values (approx)\n"
  putStrLn " -3.6303\n  1.1664"
  let yOptPoints = (toList (flatten $ x <> thetaOpt))
  putStrLn " -3.6303\n  1.1664"
  let yOptPoints = (toList (flatten $ x <> thetaOpt))
  GP.plotSync DefaultTerm.cons $ plotLine xPoints yOptPoints `mappend` plotPoints xPoints yPoints

  -- Predicting values
  let predict1 = row [1, 3.5] <> thetaOpt
  putStrLn "For population = 35,000 we predict a profit of "
  print ((predict1 * 1000) `atIndex` (0,0))

  let predict2 = row [1, 20] <> thetaOpt
  putStrLn "For population = 70,000 we predict a profit of "
  print ((predict2 * 1000) `atIndex` (0,0))

  -- Visualising J(theta_0, theta_1)
  let theta0 = toList (linspace 100 (-10 :: Double, 10 :: Double))
  let theta1 = toList (linspace 100 (-1  :: Double, 4  :: Double))
  let jVals = [computeCost x y (col [i, j]) | i <- theta0, j <- theta1]
  -- TODO(mchav): Draw contour plot as well.
  GP.plotSync DefaultTerm.cons $ Plot3D.surface theta0 theta1 (\i j -> computeCost x y (col [i, j]))

exerciseOneOptional = do
  dataTwo <- exerciseOneDataTwo
  columns <- fmap toColumns (loadMatrix dataTwo)
  let y' = last columns
  let y = asColumn y'
  let x' = fromColumns (init columns)
  let m = size y'

  let (xNorm, mu, sigma) = featureNormalize x'
  let x = ones m 1 ||| xNorm

  let theta = zeros 3 1
  let iterations = 50 :: Int
  let alpha = 1

  -- TODO(mchav): Include cost history in gradient descent function and then plot it.
  let (costHist, thetaOpt) = gradientDescent x y theta alpha iterations
  putStrLn ("Theta found by gradient descent: " ++ show thetaOpt)
  -- Plot the graph ingoring the first (zeroth) iteration.
  GP.plotSync DefaultTerm.cons $ plotLine [1..(fromIntegral iterations)] ((tail . reverse) costHist)
  putStrLn "Predicted price of a 1650 sq-ft, 3 br house (using gradient descent)"
  let xTest = ones 1 1 ||| ((row [1650, 3] - mu) / sigma)
  print (xTest <> thetaOpt)
  let x'' = ones m 1 ||| x'
  let normTheta = normalEquation x'' y
  print normTheta
  let xTest' = ones 1 1 ||| row [1650, 3]
  print (xTest' <> normTheta)

max_stars_repo_name: mchav/CourseraMachineLearning

max_issues_repo_name: mchav/CourseraMachineLearning

max_forks_repo_name: mchav/CourseraMachineLearning


text:
{-# LANGUAGE BangPatterns                             #-}
{-# LANGUAGE DeriveFunctor                            #-}
{-# LANGUAGE DeriveGeneric                            #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE KindSignatures                           #-}
{-# LANGUAGE LambdaCase                               #-}
{-# LANGUAGE NoStarIsType                             #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE StandaloneDeriving                       #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilyDependencies                   #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# LANGUAGE ViewPatterns                             #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module Numeric.Wavelet.Continuous (
  -- * Continuous Wavelet Transform
    CWD(..), mapCWD
  , CWDLine(..), mapCWDLine
  , LNorm(..), CWDOpts(..), defaultCWDO
  , cwd
  , cwdReal
  -- * Wavelets
  , AWavelet(..), mapAWavelet
  , morlet, morletFunc
  , meyer, meyerFunc
  , fbsp, fbspFunc
  ) where

import           Data.Complex
import           Data.Finite
import           Data.Foldable
import           Data.Maybe
import           Data.Ord
import           Data.Proxy
import           Data.Type.Equality
import           Data.Vector.Generic.Sized    (Vector)
import           GHC.Generics                 (Generic)
import           GHC.TypeLits.Compare
import           GHC.TypeNats
import           Math.FFT.Base
import           Numeric.Wavelet.Internal.FFT
import qualified Data.Vector.Generic          as UVG
import qualified Data.Vector.Generic.Sized    as VG
import qualified Data.Vector.Sized            as V

newtype CWD v n m a b = CWD { cwdLines :: V.Vector m (CWDLine v n a b) }
  deriving (Show, Functor)

data CWDLine v n a b = CWDLine
    { cwdlData  :: Vector v n b
    , cwdlScale :: Finite (n `Div` 2 + 1) -- ^ Scale factor, in number of ticks.
    , cwdlFreq  :: a                      -- ^ The frequency associated with this scale, in inverse tick
    , cwdlCoI   :: Finite (n `Div` 2 + 1) -- ^ How many items are /outside/ of the Cone of Influence, on each side.
    }
  deriving (Show, Functor)

mapCWDLine
    :: (UVG.Vector v b, UVG.Vector v c)
    => (b -> c)
    -> CWDLine v n a b
    -> CWDLine v n a c
mapCWDLine f (CWDLine d s q c) = CWDLine (VG.map f d) s q c

mapCWD
    :: (UVG.Vector v b, UVG.Vector v c)
    => (b -> c)
    -> CWD v n m a b
    -> CWD v n m a c
mapCWD f (CWD l) = CWD (mapCWDLine f <$> l)

data LNorm = L1 | L2
  deriving (Show, Eq, Ord, Enum, Bounded, Generic)

data CWDOpts n = CWDO
    { cwdoNorm     :: LNorm                     -- ^ wavelet normalization
    , cwdoMinScale :: Finite (n `Div` 2 + 1)    -- ^ min scale (period)
    , cwdoMaxScale :: Finite (n `Div` 2 + 1)    -- ^ max scale (period)
    }
  deriving (Show, Eq, Ord, Generic)

defaultCWDO :: KnownNat n => CWDOpts n
defaultCWDO = CWDO
    { cwdoNorm     = L2
    , cwdoMinScale = minBound
    , cwdoMaxScale = maxBound
    }

-- | 'cwd' for complex-valued analytic wavelets.
cwd :: forall v n m a b.
     ( UVG.Vector v (Complex b)
     , KnownNat n
     , KnownNat m
     , FFTWReal b
     , 1 <= n
     , RealFloat a
     )
    => AWavelet v a (Complex b)
    -> CWDOpts n
    -> Vector v n (Complex b)
    -> CWD v n m a (Complex b)
cwd AW{..} CWDO{..} xs = CWD . VG.generate $ \i ->
      let s   = scaleOf i
          dt  = 1/s
      in  case awVector dt of
            VG.SomeSized (wv :: Vector v q (Complex b))
              | Just Refl <- isLE (Proxy @1) (Proxy @q)
              , Just Refl <- isLE (Proxy @((q-1)`Div`2)) (Proxy @(q-1))
              -> let ys :: Vector v (n + q - 1) (Complex b)
                     ys  = (* (realToFrac (normie dt) :+ 0)) `VG.map` convolve xs wv
                     coi = fromMaybe maxBound . packFinite . round @a @Integer $ sqrt 2 * s
                     s'  = fromMaybe maxBound . packFinite . round @a @Integer $ s
                     ys' :: Vector v n (Complex b)
                     ys' = VG.slice @_ @((q - 1)`Div`2) @n @((q-1)-((q-1)`Div`2)) Proxy ys
                 in  CWDLine ys' s' (awFreq / s) coi
            _ -> error "Bad scale: wavelet vector is empty?"
  where
    n = natVal (Proxy @n)
    m = natVal (Proxy @m)
    normie :: a -> a
    normie = case cwdoNorm of
      L1 -> sqrt
      L2 -> id
    minScale = fromIntegral cwdoMinScale `max` 1
    maxScale = (fromIntegral cwdoMaxScale `min` (fromIntegral n / (2 * sqrt 2)))
         `max` (minScale + 1)
    scaleStep = (log maxScale - log minScale) / (fromIntegral m - 1)
    scaleOf :: Finite m -> a
    scaleOf i = exp $ log minScale + fromIntegral i * scaleStep

-- | 'cwd' for real-valued analytic wavelets.
cwdReal
    :: forall v n m a b.
     ( UVG.Vector v b
     , UVG.Vector v (Complex b)
     , KnownNat n
     , KnownNat m
     , FFTWReal b
     , 1 <= n
     , RealFloat a
     )
    => AWavelet v a b
    -> CWDOpts n
    -> Vector v n b
    -> CWD v n m a b
cwdReal aw cwdo = mapCWD realPart
                . cwd (mapAWavelet (:+ 0) aw) cwdo
                . VG.map (:+ 0)

-- | Analytical Wavelet
data AWavelet v a b = AW
    { awVector :: a -> v b    -- ^ generate a vector within awRange with a given dt
    , awFreq   :: a           -- ^ Dominant frequency component
    , awRange  :: a           -- ^ range away from zero outside of which wavelet can be considered negligible
    }
  deriving Functor

mapAWavelet
    :: (UVG.Vector v b, UVG.Vector v c)
    => (b -> c)
    -> AWavelet v a b
    -> AWavelet v a c
mapAWavelet f (AW v q r) = AW (UVG.map f . v) q r

morlet
    :: (UVG.Vector v (Complex a), RealFloat a)
    => a
    -> AWavelet v a (Complex a)
morlet σ = AW{..}
  where
    awRange  = 4
    (!awFreq, mf) = morletFunc_ σ
    awVector = renderFunc awRange mf

morletFunc :: RealFloat a => a -> a -> Complex a
morletFunc = snd . morletFunc_

morletFunc_ :: RealFloat a => a -> (a, a -> Complex a)
morletFunc_ σ = (q, f)
  where
    f t = (c * exp(-t*t/2) :+ 0) * (exp (0 :+ (σ * t)) - (exp (-σ2/2) :+ 0))
    !c   = pi ** (-1/4) * (1 + exp (-σ2) - 2 * exp (-3/4*σ2)) ** (-1/2)
    !σ2  = σ * σ
    !q   = converge 20 iter σ / (2 * pi)
    iter w = σ / (1 - exp (-σ * w))

meyer
    :: (UVG.Vector v a, RealFloat a)
    => AWavelet v a a
meyer = AW{..}
  where
   awRange  = 6
   awVector = renderFunc awRange meyerFunc
   awFreq   = 4 * pi / 3

meyerFunc :: RealFloat a => a -> a
meyerFunc t
    | isNaN ψ || isInfinite ψ = 0
    | otherwise               = ψ
  where
    t' = t - 0.5
    t'3 = t'**3
    sinTerm = sin(4*pi/3*t') / pi
    ψ1 = (4/3/pi*t'*cos(2*pi/3*t') - sinTerm)
       / (t' - 16/9 * t'3)
    ψ2 = (8/3/pi*t'*cos(8*pi/3*t') + sinTerm)
       / (t' - 64/9 * t'3)
    ψ = ψ1 + ψ2

fbsp
    :: (UVG.Vector v (Complex a), FFTWReal a)
    => Int      -- ^ m, >= 1
    -> a        -- ^ f_b, bandwidth
    -> a        -- ^ f_c, wavelet center frequency
    -> AWavelet v a (Complex a)
fbsp m fb fc = AW{..}
  where
    awRange  = 4/fb
    awVector = renderFunc awRange (fbspFunc m fb fc)
    awFreq   = autoDeriveFreq awRange awVector

autoDeriveFreq
    :: (UVG.Vector v (Complex a), FFTWReal a)
    => a
    -> (a -> v (Complex a))
    -> a
autoDeriveFreq r fv = case fv 0.001 of
    VG.SomeSized v ->
      let vv    = zip [1..] . map magnitude . VG.toList $ fft v
          (i,_) = maximumBy (comparing snd) vv
      in  fromInteger i / (r * 2)
    _ -> error "bad vector"

fbspFunc :: RealFloat a => Int -> a -> a -> a -> Complex a
fbspFunc m fb fc t =
    ((sqrt fb * sinc (fb * t / fromIntegral m))^m :+ 0) * exp (0 :+ (2 * pi * fc * t))
  where
    sinc x = sin x / x

-- | Render the effective range of a wavelet (based on 'awRange'), centered
-- around zero.  Takes a timestep.
renderFunc
    :: (UVG.Vector v b, RealFrac a)
    => a              -- ^ range about zero
    -> (a -> b)       -- ^ func
    -> a              -- ^ dt
    -> v b
renderFunc r f dt = UVG.generate (round n) $ \i ->
    f (fromIntegral i * dt - r)
  where
    n  = r * 2 / dt

converge
    :: (Fractional a, Ord a)
    => Int      -- ^ maximum iterations
    -> (a -> a) -- ^ function to find the fixed point convergence
    -> a        -- ^ starting value
    -> a
converge n f = go 0
  where
    go !i !x
        | i >= n                = 0
        | abs (x - y) < 0.00001 = x
        | otherwise             = go (i + 1) y
      where
        !y = f x

max_stars_repo_name: mstksg/wavelets

max_issues_repo_name: mstksg/pure-wavelets

max_forks_repo_name: mstksg/pure-wavelets


text:
-- | Helpers for testing
module Tests.Helpers (
    -- * helpers
    T(..)
  , typeName
  , eq
  , eqC
    -- * Generic QC tests
  , monotonicallyIncreases
    -- * HUnit helpers
  , testAssertion
  , testEquality
  ) where

import Data.Complex
import Data.Typeable

import qualified Test.HUnit      as HU
import Test.Framework
import Test.Framework.Providers.HUnit

import Numeric.MathFunctions.Comparison

----------------------------------------------------------------
-- Helpers
----------------------------------------------------------------

-- | Phantom typed value used to select right instance in QC tests
data T a = T

-- | String representation of type name
typeName :: Typeable a => T a -> String
typeName = show . typeOf . typeParam
  where
    typeParam :: T a -> a
    typeParam _ = undefined

-- | Approximate equality for 'Double'. Doesn't work well for numbers
--   which are almost zero.
eq :: Double                    -- ^ Relative error
   -> Double -> Double -> Bool
eq = eqRelErr

-- | Approximate equality for 'Complex Double'
eqC :: Double                   -- ^ Relative error
    -> Complex Double
    -> Complex Double
    -> Bool
eqC eps a@(ar :+ ai) b@(br :+ bi)
  | a == 0 && b == 0 = True
  | otherwise        = abs (ar - br) <= eps * d
                    && abs (ai - bi) <= eps * d
  where
    d = max (realPart $ abs a) (realPart $ abs b)

----------------------------------------------------------------
-- Generic QC
----------------------------------------------------------------

-- Check that function is nondecreasing
monotonicallyIncreases :: (Ord a, Ord b) => (a -> b) -> a -> a -> Bool
monotonicallyIncreases f x1 x2 = f (min x1 x2) <= f (max x1 x2)

----------------------------------------------------------------
-- HUnit helpers
----------------------------------------------------------------

testAssertion :: String -> Bool -> Test
testAssertion str cont = testCase str $ HU.assertBool str cont

testEquality :: (Show a, Eq a) => String -> a -> a -> Test
testEquality msg a b = testCase msg $ HU.assertEqual msg a b

max_stars_repo_name: Javran/math-functions

max_issues_repo_name: Javran/math-functions

max_forks_repo_name: Javran/math-functions


text:
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeOperators       #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE TypeApplications    #-}

{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Test.Grenade.Batch where

import           Grenade.Types
import           Grenade.Core.Shape
import           Grenade.Core.Network
import           Grenade.Core.Layer
import           Grenade.Layers.FullyConnected
import           Grenade.Layers.Convolution
import           Grenade.Utils.ListStore

import           Numeric.LinearAlgebra hiding (uniformSample, konst, (===))
import           Numeric.LinearAlgebra.Static as H hiding ((===))
import           Numeric.LinearAlgebra.Data as D hiding ((===))

import           Hedgehog
import           GHC.TypeLits

type FFNetwork = Network '[ FullyConnected 3 5, FullyConnected 5 4 ] '[ 'D1 3, 'D1 5, 'D1 4 ]

prop_networkBatchFeedforward = property $ do
  let bias  :: H.R 5              = H.fromList [1..5]
      bias' :: H.R 4              = H.fromList [1..4]
      acts  :: H.L 5 3            = H.fromList [1..15]
      acts' :: H.L 4 5            = H.fromList [1..20]
      fc    :: FullyConnected 3 5 = FullyConnected (FullyConnected' bias  acts)  mkListStore
      fc'   :: FullyConnected 5 4 = FullyConnected (FullyConnected' bias' acts') mkListStore
      ins   :: [S ('D1 3)]        = [S1D (H.fromList [1, 2, 3]), S1D (H.fromList [4, 5, 6])]
      net   :: FFNetwork          = fc :~> fc' :~> NNil
      (_, outs :: [S ('D1 4)]) = runBatchForwards net ins
      outs' = map (\(S1D v) -> (D.toList . H.extract) v) outs
  outs' === [[986, 2312, 3638, 4964], [2336, 5462, 8588, 11714]]

extractFCGrads :: Gradients '[FullyConnected 3 5, FullyConnected 5 4] -> [(Vector RealNum, Matrix RealNum)]
extractFCGrads ((FullyConnected' wB wN) :/> ((FullyConnected' wB' wN') :/> GNil))
  = [(H.extract wB, H.extract wN), (H.extract wB', H.extract wN')]

prop_networkBackpropCalculatesGradients = property $ do
  let bias  :: H.R 5              = H.fromList [1..5]
      bias' :: H.R 4              = H.fromList [1..4]
      acts  :: H.L 5 3            = H.fromList [1..15]
      acts' :: H.L 4 5            = H.fromList [1..20]
      fc    :: FullyConnected 3 5 = FullyConnected (FullyConnected' bias acts) mkListStore
      fc'   :: FullyConnected 5 4 = FullyConnected (FullyConnected' bias' acts') mkListStore
      ins   :: [S ('D1 3)]        = [S1D (H.fromList [1, 2, 3]), S1D (H.fromList [4, 5, 6])]
      net   :: FFNetwork          = fc :~> fc' :~> NNil
      (tapes, outs :: [S ('D1 4)]) = runBatchForwards net ins
      (grads, vs   :: [S ('D1 3)]) = runBatchBackwards net tapes outs
      (grad,  v    :: S ('D1 3))   = runBackwards net (tapes!!0) (outs!!0)
      (grad', v'   :: S ('D1 3))   = runBackwards net (tapes!!1) (outs!!1)
      grads'                       = map extractFCGrads grads
      grads''                      = map extractFCGrads [grad, grad']
      vs'                          = map (\(S1D vec) -> (D.toList . H.extract) vec) vs
      vs''                         = map (\(S1D vec) -> (D.toList . H.extract) vec) [v, v']

  grads' === grads''
  vs'    === vs''

prop_networkAveragesGradients = property $ do
  let bias  :: H.R 5                 = H.fromList [1..5]
      bias' :: H.R 4                 = H.fromList [1..4]
      acts  :: H.L 5 3               = H.fromList [1..15]
      acts' :: H.L 4 5               = H.fromList [1..20]
      fc    :: FullyConnected 3 5    = FullyConnected (FullyConnected' bias acts) mkListStore
      fc'   :: FullyConnected 5 4    = FullyConnected (FullyConnected' bias' acts') mkListStore
      ins   :: [S ('D1 3)]           = [S1D (H.fromList [1, 2, 3]), S1D (H.fromList [4, 5, 6])]
      net   :: FFNetwork             = fc :~> fc' :~> NNil
      (tapes, outs :: [S ('D1 4)])   = runBatchForwards net ins
      (grads, _    :: [S ('D1 3)])   = runBatchBackwards net tapes outs
      (grad,  _    :: S ('D1 3))     = runBackwards net (tapes!!0) (outs!!0)
      (grad', _    :: S ('D1 3))     = runBackwards net (tapes!!1) (outs!!1)
      rgrad                          = extractFCGrads (reduceGradient @FFNetwork grads)
      [(wB, wN), (wB', wN')]         = extractFCGrads grad
      [(wB'', wN''), (wB''', wN''')] = extractFCGrads grad'
      rgrad'                         = (0.5 * (wB + wB''), 0.5 * (wN + wN''))
      rgrad''                        = (0.5 * (wB' + wB'''), 0.5 * (wN' + wN'''))

  rgrad === [rgrad', rgrad'']

prop_convolutionCalculatesOutputOfBatches = property $ do
  let weights :: H.L 25 1 = H.fromList [1..25]
      convLayer :: Convolution 'WithoutBias 'NoPadding 1 1 5 5 2 2 = Convolution weights mkListStore
      ins :: [S ('D2 11 11)] = [S2D (H.fromList [1..121]), S2D (H.fromList [2..122])]
      (_, outs :: [S ('D2 4 4)]) = runBatchForwards convLayer ins
      outs' = map (\(S2D v) -> (concat . D.toLists . H.extract) v) outs
  (take 9 (concat outs')) === [10925, 11575, 12225, 12875, 18075, 18725, 19375, 20025, 25225]

unwrapGradConv :: ( KnownNat c
              , KnownNat f
              , KnownNat kR
              , KnownNat kC
              , KnownNat sR
              , KnownNat sC
              , KnownNat kF
              , kF ~ (kR * kC * c)) => Convolution' 'WithoutBias 'NoPadding c f kR kC sR sC -> Matrix RealNum
unwrapGradConv (Convolution' mat) = H.extract mat

prop_convolutionBackprop = property $ do
  let weights :: H.L 25 1 = H.fromList [1..25]
      convLayer :: Convolution 'WithoutBias 'NoPadding 1 1 5 5 2 2 = Convolution weights mkListStore
      ins :: [S ('D2 11 11)] = [S2D (H.fromList [1..121]), S2D (H.fromList [2..122])]
      (tapes, outs :: [S ('D2 4 4)]) = runBatchForwards convLayer ins
      (grads, vs   :: [S ('D2 11 11)]) = runBatchBackwards convLayer tapes outs
      (grad,  v    :: S ('D2 11 11))   = runBackwards convLayer (tapes!!0) (outs!!0)
      (grad', v'   :: S ('D2 11 11))   = runBackwards convLayer (tapes!!1) (outs!!1)
      grads'                           = map unwrapGradConv grads
      grads''                          = map unwrapGradConv [grad, grad']
      vs'                              = map (\(S2D u) -> (D.toLists . H.extract) u) vs
      vs''                             = map (\(S2D u) -> (D.toLists . H.extract) u) [v, v']

  grads' === grads''
  vs'    === vs''

prop_convolutionAveragesGradients = property $ do
  let weights :: H.L 25 1                  = H.fromList (concat (replicate 5 [1, 2, 3, 4, 5]))
      convLayer :: Convolution 'WithoutBias 'NoPadding 1 1 5 5 2 2 = Convolution weights mkListStore
      ins :: [S ('D2 11 11)]               = [S2D (H.fromList [1..121]), S2D (H.fromList [2..122])]
      (tapes, outs :: [S ('D2 4 4)])       = runBatchForwards convLayer ins
      (grads, _    :: [S ('D2 11 11)])     = runBatchBackwards convLayer tapes outs
      (grad,  _    :: S ('D2 11 11))       = runBackwards convLayer (tapes!!0) (outs!!0)
      (grad', _    :: S ('D2 11 11))       = runBackwards convLayer (tapes!!1) (outs!!1)
      rgrad                                = unwrapGradConv (reduceGradient @(Convolution 'WithoutBias 'NoPadding 1 1 5 5 2 2) grads)
      w                                    = unwrapGradConv grad
      w'                                   = unwrapGradConv grad'
      rgrad'                       = 0.5 * (w + w')

  rgrad === rgrad'

prop_fullyConnectedCalculatesOutputOfBatches = property $ do
  let bias :: H.R 5 = H.fromList [1..5]
      acts :: H.L 5 3 = H.fromList [1..15]
      fc :: FullyConnected 3 5 = FullyConnected (FullyConnected' bias acts) mkListStore
      ins :: [S ('D1 3)] = [S1D (H.fromList [1, 2, 3]), S1D (H.fromList [4, 5, 6])]
      (_, outs :: [S ('D1 5)]) = runBatchForwards fc ins
      outs' = map (\(S1D v) -> (D.toList . H.extract) v) outs
  outs' === [[15, 34, 53, 72, 91], [33, 79, 125, 171, 217]]

prop_fullyConnectedBackprop = property $ do
  let bias :: H.R 5 = H.fromList [1..5]
      acts :: H.L 5 3 = H.fromList [1..15]
      fc :: FullyConnected 3 5 = FullyConnected (FullyConnected' bias acts) mkListStore
      ins :: [S ('D1 3)] = [S1D (H.fromList [1, 2, 3]), S1D (H.fromList [4, 5, 6])]
      (tapes, outs :: [S ('D1 5)]) = runBatchForwards fc ins
      (grads, vs   :: [S ('D1 3)]) = runBatchBackwards fc tapes outs
      (grad,  v    :: S ('D1 3))   = runBackwards fc (tapes!!0) (outs!!0)
      (grad', v'   :: S ('D1 3))   = runBackwards fc (tapes!!1) (outs!!1)
      unwrapGrad                   = \(FullyConnected' wB wN) -> (H.extract wB, H.extract wN)
      grads'                       = map unwrapGrad grads
      grads''                      = map unwrapGrad [grad, grad']
      vs'                          = map (\(S1D vec) -> (D.toList . H.extract) vec) vs
      vs''                         = map (\(S1D vec) -> (D.toList . H.extract) vec) [v, v']

  grads' === grads''
  vs'    === vs''

prop_fullyConnectedAveragesGradients = property $ do
  let bias :: H.R 5 = H.fromList [1..5]
      acts :: H.L 5 3 = H.fromList [1..15]
      fc :: FullyConnected 3 5 = FullyConnected (FullyConnected' bias acts) mkListStore
      ins :: [S ('D1 3)] = [S1D (H.fromList [1, 2, 3]), S1D (H.fromList [4, 5, 6])]
      (tapes, outs :: [S ('D1 5)]) = runBatchForwards fc ins
      (grads, _    :: [S ('D1 3)]) = runBatchBackwards fc tapes outs
      (grad,  _    :: S ('D1 3))   = runBackwards fc (tapes!!0) (outs!!0)
      (grad', _    :: S ('D1 3))   = runBackwards fc (tapes!!1) (outs!!1)
      f                            = \(FullyConnected' bs ns) -> (H.extract bs, H.extract ns)
      rgrad                        = f (reduceGradient @(FullyConnected 3 5) grads)
      (wB, wN)                     = f grad
      (wB', wN')                   = f grad'
      rgrad'                       = (0.5 * (wB + wB'), 0.5 * (wN + wN'))

  rgrad === rgrad'

tests :: IO Bool
tests = checkParallel $$(discover)

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-|
Module: MachineLearning.Regression
Description: Regression
Copyright: (c) Alexander Ignatyev, 2016-2018.
License: BSD-3
Stability: experimental
Portability: POSIX
-}

module MachineLearning.Regression
(
  Model.Model(..)
  , LeastSquares.LeastSquaresModel(..)
  , Optimization.MinimizeMethod(..)
  , Optimization.minimize
  , normalEquation
  , normalEquation_p
  , Regularization(..)
)

where

import Prelude hiding ((<>))
import MachineLearning.Types (Vector, Matrix)
import MachineLearning.Optimization as Optimization
import MachineLearning.Model as Model
import MachineLearning.LeastSquaresModel as LeastSquares
import MachineLearning.Regularization (Regularization(..))

import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra ((<>), (#>))

-- | Normal equation using inverse, does not require feature normalization
-- It takes X and y, returns theta.
normalEquation :: Matrix -> Vector -> Vector
normalEquation x y =
  let trX = LA.tr x
  in (LA.inv (trX <> x) <> trX) #> y

-- | Normal equation using pseudo inverse, requires feature normalization
-- It takes X and y, returns theta.
normalEquation_p :: Matrix -> Vector -> Vector
normalEquation_p x y =
  let trX = LA.tr x
  in (LA.pinv (trX <> x) <> trX) #> y

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE BangPatterns      #-}
{-# LANGUAGE OverloadedStrings #-}
module ML.Gym.Gym
    ( Gym (..)
    , GymResult (..)
    , GymObservation
    , initGym
    , resetGym
    , stepGymRandom
    , stepGymRenderRandom
    , stepGym
    , stepGymRender
    , setMaxEpisodeSteps
    , getMaxEpisodeSteps
    , getElapsedSteps
    ) where

import           ML.Gym.Data
import           ML.Gym.DType
import           ML.Gym.Range
import           ML.Gym.Shape
import           ML.Gym.Space
import           ML.Gym.Util

import qualified Control.Exception        as E
import           Control.Monad            (void, when)
import qualified CPython                  as Py
import qualified CPython.Constants        as Py
import qualified CPython.Protocols.Object as Py
import qualified CPython.Types            as Py
import qualified CPython.Types.Exception  as Py
import qualified CPython.Types.Module     as Py
import qualified Data.ByteString          as BS
import           Data.Complex             (Complex)
import           Data.Maybe               (fromMaybe)
import qualified Data.Text                as T
import           System.Environment       (getArgs, getProgName)
import           System.Exit
import           System.IO
import           System.Random            (randomRIO)

type GymEnv = Py.SomeObject
type GymModule = Py.Module
type GymObservation = GymData
type GymActionSpace = GymSpace
type GymObservationSpace = GymSpace

data Gym = Gym
  { name             :: !T.Text
  , gym              :: !GymModule
  , env              :: !GymEnv
  , actionSpace      :: !GymActionSpace
  , observationSpace :: !GymObservationSpace
  }

instance Show Gym where
  show (Gym n _ _ acts obs) = T.unpack n ++  "; Action Space: " ++ show acts ++ "; Observation Space: " ++ show obs

data GymResult = GymResult
  { observation :: !GymObservation
  , reward      :: !Double
  , episodeDone :: !Bool
  } deriving (Show, Eq)

-------------------- Init --------------------

getActionSpace :: GymEnv -> IO GymActionSpace
getActionSpace gEnv = do
  actSpace <- python $ Py.getAttribute gEnv =<< Py.toUnicode "action_space"
  toGymSpace actSpace

getObservationSpace :: GymEnv -> IO GymObservationSpace
getObservationSpace gEnv = do
  obsSpace <- python $ Py.getAttribute gEnv =<< Py.toUnicode "observation_space"
  toGymSpace obsSpace

-- | Initializes the gym environment. It also resets the environment, s.t. `stepGym` can be called immediately.
initGym :: T.Text -> IO (GymObservation, Gym)
initGym envName = do
  name <- getProgName
  args <- getArgs
  -- putStrLn $ "Called: " ++ name  ++ " " ++ unwords args
  Py.initialize
  res <- Py.isInitialized
  version <- Py.getVersion
  -- putStrLn $ "Initialized python: " ++ show res ++ "; Version: " ++ show version
  Py.setArgv (T.pack name <> "-" <> envName) (map T.pack args)
  -- putStrLn "Arguments set"
  gGym <- python $ Py.importModule "gym"
  -- putStrLn "Module imported"
  pyName <- Py.toUnicode envName
  -- putStrLn "Calling make"
  gEnv <- python $ Py.callMethodArgs gGym "make" [Py.toObject pyName]
  -- putStrLn "Createing action and observation space make"
  gActSpace <- getActionSpace gEnv
  gObsSpace <- getObservationSpace gEnv
  -- putStrLn "Resetting environment"
  gObservation <- python $ Py.callMethodArgs gEnv "reset" []
  obs <- fromMaybe (error "could not convert observation to GymData") <$> python (getGymData gObsSpace gObservation)
  -- putStrLn "Done"
  return (obs, Gym envName gGym gEnv gActSpace gObsSpace)

-------------------- Functions --------------------

-- | Set the maximum number of steps before an episode is finished.
setMaxEpisodeSteps :: Gym -> Integer -> IO ()
setMaxEpisodeSteps gym val = do
  pyVal <- Py.toInteger val
  python $ Py.toUnicode "_max_episode_steps" >>= \attr -> Py.setAttribute (env gym) attr (Py.toObject pyVal)

-- | Get the maximum number of steps before an episode is finished.
getMaxEpisodeSteps :: Gym -> IO (Maybe Integer)
getMaxEpisodeSteps gym = python $ Py.toUnicode "_max_episode_steps" >>= Py.getAttribute (env gym) >>= pyToInteger

-- | Get the number of elapsed steps in this episode.
getElapsedSteps :: Gym -> IO (Maybe Integer)
getElapsedSteps gym = python $ Py.toUnicode "_elapsed_steps" >>= Py.getAttribute (env gym) >>= pyToInteger

-- | Reset the environment.
resetGym :: Gym -> IO GymObservation
resetGym gym = do
  gObs <- python $ Py.callMethodArgs (env gym) "reset" []
  fromMaybe (error "could not convert observation to GymData") <$> python (getGymData (observationSpace gym) gObs)

-- | Take a random step and renders the environment. Also see @stepGymRenderRandom@.
stepGymRandom :: Gym -> IO GymResult
stepGymRandom = stepGymRenderRandom True

-- | Take a random step and specify wether to render the environment or not.
stepGymRenderRandom :: Bool -> Gym -> IO GymResult
stepGymRenderRandom render gym = do
  idx <- randomRIO (0, spaceSize (actionSpace gym) - 1)
  stepGymRender render gym idx

-- | Take a step as specified by the given action index. Renders the environment. See also @stepGymRender@.
stepGym :: Gym -> Integer -> IO GymResult
stepGym = stepGymRender True

-- | Take a step as specified by the given action index. Choose wether to render the environment.
stepGymRender :: Bool -> Gym -> Integer -> IO GymResult
stepGymRender render gym actIdx = do
  when render $ void $ python $ Py.callMethodArgs (env gym) "render" []
  act <- toAction (actionSpace gym) actIdx
  Just tuple <- python $ Py.callMethodArgs (env gym) "step" [Py.toObject act] >>= Py.cast
  [gObs, gReward, gDone, gInfo] <- python $ Py.fromTuple tuple
  obs <- fromMaybe (error "could not convert observation to GymData") <$> python (getGymData (observationSpace gym) gObs)
  rew <- fromMaybe (error "could not convert reward to Double") <$> python (pyToDouble gReward)
  done <- fromMaybe (error "could not convert done falg to Bool") <$> python (pyToBool gDone)
  -- Py.print gDone
  -- void $ python $ Py.print reward stdout
  -- void $ python $ Py.print gDone stdout
  -- void $ python $ Py.print info stdout
  return $ GymResult obs rew done

max_stars_repo_name: schnecki/gym-haskell

max_issues_repo_name: schnecki/gym-haskell

max_forks_repo_name: schnecki/gym-haskell


text:
{-# LANGUAGE DeriveGeneric     #-}
{-# LANGUAGE OverloadedLists   #-}
{-# LANGUAGE OverloadedStrings #-}

import qualified Data.ByteString.Char8      as S
import qualified Data.ByteString.Lazy.Char8 as L
import           Data.Complex
import qualified Data.IntMap                as IntMap
import qualified Data.IntSet                as IntSet
import qualified Data.Map                   as Map
import           Data.Ratio
import qualified Data.Sequence              as Seq
import qualified Data.Set                   as Set
import qualified Data.Text                  as T
import qualified Data.Text.Encoding         as T
import qualified Data.Text.Lazy             as TL
import qualified Data.Vector                as V
import qualified Data.Vector.Primitive      as VP
import qualified Data.Vector.Storable       as VS
import qualified Data.Vector.Unboxed        as VU
import           GHC.Generics
import           Test.Tasty
import           Test.Tasty.HUnit

import           Text.PrettyPrint.Generic

data Foo = Foo { fooA :: Int, fooB :: String } deriving Generic
instance Pretty Foo

data Bar a = Bar { barA :: Foo, barB :: a } deriving Generic
instance Pretty a => Pretty (Bar a)

data Baz = Baz1 | Baz2 Int | Baz3 Foo (Bar Baz) deriving Generic
instance Pretty Baz

main :: IO ()
main =
  defaultMain $ testGroup "test cases"
    [ testCase "()" $
      prettyShow' () @?= "()"
    , testCase "Char" $
      prettyShow' 'a' @?= "'a'"
    , testCase "Int" $
      prettyShow' (123 :: Int) @?= "123"
    , testCase "Integer" $
      prettyShow' (2^100 :: Integer) @?= "1267650600228229401496703205376"
    , testCase "Float" $
      prettyShow' (pi :: Float) @?= "3.1415927"
    , testCase "Double" $
      prettyShow' (pi :: Double) @?= "3.141592653589793"
    , testCase "String" $
      prettyShow' ("Hello" :: String) @?= "\"Hello\""
    , testCase "Bool" $
      prettyShow' True @?= "True"

    , testCase "Rational" $
      prettyShow' (123 % 456 :: Rational) @?= "41 % 152"
    , testCase "Complex" $
      prettyShow' (123 :+ 456 :: Complex Double) @?= "123.0 :+ 456.0"

    , testCase "List" $
      prettyShow' ([1..5] :: [Int]) @?= "[ 1, 2, 3, 4, 5 ]"

    , testCase "2-tuple" $
      prettyShow' ('a', 'b') @?= "('a', 'b')"
    , testCase "3-tuple" $
      prettyShow' ('a', 'b', 'c') @?= "('a', 'b', 'c')"
    , testCase "4-tuple" $
      prettyShow' ('a', 'b', 'c', 'd') @?= "('a', 'b', 'c', 'd')"
    , testCase "5-tuple" $
      prettyShow' ('a', 'b', 'c', 'd', 'e') @?= "('a', 'b', 'c', 'd', 'e')"
    , testCase "6-tuple" $
      prettyShow' ('a', 'b', 'c', 'd', 'e', 'f') @?= "('a', 'b', 'c', 'd', 'e', 'f')"
    , testCase "7-tuple" $
      prettyShow' ('a', 'b', 'c', 'd', 'e', 'f', 'g') @?= "('a', 'b', 'c', 'd', 'e', 'f', 'g')"

    , testCase "Maybe 1" $
      prettyShow' (Nothing :: Maybe Int) @?= "Nothing"
    , testCase "Maybe 2" $
      prettyShow' (Just 123 :: Maybe Int) @?= "Just 123"
    , testCase "Maybe 3" $
      prettyShow' (Just (Just 123) :: Maybe (Maybe Int)) @?= "Just (Just 123)"
    , testCase "Either 1" $
      prettyShow' (Left "Left" :: Either String Double) @?= "Left \"Left\""
    , testCase "Either 2" $
      prettyShow' (Right pi :: Either String Double) @?= "Right 3.141592653589793"

    , testCase "Strict ByteString" $
      prettyShow' (T.encodeUtf8 "日本語" :: S.ByteString) @?= "\"日本語\""
    , testCase "Lazy ByteString" $
      prettyShow' (L.fromStrict $ T.encodeUtf8 "日本語" :: L.ByteString) @?= "\"日本語\""
    , testCase "Strict Text" $
      prettyShow' ("日本語" :: T.Text) @?= "\"日本語\""
    , testCase "Lazy Text" $
      prettyShow' ("日本語" :: TL.Text) @?= "\"日本語\""

    , testCase "Set" $
      prettyShow' (["foo", "bar", "baz"] :: Set.Set String) @?= "[ \"bar\", \"baz\", \"foo\" ]"
    , testCase "IntSet" $
      prettyShow' ([1..5] :: IntSet.IntSet) @?= "[ 1, 2, 3, 4, 5 ]"
    , testCase "Map" $
      prettyShow' ([("foo", 123), ("bar", 456)] :: Map.Map String Int) @?= "{ \"bar\": 456, \"foo\": 123 }"
    , testCase "IntMap" $
      prettyShow' ([(123, "foo"), (456, "bar")] :: IntMap.IntMap String) @?= "{ 123: \"foo\", 456: \"bar\" }"
    , testCase "Seq" $
      prettyShow' ([1..5] :: Seq.Seq Int) @?= "[ 1, 2, 3, 4, 5 ]"

    , testCase "Vector" $
      prettyShow' ([1..5] :: V.Vector Int) @?= "[ 1, 2, 3, 4, 5 ]"
    , testCase "Primitive Vector" $
      prettyShow' ([1..5] :: VP.Vector Int) @?= "[ 1, 2, 3, 4, 5 ]"
    , testCase "Storable Vector" $
      prettyShow' ([1..5] :: VS.Vector Int) @?= "[ 1, 2, 3, 4, 5 ]"
    , testCase "Unboxed Vector" $
      prettyShow' ([1..5] :: VU.Vector Int) @?= "[ 1, 2, 3, 4, 5 ]"

    , testCase "User defined 1" $
      prettyShow' (Foo 123 "foo") @?= "Foo { fooA = 123, fooB = \"foo\" }"
    , testCase "User defined 2" $
      prettyShow' (Bar (Foo 123 "foo") (Just True)) @?= "Bar { barA = Foo { fooA = 123\n                 , fooB = \"foo\" }\n    , barB = Just True }"
    , testCase "User defined 3" $
      prettyShow' Baz1 @?= "Baz1"
    , testCase "User defined 4" $
      prettyShow' (Baz2 123) @?= "Baz2 123"
    , testCase "User defined 5" $
      prettyShow' (Baz3 (Foo 123 "foo") (Bar (Foo 456 "bar") Baz1)) @?= "Baz3 Foo { fooA = 123\n         , fooB = \"foo\" }\n     Bar { barA = Foo { fooA = 456\n                      , fooB = \"bar\" }\n         , barB = Baz1 }"
    ]

max_stars_repo_name: neilmayhew/generic-pretty

max_issues_repo_name: neilmayhew/generic-pretty

max_forks_repo_name: neilmayhew/generic-pretty


text:
module Statistics.GModeling.Models.LDA
  (
  ) where

import Statistics.GModeling.DSL

data LDALabels = Alpha | Beta | Topics | Topic
               | Doc | Symbols | Symbol

lda :: Network LDALabels
lda =
  [
    Only Alpha :-> Topics
  , Only Beta :-> Symbols
  , (Topics :@ Doc) :-> Topic
  , (Symbols :@ Topic) :-> Symbol
  ]

max_stars_repo_name: nanonaren/gmodeling

max_issues_repo_name: nanonaren/gmodeling

max_forks_repo_name: nanonaren/gmodeling


text:
{-# LANGUAGE ScopedTypeVariables #-}
module Main where

import Numeric.LinearAlgebra
import Control.DeepSeq
import Control.Concurrent
import Control.Monad
import Control.Exception
import System.Mem
import Control.Concurrent.Async
import System.Environment

main :: IO ()
main = do
  [matSizeS] <- getArgs
  let matSize = read matSizeS
  let gcs = do
        performMinorGC
        performMajorGC
        threadDelay 1000
  caps <- getNumCapabilities
  void $ concurrently gcs $ flip mapConcurrently [1..caps] $ \i0 -> do
    let matrices :: Int -> [Matrix Double]
        matrices i = if i == 10000
          then matrices (-10000)
          else let
            a :: Matrix Double = (matSize >< matSize) (take (matSize * matSize) [fromIntegral i..])
            in a : matrices (i+1)
    forM_ (matrices i0) $ \a -> do
      let (l, m, r) = svd a
      evaluate (force (l <> diagRect 0 m matSize matSize <> tr r))

max_stars_repo_name: bitonic/hmatrix-nasal-demons

max_issues_repo_name: bitonic/hmatrix-nasal-demons

max_forks_repo_name: bitonic/hmatrix-nasal-demons


text:
{-|
 Module      : Data.Boltzmann.System
 Description : System utilities for combinatorial specifications.
 Copyright   : (c) Maciej Bendkowski, 2017-2020

 License     : BSD3
 Maintainer  : maciej.bendkowski@tcs.uj.edu.pl
 Stability   : experimental

 General utilities for combinatorial system of algebraic and rational systems.
 -}
{-# LANGUAGE OverloadedStrings #-}
module Data.Boltzmann.System
  ( Format(..)
  , isAlgebraicF
  , isRationalF
  , System(..)
  , size
  , constructors
  , Cons(..)
  , Arg(..)
  , argName
  , types
  , PSystem(..)
  , typeList
  , initType
  , paramTypes
  , paramTypesW
  , typeWeight
  , seqTypes
  , SystemType(..)
  , systemType
  , isAlgebraic
  , isRational
  , hasAtoms
  , isAtomic
  , evalT
  , evalC
  , evalA
  , getIdx
  , value
  , eval
  , structureSize
  , SystemT
  , toSystemT
  )
where

import           Data.Set                       ( Set )
import qualified Data.Set                      as S

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

import           Numeric.LinearAlgebra   hiding ( size )

import           Data.Maybe                     ( mapMaybe )

import           Data.List                      ( nub )
import           Data.Graph

import           Data.Aeson

import           Data.Boltzmann.System.Annotations
import           Data.Boltzmann.Internal.Annotations
                                                ( withInt )

-- | Input specification format.
data Format = AlgebraicF
            | RationalF

instance Show Format where
  show AlgebraicF = "algebraic"
  show RationalF  = "rational"

isAlgebraicF :: Format -> Bool
isAlgebraicF AlgebraicF = True
isAlgebraicF _          = False

isRationalF :: Format -> Bool
isRationalF = not . isAlgebraicF

-- | System of combinatorial structures.
data System a = System { defs        :: Map String [Cons a] -- ^ Type definitions.
                       , annotations :: Annotations         -- ^ System annotations.
                       } deriving (Show)

data SystemT a = SystemT { systemTypes       :: [TypeT a]
                         , systemAnnotations :: Annotations
                         } deriving (Show)

instance ToJSON a => ToJSON (SystemT a) where
  toJSON sys =
    object ["types" .= systemTypes sys, "annotations" .= systemAnnotations sys]

data TypeT a = TypeT { typeName :: String
                     , constrs  :: [ConsT a]
                     } deriving (Show)

instance ToJSON a => ToJSON (TypeT a) where
  toJSON t = object ["name" .= typeName t, "constructors" .= constrs t]

data ConsT a = ConsT { constrName   :: String
                     , arguments    :: [ArgT]
                     , constrWeight :: a
                     , constrFreq   :: Maybe Int
                     } deriving (Show)

instance ToJSON a => ToJSON (ConsT a) where
  toJSON t = object
    [ "name" .= constrName t
    , "args" .= arguments t
    , "weight" .= constrWeight t
    , "freq" .= constrFreq t
    ]

data ArgT = ArgT { argumentName :: String
                 , argumentType :: String
                 } deriving (Show)

instance ToJSON ArgT where
  toJSON t = object ["name" .= argumentName t, "type" .= argumentType t]

-- | Converts a given system to an output format.
toSystemT :: System a -> SystemT a
toSystemT sys = SystemT { systemTypes       = map toTypeT (M.toList $ defs sys)
                        , systemAnnotations = annotations sys
                        }

toTypeT :: (String, [Cons a]) -> TypeT a
toTypeT (t, cons) = TypeT { typeName = t, constrs = map toConsT cons }

toConsT :: Cons a -> ConsT a
toConsT con = ConsT { constrName   = func con
                    , arguments    = map toArgT (args con)
                    , constrWeight = weight con
                    , constrFreq   = frequency con
                    }

toArgT :: Arg -> ArgT
toArgT (Type arg) = ArgT { argumentName = arg, argumentType = "type" }
toArgT (List arg) = ArgT { argumentName = arg, argumentType = "list" }

-- | Size of a combinatorial system.
size :: System a -> Int
size = M.size . defs

-- | First (in alphabetical order) type in system.
initType :: System a -> String
initType = head . M.keys . defs

-- | Constructors of a combinatorial system.
constructors :: System a -> Int
constructors = length . concat . M.elems . defs

-- | Type constructor.
data Cons a = Cons { func      :: String        -- ^ Constructor name.
                   , args      :: [Arg]         -- ^ Argument list.
                   , weight    :: a             -- ^ Constructor weight.
                   , frequency :: Maybe Int     -- ^ Marking parameter.
                   } deriving (Eq,Show)

-- | Type constructor arguments.
data Arg = Type String                       -- ^ Regular type reference.
         | List String                       -- ^ Type list reference.
           deriving (Eq,Show)

-- | The name of an argument.
argName :: Arg -> String
argName (Type s) = s
argName (List s) = s

-- | Type set of the given system.
types :: System a -> Set String
types = M.keysSet . defs

-- | Parametrised system of combinatorial structures.
data PSystem a = PSystem { system  :: System a      -- ^ System with probability weights.
                         , values  :: Vector a      -- ^ Numerical values of corresponding types.
                         , param   :: a             -- ^ Evaluation parameter.
                         , weights :: System Int    -- ^ System with input weights.
                         }

-- | Type list of the given parametrised system.
typeList :: PSystem a -> [String]
typeList = S.toList . M.keysSet . defs . system

-- | List of types with corresponding constructors.
paramTypes :: PSystem a -> [(String, [Cons a])]
paramTypes = M.toList . defs . system

-- | List of types with corresponding constructors and input weights.
paramTypesW :: PSystem a -> [(String, [(Cons a, Int)])]
paramTypesW sys = map (addW $ weights sys) xs where xs = paramTypes sys

addW :: System Int -> (String, [a]) -> (String, [(a, Int)])
addW sys (s, cons) = (s, zip cons ws) where ws = typeW sys s

typeW :: System Int -> String -> [Int]
typeW sys s = case s `M.lookup` defs sys of
  Just cons -> map weight cons
  Nothing   -> []

-- | Type weight of the given parametrised system.
typeWeight :: PSystem Double -> String -> Double
typeWeight sys t = vec ! idx
 where
  m   = defs $ system sys
  vec = values sys
  idx = M.findIndex t m

-- | List of sequence types.
seqTypes :: System a -> [String]
seqTypes =
  S.elems . S.fromList . concatMap seqTypesCons . concat . M.elems . defs

seqTypesCons :: Cons a -> [String]
seqTypesCons = mapMaybe listN . args
 where
  listN (List s) = Just s
  listN _        = Nothing

-- | Checks it the argument is a list.
isListArg :: Arg -> Bool
isListArg (List _) = True
isListArg _        = False

-- | Type of a combinatorial system.
--   Note: System other than rational or algebraic are not yet supported.
data SystemType = Rational
                | Algebraic
                | Unsupported String   -- ^ error message

instance Show SystemType where
  show Rational        = "rational"
  show Algebraic       = "algebraic"
  show (Unsupported _) = "unsupported"

isAlgebraic :: SystemType -> Bool
isAlgebraic Algebraic = True
isAlgebraic _         = False

isRational :: SystemType -> Bool
isRational Rational = True
isRational _        = False

-- | Determines the system type.
systemType :: System a -> SystemType
systemType sys
  | not (isLinear sys)
  = Algebraic
  | not (isInterruptible sys)
  = Unsupported "Given rational system is not interruptible."
  | otherwise
  = let depGraph = dependencyGraph sys
    in  case scc depGraph of
          [_] -> Rational
          xs ->
            Unsupported
              $  "Given rational system has "
              ++ show (length xs)
              ++ " strongly connected components."

-- | Constructs a dependency graph for the given system.
dependencyGraph :: System a -> Graph
dependencyGraph sys = buildG (0, n + d - 1) (edgs ++ edgs')
 where
  idx s = M.findIndex s (defs sys)
  idxSeq s = n + S.findIndex s seqsSet
  edgs    = concatMap (edges' atomicT idx idxSeq) $ M.toList (defs sys)
  edgs'   = concatMap (\t -> [(idxSeq t, idxSeq t), (idxSeq t, idx t)]) seqs
  atomicT = atomicTypes sys
  seqsSet = S.fromAscList seqs
  seqs    = seqTypes sys
  d       = S.size seqsSet
  n       = size sys

edges'
  :: Set String
  -> (String -> Int)
  -> (String -> Int)
  -> (String, [Cons b])
  -> [(Vertex, Vertex)]
edges' atomicT idx idxSeq (t, cons) = concatMap edge' $ neighbours cons
 where
  tidx       = idx t
  neighbours = nub . concatMap args
  edge' (List s) = [(tidx, idxSeq s)]
  edge' (Type s) | s `S.member` atomicT = [(tidx, idx s), (idx s, tidx)]
                 | -- double edge
                   otherwise            = [(tidx, idx s)]

-- | Checks whether the system is linear, i.e.
--   each constructor references at most one type.
isLinear :: System a -> Bool
isLinear sys = all (all linear) (M.elems $ defs sys)
 where
  atomicT = atomicTypes sys
  linear cons =
    not (any isListArg $ args cons)
      && length (compoundArgs atomicT $ args cons)
      <= 1

-- | Determines whether each constructor in the system has at most one atom.
--   Note: the system is assumed to contain some atoms (see hasAtoms).
isInterruptible :: System a -> Bool
isInterruptible sys = all interruptible' $ M.elems (defs sys)
  where interruptible' cons = length (filter isAtomic cons) <= 1

compoundArgs :: Set String -> [Arg] -> [Arg]
compoundArgs atomicT = filter (\x -> argName x `S.notMember` atomicT)

-- | Determines the set of "atomic" types.
atomicTypes :: System a -> Set String
atomicTypes sys = S.fromList $ map fst ts
 where
  ts = filter isAtomic' $ M.toList (defs sys)
  isAtomic' (_, cons) = all isAtomic cons

isAtomic :: Cons a -> Bool
isAtomic = null . args

-- | Determines whether the system has atoms.
hasAtoms :: System a -> Bool
hasAtoms sys = any (any isAtomic) $ M.elems (defs sys)

-- | Evaluates the type in the given coordinates.
evalT :: System Int -> Double -> Vector Double -> [Cons Int] -> Double
evalT sys z ys cons = sum $ map (evalC sys z ys) cons

-- | Evaluates the constructor in the given coordinates.
evalC :: System Int -> Double -> Vector Double -> Cons Int -> Double
evalC sys z ys con = foldl (*) start $ map (evalA sys ys) (args con)
 where
  w     = weight con
  start = if w > 0 then z ^^ w else 1

-- | Evaluates the argument in the given coordinates.
evalA :: System Int -> Vector Double -> Arg -> Double
evalA sys ys (Type t) = ys ! getIdx sys t
evalA sys ys (List t) = recip $ 1 - ys ! getIdx sys t

getIdx :: System Int -> String -> Int
getIdx sys x = x `M.findIndex` defs sys

value :: String -> System b -> Vector Double -> Double
value t sys vec = vec ! M.findIndex t (defs sys)

-- | Evaluates the system at the given coordinates.
eval :: System Int -> Vector Double -> Double -> Vector Double
eval sys ys z = n |> map update [0 .. n]
 where
  n = size sys
  f k = snd $ M.elemAt k (defs sys)
  update idx = evalT sys z ys $ f idx

structureSize :: System a -> Int
structureSize sys = withInt (annotations sys) "size" 1000

max_stars_repo_name: timsears/boltzmann-brain

max_issues_repo_name: timsears/boltzmann-brain

max_forks_repo_name: timsears/boltzmann-brain


text:
-----------------------------------------------------------------------------
{- |
Module      :  Data.Packed
Copyright   :  (c) Alberto Ruiz 2006-2010
License     :  GPL-style

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional
Portability :  uses ffi

Types for dense 'Vector' and 'Matrix' of 'Storable' elements.

-}
-----------------------------------------------------------------------------

module Data.Packed (
    module Data.Packed.Vector,
    module Data.Packed.Matrix,
--    module Numeric.Conversion,
--    module Data.Packed.Random,
--    module Data.Complex
) where

import Data.Packed.Vector
import Data.Packed.Matrix
--import Data.Packed.Random
--import Data.Complex
--import Numeric.Conversion

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
module Trade.Statistics.Summary where

import Control.Applicative (liftA2)

import qualified Data.Map as Map

import qualified Data.Vector as Vec

import qualified Statistics.Sample as Sample

import qualified Trade.Type.DeltaSignal.Algorithm as DSA
import Trade.Type.DeltaSignal (DeltaSignal(..))

import Trade.Type.NestedMap (NestedMap(..))

import Trade.Type.Position (Position(..))
import Trade.Type.WinningLosing (WinningLosing(..))

import Trade.Type.Yield (LogYield(..))

import Trade.Report.Pretty (pretty)
import qualified Trade.Report.Table as Table
import Trade.Report.ToReport (ToReport, toReport)

data Summary = Summary {
  totalTrades :: Int
  , profitFactor :: Double
  , avgProfit :: Double
  , stdDevProfit :: Double
  , winners :: Int
  , largestWinner :: Double
  , avgWinner :: Double
  , losers :: Int
  , largestLoser :: Double
  , avgLoser :: Double
  , cumulativeProfit :: Double
  }

summary :: [LogYield ohlc] -> [LogYield ohlc] -> Summary
summary wins losses =
  let ws = Vec.fromList (map logYield wins)
      ls = Vec.fromList (map logYield losses)
  in Summary {
    totalTrades = length wins + length losses
    , profitFactor = exp (sum ws) / exp (sum ls)
    , avgProfit = exp (Sample.mean (Vec.concat [ws, ls]))
    , stdDevProfit = exp (Sample.stdDev ws)

    , winners = Vec.length ws
    , largestWinner = exp (Vec.maximum ws)
    , avgWinner = exp (Sample.mean ws)

    , losers = Vec.length ls
    , largestLoser = exp (Vec.minimum ls)
    , avgLoser = exp (Sample.mean ls)

    , cumulativeProfit = exp (sum ws + sum ls)
    }

summary2table :: Summary -> [[String]]
summary2table ys =
  [ [ "Summary" ]
  , []
  , [ "Total trades", pretty (totalTrades ys) ]
  , [ "Profit factor", pretty (profitFactor ys) ]
  , [ "Avg. profit", pretty (avgProfit ys) ]
  , [ "StdDev. profit", pretty (stdDevProfit ys) ]
  , []
  , [ "Winners", pretty (winners ys) ]
  , [ "Largest winner", pretty (largestWinner ys) ]
  , [ "Avg. winner", pretty (avgWinner ys) ]
  , []
  , [ "Losers", pretty (losers ys) ]
  , [ "Largest loser", pretty (largestLoser ys) ]
  , [ "Avg. loser", pretty (avgLoser ys) ]
  , []
  , [ "Cumulativ profit", pretty (cumulativeProfit ys) ]
  ]

instance ToReport Summary where
  toReport = Table.table . summary2table

toSummary ::
  NestedMap Position WinningLosing [DeltaSignal ohlc] -> Maybe Summary
toSummary (NestedMap m) =
  let nm = fmap (fmap (map DSA.yield)) m
      win = Map.lookup Invested nm >>= Map.lookup Winning
      lose = Map.lookup Invested nm >>= Map.lookup Losing
  in liftA2 summary win lose

max_stars_repo_name: fphh/trade

max_issues_repo_name: fphh/trade

max_forks_repo_name: fphh/trade


text:
module Statistics.Information.Discrete.Entropy where

import Data.Function.Memoize
import Data.Matrix
import Statistics.Information.Utils.List

entropy :: Eq a => Int -> Matrix a -> Double
entropy = memoize $ \base xs ->
  -(sum (map elog . hist . toLists $ xs)) / log (fromIntegral base)
  where
    elog x | x <= 0.0 = 0
    elog x | x >= 1.0 = 0
    elog x = x * log x

centropy :: Eq a => Int -> Matrix a -> Matrix a -> Double
centropy base xs ys = entropy base (xs <|> ys) - entropy base ys

max_stars_repo_name: eligottlieb/Shannon

max_issues_repo_name: eligottlieb/Shannon

max_forks_repo_name: eligottlieb/Shannon


text:
import Criterion.Main
import Data.Default.Class
import qualified Data.Vector.Unboxed as U
import Text.Printf

import Numeric.SpecFunctions
import Numeric.Polynomial
import Numeric.RootFinding

-- Uniformly sample logGamma performance between 10^-6 to 10^6
benchmarkLogGamma logG =
  [ bench (printf "%.3g" x) $ nf logG x
  | x <- [ m * 10**n | n <- [ -8 .. 8 ]
                     , m <- [ 10**(i / tics) | i <- [0 .. tics-1] ]
         ]
  ]
  where tics = 3
{-# INLINE benchmarkLogGamma #-}

-- Power of polynomial to be evaluated (In other words length of coefficients vector)
coef_size :: [Int]
coef_size = [ 1,2,3,4,5,6,7,8,9
            , 10,    30
            , 100,   300
            , 1000,  3000
            , 10000, 30000
            ]
{-# INLINE coef_size #-}

-- Precalculated coefficients
coef_list :: [U.Vector Double]
coef_list = [ U.replicate n 1.2 | n <- coef_size]
{-# NOINLINE coef_list #-}

main :: IO ()
main = defaultMain
  [ bgroup "logGamma" $
    benchmarkLogGamma logGamma
  , bgroup "logGammaL" $
    benchmarkLogGamma logGammaL
  , bgroup "incompleteGamma" $
      [ bench (show p) $ nf (incompleteGamma p) p
      | p <- [ 0.1
             , 1,   3
             , 10,  30
             , 100, 300
             , 999, 1000
             ]
      ]
  , bgroup "factorial"
    [ bench (show n) $ nf factorial n
    | n <- [ 0, 1, 3, 6, 9, 11, 15
           , 20, 30, 40, 50, 60, 70, 80, 90, 100
           ]
    ]
  , bgroup "incompleteBeta"
    [ bench (show (p,q,x)) $ nf (incompleteBeta p q) x
    | (p,q,x) <- [ (10,      10,      0.5)
                 , (101,     101,     0.5)
                 , (1010,    1010,    0.5)
                 , (10100,   10100,   0.5)
                 , (100100,  100100,  0.5)
                 , (1001000, 1001000, 0.5)
                 , (10010000,10010000,0.5)
                 ]
    ]
  , bgroup "log1p"
      [ bench (show x) $ nf log1p x
      | x <- [ -0.9
             , -0.5
             , -0.1
             ,  0.1
             ,  0.5
             ,  1
             ,  10
             ,  100
             ] :: [Double]
      ]
  , bgroup "sinc" $
        bench "sin" (nf sin (0.55 :: Double))
      : [ bench (show x) $ nf sinc x
        | x <- [0, 1e-6, 1e-3,  0.5]
        ]
  , bgroup "poly"
      $  [ bench ("vector_"++show (U.length coefs)) $ nf (\x -> evaluatePolynomial x coefs) (1 :: Double)
         | coefs <- coef_list ]
      ++ [ bench ("unpacked_"++show n) $ nf (\x -> evaluatePolynomialL x (map fromIntegral [1..n])) (1 :: Double)
         | n <- coef_size ]
  , bgroup "RootFinding"
    [ bench "ridders sin" $ nf (ridders       def (0,pi/2))     (\x -> sin x - 0.525)
    , bench "newton sin"  $ nf (newtonRaphson def (0,1.2,pi/2)) (\x -> (sin x - 0.525,cos x))
    ]
  ]