text:
{-# LANGUAGE ScopedTypeVariables, TypeFamilies, TupleSections #-}

module Fuml.Supervised.Accuracy where

import Fuml.Core
import Numeric.LinearAlgebra
import Data.List (partition, transpose)

accuracy :: (Eq o) => [(Vector Double, o)] -> Predict p o -> Double
accuracy theData model =
  let isCorrect (v,y) = predict model v == y
      corrects = map isCorrect theData
  in realToFrac (length $ filter id corrects) / realToFrac (length corrects)

rmse :: [(Vector Double, Double)] -> Predict p Double -> Double
rmse theData model =
  let getDev (v,y) = let d = y - predict model v in d*d
      devs = map getDev theData
  in sqrt $ sum devs / realToFrac (length devs)

splitData :: Int -> Int -> [a] -> ([a], [a])
splitData nfolds fold xs
  = let xfs = zip xs (cycle $ [0..(nfolds-1)])
        (zins, zouts) = partition ((==fold) . snd) xfs
    in (map fst zins,map fst zouts)

crossValidate :: Monad m
               => Int
               -> [(Vector Double, o)]
               -> ( [(Vector Double, o)] -> Predict p o1 -> Double )
               -> Supervisor m o p o1
               -> Maybe Int
               -> m Double
crossValidate nfolds theData metric super (Just fold) =
  let (testD, trainD) = splitData nfolds fold theData
      fit = runSupervisor super Nothing trainD
  in fmap (metric testD) fit
crossValidate nfolds theData metric super Nothing =
  let folds = [0..(nfolds-1)]
      accs = mapM (crossValidate nfolds theData metric super . Just) folds
  in fmap avg accs

gridSearch  ::  Monad m
               => Int
               -> [(Vector Double, o)]
               -> ( [(Vector Double, o)] -> Predict p o1 -> Double )
               -> (hyper -> Supervisor m o p o1)
               -> [hyper]
               -> Maybe Int
               -> m [(hyper,Double)]
gridSearch nfolds theData metric hsuper hypers (Just fold) =
  let (testD, trainD) = splitData nfolds fold theData
      getAcc h = fmap ((h,) . metric testD) $ runSupervisor (hsuper h) Nothing $ trainD
  in mapM getAcc hypers
gridSearch nfolds theData metric hsuper hypers Nothing = do
  let folds = [0..(nfolds-1)]
  accs <- mapM (gridSearch nfolds theData metric hsuper hypers . Just) folds
  let f accs = (fst $ head accs, avg $ map snd accs)
  return $ map f $ transpose accs

avg xs = sum xs / realToFrac (length xs)

max_stars_repo_name: futurice/diffusionkinetics-open

max_issues_repo_name: futurice/diffusionkinetics-open

max_forks_repo_name: futurice/diffusionkinetics-open


text:
{-# OPTIONS -Wall #-}

module Types where

import Control.Monad.Except
import Data.Array (Array)
import Data.Complex
import Data.Foldable
import Data.IORef
import Text.ParserCombinators.Parsec
import System.IO (Handle)

data LispVal
  = Atom String
  | List [LispVal]
  | DottedList [LispVal] LispVal
  | Number Integer
  | String String
  | Bool Bool
  | Character Char
  | Float Double
  | Ratio Rational
  | Complex (Complex Double)
  | Vector (Array Int LispVal)
  | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
  | Func
      { params :: [String]
      , vararg :: Maybe String
      , body :: [LispVal]
      , closure :: Env
      }
  | IOFunc ([LispVal] -> IOThrowsError LispVal)
  | Port Handle

instance Show LispVal where
  show = showVal

showVal :: LispVal -> String
showVal (Atom x) = x
showVal (List xs) = "(" ++ unwords (map showVal xs) ++ ")"
showVal (DottedList xs y) = "(" ++ unwords (map showVal xs) ++ " . " ++ showVal y ++ ")"
showVal (Number x) = show x
showVal (String x) = x
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (Character c) = [c]
showVal (Float x) = show x
showVal (Ratio x) = show x
showVal (Complex x) = show x
showVal (Vector xs) = "#(" ++ unwords (map showVal (toList xs)) ++ ")"
showVal (PrimitiveFunc _) = "<primitive>"
showVal Func {params = args, vararg = varargs} =
  "(lambda (" ++
  unwords (map show args) ++
  (case varargs of
     Nothing -> ""
     Just arg -> " . " ++ arg) ++
  ") ...)"
showVal (Port _) = "<IO port>"
showVal (IOFunc _) = "<IO primitive>"

data LispError
  = NumArgs Integer [LispVal]
  | TypeMismatch String LispVal
  | Parser ParseError
  | BadSpecialForm String LispVal
  | NotFunction String String
  | UnboundVar String String
  | Default String

instance Show LispError where
  show = showError

showError :: LispError -> String
showError (UnboundVar message varname) = message ++ ": " ++ varname
showError (BadSpecialForm message form) = message ++ ": " ++ show form
showError (NotFunction message func) = message ++ ": " ++ show func
showError (NumArgs expected found) = "Expected " ++ show expected ++ " args: found values " ++ unwords (map show found)
showError (TypeMismatch expected found) = "Invalid type: expected " ++ expected ++ ", found " ++ show found
showError (Parser parseErr) = "Parse error at " ++ show parseErr
showError (Default s) = s

type IOThrowsError = ExceptT LispError IO

type ThrowsError = Either LispError

type Env = IORef [(String, IORef LispVal)]

ext: hs

max_stars_repo_name: willGuimont/write-yourself-a-scheme

max_issues_repo_name: willGuimont/write-yourself-a-scheme

max_forks_repo_name: willGuimont/write-yourself-a-scheme


text:
{-# LANGUAGE EmptyCase #-}
{-# LANGUAGE TemplateHaskell #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# OPTIONS_GHC -freduction-depth=100 #-}
{-# OPTIONS_GHC -fno-warn-deprecations #-}
-- | Instances for 'Generic' and 'HasMetadata'.
--
-- We define instances for datatypes from @generics-sop@ and
-- @base@ that are supported.
--
-- (There are only instances defined in this module, so the
-- documentation is empty.)
--
module Generics.SOP.Instances () where

-- GHC versions and base versions:
--
-- 7.6.3:  4.6.0.1
-- 7.8.3:  4.7.0.1
-- 7.8.4:  4.7.0.2
-- 7.10.3: 4.8.2.0
-- 8.0.2:  4.9.1.0
-- 8.2.2:  4.10.1.0
-- 8.4.3:  4.11.1.0
-- 8.6.1:  4.12.0.0

import Control.Exception
import Data.Char
import Data.Complex
import Data.Data
import Data.Fixed
import Data.Functor.Compose -- new
import qualified Data.Functor.Const -- new
import Data.Functor.Identity -- new
import Data.Functor.Product -- new
import Data.Functor.Sum -- new
import Data.List.NonEmpty -- new
import qualified Data.Monoid
import Data.Ord
import qualified Data.Semigroup -- new
import Data.Version
import Data.Void -- new
import Foreign.C.Error
import Foreign.C.Types
#if MIN_VERSION_base(4,11,0)
import GHC.ByteOrder -- new
#endif
import GHC.Conc -- new
import GHC.ExecutionStack -- new
import GHC.Exts -- new
-- import GHC.Events -- platform-specific, omitted
import GHC.Fingerprint -- new
import GHC.Float -- new
import qualified GHC.Generics -- new
import GHC.IO.Buffer -- new
import GHC.IO.Device -- new
import GHC.IO.Encoding -- new
import GHC.IO.Encoding.Failure -- new
import GHC.IO.Exception -- new
import GHC.IO.Handle -- new
import GHC.RTS.Flags -- new
import qualified GHC.Stack -- new
import GHC.StaticPtr -- new
import GHC.Stats -- new
import System.Console.GetOpt
import System.IO
import Text.Printf
import Text.Read.Lex

import Generics.SOP.BasicFunctors
import Generics.SOP.Classes
import Generics.SOP.TH

-- Types from Generics.SOP:

deriveGeneric ''I
deriveGeneric ''K
deriveGeneric ''(:.:)
deriveGeneric ''(-.->) -- new

-- Cannot derive instances for Sing
-- Cannot derive instances for Shape
-- Cannot derive instances for NP, NS, POP, SOP
-- Cannot derive instances for metadata types

-- Types from the Prelude:

deriveGeneric ''Bool
deriveGeneric ''Ordering
deriveGeneric ''Maybe
deriveGeneric ''Either
deriveGeneric ''()
deriveGeneric ''(,)              -- 2
deriveGeneric ''(,,)
deriveGeneric ''(,,,)
deriveGeneric ''(,,,,)           -- 5
deriveGeneric ''(,,,,,)
deriveGeneric ''(,,,,,,)
deriveGeneric ''(,,,,,,,)
deriveGeneric ''(,,,,,,,,)
deriveGeneric ''(,,,,,,,,,)      -- 10
deriveGeneric ''(,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,) -- 15
deriveGeneric ''(,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,) -- 20
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,) -- 25
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,,,,,) -- 30
deriveGeneric ''[]

-- Other types from base:

-- From Control.Exception:
deriveGeneric ''IOException
deriveGeneric ''ArithException
deriveGeneric ''ArrayException
deriveGeneric ''AssertionFailed
deriveGeneric ''AsyncException
deriveGeneric ''NonTermination
deriveGeneric ''NestedAtomically
deriveGeneric ''BlockedIndefinitelyOnMVar
deriveGeneric ''BlockedIndefinitelyOnSTM
deriveGeneric ''AllocationLimitExceeded -- new
deriveGeneric ''Deadlock
deriveGeneric ''NoMethodError
deriveGeneric ''PatternMatchFail
deriveGeneric ''RecConError
deriveGeneric ''RecSelError
deriveGeneric ''RecUpdError
deriveGeneric ''ErrorCall
deriveGeneric ''TypeError -- new
deriveGeneric ''MaskingState

-- From Data.Char:
deriveGeneric ''GeneralCategory

-- From Data.Complex:
deriveGeneric ''Complex

-- From Data.Data:
deriveGeneric ''DataRep
deriveGeneric ''Fixity
deriveGeneric ''ConstrRep

-- From Data.Fixed:
deriveGeneric ''Fixed
deriveGeneric ''E0
deriveGeneric ''E1
deriveGeneric ''E2
deriveGeneric ''E3
deriveGeneric ''E6
deriveGeneric ''E9
deriveGeneric ''E12

-- From Data.Functor.Compose
deriveGeneric ''Compose -- new

-- From Data.Functor.Const
deriveGeneric ''Data.Functor.Const.Const -- new

-- From Data.Functor.Identity
deriveGeneric ''Identity -- new

-- From Data.Functor.Product
deriveGeneric ''Product -- new

-- From Data.Functor.Sum
deriveGeneric ''Sum -- new

-- From Data.List.NonEmpty
deriveGeneric ''NonEmpty -- new

-- From Data.Monoid:
deriveGeneric ''Data.Monoid.Dual
deriveGeneric ''Data.Monoid.Endo
deriveGeneric ''Data.Monoid.All
deriveGeneric ''Data.Monoid.Any
deriveGeneric ''Data.Monoid.Sum
deriveGeneric ''Data.Monoid.Product
deriveGeneric ''Data.Monoid.First
deriveGeneric ''Data.Monoid.Last
deriveGeneric ''Data.Monoid.Alt -- new

-- From Data.Ord:
deriveGeneric ''Down

-- From Data.Proxy:
deriveGeneric ''Proxy

-- From Data.Semigroup:
deriveGeneric ''Data.Semigroup.Min -- new
deriveGeneric ''Data.Semigroup.Max -- new
deriveGeneric ''Data.Semigroup.First -- new
deriveGeneric ''Data.Semigroup.Last -- new
deriveGeneric ''Data.Semigroup.WrappedMonoid -- new
deriveGeneric ''Data.Semigroup.Option -- new
deriveGeneric ''Data.Semigroup.Arg -- new

-- From Data.Version:
deriveGeneric ''Version

-- From Data.Void:
deriveGeneric ''Void -- new

-- From Foreign.C.Error:
deriveGeneric ''Errno

-- From Foreign.C.Types:
deriveGeneric ''CChar
deriveGeneric ''CSChar
deriveGeneric ''CUChar
deriveGeneric ''CShort
deriveGeneric ''CUShort
deriveGeneric ''CInt
deriveGeneric ''CUInt
deriveGeneric ''CLong
deriveGeneric ''CULong
deriveGeneric ''CPtrdiff
deriveGeneric ''CSize
deriveGeneric ''CWchar
deriveGeneric ''CSigAtomic
deriveGeneric ''CLLong
deriveGeneric ''CULLong
deriveGeneric ''CIntPtr
deriveGeneric ''CUIntPtr
deriveGeneric ''CIntMax
deriveGeneric ''CUIntMax
deriveGeneric ''CClock
deriveGeneric ''CTime
deriveGeneric ''CUSeconds
deriveGeneric ''CSUSeconds
deriveGeneric ''CFloat
deriveGeneric ''CDouble

#if MIN_VERSION_base(4,11,0)
-- From GHC.ByteOrder:
deriveGeneric ''ByteOrder -- new
#endif

-- From GHC.Conc:
deriveGeneric ''ThreadStatus -- new
deriveGeneric ''BlockReason -- new

-- From GHC.ExecutionStack:
deriveGeneric ''Location -- new
deriveGeneric ''SrcLoc -- new

-- From GHC.Exts:
deriveGeneric ''RuntimeRep -- new
deriveGeneric ''VecCount -- new
deriveGeneric ''VecElem -- new
deriveGeneric ''SpecConstrAnnotation -- new

-- From GHC.Generics:
deriveGeneric ''GHC.Generics.K1 -- new
deriveGeneric ''GHC.Generics.U1 -- new
deriveGeneric ''GHC.Generics.V1 -- new
deriveGeneric ''GHC.Generics.Par1 -- new
deriveGeneric ''GHC.Generics.M1 -- new
deriveGeneric ''GHC.Generics.R -- new
deriveGeneric ''GHC.Generics.S -- new
deriveGeneric ''GHC.Generics.D -- new
deriveGeneric ''GHC.Generics.C -- new
deriveGeneric ''(GHC.Generics.:*:) -- new
deriveGeneric ''(GHC.Generics.:+:) -- new
deriveGeneric ''(GHC.Generics.:.:) -- new
deriveGeneric ''GHC.Generics.Associativity -- new
deriveGeneric ''GHC.Generics.DecidedStrictness -- new
deriveGeneric ''GHC.Generics.SourceStrictness -- new
deriveGeneric ''GHC.Generics.SourceUnpackedness -- new
deriveGeneric ''GHC.Generics.Fixity -- new

-- From GHC.IO.Buffer:
deriveGeneric ''Buffer -- new
deriveGeneric ''BufferState -- new

-- From GHC.IO.Device:
deriveGeneric ''IODeviceType -- new

-- From GHC.IO.Encoding:
deriveGeneric ''BufferCodec -- new
deriveGeneric ''CodingProgress -- new

-- From GHC.IO.Encoding.Failure:
deriveGeneric ''CodingFailureMode -- new

-- From GHC.Fingerprint
deriveGeneric ''Fingerprint -- new

-- From GHC.Float
deriveGeneric ''FFFormat -- new

-- From GHC.IO.Exception:
#if MIN_VERSION_base(4,11,0)
deriveGeneric ''FixIOException -- new
deriveGeneric ''IOErrorType -- new
#endif

-- From GHC.IO.Handle:
deriveGeneric ''HandlePosn -- new
#if MIN_VERSION_base(4,10,0)
deriveGeneric ''LockMode -- new
#endif

-- From GHC.RTS.Flags:
deriveGeneric ''RTSFlags -- new
deriveGeneric ''GiveGCStats -- new
deriveGeneric ''GCFlags -- new
deriveGeneric ''ConcFlags -- new
deriveGeneric ''MiscFlags -- new
deriveGeneric ''DebugFlags -- new
deriveGeneric ''DoCostCentres -- new
deriveGeneric ''CCFlags -- new
deriveGeneric ''DoHeapProfile -- new
deriveGeneric ''ProfFlags -- new
deriveGeneric ''DoTrace -- new
deriveGeneric ''TraceFlags -- new
deriveGeneric ''TickyFlags -- new
#if MIN_VERSION_base(4,10,0)
deriveGeneric ''ParFlags -- new
#endif

-- From GHC.Stack:
deriveGeneric ''GHC.Stack.SrcLoc -- new
deriveGeneric ''GHC.Stack.CallStack -- new

-- From GHC.StaticPtr:
deriveGeneric ''StaticPtrInfo -- new

-- From GHC.Stats:
#if MIN_VERSION_base(4,10,0)
deriveGeneric ''RTSStats -- new
deriveGeneric ''GCDetails -- new
#endif
#if !MIN_VERSION_base(4,11,0)
deriveGeneric ''GCStats -- new
#endif

-- From System.Console.GetOpt:

deriveGeneric ''ArgOrder
deriveGeneric ''OptDescr
deriveGeneric ''ArgDescr

-- From System.Exit:

deriveGeneric ''ExitCode

-- From System.IO:

deriveGeneric ''IOMode
deriveGeneric ''BufferMode
deriveGeneric ''SeekMode
deriveGeneric ''Newline
deriveGeneric ''NewlineMode

-- From Text.Printf:

deriveGeneric ''FieldFormat
deriveGeneric ''FormatAdjustment
deriveGeneric ''FormatSign
deriveGeneric ''FormatParse

-- From Text.Read.Lex:

deriveGeneric ''Lexeme
deriveGeneric ''Number

-- Abstract / primitive datatypes (we don't derive Generic for these):
--
-- Ratio
-- Integer
-- ThreadId
-- Chan
-- MVar
-- QSem
-- QSemN
-- DataType
-- Dynamic
-- IORef
-- TypeRep
-- TyCon
-- TypeRepKey
-- KProxy -- not abstract, but intended for kind-level use
-- STRef
-- Unique
-- ForeignPtr
-- CFile
-- CFpos
-- CJmpBuf
-- Pool
-- Ptr
-- FunPtr
-- IntPtr
-- WordPtr
-- StablePtr
-- Char
-- Double
-- Float
-- Int
-- Int8
-- Int16
-- Int32
-- Int64
-- Word
-- Word8
-- Word16
-- Word32
-- Word64
-- IO
-- ST
-- (->)
-- RealWorld
-- Handle
-- HandlePosn
-- TextEncoding
-- StableName
-- Weak
-- ReadP
-- ReadPrec
-- STM
-- TVar
-- Natural
-- Event
-- EventManager
-- CostCentre
-- CostCentreStack
--
-- Datatypes we cannot currently handle:
--
-- SomeException
-- SomeAsyncException
-- Handler
-- Coercion
-- (:~:)


text:
module ActivFunc (
    sigmoidC,
    sigmoids,
    rampC,
    ramps,
    softmaxC
  ) where

import Numeric.LinearAlgebra

sigmoid x = 1 / (1 + exp (-x))
dsigmoid x = sigmoid x * (1 - sigmoid x)

sigmoidC :: Matrix R -> Matrix R
sigmoidC = cmap sigmoid

dsigmoidC :: Matrix R -> Matrix R
dsigmoidC = cmap dsigmoid

sigmoids = (sigmoidC, dsigmoidC)

ramp x = if x > 0 then x else 0
dramp x = if x > 0 then 1 else 0 -- as a matter of convenience

rampC :: Matrix R -> Matrix R
rampC = cmap ramp

drampC :: Matrix R -> Matrix R
drampC = cmap dramp

ramps = (rampC, drampC)

softmax :: Vector R -> Vector R
softmax us = cmap (/ (sum . toList) es) es
  where es = exp us

softmaxC :: Matrix R -> Matrix R
softmaxC = fromColumns . fmap softmax . toColumns

max_stars_repo_name: pupuu/deep-neuralnet

max_issues_repo_name: pupuu/deep-neuralnet

max_forks_repo_name: pupuu/deep-neuralnet


text:
{-# LANGUAGE NamedFieldPuns, RecordWildCards, RankNTypes #-}

module Train ( Configuration(..)
             , train
             , buildMatrices
             , sigmoid
             , skipgrams
             ) where

import Huffman
import Vocabulary

import Prelude hiding (lookup)
import Control.Monad.ST (ST, runST)
import Control.Monad.ST.Unsafe (unsafeIOToST)
import Control.Arrow ((***))
import Data.List (foldl')
import Data.Maybe (fromJust)
import qualified Data.HashMap.Strict as HM
import qualified Data.IntMap.Strict as IM
import Foreign.Storable (Storable)
import Numeric.LinearAlgebra.Data (Matrix, Vector, (><), (!), (?), toList, asRow, cmap, size)
import Numeric.LinearAlgebra.HMatrix (app, mul, fromList, scale, outer)
import Numeric.LinearAlgebra.Devel (STMatrix (..), runSTMatrix, thawMatrix, writeMatrix, modifyMatrix, mapVectorWithIndexM_, mapMatrixWithIndexM_, unsafeFreezeMatrix, freezeMatrix)
import System.Random (RandomGen, randomRs)

data Configuration = Configuration
                   { vocabulary :: Vocabulary
                   , huffman :: IM.IntMap Encoded
                   , learningRate :: Double
                   , vectorSize :: Int
                   } deriving (Show)

-- | Weight matrix of input to hidden and hidden to output.
type WeightMatrix = (Matrix Double, Matrix Double)

-- TODO: Bias node/row?
-- | Build a weight motrix with a model.
-- | - input to hidden: Voc x Vec
-- | - hidden to output: Vec x (Voc - 1)
buildMatrices :: RandomGen g => g -> Int -> Int -> WeightMatrix
buildMatrices g vocSize vecSize =
  let i2h = (vocSize >< vecSize) [0, 0..]
      -- h2o is transposed.
      bound = 0.5 / fromIntegral vocSize
      h2o = ((vocSize - 1) >< vecSize) $ randomRs (-bound, bound) g
  in (i2h, h2o)

-- | Train a weight matrix with sentences.
train :: Configuration
      -> WeightMatrix
      -> [Sentence]
      -> WeightMatrix
train config (i2h, h2o) sentences = runSTMatrices $ do
  si2h <- thawMatrix i2h
  sh2o <- thawMatrix h2o
  mapM_ (trainSentence config (si2h, sh2o)) sentences
  return (si2h, sh2o)

-- | Train a weight matrix with a sentence.
trainSentence :: Configuration
              -> (STMatrix s Double, STMatrix s Double)
              -> Sentence
              -> ST s ()
trainSentence config@(Configuration {vocabulary, ..}) matrices sentence = do
  let indices = map (fst . fromJust . (`HM.lookup` vocabulary)) sentence
  mapM_ (trainSkipgram config matrices) $ skipgrams 5 indices

trainSkipgram :: Configuration
              -> (STMatrix s Double, STMatrix s Double)
              -> (Int, [Int])
              -> ST s ()
trainSkipgram config matrices (inputIndex, outputIndices) =
  mapM_ (trainPair config matrices inputIndex) outputIndices

trainPair :: Configuration  -- ^ Configuration for training.
          -> (STMatrix s Double, STMatrix s Double) -- ^ Old matrices.
          -> Int            -- ^ Input word's index in the input layer.
          -> Int            -- ^ Output word's index in the output layer.
          -> ST s ()
trainPair Configuration {huffman, learningRate, ..} (si2h, sh2o) inputIndex outputIndex = do
  -- HACK: Using `lookup` because `!` is reserved for `Matrix`.
  let encoded = fromJust $ IM.lookup outputIndex huffman
  -- TODO: Guessing `unsafeFreezeMatrix` doesn't create a new matrix...
  i2h <- unsafeFreezeMatrix si2h
  h2o <- unsafeFreezeMatrix sh2o
      -- TODO: No activate function for the hidden layer?
      -- Vector: vectorSize
  let hidden = i2h ! inputIndex
      -- TODO: This seems to be slow because it's O(vocSize) and copies a new matrix.
      -- TODO: Create `unsafe?` or something.
      -- Matrix: (L(w) - 1 >< vectorSize)
      wo = h2o ? point encoded
      -- Vector: L(w) - 1
      idealOutputs = fromList $ map direction $ code encoded
      -- Vector: L(w) - 1
      errors = cmap sigmoid (app wo hidden) - idealOutputs
      -- Calculate deltas to update weight matrices.
  -- TODO: Does separated `let` in `do` strictly calculate the result?
  let deltaH2O = - learningRate `scale` (errors `outer` hidden)
  let deltaI2H = - learningRate `scale` (asRow errors `mul` wo)

  updateRows sh2o deltaH2O $ point encoded
  updateRows si2h deltaI2H [inputIndex]

updateRows :: STMatrix s Double
           -> Matrix Double
           -> [Int]
           -> ST s ()
updateRows mat delta rowIndices = mapMatrixWithIndexM_ modify delta
  where modify (i, j) d = modifyMatrix mat  (rowIndices !! i) j (+ d)

-- | Sigmoid function.
sigmoid :: Double -> Double
sigmoid x = 1 / (1 + exp (-x))

-- TODO: Randomize the window size (1..size).
-- | Return skipgram pairs from a list.
skipgrams :: Int -> [a] -> [(a, [a])]
skipgrams size = skipgrams' size []
  -- `prev` is a reversed list of previous items.
  where skipgrams' size prev [] = []
        skipgrams' size prev (x:xs) =
          (x, reverse (take size prev) ++ take size xs) : skipgrams' size (x:prev) xs

runSTMatrices :: (Storable t1, Storable t2)
              => (forall s. ST s (STMatrix s t1, STMatrix s t2))
              -> (Matrix t1, Matrix t2)
runSTMatrices st = runST $ do
  (sm1, sm2) <- st
  m1 <- unsafeFreezeMatrix sm1
  m2 <- unsafeFreezeMatrix sm2
  return (m1, m2)

max_stars_repo_name: shuhei/wv

max_issues_repo_name: shuhei/wv

max_forks_repo_name: shuhei/wv


text:
module Data.MyPolynomial (
	Monomial ((:*^:))
	, (*^)
	, Equation(Eq)
	, getEq
	, getL
	, getR
	, Polynomial
	, Roots
	, mbrCoeff
	, mbrPower
	, discriminantQuadratic
	, roots
	, canonify
	, degree
	, handleNegativePowers
	, getABC
	, toMap
	, absurd
	, divergent
	, BroadSol(Absurd, Real)
	, Solution(Quadratic, Simple, Broad)
	, solveEquation
	) where

import Data.Complex (Complex((:+)))
import Data.List
import qualified Data.IntMap.Lazy as M
import Control.Monad.Trans.Reader
import Control.Monad.Trans.State

import Data.MyPolynomial.Type
import Data.MyPolynomial.Parser
import Data.MyPolynomial.Print

type ComplexF = Complex Float
type Roots = (ComplexF, ComplexF)
type NonRoots = [ComplexF]

 -- Sums every Monomial of equal power.
squeeze :: Polynomial -> Polynomial
squeeze poly = runReader (doSqueeze poly) (powers poly)
	where
		 -- Build an unique list of represented powers in @pl@
		powers pl = nub $ fmap mbrPower pl

		-- Every @mn@ in @pl@ whose power is equal to the provided @pw@.
		samePower pw pl = partition (\ (_ :*^: mp) -> mp == pw) pl

		-- For each the provided power @pw@, partition out bros of power p, sum them, add them back to the other partition.
		foldPower pw pl = let	(bros, others) = samePower pw pl
					brogogo = foldr (\ m macc -> head $ macc +^+ m) (zeroP pw) bros
					in (brogogo:others)

		-- Pseudo type annotation :
		-- doSequeeze :: Polynomial -> Reader Env@[Power] Polynomial
		doSqueeze pl = do
			e:es <- ask
			let pl' = foldPower e pl
			  in case es of	[] -> return pl'
					_ -> local (\ (_:en) -> en ) (doSqueeze pl')

 -- Returns a, b, c coefficients as in the canonical expression aX^2 + bX + c
getABC :: M.IntMap Float -> (Float, Float, Float)
getABC mmap = (mmap M.! 2, mmap M.! 1, mmap M.! 0)

 -- Takes a condensed quadratic list polynomial; calls error otherwise.
toMap :: Polynomial -> M.IntMap Float
toMap pl = mapIt (squeeze pl) (M.singleton 0 0)
	where
	  mapIt ((c :*^: p):mn) map = mapIt mn (M.insert p c map)
	  mapIt [] map = map

 -- Calls error if the polynomial degree is out of bounds ([0 <-> 2]).
assertMap :: M.IntMap Float -> ()
assertMap mmap = let d = degree mmap
		   in if d > 2 || d < 0
		   	then	error "Uncondensed or non-quadratic list polynomial."
		   	else	()

 -- Expects a canonical quadratic representation
discriminantQuadratic :: M.IntMap Float -> Float
discriminantQuadratic mmap = let  (a, b, c) = getABC mmap
				in (b ^ 2 - 4.0 * a * c)

data BroadSol = Absurd | Real
data Solution = Quadratic Roots | Simple Float | Broad BroadSol

solveEquation :: Equation -> Solution
solveEquation eq = (doSolve . toMap . getL . canonify) eq
  where
    doSolve mmap = let 	a = mmap M.! 2
    			b = mmap M.! 1
			c = mmap M.! 0
		    in dispatch a b c mmap
    quadratic m = Quadratic ( roots m )
    simple m = Simple ( solveDeg1 m )
    broad m = Broad ( broadSolve m )
    dispatch a b c m	| a /= 0 = quadratic m
			| a == 0 && b /= 0 = simple m
			| a == 0 && b == 0 = broad m

broadSolve :: M.IntMap Float -> BroadSol
broadSolve mmap	= case mmap M.! 0 of	0 -> Real
					_ -> Absurd

solveDeg1 :: M.IntMap Float -> Float
solveDeg1 mmap = let	b = mmap M.! 1
			c = mmap M.! 0
		   in (-c) / b

roots :: M.IntMap Float -> Roots
roots mmap
	| delta < 0 = (cr1, cr2)
	| delta == 0 = (rr0, rr0)
	| delta > 0 = (rr1, rr2)
	where
	 (a, b, c) = getABC mmap
	 delta = discriminantQuadratic mmap
	 cr1 = (-b / (2 * a)) :+ (sqrt(-delta) / (2 * a))
	 cr2 = (-b / (2 * a)) :+ (-sqrt(-delta) / (2 * a))
	 rr1 = ((-b - sqrt(delta)) / (2 * a)) :+ 0
	 rr2 = ((-b + sqrt(delta)) / (2 * a)) :+ 0
	 rr0 = rr1

yankLeft :: Equation -> Equation
yankLeft (Eq (lp, rp)) = Eq (lp ++ (move rp), [zero])
  where
    move [] = []
    move ((c :*^: p):rpn) = ((-c) :*^: p):move(rpn)

canonify :: Equation -> Equation
canonify (Eq (l,r)) = let l' = zero : zeroP 1 : zeroP 2 : l
			  Eq (l'', r') = yankLeft $ Eq (l', r)
			  l''' = squeeze l''
			 in Eq (l''', r')

canonifyP :: Polynomial -> Polynomial
canonifyP pl = let pl' = zero : zeroP 1 : zeroP 2 : pl
		 in squeeze pl'

degree :: M.IntMap Float -> Int
degree mmap = highest 0
	where
	  highest k = case M.lookupGT k mmap of { Nothing -> smaller k;	Just (kn, _) -> highest kn; }
	  smaller k = case M.lookupLE k mmap of { Nothing -> k;		Just (kn, v) -> check kn v; }
	  check kn v = if ( v == 0 && kn > 0 )	then	smaller (kn - 1)
				   		else	kn

cutHighZeroes :: Equation -> Equation
cutHighZeroes (Eq (l, r)) = Eq (skim l, skim r)
	where skim = filter (\(c :*^: p) -> c /= 0 || (c == 0 && p == 0))

absurd :: Equation -> Bool
absurd eq	| Eq ([cl :*^: 0], [cr :*^: 0]) <- eq'  = cl /= cr
		| otherwise = False
	where
	  Eq (l, r) = yankLeft eq
	  eq' = cutHighZeroes $ Eq (squeeze l, squeeze r)

divergent :: Equation -> Bool
divergent eq	| Eq ([cl :*^: 0], [cr :*^: 0]) <- eq'  = cl == cr
		| otherwise = False
	where
	  Eq (l, r) = yankLeft eq
	  eq' = cutHighZeroes $ Eq (squeeze l, squeeze r)

-- Monomial multiplication
infix 5 *^*
(*^*) :: Monomial -> Monomial -> Monomial
c1 :*^: p1 *^* c2 :*^: p2 = (c1 * c2) :*^: (p1 + p2)

infixr 5 /^/
(/^/) :: Monomial -> Monomial -> Monomial
c1 :*^: p1 /^/ c2 :*^: p2 = (c1 * (1.0 / c2)) :*^: (p1 - p2)

-- Monomial sum
infix 4 +^+
(+^+) :: Monomial -> Monomial -> Polynomial
m1@(c1 :*^: p1) +^+ m2@(c2 :*^: p2) | p1 == p2 = [(c1 + c2) :*^: p1]
				    | otherwise = m1:m2:[]

infix 4 -^-
(-^-) :: Monomial -> Monomial -> Polynomial
m1 -^- (c2 :*^: p2) = m1 +^+ ((-c2) :*^: p2)

-- Transforms negative power monomes away and adds forbidden roots to the State.
handleNegativePowers :: Polynomial -> State NonRoots Polynomial
handleNegativePowers p = case (findNeg p) of	Nothing -> get >>= (\s -> state (\_ -> (p, s)))
						Just mn -> handleNext mn p
  where
    bundle (a, b) = nub [a, b]
    addForbidden st nr = nub $ st ++ nr

    handleNext :: Monomial -> Polynomial -> State NonRoots Polynomial
    handleNext mx pl = do
    	get >>= (\s -> put $ addForbidden s (badRoots mx))
	handleNegativePowers ( applyRcp mx pl )

    applyRcp :: Monomial -> Polynomial -> Polynomial
    applyRcp negmn pl = canonifyP $ fmap (/^/ negmn) pl

    badRoots :: Monomial -> [ComplexF]
    badRoots (_ :*^: 0) = [] -- Should not happen.
    badRoots (c :*^: p)	| (-p) == 1 = ((:[]) . (:+ 0) . solveDeg1) $ toMap [zero, (c :*^: 1), zeroP 2]
			| (-p) == 2 = bundle $ roots $ toMap [(c :*^: 2), zeroP 1, zeroP 0]
			| otherwise = error $ "Will not solve: monomial of a higher absolute degree than 2:\n" ++ prettyMonomial (c :*^: p)

    findNeg = find (\(c :*^: p) -> p < 0 && c /= 0)


text:
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE UnicodeSyntax #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DeriveFoldable #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE NoMonomorphismRestriction #-}

module Main where

import Control.Monad ( zipWithM_ )
import Data.Foldable ( toList )
import Data.Monoid ( Monoid(..) )
import Data.SBV
import Data.SBV.Tools.CodeGen
import Data.SBV.Internals ( SolverContext )
import Data.Semigroup ( Semigroup(..) )

import qualified Numeric.GSL.ODE as ODE
import qualified Numeric.LinearAlgebra as LA

import Graphics.Rendering.Chart.Easy hiding ( (.>) )
import Graphics.Rendering.Chart.Backend.Cairo

data Pendulum a = Pendulum
  { pendulumLength :: a
  , pendulumDampingConstant :: a
  , pendulumMass :: a
  , pendulumGravity :: a
  } deriving (Eq, Show, Functor, Foldable, Traversable)

data State a = State
  { stateθ :: a
  , stateω :: a
  } deriving (Eq, Show, Functor, Foldable, Traversable)

newtype Controller a = Controller
  { controllerDamping :: a
  } deriving (Eq, Show, Functor, Foldable, Traversable)

instance EqSymbolic a => EqSymbolic (State a) where
  (State x a) .== (State y b) = x .== y &&& a .== b

stateLabels :: State String
stateLabels = State "θ" "ω"

pendulum :: Fractional a => Pendulum a -> a -> State a -> State a
pendulum sys@(Pendulum len damping _ grav) τ (State θ ω) =
  State ω $
  (grav * taylorSin θ) / len + (damping * ω) / inertia + τ / inertia
  where
    inertia = pendulumInertia sys

pendulumInertia :: Fractional a => Pendulum a -> a
pendulumInertia (Pendulum len _ mass _) = mass * len * len

kineticEnergy :: Fractional a => Pendulum a -> State a -> a
kineticEnergy system (State _ ω) =
  0.5 * pendulumInertia system * ω * ω

dkedt ::
  Fractional a => Controller a -> Pendulum a -> State a -> a
dkedt ctrl system state@(State _ ω) =
  pendulumInertia system * ω * stateω (closedLoop ctrl system state)

-- | Potential energy spans [-2, 0] * mg.
potentialEnergy :: Fractional a => Pendulum a -> State a -> a
potentialEnergy (Pendulum len _ mass grav) (State θ _) =
  len * mass * grav * (taylorCos θ - 1)

dpedt :: Fractional a => Pendulum a -> State a -> a
dpedt (Pendulum len _ mass grav) (State θ ω) =
  len * mass * grav * (- taylorSin θ) * ω

lyapunovController :: Fractional a => Controller a -> Pendulum a -> State a -> a
lyapunovController (Controller kd) (Pendulum len _ mass grav) (State θ ω) =
  (-2) * mass * grav * len * taylorSin θ - kd * ω

closedLoop ::
  Fractional a => Controller a -> Pendulum a -> State a -> State a
closedLoop ctrl system initialState =
  pendulum system torque initialState
  where
    torque = lyapunovController ctrl system initialState

{- Proofs -}

v :: Fractional a => Pendulum a -> State a -> a
v system st =
  kineticEnergy system st - potentialEnergy system st

dvdt :: Fractional a => Controller a -> Pendulum a -> State a -> a
dvdt ctrl system st =
  dkedt ctrl system st - dpedt system st

newtype SAll a = SAll { getSAll :: a }

instance Boolean a => Semigroup (SAll a) where
  (SAll x) <> (SAll y) = SAll $ x &&& y

instance Boolean a => Monoid (SAll a) where
  mempty = SAll true
  mappend = (<>)

allIsPoint :: (IEEEFloating a, Foldable t) => t (SBV a) -> SBool
allIsPoint = getSAll . foldMap (SAll . fpIsPoint)

nanFree ::
  (IEEEFloating a, SolverContext m, Monad m) =>
  (String -> m (SBV a)) -> (State (SBV a) -> SBV a) -> m SBool
nanFree gen f = do
  st <- traverse gen stateLabels
  constrainPi st
  constrainFP st
  pure . fpIsPoint $ f st
  where
    constrainFP = constrain . allIsPoint
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751

lyapunov'sTheorem ::
  ( SymWord a, Fractional a
  , SolverContext m, Monad m) =>
  (String -> m (SBV a)) -> (State (SBV a) -> SBV a) -> (State (SBV a) -> SBV a) -> m SBool
lyapunov'sTheorem gen f dfdt = do
  st <- traverse gen stateLabels
  constrainPi st
--   constrainFP st
--   constrainFP [f st, dfdt st]
  eq <- lyapunovEquilibrium st
  nn <- lyapunovNonNegative st
  gn <- lyapunovGradNegative st
  pure $ eq &&& nn &&& gn
  where
    constrainPi (State θ _) = constrain $ θ .<= π &&& θ .> -π
    π = 3.1415926535897932384626433832795028841971693993751
    lyapunovEquilibrium _ = pure $
      f (State 0 0) .== 0.0

    lyapunovNonNegative st = do
      constrain $ st ./= State 0 0
      pure $ f st .> 0.0

    lyapunovGradNegative st = pure $
      dfdt st .<= 0.0 &&& dfdt (State 0 0) .<= 0.0

nominalController :: Fractional a => Controller a
nominalController = Controller 0.3

nominalSystem :: Fractional a => Pendulum a
nominalSystem = Pendulum 0.5 (-0.03) 0.1 9.81

systemLabels :: Pendulum String
systemLabels = Pendulum "length" "damping" "mass" "gravity"

controllerLabels :: Controller String
controllerLabels = Controller "kd"

proveStability :: IO ThmResult
proveStability =
  prove $ lyapunov'sTheorem sReal v' dvdt'
  where
    v' = v nominalSystem
    dvdt' = dvdt nominalController nominalSystem

proveNanSafety :: IO ThmResult
proveNanSafety =
  prove $ nanFree sFloat controller
  where
    controller = lyapunovController nominalController nominalSystem

-- Simulation

dxdt :: Floating a => p -> [a] -> [a]
dxdt _ [θ, ω] = toList $
  closedLoop nominalController nominalSystem
  (State θ ω)
dxdt _ _ = error "Invalid arguments to 'dxdt'"

dxdtOpenLoop :: Fractional a => p -> [a] -> [a]
dxdtOpenLoop _ [θ, ω] =
  toList $ pendulum nominalSystem 0 (State θ ω)
dxdtOpenLoop _ _ = error "Invalid arguments to 'dxdtOpenLoop'"

solution :: State Double
         -> (Double -> [Double] -> [Double])
         -> LA.Vector Double
         -> LA.Matrix Double
solution state0 f = ODE.odeSolve f $ toList state0

listSolution :: State Double
             -> (Double -> [Double] -> [Double])
             -> [Double]
             -> [(Double, Double)]
listSolution state0 f =
  fmap assoc . LA.toRows . solution state0 f . LA.fromList
  where
    assoc vec = let [t, x] = LA.toList vec in (t, x)

initialStates :: Fractional a => [State a]
initialStates = zipWith State [1e-3, -0.5, 0.3] [1e-3, 0.1, 0.3]

sampleTs :: (Enum a, Fractional a) => [a]
sampleTs = [0, 0.01 .. 7]

makePlot :: PlotValue a
         => String -> String -> String -> [[(a, a)]] -> IO ()
makePlot nm title lbl trajectories =
  toFile opts (nm <> ".png") $ do
  layout_title .= title
  setColors $ fmap opaque [red, blue, green]
  zipWithM_ mkPlot [0 :: Int ..] trajectories
  where
    mkPlot num = plot . line (lbl <> show num) . pure
    opts = def { _fo_size = (1280, 720) }

plotStates :: String -> (Double -> [Double] -> [Double]) -> IO ()
plotStates prefix dynamics = do
  makePlot (prefix <> "_theta") (prefix <> " pendulum angle") "θ" thetas
  makePlot (prefix <> "_omega") (prefix <> " pendulum velocity") "ω" omegas
  where
    trajs = fmap (\st -> listSolution st dynamics sampleTs) initialStates
    withTime sel = zip sampleTs . fmap sel
    thetas, omegas :: [[(Double, Double)]]
    thetas = fmap (withTime fst) trajs
    omegas = fmap (withTime snd) trajs

unstabilized :: IO ()
unstabilized = plotStates "Unstabilized" dxdtOpenLoop

stabilized :: IO ()
stabilized = plotStates "Stabilized" dxdt

main :: IO ()
main = do
  unstabilized
  stabilized
  genCCode

-- Trigonometry

taylorCos :: Fractional a => a -> a
taylorCos x = 1 + sum (take 10 series)
  where
    inc num old =
      let new = old * x * x / (num * (num + 1))
      in new : inc (num + 2) new
    signs = cycle [negate, id]
    series = zipWith ($) signs (inc 1 1)

taylorSin :: Fractional a => a -> a
taylorSin x = x + sum (take 10 series)
  where
    inc num old =
      let new = old * x * x / (num * (num + 1))
      in new : inc (num + 2) new
    signs = cycle [negate, id]
    series = zipWith ($) signs (inc 2 x)

-- C code generation

emitController ::
  (Fractional a, SymWord a) =>
  (String -> SBVCodeGen (SBV a)) -> IO ()
emitController gen = compileToC Nothing "lyapunovController" $ do
  system <- traverse gen systemLabels
  controller <- traverse gen controllerLabels
  state <- traverse gen stateLabels
  cgReturn $ lyapunovController controller system state

genCCode :: IO ()
genCCode = do
  emitController cgGen
  emitTaylor taylorSin "taylorSin" cgGen
  emitTaylor taylorCos "taylorCos" cgGen
  emitCalloutController cgGen
  where
    cgGen :: String -> SBVCodeGen SDouble
    cgGen = cgInput

emitTaylor
  :: (Fractional a, SymWord a)
  => (SBV a -> SBV a)
  -> String
  -> (String -> SBVCodeGen (SBV a))
  -> IO ()
emitTaylor f funName gen = compileToC Nothing funName $
  gen "x" >>= cgReturn . f

taylorSin' :: (Fractional a, SymWord a) => SBV a -> SBV a
taylorSin' = cgUninterpret "taylorSin" mempty taylorSin

taylorCos' :: (Fractional a, SymWord a) => SBV a -> SBV a
taylorCos' = cgUninterpret "taylorCos" mempty taylorCos

lyapunovController'
  :: (SymWord a, Fractional a) =>
     Controller (SBV a)
     -> Pendulum (SBV a) -> State (SBV a) -> SBV a
lyapunovController' (Controller kd) (Pendulum len _ mass grav) (State θ ω) =
  -2 * mass * grav * len * taylorSin' θ + kd * (-ω)

emitCalloutController
  :: (SymWord a, Fractional a) =>
     (String -> SBVCodeGen (SBV a)) -> IO ()
emitCalloutController gen = compileToC Nothing "lyapunovController2" $ do
  system <- traverse gen systemLabels
  controller <- traverse gen controllerLabels
  state <- traverse gen stateLabels
  cgReturn $ lyapunovController' controller system state

max_stars_repo_name: peddie/smt-solving

max_issues_repo_name: peddie/smt-solving

max_forks_repo_name: peddie/smt-solving


text:
{-# language TemplateHaskellQuotes #-}
{-# language FlexibleContexts #-}
{-# language FlexibleInstances #-}
{-# language ScopedTypeVariables #-}
{-# language DataKinds #-}
{-# language TypeOperators #-}
{-# language EmptyCase #-}
{-# language DefaultSignatures #-}
{-# language BangPatterns #-}

module Language.Haskell.TH.TraverseCode
  ( TraverseCode (..)
  , sequenceCode
  , genericTraverseCode
  , genericSequenceCode
  ) where
import Generics.Linear
import Language.Haskell.TH.Syntax (Code, Lift (..), Exp (..), Quote, Name)
import qualified Language.Haskell.TH.Syntax as TH
import Language.Haskell.TH.Lib (conE)

-- for instances
import qualified Data.Functor.Product as FProd
import qualified Data.Functor.Sum as FSum
import Data.Functor.Identity
import qualified Data.Sequence.Internal as Seq
import Data.Coerce
import Control.Applicative
import qualified Data.Semigroup as S
import qualified Data.Monoid as M
import qualified Data.Complex as Complex
import qualified Data.Ord as Ord
import qualified Data.Functor.Compose as Compose
import qualified Data.Tree as Tree
import qualified Data.List.NonEmpty as NonEmpty
import GHC.Tuple (Solo)
import Data.Proxy (Proxy)
import qualified GHC.Generics as GHCGenerics
import qualified Data.Array as Ar
import qualified Data.Primitive.Array as PAr
import qualified Data.Primitive.SmallArray as PSmAr
import Data.Foldable (toList)

-- | Containers supporting \"traversal\" in 'Code'.
class TraverseCode t where
  -- | Given a container and a function to fill it with splices,
  -- produce a splice that will generate a container of their results.
  traverseCode :: Quote m => (a -> Code m b) -> t a -> Code m (t b)

  default traverseCode :: (Quote m, GTraverseCode (Rep1 t), Generic1 t) => (a -> Code m b) -> t a -> Code m (t b)
  traverseCode = genericTraverseCode

-- | Given a container of splices, produce a splice that will generate a
-- container of their results.
sequenceCode :: (TraverseCode t, Quote m) => t (Code m a) -> Code m (t a)
sequenceCode = traverseCode id

-- | Like 'sequenceCode', but using the @"Generics.Linear".'Generic1'@ instance.
genericSequenceCode :: (Quote m, GTraverseCode (Rep1 t), Generic1 t) => t (Code m a) -> Code m (t a)
genericSequenceCode = TH.unsafeCodeCoerce . gtraverseCode id . from1

-- | Like 'traverseCode', but using the @"Generics.Linear".'Generic1'@ instance.
genericTraverseCode :: (Quote m, GTraverseCode (Rep1 t), Generic1 t) => (a -> Code m b) -> t a -> Code m (t b)
genericTraverseCode f = TH.unsafeCodeCoerce . gtraverseCode f . from1

class GTraverseCode f where
  gtraverseCode :: Quote m => (a -> Code m b) -> f a -> m Exp

instance (Datatype c, GTraverseCodeCon f) => GTraverseCode (D1 c f) where
  gtraverseCode f m@(M1 x) = gtraverseCodeCon pkg modl f x
    where
      pkg = packageName m
      modl = moduleName m

class GTraverseCodeCon f where
  gtraverseCodeCon :: Quote m => String -> String -> (a -> Code m b) -> f a -> m Exp

-- This instance seems totally useless, but it's obviously valid.
instance GTraverseCodeCon V1 where
  gtraverseCodeCon _pkg _modl _f x = case x of

instance (GTraverseCodeCon f, GTraverseCodeCon g) => GTraverseCodeCon (f :+: g) where
  gtraverseCodeCon pkg modl f (L1 x) = gtraverseCodeCon pkg modl f x
  gtraverseCodeCon pkg modl f (R1 y) = gtraverseCodeCon pkg modl f y

instance (Constructor c, GTraverseCodeFields f) => GTraverseCodeCon (C1 c f) where
  gtraverseCodeCon pkg modl f m@(M1 x) = gtraverseCodeFields (conE conN) f x
    where
      conBase = conName m
      conN :: Name
      conN = TH.mkNameG_d pkg modl conBase

class GTraverseCodeFields f where
  gtraverseCodeFields :: Quote m => m Exp -> (a -> Code m b) -> f a -> m Exp

instance GTraverseCodeFields f => GTraverseCodeFields (S1 c f) where
  gtraverseCodeFields c f (M1 x) = gtraverseCodeFields c f x

instance (GTraverseCodeFields f, GTraverseCodeFields g) => GTraverseCodeFields (f :*: g) where
  gtraverseCodeFields c f (x :*: y) =
    gtraverseCodeFields (gtraverseCodeFields c f x) f y

instance Lift p => GTraverseCodeFields (K1 i p) where
  gtraverseCodeFields c _f (K1 x) = [| $c x |]

instance GTraverseCodeFields Par1 where
  gtraverseCodeFields cc f (Par1 ca) = [| $cc $(TH.unTypeCode (f ca)) |]

instance GTraverseCodeFields U1 where
  gtraverseCodeFields cc _f U1 = cc

instance (GTraverseCodeFields f, TraverseCode g) => GTraverseCodeFields (f :.: g) where
  gtraverseCodeFields cc f (Comp1 x) =
    gtraverseCodeFields cc (traverseCode f) x

instance GTraverseCodeFields (URec Char) where
  gtraverseCodeFields c _f (UChar ch) = [| $c (UChar ch) |]

instance GTraverseCodeFields (URec Float) where
  gtraverseCodeFields c _f (UFloat ch) = [| $c (UFloat ch) |]

instance GTraverseCodeFields (URec Double) where
  gtraverseCodeFields c _f (UDouble ch) = [| $c (UDouble ch) |]

instance GTraverseCodeFields (URec Int) where
  gtraverseCodeFields c _f (UInt ch) = [| $c (UInt ch) |]

instance GTraverseCodeFields (URec Word) where
  gtraverseCodeFields c _f (UWord ch) = [| $c (UWord ch) |]

instance GTraverseCodeFields f => GTraverseCodeFields (MP1 m f) where
  gtraverseCodeFields c f (MP1 x) = gtraverseCodeFields c f x

-- Would an instance for URec (Ptr ()) make any sense?

-- TraverseCode instances

instance TraverseCode Maybe
instance TraverseCode Identity
instance TraverseCode []
instance TH.Lift a => TraverseCode (Either a)

instance TH.Lift a => TraverseCode ((,) a)
instance (TH.Lift a, TH.Lift b) => TraverseCode ((,,) a b)
instance (TH.Lift a, TH.Lift b, TH.Lift c) => TraverseCode ((,,,) a b c)
instance (TH.Lift a, TH.Lift b, TH.Lift c, TH.Lift d) => TraverseCode ((,,,,) a b c d)
instance (TH.Lift a, TH.Lift b, TH.Lift c, TH.Lift d, TH.Lift e) => TraverseCode ((,,,,,) a b c d e)

instance (TraverseCode f, TraverseCode g) => TraverseCode (FProd.Product f g)
instance (TraverseCode f, TraverseCode g) => TraverseCode (FSum.Sum f g)
instance Lift a => TraverseCode (Const a)

instance TraverseCode V1
instance TraverseCode U1
instance (TraverseCode f, TraverseCode g) => TraverseCode (f :*: g)
instance (TraverseCode f, TraverseCode g) => TraverseCode (f :+: g)
instance TraverseCode f => TraverseCode (M1 i c f)
instance TraverseCode Par1
instance Lift a => TraverseCode (K1 i a)
instance (TraverseCode f, TraverseCode g) => TraverseCode (f :.: g)
instance (TraverseCode f, TraverseCode g) => TraverseCode (Compose.Compose f g)
instance TraverseCode f => TraverseCode (MP1 m f)
instance TraverseCode f => TraverseCode (GHCGenerics.Rec1 f)
instance (TraverseCode f, TraverseCode g) => TraverseCode (f GHCGenerics.:.: g)

-- The Elem instance isn't needed for the Seq instance
instance TraverseCode Seq.Elem
instance TraverseCode Seq.Digit
instance TraverseCode Seq.Node
instance TraverseCode Seq.FingerTree
instance TraverseCode Seq.Seq where
  -- This wonky way of doing it makes for a more compact
  -- splice.
  traverseCode f s = [|| coerceFT $$(traverseCode f ft') ||]
    where
      ft' = coerceSeq (Seq.zipWith (flip const) (Seq.replicate (Seq.length s) ()) s)

coerceFT :: Seq.FingerTree a -> Seq.Seq a
coerceFT = coerce
coerceSeq :: Seq.Seq a -> Seq.FingerTree a
coerceSeq = coerce

instance TraverseCode Seq.ViewL
instance TraverseCode Seq.ViewR

instance TraverseCode ZipList
instance TraverseCode Complex.Complex
instance TraverseCode S.First
instance TraverseCode M.First
instance TraverseCode S.Last
instance TraverseCode M.Last
instance TraverseCode S.Min
instance TraverseCode S.Max
instance TraverseCode Ord.Down
instance TraverseCode S.WrappedMonoid
instance TraverseCode S.Dual
instance TraverseCode S.Product
instance TraverseCode S.Sum
instance TraverseCode Solo
instance TraverseCode Tree.Tree
instance TraverseCode NonEmpty.NonEmpty
-- instance TraverseCode TH.TyVarBndr  (Lift Name isn't in base)
instance TraverseCode m => TraverseCode (WrappedMonad m)
instance TraverseCode (p a) => TraverseCode (WrappedArrow p a)
instance TH.Lift a => TraverseCode (S.Arg a)
instance TraverseCode Proxy
instance TraverseCode f => TraverseCode (M.Ap f)
instance TraverseCode f => TraverseCode (M.Alt f)
instance TraverseCode (URec Char)
instance TraverseCode (URec Double)
instance TraverseCode (URec Float)
instance TraverseCode (URec Int)
instance TraverseCode (URec Word)

instance (Ar.Ix i, TH.Lift i) => TraverseCode (Ar.Array i) where
  traverseCode f xs = [|| Ar.listArray bnds $$csc ||]
    where
      csc = traverseCode f (Ar.elems xs)
      bnds = Ar.bounds xs

instance TraverseCode PAr.Array where
  traverseCode f ary = case lst of
      [] -> [|| PAr.emptyArray ||]
      [x] -> [|| pure $$(f x) ||]
      x : xs -> [|| unsafeArrayFromListN len $$(f x) $$(traverseCode f xs) ||]
    where
      len = length ary
      lst = toList ary

instance TraverseCode PSmAr.SmallArray where
  traverseCode f ary = case lst of
      [] -> [|| PSmAr.emptySmallArray ||]
      [x] -> [|| pure $$(f x) ||]
      x : xs -> [|| unsafeSmallArrayFromListN len $$(f x) $$(traverseCode f xs) ||]
    where
      len = length ary
      lst = toList ary

-- | Strictly create an array from a nonempty list (represented as
-- a first element and a list of the rest) of a known length. If the length
-- of the list does not match the given length, this makes demons fly
-- out of your nose.
unsafeArrayFromListN :: Int -> a -> [a] -> PAr.Array a
unsafeArrayFromListN n y ys =
  PAr.createArray n y $ \ma ->
    let go !_ix [] = pure ()
        go !ix (x : xs) = do
            PAr.writeArray ma ix x
            go (ix+1) xs
    in go 1 ys

-- | Strictly create an array from a nonempty list (represented as
-- a first element and a list of the rest) of a known length. If the length
-- of the list does not match the given length, this makes demons fly
-- out of your nose.
unsafeSmallArrayFromListN :: Int -> a -> [a] -> PSmAr.SmallArray a
unsafeSmallArrayFromListN n y ys =
  PSmAr.createSmallArray n y $ \ma ->
    let go !_ix [] = pure ()
        go !ix (x : xs) = do
            PSmAr.writeSmallArray ma ix x
            go (ix+1) xs
    in go 1 ys

max_stars_repo_name: treeowl/traverse-code

max_issues_repo_name: treeowl/traverse-code

max_forks_repo_name: treeowl/traverse-code


text:
{-# LANGUAGE NoMonomorphismRestriction, QuasiQuotes #-}
import Control.Arrow
import Control.Monad
import NeuralNetwork
import Numeric.LinearAlgebra.HMatrix
import Text.Printf.TH
import VisualizeFunction

andDataSet = map (vector *** vector) [
    ([0,0], [0]),
    ([0,1], [0]),
    ([1,0], [0]),
    ([1,1], [1])
    ]

xorDataSet = map (vector *** vector) [
    ([0,0], [0]),
    ([0,1], [1]),
    ([1,0], [1]),
    ([1,1], [0])
    ]

nn = randomlyWeightedNetwork 0 [2, 9, 1] logisticAF
trainOn dataset = iterate (stochasticGradientDescent dataset 0.25) nn

andNNs = trainOn andDataSet
xorNNs = trainOn xorDataSet

outputs dataset nns = map (\nn -> map (applyNN nn) (map fst dataset)) nns
andOutputs = outputs andDataSet andNNs
xorOutputs = outputs xorDataSet xorNNs

image = visualizeFunctionBMP (0,0) (1,1) (1024,1024) (0,1) . ((!0).) . applyNN

main = forM_ [0..10] $ \i -> do
    saveImage ([s|neuralnetwork_and_%04d.bmp|] (2^i)) (image (andNNs !! (2^i)))
    saveImage ([s|neuralnetwork_xor_%04d.bmp|] (2^i)) (image (xorNNs !! (2^i)))


text:
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.LAPACK

m :: Matrix Double
m = (3><3)
  [7.589183,1.703609,-4.477162,
    -4.597851,9.434889,-6.543450,
    0.4588202,-6.115153,1.331191]

v :: Matrix Double
v = (3><1)
  [1.745005,-4.448092,-4.160842]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
{-# LANGUAGE DataKinds            #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE GADTs                #-}
{-# LANGUAGE PatternSynonyms      #-}
{-# LANGUAGE RankNTypes           #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE StandaloneDeriving   #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE TypeOperators        #-}
{-# LANGUAGE UndecidableInstances #-}
{-# OPTIONS_GHC -Wall #-}
-- | This module is wrongly named.
module VectorSpace (
    LinMap (.., LI),
    lmul,
    HasDim(Dim, dimDict),
    toRawMatrix,
    evalL,
    L (..),
    linear,
    VectorSpace (..),
    toVector,
    fromVector,
) where

import Data.Constraint       (Dict (..), withDict, (:-))
import Data.Proxy            (Proxy (..))
import GHC.TypeLits
import Overloaded.Categories

import qualified Control.Category      as C
import qualified Data.Constraint.Nat   as C
import qualified Numeric.LinearAlgebra as L

-- import qualified Numeric.LinearAlgebra.Static as LS

data LinMap a b where
    LZ :: LinMap a b
    LD :: Double -> LinMap a a
    LH :: LinMap a b -> LinMap a c -> LinMap a (b, c)
    LV :: LinMap a c -> LinMap b c -> LinMap (a, b) c
    LA :: LinMap a b -> LinMap a b -> LinMap a b

deriving instance Show (LinMap a b)

pattern LI :: forall a b. () => (b ~ a) => LinMap a b
pattern LI = LD 1

lmul :: Double -> LinMap a b -> LinMap a b
lmul _ LZ       = LZ
lmul k (LD x)   = LD (k * x)
lmul k (LH f g) = LH (lmul k f) (lmul k g)
lmul k (LV f g) = LV (lmul k f) (lmul k g)
lmul k (LA f g) = LA (lmul k f) (lmul k g)

lcomp :: LinMap b c -> LinMap a b -> LinMap a c
lcomp LZ       _        = LZ
lcomp _        LZ       = LZ
lcomp (LD k)   h        = lmul k h
lcomp h        (LD k)   = lmul k h
lcomp (LA f g) h        = LA (lcomp f h) (lcomp g h)
lcomp f        (LA g h) = LA (lcomp f g) (lcomp f h)
lcomp (LH f g) h        = LH (lcomp f h) (lcomp g h)
lcomp h        (LV f g) = LV (lcomp h f) (lcomp h g)
lcomp (LV f g) (LH u v) = LA (lcomp f u) (lcomp g v)

instance Category LinMap where
    id  = LI
    (.) = lcomp

instance CategoryWith1 LinMap where
    type Terminal LinMap = ()
    terminal = LZ

instance CartesianCategory LinMap where
    type Product LinMap = (,)
    proj1  = LV C.id LZ
    proj2  = LV LZ C.id
    fanout = LH

instance CategoryWith0 LinMap where
    type Initial LinMap = ()
    initial = LZ

instance CocartesianCategory LinMap where
    type Coproduct LinMap = (,)
    inl   = LH C.id LZ
    inr   = LH LZ C.id
    fanin = LV

instance BicartesianCategory LinMap where
    distr = LH
        (LH (LV (LV LI LZ) LZ) (LV LZ LI))
        (LH (LV (LV LZ LI) LZ) (LV LZ LI))

newtype L a b = L (forall r. LinMap r a -> LinMap r b)

lfst :: LinMap a (b, c) -> LinMap a b
lfst (LA f g) = LA (lfst f) (lfst g)
lfst (LH f _) = f
lfst (LV f g) = LV (lfst f) (lfst g)
lfst LZ       = LZ
lfst (LD k)   = LV (LD k) LZ

lsnd :: LinMap a (b, c) -> LinMap a c
lsnd (LH _ g) = g
lsnd (LA f g) = LA (lsnd f) (lsnd g)
lsnd (LV f g) = LV (lsnd f) (lsnd g)
lsnd LZ       = LZ
lsnd (LD k)   = LV LZ (LD k)

linitial :: LinMap r () -> LinMap r a
linitial _ = LZ

linear :: Double -> L a a
linear k = L $ lmul k

-- lmult :: Double -> Double -> LinMap r (a, a) -> LinMap r a
-- lmult x y (LH f g) = LA (LK y f) (LK x g)
-- lmult x y (LV f g) = LV (lmult x y f) (lmult x y g)
-- lmult x y (LA f g) = LA (lmult x y f) (lmult x y g)
-- lmult x y (LK k f) = LK k (lmult x y f)
-- lmult _ _ LZ       = LZ
-- lmult x y LI       = LV (LK y LI) (LK x LI)

instance Category L where
    id = L id
    L f . L g = L (f . g)

instance CategoryWith1 L where
    type Terminal L = ()

    terminal = L (\_ -> LZ)

instance CartesianCategory L where
    type Product L = (,)

    proj1 = L lfst
    proj2 = L lsnd

    fanout (L f) (L g) = L $ \x -> LH (f x) (g x)

instance CategoryWith0 L where
    type Initial L = ()

    initial = L linitial

-- Is this correct?
instance CocartesianCategory L where
    type Coproduct L = (,)

    inl = L $ \f -> LH f LZ
    inr = L $ \g -> LH LZ g

    fanin (L f) (L g) = L $ \x -> LA (f (lfst x)) (g (lsnd x))

class HasDim a where
    type Dim a :: Nat

    dimDict :: Proxy a -> Dict (KnownNat (Dim a))

    splitPair :: (a ~ (b, c)) => (Dict (HasDim b), Dict (HasDim c))
    splitPair = error "impossible: splitPair"

instance HasDim () where
    type Dim () = 0
    dimDict _ = Dict

instance HasDim Double where
    type Dim Double = 1
    dimDict _ = Dict

instance (HasDim a, HasDim b) => HasDim (a, b) where
    type Dim (a, b) = Dim a + Dim b

    dimDict _ =
        withDimDict (Proxy :: Proxy a) $
        withDimDict (Proxy :: Proxy b) $
        withDict (C.plusNat :: (KnownNat (Dim a), KnownNat (Dim b)) :- KnownNat (Dim a + Dim b))
        Dict

    splitPair = (Dict, Dict)

withDimDict :: HasDim a => Proxy a -> (KnownNat (Dim a) => r) -> r
withDimDict p = withDict (dimDict p)

dim :: forall a. HasDim a => Proxy a -> Int
dim p = withDimDict p $ fromInteger $ natVal (Proxy :: Proxy (Dim a))

toRawMatrix :: forall a b. (HasDim a, HasDim b) => LinMap a b -> L.Matrix Double
toRawMatrix LZ       = (dim (Proxy :: Proxy a) L.>< dim (Proxy :: Proxy b)) (repeat 0)
toRawMatrix (LD k)   = L.scale k (L.ident (dim (Proxy :: Proxy a)))
toRawMatrix (LA f g) = L.add (toRawMatrix f) (toRawMatrix g)
toRawMatrix (LH f g) = go splitPair f g where
    go :: (Dict (HasDim x), Dict (HasDim y)) -> LinMap a x -> LinMap a y -> L.Matrix Double
    go (Dict, Dict) f' g' = toRawMatrix f' L.||| toRawMatrix g'
toRawMatrix (LV f g) = go splitPair f g where
    go :: (Dict (HasDim x), Dict (HasDim y)) -> LinMap x b -> LinMap y b -> L.Matrix Double
    go (Dict, Dict) f' g' = toRawMatrix f' L.=== toRawMatrix g'

evalL :: (HasDim a, HasDim b) => L a b -> L.Matrix Double
evalL (L f) = toRawMatrix (f (LD 1))

-- toStaticMatrix :: forall a b. (HasDim a, HasDim b) => LinMap a b -> LS.L (Dim a) (Dim b)
-- toStaticMatrix LZ =
--     withDimDict (Proxy :: Proxy a) $
--     withDimDict (Proxy :: Proxy b) 0
-- toStaticMatrix LI =
--     withDimDict (Proxy :: Proxy a) LS.eye
-- toStaticMatrix (LA f g) =
--     withDimDict (Proxy :: Proxy a) $
--     withDimDict (Proxy :: Proxy b) $
--     L.add (toStaticMatrix f) (toStaticMatrix g)
-- toStaticMatrix (LK k f) =
--     withDimDict (Proxy :: Proxy a) $
--     withDimDict (Proxy :: Proxy b) $
--     toStaticMatrix f LS.<> LS.diag (LS.konst k)
-- toStaticMatrix (LH f g) = go splitPair f g where
--     go :: forall x y. (x,y) ~ b => (Dict (HasDim x), Dict (HasDim y)) -> LinMap a x -> LinMap a y -> LS.L (Dim a) (Dim x + Dim y)
--     go (Dict, Dict) f' g' =
--         withDimDict (Proxy :: Proxy a) $
--         withDimDict (Proxy :: Proxy b) $
--         withDimDict (Proxy :: Proxy x) $
--         withDimDict (Proxy :: Proxy y) $
--         toStaticMatrix f' LS.||| toStaticMatrix g'
-- toStaticMatrix (LV f g) = go splitPair f g where
--     go :: forall x y. (x,y) ~ a => (Dict (HasDim x), Dict (HasDim y)) -> LinMap x b -> LinMap y b -> LS.L (Dim x + Dim y) (Dim b)
--     go (Dict, Dict) f' g' =
--         withDimDict (Proxy :: Proxy a) $
--         withDimDict (Proxy :: Proxy b) $
--         withDimDict (Proxy :: Proxy x) $
--         withDimDict (Proxy :: Proxy y) $
--         toStaticMatrix f' LS.=== toStaticMatrix g'

-------------------------------------------------------------------------------
-- Vector space
-------------------------------------------------------------------------------

class HasDim a => VectorSpace a where
    toVector' :: a -> [Double] -> [Double]

    fromVector' :: [Double] -> (a -> [Double] -> r) -> r

toVector :: VectorSpace a => a -> [Double]
toVector x = toVector' x []

fromVector :: VectorSpace a => [Double] -> a
fromVector ds = fromVector' ds const

instance VectorSpace Double where
    toVector' d = (d :)

    fromVector' []     k = k 0 []
    fromVector' (d:ds) k = k d ds

instance (VectorSpace a, VectorSpace b) => VectorSpace (a, b) where
    toVector' (a, b) = toVector' a . toVector' b

    fromVector' xs k =
        fromVector' xs $ \a ys ->
        fromVector' ys $ \b zs ->
        k (a, b) zs

max_stars_repo_name: Mikolaj/overloaded

max_issues_repo_name: Mikolaj/overloaded

max_forks_repo_name: Mikolaj/overloaded


text:
module MachineLearning.Optimization.GradientDescentTest
(
  tests
)

where

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

import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA

import qualified MachineLearning as ML
import MachineLearning.Regularization (Regularization(..))
import MachineLearning.LeastSquaresModel (LeastSquaresModel(..))
import MachineLearning.Optimization.GradientDescent

import MachineLearning.DataSets (dataset1)

(x, y) = ML.splitToXY dataset1

muSigma = ML.meanStddev x
xNorm = ML.featureNormalization muSigma x
x1 = ML.addBiasDimension xNorm
initialTheta = LA.konst 0 (LA.cols x1)
-- Normal Equation's Result: 340412.660,110631.050,-6649.474
lsExpectedTheta = LA.vector [340412.660, 110630.886, -6649.310]
eps = 1e-3

isInDescendingOrder :: [Double] -> Bool
isInDescendingOrder lst = and . snd . unzip $ scanl (\(prev, _) current -> (current, prev >= current)) (1/0, True) lst

testGradientDescent model expectedTheta = do
  let (theta, optPath) = gradientDescent 0.01 model eps 5000 RegNone x1 y initialTheta
      js = V.toList $ (LA.toColumns optPath) !! 1
  assertVector "theta" 0.01 expectedTheta theta
  assertBool "non-increasing errors" $ isInDescendingOrder js

tests = [testGroup "gradientDescent" [
            testCase "leastSquares" $ testGradientDescent LeastSquares lsExpectedTheta
            ]
        ]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module VoiceLeading.Theory where

import           VoiceLeading.Base
import           VoiceLeading.Helpers           ( rotate )

import           Data.List                      ( elemIndex )
import           Data.Maybe                     ( mapMaybe )
import           Data.Ratio
import           Data.Function.Memoize          ( memoize )

import           Data.Aeson
import qualified Data.Map                      as M
import qualified Data.Vector                   as V
import qualified Data.Vector.Unboxed           as VU
import qualified Data.Vector.Unboxed.Mutable   as VUM

import           Numeric.SpecFunctions          ( logBeta )
import qualified Control.Loop                  as L

-- simple ad-hoc heuristic
--------------------------

modal :: Int -> Int -> Int
modal i m = [0, 2, 4, 5, 7, 9, 11] !! mod (i + m) 7

cMajor :: [Int]
cMajor = [2, 2, 1, 2, 2, 2, 1]

scale :: KeySig -> [Int]
scale = memoize s
 where
  s (KeySig root mode) =
    scanl (\a b -> (a + b) `mod` 12) root (init $ rotate mode cMajor)

leadingTone :: KeySig -> Int
leadingTone = memoize lt where lt (KeySig root _) = (root - 1) `mod` 12

majorChord :: [Int]
majorChord = [0, 4, 7, 10, 11, 2, 1, 5, 6, 9, 8, 3]

minorChord :: [Int]
minorChord = [0, 3, 7, 10, 2, 11, 1, 5, 6, 8, 9, 4]

-- | Returns a pair of the estimated root and the "chordness" of the given set of pitches
findHarm :: [Int] -> (Int, Double)
findHarm []      = (0, 0)
findHarm pitches = testRoot [0 .. 11]
 where
  ps = map (`mod` 12) pitches
  indexVal :: Int -> Rational
  indexVal i = 1 % (fromIntegral i + 1)
  testChord chord tps = sum $ map val [0 .. maximum is]
   where
    is = mapMaybe (`elemIndex` chord) tps
    val i = if i `elem` is then indexVal i else negate (indexVal i)
  testRoot [] = (0, 0) -- just a catchall
  testRoot (r : rs) =
    let tps   = map ((`mod` 12) . (\x -> x - r)) ps
        score = max (testChord majorChord tps) (testChord minorChord tps)
    in  testRoot' rs r score
  testRoot' []       best score = (best, fromRational score)
  testRoot' (r : rs) best score = if newScore > score
    then testRoot' rs r newScore
    else testRoot' rs best score
   where
    tps      = map ((`mod` 12) . (\x -> x - r)) ps
    newScore = max (testChord majorChord tps) (testChord minorChord tps)

-- chord profiles
-----------------

type Profile = (Double, VU.Vector Double)
type Profiles = V.Vector Profile

loadProfiles :: FilePath -> IO (Maybe (M.Map String (M.Map Int Double)))
loadProfiles = decodeFileStrict'

-- | Returns a vector of profiles.
-- Each profile is a pair of the total counts in the profile and the smoothed count vector.
vectorizeProfiles :: M.Map String (M.Map Int Double) -> Profiles
vectorizeProfiles profmap = V.fromList $ profList <$> M.elems profmap
 where
  profList pm = (VU.sum lst, lst)
    where lst = VU.fromList [ 1 + pm M.! i | i <- [0 .. 11] ]

logDirichletMultinomial :: Profile -> VU.Vector Double -> Double -> Double
logDirichletMultinomial (a0, as) xs n = num - denom
 where
  l = VU.length xs
  entry acc i = if xi > 0 then acc + log xi + logBeta (as VU.! i) xi else acc
    where xi = xs VU.! i
  num   = log n + logBeta a0 n
  denom = L.forLoopFold 0 (< l) (+ 1) 0 entry

mkChord :: Int -> [Int] -> VU.Vector Double
mkChord r pitches = VU.create $ do
  chord <- VUM.replicate 12 0
  mapM_ (VUM.unsafeModify chord (+ 1) . (`mod` 12) . subtract r) pitches
  pure chord

matchChordProfiles :: Profiles -> [Int] -> (Int, Double)
matchChordProfiles _        []      = (0, 0)
matchChordProfiles profiles pitches = (maxRoot, chordness / 0.3) -- normalize by maximum value
 where
  roots = V.fromListN 12 [0 .. 11]
  n     = fromIntegral $ length pitches
  match chord prof = logDirichletMultinomial prof chord n
  bestMatch r = V.maximum $ match (mkChord r pitches) <$> profiles
  bestMatches = bestMatch <$> roots
  maxRoot     = V.maxIndex bestMatches
  chordness   = exp $ bestMatches V.! maxRoot

max_stars_repo_name: DCMLab/part-writing

max_issues_repo_name: DCMLab/part-writing

max_forks_repo_name: DCMLab/part-writing


text:
{-# LANGUAGE BangPatterns   #-}
{-# LANGUAGE NamedFieldPuns #-}
module Network
    ( NeuralNet (..),
      train,
      test,
      initWB,
      loadData,
      loadTestData,
    ) where

import           Codec.Compression.GZip (decompress)
import           Control.Monad
import qualified Data.Array.IO          as A
import qualified Data.ByteString        as BS
import qualified Data.ByteString.Lazy   as BL
import           Data.Int               (Int64)
import           Data.List
import           GHC.Base               (build)
import           Numeric.LinearAlgebra  hiding (build, normalize)
import           Prelude                hiding ((<>))
import           System.Directory       (getCurrentDirectory)
import qualified System.Random          as R
import           Util

type Image        = Vector Double
type Label        = Vector Double
type Bias         = Vector Double
type Weight       = Matrix Double
type Activation   = Matrix Double
type Nabla        = Matrix Double
type ImagesMat    = Matrix Double
type LabelsMat    = Matrix Double

data NeuralNet = NN
    { weights   :: ![Weight]   -- ^ List of all weights
    , biases    :: ![Bias]    -- ^ List of all biases
    , eta       :: !Double    -- ^ Learning rate
    , epochs    :: !Int       -- ^ No. of epochs
    , layers    :: !Int       -- ^ No. of layers
    , layerSize :: !Int       -- ^ Size of hidden layers
    , batchSize :: !Int       -- ^ Size of mini-batch
    , trainData :: ![(Image, Label)]
    , testData  :: ![(ImagesMat, Int)]
    } deriving (Show, Eq)

test :: NeuralNet -> IO ()
test net@NN{epochs, biases, weights, testData} = do
    let (imgs, labels) = unzip testData
    let guesses = [guess img (zip weights $ toMatrix <$> biases) | img <- imgs]
    let c = sum $ (\(a,b) -> if a == b then 1 else 0) <$> zip guesses labels
    putStrLn $ "Epoch #" ++ (show $ 30-epochs) ++ ": " ++ (show $ c/100) ++ "%"
  where
    guess :: ImagesMat -> [(Weight, Matrix Double)] -> Int
    guess i wb = maxIndex . flatten . head $ feedforward i wb

-- | Fully matrix-based approach to backpropagation.
train :: NeuralNet -> IO NeuralNet
train net@NN{weights, batchSize, epochs, trainData} = do
    shuffledData <- shuffle trainData
    let miniBatches = chunkList batchSize shuffledData
    let batches = fmap getMiniBatch miniBatches
    let !(newWs, newBs) = (recurTest batchSize (weights, biases net) batches) :: ([Weight], [Bias])
    let newNet = net { weights = newWs, biases = newBs, epochs = epochs-1 }
    test newNet
    case epochs of
        1 -> return newNet
        _ -> train newNet

recurTest :: Int -> ([Weight], [Bias]) -> [(ImagesMat, LabelsMat)] -> ([Weight], [Bias])
recurTest batchSize !wb [batch] = trainBatch batchSize wb batch
recurTest batchSize !wb (firstBatch:batches) = recurTest batchSize (trainBatch batchSize wb firstBatch) batches

trainBatch :: Int -> ([Weight], [Bias]) -> (ImagesMat, LabelsMat) -> ([Weight], [Bias])
trainBatch batchSize !(ws, bs) (mX, mY) = backprop mY ws bs $ feedforward mX (zip ws $ bTm bs batchSize)

-- | Make biases suitable for full-matrix backprop.
bTm :: [Bias] -> Int -> [Matrix Double]
bTm biases batchSize = fmap (fromColumns . replicate batchSize) biases

-- | Turn mini-batch matrices into one matrix.
getMiniBatch :: [(Image, Label)] -> (ImagesMat, LabelsMat)
getMiniBatch mBatchData = (fromColumns imgs, fromColumns labels)
  where
    (imgs, labels) = unzip mBatchData

feedforward :: ImagesMat -> [(Weight, Matrix Double)] -> [Activation]
feedforward = scanr (\(w, b) a -> sigmoid $ w <> a + b)

backprop :: LabelsMat -> [Weight] -> [Bias] -> [Activation] -> ([Weight], [Bias])
backprop y ws bs (aL:as) = (zipWithSafe (-) ws (zipWithSafe (<>) deltaList (tr <$> as)), zipWithSafe (-) bs $ meanRows <$> deltaList)
    where
        deltaList :: [Matrix Double]
        deltaList = init $ deltas [edgeDelta y aL] y as ws

deltas :: [Matrix Double]
       -> LabelsMat
       -> [Activation]
       -> [Weight]
       -> [Matrix Double]
deltas deltal _ [] [] = deltal
deltas deltal y (a:as) (wlp1:ws) = deltas (deltal++[(tr wlp1 <> last deltal) * sigmoid' a]) y as ws

edgeDelta :: LabelsMat -- ^ Correct output Y for input X
           -> Activation -- ^ Predicted ouput for input X
           -> Matrix Double
edgeDelta y aL = (aL - y) * sigmoid' aL

{-- Arithmetic Functions --}
cost :: Matrix Double -> Matrix Double -> Matrix Double
cost a y = ((a-y)**2) / 2

sigmoid :: Floating a => a -> a
sigmoid x = 1 / (1 + (exp (-x)))

sigmoid' :: Floating a => a -> a
sigmoid' x = x * (1 - x)

normalize :: (Integral a, Floating b) => a -> b
normalize x = (fromIntegral x) / 255

gauss :: IO Double
gauss = do
    u1 <- R.randomRIO (0 :: Double, 1 :: Double)
    u2 <- R.randomRIO (0 :: Double, 1 :: Double)
    return $ (sqrt (-2*log u1)) * (cos (2*pi*u2))

{-- Data Processing Functions --}
initWB :: NeuralNet -> IO NeuralNet
initWB net@NN{batchSize = bs, layerSize = ls, layers = n} = do
    headW <- randn ls 784
    midW  <- replicateM (n-2) (randn ls ls)
    lastW <- randn 10 ls

    initB <- replicateM (n-1) (randn ls 1)
    lastB <- randn 10 1

    let weights = reverse $ headW : midW ++ [lastW]
    let biases = reverse . fmap flatten $ initB ++ [lastB]

    return net { weights = weights
               , biases  = biases }

randomMatrix :: Int -> Int -> IO (Matrix Double)
randomMatrix nrows ncols = do
    gaussList <- replicateM (ncols*nrows) gauss
    return $ (nrows><ncols) gaussList

loadData :: IO [(Vector Double, Vector Double)]
loadData = do
    trainImgs <- getData "train-images-idx3-ubyte.gz"
    trainLabels <- getData "train-labels-idx1-ubyte.gz"
    let !labels = BL.toStrict trainLabels
    let !imgs = BL.toStrict trainImgs
    let !l = [vectorizeLabel $ getLabel n labels | n <- [0..49999]] :: [Label]
    let !i = [getImage n imgs | n <- [0..49999]] :: [Image]
    return $ zip i l

loadTestData :: IO [(Matrix Double, Int)]
loadTestData = do
    testImgs <- getData "t10k-images-idx3-ubyte.gz"
    testLabels <- getData "t10k-labels-idx1-ubyte.gz"
    let !labels = BL.toStrict testLabels
    let !imgs = BL.toStrict testImgs
    let !l = [getLabel n labels | n <- [0..9999]] :: [Int]
    let !i = asColumn <$> [getImage n imgs | n <- [0..9999]] :: [Matrix Double]
    return $ zip i l

getData :: FilePath -> IO BL.ByteString
getData path = do
    currentDir <- getCurrentDirectory
    fileData <- decompress <$> BL.readFile (currentDir ++ "/data/mnist_dataset/" ++ path)
    return fileData

getImage :: Int -> BS.ByteString -> Image
getImage n imgs = fromList [normalize $ BS.index imgs (16 + n*784 + s) | s <- [0..783]]

-- | `getImage` but lazy
getImage' :: Int64 -> BL.ByteString -> Image
getImage' n imgs = fromList [normalize $ BL.index imgs (16 + n*784 + s) | s <- [0..783]]

getLabel :: Num a => Int -> BS.ByteString -> a
getLabel n labels = fromIntegral $ BS.index labels (n+8)

getGuess :: [Activation] -> Int
getGuess = maxIndex . flatten . head

vectorizeLabel :: Int -> Vector Double
vectorizeLabel l = fromList $ x ++ 1 : y
    where (x,y) = splitAt l $ replicate 9 0

{-- HELPER FUNCTIONS --}
toMatrix :: Vector Double -> Matrix Double
toMatrix = reshape 1

-- | Means of rows in matrix as a column vector
meanRows :: Matrix Double -> Vector Double
meanRows = fst . meanCov . tr

-- | zipWith, but throws error if lengths don't match
zipWithSafe :: (a -> b -> c) -> [a] -> [b] -> [c]
zipWithSafe f as bs
    | length as == length bs = zipWith f as bs
    | otherwise = error $ "zipWith: lengths don't match (" ++ (show $ length as) ++ " " ++ (show $ length bs) ++ ")"

-- | Randomly shuffle a list
--   /O(N)/
shuffle :: [a] -> IO [a]
shuffle xs = do
        ar <- newArray n xs
        forM [1..n] $ \i -> do
            j <- R.randomRIO (i,n)
            vi <- A.readArray ar i
            vj <- A.readArray ar j
            A.writeArray ar j vi
            return vj
  where
    n = length xs
    newArray :: Int -> [a] -> IO (A.IOArray Int a)
    newArray n xs =  A.newListArray (1,n) xs

max_stars_repo_name: WhateverHappns/haskell-mnist

max_issues_repo_name: WhateverHappns/haskell-mnist

max_forks_repo_name: WhateverHappns/haskell-mnist


text:
-- (C) Copyright Chris Banks 2011-2012

-- This file is part of The Continuous Pi-calculus Workbench (CPiWB).

--     CPiWB is free software: you can redistribute it and/or modify
--     it under the terms of the GNU General Public License as published by
--     the Free Software Foundation, either version 3 of the License, or
--     (at your option) any later version.

--     CPiWB is distributed in the hope that it will be useful,
--     but WITHOUT ANY WARRANTY; without even the implied warranty of
--     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
--     GNU General Public License for more details.

--     You should have received a copy of the GNU General Public License
--     along with CPiWB.  If not, see <http://www.gnu.org/licenses/>.

{-# OPTIONS_GHC -w #-}

module CPi.Matlab
    (solveODEoctave,
     matlabScript
    ) where

import BioCalcLib.Lib
import CPi.Lib
import CPi.Semantics
import CPi.ODE

import qualified Data.List as L
import qualified Control.Exception as X
import qualified Data.Map as Map
import Data.Map (Map)
--import Data.Maybe

--import qualified Numeric.LinearAlgebra as LA

import System.IO.Unsafe
import qualified System.Process as OS

-------------------------------
-- MATLAB script output:
-------------------------------

matlabODE :: Env -> Process -> P' -> (Int,(Double,Double)) -> String
matlabODE env p dpdt (n,(t0,tn))
    = "init = [" ++ (L.concat (L.intersperse ";" (map show (initials env p p')))) ++ "];\n" ++
      "t = linspace(" ++ show t0 ++ "," ++ show tn ++ "," ++ show n ++ ");\n" ++
      "function xdot = f(x,t) \n\n" ++
      (L.concat (L.intersperse ";\n"
                 ((map (\(k,v) ->
                        (maybe "XXX" id (Map.lookup k (svxd)) ++ " = " ++
                         matlabExpr env (svx) v)))
                  (Map.toList p')))) ++ ";\nendfunction;\n\n" ++
      "function jac = jj(x,t)\n\n" ++
      matlabJac env p' ++ ";\nendfunction;\n\n" ++
      "x = lsode({@f,@jj}, init, t);\n" ++
      "save (\"-ascii\", \"-\", \"x\");"
          where
            xdots = ["xdot("++ show i ++")" | i<-[1..]]
            svxd = speciesVars env xdots p'
            xs = ["x("++ show i ++")" | i<-[1..]]
            svx = speciesVars env xs p'
            p' = simpP' env dpdt

matlabExpr :: Env -> Map Species String -> Expr -> String
matlabExpr env vs (Var x)
    = maybe (ex x) id (Map.lookup x vs)
      where
        ex x = X.throw $ CpiException
               ("Bug: failed var lookup in CPi.ODE.matlabExpr: " ++ show x)
matlabExpr env vs (Plus x y) = matlabExpr env vs x ++ " .+ " ++ matlabExpr env vs y
matlabExpr env vs (Times (Plus x y) (Plus x' y'))
    = "(" ++ matlabExpr env vs (Plus x y) ++ ").*(" ++ matlabExpr env vs (Plus x' y') ++ ")"
matlabExpr env vs (Times (Plus x y) z)
    = "(" ++ matlabExpr env vs (Plus x y) ++ ").*" ++ matlabExpr env vs z
matlabExpr env vs (Times x (Plus y z))
    = matlabExpr env vs x ++ ".*(" ++ matlabExpr env vs (Plus y z) ++ ")"
matlabExpr env vs (Times x y) = matlabExpr env vs x ++ ".*" ++ matlabExpr env vs y
matlabExpr env vs (Num k) = show k

matlabJac :: Env -> P' -> String
matlabJac env p' = L.concat $ L.intersperse ";\n" $
                   map (\(a,b) -> a ++ " = " ++ b) $
                   zip rhss $
                   map (matlabExpr env (speciesVars env xs p')) $
                   map (simp env) $
                   map diff' $ cp (unzip (Map.toList (sp')))
    where
      cp (x,y) = [(a,b)|b<-y,a<-x]
      diff' (x,e) = diff x e
      xs = ["x("++ show i ++")" | i<-[1..]]
      sp' = simpP' env p'
      len = Map.size p'
      rhss = ["jac("++ show x ++ "," ++ show y ++ ")" | x<-[1..len], y<-[1..len]]

---------------------------------------
-- Using Octave to execute the scripts
---------------------------------------

callOctave env p mts p' ts = let
    script = matlabODE env (wholeProc env p mts) p' ts
  in do
    putStrLn $ "Octave Script: \n" ++ script
    OS.readProcess
      "octave" ["-q", "--eval", script] []

-- | Solver which calculates the symbolic Jacobian, writes MATLAB code, and executes it with GNU Octave. (General purpose, deals with stiff systems, uses LSODE.)
solveODEoctave :: Solver
solveODEoctave env p mts p' ts@(n,(t0,tn))
    = let raw = unsafePerformIO (callOctave env p mts p' ts)
      in (n><(Map.size p')) $ map s2d $ words raw

-- Return the MATLAB script for ODEs
matlabScript :: Env
             -> Process
             -> MTS
             -> P'
             -> (Int, (Double, Double))
             -> String
matlabScript env p mts p' ts = matlabODE env (wholeProc env p mts) p' ts

max_stars_repo_name: twright/bondwb

max_issues_repo_name: twright/bondwb

max_forks_repo_name: twright/bondwb


text:
{-# LANGUAGE ForeignFunctionInterface #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Foreign.BLAS.Zomplex
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--

module Foreign.BLAS.Zomplex
    where

import Data.Complex( Complex )
import Foreign
import Foreign.BLAS.Types
import Foreign.C.Types

#include "config.h"
#include "f77_func-hsc.h"
#define la_int int

---------------------------- Level 1 Routines -------------------------------

foreign import ccall unsafe #f77_func zdotu
    zdotu :: Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zdotc
    zdotc :: Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func dznrm2
    znrm2  :: Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO Double

foreign import ccall unsafe #f77_func dzasum
    zasum  :: Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO Double

foreign import ccall unsafe #f77_func izamax
    izamax_hidden :: Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO CInt
izamax :: Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO LAInt
izamax a b c = do
                res <- izamax_hidden a b c
                return $! LAInt res

foreign import ccall unsafe #f77_func zscal
    zscal  :: Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zswap
    zswap  :: Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zcopy
    zcopy  :: Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zaxpy
    zaxpy  :: Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zrotg
    zrotg  :: Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe #f77_func zdrot
    zdrot :: Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr Double -> Ptr Double -> IO ()

---------------------------- Level 2 Routines -------------------------------

foreign import ccall unsafe #f77_func zgemv
    zgemv :: BLASTrans -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zgbmv
    zgbmv ::  BLASTrans -> Ptr LAInt -> Ptr LAInt -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func ztrmv
    ztrmv ::  BLASUplo -> BLASTrans -> BLASDiag -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func ztpmv
    ztpmv ::  BLASUplo -> BLASTrans -> BLASDiag -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func ztpsv
    ztpsv ::  BLASUplo -> BLASTrans -> BLASDiag -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func ztbmv
    ztbmv ::  BLASUplo -> BLASTrans -> BLASDiag -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func ztrsv
    ztrsv ::  BLASUplo -> BLASTrans -> BLASDiag -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func ztbsv
    ztbsv ::  BLASUplo -> BLASTrans -> BLASDiag -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zhemv
    zhemv ::  BLASUplo -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zhbmv
    zhbmv ::  BLASUplo -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zgeru
    zgeru  ::  Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zgerc
    zgerc  ::  Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zher
    zher  ::  BLASUplo -> Ptr LAInt -> Ptr Double -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zher2
    zher2 ::  BLASUplo -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zhpmv
    zhpmv ::  BLASUplo -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zhpr
    zhpr  ::  BLASUplo -> Ptr LAInt -> Ptr Double -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe #f77_func zhpr2
    zhpr2 ::  BLASUplo -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> IO ()

---------------------------- Level 3 Routines -------------------------------

foreign import ccall unsafe #f77_func zgemm
    zgemm  ::  BLASTrans -> BLASTrans -> Ptr LAInt -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zsymm
    zsymm  ::  BLASSide -> BLASUplo -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zhemm
    zhemm  ::  BLASSide -> BLASUplo -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func ztrmm
    ztrmm  ::  BLASSide -> BLASUplo -> BLASTrans -> BLASDiag -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func ztrsm
    ztrsm  ::  BLASSide -> BLASUplo -> BLASTrans -> BLASDiag -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zsyrk
    zsyrk  ::  BLASUplo -> BLASTrans -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zsyr2k
    zsyr2k ::  BLASUplo -> BLASTrans -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zherk
    zherk  ::  BLASUplo -> BLASTrans -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

foreign import ccall unsafe #f77_func zher2k
    zher2k ::  BLASUplo -> BLASTrans -> Ptr LAInt -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr LAInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr LAInt -> IO ()

ext: hsc

max_stars_repo_name: cartazio/hs-cblas

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE OverloadedStrings, ScopedTypeVariables #-}

module Main where

import Control.Monad (forM, forM_)
import Data.List (intercalate)
import Statistics.Sample (mean, stdDev)
import Turtle
import qualified Data.Text as T
import qualified Data.Vector as V

data Platform = PlatformLinux | PlatformMacOS deriving (Eq)

main :: IO ()
main = do
    uname <- strict $ inshell "uname" empty
    platform <-
        if "Linux" `T.isInfixOf` uname then return PlatformLinux
        else if "Darwin" `T.isInfixOf` uname then return PlatformMacOS
        else failWithMsg "Currently only Linux and macOS are supported."

    -- Required for pcregrep on Linux.
    when (platform == PlatformLinux) $ export "LC_ALL" "en_US.UTF-8"

    let logFile :: (IsString a) => a
        logFile = "./tmp/bench.sh.log"

    export "logfile" (fromString logFile)
    testfile logFile >>= flip when (rm logFile)

    echo "Compiling C program... "
    shells "mkdir -p ./tmp" empty
    shells ("gcc -O2 -Wall bench-archive.c -I/usr/local/opt/libarchive/include -L/usr/local/opt/libarchive/lib "
            `T.append` "-larchive -o bench-archive >> $logfile 2>&1") empty

    echo "Compiling Haskell programs... "
    shells "stack build >> $logfile 2>&1" empty

    let c_executable = "./bench-archive"
    hs_plain_executable <- T.strip <$> strict (inshell "stack exec -- which bench-archive-plain 2> /dev/null" empty)
    hs_streamly_executable <- T.strip <$> strict (inshell "stack exec -- which bench-archive-streamly 2> /dev/null" empty)

    let fileCountsAndFileSizes :: [(Int,Int)] =
            [(2500000,1), (5000000,1), (20,250000000), (40,250000000), (80,250000000)]

    echoHrule
    answerIsYes "Create archives?" >>= flip when (do
        echo "Creating archives using C program (unless they already exist)..."
        let tmpDir :: (IsString a) => a
            tmpDir = "./tmp"
        mktree tmpDir

        forM_ fileCountsAndFileSizes $ \(fileCount, fileSize') -> do
            let archivePath = format (""%fp%"/test_"%d%"_"%d%".tar") tmpDir fileCount fileSize'
            (_, existingFileCount) <-
                shellStrict (format ("tar -tf "%s%" 2> /dev/null | wc -l | bc") archivePath) empty
            if T.strip existingFileCount == format (""%d%"") fileCount then
                procs "echo" [format ("    Archive "%s%" already exists.") archivePath] empty
            else do
                procs "echo" ["-n", format ("    Creating archive "%s%"... ") archivePath] empty
                shells (format ("rm -rf "%s%" 2> /dev/null") archivePath) empty
                shells (format (""%s%" write "%s%" "%d%" "%d%"") c_executable archivePath fileCount fileSize') empty
                echo "success.")

    echoHrule
    answerIsYes "Measure read performance?" >>= flip when (do
        let tmpDir :: (IsString a) => a
            tmpDir = "./tmp"
        mktree tmpDir

        let csvFile :: (IsString a) => a
            csvFile = "./tmp/read.csv"
        shells (format ("rm -f "%fp%"") csvFile) empty

        output csvFile . return . unsafeTextToLine . T.pack $
            "fileCount,fileSize,"
            ++ "c_mean,c_std,"
            ++ "hsPlain_mean,hsPlain_std,"
            ++ "hsStreamly_mean,hsStreamly_std,"
            ++ "hsPlain_div_c_mean,hsPlain_div_c_std,"
            ++ "hsStreamly_div_c_mean,hsStreamly_div_c_std,"
            ++ "hsStreamly_div_hsPlain_mean,hsStreamly_div_hsPlain_std"

        forM_ fileCountsAndFileSizes $ \(fileCount, fileSize') -> do
            let archivePath = format (""%fp%"/test_"%d%"_"%d%".tar") tmpDir fileCount fileSize'
            let args  = format ("read "%s%"") archivePath
            let expectedInOutput1 = format ("Total filesize: "%d%"") (fileSize' * fileCount)
            let expectedInOutput2 = format ("File count:     "%d%"") fileCount
            let warmCount = 1
            let timeCount = 10
            procs "echo" [format ("    Reading archive "%s%" with C...") archivePath] empty
            cTimes <- timeCommand Nothing (format (""%s%" "%s%"") c_executable args)
                        [expectedInOutput1, expectedInOutput2] warmCount timeCount
            procs "echo" [format ("    Reading archive "%s%" with Haskell (plain)...") archivePath] empty
            hsPlainTimes <- timeCommand Nothing (format (""%s%" "%s%"") hs_plain_executable args)
                        [expectedInOutput1, expectedInOutput2] warmCount timeCount
            procs "echo" [format ("    Reading archive "%s%" with Haskell (streamly)...") archivePath] empty
            hsStreamlyTimes <- timeCommand Nothing (format (""%s%" "%s%"") hs_streamly_executable args)
                        [expectedInOutput1, expectedInOutput2] warmCount timeCount

            let hsPlainDivCTimes = [a / b | a <- hsPlainTimes, b <- cTimes]
            let hsStreamlyDivCTimes = [a / b | a <- hsStreamlyTimes, b <- cTimes]
            let hsStreamlyDivHsPlainTimes = [a / b | a <- hsStreamlyTimes, b <- hsPlainTimes]

            let mean' :: [NominalDiffTime] -> Double
                mean' xs = mean . V.fromList $ map realToFrac xs

            let std' :: [NominalDiffTime] -> Double
                std' xs = stdDev . V.fromList $ map realToFrac xs

            let line = intercalate "," $ [show fileCount, show fileSize'] ++ map show [
                    mean' cTimes, std' cTimes,
                    mean' hsPlainTimes, std' hsPlainTimes,
                    mean' hsStreamlyTimes, std' hsStreamlyTimes,
                    mean' hsPlainDivCTimes, std' hsPlainDivCTimes,
                    mean' hsStreamlyDivCTimes, std' hsStreamlyDivCTimes,
                    mean' hsStreamlyDivHsPlainTimes, std' hsStreamlyDivHsPlainTimes ]
            append csvFile . return . unsafeTextToLine . T.pack $ line
            procs "echo" [format ("    Appended results to "%s%"") csvFile] empty
            return ())

failWithMsg :: Line -> IO a
failWithMsg msg = echo msg >> exit (ExitFailure 1)

echoHrule :: IO ()
echoHrule = echo "----------------------------------------"

answerIsYes :: String -> IO Bool
answerIsYes question = do
    procs "echo" ["-n", T.pack $ question ++ " (y/N) "] empty
    line <- readline
    case line of
        Just answer | answer == "y" -> return True
        _ -> return False

timeCommand' :: Text -> [Text] -> IO NominalDiffTime
timeCommand' cmdAndArgs expectedInOut = do
    (out, t) <- time . strict $ inshell cmdAndArgs empty
    if all (`T.isInfixOf` out) expectedInOut then
        return t
    else
        error $ "all of " ++ show expectedInOut ++ " were not found in output "
            ++ show out ++ " of command " ++ show cmdAndArgs

timeCommand :: Maybe T.Text -> Text -> [Text] -> Int -> Int -> IO [NominalDiffTime]
timeCommand forceRemovalPath cmdAndArgs expectedInOut warmCount timeCount = do
    let removePath = maybe (return ()) (\fp' -> shells (format ("rm -rf "%s%"") fp') empty) forceRemovalPath
    let echoCount c = procs "echo" ["-n", format (""%d%" ") c] empty
    when (warmCount > 0) $ do
        procs "echo" ["-n", "        Warming up... "] empty
        forM_ [1..warmCount] $ \c ->
            echoCount c >> removePath >> timeCommand' cmdAndArgs expectedInOut
        echo "done."
    procs "echo" ["-n", "        Timing... "] empty
    times <- forM [1..timeCount] $ \c ->
        echoCount c >> removePath >> timeCommand' cmdAndArgs expectedInOut
    removePath
    echo "done."
    return times

max_stars_repo_name: shlok/streamly-archive

max_issues_repo_name: shlok/streamly-archive

max_forks_repo_name: shlok/streamly-archive


text:
{-# LANGUAGE UndecidableInstances #-} --For the BitPack instance
module Clash.DSP.Complex where

import Clash.Prelude

import GHC.Generics
import Test.QuickCheck

import qualified Data.Complex as C

{-| I defined my own complex type so that I can write a Num instance without the RealFloat constraint. TODO: think about whether this is really a good idea. -}
data Complex a = a :+ a deriving (Show, Lift, Generic, ShowX, NFDataX, Functor, Foldable, Traversable, Eq, Ord)

deriving instance (BitPack a, KnownNat (BitSize a)) => BitPack (Complex a)

instance Num a => Num (Complex a) where
    (a :+ b) + (c :+ d) = (a + c) :+ (b + d)
    (a :+ b) - (c :+ d) = (a - c) :+ (b - d)
    (a :+ b) * (c :+ d) = (a * c - b * d) :+ (a * d + b * c)
    fromInteger x       = (fromInteger x :+ 0)

instance Applicative Complex where
    pure a = a :+ a
    f :+ g <*> a :+ b = f a :+ g b

instance Arbitrary a => Arbitrary (Complex a)
    where
    arbitrary = liftA2 (:+) arbitrary arbitrary

instance Default a => Default (Complex a)
    where
    def = def :+ def

realPart :: Complex a -> a
realPart (x :+ _) = x

imagPart :: Complex a -> a
imagPart (_ :+ y) = y

fromComplex (a C.:+ b) = a :+ b
toComplex   (a :+ b)   = a C.:+ b

conjugate :: Num a => Complex a -> Complex a
conjugate (x :+ y) = x :+ (-y)

cMul
    :: ExtendingNum a b
    => ExtendingNum (MResult a b) (MResult a b)
    => Complex a
    -> Complex b
    -> Complex (AResult (MResult a b) (MResult a b))
cMul (xr :+ xi) (yr :+ yi)
    =  ((xr `mul` yr) `sub` (xi `mul` yi))
    :+ ((xr `mul` yi) `add` (xi `mul` yr))

cMul3
    :: ExtendingNum b b
    => ExtendingNum a a
    => ExtendingNum a (AResult b b)
    => ExtendingNum b (AResult a a)
    => ExtendingNum (MResult b (AResult a a)) (MResult a (AResult b b))
    => Complex a
    -> Complex b
    -> Complex (AResult (MResult b (AResult a a)) (MResult a (AResult b b)))
cMul3 (xr :+ xi) (yr :+ yi)
    =  (k1 `sub` k3)
    :+ (k1 `add` k2)
    where
    k1 = yr `mul` (xr `add` xi)
    k2 = xr `mul` (yi `sub` yr)
    k3 = xi `mul` (yr `add` yi)

cMulPipe
    :: HiddenClockResetEnable dom
    => ExtendingNum a b
    => ExtendingNum (MResult a b) (MResult a b)
    => NFDataX (MResult a b)
    => NFDataX (AResult (MResult a b) (MResult a b))
    => Signal dom Bool
    -> Signal dom (Complex a)
    -> Signal dom (Complex b)
    -> Signal dom (Complex (AResult (MResult a b) (MResult a b)))
cMulPipe en x y
    = liftA2 (:+) rr ri
    where
    --Products
    xryr = regEn (errorX "initial xryr") en $ liftA2 mul (realPart <$> x) (realPart <$> y)
    xiyi = regEn (errorX "initial xiyi") en $ liftA2 mul (imagPart <$> x) (imagPart <$> y)
    xryi = regEn (errorX "initial xryi") en $ liftA2 mul (realPart <$> x) (imagPart <$> y)
    xiyr = regEn (errorX "initial xiyr") en $ liftA2 mul (imagPart <$> x) (realPart <$> y)
    --Sums
    rr   = regEn (errorX "initial rr")   en $ liftA2 sub xryr xiyi
    ri   = regEn (errorX "initial ri")   en $ liftA2 add xryi xiyr

cMul3Pipe
    :: HiddenClockResetEnable dom
    => ExtendingNum b b
    => ExtendingNum a a
    => ExtendingNum a (AResult b b)
    => ExtendingNum b (AResult a a)
    => ExtendingNum (MResult b (AResult a a)) (MResult a (AResult b b))
    => NFDataX (AResult a a)
    => NFDataX (AResult b b)
    => NFDataX (MResult b (AResult a a))
    => NFDataX (MResult a (AResult b b))
    => NFDataX (AResult (MResult b (AResult a a)) (MResult a (AResult b b)))
    => NFDataX a
    => NFDataX b
    => Signal dom Bool
    -> Signal dom (Complex a)
    -> Signal dom (Complex b)
    -> Signal dom (Complex (AResult (MResult b (AResult a a)) (MResult a (AResult b b))))
cMul3Pipe en x y
    =  liftA2 (:+) rr ri
    where
    --Sums
    xSum    = regEn (errorX "initial xSum")    en $ liftA2 add (realPart <$> x) (imagPart <$> x)
    yDiff   = regEn (errorX "initial yDiff")   en $ liftA2 sub (imagPart <$> y) (realPart <$> y)
    ySum    = regEn (errorX "initial ySum")    en $ liftA2 add (realPart <$> y) (imagPart <$> y)
    --Delays
    xD      = regEn (errorX "initial xD") en x
    yD      = regEn (errorX "initial yD") en y
    --Products
    yrXSum  = regEn (errorX "initial yrXSum")  en $ liftA2 mul (realPart <$> yD) xSum
    xrYDiff = regEn (errorX "initial xrYDiff") en $ liftA2 mul (realPart <$> xD) yDiff
    xiYSum  = regEn (errorX "initial xiYSum")  en $ liftA2 mul (imagPart <$> xD) ySum
    --Sums
    rr      = regEn (errorX "initial rr")      en $ liftA2 sub yrXSum xiYSum
    ri      = regEn (errorX "initial ri")      en $ liftA2 add yrXSum xrYDiff

max_stars_repo_name: adamwalker/clash-utils

max_issues_repo_name: adamwalker/clash-utils

max_forks_repo_name: adamwalker/clash-utils