text:
module SchemeParser
  (readExpr,
   readExprList) where

import Text.ParserCombinators.Parsec hiding (spaces)
import Control.Monad
import Numeric
import Data.Ratio
import Data.Complex
import qualified Data.Vector as V
import CommonTypes
import Control.Monad.Except

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

spaces :: Parser ()
spaces = skipMany space

escapedCharacters :: Parser Char
escapedCharacters = do
                      char '\\'
                      x <- oneOf "\\\"nrt"
                      return $ case x of
                                 '\\' -> x
                                 '"'  -> x
                                 'n'  -> '\n'
                                 't'  -> '\t'
                                 'r'  -> '\r'

parseString :: Parser LispVal
parseString = do
               char '"'
               x <- many $ escapedCharacters <|> noneOf "\\\""
               char '"'
               return $ String x

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

parseAtom :: Parser LispVal
parseAtom = do
             first <- letter <|> symbol
             rest <- many $ letter <|> digit <|> symbol
             let atom = first:rest
             return $ Atom atom

parseDecimal :: Parser LispVal
parseDecimal = Number . read <$> many1 digit

myParseDecimal :: Parser LispVal
myParseDecimal = many1 digit >>= (return . Number . read)

parseDecimal2 :: Parser LispVal
parseDecimal2 = do
                  try $ string "#d"
                  x <- many1 digit
                  (return . Number . read) x

parseOct :: Parser LispVal
parseOct = do
              try $ string "#o"
              x <- many1 octDigit
              (return . Number . oct2dig) x

parseHex :: Parser LispVal
parseHex = do
            try $ string "#x"
            x <- many1 hexDigit
            (return . Number . hex2dig) x

parseBin :: Parser LispVal
parseBin = do
            try $ string "#b"
            x <- many1 $ oneOf "10"
            (return . Number . bin2dig) x

oct2dig x = fst $ head $ readOct x

hex2dig x = fst $ head $ readHex x

bin2dig  = bin2dig' 0
bin2dig' digint "" = digint
bin2dig' digint (x:xs) = let old = 2 * digint + (if x == '0' then 0 else 1)
                            in bin2dig' old xs

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

parseCharacter :: Parser LispVal
parseCharacter = do
                  try $ string "#\\"
                  value <- try (string "newline" <|> string "space")
                    <|> do { x <- anyChar; notFollowedBy alphaNum ; return [x]}
                  return $ Character $ case value of
                                         "space"    -> ' '
                                         "newline"  -> '\n'
                                         otherwise  -> (head value)

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

parseRatio :: Parser LispVal
parseRatio = do
              x <- many1 digit
              char '/'
              y <- many1 digit
              return $ Ratio (read x % read y)

toDouble :: LispVal -> Double
toDouble (Float d) = d
toDouble (Number n) = fromInteger n

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

parseList :: Parser LispVal
parseList = liftM List $ sepBy parseExpr spaces

parseDottedList :: Parser LispVal
parseDottedList = do
                   head <- endBy parseExpr spaces
                   tail <- char '.' >> spaces >> parseExpr
                   return $ DottedList head tail

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

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

parseUnQuote :: Parser LispVal
parseUnQuote = do
                char ','
                x <- parseExpr
                return $ List [Atom "unquote", x]

parseVector :: Parser LispVal
parseVector = do
                vectorValues <- sepBy parseExpr spaces
                return $ Vector (V.fromList vectorValues)

parseExpr :: Parser LispVal
parseExpr = parseAtom
            <|> parseString
            <|> try parseComplex
            <|> try parseFloat
            <|> try parseRatio
            <|> try parseNumber
            <|> try parseBool
            <|> try parseCharacter
            <|> parseQuoted
            <|> parseQuasiQuoted
            <|> parseUnQuote
            <|> try  (do
                        string "#("
                        x <- parseVector
                        char ')'
                        return x)
            <|> do
                  char '('
                  x <- try parseList <|> parseDottedList
                  char ')'
                  return x

readOrThrow :: Parser a -> String -> ThrowsError a
readOrThrow parser input = case parse parser "lisp" input of
                   Left err  -> throwError $ Parser err
                   Right val -> return val

readExpr = readOrThrow parseExpr
readExprList = readOrThrow (endBy parseExpr spaces)


text:
{-# OPTIONS_GHC -Wno-incomplete-patterns #-}

module VectorSpec where

import Data.Complex

import Test.QuickCheck
import Test.QuickCheck.Poly

import Category
import Comonad
import Functor
import Unboxed
import Vector

prop_Vector_Functor_id :: FnProp (Vector A)
prop_Vector_Functor_id = checkFnEqual law_Functor_id

prop_Vector_Functor_comp :: Fun B C -> Fun A B -> FnProp (Vector A)
prop_Vector_Functor_comp f g = checkFnEqual (law_Functor_comp f g)

prop_Vector_Semicomonad_comm ::
    Fun (Vector B) C -> Fun (Vector A) B -> FnProp (Vector A)
prop_Vector_Semicomonad_comm f g = checkFnEqual (law_Semicomonad_comm f g)

prop_Vector_Semicomonad1_comm ::
    Fun (Vector B) C -> Fun (Vector A) B -> FnProp (Vector A)
prop_Vector_Semicomonad1_comm f g = checkFnEqual (law_Semicomonad1_comm f g)

prop_Vector_Semicomonad1_comm' ::
    Fun (Vector B) C -> Fun (Vector A) B -> FnProp (Vector A)
prop_Vector_Semicomonad1_comm' (Fn f) (Fn g) =
    checkFnEqual (law_Semicomonad1_comm' f g)

type UA = Int
type UB = Double
type UC = Complex Double

prop_UVector_Functor_id :: FnProp (UVector UA)
prop_UVector_Functor_id = checkFnEqual law_Functor_id

prop_UVector_Functor_comp ::
    (UB -#> UC) -> (UA -#> UB) -> FnProp (UVector UA)
prop_UVector_Functor_comp f g = checkFnEqual (law_Functor_comp f g)

prop_UVector_Semicomonad1_comm' ::
    Fun (UVector UB) UC -> Fun (UVector UA) UB -> FnProp (UVector UA)
prop_UVector_Semicomonad1_comm' (Fn f) (Fn g) =
    checkFnEqual (law_Semicomonad1_comm' f g)

ext: hs


text:
{-# LANGUAGE FlexibleContexts #-}

module STC.OrientationScaleAnalysis where

import           Control.Monad.Parallel    as MP
import           Data.Array.Repa           as R
import           Data.Complex
import           Data.List                 as L
import           Data.Vector.Unboxed       as VU
import           DFT.Plan
import           FokkerPlanck.DomainChange (r2z1Tor2s1)
import           Graphics.Gnuplot.Simple
import           Image.IO                  (ImageRepa (..), plotImageRepa,
                                            plotImageRepaComplex)
import           STC.Convolution
import           System.FilePath
import           Types
import           Utils.Array
import           Utils.Parallel
import Text.Printf

-- {-# INLINE analyzeOrientation #-}
-- analyzeOrientation :: Int -> [Double] -> R2T0Array -> (R2T0Array, R2T0Array)
-- analyzeOrientation numOrientation thetaFreqs arr =
--   let arrR2S1 = R.map magnitude . r2z1Tor2s1 numOrientation thetaFreqs $ arr
--       deltaTheta = 2 * pi / fromIntegral numOrientation :: Double
--       (Z :. _ :. cols :. rows) = extent arrR2S1
--       orientationArr =
--         fromListUnboxed (Z :. cols :. rows) .
--         parMap
--           rdeepseq
--           (\(i, j) ->
--              deltaTheta *
--              (fromIntegral .
--               VU.maxIndex . toUnboxed . computeS . R.slice arrR2S1 $
--               (Z :. All :. i :. j))) $
--         [(i, j) | i <- [0 .. cols - 1], j <- [0 .. rows - 1]]
--       freqArr = fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs
--       magArr =
--         R.map sqrt . sumS . R.map (\x -> (magnitude x) ^ 2) . rotate3D $ arr
--       func theta =
--         computeS .
--         R.traverse3 orientationArr freqArr magArr (\_ _ _ -> extent arr) $ \fOri fFreq fMag (Z :. k :. i :. j) ->
--           (fMag (Z :. i :. j) :+ 0) *
--           (exp $ 0 :+ (1) * fFreq (Z :. k) * (fOri (Z :. i :. j) + theta))
--    in (func 0, func 0)

-- {-# INLINE analyzeOrientationR2Z1T0 #-}
-- analyzeOrientationR2Z1T0 ::
--      Int -> [Double] -> [Double] -> R2T0Array -> (R2Z1T0Array, R2Z1T0Array)
-- analyzeOrientationR2Z1T0 numOrientation thetaFreqs theta0Freqs arr =
--   let arrR2S1 = R.map magnitude . r2z1Tor2s1 numOrientation thetaFreqs $ arr
--       deltaTheta = 2 * pi / fromIntegral numOrientation :: Double
--       (Z :. _ :. cols :. rows) = extent arrR2S1
--       orientationArr =
--         fromListUnboxed (Z :. cols :. rows) .
--         parMap
--           rdeepseq
--           (\(i, j) ->
--              deltaTheta *
--              (fromIntegral .
--               VU.maxIndex . toUnboxed . computeS . R.slice arrR2S1 $
--               (Z :. All :. i :. j))) $
--         [(i, j) | i <- [0 .. cols - 1], j <- [0 .. rows - 1]]
--       freqArr = fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs
--       freq0Arr = fromListUnboxed (Z :. (L.length theta0Freqs)) theta0Freqs
--       magArr =
--         R.map sqrt . sumS . R.map (\x -> (magnitude x) ^ 2) . rotate3D $ arr
--       func theta =
--         computeS .
--         R.traverse4
--           orientationArr
--           freqArr
--           freq0Arr
--           magArr
--           (\_ _ _ _ ->
--              (Z :. (L.length thetaFreqs) :. (L.length theta0Freqs) :. cols :.
--               rows)) $ \fOri fFreq fFreq0 fMag (Z :. k :. l :. i :. j) ->
--           (fMag (Z :. i :. j) :+ 0) *
--           (exp $ 0 :+ (-1) * (fFreq (Z :. k) + fFreq0 (Z :. l)) * (fOri (Z :. i :. j) + theta))
--    in (func (pi / 2),func 0)

-- {-# INLINE normalizeList #-}
-- normalizeList :: (Ord e, Fractional e) => [e] -> [e]
-- normalizeList xs = L.map (/ L.maximum xs) xs

-- plotMagnitudeOrientation ::
--      (R.Source s (Complex Double))
--   => FilePath
--   -> Int
--   -> [Double]
--   -> R.Array s DIM3 (Complex Double)
--   -> (Int, Int)
--   -> IO (Int,Int)
-- plotMagnitudeOrientation folderPath numOrientationSample thetaFreqs arr (i', j') = do
--   let orientationSampleRad =
--         [ 2 * pi / fromIntegral numOrientationSample * fromIntegral i
--         | i <- [0 .. numOrientationSample - 1]
--         ]
--       orientationSampleDeg =
--         [ 360 * fromIntegral i / fromIntegral numOrientationSample
--         | i <- [0 .. numOrientationSample - 1]
--         ]
--       freqArr = fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs
--       xsFreqDomain
--         -- normalizeList $
--        =
--         parMap
--           rdeepseq
--           (\theta ->
--              magnitude .
--              R.sumAllS .
--              R.zipWith (\freq x -> x * exp (0 :+ theta * freq)) freqArr .
--              R.slice arr $
--              (Z :. All :. i :. j))
--           orientationSampleRad
--       xs =
--         normalizeList $
--         R.toList .
--         R.map magnitude .
--         r2z1Tor2s1 numOrientationSample thetaFreqs .
--         extend (Z :. All :. (1 :: Int) :. (1 :: Int)) . R.slice arr $
--         (Z :. All :. i :. j)
--       (Z :. _ :. cols :. rows) = extent arr
--       magVec =
--         VU.concat $
--         parMap
--           rdeepseq
--           (\theta ->
--              toUnboxed .
--              computeS .
--              R.map magnitude . R.sumS . rotate3D . R.traverse2 arr freqArr const $ \f1 f2 idx@(Z :. k :. _ :. _) ->
--                f1 idx * exp (0 :+ theta * f2 (Z :. k)))
--           orientationSampleRad
--       (Z :. c :. a :. b) =
--         fromIndex (Z :. (L.length thetaFreqs) :. cols :. rows) . VU.maxIndex $
--         magVec
--       maxMag = VU.maximum magVec
--       (i,j) = (a,b)
--   printf
--     "Max magnitude: %0.5f at (%d,%d) %f degree.\n"
--     maxMag
--     a
--     b
--     ((fromIntegral c :: Double) / (fromIntegral numOrientationSample) * 360)
--   plotPathsStyle
--     [ PNG (folderPath </> "Magnitude.png")
--     , Title ("Magnitude at " L.++ show (i, j))
--     ] $
--     L.zip
--         -- defaultStyle
--       --     { plotType = LinesPoints
--       --     , lineSpec = CustomStyle [LineTitle "Spatial Domain", PointType 1]
--       --     }
--       -- ,
--       [ defaultStyle
--           { plotType = Lines
--           , lineSpec =
--               CustomStyle
--                 [ LineTitle "Frequency Domain" -- , PointType 0
--                 ]
--           }
--       ]
--        -- L.zip orientationSampleDeg xs,
--       [L.zip orientationSampleDeg xsFreqDomain]
--   return (a,b)

-- plotMagnitudeOrientationSource ::
--      (R.Source s (Complex Double))
--   => DFTPlan
--   -> FilePath
--   -> Int
--   -> Int
--   -> [Double]
--   -> R2Z1T0Array
--   -> R.Array s DIM3 (Complex Double)
--   -> (Int, Int)
--   -> IO ()
-- plotMagnitudeOrientationSource plan folderPath numOrientationSample numOrientation thetaFreqs filter input (i, j) = do
--   let orientationSampleRad =
--         [ 2 * pi / fromIntegral numOrientationSample * fromIntegral i
--         | i <- [0 .. numOrientationSample - 1]
--         ]
--       orientationSampleDeg =
--         [ 360 * fromIntegral i / fromIntegral numOrientationSample
--         | i <- [0 .. numOrientationSample - 1]
--         ]
--       freqArr = fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs
--       (Z :. numThetaFreq :. cols :. rows) = extent input
--   xs <-
--     MP.mapM
--       (\theta -> do
--          initialDistF <-
--            fmap (fromUnboxed (Z :. numThetaFreq :. cols :. rows) . VU.convert) .
--            dftExecute plan (DFTPlanID DFT1DG [numThetaFreq, cols, rows] [1, 2]) .
--            VU.convert . toUnboxed . computeS . R.traverse2 input freqArr const $
--            (\f1 f2 idx@(Z :. k :. i :. j) ->
--               f1 idx * exp (0 :+ theta * (1) * f2 (Z :. k)))
--          sourceArr <- convolveR2T0 plan filter initialDistF
--          let sourceR2Z1 = R.sumS . rotateR2Z1T0Array $ sourceArr
--              mag =
--                magnitude .
--                R.sumAllS .
--                r2z1Tor2s1 numOrientation thetaFreqs .
--                extend (Z :. All :. (1 :: Int) :. (1 :: Int)) .
--                R.slice sourceR2Z1 $
--                (Z :. All :. i :. j)
--          return mag)
--       orientationSampleRad
--   plotPath
--     [ PNG (folderPath </> "SourceMagnitude.png")
--     , Title ("Source Magnitude at " L.++ show (i, j))
--     ] .
--     L.zip orientationSampleDeg . normalizeList $
--     xs

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

import Data.Complex
import Data.List (genericLength)

-- |compute x to the n-th power
power :: Complex Double -> Integer -> Complex Double
power x 0 = 1
power x n | n > 0 = x * power x (n - 1)

examplePair :: (Double, Bool)
examplePair = (3.14, True)

exampleTriple :: (Bool, Integer, String)
exampleTriple = (False, 42, "hello, world")

exampleFunction :: (Bool, Integer, String ) -> Bool
exampleFunction (b, i, s) = b || fromInteger i > length s
-- alternative: ... = b || i > genericLength s

exampleTriple' :: (Bool, Int, String)
exampleTriple' = (False, 42, "hello, world")

exampleFunction' :: (Bool, Int, String ) -> Bool
exampleFunction' (b, i, s) = b || i > length s

summerize :: [String] -> String
summerize [] = "None"
summerize [x] = "Only " ++ x
summerize [x,y] = "Two things: " ++ x ++ " and " ++ y
summerize [_ ,_ ,_ ] = "Three things: ???"
summerize _ = "Several things."

-- |doubles every element of a list of numbers
doubles :: [Integer] -> [Integer]
doubles [] = []
doubles (x:xs) = (2 * x) : doubles xs

-- |map a function on all elements of a list
mymap :: (a -> b) -> [a] -> [b]
mymap f [] = []
mymap f (x:xs) = f x : mymap f xs

double x = 2 * x

doubles' xs = mymap double xs

-- task: use quickcheck to compare doubles and doubles'

max_stars_repo_name: proglang/FunctionalProgramming

max_issues_repo_name: proglang/FunctionalProgramming

max_forks_repo_name: proglang/FunctionalProgramming


text:
-- |
-- Module    : GRN.Uncertainty
-- Copyright : (c) 2011 Jason Knight
-- License   : BSD3
--
-- Maintainer  : jason@jasonknight.us
-- Stability   : experimental
-- Portability : portable
--
-- Print out a set of boolean networks that represent the uncertainty class.
--

module GRN.Uncertainty where

import GRN.StateTransition
import GRN.Render
import GRN.Sparse
import GRN.Types
import GRN.Utils
import GRN.EM

import qualified Data.Map as M
import Data.List
import Data.Char
import Data.Bits
import Control.Monad
import Text.Printf
import System.Console.ParseArgs

import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Generic as G

import Control.Parallel.Strategies
import Data.Vector.Strategies

import Data.Graph.Inductive
import Data.Graph.Analysis.Algorithms.Common (componentsOf)

import Statistics.Distribution
import Statistics.Distribution.Normal

import Graphics.Gnuplot.Simple
import Text.PrettyPrint

-- | Collapses SSD onto space of non-input nodes
collapseSSD :: ParseData -> SSD -> SSD
collapseSSD pdata orig = snd $ foldl' splitAndSum (pdata,orig) inputSet where
  inputSet = inputGeneNames pdata

-- | Partition an SSD based on the value of a gene and sum the SSD (reducing
-- its size by 2 in the process).
splitAndSum :: (ParseData,SSD) -> Gene -> (ParseData,SSD)
splitAndSum (pd, sd) gene = (M.delete gene pd, U.zipWith (filtSum) sp1 sp2) where
  filtSum (_,x1) (_,x2) = x1+x2 -- Have to get rid of the indices
  (sp1,sp2) = U.partition evalGene (U.indexed sd) -- Wish there was a U.ipartition
  evalGene (ind,val) = testBit ind (numGenes-geneIndex-1)
  geneIndex = M.findIndex gene pd
  numGenes = M.size pd

-- | Find the names of the input genes
inputGeneNames :: ParseData -> [Gene]
inputGeneNames pdata = filter (inputGene pdata) (M.keys pdata)

-- | Use this to find input genes by the parse data.
inputGene :: ParseData -> Gene -> Bool
inputGene pdata g = null $ depends gi
  where
    Just gi = M.lookup g pdata

computeSSD :: IO ()
computeSSD = do
    let r file = liftM (U.fromList . head . read) $ readFile file :: IO (U.Vector Double)
    con1a <- r "moham_100.dat"
    con1b <- r "moham_111.dat"

    con2a <- r "moham_111.dat"
    con2b <- r "moham_101.dat"

    con3a <- r "moham_x1x.dat"
    con3b <- r "moham_1xx.dat"

    let ssd1 = U.sum . U.zipWith (\x y->(min x y)*0.5) con1a $ con1b
        ssd2 = U.sum . U.zipWith (\x y->(min x y)*0.5) con2a $ con2b
        ssd3 = U.sum . U.zipWith (\x y->(min x y)*0.5) con3a $ con3b

    print ssd1
    print ssd2
    print ssd3

uncertaintyPrint :: Args String -> ParseData -> IO ()
uncertaintyPrint args p = do
  let
    {-gen = gotArg args "generate"-}
    ks = buildKmaps p
    ms = measureMap p
    n1 = getRequiredArg args "n1" :: Int
    uvals = uncertainVals ks

    ng = M.size ks
    n = length uvals

    {-hypercorners = map (U.fromList.map fromIntegral.dec2bin n) [0..2^n-1]-}

    xs = clamped [0,0.5,1]
    center = replicate n 0.5
    clamped xs = map (\x-> if (x>1.0) then 1.0 else (if (x<0.0) then 0.0 else x)) xs
    grid = [[]] >>= foldr (>=>) return (replicate n (\x-> [x++[y] | y <- xs]))
    specialgrid = center : (grid \\ [center])
    hypergrid = map (U.fromList) specialgrid

    kmaps = map (fillKentries uvals ks) hypergrid
    ssds = parMap rdeepseq (collapseSSD p . simulateDOKUnif args . flip kmapToDOK 0) kmaps

  putStr (pprint ssds)
  {-plotList [] (U.toList . head $ ssds)-}

  {-print $ map U.length ssds-}
  {-putStrLn ""-}
  {-print $ length ssds-}
  {-print $ map (fst . normalizeSSD) ssds-}
  {-print $ "Optimum Theta: " ++ show optimumTheta-}

  return ()

class Pretty a where
  prettyPrint :: a -> String
  prettyPrint = renderStyle (Style LeftMode 80 0) . prettyShow
  prettyShow :: a -> Doc

pprint :: Pretty a => a -> String
pprint = prettyPrint

instance Pretty Int where
  prettyShow = text . show
instance Pretty Double where
  prettyShow = text . printf "%f"
instance Pretty a => Pretty [a] where
  prettyShow = brackets . hcat . punctuate comma . map prettyShow
instance (Pretty a, U.Unbox a) => Pretty (U.Vector a) where
  prettyShow = brackets . hcat . punctuate comma . map prettyShow . U.toList
instance (Pretty a, Pretty b) => Pretty (a,b) where
  prettyShow (a,b) = prettyShow a <> colon <> prettyShow b

max_stars_repo_name: binarybana/grn-pathways

max_issues_repo_name: binarybana/grn-pathways

max_forks_repo_name: binarybana/grn-pathways


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE InstanceSigs          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE Strict                #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}

module Grenade.Layers.FullyConnected (
    FullyConnected (..)
  , FullyConnected' (..)
  , randomFullyConnected
  , SpecFullyConnected (..)
  , specFullyConnected
  , fullyConnected
  ) where

import           Control.DeepSeq
import           Control.Monad
import           Control.Monad.Primitive        (PrimBase, PrimState)
import           Control.Parallel.Strategies
import           Data.Maybe                     (fromMaybe)
import           Data.Reflection                (reifyNat)
import qualified Data.Vector.Storable           as V
import           GHC.Generics                   (Generic)
import           GHC.TypeLits
import           System.Random.MWC              hiding (create)
#if MIN_VERSION_singletons(2,6,0)
import           Data.Singletons.TypeLits       (SNat (..))
#endif
import           Data.List                      (foldl')
import           Data.Proxy
import           Data.Serialize
import           Data.Singletons
import           Data.Singletons.Prelude.Num    ((%*))
import qualified Numeric.LinearAlgebra          as LA
import           Numeric.LinearAlgebra.Static   hiding (zipWithVector)
import           Text.Printf

import           Control.Monad                  (void)
import           Foreign.Storable               (peekElemOff, pokeElemOff, sizeOf)
import           System.IO.Unsafe               (unsafePerformIO)

import           Grenade.Core
import           Grenade.Dynamic
import           Grenade.Dynamic.Internal.Build
import           Grenade.Layers.Internal.BLAS
import           Grenade.Layers.Internal.CUDA
import           Grenade.Layers.Internal.Update
import           Grenade.Types
import           Grenade.Utils.Conversion
import           Grenade.Utils.LinearAlgebra
import           Grenade.Utils.ListStore
import           Grenade.Utils.Vector

import           Debug.Trace

-- | A basic fully connected (or inner product) neural network layer.
data FullyConnected i o = FullyConnected
                        !(FullyConnected' i o)             -- Neuron weights
                        !(ListStore (FullyConnected' i o)) -- momentum store
                        deriving (Generic)

instance NFData (FullyConnected i o) where
  rnf (FullyConnected w store) = rnf w `seq` rnf store

instance Show (FullyConnected i o) where
  show FullyConnected {} = "FullyConnected"

-- | How to store the data .
data FullyConnected' i o
  = FullyConnectedHMatrix
      !(R o)   -- ^ Bias
      !(L o i) -- ^ Activations
  | FullyConnectedBLAS
    !(Int, Int)         -- ^ Input, output
    !(V.Vector RealNum) -- ^ Bias, Temporary vector of same size
    !(V.Vector RealNum) -- ^ Activations
  deriving (Generic)

instance Show (FullyConnected' i o) where
  show FullyConnectedBLAS{}    = "FullyConnectedBLAS"
  show FullyConnectedHMatrix{} = "FullyConnectedHMatrix"

instance NFData (FullyConnected' i o) where
  rnf (FullyConnectedHMatrix b w)    = rnf b `seq` rnf w
  rnf (FullyConnectedBLAS !io !b !w) = rnf io `seq` rnf b `seq` rnf w

instance (KnownNat i, KnownNat o, KnownNat (i * o)) => UpdateLayer (FullyConnected i o) where
  type Gradient (FullyConnected i o) = (FullyConnected' i o)
  type MomentumStore (FullyConnected i o) = ListStore (FullyConnected' i o)
  runUpdate opt@OptSGD {} x@(FullyConnected (FullyConnectedHMatrix oldBias oldActivations) store) (FullyConnectedHMatrix biasGradient activationGradient) =
    let (FullyConnectedHMatrix oldBiasMomentum oldMomentum) = getData opt x store
        VectorResultSGD newBias newBiasMomentum = descendVector opt (VectorValuesSGD oldBias biasGradient oldBiasMomentum)
        MatrixResultSGD newActivations newMomentum = descendMatrix opt (MatrixValuesSGD oldActivations activationGradient oldMomentum)
        newStore = setData opt x store (FullyConnectedHMatrix newBiasMomentum newMomentum)
     in FullyConnected (FullyConnectedHMatrix newBias newActivations) newStore
  runUpdate opt@OptAdam {} x@(FullyConnected (FullyConnectedHMatrix oldBias oldActivations) store) (FullyConnectedHMatrix biasGradient activationGradient) =
    let [FullyConnectedHMatrix oldMBias oldMActivations, FullyConnectedHMatrix oldVBias oldVActivations] = getData opt x store
        VectorResultAdam newBias newMBias newVBias = descendVector opt (VectorValuesAdam (getStep store) oldBias biasGradient oldMBias oldVBias)
        MatrixResultAdam newActivations newMActivations newVActivations =
          descendMatrix opt (MatrixValuesAdam (getStep store) oldActivations activationGradient oldMActivations oldVActivations)
        newStore = setData opt x store [FullyConnectedHMatrix newMBias newMActivations, FullyConnectedHMatrix newVBias newVActivations]
     in FullyConnected (FullyConnectedHMatrix newBias newActivations) newStore
  runUpdate opt@OptSGD {} x@(FullyConnected (FullyConnectedBLAS io@(i,o) oldBiasH oldActivationsH) store) (FullyConnectedBLAS _ biasGradient activationGradient) =
      let oldBias = oldBiasH
          oldActivations = oldActivationsH
          (oldMBias, oldMActivations) = case getData opt x store of -- In the first periods until the store is filled newData is called, which will generate FullyConnectedHMatrix instances!
            FullyConnectedBLAS _ oldMBias' oldMActivations' -> (oldMBias', oldMActivations')
            FullyConnectedHMatrix oldMBias' oldMActivations' -> (extract oldMBias', extractM oldMActivations')
          VectorResultSGDV newBias newMBias               = descendVectorV opt (VectorValuesSGDV oldBias biasGradient oldMBias)
          MatrixResultSGDV newActivations newMActivations = descendMatrixV opt (MatrixValuesSGDV oldActivations activationGradient oldMActivations)
          newStore = setData opt x store (FullyConnectedBLAS io newMBias newMActivations)
      in FullyConnected (FullyConnectedBLAS io newBias newActivations) newStore
    where extractM mat = (\(S2DV vec) -> vec) . fromS2D $ S2D mat
  runUpdate opt@OptAdam {} x@(FullyConnected (FullyConnectedBLAS io@(i,o) oldBiasH oldActivationsH) store) (FullyConnectedBLAS _ biasGradient activationGradient) =
      let oldBias = oldBiasH
          oldActivations = oldActivationsH
          (oldMBias, oldMActivations, oldVBias, oldVActivations) = case getData opt x store of -- In the first periods until the store is filled newData is called, which will generate FullyConnectedHMatrix instances!
            [FullyConnectedBLAS _ oldMBias' oldMActivations', FullyConnectedBLAS _ oldVBias' oldVActivations'] -> (oldMBias', oldMActivations', oldVBias', oldVActivations')
            [FullyConnectedHMatrix oldMBias' oldMActivations', FullyConnectedHMatrix oldVBias' oldVActivations'] -> (extract oldMBias', extractM oldMActivations', extract oldVBias', extractM oldVActivations')
            [FullyConnectedBLAS _ oldMBias' oldMActivations', FullyConnectedHMatrix oldVBias' oldVActivations'] -> (oldMBias', oldMActivations', extract oldVBias', extractM oldVActivations')
            xs -> error $ "unexpected data in ListStore in FullyConnected BLAS implementation: " ++ show xs
          VectorResultAdamV newBias newMBias newVBias                      = descendVectorV opt (VectorValuesAdamV (getStep store) oldBias biasGradient oldMBias oldVBias)
          MatrixResultAdamV newActivations newMActivations newVActivations = descendMatrixV opt (MatrixValuesAdamV (getStep store) oldActivations activationGradient oldMActivations oldVActivations)
          newStore = setData opt x store [FullyConnectedBLAS io newMBias newMActivations, FullyConnectedBLAS io newVBias newVActivations]
      in FullyConnected (FullyConnectedBLAS io newBias newActivations) newStore
    where extractM mat = (\(S2DV vec) -> vec) . fromS2D $ S2D mat
  runUpdate opt (FullyConnected layer _) _ = error $ "Unexpected input in runUpdate in FullyConnected layer. Optimizer" ++ show opt ++ ". Layer: " ++ show layer

instance (KnownNat i, KnownNat o, KnownNat (i * o)) => LayerOptimizerData (FullyConnected i o) (Optimizer 'SGD) where
  type MomentumDataType (FullyConnected i o) (Optimizer 'SGD) = FullyConnected' i o
  getData opt x store = head $ getListStore opt x store
  setData opt x store = setListStore opt x store . return
  newData _ _ = FullyConnectedHMatrix (konst 0) (konst 0)

instance (KnownNat i, KnownNat o, KnownNat (i * o)) => LayerOptimizerData (FullyConnected i o) (Optimizer 'Adam) where
  type MomentumDataType (FullyConnected i o) (Optimizer 'Adam) = FullyConnected' i o
  type MomentumExpOptResult (FullyConnected i o) (Optimizer 'Adam) = [FullyConnected' i o]
  getData = getListStore
  setData = setListStore
  newData _ _ = FullyConnectedHMatrix (konst 0) (konst 0)

instance (KnownNat i, KnownNat o) => FoldableGradient (FullyConnected' i o) where
  mapGradient f (FullyConnectedHMatrix bias activations) = FullyConnectedHMatrix (dvmap f bias) (dmmap f activations)
  mapGradient f (FullyConnectedBLAS io bias activations) = FullyConnectedBLAS io (mapVector f bias) (mapVector f activations)
  squaredSums (FullyConnectedHMatrix bias activations) = [sumV . squareV $ bias, sumM . squareM $ activations]
  squaredSums (FullyConnectedBLAS _ bias activations) = [V.sum . mapVector (^(2::Int)) $ bias, V.sum . mapVector (^(2::Int)) $ activations]

runForward :: forall i o. (KnownNat i, KnownNat o) => FullyConnected i o -> S ('D1 i) -> (Tape (FullyConnected i o) ('D1 i) ('D1 o), S ('D1 o))
runForward (FullyConnected (FullyConnectedHMatrix wB wN) _) (S1D v) = (S1D v, S1D (wB + wN #> v))
runForward (FullyConnected (FullyConnectedBLAS _ wB wN) _) (S1DV v) =
  let inp = unsafeMemCopyVectorFromTo wB (createVectorUnsafe (V.length wB))
      !out' = unsafePerformIO $ matXVec BlasNoTranspose wN v 1 inp
   in force out' `seq` (S1DV v, S1DV out')
runForward lay@(FullyConnected (FullyConnectedHMatrix _ _) _ ) v@S1DV{} = runForward lay (toS1D v)
runForward lay@(FullyConnected FullyConnectedBLAS{} _) v@S1D{} = runForward lay (fromS1D v)

runBackward :: forall i o . (KnownNat i, KnownNat o) => FullyConnected i o -> Tape (FullyConnected i o) ('D1 i) ('D1 o) -> S ('D1 o) -> (Gradient (FullyConnected i o), S ('D1 i))
runBackward (FullyConnected (FullyConnectedHMatrix _ wN) _) (S1D x) (S1D dEdy) =
  let wB' = dEdy
      mm' = dEdy `outer` x
            -- calcluate derivatives for next step
      dWs = tr wN #> dEdy
   in (FullyConnectedHMatrix wB' mm', S1D dWs)
runBackward (FullyConnected (FullyConnectedBLAS io@(i, o) _ wN) _) (S1DV x) (S1DV dEdy) =
  let !mm' = unsafePerformIO $ outerV dEdy x
      dWsInp = createVectorUnsafe i
      !dWs' = unsafePerformIO $ matXVec BlasTranspose wN dEdy 0 dWsInp
   in force mm' `seq` force dWs' `seq` (FullyConnectedBLAS io dEdy mm', S1DV dWs')
runBackward l x dEdy = runBackward l x (toLayerShape x dEdy)

instance (KnownNat i, KnownNat o, KnownNat (i * o)) => Layer (FullyConnected i o) ('D1 i) ('D1 o) where
  type Tape (FullyConnected i o) ('D1 i) ('D1 o) = S ('D1 i)
  runForwards = runForward   -- Do a matrix vector multiplication and return the result.
  runBackwards = runBackward -- Run a backpropogation step for a full connected layer.

instance (KnownNat i, KnownNat o) => Serialize (FullyConnected i o) where
  put (FullyConnected w ms) = put w >> put ms
  get = FullyConnected <$> get <*> get

instance (KnownNat i, KnownNat o) => Serialize (FullyConnected' i o) where
  put (FullyConnectedHMatrix b w) = do
    put (0 :: Int)
    putListOf put . LA.toList . extract $ b
    putListOf put . LA.toList . LA.flatten . extract $ w
  put (FullyConnectedBLAS io b w) = do
    put (1 :: Int)
    put io
    putListOf put . V.toList $ b
    putListOf put . V.toList $ w
  get = do
    (nr :: Int) <- get
    case nr of
      0 -> do
        let f = fromIntegral $ natVal (Proxy :: Proxy i)
        b <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get
        k <- maybe (fail "Vector of incorrect size") return . create . LA.reshape f . LA.fromList =<< getListOf get
        return $ FullyConnectedHMatrix b k
      1 -> do
        io <- get
        b <- V.fromList <$> getListOf get
        w <- V.fromList <$> getListOf get
        return $ FullyConnectedBLAS io b w
      _ -> error $ "Unexpected nr in get in Serialize of FullyConnected' " ++ show nr

instance (KnownNat i, KnownNat o, KnownNat (i*o)) => RandomLayer (FullyConnected i o) where
  createRandomWith = randomFullyConnected

randomFullyConnected ::
     forall m i o. (PrimBase m, KnownNat i, KnownNat o, KnownNat (i * o))
  => NetworkInitSettings
  -> Gen (PrimState m)
  -> m (FullyConnected i o)
randomFullyConnected (NetworkInitSettings m HMatrix _) gen = do
  wN <- getRandomMatrix i o m gen
  wB <- getRandomVector i o m gen
  return $!! FullyConnected (FullyConnectedHMatrix wB wN) mkListStore
  where i = natVal (Proxy :: Proxy i)
        o = natVal (Proxy :: Proxy o)
randomFullyConnected (NetworkInitSettings m BLAS _) gen = do

  wB <- getRandomVectorV i o o' m gen
  wN <- getRandomVectorV i o (i' * o') m gen
  return $!! FullyConnected (FullyConnectedBLAS (i', o') wB wN) mkListStore
  where
    i = natVal (Proxy :: Proxy i)
    i' = fromIntegral i
    o = natVal (Proxy :: Proxy o)
    o' = fromIntegral o

-------------------- DynamicNetwork instance --------------------

instance (KnownNat i, KnownNat o) => FromDynamicLayer (FullyConnected i o) where
  fromDynamicLayer _ _ _ = SpecNetLayer $ SpecFullyConnected (natVal (Proxy :: Proxy i)) (natVal (Proxy :: Proxy o))

instance ToDynamicLayer SpecFullyConnected where
  toDynamicLayer wInit gen (SpecFullyConnected nrI nrO) =
    reifyNat nrI $ \(pxInp :: (KnownNat i) => Proxy i) ->
      reifyNat nrO $ \(pxOut :: (KnownNat o') => Proxy o') ->
        case singByProxy pxInp %* singByProxy pxOut of
          SNat -> do
            (layer :: FullyConnected i o') <- randomFullyConnected wInit gen
            return $ SpecLayer layer (sing :: Sing ('D1 i)) (sing :: Sing ('D1 o'))

-- | Make a specification of a fully connected layer (see Grenade.Dynamic.Build for a user-interface to specifications).
specFullyConnected :: Integer -> Integer -> SpecNet
specFullyConnected nrI nrO = SpecNetLayer $ SpecFullyConnected nrI nrO

-- | A Fully-connected layer with input dimensions as given in last output layer and output dimensions specified. 1D only!
fullyConnected :: Integer -> BuildM ()
fullyConnected rows = do
  (inRows, _, _) <- buildRequireLastLayerOut Is1D
  buildAddSpec (SpecNetLayer $ SpecFullyConnected inRows rows)
  buildSetLastLayer (rows, 1, 1)

-------------------- GNum instances --------------------

instance (KnownNat i, KnownNat o) => GNum (FullyConnected i o) where
  s |* FullyConnected w store = FullyConnected (s |* w) (s |* store)
  FullyConnected w1 store1 |+ FullyConnected w2 store2 = FullyConnected (w1 |+ w2) (store1 |+ store2)

instance (KnownNat i, KnownNat o) => GNum (FullyConnected' i o) where
  s |* FullyConnectedHMatrix b w = FullyConnectedHMatrix (dvmap (fromRational s *) b) (dmmap (fromRational s *) w)
  s |* FullyConnectedBLAS io b w = FullyConnectedBLAS io (mapVector (fromRational s *) b) (mapVector (fromRational s *) w)
  FullyConnectedHMatrix b1 w1 |+ FullyConnectedHMatrix b2 w2 = FullyConnectedHMatrix (b1 + b2) (w1 + w2)
  FullyConnectedBLAS io b1 w1 |+ FullyConnectedBLAS _ b2 w2 = FullyConnectedBLAS io (zipWithVector (+) b2 b1) (zipWithVector (+) w2 w1)
  x |+ y = error $ "Cannot add different network types in |+ in FullyConnected: " ++ show (x, y)
  sumG xs@(FullyConnectedBLAS io _ _:_) = FullyConnectedBLAS io bs' ws'
    where
      (bs, ws) = unzip $ map (\(FullyConnectedBLAS _ b w) -> (b, w)) xs
      bs' = sumVectors bs `using` rparWith rdeepseq
      ws' = sumVectors ws `using` rparWith rdeepseq
  sumG xs = foldl1 (|+) xs

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE OverloadedStrings #-}

module StatReport
  ( statReport
  , showStatEntryValue
  )
where

import           Data.Fixed                     ( showFixed )
import           Data.Text                      ( Text )
import           Fmt
import           QuoteData
import           Statistics

instance Buildable Statistic where
  build Mean = "Mean"
  build Min  = "Minimum"
  build Max  = "Maximum"
  build Days = "Days between Min/Max"

showStatEntryValue :: StatEntry -> String
showStatEntryValue StatEntry {..} = showFixed (removeTrailing stat qfield)
                                              value
 where
  removeTrailing Days _      = True
  removeTrailing Min  Volume = True
  removeTrailing Max  Volume = True
  removeTrailing _    _      = False

instance Buildable StatEntry where
  build se@StatEntry {..} = "" +| stat |+ ": " +| showStatEntryValue se |+ ""

instance Buildable StatQFieldData where
  build (qf, stats) = nameF ("Statistics for " +|| qf ||+ "") $ unlinesF stats

statReport :: StatInfo -> Text
statReport = fmt . unlinesF


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
module Math.HiddenMarkovModel (
   T(..), Distr.State, state,
   Discrete, DiscreteTrained,
   Gaussian, GaussianTrained,
   uniform,
   generate,
   probabilitySequence,
   Normalized.logLikelihood,
   Normalized.reveal,
   Normalized.posterior,

   Trained(..),
   trainSupervised,
   Normalized.trainUnsupervised,
   mergeTrained, finishTraining, trainMany,
   deviation,

   toCSV,
   fromCSV,
   ) where

import qualified Math.HiddenMarkovModel.Distribution as Distr
import qualified Math.HiddenMarkovModel.Normalized as Normalized
import qualified Math.HiddenMarkovModel.CSV as HMMCSV
import Math.HiddenMarkovModel.Private
          (T(..), Trained(..), mergeTrained, toCells, parseCSV)
import Math.HiddenMarkovModel.Distribution (State(State))
import Math.HiddenMarkovModel.Utility
          (randomItemProp, normalizeProb, attachOnes)

import qualified Numeric.LinearAlgebra.Algorithms as Algo
import qualified Numeric.Container as NC
import qualified Data.Packed.Matrix as Matrix
import qualified Data.Packed.Vector as Vector
import Data.Packed.Matrix (Matrix)
import Data.Packed.Vector (Vector)

import qualified Text.CSV.Lazy.String as CSV

import qualified System.Random as Rnd

import qualified Control.Monad.Exception.Synchronous as ME
import qualified Control.Monad.Trans.State as MS
import qualified Control.Monad.HT as Monad

import qualified Data.NonEmpty as NonEmpty
import qualified Data.Array as Array
import Data.Traversable (Traversable, mapAccumL)
import Data.Foldable (Foldable)
import Data.Array (accumArray)

state :: Int -> State
state = State

type DiscreteTrained prob symbol = Trained (Distr.DiscreteTrained prob symbol) prob
type Discrete prob symbol = T (Distr.Discrete prob symbol) prob

type GaussianTrained a = Trained (Distr.GaussianTrained a) a
type Gaussian a = T (Distr.Gaussian a) a

{- |
Create a model with uniform probabilities
for initial vector and transition matrix
given a distribution for the emissions.
You can use this as a starting point for 'Normalized.trainUnsupervised'.
-}
uniform ::
   (Distr.Info distr, Distr.Probability distr ~ prob) =>
   distr -> T distr prob
uniform distr =
   let n = Distr.numberOfStates distr
       c = recip $ fromIntegral n
   in  Cons {
          initial = NC.constant c n,
          transition = NC.konst c (n,n),
          distribution = distr
       }

probabilitySequence ::
   (Traversable f, Distr.EmissionProb distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
   T distr prob -> f (State, emission) -> f prob
probabilitySequence hmm =
   snd
   .
   mapAccumL
      (\index (State s, e) ->
         (NC.atIndex (transition hmm) . flip (,) s,
          index s * Distr.emissionStateProb (distribution hmm) e (State s)))
      (NC.atIndex (initial hmm))

generate ::
   (Rnd.RandomGen g, Ord prob, Rnd.Random prob,
    Distr.Generate distr, Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
   T distr prob -> g -> [emission]
generate hmm =
   MS.evalState $
   flip MS.evalStateT (initial hmm) $
   Monad.repeat $ MS.StateT $ \v0 -> do
      s <- randomItemProp $ zip [0..] (Vector.toList v0)
      x <- Distr.generate (distribution hmm) (State s)
      return (x, takeColumn s $ transition hmm)

takeColumn :: (Matrix.Element a) => Int -> Matrix a -> Vector a
takeColumn n  =  Matrix.flatten . Matrix.extractColumns [n]

{- |
Contribute a manually labeled emission sequence to a HMM training.
-}
trainSupervised ::
   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
   Int -> NonEmpty.T [] (State, emission) -> Trained tdistr prob
trainSupervised n xs =
   let getState (State s, _x) = s
   in  Trained {
          trainedInitial = NC.assoc n 0 [(getState (NonEmpty.head xs), 1)],
          trainedTransition =
             Matrix.trans $ NC.accum (NC.konst 0 (n,n)) (+) $
             attachOnes $ NonEmpty.mapAdjacent (,) $ fmap getState xs,
          trainedDistribution =
             Distr.accumulateEmissions $ map attachOnes $ Array.elems $
             accumArray (flip (:)) [] (State 0, State (n-1)) $ NonEmpty.flatten xs
       }

finishTraining ::
   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
    Distr.Probability distr ~ prob) =>
   Trained tdistr prob -> T distr prob
finishTraining hmm =
   Cons {
      initial = normalizeProb $ trainedInitial hmm,
      transition =
         Matrix.fromColumns $ map normalizeProb $
         Matrix.toColumns $ trainedTransition hmm,
      distribution = Distr.normalize $ trainedDistribution hmm
   }

trainMany ::
   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
    Distr.Probability distr ~ prob,
    Foldable f) =>
   (trainingData -> Trained tdistr prob) ->
   NonEmpty.T f trainingData -> T distr prob
trainMany train =
   finishTraining . NonEmpty.foldl1Map mergeTrained train

{- |
Compute maximum deviation between initial and transition probabilities.
You can use this as abort criterion for unsupervised training.
We omit computation of differences between the emission probabilities.
This simplifies matters a lot and
should suffice for defining an abort criterion.
-}
deviation ::
   (Algo.Field prob, Ord prob) => T distr prob -> T distr prob -> prob
deviation hmm0 hmm1 =
   deviationVec (initial hmm0) (initial hmm1)
   `max`
   deviationVec (transition hmm0) (transition hmm1)

deviationVec ::
   (Ord a, NC.Container c a) =>
   c a -> c a -> a
deviationVec x y =
   let d = NC.sub x y
   in  NC.maxElement d `max` negate (NC.minElement d)

toCSV ::
   (Distr.CSV distr, Algo.Field prob, Show prob) =>
   T distr prob -> String
toCSV hmm =
   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $
   toCells hmm

fromCSV ::
   (Distr.CSV distr, Algo.Field prob, Read prob) =>
   String -> ME.Exceptional String (T distr prob)
fromCSV =
   MS.evalStateT parseCSV . map HMMCSV.fixShortRow . CSV.parseCSV

max_stars_repo_name: rybern/hmm-hmatrix

max_issues_repo_name: rybern/hmm-hmatrix

max_forks_repo_name: rybern/hmm-hmatrix


text:
import qualified Prelude as P
import Control.Constrained.Prelude

import Control.Constrained.Applicative
import Control.Constrained.Category
import Control.Constrained.Comonad
import Control.Constrained.Foldable
import Control.Constrained.Functor
import Control.Constrained.Plain
import qualified Control.Monad as P
import Data.Aeson (genericToEncoding, defaultOptions)
import Data.Complex
import qualified Data.Monoid as M
import Data.Map (Map)
import qualified Data.Map as Map
import Data.Monoid.Instances ()
import qualified Data.Vector.Unboxed as U
import Data.Text (Text)
import Data.Yaml
import Formatting
import System.Directory

sum :: Foldable f => Dom f ~ (-#>) => Num a => PCon a => PCon (M.Sum a)
    => f a -> a
sum xs = M.getSum (foldMap (PFun M.Sum) xs)

sum2 :: Foldable f => Dom f ~ (-#>) => Num a => PCon a => PCon (M.Sum a)
     => f a -> a
sum2 xs = M.getSum (foldMap (PFun M.Sum . PFun \x -> x^(2::Int)) xs)

integral :: PCon a => Floating a => UIVector a -> a
integral xs = sum xs / fromIntegral (length xs)

norm2 :: PCon a => Floating a => UIVector a -> a
norm2 xs = rms (sum2 xs) (length xs)
  where rms :: Floating a => Integral i => a -> i -> a
        rms s2 c = sqrt (s2 / fromIntegral c)

-- physical boundaries
xmin :: Fractional a => a
xmin = 0
xmax :: Fractional a => a
xmax = 1

-- discretization
npts :: Int
npts = 100
hgsp :: Fractional a => a
hgsp = (xmax - xmin) / fromIntegral npts

cfl :: Fractional a => a
cfl = 0.5
nsteps :: Int
nsteps = 200

coord :: Fractional a => Int -> a
coord i = xavg + xrad / irad * (fromIntegral i - iavg)
  where -- discrete boundaries
        imin = -1/2
        imax = fromIntegral npts + 1/2
        -- auxiliary functions
        iavg = (imin + imax) / 2
        irad = (imax - imin) / 2
        xavg = (xmin + xmax) / 2
        xrad = (xmax - xmin) / 2

-- TODO: implement 'ToList UIVector', use it
coords :: PCon a => Fractional a => UIVector a
coords = uivector . U.fromListN npts . fmap coord $ [0..npts-1]

setup :: Floating a => a -#> Complex a
setup = PFun \x -> sin (2*pi*x) :+ 0

setups :: Functor f => Dom f ~ (-#>) => Cod f ~ (->)
       => PCon a => Floating a
       => f a -> f (Complex a)
setups = fmap setup

deriv :: ComonadStore f
      => Cod f ~ (->) => Eq (Index f) => Num (Index f)
      => Ok (Dom f) a => Fractional a
      => f a -> a
deriv xs = if | n == 0 -> 0
              | i == 0 -> (peek (i+1) xs - peek i xs) / hgsp
              | i == n-1 -> (peek i xs - peek (i-1) xs) / hgsp
              | otherwise -> (peek (i+1) xs - peek (i-1) xs) / (2 * hgsp)
  where i = pos xs
        n = fromIntegral npts

derivs :: ComonadStore f
       => Cod f ~ (->) => Eq (Index f) => Num (Index f)
       => Ok (Dom f) a => Fractional a
       => f a -> f a
derivs = extend deriv

energy :: ComonadStore f
       => Cod f ~ (->) => Eq (Index f) => Num (Index f)
       => Ok (Dom f) (Complex a) => RealFloat a
       => f (Complex a) -> a
energy us = (realPart du ^^ (2::Int) + imagPart du ^^ (2::Int)) / 2
   where du = deriv us

energies :: ComonadStore f
         => Cod f ~ (->) => Eq (Index f) => Num (Index f)
         => Ok (Dom f) a => Ok (Dom f) (Complex a) => RealFloat a
         => f (Complex a) -> f a
energies = extend energy

rhs :: ComonadStore f
    => Cod f ~ (->) => Eq (Index f) => Num (Index f)
    => Ok (Dom f) (Complex a) => RealFloat a
    => f (Complex a) -> Complex a
rhs us = du
   where du = deriv us

rhss :: ComonadStore f
     => Cod f ~ (->) => Eq (Index f) => Num (Index f)
     => Ok (Dom f) (Complex a) => RealFloat a
     => f (Complex a) -> f (Complex a)
rhss = extend rhs

rk2 :: Apply f
    => Dom f ~ (-#>) => Cod f ~ (->)
    => Ok (Dom f) a => Fractional a
    => (f a -> f a) -> a -> f a -> f a
rk2 f dx y0 = let k0 = f y0
                  y1 = liftA2u (PFun \(y, k) -> y + 0.5 * dx * k) y0 k0
                  k1 = f y1
                  y2 = liftA2u (PFun \(y, k) -> y + dx * k) y0 k1
              in y2

outputFilePath :: Int -> FilePath
outputFilePath i =
  formatToString ("output/wavetoy.it" % (left 6 '0' %. int) % ".yaml") i

instance (ToJSON a, PCon a) => ToJSON (Complex a) where
  toEncoding = genericToEncoding defaultOptions

instance (ToJSON a, PCon a) => ToJSON (UIVector a) where
  toJSON = toJSON . getUIVector
  toEncoding = toEncoding . getUIVector

data JSONData where
  JD :: ToJSON a => a -> JSONData

instance ToJSON JSONData where
  toJSON (JD x) = toJSON x
  toEncoding (JD x) = toEncoding x

output :: FilePath -> Map Text JSONData -> IO ()
output = encodeFile

main:: IO ()
main = do putStrLn "Constrained functors: Numeric functions"
          let x = coords @Double
          putStrLn $ "L2[x] = " ++ show (norm2 x)
          let u0 = setups x
          putStrLn $ "L2[u0] = " ++ show (realPart (norm2 u0))
          let e0 = energies u0
          putStrLn $ "E0 = " ++ show (integral e0)
          -- let du = rhss u0
          -- putStrLn $ "L2[du] = " ++ show (realPart (norm2 du))
          -- let u1 = rk2 rhss (cfl * hgsp) u0
          -- putStrLn $ "L2[u1] = " ++ show (realPart (norm2 u1))
          let step :: PCon a => RealFloat a
                   => UIVector (Complex a) -> UIVector (Complex a)
              step = rk2 rhss (cfl * hgsp)
          let uis = take (nsteps + 1) (P.iterate step u0)
          createDirectoryIfMissing True "output"
          P.forM (zip [0..] uis) \(i, u) ->
            do putStrLn $ "Iteration " ++ show i
               putStrLn $ "  L2[u] = " ++ show (realPart (norm2 u))
               let e = energies u
               putStrLn $ "  E = " ++ show (integral e)
               let vars = Map.fromList [ ("iteration", JD i)
                                       , ("x", JD x)
                                       , ("u", JD u0)
                                       , ("e", JD e0)
                                       ]
               P.when (i `mod` 10 == 0) $
                 output (outputFilePath i) vars
          putStrLn "Done."

max_stars_repo_name: eschnett/constrained-categories

max_issues_repo_name: eschnett/constrained-functors

max_forks_repo_name: eschnett/constrained-functors


text:
{-# LANGUAGE AllowAmbiguousTypes                      #-}
{-# LANGUAGE ApplicativeDo                            #-}
{-# LANGUAGE DeriveAnyClass                           #-}
{-# LANGUAGE DeriveDataTypeable                       #-}
{-# LANGUAGE DeriveGeneric                            #-}
{-# LANGUAGE DerivingVia                              #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE InstanceSigs                             #-}
{-# LANGUAGE KindSignatures                           #-}
{-# LANGUAGE MultiParamTypeClasses                    #-}
{-# LANGUAGE PatternSynonyms                          #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE StandaloneDeriving                       #-}
{-# LANGUAGE TemplateHaskell                          #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module Backprop.Learn.Model.Regression (
  -- * Linear Regression
    linReg, logReg
  , LRp(..), lrAlpha, lrBeta, runLRp
  -- * Reshape
  , reshapeLRpInput, reshapeLRpOutput
  , expandLRpInput, expandLRpOutput
  , premuteLRpInput, permuteLRpOutput
  -- * ARIMA
  , arima, autoregressive, movingAverage, arma
  , ARIMAp(..), ARIMAs(..)
  , arimaPhi, arimaTheta, arimaConstant, arimaYPred, arimaYHist, arimaEHist
  ) where

import           Backprop.Learn.Initialize
import           Backprop.Learn.Model.Combinator
import           Backprop.Learn.Model.Function
import           Backprop.Learn.Model.Types
import           Backprop.Learn.Regularize
import           Control.DeepSeq
import           Control.Monad.Primitive
import           Data.Finite
import           Data.Kind
import           Data.List
import           Data.Maybe
import           Data.Proxy
import           Data.Type.Equality
import           Data.Typeable
import           GHC.Generics                          (Generic)
import           GHC.TypeLits.Compare
import           GHC.TypeNats
import           Lens.Micro.TH
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import           Numeric.LinearAlgebra.Static.Vector
import           Numeric.OneLiner
import           Numeric.Opto.Ref
import           Numeric.Opto.Update hiding            ((<.>))
import           Statistics.Distribution
import           Unsafe.Coerce
import qualified Data.Binary                           as Bi
import qualified Data.Vector.Generic.Sized             as SVG
import qualified Data.Vector.Sized                     as SV
import qualified Data.Vector.Storable.Sized            as SVS
import qualified Numeric.LinearAlgebra                 as HU
import qualified Numeric.LinearAlgebra.Static          as H
import qualified System.Random.MWC                     as MWC

-- | Linear Regression parameter
data LRp (i :: Nat) (o :: Nat) = LRp
    { _lrAlpha :: !(R o)
    , _lrBeta  :: !(L o i)
    }
  deriving stock     (Generic, Typeable, Show)
  deriving anyclass  (NFData, Linear Double, Metric Double, Bi.Binary, Initialize, Backprop)

deriving via (GNum (LRp i o)) instance (KnownNat i, KnownNat o) => Num (LRp i o)
deriving via (GNum (LRp i o)) instance (KnownNat i, KnownNat o) => Fractional (LRp i o)
deriving via (GNum (LRp i o)) instance (KnownNat i, KnownNat o) => Floating (LRp i o)

makeLenses ''LRp

instance (PrimMonad m, KnownNat i, KnownNat o) => Mutable m (LRp i o) where
    type Ref m (LRp i o) = GRef m (LRp i o)
    thawRef   = gThawRef
    freezeRef = gFreezeRef
    copyRef   = gCopyRef
instance (PrimMonad m, KnownNat i, KnownNat o) => LinearInPlace m Double (LRp i o)

instance (KnownNat i, KnownNat o) => Regularize (LRp i o) where
    rnorm_1 = rnorm_1 . _lrBeta
    rnorm_2 = rnorm_2 . _lrBeta
    lasso r LRp{..} = LRp { _lrAlpha = 0
                          , _lrBeta  = lasso r _lrBeta
                          }
    ridge r LRp{..} = LRp { _lrAlpha = 0
                          , _lrBeta  = ridge r _lrBeta
                          }

instance (KnownNat i, KnownNat o, PrimMonad m) => Learnable m (LRp i o)

runLRp
    :: (KnownNat i, KnownNat o, Reifies s W)
    => BVar s (LRp i o)
    -> BVar s (R i)
    -> BVar s (R o)
runLRp lrp x = (lrp ^^. lrBeta) #> x + (lrp ^^. lrAlpha)

-- | Reshape an 'LRp' to take more or less inputs.  If more, new parameters
-- are initialized randomly according to the given distribution.
reshapeLRpInput
    :: (ContGen d, PrimMonad m, KnownNat i, KnownNat i', KnownNat o)
    => d
    -> MWC.Gen (PrimState m)
    -> LRp i o
    -> m (LRp i' o)
reshapeLRpInput d g (LRp α β) =
    LRp α <$> reshapeLCols d g β

-- | Reshape an 'LRp' to return more or less outputs  If more, new
-- parameters are initialized randomly according to the given distribution.
reshapeLRpOutput
    :: (ContGen d, PrimMonad m, KnownNat i, KnownNat o, KnownNat o')
    => d
    -> MWC.Gen (PrimState m)
    -> LRp i o
    -> m (LRp i o')
reshapeLRpOutput d g (LRp α β) =
    LRp <$> reshapeR d g α
        <*> reshapeLRows d g β

linReg
    :: (KnownNat i, KnownNat o)
    => Model ('Just (LRp i o)) 'Nothing (R i) (R o)
linReg = modelStatelessD (\(PJust p) -> runLRp p)

logReg
    :: (KnownNat i, KnownNat o)
    => Model ('Just (LRp i o)) 'Nothing (R i) (R o)
logReg = funcD logistic <~ linReg

-- | Adjust an 'LRp' to take extra inputs, initialized randomly.
--
-- Initial contributions to each output is randomized.
expandLRpInput
    :: (PrimMonad m, ContGen d, KnownNat i, KnownNat j, KnownNat o)
    => LRp i o
    -> d
    -> MWC.Gen (PrimState m)
    -> m (LRp (i + j) o)
expandLRpInput LRp{..} d g = LRp _lrAlpha . (_lrBeta H.|||) <$> initialize d g

-- | Adjust an 'LRp' to return extra ouputs, initialized randomly
expandLRpOutput
    :: (PrimMonad m, ContGen d, KnownNat i, KnownNat o, KnownNat p)
    => LRp i o
    -> d
    -> MWC.Gen (PrimState m)
    -> m (LRp i (o + p))
expandLRpOutput LRp{..} d g = do
    newAlpha <- initialize d g
    newBeta  <- initialize d g
    pure (LRp (_lrAlpha H.# newAlpha) (_lrBeta H.=== newBeta))

-- | Premute (or remove) inputs
--
-- Removed inputs will simply have their contributions removed from each
-- output.
premuteLRpInput
    :: (KnownNat i, KnownNat o)
    => SV.Vector i' (Finite i)
    -> LRp i o
    -> LRp i' o
premuteLRpInput is p = p { _lrBeta = colsL . fmap (β `SV.index`) $ is }
  where
    β = lCols (_lrBeta p)

-- | Premute (or remove) outputs
permuteLRpOutput
    :: (KnownNat i, KnownNat o)
    => SV.Vector o' (Finite o)
    -> LRp i o
    -> LRp i o'
permuteLRpOutput is LRp{..} =
    LRp { _lrAlpha = vecR . SVG.convert . fmap (α `SVS.index`) $ is
        , _lrBeta  = rowsL . fmap (β `SV.index`) $ is
        }
  where
    α = rVec _lrAlpha
    β = lRows _lrBeta

-- | 'ARIMA' parmaeters
data ARIMAp :: Nat -> Nat -> Type where
    ARIMAp :: { _arimaPhi      :: !(R p)
              , _arimaTheta    :: !(R q)
              , _arimaConstant :: !Double
              }
           -> ARIMAp p q
  deriving stock     (Generic, Typeable, Show)
  deriving anyclass  (NFData, Linear Double, Metric Double, Initialize, Backprop, Bi.Binary)

deriving via (GNum (ARIMAp p q)) instance Num (ARIMAp p q)
deriving via (GNum (ARIMAp p q)) instance Fractional (ARIMAp p q)
deriving via (GNum (ARIMAp p q)) instance Floating (ARIMAp p q)

makeLenses ''ARIMAp

-- | 'ARIMA' state
data ARIMAs :: Nat -> Nat -> Nat -> Type where
    ARIMAs :: { _arimaYPred :: !Double
              , _arimaYHist :: !(R (p + d))
              , _arimaEHist :: !(R q)
              }
          -> ARIMAs p d q
  deriving stock     (Generic, Typeable, Show)
  deriving anyclass  (NFData, Linear Double, Metric Double, Initialize, Backprop, Bi.Binary)

deriving via (NoRegularize (ARIMAs p d q)) instance Regularize (ARIMAs p d q)
deriving via (GNum (ARIMAs p d q)) instance Num (ARIMAs p d q)
deriving via (GNum (ARIMAs p d q)) instance Fractional (ARIMAs p d q)
deriving via (GNum (ARIMAs p d q)) instance Floating (ARIMAs p d q)

makeLenses ''ARIMAs

arima
    :: forall p d q. (KnownNat p, KnownNat d, KnownNat q)
    => Model ('Just (ARIMAp p q)) ('Just (ARIMAs p d q)) Double Double
arima = modelD $ \(PJust p) x (PJust s) ->
    let d :: L p (p + d)
        d = difference
        e  = x - (s ^^. arimaYPred)
        y  = (p ^^. arimaConstant)
           + (p ^^. arimaPhi  ) <.> (constVar d #> (s ^^. arimaYHist))
           + (p ^^. arimaTheta) <.> (s ^^. arimaEHist)
        yHist' = case Proxy @1 %<=? Proxy @(p + d) of
          LE Refl -> single y # constVar dropLast #> (s ^^. arimaYHist)
          NLE _ _ -> 0
        eHist' = case Proxy @1 %<=? Proxy @q of
          LE Refl -> single e # constVar dropLast #> (s ^^. arimaEHist)
          NLE _ _ -> 0
        s' = isoVar3 ARIMAs (\(ARIMAs pr yh eh) -> (pr,yh,eh))
                y
                yHist'
                eHist'
    in  (y, PJust s')

autoregressive
    :: KnownNat p
    => Model ('Just (ARIMAp p 0)) ('Just (ARIMAs p 0 0)) Double Double
autoregressive = arima

movingAverage
    :: KnownNat q
    => Model ('Just (ARIMAp 0 q)) ('Just (ARIMAs 0 0 q)) Double Double
movingAverage = arima

arma
    :: (KnownNat p, KnownNat q)
    => Model ('Just (ARIMAp p q)) ('Just (ARIMAs p 0 q)) Double Double
arma = arima

monosquare :: forall n. (n <=? (n ^ 2)) :~: 'True
monosquare = unsafeCoerce Refl

dropLast :: forall n. (KnownNat n, 1 <= n) => L (n - 1) n
dropLast = case monosquare @n of
    Refl -> vecL . SVS.generate $ \ij ->
      let i :: Finite n
          j :: Finite (n - 1)
          (i, j) = separateProduct ij
      in  if fromIntegral @_ @Int i == fromIntegral j
            then 1
            else 0

single :: Reifies s W => BVar s Double -> BVar s (R 1)
single = konst

difference'
    :: Int                  -- ^ initial
    -> Int                  -- ^ target
    -> HU.Matrix Double     -- ^ target x initial
difference' n m = foldl' go (HU.ident m) [m + 1 .. n]
  where
    go x k = x HU.<> d k
    d k = HU.build (k-1, k) $ \i j ->
        case round @_ @Int (j - i) of
          0 -> 1
          1 -> -1
          _ -> 0

difference :: forall n m. (KnownNat n, KnownNat m) => L n (n + m)
difference = fromJust . H.create $ difference' (n + m) n
  where
    n = fromIntegral $ natVal (Proxy @n)
    m = fromIntegral $ natVal (Proxy @m)

instance (PrimMonad m, KnownNat p, KnownNat q) => Mutable m (ARIMAp p q) where
    type Ref m (ARIMAp p q) = GRef m (ARIMAp p q)
    thawRef = gThawRef
    freezeRef = gFreezeRef
    copyRef = gCopyRef

instance (PrimMonad m, KnownNat p, KnownNat d, KnownNat q) => Mutable m (ARIMAs p d q) where
    type Ref m (ARIMAs p d q) = GRef m (ARIMAs p d q)
    thawRef = gThawRef
    freezeRef = gFreezeRef
    copyRef = gCopyRef

instance (KnownNat p, KnownNat q, PrimMonad m)  => LinearInPlace m Double (ARIMAp p q)
instance (KnownNat p, KnownNat d, KnownNat q, PrimMonad m) => LinearInPlace m Double (ARIMAs p d q)

instance (KnownNat p, KnownNat q) => Regularize (ARIMAp p q) where
    rnorm_1 ARIMAp{..} = rnorm_1 _arimaPhi + rnorm_1 _arimaTheta
    rnorm_2 ARIMAp{..} = rnorm_2 _arimaPhi + rnorm_2 _arimaTheta
    lasso r ARIMAp{..} = ARIMAp { _arimaPhi      = lasso r _arimaPhi
                                , _arimaTheta    = lasso r _arimaTheta
                                , _arimaConstant = 0
                                }
    ridge r ARIMAp{..} = ARIMAp { _arimaPhi      = ridge r _arimaPhi
                                , _arimaTheta    = ridge r _arimaTheta
                                , _arimaConstant = 0
                                }

instance (KnownNat p, KnownNat q, PrimMonad m) => Learnable m (ARIMAp p q)
instance (KnownNat p, KnownNat d, KnownNat q, PrimMonad m) => Learnable m (ARIMAs p d q)

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# OPTIONS_GHC -Wall #-}

-- | 'Distribution' avoids a name clash with 'Data.Distributive'
module Plankton.Distribution
  ( Distribution
  ) where

import Data.Complex (Complex(..))
import Plankton.Additive
import Plankton.Multiplicative
import Protolude (Bool(..), Double, Float, Int, Integer)

-- | Distribution (and annihilation) laws
--
-- > a * (b + c) == a * b + a * c
-- > (a + b) * c == a * c + b * c
-- > a * zero == zero
-- > zero * a == zero
class (Additive a, MultiplicativeMagma a) =>
      Distribution a

instance Distribution Double

instance Distribution Float

instance Distribution Int

instance Distribution Integer

instance Distribution Bool

instance (AdditiveGroup a, Distribution a) => Distribution (Complex a)

max_stars_repo_name: chessai/plankton

max_issues_repo_name: chessai/plankton

max_forks_repo_name: chessai/plankton


text:
{-#OPTIONS_GHC -w#-}
module Plots (plotAndSave2D') where

-- import Data.Colour
-- import Data.Colour.Names
import           Graphics.Rendering.Chart.Backend.Cairo
import           Graphics.Rendering.Chart.Easy
import           Numeric.LinearAlgebra                  as L (Vector (..), size,
                                                              toList)
import           System.FilePath                        (takeBaseName)

plotAndSave2D' :: FilePath
               -> Int
               -> Int
               -> (Int, Int)
               -> [L.Vector Double]
               -> IO ()
plotAndSave2D' file len basePoint (gn, pn) results =
  if (L.size . head $ results) /= 2
    then error "This function only works for 2D points!"
    else do
      putStrLn "Plotting..."
      toFile def plotFileName $ chartPlot2D results'
      putStrLn $ "Writing data points to file: " ++ dataFileName
      pp2DResults' dataFileName results'
  where
    actualFileName = takeBaseName file
    baseOutputName = actualFileName ++ "-"
                                       ++ "base" ++ show basePoint ++
                                                    "-" ++ show gn ++ "-" ++ show pn
    plotFileName = baseOutputName ++ ".png"
    dataFileName = baseOutputName ++ ".txt"
    results' = map ((\[a, b] -> (a, b)) . L.toList) results

chartPlot2D results = do
  layout_title .= "Normal Coordinates"
  setColors [opaque red]
  plot (points "points" results)

pp2DResults' :: FilePath -> [(Double,Double)] -> IO ()
pp2DResults' file xs = writeFile file . transPairLine $ xs
  where
    transPairLine = concatMap (\x -> show (fst x) ++ " " ++ show (snd x) ++ "\n")

max_stars_repo_name: emmanueldenloye/manifoldRNC

max_issues_repo_name: emmanueldenloye/manifoldRNC

max_forks_repo_name: emmanueldenloye/manifoldRNC


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE DeriveDataTypeable    #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TemplateHaskell       #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE ViewPatterns          #-}
-- |
-- Monoids for calculating various statistics in constant space
module Statistics.Monoid.Numeric (
    -- * Mean & Variance
    -- ** Number of elements
    CountG(..)
  , Count
  , asCount
    -- ** Mean
  , Mean
  , asMean
  , MeanKBN(..)
  , asMeanKBN
    -- ** Variance
  , Variance(..)
  , asVariance
    -- * Maximum and minimum
  , Max(..)
  , Min(..)
  , MaxD(..)
  , MinD(..)
    -- * Binomial trials
  , BinomAcc(..)
  , asBinomAcc
    -- * Accessors
  , CalcCount(..)
  , CalcMean(..)
  , CalcVariance(..)
  , calcStddev
  , calcStddevML
    -- * References
    -- $references
  ) where

import Data.Semigroup               (Semigroup(..))
import Data.Monoid                  (Monoid(..))
import Data.Data                    (Typeable,Data)
import Data.Vector.Unboxed          (Unbox)
import Data.Vector.Unboxed.Deriving (derivingUnbox)
import Numeric.Sum
import GHC.Generics                 (Generic)

import Statistics.Monoid.Class
-- COMPAT
import qualified Data.Vector.Generic         -- Needed for GHC7.4
import qualified Data.Vector.Generic.Mutable -- Needed for GHC7.4

----------------------------------------------------------------
-- Statistical monoids
----------------------------------------------------------------

-- | Calculate number of elements in the sample.
newtype CountG a = CountG { calcCountN :: a }
                  deriving (Show,Eq,Ord,Typeable)

type Count = CountG Int

-- | Type restricted 'id'
asCount :: CountG a -> CountG a
asCount = id

instance Integral a => Semigroup (CountG a) where
  (<>) = mappend
  {-# INLINE (<>) #-}

instance Integral a => Monoid (CountG a) where
  mempty                      = CountG 0
  CountG i `mappend` CountG j = CountG (i + j)
  {-# INLINE mempty  #-}
  {-# INLINE mappend #-}

instance (Integral a) => StatMonoid (CountG a) b where
  singletonMonoid _            = CountG 1
  addValue        (CountG n) _ = CountG (n + 1)
  {-# INLINE singletonMonoid #-}
  {-# INLINE addValue        #-}

instance CalcCount (CountG Int) where
  calcCount = calcCountN
  {-# INLINE calcCount #-}

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

-- | Type alias for currently recommended algorithms for calculation
--   of mean. It should be default choice
type Mean = MeanKBN

asMean :: Mean -> Mean
asMean = id

-- | Incremental calculation of mean. Sum of elements is calculated
--   using Kahan-Babuška-Neumaier summation.
data MeanKBN = MeanKBN !Int !KBNSum
             deriving (Show,Eq,Typeable,Data,Generic)

asMeanKBN :: MeanKBN -> MeanKBN
asMeanKBN = id

instance Semigroup MeanKBN where
  (<>) = mappend
  {-# INLINE (<>) #-}

instance Monoid MeanKBN where
  mempty = MeanKBN 0 mempty
  MeanKBN 0  _  `mappend` m             = m
  m             `mappend` MeanKBN 0  _  = m
  MeanKBN n1 s1 `mappend` MeanKBN n2 s2 = MeanKBN (n1+n2) (s1 `mappend` s2)

instance Real a => StatMonoid MeanKBN a where
  addValue (MeanKBN n m) x = MeanKBN (n+1) (addValue m x)

instance CalcCount MeanKBN where
  calcCount (MeanKBN n _) = n
instance CalcMean MeanKBN where
  calcMean (MeanKBN 0 _) = Nothing
  calcMean (MeanKBN n s) = Just (kbn s / fromIntegral n)

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

-- | This is algorithm for estimation of mean and variance of sample
--   which uses modified Welford algorithm. It uses KBN summation and
--   provides approximately 2 additional decimal digits
data VarWelfordKBN = VarWelfordKBN
  {-# UNPACK #-} !Int    --  Number of elements in the sample
  {-# UNPACK #-} !KBNSum -- Current sum of elements of sample
  {-# UNPACK #-} !KBNSum -- Current sum of squares of deviations from current mean

asVarWelfordKBN :: VarWelfordKBN -> VarWelfordKBN
asVarWelfordKBN = id

-- | Incremental algorithms for calculation the standard deviation.
data Variance = Variance {-# UNPACK #-} !Int    --  Number of elements in the sample
                         {-# UNPACK #-} !Double -- Current sum of elements of sample
                         {-# UNPACK #-} !Double -- Current sum of squares of deviations from current mean
                deriving (Show,Eq,Typeable)

-- | Type restricted 'id '
asVariance :: Variance -> Variance
asVariance = id
{-# INLINE asVariance #-}

instance Semigroup Variance where
  (<>) = mappend
  {-# INLINE (<>) #-}

-- | Iterative algorithm for calculation of variance [Chan1979]
instance Monoid Variance where
  mempty = Variance 0 0 0
  mappend (Variance n1 ta sa) (Variance n2 tb sb)
    = Variance (n1+n2) (ta+tb) sumsq
    where
      na = fromIntegral n1
      nb = fromIntegral n2
      nom = sqr (ta * nb - tb * na)
      sumsq | n1 == 0   = sb
            | n2 == 0   = sa
            | otherwise = sa + sb + nom / ((na + nb) * na * nb)
  {-# INLINE mempty #-}
  {-# INLINE mappend #-}

instance Real a => StatMonoid Variance a where
  addValue (Variance 0 _ _) x = singletonMonoid x
  addValue (Variance n t s) (realToFrac -> x)
    = Variance (n + 1) (t + x) (s + sqr (t  - n' * x) / (n' * (n'+1)))
    where
      n' = fromIntegral n
  singletonMonoid x = Variance 1 (realToFrac x) 0
  {-# INLINE singletonMonoid #-}

instance CalcCount Variance where
  calcCount (Variance n _ _) = n

instance CalcMean Variance where
  calcMean (Variance 0 _ _) = Nothing
  calcMean (Variance n s _) = Just (s / fromIntegral n)

instance CalcVariance Variance where
  calcVariance (Variance n _ s)
    | n < 2     = Nothing
    | otherwise = Just $! s / fromIntegral (n - 1)
  calcVarianceML (Variance n _ s)
    | n < 1     = Nothing
    | otherwise = Just $! s / fromIntegral n

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

-- | Calculate minimum of sample
newtype Min a = Min { calcMin :: Maybe a }
              deriving (Show,Eq,Ord,Typeable,Data,Generic)

instance Ord a => Semigroup (Min a) where
  (<>) = mappend
  {-# INLINE (<>) #-}

instance Ord a => Monoid (Min a) where
  mempty = Min Nothing
  Min (Just a) `mappend` Min (Just b) = Min (Just $! min a b)
  Min a        `mappend` Min Nothing  = Min a
  Min Nothing  `mappend` Min b        = Min b

instance (Ord a, a ~ a') => StatMonoid (Min a) a' where
  singletonMonoid a = Min (Just a)

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

-- | Calculate maximum of sample
newtype Max a = Max { calcMax :: Maybe a }
              deriving (Show,Eq,Ord,Typeable,Data,Generic)

instance Ord a => Semigroup (Max a) where
  (<>) = mappend
  {-# INLINE (<>) #-}

instance Ord a => Monoid (Max a) where
  mempty = Max Nothing
  Max (Just a) `mappend` Max (Just b) = Max (Just $! min a b)
  Max a        `mappend` Max Nothing  = Max a
  Max Nothing  `mappend` Max b        = Max b

instance (Ord a, a ~ a') => StatMonoid (Max a) a' where
  singletonMonoid a = Max (Just a)

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

-- | Calculate minimum of sample of Doubles. For empty sample returns NaN. Any
--   NaN encountered will be ignored.
newtype MinD = MinD { calcMinD :: Double }
              deriving (Show,Typeable,Data,Generic)

instance Eq MinD where
  MinD a == MinD b
    | isNaN a && isNaN b = True
    | otherwise          = a == b

instance Semigroup MinD where
  (<>) = mappend
  {-# INLINE (<>) #-}

-- N.B. forall (x :: Double) (x <= NaN) == False
instance Monoid MinD where
  mempty = MinD (0/0)
  mappend (MinD x) (MinD y)
    | isNaN x   = MinD y
    | isNaN y   = MinD x
    | otherwise = MinD (min x y)
  {-# INLINE mempty  #-}
  {-# INLINE mappend #-}

instance a ~ Double => StatMonoid MinD a where
  singletonMonoid = MinD

-- | Calculate maximum of sample. For empty sample returns NaN. Any
--   NaN encountered will be ignored.
newtype MaxD = MaxD { calcMaxD :: Double }
              deriving (Show,Typeable,Data,Generic)

instance Eq MaxD where
  MaxD a == MaxD b
    | isNaN a && isNaN b = True
    | otherwise          = a == b

instance Semigroup MaxD where
  (<>) = mappend
  {-# INLINE (<>) #-}

instance Monoid MaxD where
  mempty = MaxD (0/0)
  mappend (MaxD x) (MaxD y)
    | isNaN x   = MaxD y
    | isNaN y   = MaxD x
    | otherwise = MaxD (max x y)
  {-# INLINE mempty  #-}
  {-# INLINE mappend #-}

instance a ~ Double => StatMonoid MaxD a where
  singletonMonoid = MaxD

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

-- | Accumulator for binomial trials.
data BinomAcc = BinomAcc { binomAccSuccess :: !Int
                         , binomAccTotal   :: !Int
                         }
  deriving (Show,Eq,Ord,Typeable,Data,Generic)

-- | Type restricted 'id'
asBinomAcc :: BinomAcc -> BinomAcc
asBinomAcc = id

instance Semigroup BinomAcc where
  (<>) = mappend
  {-# INLINE (<>) #-}

instance Monoid BinomAcc where
  mempty = BinomAcc 0 0
  mappend (BinomAcc n1 m1) (BinomAcc n2 m2) = BinomAcc (n1+n2) (m1+m2)

instance StatMonoid BinomAcc Bool where
  addValue (BinomAcc nS nT) True  = BinomAcc (nS+1) (nT+1)
  addValue (BinomAcc nS nT) False = BinomAcc  nS    (nT+1)

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

sqr :: Double -> Double
sqr x = x * x
{-# INLINE sqr #-}

----------------------------------------------------------------
-- Unboxed instances
----------------------------------------------------------------

derivingUnbox "CountG"
  [t| forall a. Unbox a => CountG a -> a |]
  [| calcCountN |]
  [| CountG     |]

derivingUnbox "MeanKBN"
  [t| MeanKBN -> (Int,Double,Double) |]
  [| \(MeanKBN a (KBNSum b c)) -> (a,b,c)   |]
  [| \(a,b,c) -> MeanKBN a (KBNSum b c) |]

derivingUnbox "Variance"
  [t| Variance -> (Int,Double,Double) |]
  [| \(Variance a b c) -> (a,b,c)  |]
  [| \(a,b,c) -> Variance a b c    |]

derivingUnbox "MinD"
  [t| MinD -> Double |]
  [| calcMinD |]
  [| MinD     |]

derivingUnbox "MaxD"
  [t| MaxD -> Double |]
  [| calcMaxD |]
  [| MaxD     |]

-- $references
--
-- * [Welford1962] Welford, B.P. (1962) Note on a method for
--   calculating corrected sums of squares and
--   products. /Technometrics/
--   4(3):419-420. <http://www.jstor.org/stable/1266577>
--
-- * [Chan1979] Chan, Tony F.; Golub, Gene H.; LeVeque, Randall
--   J. (1979), Updating Formulae and a Pairwise Algorithm for
--   Computing Sample Variances., Technical Report STAN-CS-79-773,
--   Department of Computer Science, Stanford University. Page 4.


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE OverloadedStrings #-}
-- | Various tests
module Main where

import AI.Network.RNN.LSTM
import AI.Network.RNN.RNN
import AI.Network.RNN.Genetic
import AI.Network.RNN.Data
import AI.Network.RNN.Types
import AI.Network.RNN.Util
import AI.Network.RNN.Expr

import Test.Tasty
import Test.Tasty.HUnit
import Test.Tasty.QuickCheck

import qualified Numeric.LinearAlgebra.HMatrix as M

import Control.Monad.Random
import qualified Data.Text as T
import qualified Data.Set as DS

import Data.Char
import Data.List
import Numeric

import Debug.Trace

import Numeric.AD
import qualified Data.IntMap as I
import Text.PrettyPrint.HughesPJClass

main :: IO()
main = defaultMain tests

tests :: TestTree
tests = testGroup "Tests"
    [  testGroup "Utils" [
          testProperty "Matrix/vector product" prop_product
        , testCase "takes" $ do
            let l=[1,2,3,4,5,6]::[Int]
            [[1],[2,3],[4,5,6]] @=? takes ([1,2,3]::[Int]) l
            [[],[],[]] @=? takes [1,2,3] ([]::[Int])
        , testProperty "takes keep all data in order" prop_takes_concat
        , testProperty "takes has proper number of lists" prop_takes_length
        , testProperty "binary digits" prop_binary_digits
        , testCase "softmax" $ do
            let l=[1,2,3,4,5,6]::[Double]
                sf=softmax l
            sum sf @?= 1
        , testProperty "euclidian" prop_euclidian
        , testProperty "equilateral size" prop_equilateral_size
        , testProperty "equilateral distance" prop_equilateral_distance
        , testProperty "equilateral decoding" prop_equilateral_decoding
       ]
     , testGroup "RNN" [
        testCase "Check Steps without Back" $ checkSteps False
        , testCase "Check Steps with Back" $ checkSteps True
        , testCase "Check array conversion without Back" $ checkArray False
        , testCase "Check array conversion with Back" $ checkArray True
        , testCase "Check vector conversion without Back" $ checkArray False
        , testCase "Check vector conversion with Back" $ checkArray True
        ]
     , testGroup "Data" [
        testCase "Text" $ do
            testTextData "hello"
            testTextData "hello world!"
        , testCase "Text Sparse" $ do
            testTextDataS "hello"
            testTextDataS "hello world!"
     , testGroup "Genetic" [
            testCase "MixVector" $ do
                let v1 = M.fromList [1,1,1,1]
                let v2 = M.fromList [2,2,2,2]
                (v3,v4) <- evalRandIO $ mixVector v1 v2 0.5
                (v3 /= v4) @? "v3 == v4"
                M.size v1 @=? M.size v3
                M.size v1 @=? M.size v4
            , testCase "crossNetworkFull" $ testCrossover crossNetworkFull
            , testCase "crossNetworkHalf" $ testCrossover crossNetworkHalf
            , testCase "pointMutation" $ testMutation pointMutation
            , testCase "swapMutation" $ testMutation swapMutation
            , testCase "insertMutation" $ testMutation insertMutation
        ]
      , testGroup "LSTM" [
            testCase "Check steps" checkLSTMSteps
          , testProperty "Property steps" prop_lstm_eval_steps
          , testCase "Check vector conversion" checkLSTMVector
          , testProperty "list and vector match" prop_lstm_step
          , testCase "Check list of LSTMs" checkLSTMList
          , testProperty "Property list" prop_lstm_list
        ]
       , testGroup "LSTMIO" [
            testCase "LSTMIO basics" checkLSTMIO
        ]
       , testGroup "Expr" [
            testCase "Expr basics" testExpr
        ]
     ]
    ]

prop_product :: MatrixVector -> Bool
prop_product (MatrixVector sz ms vs) =
    let l=listMProd ms vs
        m=M.toList (M.matrix sz ms M.#> M.vector vs)
    in l == m -- (trace (show l) l)  == (trace (show m) m)

data MatrixVector = MatrixVector Int [Double] [Double]
    deriving (Show,Read,Eq,Ord)

instance Arbitrary MatrixVector where
    arbitrary = do
        Positive rows<-arbitrary
        Positive cols<-arbitrary
        ms <- vector (rows*cols)
        vs <- vector cols
        return $ MatrixVector cols ms vs

prop_takes_concat ::  String -> [Positive Int] -> Bool
prop_takes_concat xs ps = let
    idxs = map getPositive ps
    tot = sum idxs
    in concat (takes idxs xs) == take tot xs

prop_takes_length ::  String -> [Positive Int] -> Bool
prop_takes_length xs ps = let
    idxs = map getPositive ps
    in length (takes idxs xs) == length ps

prop_binary_digits :: Positive Int -> Bool
prop_binary_digits (Positive a) = binaryDigits a == length (showIntAtBase 2 intToDigit a "")

data EuclidianData = EuclidianData [Double] [Double]
    deriving Show

instance Arbitrary EuclidianData where
    arbitrary = do
        Positive sz<-arbitrary
        is1 <- vector sz
        is2 <- vector sz
        return $ EuclidianData is1 is2

prop_euclidian :: EuclidianData -> Bool
prop_euclidian (EuclidianData ds1 ds2) =
    euclidian (M.fromList ds1) (M.fromList ds2)
        == sum (map (\x->x*x) $ zipWith (-) ds2 ds1)

prop_equilateral_size :: Positive Int -> Bool
prop_equilateral_size (Positive n) =
    let m = equilateralEncoding n
    in (n,n-1) == (M.rows m,M.cols m)

prop_equilateral_distance :: Positive Int -> Bool
prop_equilateral_distance (Positive n) =
    let m = equilateralEncoding (n+1)
        ts = M.toRows m
        pairs = concatMap (\(x:ys)->map (\y->(x,y)) ys) $ init $ tails ts
        dis = map (roundTo 5 . uncurry euclidian) pairs
    in 1 == length (ordNub dis)

prop_equilateral_decoding :: Positive Int -> Bool
prop_equilateral_decoding (Positive n) =
    let rn = n + 1
        m  = equilateralEncoding rn
        vs= M.toRows m
        rs = map (equilateralDecoding m) vs
    in rs == [0..n]

testTextData :: T.Text -> IO()
testTextData t = do
    let (TrainData isSt is sz m) = textToTrainData t
    let charSet = T.foldl (flip DS.insert) DS.empty t
    DS.size charSet @=? sz+1
    t @=? dataToText m is
    T.init t @=? dataToText m (tail isSt)

testTextDataS :: T.Text -> IO()
testTextDataS t = do
    let (TrainData isSt is sz m) = textToTrainDataS t
    -- let charSet = T.foldl (flip DS.insert) DS.empty t
    -- sparseSize (DS.size charSet) @=? sz
    t @=? dataToTextS m is
    T.init t @=? dataToTextS m (tail isSt)

checkSteps :: Bool -> IO ()
checkSteps back = do
    (n::RNNetwork) <- evalRandIO $ randomNetwork (RNNDimensions 1 2 3 back) totalDataLength
    let (n2,out)=evalStep n $ M.fromList [1::Double]
        (n3,out1)=evalStep n2 $ M.fromList [3]
        (n4,out2)=evalSteps n [M.fromList [1],M.fromList [3]]
    n3 @=? n4
    out @=? head out2
    out1 @=? last out2

checkLSTMSteps :: IO ()
checkLSTMSteps = do
    (n::LSTMNetwork) <- evalRandIO $ randomNetwork 2 lstmFullSize
    let (n2,out)=evalStep n $ M.fromList [1::Double,2]
        (n3,out1)=evalStep n2 $ M.fromList [3,4]
        (n4,out2)=evalSteps n [M.fromList [1,2],M.fromList [3,4]]
    n3 @=? n4
    out @=? head out2
    out1 @=? last out2

prop_lstm_eval_steps :: LSTMData2 -> Bool
prop_lstm_eval_steps (LSTMData2 n1 is1 _ is2) =
    let
        (n2,out)=evalStep n1 $ M.fromList is1
        (n3,out1)=evalStep n2 $ M.fromList is2
        (n4,out2)=evalSteps n1 [M.fromList is1,M.fromList is2]
    in n3==n4 && out==head out2 && out1 == last out2

checkLSTMList :: IO ()
checkLSTMList = do
    (n1::LSTMNetwork) <- evalRandIO $ randomNetwork 2 lstmFullSize
    (n2::LSTMNetwork) <- evalRandIO $ randomNetwork 2 lstmFullSize
    let (n1_2,out)=evalStep n1 $ M.fromList [1::Double,2]
        (n2_2,out1)=evalStep n2 out
        (n3,out2)=evalStep [n1,n2] $ M.fromList [1,2]
    n3 @=? [n1_2,n2_2]
    out1 @=? out2

prop_lstm_list :: LSTMData2 -> Bool
prop_lstm_list (LSTMData2 n1 is1 n2 _) =
    let (n1_2,out)=evalStep n1 $ M.fromList is1
        (n2_2,out1)=evalStep n2 out
        (n3,out2)=evalStep [n1,n2] $ M.fromList is1
    in n3 == [n1_2,n2_2] && out1 == out2

checkLSTMVector :: IO ()
checkLSTMVector = do
    (n::LSTMNetwork) <- evalRandIO $ randomNetwork 2 lstmFullSize
    let arr = toVector n
        n2  = fromVector 2 arr
    n @=? n2

prop_lstm_step :: LSTMData -> Bool
prop_lstm_step (LSTMData n is)= snd (evalStep n (M.fromList is)) == M.fromList (snd $ lstmList 10 (M.toList $ toVector n) is)

checkLSTMIO :: IO()
checkLSTMIO = do
    print $ lstmFullSize 4
    print $ lstmioFullSize (2,4,3,1)
    (n1::LSTMIO) <- evalRandIO $ randomNetwork (2,4,3,1) lstmioFullSize
    let sz=rnnsize n1
        v1=toVector n1
        n2=fromVector sz v1
    n1 @=? n2
    let is=M.fromList [0.3,0.4]
        r1@(_,o1) = evalStep n1 is
    M.size o1 @?= 1
    let r2 = evalStep n2 is
    r1 @=? r2

data LSTMData = LSTMData LSTMNetwork [Double]
    deriving Show

instance Arbitrary LSTMData where
    arbitrary = do
        let sz=10
        ls <- vector (lstmFullSize sz)
        is <- vector sz
        let n = fromVector sz $ M.fromList ls
        return $ LSTMData n is

data LSTMData2 = LSTMData2 LSTMNetwork [Double] LSTMNetwork [Double]
    deriving Show

instance Arbitrary LSTMData2 where
    arbitrary = do
        let sz=10
        ls1 <- vector (lstmFullSize sz)
        is1 <- vector sz
        let n1 = fromVector sz $ M.fromList ls1
        ls2 <- vector (lstmFullSize sz)
        is2 <- vector sz
        let n2 = fromVector sz $ M.fromList ls2
        return $ LSTMData2 n1 is1 n2 is2

checkArray :: Bool -> IO ()
checkArray back = do
    let dim = RNNDimensions 1 2 3 back
    (n::RNNetwork) <- evalRandIO $ randomNetwork dim totalDataLength
    let arr = networkToArray n
        n2  = createNetworkFromArray dim arr
    case n2 of
        Right rn2 -> n @=? rn2
        Left err  -> assertFailure err

checkVector :: Bool -> IO ()
checkVector back = do
    let dim = RNNDimensions 1 2 3 back
    (n::RNNetwork) <- evalRandIO $ randomNetwork dim totalDataLength
    let arr = toVector n
        n2  = createNetworkFromVector dim arr
    case n2 of
        Right rn2 -> n @=? rn2
        Left err  -> assertFailure err

testCrossover :: (RNNetwork
                        -> RNNetwork -> Rand StdGen [RNNetwork])
                       -> IO ()
testCrossover f = do
    let dim = RNNDimensions 1 2 3 True
    (n1::RNNetwork) <- evalRandIO $ randomNetwork dim totalDataLength
    (n2::RNNetwork) <- evalRandIO $ randomNetwork dim totalDataLength
    rnns <- evalRandIO $ f n1 n2
    2 @=? length rnns
    notElem n1 rnns @? "n1 in result"
    notElem n2 rnns @? "n2 in result"

testMutation :: (LSTMNetwork -> Rand StdGen LSTMNetwork) -> IO()
testMutation f = do
    (n1::LSTMNetwork) <- evalRandIO $ randomNetwork 5 lstmFullSize
    n2 <- evalRandIO $ f n1
    toVector n1 /= toVector n2 @? "n1==n2"

testExpr :: IO()
testExpr = do
    let
        f :: (Num a,Floating a)=> [a] -> a
        f= \[x,y,_] -> x * sin (x + log y)
        i = I.fromList [(1,2::Double),(2,3),(3,4)]
        gs0 = grad f [Var 1,Var 2,Var 3]
        gs02 = grad f [2::Double,3,4]
        cexpr = map (\g-> close g i) gs0
        eexpr = map eval cexpr
    gs02 @=? eexpr
    let
        f1 :: (Num a,Floating a)=> a -> [a] -> a
        f1 = \x [y,_] -> x * sin (x + log y)
        i1 = I.fromList [(2,3),(3,4)]
        x1 :: Double
        x1 = 2
        gs1 = grad (f1 (autoEval (Lit $ Lit x1) $ I.empty)) [Var 2,Var 3]
        gs12 = grad (f1 2) [3,4]
        cexpr1 = map (\g-> fullSimplify $ close g i1) gs1
        eexpr1 = map eval cexpr1
    gs12 @=? eexpr1
    let
        f2 :: (Num a,Floating a)=> [a] -> [a] -> a
        f2= \[x] [y,_] -> x * sin (x + log y)
        i2 = I.fromList [(2,3),(3,4)]
        --ex :: [Expr (Expr Double)]
        --ex = map (\x->autoEval (Lit x) I.empty) [2::Double]
        x2 :: [Double]
        x2 = [2]
        gs2 = grad (f2 (map (\x->autoEval (Lit $ Lit x) I.empty) x2)) [Var 2,Var 3]
        gs22 = grad (f2 [2]) [3,4]
        cexpr2 = map (\g-> fullSimplify $ close g i2) gs2
        eexpr2 = map eval cexpr2
    gs22 @=? eexpr2
    let
        fs = lstmFullSize 2
        ixs = [1..fs]
        ixsF = map fromIntegral ixs
        gs3::[Expr Double] = grad costT $ map Var ixs
        i3 = I.fromList $ zip ixs ixsF
        gs32 = grad (costT) ixsF
        cexpr3 = map (\g-> fullSimplify $ close g i3) gs3
        eexpr3 = map eval cexpr3
    (map (rnd 6) gs32) @=? (map (rnd 6) eexpr3)
    where rnd n f = (fromInteger $ round $ f * (10^n)) / (10.0^^n)
--  print gs3
--  mapM_ (print . prettyShow . fullSimplify) (take 1 gs3)
--  print $ length gs3
--  print fs

costT :: (Num b,Floating b,Fractional b) => [b] -> b
costT lstm = let
    res = snd $ mapAccumL (lstmList 2) lstm [[1,2],[2,3]]
    in sum $ zipWith err [[3,4],[4,5]] res
    where
      err :: (Num b,Floating b) => [b] -> [b] -> b
      err a b  = sum (zipWith (\c d -> (c- d)**2 ) a b)

max_stars_repo_name: JPMoresmau/rnn

max_issues_repo_name: JPMoresmau/rnn

max_forks_repo_name: JPMoresmau/rnn


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards       #-}

module Q.Options.BlackScholes (
    BlackScholes(..)
  , atmf
  , euOption
  , eucall
  , euput
  , module Q.Options
) where

import           Control.Monad.State
import           Data.Random                    hiding (Gamma)
import           Data.Time
import           Numeric.RootFinding
import           Q.ContingentClaim.Options
import           Q.MonteCarlo
import           Q.Options
import           Q.Stochastic.Discretize
import           Q.Stochastic.Process
import           Q.Time
import           Q.Types
import           Statistics.Distribution        (cumulative, density)
import           Statistics.Distribution.Normal (standard)
import qualified Q.Options.Black76 as B76

dcf = dcYearFraction ThirtyUSA

-- | Parameters for a simplified black scholes equation.
data BlackScholes = BlackScholes {
    bsSpot :: Spot -- ^ The asset's spot on the valuation date.
  , bsRate :: Rate   -- ^ Risk free rate.
  , bsVol  :: Vol    -- ^ Volatility.
} deriving Show

instance Model BlackScholes Double where
  discountFactor BlackScholes{..} t1 t2 = return $ exp (scale dt bsRate)
    where dt = t2 - t1

  evolve (BlackScholes spot (Rate r) (Vol sigma)) (YearFrac t) = do
    (YearFrac t0, s0) <- get
    let dt = t - t0
    dw <- (lift stdNormal)::StateT (YearFrac, Double) RVar Double
    let st = s0 * exp ((r - 0.5 * sigma * sigma) * dt + sigma * dw * sqrt dt)
    put (YearFrac t, st)
    return st

atmf :: BlackScholes -> YearFrac -> Strike
atmf BlackScholes{..} t = Strike $ s / d where
  (Rate d) = exp (scale t (-bsRate))
  (Spot s) = bsSpot

-- | European option valuation with black scholes.
euOption ::  BlackScholes ->  YearFrac -> OptionType -> Strike ->Valuation
euOption bs@BlackScholes{..} t cp k =
  let b76 = B76.Black76 {
          b76F  = forward bs t
        , b76DF = Q.Types.discountFactor t bsRate
        , b76T  = t
        , b76Vol = bsVol
        }
  in B76.euOption b76 cp k

-- | see 'euOption'
euput bs t = euOption bs t Put

-- | see 'euOption'
eucall bs t = euOption bs t Call

forward BlackScholes{..} (YearFrac t) = Forward $ s * exp (r * t)
  where (Spot s) = bsSpot
        (Rate r) = bsRate

corradoMillerIniitalGuess bs@BlackScholes{..} cp (Strike k) (YearFrac t) (Premium premium) =
  (recip $ sqrt t) * ((sqrt (2 * pi)/ (s + discountedStrike)) + (premium - (s - discountedStrike)/2) + sqrt ((premium - (s - discountedStrike)/2)**2 - ((s - discountedStrike)**2/pi))) where
    discountedStrike = k * (exp $ (-r) * t)
    (Rate r) = bsRate
    (Spot s) = bsSpot

max_stars_repo_name: ghais/lowq

max_issues_repo_name: ghais/lowq

max_forks_repo_name: ghais/lowq


text:
{-# LANGUAGE BangPatterns     #-}
{-# LANGUAGE CPP              #-}
{-# LANGUAGE DataKinds        #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs            #-}
{-# LANGUAGE TemplateHaskell  #-}
{-# LANGUAGE TypeOperators    #-}

module Test.Grenade.Sys.Training where

-- TODO: These tests are unreliable and are not included.
-- In particular, there needs to be more investigation
-- into how to fix the seed that fit uses for random generation
-- or otherwise expose the combineTraining functions from the
-- Training module.

import qualified Numeric.LinearAlgebra.Data   as D
import qualified Numeric.LinearAlgebra.Static as H
import           GHC.TypeLits
import           Control.Monad.Random

import           Grenade

import           Test.Grenade.Sys.Utils

import           Hedgehog
import qualified Hedgehog.Gen                 as Gen
import qualified Hedgehog.Range               as Range

type FFNet
  = Network
    '[ FullyConnected 2 40, Tanh,
       FullyConnected 40 10, Relu,
       FullyConnected 10 1, Logit
     ]
    '[ 'D1 2,          -- Input
       'D1 40, 'D1 40, -- Fully connected, Tanh
       'D1 10, 'D1 10, -- Fully connected, Relu
       'D1 1, 'D1 1    -- Fully connected, Logit
     ]

circle :: S ('D1 2) -> S ('D1 1)
circle (S1D v) =
  if v `inCircle` (fromRational 0, 0.5)
     then S1D $ fromRational 1
     else S1D $ fromRational 0
  where
    inCircle :: KnownNat n => H.R n -> (H.R n, RealNum) -> Bool
    u `inCircle` (o, r) = H.norm_2 (u - o) <= r

prop_training_decreases_loss :: Property
prop_training_decreases_loss = withTests 10 $ property $ do
    let n = 1000
    rand <- forAll $ Gen.int $ Range.linear 1 10000
    evalIO $ setStdGen (mkStdGen rand)
    inps <- evalIO $ replicateM n randomOfShape
    let epochs = 10
        ts     = zip inps (map circle inps)
        trainOptsSgd = TrainingOptions { optimizer = defSGD
                                       , batchSize = 10
                                       , validationFreq = 1
                                       , verbose = Silent
                                       , metrics = []
                                       }
        trainOptsAdam = TrainingOptions { optimizer = defAdam
                                        , batchSize = 10
                                        , validationFreq = 1
                                        , verbose = Silent
                                        , metrics = []
                                        }
    losses  <- evalIO $ mapM (getLosses trainOptsSgd ts rand)  [1..epochs]
    losses' <- evalIO $ mapM (getLosses trainOptsAdam ts rand) [1..epochs]
    -- Perform linear regression on the log of the losses
    let m         = linearRegression epochs losses
        m'        = linearRegression epochs losses'
    -- We expect the loss to decrease over time
    assert $ m < 0
    assert $ m' < 0
  where
    getLosses opts trainData s i = do
      setStdGen (mkStdGen s)
      net <- fit trainData [] opts i quadratic' :: IO FFNet
      let total = foldr (\(x, y) a -> let out = runNet net x in a + (quadratic out y)) 0.0 trainData
          loss  = total / (fromIntegral $ length trainData)
      return loss

prop_training_improves_accuracy :: Property
prop_training_improves_accuracy = withTests 10 $ property $ do
  let n = 1000
  rand <- forAll $ Gen.int $ Range.linear 1 10000
  evalIO $ setStdGen (mkStdGen rand)
  inps <- evalIO $ replicateM n randomOfShape
  let epochs = 20
      ts     = zip inps (map circle inps)
      trainOptsSgd = TrainingOptions { optimizer = defSGD
                                     , batchSize = 10
                                     , validationFreq = 1
                                     , verbose = Silent
                                     , metrics = []
                                     }
      trainOptsAdam = TrainingOptions { optimizer = defAdam
                                      , batchSize = 10
                                      , validationFreq = 1
                                      , verbose = Silent
                                      , metrics = []
                                      }
  accuracies  <- evalIO (mapM (getAccuracy trainOptsSgd ts rand) [1..epochs]  :: IO [Int])
  accuracies' <- evalIO (mapM (getAccuracy trainOptsAdam ts rand) [1..epochs] :: IO [Int])
  -- Perform linear regression on the log of the accuracies
  let acc = map fromIntegral accuracies   :: [Double]
      acc' = map fromIntegral accuracies' :: [Double]
      m  = linearRegression epochs acc
      m' = linearRegression epochs acc'
  -- We expect the accuracy to increase over time
  assert $ m > 0
  assert $ m' > 0
  where
    getAccuracy opts trainData s i = do
      setStdGen (mkStdGen s)
      net <- fit trainData [] opts i quadratic' :: IO FFNet
      inps <- replicateM 100 randomOfShape
      let unshape :: S ('D1 1) -> RealNum
          unshape  = \(S1D x) -> (H.unwrap x) `D.atIndex` 0
          expected = map (unshape . circle) inps
          actual   = map (unshape . runNet net) inps
      let accuracy = length $ filter (<= 0.2) $ zipWith (\x y -> abs (x - y)) expected actual
      return accuracy

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:
import Data.Array.CArray
import Data.Complex
import Math.FFT (dft, idft) -- Binding to fftw

type Vector = CArray Int (Complex Double)

calculateEnergy :: Double -> Vector -> Vector -> Vector -> Double
calculateEnergy dx kin pot wfc = (* dx) . sum . map realPart $ elems total
  where
    total = liftArray2 (+) kineticE potentialE
    potentialE = wfcConj .* pot .* wfc
    kineticE = wfcConj .* idft (kin .* dft wfc)
    wfcConj = liftArray conjugate wfc
    a .* b = liftArray2 (*) a b

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
module Statistics.Classification.ConfusionMatrix (
    module ConfusionMatrix
  ) where

import Statistics.Classification.ConfusionMatrix.Binary as ConfusionMatrix
import Statistics.Classification.ConfusionMatrix.Multi as ConfusionMatrix

max_stars_repo_name: tsbattman/rochs

max_issues_repo_name: tsbattman/rochs

max_forks_repo_name: tsbattman/rochs


text:
import Statistics.Distribution
import Statistics.Distribution.Normal

nd :: NormalDistribution
nd = normalDistr mean sd
     where
       mean = 75
       sd   = 10

main = do
  putStrLn $ "P(80<=Z<=90)= " ++ show (cumulative nd 90 - cumulative nd 80)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat