text:
-- Copyright 2021 Google LLC.
-- SPDX-License-Identifier: Apache-2.0

-- #region language extensions
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DuplicateRecordFields #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE NumericUnderscores #-}
{-# LANGUAGE OverloadedLabels #-}
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE ViewPatterns #-}

-- #endregion

module Main where

-- #region imports

import Control.DeepSeq
import Control.Lens (at, each, folded, isn't, ix, traversed, (%~), (&), (.~), (?~), (^.), (^..), (^?), _1, _2, _3, _4)
import Control.Lens.Extras (is)
import Control.Monad (forM_)
import qualified Data.Array as A
import qualified Data.Array.Unboxed as UA
import qualified Data.ByteString as BS
import qualified Data.ByteString.Char8 as BS8
import qualified Data.ByteString.Lazy as LBS
import qualified Data.ByteString.Lazy.Char8 as LBS8
import Data.Complex (Complex (..), cis, imagPart, realPart)
import Data.Functor ((<&>))
import Data.Generics.Labels ()
import qualified Data.Graph.Inductive as G
import Data.Graph.Inductive.PatriciaTree (Gr)
import Data.List (elemIndex, scanl', sortOn, transpose, foldl', tails)
import Data.List.Split (chunksOf, splitOn)
import Data.Map.Strict (Map)
import qualified Data.Map.Strict as M
import Data.Maybe (catMaybes, fromJust, mapMaybe)
import Data.MemoTrie
import Data.Sequence (Seq (..))
import qualified Data.Sequence as Seq
import Data.Set (Set)
import qualified Data.Set as S
import qualified Data.Text as T
import qualified Data.Text.Lazy as LT
import Data.Vector (Vector)
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as UV
import qualified Data.Vector.Unboxed.Mutable as MUV
import Debug.Trace (trace, traceShow, traceShowId)
import GHC.Generics (Generic)
import GHC.IO (unsafePerformIO)
import Linear (V2 (..), angle, (*^), _x, _y)
import QuickSpec
import Safe hiding (at)
import Test.QuickCheck (Arbitrary (..), NonNegative (NonNegative), chooseInt, vectorOf, (==>), quickCheckResult)
import Test.QuickCheck.Modifiers (Positive (Positive, getPositive))
import qualified Text.Megaparsec as P
import qualified Text.Megaparsec.Char as P
import qualified Text.Megaparsec.Char.Lexer as L
import Text.Pretty.Simple (pPrint, pShow)
import Text.Regex.PCRE ((=~))
import Util

-- #endregion

type Input = UV.Vector Int

parse :: String -> Input
parse = UV.fromList . map (read . (: [])) . head . lines

input :: Input
input = parse (unsafePerformIO (readFile "/home/polux/projects/aoc/2019/day16-2/input.txt"))

vecSize :: Int
vecSize = UV.length input * 10_000

offset :: Int
offset = UV.take 7 input & UV.toList & concatMap show & read

inputFromOffset :: [Int]
inputFromOffset = [input UV.! (i `mod` UV.length input) | i <- [offset..vecSize-1]]

multFromOffset :: [Int] -> [Int]
multFromOffset vec = scanr (+) 0 vec & map (`mod` 10)

result :: [Int]
result = take 8 (iterate multFromOffset inputFromOffset !! 100)

main :: IO ()
main = print result

ext: hs

max_stars_repo_name: polux/advent-of-code

max_issues_repo_name: polux/advent-of-code

max_forks_repo_name: polux/advent-of-code


text:
module Scheme.Parser where

import Control.Monad (unless, void)
import Control.Monad.Except (MonadError (throwError))
import Data.Array (listArray)
import qualified Data.ByteString as BS
import qualified Data.Char as C
import Data.Complex (Complex ((:+)))
import Data.Ratio ((%))
import Data.Text (Text)
import qualified Data.Text as T
import Data.Text.Read (hexadecimal)
import Data.Word (Word8)
import Scheme.Types (Number (..), SchemeError (..), SchemeResult, SchemeVal (..))
import Text.Megaparsec
import Text.Megaparsec.Char
import qualified Text.Megaparsec.Char.Lexer as L

newtype ParserError = Unimplemented Text deriving (Show, Eq, Ord)

instance ShowErrorComponent ParserError where
  showErrorComponent (Unimplemented text) = T.unpack text <> " is not implemented"

type Parser = Parsec ParserError Text

sc :: Parser ()
sc =
  L.space
    space1
    (L.skipLineComment ";")
    (L.skipBlockCommentNested "#|" "|#")

lexeme :: Parser a -> Parser a
lexeme = L.lexeme sc

symbol :: Text -> Parser Text
symbol = L.symbol sc

parens :: Parser a -> Parser a
parens = between (symbol "(") (symbol ")")

brackets :: Parser a -> Parser a
brackets = between (symbol "[") (symbol "]")

braces :: Parser a -> Parser a
braces = between (symbol "{") (symbol "}")

pInteger :: Parser Integer
pInteger = L.signed (return ()) L.decimal

pBinaryInteger :: Parser Integer
pBinaryInteger = chunk "#b" >> L.binary

pOctalInteger :: Parser Integer
pOctalInteger = chunk "#o" >> L.octal

pDecimalInteger :: Parser Integer
pDecimalInteger = chunk "#d" >> L.decimal

pHexadecimalInteger :: Parser Integer
pHexadecimalInteger = chunk "#x" >> L.hexadecimal

integer :: Parser Number
integer =
  Integer
    <$> ( choice
            [ pInteger,
              pBinaryInteger,
              pOctalInteger,
              pDecimalInteger,
              pHexadecimalInteger
            ]
            <?> "integer"
        )

pReal :: Parser Double
pReal = L.signed (return ()) L.float

pDouble :: Parser Number
pDouble = Real <$> (pReal <?> "double")

pRational :: Parser Number
pRational = do
  numerator <- pInteger
  void (char '/')
  denominator <- pInteger
  -- TODO: Fix error if denominator is 0
  return $ Rational (numerator % denominator)

pComplex :: Parser Number
pComplex = do
  real <- pNan' <|> pInfinity' <|> try pReal <|> (fromInteger <$> pInteger)
  void (char '+')
  imag <- try pReal <|> fromInteger <$> pInteger <|> read "NaN" <$ chunk "nan.0" <|> read "Infinity" <$ chunk "inf.0"
  void (char 'i')
  return $ Complex (real :+ imag)

pInfinity' :: Parser Double
pInfinity' =
  read "Infinity" <$ chunk "+inf.0"
    <|> read "-Infinity" <$ chunk "-inf.0"

pInfinity :: Parser Number
pInfinity = Real <$> pInfinity'

pNan' :: Parser Double
pNan' = read "NaN" <$ chunk "+nan.0"

pNan :: Parser Number
pNan = Real <$> pNan'

pExact :: Parser Number
pExact = do
  void (try $ chunk "#e")
  num <- number'
  case num of
    x@(Integer _) -> return x
    x@(Rational _) -> return x
    (Real x) -> return $ Rational (toRational x)
    _ -> customFailure $ Unimplemented "Exactness for complex numbers"

pInexact :: Parser Number
pInexact = do
  void (try $ chunk "#i")
  num <- number'
  case num of
    x@(Real _) -> return x
    (Integer x) -> return $ Real (fromInteger x)
    (Rational x) -> return $ Real (fromRational x)
    _ -> customFailure $ Unimplemented "Exactness for complex numbers"

number' :: Parser Number
number' = choice (map (try . lexeme) [pComplex, pRational, pDouble, integer, pInfinity, pNan, pExact, pInexact])

number :: Parser SchemeVal
number = Number <$> number'

pSymbol :: Parser SchemeVal
pSymbol = try $ do
  ident <- lexeme (T.pack <$> start <> rest <?> "identifier")
  if ident == "."
    then empty
    else return $ Symbol ident
  where
    extendedSymbols = satisfy (`elem` ['!', '$', '%', '&', '*', '+', '-', '.', '/', ':', '<', '=', '>', '?', '@', '^', '_', '~'])
    start = some (letterChar <|> extendedSymbols)
    rest = many (alphaNumChar <|> extendedSymbols)

pString :: Parser SchemeVal
pString = lexeme $ String . T.pack <$> (char '"' *> manyTill L.charLiteral (char '"'))

pChar :: Parser SchemeVal
pChar = do
  void (try $ chunk "#\\")
  chr <- many asciiChar
  case chr of
    "alarm" -> return $ Character '\BEL'
    "backspace" -> return $ Character '\BS'
    "delete" -> return $ Character '\DEL'
    "escape" -> return $ Character '\ESC'
    "newline" -> return $ Character '\LF'
    "null" -> return $ Character '\NUL'
    "return" -> return $ Character '\CR'
    "space" -> return $ Character ' '
    "tab" -> return $ Character '\HT'
    _ -> pChar' chr

pChar' :: [Char] -> Parser SchemeVal
pChar' [c] = pure (Character c) <?> "char"
pChar' ('x' : hex) = case hexadecimal (T.pack hex) of
  Right n -> return (Character <$> C.chr $ fst n) <?> "hex char"
  _ -> empty
pChar' _ = empty

boolean :: Parser SchemeVal
boolean =
  choice $
    map
      lexeme
      [ Boolean True <$ chunk "#t",
        Boolean False <$ chunk "#f"
      ]

pPairList :: Parser SchemeVal
pPairList = try $
  parens $ do
    ls <- pExpr `sepEndBy` sc <?> "dotted pair car"
    dot <- optional (symbol "." >> pExpr) <?> "pair"
    pure $ case dot of
      Nothing -> List ls
      Just (List c) -> List (ls ++ c)
      Just (PairList car cdr) -> PairList (ls ++ car) cdr
      Just val -> PairList ls val

pList :: Parser SchemeVal
pList = List <$> try (parens (pExpr `sepEndBy` sc)) <?> "list"

pVector :: Parser SchemeVal
pVector = do
  void (try $ symbol "#(")
  ls <- pExpr `sepEndBy` sc
  void (symbol ")")
  return $ Vector (listArray (0, length ls - 1) ls)

pBytevector :: Parser SchemeVal
pBytevector = do
  void (try $ symbol "#u8(")
  ls <- pInteger `sepEndBy` sc
  void (symbol ")")
  return $ Bytevector (BS.pack $ map toByte ls)
  where
    toByte num = fromInteger num :: Word8

pQuote :: Parser SchemeVal
pQuote = do
  void (try $ lexeme $ char '\'')
  expr <- pExpr
  return $ List [Symbol "quote", expr]

pQuasiquote :: Parser SchemeVal
pQuasiquote = do
  void (try $ lexeme $ char '`')
  expr <- pExpr
  return $ List [Symbol "quasiquote", expr]

pUnquote :: Parser SchemeVal
pUnquote = do
  void (try $ lexeme $ char ',')
  expr <- pExpr
  return $ List [Symbol "unquote", expr]

pUnquoteSplicing :: Parser SchemeVal
pUnquoteSplicing = do
  void (try $ lexeme $ chunk ",@")
  expr <- pExpr
  return $ List [Symbol "unquote-splicing", expr]

pExpr :: Parser SchemeVal
pExpr =
  number
    <|> pChar
    <|> pSymbol
    <|> pString
    <|> boolean
    <|> pList
    <|> pPairList
    <|> pVector
    <|> pBytevector
    <|> pQuote
    <|> pQuasiquote
    <|> pUnquoteSplicing
    <|> pUnquote
    <?> "expression"

parseExpr :: FilePath -> Text -> Either SchemeError SchemeVal
parseExpr file input = case runParser pExpr file input of
  Right out -> Right out
  Left err -> Left $ ParserError (T.pack $ errorBundlePretty err)

parseInput :: Text -> SchemeResult SchemeVal
parseInput = parseExpr "file"

parseFile :: FilePath -> Text -> Either SchemeError SchemeVal
parseFile = parseExpr

parseOrThrow :: Parser a -> Text -> SchemeResult a
parseOrThrow parser input =
  let checkEmpty = do
        rest <- getInput
        unless (T.null rest) $ fail $ "Non-empty: " ++ T.unpack rest
      parser' = parser <* checkEmpty
   in case runParser parser' "scheme" input of
        Right val -> return val
        Left err -> throwError $ ParserError (T.pack $ errorBundlePretty err)

readExpr :: Text -> SchemeResult SchemeVal
readExpr = parseOrThrow pExpr

readManyExpr :: Text -> SchemeResult [SchemeVal]
readManyExpr = parseOrThrow $ some pExpr


text:
module CORDICSpec where

import qualified Clash.Prelude as Clash
import Clash.Prelude (Signal, Vec(..), BitVector, Index, Signed, Unsigned, SFixed, Bit, SNat(..),
                      simulate, simulate_lazy, listToVecTH, KnownNat, pack, unpack, (++#), mealy, mux, bundle, unbundle,
                      HiddenClockResetEnable, sample, System)
import Test.Hspec
import Test.QuickCheck hiding (sample)

import qualified Data.Complex as C
import Clash.DSP.CORDIC

approxEqual :: Double -> Double -> Bool
approxEqual x y = abs (x - y) < 0.0001

--CORDIC testing
spec = describe "CORDIC" $ do
    specify "vector mode"    $ property prop_CORDICVectorMode
    specify "rotation mode " $ property prop_CORDICRotationMode
    specify "pipelined "     $ property prop_CORDICPipelined

consts :: Vec 100 (SFixed 32 32)
consts = $(listToVecTH (Prelude.take 100 arctans))

--Test CORDIC vector mode by calculating the magnitude and phase of a complex number
prop_CORDICVectorMode =
    forAll (choose (0, 500000000)) $ \(x :: Double) -> --Restrict the input ranges so that the error does not grow beyond what approxEqual considers acceptable
        forAll (choose (0, 500000000)) $ \(y :: Double) ->
            let res     = doIt $ CordicState (fromRational (toRational x) :+ fromRational (toRational y)) 0
                cplxNum = x C.:+ y
            in     approxEqual (kValue 100 * realToFrac (realPart (cplx res))) (C.magnitude cplxNum)
                && approxEqual (realToFrac (arg res)) (C.phase cplxNum)
    where

    doIt :: CordicState (SFixed 32 32) (SFixed 32 32) -> CordicState (SFixed 32 32) (SFixed 32 32)
    doIt = cordicSteps (\(CordicState (_ :+ y) _) -> y < 0) (0 :: Index 100) consts

--Test CORDIC rotation mode by calculating the real and imaginary part of a complex number given in polar form
prop_CORDICRotationMode =
    forAll (choose (0, 5000)) $ \(x :: Double) ->  --Restrict the input ranges so that the error does not grow beyond what approxEqual considers acceptable
        forAll (choose (0, 0.75)) $ \(y :: Double) ->
            let res     = doIt $ CordicState (fromRational (toRational x) :+ 0) (fromRational (toRational y))
            in     approxEqual (kValue 100 * realToFrac (realPart (cplx res))) (x * Prelude.cos y)
                && approxEqual (kValue 100 * realToFrac (imagPart (cplx res))) (x * Prelude.sin y)
    where

    doIt :: CordicState (SFixed 32 32) (SFixed 32 32) -> CordicState (SFixed 32 32) (SFixed 32 32)
    doIt = cordicSteps (\(CordicState _ a) -> a > 0) (0 :: Index 100) consts

--Test CORDIC vector mode by calculating the magnitude and phase of a complex number
prop_CORDICPipelined =
    forAll (choose (0, 500000000)) $ \(x :: Double) -> --Restrict the input ranges so that the error does not grow beyond what approxEqual considers acceptable
        forAll (choose (0, 500000000)) $ \(y :: Double) ->
            let res     = doIt $ CordicState (fromRational (toRational x) :+ fromRational (toRational y)) 0
                cplxNum = x C.:+ y
            in     approxEqual (kValue 100 * realToFrac (realPart (cplx res))) (C.magnitude cplxNum)
                && approxEqual (realToFrac (arg res)) (C.phase cplxNum)
    where

    doIt :: CordicState (SFixed 32 32) (SFixed 32 32) -> CordicState (SFixed 32 32) (SFixed 32 32)
    doIt input
        = last
        $ take 12
        $ sample @System
        $ cordicPipeline (\(CordicState (_ :+ y) _) -> y < 0) (0 :: Index 100) (Clash.unconcat (SNat @ 10) consts) (pure True)
        $ pure input

max_stars_repo_name: adamwalker/clash-utils

max_issues_repo_name: adamwalker/clash-utils

max_forks_repo_name: adamwalker/clash-utils


text:
module MatrixPlus where

import Numeric.LinearAlgebra as N
import Test.QuickCheck

dummyMatrix :: (Field a) => Matrix a
dummyMatrix = (0 N.><0) []

mapMatrix :: (Field a)
  => (a -> a)
  -> Matrix a
  -> Matrix a
mapMatrix f x = (r N.><c) y
  where x' = toList (flatten x)
        y = map f x'
        r = rows x
        c = cols x

zipMatricesWith :: (Field a)
  => (a -> a -> a)
  -> Matrix a
  -> Matrix a
  -> Matrix a
zipMatricesWith f x y = (r N.><c) z
  where x' = toList (flatten x)
        y' = toList (flatten y)
        z = zipWith f x' y'
        r = rows x
        c = cols x

hadamardProduct :: (Field a)
  => Matrix a
  -> Matrix a
  -> Matrix a
hadamardProduct = zipMatricesWith (*)

average :: (Field a)
  => Matrix a
  -> a
average m = sum ms / fromIntegral (length ms)
  where ms = toList (flatten m)

magnitude :: (Field a, Floating a)
  => Matrix a
  -> a
magnitude x =
  if cols x == 1
  then sqrt (sum xsxs)
  else error "not a column vector"
    where xs = toList (flatten x)
          xsxs = zipWith (*) xs xs

pseudoMagnitude :: (Field a, Floating a)
  => Matrix a
  -> a
pseudoMagnitude m = sqrt (sum msms)
    where ms = toList (flatten m)
          msms = zipWith (*) ms ms

-- | Inputs, outputs and targets are represented as column vectors instead of lists
type ColumnVector a = Matrix a

-- | Convert a list to a column vector
listToColumnVector :: (Ord a, Field a)
    -- | the list to convert
    => [a]
    -- | the resulting column vector
    -> ColumnVector a
listToColumnVector x = (len N.><1 ) x
    where len = length x

-- | Convert a column vector to a list
columnVectorToList :: (Ord a, Field a)
    -- | The column vector to convert
    => ColumnVector a
    -- | The resulting list
    -> [a]
columnVectorToList = toList . flatten

-- | Testable property:
-- | If we take a non-empty list, we should be able to convert it to a column vector
-- | and back, and get the result we started with.
-- | To test this property, in GHCi type: quickCheck prop_tocolumnVectorToListRoundtrippable
prop_tocolumnVectorToListRoundtrippable :: [Double] -> Property
prop_tocolumnVectorToListRoundtrippable x =
    (length x > 0) ==>
        -- collect x $ -- uncomment to view test data
        (columnVectorToList (listToColumnVector x) == x)

max_stars_repo_name: mhwombat/backprop-example

max_issues_repo_name: mhwombat/backprop-example

max_forks_repo_name: mhwombat/backprop-example


text:
module Lambda.Parse where

-- Local Imports
import Lambda.Types
-- Global Imports
import Control.Monad
import Control.Monad.Except
import System.Environment
import Text.ParserCombinators.Parsec hiding (spaces)
import Numeric
import Data.Ratio
import Data.Complex

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

spaces :: Parser ()
spaces = skipMany1 space

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

-- Parsing
parseAtom :: Parser LispVal
parseAtom = do
  first <- letter <|> symbol
  rest  <- many (letter <|> digit <|> symbol)
  let atom = first:rest
  return $
    case atom of
      "#t" -> Bool True
      "#f" -> Bool False
      _    -> Atom atom

parseNumber :: Parser LispVal
parseNumber = many1 digit >>= return . Number . read

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

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

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

hex2dig x = fst $ readHex x !! 0

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

oct2dig x = fst $ readOct x !! 0

parseBin :: Parser LispVal
parseBin = do try $ string "#b"
              x <- many1 (oneOf "10")
              return $ Number (bin2dig 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

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))

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

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

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

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

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 -> (value !! 0)

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

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

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

parseBracketed :: Parser LispVal -> Parser LispVal
parseBracketed p = try $ do
  char '('
  many space
  x <- p
  char ')'
  return x

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

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

parseVector :: Parser LispVal
parseVector = do
  string "#("
  many space
  x <- liftM Vector $ sepEndBy parseExpr spaces
  char ')'
  return x

parseExpr :: Parser LispVal
parseExpr = parseAtom
  <|> parseString
  <|> try parseFloat
  <|> try parseRatio
  <|> try parseComplex
  <|> try parseNumber
  <|> try parseBool
  <|> try parseCharacter
  <|> try parseQuoted
  <|> try parseUnQuoted
  <|> try parseQuasiQuoted
  <|> parseBracketed parseList
  <|> parseBracketed parseDottedList
  <|> try parseVector

max_stars_repo_name: tdfirth/lambda

max_issues_repo_name: tdfirth/lambda

max_forks_repo_name: tdfirth/lambda


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

This module provides the Convolution layer, which is critical in many computer vision tasks.

-}
module Grenade.Layers.Convolution (
    Convolution (..)
  , Convolution' (..)
  , randomConvolution
  ) where

import           Control.Monad.Random hiding ( fromList )
import           Data.Maybe
import           Data.Proxy
import           Data.Serialize
import           Data.Singletons.TypeLits

#if MIN_VERSION_base(4,11,0)
import           GHC.TypeLits hiding (natVal)
#else
import           GHC.TypeLits
#endif
#if MIN_VERSION_base(4,9,0)
import           Data.Kind (Type)
#endif

import           Numeric.LinearAlgebra hiding ( uniformSample, konst )
import qualified Numeric.LinearAlgebra as LA
import           Numeric.LinearAlgebra.Static hiding ((|||), build, toRows)

import           Grenade.Core
import           Grenade.Layers.Internal.Convolution
import           Grenade.Layers.Internal.Update

-- | A convolution layer for a neural network.
--   This uses the im2col convolution trick popularised by Caffe, which essentially turns the
--   many, many, many, many loop convolution into a single matrix multiplication.
--
--   The convolution layer takes all of the kernels for the convolution, which are flattened
--   and then put into columns in the matrix.
--
--   The kernel size dictates which input and output sizes will "fit". Fitting the equation:
--   `out = (in - kernel) / stride + 1` for both dimensions.
--
--   One probably shouldn't build their own layer, but rather use the randomConvolution function.
data Convolution :: Nat -- Number of channels, for the first layer this could be RGB for instance.
                 -> Nat -- Number of filters, this is the number of channels output by the layer.
                 -> Nat -- The number of rows in the kernel filter
                 -> Nat -- The number of column in the kernel filter
                 -> Nat -- The row stride of the convolution filter
                 -> Nat -- The columns stride of the convolution filter
                 -> Type where
  Convolution :: ( KnownNat channels
                 , KnownNat filters
                 , KnownNat kernelRows
                 , KnownNat kernelColumns
                 , KnownNat strideRows
                 , KnownNat strideColumns
                 , KnownNat kernelFlattened
                 , kernelFlattened ~ (kernelRows * kernelColumns * channels))
              => !(L kernelFlattened filters) -- The kernel filter weights
              -> !(L kernelFlattened filters) -- The last kernel update (or momentum)
              -> Convolution channels filters kernelRows kernelColumns strideRows strideColumns

data Convolution' :: Nat -- Number of channels, for the first layer this could be RGB for instance.
                  -> Nat -- Number of filters, this is the number of channels output by the layer.
                  -> Nat -- The number of rows in the kernel filter
                  -> Nat -- The number of column in the kernel filter
                  -> Nat -- The row stride of the convolution filter
                  -> Nat -- The columns stride of the convolution filter
                  -> Type where
  Convolution' :: ( KnownNat channels
                  , KnownNat filters
                  , KnownNat kernelRows
                  , KnownNat kernelColumns
                  , KnownNat strideRows
                  , KnownNat strideColumns
                  , KnownNat kernelFlattened
                  , kernelFlattened ~ (kernelRows * kernelColumns * channels))
               => !(L kernelFlattened filters) -- The kernel filter gradient
               -> Convolution' channels filters kernelRows kernelColumns strideRows strideColumns

instance Show (Convolution c f k k' s s') where
  show (Convolution a _) = renderConv a
    where
      renderConv mm =
        let m  = extract mm
            ky = fromIntegral $ natVal (Proxy :: Proxy k)
            rs = LA.toColumns m
            ms = map (take ky) $ toLists . reshape ky <$> rs

            render n'  | n' <= 0.2  = ' '
                       | n' <= 0.4  = '.'
                       | n' <= 0.6  = '-'
                       | n' <= 0.8  = '='
                       | otherwise =  '#'

            px = (fmap . fmap . fmap) render ms
        in unlines $ foldl1 (zipWith (\a' b' -> a' ++ "   |   " ++ b')) $ px

randomConvolution :: ( MonadRandom m
                     , KnownNat channels
                     , KnownNat filters
                     , KnownNat kernelRows
                     , KnownNat kernelColumns
                     , KnownNat strideRows
                     , KnownNat strideColumns
                     , KnownNat kernelFlattened
                     , kernelFlattened ~ (kernelRows * kernelColumns * channels))
                  => m (Convolution channels filters kernelRows kernelColumns strideRows strideColumns)
randomConvolution = do
    s     <- getRandom
    let wN = uniformSample s (-1) 1
        mm = konst 0
    return $ Convolution wN mm

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * channels)
         ) => UpdateLayer (Convolution channels filters kernelRows kernelColumns strideRows strideColumns) where
  type Gradient (Convolution channels filters kernelRows kernelColumns strideRows strideColumns) = (Convolution' channels filters kernelRows kernelColumns strideRows strideColumns)
  runUpdate LearningParameters {..} (Convolution oldKernel oldMomentum) (Convolution' kernelGradient) =
    let (newKernel, newMomentum) = descendMatrix learningRate learningMomentum learningRegulariser oldKernel kernelGradient oldMomentum
    in Convolution newKernel newMomentum

  createRandom = randomConvolution

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * channels)
         ) => Serialize (Convolution channels filters kernelRows kernelColumns strideRows strideColumns) where
  put (Convolution w _) = putListOf put . toList . flatten . extract $ w
  get = do
      let f  = fromIntegral $ natVal (Proxy :: Proxy filters)
      wN    <- maybe (fail "Vector of incorrect size") return . create . reshape f . LA.fromList =<< getListOf get
      let mm = konst 0
      return $ Convolution wN mm

-- | A three dimensional image (or 2d with many channels) can have
--   an appropriately sized convolution filter run across it.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , KnownNat channels
         , ((outputRows - 1) * strideRows) ~ (inputRows - kernelRows)
         , ((outputCols - 1) * strideCols) ~ (inputCols - kernelCols)
         , KnownNat (kernelRows * kernelCols * channels)
         , KnownNat (outputRows * filters)
         ) => Layer (Convolution channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) where

  type Tape (Convolution channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (Convolution kernel _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)

        c  = vid2col kx ky sx sy ix iy ex
        mt = c LA.<> ek
        r  = col2vid 1 1 1 1 ox oy mt
        rs = fromJust . create $ r
    in  (S3D input, S3D rs)
  runBackwards (Convolution kernel _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)

        c  = vid2col kx ky sx sy ix iy ex

        eo = extract dEdy
        ek = extract kernel

        vs = vid2col 1 1 1 1 ox oy eo

        kN = fromJust . create $ tr c LA.<> vs

        dW = vs LA.<> tr ek

        xW = col2vid kx ky sx sy ix iy dW
    in  (Convolution' kN, S3D . fromJust . create $ xW)

-- | A two dimentional image may have a convolution filter applied to it
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , ((outputRows - 1) * strideRows) ~ (inputRows - kernelRows)
         , ((outputCols - 1) * strideCols) ~ (inputCols - kernelCols)
         , KnownNat (kernelRows * kernelCols * 1)
         , KnownNat (outputRows * filters)
         ) => Layer (Convolution 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 outputRows outputCols filters) where
  type Tape (Convolution 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    runForwards c (S3D input :: S ('D3 inputRows inputCols 1))

  runBackwards c tape grads =
    case runBackwards c tape grads of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) ->  (c', S2D back)

-- | A two dimensional image may have a convolution filter applied to it producing
--   a two dimensional image if both channels and filters is 1.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , ((outputRows - 1) * strideRows) ~ (inputRows - kernelRows)
         , ((outputCols - 1) * strideCols) ~ (inputCols - kernelCols)
         , KnownNat (kernelRows * kernelCols * 1)
         , KnownNat (outputRows * 1)
         ) => Layer (Convolution 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 outputRows outputCols) where
  type Tape (Convolution 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 outputRows outputCols) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    case runForwards c (S3D input :: S ('D3 inputRows inputCols 1)) of
      (tps, S3D back :: S ('D3 outputRows outputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    case runBackwards c tape (S3D grads :: S ('D3 outputRows outputCols 1)) of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) -> (c', S2D back)

-- | A three dimensional image can produce a 2D image from a convolution with 1 filter
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , KnownNat channels
         , ((outputRows - 1) * strideRows) ~ (inputRows - kernelRows)
         , ((outputCols - 1) * strideCols) ~ (inputCols - kernelCols)
         , KnownNat (kernelRows * kernelCols * channels)
         , KnownNat (outputRows * 1)
         ) => Layer (Convolution channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 outputRows outputCols) where
  type Tape (Convolution channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 outputRows outputCols) = S ('D3 inputRows inputCols channels)
  runForwards c input =
    case runForwards c input of
      (tps, S3D back :: S ('D3 outputRows outputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    runBackwards c tape (S3D grads :: S ('D3 outputRows outputCols 1))

max_issues_repo_name: Alien-Inc/grenade

max_forks_repo_name: Alien-Inc/grenade


text:
module Language.Parser.Number where

import qualified Text.ParserCombinators.Parsec as P
import Language.Definition
import qualified Language.Number as N
import qualified Numeric
import qualified Data.Ratio as Ratio
import qualified Data.Complex as C

parse :: P.Parser LispVal
parse = P.try complex P.<|> P.try rational P.<|> P.try float P.<|> decimal1 P.<|> (P.char '#' >> (decimal2 P.<|> hex P.<|> oct P.<|> bin))
-- simple implementation
-- parse = M.liftM (Number . read) $ P.many1 P.digit
--
-- 'do' notation
-- parse :: P.Parser LispVal
-- parse = do
--                 x <- P.many1 P.digit
--                 return $ Number . read $ x
-- >>=
-- parse :: P.Parser LispVal
-- parse =  P.many1 P.digit >>= return . Number . read

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

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

hex :: P.Parser LispVal
-- hex = do try $ string "#x"
hex = do P.string "x"
         x <- P.many1 P.hexDigit
         return $ Number (N.hex2dec x)

oct :: P.Parser LispVal
-- hex = do try $ string "#x"
oct = do P.string "o"
         x <- P.many1 P.octDigit
         return $ Number (N.oct2dec x)

bin :: P.Parser LispVal
-- hex = do try $ string "#x"
bin = do P.string "b"
         x <- P.many1 (P.oneOf "10")
         return $ Number (N.bin2dec x)

-- @todo: add support for precision
float :: P.Parser LispVal
float = do x <- P.many1 P.digit
           P.char '.'
           y <- P.many1 P.digit
           return $ Float (fst . head $ Numeric.readFloat ( x ++ "." ++ y ) )

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

complex :: P.Parser LispVal
complex = do x <- (P.try float P.<|> decimal1)
             P.char '+'
             y <- (P.try float P.<|> decimal1)
             P.char 'i'
             return $ Complex (N.toDouble x C.:+ N.toDouble y)

max_stars_repo_name: ahmadnazir/scheme

max_issues_repo_name: ahmadnazir/scheme

max_forks_repo_name: ahmadnazir/scheme


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

{-# LANGUAGE AllowAmbiguousTypes   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLabels      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-|
Module      : Grenade.Layers.Mul
Description : Scalar multiplication with some constant, supporting scaling
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Mul
  (
  -- * Layer Definition
    Mul (..)

  -- * Helper function
  , initMul
  )
where

import           Data.Kind                    (Type)
import           Data.Proxy
import           Data.Serialize
import           GHC.TypeLits

import qualified Numeric.LinearAlgebra        as LA
import           Numeric.LinearAlgebra.Static (R)
import qualified Numeric.LinearAlgebra.Static as H

import           Lens.Micro                   ((^.))

import           Grenade.Core
import           Grenade.Onnx

-- | A layer allowing for element wise multiplication with broadcasting
data Mul :: Nat -- The number of channels of the bias
         -> Nat -- The number of rows of the bias
         -> Nat -- The number of columns of the bias
         -> Type where
  Mul  :: ( KnownNat channels
          , KnownNat rows
          , KnownNat columns )
          => R (channels * rows * columns)
          -> Mul channels rows columns

instance Show (Mul c h w) where
  show (Mul mat) = "Mul " ++ show mat

instance UpdateLayer (Mul c h w) where
  type Gradient (Mul c h w) = ()
  runUpdate _ x _  = x
  reduceGradient _ = ()

instance (KnownNat c, KnownNat h, KnownNat w ) => RandomLayer (Mul c h w) where
  createRandomWith _ _ = pure initMul

-- ^ Initialize a Mul layer with scale 1
initMul :: forall c h w. ( KnownNat c, KnownNat h, KnownNat w )
        => Mul c h w
initMul =
  let c'    = fromIntegral $ natVal (Proxy :: Proxy c)
      h'    = fromIntegral $ natVal (Proxy :: Proxy h)
      w'    = fromIntegral $ natVal (Proxy :: Proxy w)
      ones  = replicate (c' * h' * w') 1
      bias  = H.fromList ones :: R (c * h * w)
  in Mul bias

instance ( KnownNat c, KnownNat h, KnownNat w ) => Serialize (Mul c h w) where
  put (Mul bias) = putListOf put . LA.toList . H.extract $ bias
  get            = do
    bias <- maybe (fail "Vector of incorrect size") return . H.create . LA.fromList =<< getListOf get
    return $ Mul bias

-- | Currently, only multiplication by a single scalar is supported.
instance ( KnownNat i, KnownNat j, KnownNat k ) => Layer (Mul 1 1 1) ('D3 i j k) ('D3 i j k) where
  type Tape (Mul 1 1 1) ('D3 i j k) ('D3 i j k) = ()

  runForwards (Mul b) (S3D m)
    = let s  = H.extract b LA.! 0
      in  ((), S3D $  H.dmmap (s *) m)

  runBackwards = error "runBackwards is not implemented for the Mul layer"

instance OnnxOperator (Mul c h w) where
  onnxOpTypeNames _ = ["Mul"]

instance (KnownNat c, KnownNat h, KnownNat w) => OnnxLoadable (Mul c h w) where
  loadOnnxNode inits node = case node ^. #input of
    [_, scale] -> do
      loadedScale <- readInitializerVector inits scale

      return $ Mul loadedScale
    _               -> onnxIncorrectNumberOfInputs

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE TemplateHaskell #-}

module Shapes where

import Control.Lens hiding (element)
import Graphics.Gloss.Data.Picture
import Data.Complex

data ShapeType = Kite | Dart deriving (Show)

data Shape = Shape {_shape :: ShapeType
                   ,_translate :: Complex Double
                   ,_rotate :: Complex Double
} deriving (Show)

newtype Vertices = Vertices {_nodes :: [Complex Double]} deriving (Show)

makeLenses ''Shape
makeLenses ''Vertices

standardKite :: Shape
standardKite = Shape {_shape = Kite, _translate = 0, _rotate = 1}

standardDart :: Shape
standardDart = Shape {_shape = Dart, _translate = 0, _rotate = 1}

z :: Complex Double
z = cis (pi/5)

zPower :: Int -> Complex Double
zPower n = cis((pi/5) * fromIntegral (n `mod` 10))

longEdge :: Vertices
longEdge = Vertices {_nodes = [0,(0.5 :+ (-0.1))*z,z]}

shortEdge :: Vertices
shortEdge = Vertices {_nodes = [0,(0.5 :+ 0.1)*(1-z),1-z]}

rotranslate :: Complex Double -> Complex Double -> Vertices -> Vertices
rotranslate z w = over nodes (map (\u -> u*z+w))

reverseVertices :: Vertices -> Vertices
reverseVertices = over nodes reverse

concatVertices :: [Vertices] -> Vertices
concatVertices p = Vertices {_nodes = concatMap (tail._nodes) p}

shapeToVertices :: Shape -> Vertices
shapeToVertices (Shape Kite z x) = rotranslate x z kiteVertices
shapeToVertices (Shape Dart z x) = rotranslate x z dartVertices

complexToPoint :: Complex Double -> Point
complexToPoint z = (realToFrac x, realToFrac y)
  where
    x = realPart z
    y = imagPart z

verticesToPicture :: Vertices -> Picture
verticesToPicture = lineLoop . map complexToPoint . view nodes

shapeToPicture :: Shape -> Picture
shapeToPicture = verticesToPicture . shapeToVertices

refine :: Shape -> [Shape]
refine s@(Shape t w z) =
  let phi = (sqrt 5 - 1)/2 :+ 0
      tr t' u v = Shape t' (z*phi*v+w) (z*u*phi)
  in case t of
  Kite -> [tr Dart 1 0
          ,tr Kite (zPower 3) (1 - zPower 3)
          ,tr Kite (zPower 7) (1 + zPower 2)
          ]
  Dart -> [tr Kite (zPower 9) 0
          ,tr Dart (zPower 6) (1+phi)
          ]

tilings :: [[Shape]]
tilings = [standardKite] : map (concatMap refine) tilings

tilingPics :: [Picture]
tilingPics = map (pictures . map shapeToPicture) tilings

kiteVertices :: Vertices
kiteVertices = concatVertices
  [longEdge
  ,rotranslate 1 z shortEdge
  ,reverseVertices $ rotranslate (zPower 4) (zPower 9) shortEdge
  ,reverseVertices $ rotranslate (zPower 8) 0 longEdge
  ]

dartVertices :: Vertices
dartVertices = concatVertices
  [reverseVertices $ rotranslate (zPower 4) 1 longEdge
  ,reverseVertices $ rotranslate (zPower 3) x shortEdge
  ,rotranslate (zPower 9) x shortEdge
  ,rotranslate (zPower 2) (zPower 8) longEdge
  ]
  where
    x = (1 + (zPower 4) - (zPower 3))

printVertices :: Vertices -> IO ()
printVertices = mapM_ print . _nodes

max_stars_repo_name: bnarnold/Penrose

max_issues_repo_name: bnarnold/Penrose

max_forks_repo_name: bnarnold/Penrose


text:
{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, FlexibleInstances #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}
module Data.NeuralNetwork.Backend.HMatrix (
  module Data.NeuralNetwork.Backend.HMatrix.Layers,
  ByHmatrix(..),
  ErrCode(..)
) where

import Data.NeuralNetwork
import Data.NeuralNetwork.Common
import Data.NeuralNetwork.Backend.HMatrix.Layers
import Numeric.LinearAlgebra (Vector, Matrix)
import Control.Monad.Except
import Data.Functor.Identity
import Data.Constraint (Dict(..))

type Err     = ExceptT ErrCode IO

-- the backend type
data ByHmatrix = ByHmatrix

instance (HeadSize z, TranslateBody Err s,
          Component (SpecToCom s),
          RunInEnv (Run (SpecToCom s)) Err)
       => Backend ByHmatrix (z :++ s) where
  type Env ByHmatrix = Err
  type ConvertFromSpec ByHmatrix (z :++ s) = SpecToCom s
  compile _ (a :++ l)= trans (hsize a) l
  witness _ _ = Dict

instance RunInEnv Identity Err where
  run = return . runIdentity

instance TranslateBody Err SpecFullConnect where
  type SpecToCom SpecFullConnect = RunLayer (S F (T (SinglC :. Vector)))
  trans (D1 s) (FullConnect n) = do u <- lift $ newFLayer s n
                                    return $ Stack u (Activation (relu, relu'))
  trans _ _ = throwError ErrMismatch

instance TranslateBody Err SpecConvolution where
  type SpecToCom SpecConvolution = RunLayer (S C (T (MultiC :. Matrix)))
  trans (D2 k s t) (Convolution n f p) = do u <- lift $ newCLayer k n f p
                                            return $ Stack u (Activation (relu, relu'))
  trans _ _ = throwError ErrMismatch

instance TranslateBody Err SpecReshape2DAs1D where
  type SpecToCom SpecReshape2DAs1D = RunLayer A
  trans (D2 _ _ _) _ = return As1D
  trans (D1 _)     _ = throwError ErrMismatch

instance TranslateBody Err SpecMaxPooling where
  type SpecToCom SpecMaxPooling = RunLayer M
  trans (D2 _ _ _) (MaxPooling n) = return (MaxP n)
  trans (D1 _)     _              = throwError ErrMismatch

max_stars_repo_name: pierric/neural-network

max_issues_repo_name: pierric/neural-network

max_forks_repo_name: pierric/neural-network


text:
-- Jacobi's Method (Nonlinear Equation System Solution)
-- May 13th, 2016
import Data.Matrix
import Data.Complex

jacobiroots :: (RealFloat a) => (Matrix (Matrix a -> a))     -- F(x) vector
                             -> (Matrix a)                   -- initial state
                             -> a                            -- maximum epsilon
                             -> ((a -> a) -> a -> a -> a)    -- given solver (f, x0, emax)
                             -> (Matrix a)                   -- result
jacobiroots f x0 e sv = if abs (maximum (x1 - x0)) > e then jacobiroots f x1 e sv else x1 where
   n  = nrows x0
   x1 = foldl (\acc x -> setElem (sv (\y -> (f ! (x, 1)) $ matrix n 1 (\(z, _) -> if z /= x then x0 ! (z, 1) else y)) (x0 ! (x, 1)) e) (x, 1) acc) (matrix n 1 (\_ -> 0)) [1..n]

-- F(x) is formatted as follows
-- fns = fromList 2 1 [(\[x, y] -> x + y) . toList, (\[x, y] -> x - 2*y) . toList]

secroot :: (RealFloat a) => (a -> a)                -- f(x)
                         -> a                       -- initial state 1
                         -> a                       -- initial state 2
                         -> a                       -- maximum epsilon
                         -> a                       -- result
secroot f x1 x2 e = if abs ((-) x3 x2) > e then secroot f x2 x3 e else x3 where
    f1   = f x1
    f2   = f x2
    fpi  = if f2 /= f1 then (/) (x2 - x1) (f2 - f1) else 0
    x3   = if f2 /= f1 then x2 - f2 * fpi           else x2 + e -- reseed

-- F(x) is formatted as follows
-- fns = fromList 2 1 [(\[x, y] -> x + y) . toList, (\[x, y] -> x - 2*y) . toList]
-- complex counterparts
jacobiroots_c :: (RealFloat a) => (Matrix (Matrix (Complex a) -> (Complex a)))     -- F(x) vector
                             -> (Matrix (Complex a))                   -- initial state
                             -> (Complex a)                            -- maximum epsilon
                             -> (((Complex a) -> (Complex a)) -> (Complex a) -> (Complex a) -> (Complex a))    -- given solver (f, x0, emax)
                             -> (Matrix (Complex a))                   -- result
jacobiroots_c f x0 e sv = if (maximum (fmap magnitude (x1 - x0))) > magnitude e then jacobiroots_c f x1 e sv else x1 where
   n  = nrows x0
   x1 = foldl (\acc x -> setElem (sv (\y -> (f ! (x, 1)) $ matrix n 1 (\(z, _) -> if z /= x then x0 ! (z, 1) else y)) (x0 ! (x, 1)) e) (x, 1) acc) (matrix n 1 (\_ -> 0)) [1..n]

secroot_c :: (RealFloat a) => (Complex a -> Complex a)                -- f(x)
                         -> Complex a                       -- initial state 1
                         -> Complex a                       -- initial state 2
                         -> Complex a                       -- maximum epsilon (real, Complex for type only)
                         -> Complex a                       -- result
secroot_c f x1 x2 e = if magnitude ((-) x3 x2) > magnitude e then secroot_c f x2 x3 e else x3 where
    f1   = f x1
    f2   = f x2
    fpi  = if f2 /= f1 then (/) (x2 - x1) (f2 - f1) else 0
    x3   = if f2 /= f1 then x2 - f2 * fpi           else x2 + e -- reseed


text:
-- | <<docs/apollonian_gasket/spaced_gasket.svg>>
module Geometry.Processes.ApollonianGasket (
    -- * The classical shape
      createGasket
    , fifthCircle
    , Circle(..)

    -- * For playing around
    , ApoCircle(..)
    , newCircle
    , toCircle
    , toApoCircle
) where

import Data.Complex
import Data.Tree

import Geometry.Core as G

-- | A circle, but in a format more suitable for Descartes’ Theorem, which the
-- gasket algorithm is based on.
--
-- It interprets the center @'Vec2' x y@ as a complex number \(x + \mathrm i\,y\),
-- and stores the curvature \(k=\frac1r\) instead of the radius \(r\).
--
-- 'Circle' and 'ApoCircle' can be freely converted between using 'toCircle' and
-- 'toApoCircle'.
data ApoCircle = ApoCircle (Complex Double) Double
    deriving (Eq, Show)

instance HasBoundingBox ApoCircle where
    boundingBox = boundingBox . toCircle

-- | Simple conversion function. Inverse of 'toApoCircle'.
toCircle :: ApoCircle -> Circle
toCircle (ApoCircle (x :+ y) k) = Circle (Vec2 x y) (1/k)

-- | Simple conversion function. Inverse of 'toCircle'.
toApoCircle :: Circle -> ApoCircle
toApoCircle (Circle (Vec2 x y) r) = ApoCircle (x :+ y) (1/r)

-- | Curvature of the new Apollonian circle.
newCurvature
    :: (Double -> Double -> Double)
    -> Double -- ^ Curvature of first circle
    -> Double -- ^ Second circle
    -> Double -- ^ Third circle
    -> Double -- ^ New circle
newCurvature (+-) k1 k2 k3 =
    -- Taken straight from Wikipedia: https://en.wikipedia.org/wiki/Descartes%27_theorem
    (k1 + k2 + k3) +- (2 * sqrt (k1*k2 + k2*k3 + k3*k1))

-- | Center of the new Apollonian circle.
newCenter
    :: (Complex Double -> Complex Double -> Complex Double) -- ^ Plus or minus to yield both circles
    -> ApoCircle
    -> ApoCircle
    -> ApoCircle
    -> Double         -- ^ Curvature of the new circle
    -> Complex Double -- ^ Center of the new circle
newCenter (+-) (ApoCircle c1 k1) (ApoCircle c2 k2) (ApoCircle c3 k3) k4 =
    -- Taken straight from Wikipedia: https://en.wikipedia.org/wiki/Descartes%27_theorem
    1/k4' * ( (c1*k1' + c2*k2' + c3*k3')
              +- (2*sqrt (k1'*k2'*c1*c2 + k2'*k3'*c2*c3 + k1'*k3'*c1*c3)))
  where
    [k1', k2', k3', k4'] = map (:+ 0) [k1, k2, k3, k4]

-- | Create a new Apollonian circle, based on three existing ones. This is the raw
-- application of Descartes’ theorem. For something with less pitfalls, use
-- 'fifthCircle'.
--
-- The choice of sign is a mystery to me. For the Apollonian Gasket, the rule is:
--
--   1. The large circle is obtained using \((-,-)\) on the three initial cycles
--   2. The bottom-left circle is obtained using \((+,-)\) on the initial top, initial bottom, and large circle
--   3. All other recursive steps use \((+,+)\)
--   4. If you can explain the \((+,-)\), please contact me, I’m dying to know
--
-- If you’re wondering, this is what you get when you don’t do the \((+,-)\) trick:
--
-- <<docs/apollonian_gasket/missing_the_minus.svg>>
newCircle
    :: (Double -> Double -> Double)
        -- ^ Plus (larger curvature, smaller circle) or minus (smaller curvature, larger circle).
    -> (Complex Double -> Complex Double -> Complex Double) -- ^ Plus or minus; see this function’s doc.
    -> ApoCircle
    -> ApoCircle
    -> ApoCircle
    -> ApoCircle
newCircle plusMinusCurvature plusMinusCenter circ1@(ApoCircle _ k1) circ2@(ApoCircle _ k2) circ3@(ApoCircle _ k3) =
    let k4 = newCurvature plusMinusCurvature k1 k2 k3
        c4 = newCenter plusMinusCenter circ1 circ2 circ3 k4
    in ApoCircle c4 k4

-- | Given four mutually tangent circles, find »the other one« than the first argument.
-- This resolves the »\((+,-)\)« ambiguity of 'newCircle'.
fifthCircle :: ApoCircle -> ApoCircle -> ApoCircle -> ApoCircle -> ApoCircle
fifthCircle (ApoCircle c1 k1) (ApoCircle c2 k2) (ApoCircle c3 k3) (ApoCircle c4 k4) =
    let [k1', k2', k3', k4', k5'] = map (:+ 0) [k1, k2, k3, k4, k5]
        k5 = 2*(k2+k3+k4) - k1
        c5 = (2*(k2'*c2+k3'*c3+k4'*c4) - k1'*c1) / k5'
    in ApoCircle c5 k5

-- | The simple workhorse function. Given three mutually touching, equally sized
-- circles, it calculates the classical Apollonian Gasket:
--
--
-- <<docs/apollonian_gasket/classical_gasket.svg>>
--
-- Be warned though, the number of circles grows exponentially with recursion depth
-- (controlled with the radius parameter). This can easily melt your RAM, CPU or
-- harddrive if you aren’t careful.
--
-- It’s quite fun to play with the parameters though! Here is what happens when the three initial circles
-- aren’t touching:
--
-- <<docs/apollonian_gasket/forgetting_gen0.svg>>
--
-- To get even more creative freedom, I like to copy this algorithm and see what I
-- get, e.g. by commenting out one of the recursions.
createGasket
    :: Double -- ^ Radius threshold: do not recurse to circles smaller than this. Smaller values yield more detail.
    -> Circle -- ^ Initial left circle
    -> Circle -- ^ Initial right circle
    -> Circle -- ^ Initial bottom circle
    -> Tree Circle
createGasket minRadius gen0LCirc gen0RCirc gen0BCirc =
    let gen0L = toApoCircle gen0LCirc
        gen0R = toApoCircle gen0RCirc
        gen0B = toApoCircle gen0BCirc

        large = newCircle (-) (-) gen0L gen0R gen0B
        small = newCircle (+) (+) gen0L gen0R gen0B

        recurse :: ApoCircle -> ApoCircle -> ApoCircle -> ApoCircle -> Tree Circle
        recurse cA c1 c2 c3 =
            let cB@(ApoCircle _ kB) = fifthCircle cA c1 c2 c3
            in Node (toCircle cB)
                (if abs kB < 1/minRadius
                    then
                        [ recurse c3 cB c1 c2
                        , recurse c2 c3 cB c1
                        , recurse c1 c2 c3 cB
                        ]
                    else [])

        circles :: Tree Circle
        circles = Node (toCircle large)
            [ recurse small gen0L gen0R gen0B
            , recurse large gen0L gen0R gen0B
            ]
    in circles

max_stars_repo_name: quchen/haskell-logo-algorithmic

max_issues_repo_name: quchen/haskell-logo-algorithmic

max_forks_repo_name: quchen/haskell-logo-algorithmic


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE OverloadedLists     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeOperators       #-}

import Data.Maybe
import Data.Number.Symbolic
import qualified Data.Number.Symbolic as Sym
import Data.Proxy

import qualified Naperian as N
import qualified Data.Foldable as F
import           Control.Applicative ( liftA2 )
import qualified GHC.TypeLits as M
import           Data.Functor
import           Data.List.Split

import           Numeric.Sundials.ARKode.ODE
import           Numeric.LinearAlgebra

kx, ky :: Floating a => a
kx = 0.5
ky = 0.75

-- spatial mesh size
nx, ny :: Int
nx = 30
ny = 60

-- x mesh spacing
-- y mesh spacing
dx :: Floating a => a
dx = 1 / (fromIntegral nx - 1)

dy :: Floating a => a
dy = 1 / (fromIntegral ny - 1)

c1, c2 :: Floating a => a
c1 = kx/dx/dx
c2 = ky/dy/dy

cc4' :: forall b m n . (M.KnownNat m, M.KnownNat n, Num b) =>
        N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] b
cc4' = N.Prism $ N.Prism $ N.Prism $ N.Prism $ N.Scalar $
      N.viota @m <&> (\(N.Fin x) ->
      N.viota @n <&> (\(N.Fin w) ->
      N.viota @m <&> (\(N.Fin v) ->
      N.viota @n <&> (\(N.Fin u) ->
      (f m n x w v u)))))
        where
          m = fromIntegral $ M.natVal (undefined :: Proxy m)
          n = fromIntegral $ M.natVal (undefined :: Proxy n)
          f m n i j k l | i == 0               = 0
                        | j == 0               = 0
                        | i == n - 1           = 0
                        | j == m - 1           = 0
                        | k == i - 1 && l == j = 1
                        | k == i     && l == j = -2
                        | k == i + 1 && l == j = 1
                        | otherwise            = 0

cc5' :: forall a m n . (M.KnownNat m, M.KnownNat n, Floating a) =>
        N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] a
cc5' = N.binary (*) (N.Scalar c2) cc4'

cc5Sym' :: forall a m n . (M.KnownNat m, M.KnownNat n, Floating a, Eq a) =>
          N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] (Sym a)
cc5Sym' = N.binary (*) (N.Scalar $ var "c2") cc4'

yy4' :: forall b m n . (M.KnownNat m, M.KnownNat n, Num b) =>
        N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] b
yy4' = N.Prism $ N.Prism $ N.Prism $ N.Prism $ N.Scalar $
      N.viota @m <&> (\(N.Fin x) ->
      N.viota @n <&> (\(N.Fin w) ->
      N.viota @m <&> (\(N.Fin v) ->
      N.viota @n <&> (\(N.Fin u) ->
      (f m n x w v u)))))
        where
          m = fromIntegral $ M.natVal (undefined :: Proxy m)
          n = fromIntegral $ M.natVal (undefined :: Proxy n)
          f :: Int -> Int -> Int -> Int -> Int -> Int -> b
          f m n i j k l | i == 0                   = 0
                        | j == 0                   = 0
                        | i == n - 1               = 0
                        | j == m - 1               = 0
                        | k == i     && l == j - 1 = 1
                        | k == i     && l == j     = -2
                        | k == i     && l == j + 1 = 1
                        | otherwise                = 0

yy5' :: forall a m n . (M.KnownNat m, M.KnownNat n, Floating a) =>
        N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] a
yy5' = N.binary (*) (N.Scalar c1) yy4'

yy5Sym' :: forall a m n . (M.KnownNat m, M.KnownNat n, Floating a, Eq a) =>
           N.Hyper '[N.Vector n, N.Vector m, N.Vector n, N.Vector m] (Sym a)
yy5Sym' = N.binary (*) (N.Scalar $ var "c1") yy4'

ccSym5 = cc5Sym' @Double @4 @5
yy5Sym = yy5Sym' @Double @4 @5

ccSym5

yy5Sym

fmap (N.elements . N.Prism . N.Prism . N.Scalar) $ N.elements $ N.crystal $ N.crystal $ N.binary (+) cc5Sym yy5Sym

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE OverloadedLists       #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}

import Data.Maybe
import Data.Number.Symbolic
import Data.Proxy

import qualified Naperian as N
import qualified Data.Foldable as F
import           Control.Applicative ( liftA2 )
import qualified GHC.TypeLits as M

import           Numeric.Sundials.ARKode.ODE
import           Numeric.LinearAlgebra

x1, a, x2 :: Double
x1 = 0
a = 1.0
x2 = a

y1, y2 :: Double
y1 = 0.0
y2 = 1.0

bigT :: Double
bigT = 1000.0

n :: Int
n = 2

dx :: Double
dx = a / (fromIntegral n + 1)

dy :: Double
dy = a / (fromIntegral n + 1)

beta, s :: Double
beta = 1.0e-5
s = beta / (dx * dx)

kx, ky :: Double
kx = 0.5
ky = 0.75

c1, c2 :: Double
c1 = kx/dx/dx
c2 = ky/dy/dy

bigAA1 :: Matrix Double
bigAA1 = assoc (n * n, n * n) 0.0 [((i, j), f (i, j)) | i <- [0 .. n * n - 1]
                                                      , j <- [i - n, i,  i + n]
                                                      , j `elem` [0 .. n * n -1]]
  where
    f (i, j) | i     == j = (-2.0) * c1
             | i - n == j = 1.0    * c1
             | i + n == j = 1.0    * c1
             | otherwise = error $ show (i, j)

bigAA2 :: Matrix Double
bigAA2 = diagBlock (replicate n bigA)
  where
    bigA :: Matrix Double
    bigA = assoc (n, n) 0.0 [((i, j), f (i, j)) | i <- [0 .. n - 1]
                                                , j <- [i-1..i+1]
                                                , j `elem` [0..n-1]]
      where
        f (i, j) | i     == j = (-2.0) * c2
                 | i - 1 == j = 1.0    * c2
                 | i + 1 == j = 1.0    * c2

bigAA :: Matrix Double
bigAA = bigAA1 + bigAA2

bigAA1

bigAA2

n

bigZZ1 :: Matrix Double
bigZZ1 = assoc (m * m, m * m) 0.0 [((i, j), f (i, j)) | i <- [0 .. m * m - 1]
                                                      , j <- [0 .. m * m - 1]]
  where
    m = n + 2
    f (i, j) | i     == 0     = 0.0
             | j     == 0     = 0.0
             | i     == j     = (-2.0) * c1
             | i - n == j     = 1.0    * c1
             | i + n == j     = 1.0    * c1
             | i     == n + 1 = 0.0
             | j     == n + 1 = 0.0
             | otherwise      = 0.0

bigZZ1

x :: forall m n . (M.KnownNat m, M.KnownNat n) => N.Vector n (N.Vector m (Sym Int))
x = (fromJust . N.fromList) $
    map (fromJust . N.fromList) ([[var $ (\(x,y) -> "A" ++ show x ++ "," ++ show y) (x,y) | y <- [1..m]] | x <- [1..n]] :: [[Sym Int]])
    where
      m = M.natVal (undefined :: Proxy m)
      n = M.natVal (undefined :: Proxy n)

u1 :: N.Hyper '[N.Vector 3, N.Vector 2] (Sym Int)
u1 = N.Prism $ N.Prism (N.Scalar x)

u1

y :: forall n . M.KnownNat n => N.Vector n (Sym Int)
y = (fromJust . N.fromList) $
    (map (var . ("v" ++) . show) [1..n ] :: [Sym Int])
    where
    n = M.natVal (undefined :: Proxy n)

u2 :: N.Hyper '[N.Vector 3] (Sym Int)
u2 = N.Prism (N.Scalar y)

N.innerH u1 u2

max_stars_repo_path: Two Dimensional Heat Equation.hs

max_stars_repo_name: idontgetoutmuch/Diffusions

max_issues_repo_path: Two Dimensional Heat Equation.hs

max_issues_repo_name: idontgetoutmuch/Diffusions

max_forks_repo_path: Two Dimensional Heat Equation.hs

max_forks_repo_name: idontgetoutmuch/Diffusions


text:
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DeriveFoldable #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE ParallelListComp #-}
-----------------------------------------------------------------------------
-- |
-- Copyright   :  (c) Edward Kmett 2010-2021
-- License     :  BSD3
-- Maintainer  :  ekmett@gmail.com
-- Stability   :  experimental
-- Portability :  GHC only
--
-----------------------------------------------------------------------------

module Numeric.AD.Newton
  (
  -- * Newton's Method (Forward AD)
    findZero
  , findZeroNoEq
  , inverse
  , inverseNoEq
  , fixedPoint
  , fixedPointNoEq
  , extremum
  , extremumNoEq
  -- * Gradient Ascent/Descent (Reverse AD)
  , gradientDescent, constrainedDescent, CC(..), eval
  , gradientAscent
  , conjugateGradientDescent
  , conjugateGradientAscent
  , stochasticGradientDescent
  ) where

import Data.Foldable (all, sum)
import Data.Reflection (Reifies)
import Data.Traversable
import Numeric.AD.Internal.Combinators
import Numeric.AD.Internal.Forward (Forward)
import Numeric.AD.Internal.On
import Numeric.AD.Internal.Or
import Numeric.AD.Internal.Reverse (Reverse, Tape)
import Numeric.AD.Internal.Type (AD(..))
import Numeric.AD.Mode
import Numeric.AD.Mode.Reverse as Reverse (gradWith, gradWith', grad')
import Numeric.AD.Rank1.Kahn as Kahn (Kahn, grad)
import qualified Numeric.AD.Rank1.Newton as Rank1
import Prelude hiding (all, mapM, sum)

-- $setup
-- >>> import Data.Complex

-- | The 'findZero' function finds a zero of a scalar function using
-- Newton's method; its output is a stream of increasingly accurate
-- results.  (Modulo the usual caveats.) If the stream becomes constant
-- ("it converges"), no further elements are returned.
--
-- Examples:
--
-- >>> take 10 $ findZero (\x->x^2-4) 1
-- [1.0,2.5,2.05,2.000609756097561,2.0000000929222947,2.000000000000002,2.0]
--
-- >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1)
-- 0.0 :+ 1.0
findZero :: (Fractional a, Eq a) => (forall s. AD s (Forward a) -> AD s (Forward a)) -> a -> [a]
findZero f = Rank1.findZero (runAD.f.AD)
{-# INLINE findZero #-}

-- | The 'findZeroNoEq' function behaves the same as 'findZero' except that it
-- doesn't truncate the list once the results become constant. This means it
-- can be used with types without an 'Eq' instance.
findZeroNoEq :: Fractional a => (forall s. AD s (Forward a) -> AD s (Forward a)) -> a -> [a]
findZeroNoEq f = Rank1.findZeroNoEq (runAD.f.AD)
{-# INLINE findZeroNoEq #-}

-- | The 'inverse' function inverts a scalar function using
-- Newton's method; its output is a stream of increasingly accurate
-- results.  (Modulo the usual caveats.) If the stream becomes
-- constant ("it converges"), no further elements are returned.
--
-- Example:
--
-- >>> last $ take 10 $ inverse sqrt 1 (sqrt 10)
-- 10.0
inverse :: (Fractional a, Eq a) => (forall s. AD s (Forward a) -> AD s (Forward a)) -> a -> a -> [a]
inverse f = Rank1.inverse (runAD.f.AD)
{-# INLINE inverse  #-}

-- | The 'inverseNoEq' function behaves the same as 'inverse' except that it
-- doesn't truncate the list once the results become constant. This means it
-- can be used with types without an 'Eq' instance.
inverseNoEq :: Fractional a => (forall s. AD s (Forward a) -> AD s (Forward a)) -> a -> a -> [a]
inverseNoEq f = Rank1.inverseNoEq (runAD.f.AD)
{-# INLINE inverseNoEq #-}

-- | The 'fixedPoint' function find a fixedpoint of a scalar
-- function using Newton's method; its output is a stream of
-- increasingly accurate results.  (Modulo the usual caveats.)
--
-- If the stream becomes constant ("it converges"), no further
-- elements are returned.
--
-- >>> last $ take 10 $ fixedPoint cos 1
-- 0.7390851332151607
fixedPoint :: (Fractional a, Eq a) => (forall s. AD s (Forward a) -> AD s (Forward a)) -> a -> [a]
fixedPoint f = Rank1.fixedPoint (runAD.f.AD)
{-# INLINE fixedPoint #-}

-- | The 'fixedPointNoEq' function behaves the same as 'fixedPoint' except that
-- it doesn't truncate the list once the results become constant. This means it
-- can be used with types without an 'Eq' instance.
fixedPointNoEq :: Fractional a => (forall s. AD s (Forward a) -> AD s (Forward a)) -> a -> [a]
fixedPointNoEq f = Rank1.fixedPointNoEq (runAD.f.AD)
{-# INLINE fixedPointNoEq #-}

-- | The 'extremum' function finds an extremum of a scalar
-- function using Newton's method; produces a stream of increasingly
-- accurate results.  (Modulo the usual caveats.) If the stream
-- becomes constant ("it converges"), no further elements are returned.
--
-- >>> last $ take 10 $ extremum cos 1
-- 0.0
extremum :: (Fractional a, Eq a) => (forall s. AD s (On (Forward (Forward a))) -> AD s (On (Forward (Forward a)))) -> a -> [a]
extremum f = Rank1.extremum (runAD.f.AD)
{-# INLINE extremum #-}

-- | The 'extremumNoEq' function behaves the same as 'extremum' except that it
-- doesn't truncate the list once the results become constant. This means it
-- can be used with types without an 'Eq' instance.
extremumNoEq :: Fractional a => (forall s. AD s (On (Forward (Forward a))) -> AD s (On (Forward (Forward a)))) -> a -> [a]
extremumNoEq f = Rank1.extremumNoEq (runAD.f.AD)
{-# INLINE extremumNoEq #-}

-- | The 'gradientDescent' function performs a multivariate
-- optimization, based on the naive-gradient-descent in the file
-- @stalingrad\/examples\/flow-tests\/pre-saddle-1a.vlad@ from the
-- VLAD compiler Stalingrad sources.  Its output is a stream of
-- increasingly accurate results.  (Modulo the usual caveats.)
--
-- It uses reverse mode automatic differentiation to compute the gradient.
gradientDescent :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> [f a]
gradientDescent f x0 = go x0 fx0 xgx0 0.1 (0 :: Int)
  where
    (fx0, xgx0) = Reverse.gradWith' (,) f x0
    go x fx xgx !eta !i
      | eta == 0     = [] -- step size is 0
      | fx1 > fx     = go x fx xgx (eta/2) 0 -- we stepped too far
      | zeroGrad xgx = [] -- gradient is 0
      | otherwise    = x1 : if i == 10
                            then go x1 fx1 xgx1 (eta*2) 0
                            else go x1 fx1 xgx1 eta (i+1)
      where
        zeroGrad = all (\(_,g) -> g == 0)
        x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx
        (fx1, xgx1) = Reverse.gradWith' (,) f x1
{-# INLINE gradientDescent #-}

data SEnv (f :: * -> *) a = SEnv { sValue :: a, origEnv :: f a }
  deriving (Functor, Foldable, Traversable)

-- | Convex constraint, CC, is a GADT wrapper that hides the existential
-- ('s') which is so prevalent in the rest of the API.  This is an
-- engineering convenience for managing the skolems.
data CC f a where
  CC :: forall f a. (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> CC f a

-- |@constrainedDescent obj fs env@ optimizes the convex function @obj@
-- subject to the convex constraints @f <= 0@ where @f `elem` fs@. This is
-- done using a log barrier to model constraints (i.e. Boyd, Chapter 11.3).
-- The returned optimal point for the objective function must satisfy @fs@,
-- but the initial environment, @env@, needn't be feasible.
constrainedDescent :: forall f a. (Traversable f, RealFloat a, Floating a, Ord a)
                   => (forall s. Reifies s Tape => f (Reverse s a)
                                                -> Reverse s a)
                   -> [CC f a]
                   -> f a
                   -> [(a,f a)]
constrainedDescent objF [] env =
  map (\x -> (eval objF x, x)) (gradientDescent objF env)
constrainedDescent objF cs env =
    let s0       = 1 + maximum [eval c env | CC c <- cs]
        -- ^ s0 = max ( f_i(0) )
        cs'      = [CC (\(SEnv sVal rest) -> c rest - sVal) | CC c <- cs]
        -- ^ f_i' = f_i - s0  and thus f_i' <= 0
        envS     = SEnv s0 env
        -- feasible point for f_i', use gd to find feasiblity for f_i
        cc       = constrainedConvex' (CC sValue) cs' envS ((<=0) . sValue)
    in case dropWhile ((0 <) . fst) (take (2^(20::Int)) cc) of
        []                  -> []
        (_,envFeasible) : _ ->
            constrainedConvex' (CC objF) cs (origEnv envFeasible) (const True)
{-# INLINE constrainedDescent #-}

eval :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> a
eval f e = fst (grad' f e)
{-# INLINE eval #-}

-- | Like 'constrainedDescent' except the initial point must be feasible.
constrainedConvex' :: forall f a. (Traversable f, RealFloat a, Floating a, Ord a)
                   => CC f a
                   -> [CC f a]
                   -> f a
                   -> (f a -> Bool)
                   -> [(a,f a)]
constrainedConvex' objF cs env term =
  -- 1. Transform cs using a log barrier with increasing t values.
  let os   = map (mkOpt objF cs) tValues
  -- 2. Iteratively run gradientDescent on each os.
      envs =  [(undefined,env)] :
              [gD (snd $ last e) o
                          | o  <- os
                          | e  <- limEnvs
                          ]
      -- Obtain a finite number of elements from the initial len tValues - 1 lists.
      limEnvs = zipWith id nrSteps envs
  in dropWhile (not . term . snd) (concat $ drop 1 limEnvs)
 where
  tValues = map realToFrac $ take 64 $ iterate (*2) (2 :: a)
  nrSteps = [take 20 | _ <- [1..length tValues]] ++ [id]
  -- | `gD f e` is gradient descent with the evaulated result
  gD e (CC f)  = (eval f e, e) :
                 map (\x -> (eval f x, x)) (gradientDescent f e)
{-# INLINE constrainedConvex' #-}

-- @mkOpt u fs t@ converts an inequality convex problem (@u,fs@) into an
-- unconstrained convex problem using log barrier @u + -(1/t)log(-f_i)@.
-- As @t@ increases the approximation is more accurate but the gradient
-- decreases, making the gradient descent more expensive.
mkOpt :: forall f a. (Traversable f, RealFloat a, Floating a, Ord a)
      => CC f a -> [CC f a]
      -> a -> CC f a
mkOpt (CC o) xs t = CC (\e -> o e + sum (map (\(CC c) -> iHat t c e) xs))
{-# INLINE mkOpt #-}

iHat :: forall a f. (Traversable f, RealFloat a, Floating a, Ord a)
     => a
     -> (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a)
     -> (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a)
iHat t c e =
   let r = c e
   in if r >= 0 || isNaN r
        then 1  / 0
        else (-1 / auto t) * log( - (c  e))
{-# INLINE iHat #-}

-- | The 'stochasticGradientDescent' function approximates
-- the true gradient of the constFunction by a gradient at
-- a single example. As the algorithm sweeps through the training
-- set, it performs the update for each training example.
--
-- It uses reverse mode automatic differentiation to compute the gradient
-- The learning rate is constant through out, and is set to 0.001
stochasticGradientDescent :: (Traversable f, Fractional a, Ord a)
  => (forall s. Reifies s Tape => e -> f (Reverse s a) -> Reverse s a)
  -> [e]
  -> f a
  -> [f a]
stochasticGradientDescent errorSingle d0 x0 = go xgx0 0.001 dLeft
  where
    dLeft = tail $ cycle d0
    xgx0 = Reverse.gradWith (,) (errorSingle (head d0)) x0
    go xgx !eta d
      | eta ==0       = []
      | otherwise     = x1 : go xgx1 eta (tail d)
      where
        x1 = fmap (\(xi, gxi) -> xi - eta * gxi) xgx
        (_, xgx1) = Reverse.gradWith' (,) (errorSingle (head d)) x1
{-# INLINE stochasticGradientDescent #-}

-- | Perform a gradient descent using reverse mode automatic differentiation to compute the gradient.
gradientAscent :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> [f a]
gradientAscent f = gradientDescent (negate . f)
{-# INLINE gradientAscent #-}

-- | Perform a conjugate gradient descent using reverse mode automatic differentiation to compute the gradient, and using forward-on-forward mode for computing extrema.
--
-- >>> let sq x = x * x
-- >>> let rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)
-- >>> rosenbrock [0,0]
-- 1
-- >>> rosenbrock (conjugateGradientDescent rosenbrock [0, 0] !! 5) < 0.1
-- True
conjugateGradientDescent
  :: (Traversable f, Ord a, Fractional a)
  => (forall s. Chosen s => f (Or s (On (Forward (Forward a))) (Kahn a)) -> Or s (On (Forward (Forward a))) (Kahn a))
  -> f a -> [f a]
conjugateGradientDescent f = conjugateGradientAscent (negate . f)
{-# INLINE conjugateGradientDescent #-}

lfu :: Functor f => (f (Or F a b) -> Or F a b) -> f a -> a
lfu f = runL . f . fmap L

rfu :: Functor f => (f (Or T a b) -> Or T a b) -> f b -> b
rfu f = runR . f . fmap R

-- | Perform a conjugate gradient ascent using reverse mode automatic differentiation to compute the gradient.
conjugateGradientAscent
  :: (Traversable f, Ord a, Fractional a)
  => (forall s. Chosen s => f (Or s (On (Forward (Forward a))) (Kahn a)) -> Or s (On (Forward (Forward a))) (Kahn a))
  -> f a -> [f a]
conjugateGradientAscent f x0 = takeWhile (all (\a -> a == a)) (go x0 d0 d0 delta0)
  where
    dot x y = sum $ zipWithT (*) x y
    d0 = Kahn.grad (rfu f) x0
    delta0 = dot d0 d0
    go xi _ri di deltai = xi : go xi1 ri1 di1 deltai1
      where
        ai = last $ take 20 $ Rank1.extremum (\a -> lfu f $ zipWithT (\x d -> auto x + a * auto d) xi di) 0
        xi1 = zipWithT (\x d -> x + ai*d) xi di
        ri1 = Kahn.grad (rfu f) xi1
        deltai1 = dot ri1 ri1
        bi1 = deltai1 / deltai
        di1 = zipWithT (\r d -> r + bi1 * d) ri1 di
{-# INLINE conjugateGradientAscent #-}

max_stars_repo_name: msakai/ad

max_issues_repo_name: msakai/ad

max_forks_repo_name: msakai/ad


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE TypeOperators       #-}
import           Control.Monad
import           Control.Monad.Random
import           GHC.TypeLits

import qualified Numeric.LinearAlgebra.Static as SA

import           Grenade

{--
  Learning circles
  ================
    This simple example shows how to create a
    network which learns the shape of a circle.
--}

{--
  The Network type
  ================
    For this simple example, we will use a network consisting of
    fully connected layers. There are two hidden layers with 40 neurons
    and 10 neurons respectively. We use the tanh, relu and logit
    activation functions.
    Given an (x, y) coordinate, the output is the probability the
    coordinate is inside the circle.
--}
type FFNet
  = Network
    '[ FullyConnected 2 40, Tanh,
       FullyConnected 40 10, Relu,
       FullyConnected 10 1, Logit
     ]
    '[ 'D1 2,          -- Input
       'D1 40, 'D1 40, -- Tanh fully connected
       'D1 10, 'D1 10, -- Relu fully connected
       'D1 1, 'D1 1    -- Logit fully connected
     ]

{--
  A helper function for displaying the circles.
--}
drawCircle :: FFNet -> IO ()
drawCircle network = do
    let testIns = [ [ (x,y)  | x <- [0..50] ] | y <- [0..20] ]
        outMat  = (map . map) renderCoord testIns
    putStrLn $ unlines outMat

  where
    render n'  | n' <= 0.2  = ' '
               | n' <= 0.4  = '.'
               | n' <= 0.6  = '-'
               | n' <= 0.8  = '='
               | otherwise = '#'

    normx :: S ('D1 1) -> RealNum
    normx (S1D r) = SA.mean r

    renderCoord (x, y) = (render . normx) $ runNet network (S1D $ SA.vector [x / 25 - 1, y / 10 - 1])

{--
  This is the function we are attempting to learn.
  Given a vector (x, y), we check if the distance from the origin
  is less than a radius of 0.5, and return 1 if the coordinate is
  inside the circle and 0 otherwise.
--}
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 => SA.R n -> (SA.R n, RealNum) -> Bool
    u `inCircle` (o, r) = SA.norm_2 (u - o) <= r

{--
  Training and Execution
  ================
    Now we train our network and run it.
--}
main :: IO ()
main = do
    -- Train on 100000 examples, randomly generated with a uniform
    -- distribution.
    let n = 100000
    inps <- replicateM n $ do
      s  <- getRandom
      return $ S1D $ SA.randomVector s SA.Uniform * 2 - 1
    -- For our network, we will use stochastic gradient descent with a batch
    -- size of 1.
    let options = TrainingOptions { optimizer      = OptSGD 0.01 0.9 0.0001
                                  , batchSize      = 1
                                  , validationFreq = 1
                                  , verbose        = Full
                                  , metrics        = [Quadratic]
                                  }
    -- Now we train the network with the fit function.
    -- We pass in a list of tuples of inputs and expected outputs,
    -- and use the quadratic loss function. We train for one epoch.
    -- There is no validation data, hence we pass in an empty list.
    net <- fit (zip inps (map circle inps)) [] options 1 quadratic'

    -- Now we visualise what our network has learned.
    drawCircle net


text:
module Numeric.Morpheus.StatisticsTest
(
  tests
)

where

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

import Numeric.LinearAlgebra
import Numeric.Morpheus.Statistics

x = fromList [35,32,213,54,81,92,63,70]
meanValue = 80
stddevValue = 57.60456

a = matrix 4 [35,32,213,54,
              81,92,63,70,
              77,14,91,20]
aT = tr a

columnMeansRM = fromList [64.33333, 46, 122.33333, 48]
rowMeansRM = fromList [83.5, 76.5, 50.5]
columnMeansCM = rowMeansRM
rowMeansCM = columnMeansRM

columnStddevsRM = fromList [25.48202, 40.84116, 79.75797, 25.53429]
rowStddevsRM = fromList [86.88114, 12.71482, 39.17908]
columnStddevsCM = rowStddevsRM
rowStddevsCM = columnStddevsRM

tests = [ testGroup "Statistics" [
            testCase "mean" $ assertApproxEqual "" 1e-5 meanValue (mean x)
          , testCase "stddev_m" $ assertApproxEqual "" 1e-5 stddevValue (stddev_m meanValue x)
          , testCase "stddev" $ assertApproxEqual "" 1e-5 stddevValue (stddev x)
          ]
        , testGroup "Matrix Statistics" [
            testCase "columnMean - row major" $ assertVector "" 1e-5 columnMeansRM (columnMean a)
          , testCase "columnMean - column major" $ assertVector "" 1e-5 columnMeansCM (columnMean aT)
          , testCase "rowMean - row major" $ assertVector "" 1e-5 rowMeansRM (rowMean a)
          , testCase "rowMean - column major" $ assertVector "" 1e-5 rowMeansCM (rowMean aT)
          , testCase "columnStddev_m - row major" $ assertVector "" 1e-5 columnStddevsRM (columnStddev_m columnMeansRM a)
          , testCase "columnStddev_m - column major" $ assertVector "" 1e-5 columnStddevsCM (columnStddev_m columnMeansCM aT)
          , testCase "rowStddev_m - row major" $ assertVector "" 1e-5 rowStddevsRM (rowStddev_m rowMeansRM a)
          , testCase "rowStddev_m - column major" $ assertVector "" 1e-5 rowStddevsCM (rowStddev_m rowMeansCM aT)
          ]
        ]

max_stars_repo_name: Alexander-Ignatyev/morpheus

max_issues_repo_name: aligusnet/morpheus

max_forks_repo_name: aligusnet/morpheus


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Data.Packed.Vector
-- Copyright   :  (c) Alberto Ruiz 2009
-- License     :  GPL
--
-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
-- Stability   :  provisional
--
-- Random vectors and matrices.
--
-----------------------------------------------------------------------------

module Data.Packed.Random (
    Seed,
    RandDist(..),
    randomVector,
    gaussianSample,
    uniformSample
) where

import Numeric.GSL.Vector
import Data.Packed
import Numeric.ContainerBoot
import Numeric.LinearAlgebra.Algorithms

type Seed = Int

-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
-- Gaussian distribution.
gaussianSample :: Seed
               -> Int -- ^ number of rows
               -> Vector Double -- ^ mean vector
               -> Matrix Double -- ^ covariance matrix
               -> Matrix Double -- ^ result
gaussianSample seed n med cov = m where
    c = dim med
    meds = konst 1 n `outer` med
    rs = reshape c $ randomVector seed Gaussian (c * n)
    m = rs `mXm` cholSH cov `add` meds

-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
-- uniform distribution.
uniformSample :: Seed
               -> Int -- ^ number of rows
               -> [(Double,Double)] -- ^ ranges for each column
               -> Matrix Double -- ^ result
uniformSample seed n rgs = m where
    (as,bs) = unzip rgs
    a = fromList as
    cs = zipWith subtract as bs
    d = dim a
    dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
    am = konst 1 n `outer` a
    m = fromColumns (zipWith scale cs dat) `add` am

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# OPTIONS_GHC -Wall #-}

-- | Field classes
module Plankton.Field
  ( Semifield
  , Field
  , ExpField(..)
  , QuotientField(..)
  , BoundedField(..)
  , infinity
  , neginfinity
  , TrigField(..)
  ) where

import Data.Complex (Complex(..))
import Plankton.Additive
import Plankton.Multiplicative
import Plankton.Ring
import Protolude (Bool, Double, Float, Integer, (||))
import qualified Protolude as P

-- | A Semifield is a Field without Commutative Multiplication.
class (MultiplicativeInvertible a, Ring a) =>
      Semifield a

instance Semifield Double

instance Semifield Float

instance (Semifield a) => Semifield (Complex a)

-- | A Field is a Ring plus additive invertible and multiplicative invertible operations.
--
-- A summary of the rules inherited from super-classes of Field
--
-- > zero + a == a
-- > a + zero == a
-- > (a + b) + c == a + (b + c)
-- > a + b == b + a
-- > a - a = zero
-- > negate a = zero - a
-- > negate a + a = zero
-- > a + negate a = zero
-- > one * a == a
-- > a * one == a
-- > (a * b) * c == a * (b * c)
-- > a * (b + c) == a * b + a * c
-- > (a + b) * c == a * c + b * c
-- > a * zero == zero
-- > zero * a == zero
-- > a * b == b * a
-- > a / a = one
-- > recip a = one / a
-- > recip a * a = one
-- > a * recip a = one
class (AdditiveGroup a, MultiplicativeGroup a, Ring a) =>
      Field a

instance Field Double

instance Field Float

instance (Field a) => Field (Complex a)

-- | A hyperbolic field class
--
-- > sqrt . (**2) == identity
-- > log . exp == identity
-- > for +ive b, a != 0,1: a ** logBase a b ≈ b
class (Field a) =>
      ExpField a where
  exp :: a -> a
  log :: a -> a
  logBase :: a -> a -> a
  logBase a b = log b / log a
  (**) :: a -> a -> a
  (**) a b = exp (log a * b)
  sqrt :: a -> a
  sqrt a = a ** (one / (one + one))

instance ExpField Double where
  exp = P.exp
  log = P.log
  (**) = (P.**)

instance ExpField Float where
  exp = P.exp
  log = P.log
  (**) = (P.**)

-- | todo: bottom is here somewhere???
instance (TrigField a, ExpField a) => ExpField (Complex a) where
  exp (rx :+ ix) = exp rx * cos ix :+ exp rx * sin ix
  log (rx :+ ix) = log (sqrt (rx * rx + ix * ix)) :+ atan2 ix rx

-- | quotient fields explode constraints if they allow for polymorphic integral types
--
-- > a - one < floor a <= a <= ceiling a < a + one
-- > round a == floor (a + one/(one+one))
class (Field a) =>
      QuotientField a where
  round :: a -> Integer
  ceiling :: a -> Integer
  floor :: a -> Integer
  (^^) :: a -> Integer -> a

instance QuotientField Float where
  round = P.round
  ceiling = P.ceiling
  floor = P.floor
  (^^) = (P.^^)

instance QuotientField Double where
  round = P.round
  ceiling = P.ceiling
  floor = P.floor
  (^^) = (P.^^)

-- | A bounded field includes the concepts of infinity and NaN, thus moving away from error throwing.
--
-- > one / zero + infinity == infinity
-- > infinity + a == infinity
-- > isNaN (infinity - infinity)
-- > isNaN (infinity / infinity)
-- > isNaN (nan + a)
-- > zero / zero != nan
--
-- Note the tricky law that, although nan is assigned to zero/zero, they are never-the-less not equal. A committee decided this.
class (Field a) =>
      BoundedField a where
  maxBound :: a
  maxBound = one / zero
  minBound :: a
  minBound = negate (one / zero)
  nan :: a
  nan = zero / zero
  isNaN :: a -> Bool

-- | prints as `Infinity`
infinity :: BoundedField a => a
infinity = maxBound

-- | prints as `-Infinity`
neginfinity :: BoundedField a => a
neginfinity = minBound

instance BoundedField Float where
  isNaN = P.isNaN

instance BoundedField Double where
  isNaN = P.isNaN

-- | todo: work out boundings for complex
-- as it stands now, complex is different eg
--
-- > one / (zero :: Complex Float) == nan
instance (BoundedField a) => BoundedField (Complex a) where
  isNaN (rx :+ ix) = isNaN rx || isNaN ix

-- | Trigonometric Field
class (P.Ord a, Field a) =>
      TrigField a where
  pi :: a
  sin :: a -> a
  cos :: a -> a
  tan :: a -> a
  tan x = sin x / cos x
  asin :: a -> a
  acos :: a -> a
  atan :: a -> a
  sinh :: a -> a
  cosh :: a -> a
  tanh :: a -> a
  tanh x = sinh x / cosh x
  asinh :: a -> a
  acosh :: a -> a
  atanh :: a -> a
  atan2 :: a -> a -> a
  atan2 y x
    | x P.> zero = atan (y / x)
    | x P.== zero P.&& y P.> zero = pi / (one + one)
    | x P.< one P.&& y P.> one = pi + atan (y / x)
    | (x P.<= zero P.&& y P.< zero) || (x P.< zero) =
      negate (atan2 (negate y) x)
    | y P.== zero = pi -- must be after the previous test on zero y
    | x P.== zero P.&& y P.== zero = y -- must be after the other double zero tests
    | P.otherwise = x + y -- x or y is a NaN, return a NaN (via +)

instance TrigField Double where
  pi = P.pi
  sin = P.sin
  cos = P.cos
  asin = P.asin
  acos = P.acos
  atan = P.atan
  sinh = P.sinh
  cosh = P.cosh
  asinh = P.sinh
  acosh = P.acosh
  atanh = P.atanh

instance TrigField Float where
  pi = P.pi
  sin = P.sin
  cos = P.cos
  asin = P.asin
  acos = P.acos
  atan = P.atan
  sinh = P.sinh
  cosh = P.cosh
  asinh = P.sinh
  acosh = P.acosh
  atanh = P.atanh

max_stars_repo_name: chessai/plankton

max_issues_repo_name: chessai/plankton

max_forks_repo_name: chessai/plankton


text:
module Wavefront.Types where

import Numeric.Vector

data WavefrontLine = Comment String
                   | Vertex Double Double Double
                   | TexCoord Double Double
                   | Normal Double Double Double
                   | Face Int Int Int Int Int Int Int Int Int
                   deriving (Show)

isComment :: WavefrontLine -> Bool
isComment (Comment _) = True
isComment _ = False

isVertex :: WavefrontLine -> Bool
isVertex (Vertex _ _ _) = True
isVertex _ = False

isTexCoord :: WavefrontLine -> Bool
isTexCoord (TexCoord _ _) = True
isTexCoord _ = False

isNormal :: WavefrontLine -> Bool
isNormal (Normal _ _ _) = True
isNormal _ = False

isFace :: WavefrontLine -> Bool
isFace (Face _ _ _ _ _ _ _ _ _) = True
isFace _ = False

data Vertex_ = Vertex_ Double Double Double deriving (Show)
data TexCoord_ = TexCoord_ Double Double deriving (Show)
data Normal_ = Normal_ Double Double Double deriving (Show)
data Face_ = Face_ Int Int Int Int Int Int Int Int Int deriving (Show)

vertexFromLine :: WavefrontLine -> Vertex_
vertexFromLine (Vertex a b c) = Vertex_ a b c
vertexFromLine _ = error "vertexFromLine is a partial function: it only accepts vertex lines"

texCoordFromLine :: WavefrontLine -> TexCoord_
texCoordFromLine (TexCoord a b) = TexCoord_ a b
texCoordFromLine _ = error "texCoordFromLine is a partial function: it only accepts tex coord lines"

normalFromLine :: WavefrontLine -> Normal_
normalFromLine (Normal a b c) = Normal_ a b c
normalFromLine _ = error "normalFromLine is a partial function: it only accepts normal lines"

faceFromLine :: WavefrontLine -> Face_
faceFromLine (Face a b c d e f g h i) = Face_ a b c d e f g h i
faceFromLine _ = error "vertexFromLine is a partial function: it only accepts face lines"

vertexToVec :: Vertex_ -> Vec3d
vertexToVec (Vertex_ a b c) = vec3 a b c

texCoordToVec :: TexCoord_ -> Vec2d
texCoordToVec (TexCoord_ a b) = vec2 a b

normalToVec :: Normal_ -> Vec3d
normalToVec (Normal_ a b c) = vec3 a b c


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}

{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE QuasiQuotes #-}

module Vectors where

import Feldspar
import Feldspar.Array.Vector
import Feldspar.Array.Buffered

import Feldspar.Software
import Feldspar.Software as Soft (icompile)
import Feldspar.Software.Compile
import Feldspar.Software.Marshal

import Feldspar.Hardware hiding (Arr, IArr, Ref)
import Feldspar.Hardware as Hard (icompile, icompileSig, icompileAXILite)

import Control.Monad (void)
import Data.Complex (Complex)

-- language-c-quote
import Language.C.Quote.GCC
import qualified Language.C.Syntax as C

-- imperative-edsl
import qualified Language.Embedded.Backend.C  as Imp

import Prelude hiding (take, drop, reverse, length, zip, zipWith, sum, tail, map)

--------------------------------------------------------------------------------
-- * ...
--------------------------------------------------------------------------------

sumLast5 :: SPull (SExp Word32) -> SExp Word32
sumLast5 = sum . take 5 . reverse

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

test1 :: IO ()
test1 = Soft.icompile $ printf "%d" $ sumLast5 inp
  where
    inp :: SPull (SExp Word32)
    inp = 0 ... 10

--------------------------------------------------------------------------------
-- * ...
--------------------------------------------------------------------------------

dot :: (Vector exp, Pully exp vec a, Num a, Syntax exp a) => vec -> vec -> a
dot a b = sum $ zipWith (*) a b

dotArr :: IArr (SExp Int32) -> IArr (SExp Int32) -> SExp Int32
dotArr = dot

dotProg :: Software ()
dotProg = connectStdIO $ return . uncurry dotArr

dotIO :: IO ()
dotIO = Soft.icompile dotProg

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

fir :: (Vector exp, Num a, Syntax exp a) => Pull exp a -> Pull exp a -> Pull exp a
fir coeff = map (dot coeff . reverse) . tail . inits

firArr :: IArr (SExp Int32) -> IArr (SExp Int32) -> Pull SExp (SExp Int32)
firArr a b = fir (toPull a) (toPull b)

firProg :: Software ()
firProg = connectStdIO $ manifestFresh . uncurry firArr

firIO :: IO ()
firIO = Soft.icompile firProg

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

type SStore  a = Store  Software a

printTime_def = [cedecl|
void printTime(typename clock_t start, typename clock_t end)
{
  printf("CPU time (sec): %f\n", (double)(end-start) / CLOCKS_PER_SEC);
}
|]

sizeOf_double_complex :: SExp Length
sizeOf_double_complex = 16

-- | Measure the time for 100 runs of 'fftCore' (excluding initialization) for
-- arrays of the given size
benchmark :: SExp Length -> Software ()
benchmark n = do
  addInclude "<stdio.h>"
  addInclude "<string.h>"
  addInclude "<time.h>"
  addDefinition printTime_def
  start <- newObject "clock_t" False
  end   <- newObject "clock_t" False
  st1 :: SStore (SExp (Complex Double)) <- newStore n
  inp1 <- unsafeFreezeStore n st1
  callProc "memset"
      [ iarrArg (manifest inp1)
      , valArg (0 :: SExp Index)
      , valArg (n*sizeOf_double_complex)
      ]
  st2 :: SStore (SExp (Complex Double)) <- newStore n
  inp2 <- unsafeFreezeStore n st2
  callProc "memset"
      [ iarrArg (manifest inp1)
      , valArg (0 :: SExp Index)
      , valArg (n*sizeOf_double_complex)
      ]
  callProcAssign start "clock" []

  for 0 1 99 $ \(_ :: SExp Index) ->
    void $ manifestFresh (fir (toPull inp1) (toPull inp2))

  callProcAssign end "clock" []
  callProc "printTime" [objArg start, objArg end]

runBenchmark n = runCompiled'
  (Imp.def :: CompilerOpts)
  (Imp.def {Imp.externalFlagsPre = ["-O3"], Imp.externalFlagsPost = ["-lm"]})
  (benchmark n)

printBenchmark n = Soft.icompile (benchmark n)

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

dot_sig :: HSig (SArr Word32 -> SArr Word32 -> Signal Word32 -> ())
dot_sig =
  inputIArr 1 $ \a ->
  inputIArr 1 $ \b ->
  ret $ pure $ dot a b
  -- `dot` is pure, so we lift its result.

test3 :: IO ()
test3 = Hard.icompileSig $ dot_sig
  -- Seems the optimizer isn't happy, might be just my installation tho.

test4 :: IO ()
test4 = Hard.icompileAXILite $ dot_sig

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

dot_mmap :: Software ()
dot_mmap =
  do dot <- mmap "0x43C00000" dot_sig
     a :: Arr (SExp Word32) <- initArr [1]
     b :: Arr (SExp Word32) <- initArr [1]
     c :: Ref (SExp Word32) <- newRef
     --
     call dot (a >>: b >>: c >: nil)
     --
     result <- getRef c
     printf "dot: %d\n" result

test5 :: IO ()
test5 = Soft.icompile dot_mmap

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

max_stars_repo_name: markus-git/co-feldspar

max_issues_repo_name: markus-git/co-feldspar

max_forks_repo_name: markus-git/co-feldspar