text:
{-# LANGUAGE BangPatterns, RecordWildCards #-}
-- |
-- Module      : Criterion
-- Copyright   : (c) 2009-2014 Bryan O'Sullivan
--
-- License     : BSD-style
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : GHC
--
-- Core benchmarking code.

module Criterion.Internal
    (
      runAndAnalyse
    , runAndAnalyseOne
    , runOne
    , runFixedIters
    ) where

import qualified Data.Aeson as Aeson
import Control.DeepSeq (rnf)
import Control.Exception (evaluate)
import Control.Monad (foldM, forM_, void, when, unless)
import Control.Monad.Catch (MonadMask, finally)
import Control.Monad.Reader (ask, asks)
import Control.Monad.Trans (MonadIO, liftIO)
import Control.Monad.Trans.Except
import qualified Data.Binary as Binary
import Data.Int (Int64)
import qualified Data.ByteString.Lazy.Char8 as L
import Criterion.Analysis (analyseSample, noteOutliers, analyseSimpleSample)
import Criterion.IO (header, headerRoot, critVersion, readJSONReports, writeJSONReports)
import Criterion.IO.Printf (note, printError, prolix, writeCsv)
import Criterion.Measurement (runBenchmark, runBenchmarkable_, secs)
import Criterion.Measurement.Types (hasFixedIters)
import Criterion.Monad (Criterion)
import Criterion.Report (report)
import Criterion.Types hiding (measure)
import Criterion.Measurement.Types.Internal (fakeEnvironment)
import qualified Data.Map as Map
import qualified Data.Vector as V
import Statistics.Types (Estimate(..),ConfInt(..),confidenceInterval,cl95,confidenceLevel)
import System.Directory (getTemporaryDirectory, removeFile)
import System.IO (IOMode(..), hClose, openTempFile, openFile, hPutStr, openBinaryFile)
import Text.Printf (printf)

-- | Run a single benchmark.
runOne :: Int -> String -> Benchmarkable -> Criterion DataRecord
runOne i desc bm = do
  Config{..} <- ask
  (meas,timeTaken) <- liftIO $ runBenchmark bm timeLimit
  when (timeTaken > timeLimit * 1.25) .
    void $ prolix "measurement took %s\n" (secs timeTaken)
  return (Measurement i desc meas)

-- | Analyse a single benchmark.
analyseOne :: Int -> Bool -> String -> V.Vector Measured -> Criterion DataRecord
analyseOne i simple desc meas = do
  Config{..} <- ask
  _ <- prolix "analysing with %d resamples\n" resamples
  erp <- runExceptT $ if simple
    then analyseSimpleSample i desc meas
    else analyseSample i desc meas
  case erp of
    Left err -> printError "*** Error: %s\n" err
    Right rpt@Report{..} -> do
      let SampleAnalysis{..} = reportAnalysis
          OutlierVariance{..} = anOutlierVar
          wibble = case ovEffect of
                     Unaffected -> "unaffected" :: String
                     Slight -> "slightly inflated"
                     Moderate -> "moderately inflated"
                     Severe -> "severely inflated"
          (builtin, others) = splitAt 1 anRegress
      let r2 n = printf "%.3f R\178" n
      forM_ builtin $ \Regression{..} ->
        case Map.lookup "iters" regCoeffs of
          Nothing -> return ()
          Just t  -> bs secs "time" t >> bs r2 "" regRSquare
      bs secs "mean" anMean
      bs secs "std dev" anStdDev
      forM_ others $ \Regression{..} -> do
        _ <- bs r2 (regResponder ++ ":") regRSquare
        forM_ (Map.toList regCoeffs) $ \(prd,val) ->
          bs (printf "%.3g") ("  " ++ prd) val
      writeCsv
        (desc,
         estPoint anMean,   fst $ confidenceInterval anMean,   snd $ confidenceInterval anMean,
         estPoint anStdDev, fst $ confidenceInterval anStdDev, snd $ confidenceInterval anStdDev
        )
      when (verbosity == Verbose || (ovEffect > Slight && verbosity > Quiet)) $ do
        when (verbosity == Verbose) $ noteOutliers reportOutliers
        _ <- note "variance introduced by outliers: %d%% (%s)\n"
             (round (ovFraction * 100) :: Int) wibble
        return ()
      _ <- note "\n"
      return (Analysed rpt)
      where bs :: (Double -> String) -> String -> Estimate ConfInt Double -> Criterion ()
            bs f metric e@Estimate{..} =
              note "%-20s %-10s (%s .. %s%s)\n" metric
                   (f estPoint) (f $ fst $ confidenceInterval e) (f $ snd $ confidenceInterval e)
                   (let cl = confIntCL estError
                        str | cl == cl95 = ""
                            | otherwise  = printf ", ci %.3f" (confidenceLevel cl)
                    in str
                   )

{-}
analyseOneSimple :: Int -> String -> V.Vector Measured
                 -> Criterion ()
analyseOneSimple i desc meas = do
    erp <- runExceptT $ analyseSimpleSample i desc meas
    case erp of
      Left err -> printError "*** Error: %s\n" err
      Right (stMean, stDev) -> do
        bs secs "mean" stMean
        bs secs "std dev" stDev
  where
    bs :: (Double -> String) -> String -> Estimate ConfInt Double -> Criterion ()
    bs f metric e@Estimate{..} =
      note "%-20s %-10s (%s .. %s%s)\n" metric
           (f estPoint) (f $ fst $ confidenceInterval e) (f $ snd $ confidenceInterval e)
           (let cl = confIntCL estError
                str | cl == cl95 = ""
                    | otherwise  = printf ", ci %.3f" (confidenceLevel cl)
            in str
           )
-}

-- | Run a single benchmark and analyse its performance.
runAndAnalyseOne :: Int -> String -> Benchmarkable -> Criterion DataRecord
runAndAnalyseOne i desc bm = do
  Measurement _ _ meas <- runOne i desc bm
  analyseOne i (hasFixedIters bm) desc meas

-- | Run, and analyse, one or more benchmarks.
runAndAnalyse :: (String -> Bool) -- ^ A predicate that chooses
                                  -- whether to run a benchmark by its
                                  -- name.
              -> Benchmark
              -> Criterion ()
runAndAnalyse select bs = do
  mbJsonFile <- asks jsonFile
  (jsonFile, handle) <- liftIO $
    case mbJsonFile of
      Nothing -> do
        tmpDir <- getTemporaryDirectory
        openTempFile tmpDir "criterion.json"
      Just file -> do
        handle <- openFile file WriteMode
        return (file, handle)
  -- The type we write to the file is ReportFileContents, a triple.
  -- But here we ASSUME that the tuple will become a JSON array.
  -- This assumption lets us stream the reports to the file incrementally:
  liftIO $ hPutStr handle $ "[ \"" ++ headerRoot ++ "\", " ++
                             "\"" ++ critVersion ++ "\", [ "

  for select bs $ \idx desc bm -> do
    _ <- note "benchmarking %s\n" desc
    Analysed rpt <- runAndAnalyseOne idx desc bm
    unless (idx == 0) $
      liftIO $ hPutStr handle ", "
    liftIO $ L.hPut handle (Aeson.encode (rpt::Report))

{-}
  for select bs $ \idx desc bm -> do
    _ <- note "benchmarking %s\n" desc
    if (hasFixedIters bm)
    then do
      measur <- runAndAnalyseOne idx desc bm
      liftIO $ print measur
    else do
      Analysed rpt <- runAndAnalyseOne idx desc bm
      unless (idx == 0) $
        liftIO $ hPutStr handle ", "
      liftIO $ L.hPut handle (Aeson.encode (rpt::Report))
-}
  liftIO $ hPutStr handle " ] ]\n"
  liftIO $ hClose handle

  rpts <- liftIO $ do
    res <- readJSONReports jsonFile
    case res of
      Left err -> error $ "error reading file "++jsonFile++":\n  "++show err
      Right (_,_,rs) ->
       case mbJsonFile of
         Just _ -> return rs
         _      -> removeFile jsonFile >> return rs

  rawReport rpts
  report rpts
  json rpts
  junit rpts

-- | Write out raw binary report files.  This has some bugs, including and not
-- limited to #68, and may be slated for deprecation.
rawReport :: [Report] -> Criterion ()
rawReport reports = do
  mbRawFile <- asks rawDataFile
  case mbRawFile of
    Nothing   -> return ()
    Just file -> liftIO $ do
      handle <- openBinaryFile file ReadWriteMode
      L.hPut handle header
      forM_ reports $ \rpt ->
        L.hPut handle (Binary.encode rpt)
      hClose handle

-- | Run a benchmark without analysing its performance.
runFixedIters :: Int64            -- ^ Number of loop iterations to run.
              -> (String -> Bool) -- ^ A predicate that chooses
                                  -- whether to run a benchmark by its
                                  -- name.
              -> Benchmark
              -> Criterion ()
runFixedIters iters select bs =
  for select bs $ \_idx desc bm -> do
    _ <- note "benchmarking %s\n" desc
    liftIO $ runBenchmarkable_ bm iters

-- | Iterate over benchmarks.
for :: (MonadMask m, MonadIO m) => (String -> Bool) -> Benchmark
    -> (Int -> String -> Benchmarkable -> m ()) -> m ()
for select bs0 handle = go (0::Int) ("", bs0) >> return ()
  where
    go !idx (pfx, Environment mkenv cleanenv mkbench)
      | shouldRun pfx mkbench = do
        e <- liftIO $ do
          ee <- mkenv
          evaluate (rnf ee)
          return ee
        go idx (pfx, mkbench e) `finally` liftIO (cleanenv e)
      | otherwise = return idx
    go idx (pfx, Benchmark desc b)
      | select desc' = do handle idx desc' b; return $! idx + 1
      | otherwise    = return idx
      where desc' = addPrefix pfx desc
    go idx (pfx, BenchGroup desc bs) =
      foldM go idx [(addPrefix pfx desc, b) | b <- bs]

    shouldRun pfx mkbench =
      any (select . addPrefix pfx) . benchNames . mkbench $ fakeEnvironment

-- | Write summary JSON file (if applicable)
json :: [Report] -> Criterion ()
json rs
  = do jsonOpt <- asks jsonFile
       case jsonOpt of
         Just fn -> liftIO $ writeJSONReports fn rs
         Nothing -> return ()

-- | Write summary JUnit file (if applicable)
junit :: [Report] -> Criterion ()
junit rs
  = do junitOpt <- asks junitFile
       case junitOpt of
         Just fn -> liftIO $ writeFile fn msg
         Nothing -> return ()
  where
    msg = "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n" ++
          printf "<testsuite name=\"Criterion benchmarks\" tests=\"%d\">\n"
          (length rs) ++
          concatMap single rs ++
          "</testsuite>\n"
    single Report{..} = printf "  <testcase name=\"%s\" time=\"%f\" />\n"
               (attrEsc reportName) (estPoint $ anMean $ reportAnalysis)
    attrEsc = concatMap esc
      where
        esc '\'' = "'"
        esc '"'  = """
        esc '<'  = "<"
        esc '>'  = ">"
        esc '&'  = "&"
        esc c    = [c]


text:
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Test.Grenade.Layers.FullyConnected where

import           System.Random.MWC             (create)

import           Data.Constraint               (Dict (..))
import           Data.Proxy
import           Data.Singletons               ()
import           Data.Kind                     (Type)
import           Data.Serialize
import           Data.Either

import           GHC.TypeLits
import           Unsafe.Coerce                 (unsafeCoerce)

import           Hedgehog
import qualified Hedgehog.Range              as Range

import           Test.Hedgehog.Compat
import           Test.Hedgehog.Hmatrix
import           Test.Utils.Rnf
import           Test.Grenade.Layers.Internal.Reference

import qualified Numeric.LinearAlgebra.Data   as D
import           Numeric.LinearAlgebra.Static (R)
import qualified Numeric.LinearAlgebra.Static as H

import           Grenade.Core
import           Grenade.Layers.FullyConnected
import           Grenade.Utils.ListStore

data OpaqueFullyConnected :: Type where
     OpaqueFullyConnected :: (KnownNat i, KnownNat o, KnownNat (i * o)) => FullyConnected i o -> OpaqueFullyConnected

instance Show OpaqueFullyConnected where
    show (OpaqueFullyConnected n) = show n

genOpaqueFullyConnected :: Gen OpaqueFullyConnected
genOpaqueFullyConnected = do
  input   :: Integer  <- choose 2 100
  output  :: Integer  <- choose 1 100
  let Just input'      = someNatVal input
  let Just output'     = someNatVal output
  case (input', output') of
    (SomeNat (Proxy :: Proxy i'), SomeNat (Proxy :: Proxy o')) ->
      case (unsafeCoerce (Dict :: Dict ()) :: Dict (KnownNat (i' * o'))) of
        Dict -> do
          wB    <- randomVector
          bM    <- randomVector
          wN    <- uniformSample
          kM    <- uniformSample
          return . OpaqueFullyConnected $ (FullyConnected (FullyConnected' wB wN) (ListStore 0 [Just $ FullyConnected' bM kM]) :: FullyConnected i' o')

prop_fully_connected_run_forwards_backwards :: Property
prop_fully_connected_run_forwards_backwards = property $ do
  OpaqueFullyConnected (fc :: FullyConnected i o) <- blindForAll genOpaqueFullyConnected
  let (FullyConnected (FullyConnected' b w) _) = fc
  input <- forAll randomVector
  let (ntape, nout)  = naiveFullyConnectedRunForwards w b input
      (nnw, nnb, nd) = naiveFullyConnectedBackprop w (ntape, nout)
      (tape, S1D out :: S ('D1 o)) = runForwards fc (S1D input)
      (grad, S1D d   :: S ('D1 i)) = runBackwards fc tape (S1D out)
      FullyConnected' nb nw = grad

  H.extract out  `isSimilarVectorTo` H.extract nout
  H.extract d    `isSimilarVectorTo` H.extract nd
  H.extract nw   `isSimilarMatrixTo` H.extract nnw
  H.extract nb   `isSimilarVectorTo` H.extract nnb
  H.extract tape `isSimilarVectorTo` H.extract ntape

prop_fully_connected_rnf :: Property
prop_fully_connected_rnf = property $ do
  gen <- evalIO create
  fc :: FullyConnected 3 5 <- evalIO $ createRandomWith UniformInit gen
  let (FullyConnected fc'@(FullyConnected' b w) _) = fc

  -- Evaluates the weights and biases
  err0 <- evalIO $ tryEvalRnf (FullyConnected (error expectedErrStr) mkListStore)
  assert $ rnfRaisedErr err0

  -- Evaluates the list store
  err1 <- evalIO $ tryEvalRnf (FullyConnected fc' (error expectedErrStr))
  assert $ rnfRaisedErr err1

  -- Evaluates the biases
  err2 <- evalIO $ tryEvalRnf (FullyConnected (FullyConnected' (error expectedErrStr) w) mkListStore)
  assert $ rnfRaisedErr err2

  -- Evaluates the weights
  err3 <- evalIO $ tryEvalRnf (FullyConnected (FullyConnected' b (error expectedErrStr)) mkListStore)
  assert $ rnfRaisedErr err3

  OpaqueFullyConnected (fcValid :: FullyConnected i o) <- blindForAll genOpaqueFullyConnected
  noErrStr <- evalIO $ tryEvalRnf fcValid
  assert $ rnfNoError noErrStr

prop_fully_connected_is_randomizable :: Property
prop_fully_connected_is_randomizable = withTests 1 $ property $ do
  gen <- evalIO create
  _ :: FullyConnected 3 5 <- evalIO $ createRandomWith UniformInit gen
  success

prop_fully_connected_serialized :: Property
prop_fully_connected_serialized = property $ do
  OpaqueFullyConnected (fc :: FullyConnected i o) <- blindForAll genOpaqueFullyConnected
  let bs = encode fc
  let dec = decode bs :: Either String (FullyConnected i o)
  assert $ isRight dec
  let Right fc' = dec
      (FullyConnected (FullyConnected' b w)   _) = fc
      (FullyConnected (FullyConnected' b' w') _) = fc'
  H.extract b === H.extract b'
  H.extract w === H.extract w'

prop_fully_connected_show :: Property
prop_fully_connected_show = withTests 1 $ property $ do
  OpaqueFullyConnected fc <- blindForAll genOpaqueFullyConnected
  (show fc) `seq` success

prop_fully_connected_update :: Property
prop_fully_connected_update = withTests 1 $ property $ do
  let fc :: FullyConnected 3 5 = createFixed
      input = createFixedInput :: R 3
      (tape, out :: S ('D1 5)) = runForwards fc (S1D input)
      (grad, _   :: S ('D1 3)) = runBackwards fc tape out
      FullyConnected (FullyConnected' b w)   _ = runUpdate defSGD  fc grad
      FullyConnected (FullyConnected' b' w') _ = runUpdate defAdam fc grad
      expW :: H.L 5 3  = H.fromList [ 0.1214, 0.2429, 0.3643
                                    , 0.5229, 0.6171, 0.7114
                                    , 0.9243, 0.9914, 1.0586
                                    , 1.3257, 1.3657, 1.4057
                                    , 1.7271, 1.7400, 1.7529 ]
      expW' :: H.L 5 3 = H.fromList [ 0.1419, 0.2847, 0.4276
                                    , 0.5704, 0.7133, 0.8561
                                    , 0.9990, 1.1419, 1.2847
                                    , 1.4276, 1.5704, 1.7133
                                    , 1.8561, 1.9990, 2.1419 ]
  assert $ allCloseV b  (H.fromList [0.1214,0.2371,0.3529,0.4686,0.5843])
  assert $ allCloseV b' (H.fromList [0.1419,0.2847,0.4276,0.5704,0.7133] )
  assert $ allCloseP (S2D w)  (S2D expW) 0.001
  assert $ allCloseP (S2D w') (S2D expW') 0.001
  where
    createFixed :: forall i o. (KnownNat i, KnownNat o) => FullyConnected i o
    createFixed = fc
      where
        i' = fromIntegral $ natVal (Proxy :: Proxy i)
        o' = fromIntegral $ natVal (Proxy :: Proxy o)
        bias :: H.R o   = H.fromList (map (/7) [1..o'])
        acts :: H.L o i = H.fromList (map (/7) [1..(o' * i')])
        fc :: FullyConnected i o = FullyConnected (FullyConnected' bias acts) mkListStore

    createFixedInput :: forall i. (KnownNat i) => R i
    createFixedInput = H.fromList [1..i']
      where
          i' = fromIntegral $ natVal (Proxy :: Proxy i)

prop_fully_connected_foldable_gradient :: Property
prop_fully_connected_foldable_gradient = property $ do
  OpaqueFullyConnected (fc :: FullyConnected i o) <- blindForAll genOpaqueFullyConnected
  m <- forAll $ genRealNum (Range.constant 0 0.9)
  let FullyConnected fc'@(FullyConnected' b w) _ = fc
      FullyConnected' b' w' = mapGradient (*m) fc'
      bl  = D.toList $ H.extract b
      wl  = concat $ D.toLists $ H.extract w
      bl' = D.toList $ H.extract b'
      wl' = concat $ D.toLists $ H.extract w'
      [sb, sw] = squaredSums fc'
      ssum =  sum . map (** 2)

  (map (*m) bl) `isSimilarListTo` bl'
  (map (*m) wl) `isSimilarListTo` wl'

  ssum bl `isWithinPrecisionOf` sb
  ssum wl `isWithinPrecisionOf` sw
  where
    isWithinPrecisionOf = isWithinOf 0.001

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

ext: hs

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE OverloadedStrings #-}

module Data.Model where

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

import Statistics.Distribution
import Statistics.Distribution.Poisson

import Data.Aeson
import Data.String (IsString)

type Hist = [Double]
type PredHist = [PoissonDistribution]
type DataHist = [Int]

newtype ProcName = ProcName { unPN :: String } deriving (Eq, Ord, FromJSON, IsString)
newtype RegName = RegName { unRN :: String } deriving (Eq, Ord, FromJSON, IsString)

instance Show ProcName where
    show = unPN

instance Show RegName where
    show = unRN

rebinH :: Num a => Int -> [a] -> [a]
rebinH _ []            = []
rebinH n h | n <= 1    = h
           | otherwise = let (hs, hs') = splitAt n h
                         in  sum hs : rebinH n hs'

zipWithLen :: (a -> b -> c) -> [a] -> [b] -> [c]
zipWithLen f (a:as) (b:bs) = f a b : zipWithLen f as bs
zipWithLen _ [] []         = []
zipWithLen _ _ _           = error "zipping lists of different lengths"

-- the log likelihood of a prediction histogram given a data histogram
predHistLLH :: DataHist -> PredHist -> Double
predHistLLH dh ph = sum $ zipWithLen logProbability ph dh

-- a process's normalization is consistent across regions
type Process = Map RegName Hist

-- a collection of named processes
type Prediction = Map ProcName Process

type TotalPrediction = Map RegName PredHist

-- data in several regions
type Dataset = Map RegName DataHist

liftP :: (Double -> Double) -> PoissonDistribution -> PoissonDistribution
liftP f = poisson . f . poissonLambda

liftP2 :: (Double -> Double -> Double) -> PoissonDistribution -> PoissonDistribution -> PoissonDistribution
liftP2 f p p' = poisson $ f (poissonLambda p) (poissonLambda p')

toPred :: Hist -> PredHist
toPred = fmap (poisson . cleanBin)

scaleH :: Double -> Hist -> Hist
scaleH n = fmap (*n)

addH :: Hist -> Hist -> Hist
addH = zipWithLen (+)

cleanBin :: Double -> Double
cleanBin x = if x <= 0 then 1e-100 else x

mulH :: Hist -> Hist -> Hist
mulH = zipWithLen (*)

divH :: Hist -> Hist -> Hist
divH = zipWithLen (/)

sumH :: [Hist] -> Hist
sumH = foldl1 addH

subH :: Hist -> Hist -> Hist
subH = zipWithLen (-)

-- a ModelParam alters a model in some particular way and has a prior
-- distribution
data ModelParam = ModelParam { mpName :: String
                             , mpPrior :: Double -> Double
                             , mpAlter :: Double -> Prediction -> Prediction
                             }

-- alter a particular process's histograms uniformly by some function
alterProc :: (Hist -> Hist) -> ProcName -> Prediction -> Prediction
alterProc f = M.adjust (fmap f)

-- process normalization ModelParam
procNormParam :: ContDistr d => d -> Double -> ProcName -> ModelParam
procNormParam prior n name = ModelParam (show name ++ "_norm") (logDensity prior)
                                $ (\x -> alterProc (scaleH (1+n*x)) name)

procShapeParam :: Double -> Map RegName Hist -> Process -> Process
procShapeParam x s p = M.differenceWith f p s
    where f p' s' = Just $ addH (scaleH x s') p'

shapeParam :: ContDistr d => d -> String -> Map ProcName (Map RegName Hist) -> ModelParam
shapeParam prior name hshapes = ModelParam name (logDensity prior) f
    where f x p = M.differenceWith (g x) p hshapes
          g x p' s' = Just $ procShapeParam x s' p'

totalPrediction :: Prediction -> TotalPrediction
totalPrediction = fmap toPred . M.foldr (M.unionWith addH) M.empty

-- TODO!!!
expectedData :: Prediction -> Dataset
expectedData = fmap (fmap $ round . poissonLambda) . totalPrediction

-- the poisson likelihood of a model given the input data
modelPoissonLLH :: Dataset -> Prediction -> Double
modelPoissonLLH ds m = M.foldr (+) 0 $ M.intersectionWith predHistLLH ds
                                     $ totalPrediction m

modelLLH :: Dataset -> Prediction -> [ModelParam] -> [Double] -> Double
modelLLH ds hpred hparams params = priorLLH + poissLLH
    where
        priorLLH = sum $ zipWithLen mpPrior hparams params
        hpred' = foldr ($) hpred (zipWithLen mpAlter hparams params)
        poissLLH = modelPoissonLLH ds hpred'

max_stars_repo_name: cspollard/HMCMC

max_issues_repo_name: cspollard/HMCMC

max_forks_repo_name: cspollard/HMCMC


text:
{-# LANGUAGE ViewPatterns #-}

module Main where

import Codec.Image.DevIL
import Data.Array.Unboxed
import Math.Probably.MCMC
import Math.Probably.RandIO
import Math.Probably.Sampler
import Math.Probably.FoldingStats
import Control.Monad.State.Strict
import System.Cmd
import System.Environment
import Data.Array.IO
import System.IO

import Data.Array.Unboxed
import Numeric.LinearAlgebra hiding (find)
import qualified Math.Probably.PDF as PDF
import Data.List
import Data.Maybe
import Control.Applicative
import Data.Ord

import Data.IORef
import System.IO.Unsafe

import CVUtils
import Edge

posteriorV :: Image -> Image -> (Int, Int) -> Vector R -> R
posteriorV bgim im (cx,cy) v =
   uniformLogPdf 0 100.0 len1 +
   uniformLogPdf 0 1.0 ecc +
   uniformLogPdf 0 10000.0 noise +
   uniformLogPdf 0 1500.0 px +
   uniformLogPdf 0 1500.0 py +
   sum [ f x y chan |
                   x <- [ (cx::Int) -20.. cx +20],
                   y <- [ cy -20.. cy +20],
                   chan <- [0..2]]

  where px = v @> 0
        py = v @> 1
        noise = v @> 2
        len1 = v@> 3
        rot = v@> 4
        ecc = v@> 5
        red = v@> 6
        green = v@> 7
        blue = v@> 8
        noiseInside = v@> 9
        len2 = v@> 10
        f1x = px+(len1*ecc)*cos rot
        f1y = py+(len1*ecc)*sin rot
        f2x = px-(len1*ecc)*cos rot
        f2y = py-(len1*ecc)*sin rot
        f1x2 = px+(len2*ecc)*cos rot
        f1y2 = py+(len2*ecc)*sin rot
        f2x2 = px-(len2*ecc)*cos rot
        f2y2 = py-(len2*ecc)*sin rot
        colVec = fromList [red,green,blue]
        f :: Int -> Int -> Int -> R
        f x y ch
         = if dist  f1x f1y   x  y  + dist  f2x f2y   x  y  < 2 * len1
              then gaussR noiseInside (colVec@>ch) $ im!(y,x,ch)
              else if dist  f1x2 f1y2   x  y  + dist  f2x2 f2y2   x  y  < 2 * len2
                      then gaussW8 noise ((bgim!(y,x,ch))`div` 2) $ im!(y,x,ch)
                      else gaussW8 noise (bgim!(y,x,ch)) $ im!(y,x,ch)

main = do
     ilInit
     bgnm : fvid :_ <- getArgs
     bgIm <-readImage bgnm
     system $ "~/cvutils/extract "++fvid++" 3200"
     frame0 <-readImage "extract.png"
     let x = 1092
         y = 518
         rot = -5.12
         posterior = posteriorV bgIm frame0 (round x,round y)
         postAndV v = (posteriorV bgIm frame0 (round x,round y) $ fromList v, v)

         initialsV = fromList [x,y,218, 4.5, rot, 0.9, 0.1,0.1,0.1, 218, 6]
     print [frame0!(round y,round x,c) | c <- [0..2]]
     print [frame0!(round y+1,round x,c) | c <- [0..2]]
     print [frame0!(round y,round x+1,c) | c <- [0..2]]
     runRIO $ do
         iniampar <- sample $ initialAdaMet 300 1.5e-4 posterior initialsV
         AMPar v _ _ _ _ _ _ <- runAndDiscard 50000 (show . ampPar) iniampar $ adaMet False posterior
         lift $ print initialsV
         lift $ print v

max_stars_repo_name: glutamate/cvutils

max_issues_repo_name: glutamate/cvutils

max_forks_repo_name: glutamate/cvutils


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE DefaultSignatures #-}

-----------------------------------------------------------------------------
-- |
-- Copyright   :  (C) 2012-2015 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  experimental
-- Portability :  non-portable
--
-- Involutive rings
----------------------------------------------------------------------------
module Linear.Conjugate
  ( Conjugate(..)
  , TrivialConjugate
  ) where

import Data.Complex hiding (conjugate)
import Data.Int
import Data.Word
import Foreign.C.Types (CFloat, CDouble)

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

-- | An involutive ring
class Num a => Conjugate a where
  -- | Conjugate a value. This defaults to the trivial involution.
  --
  -- >>> conjugate (1 :+ 2)
  -- 1.0 :+ (-2.0)
  --
  -- >>> conjugate 1
  -- 1
  conjugate :: a -> a
#ifndef HLINT
  default conjugate :: TrivialConjugate a => a -> a
  conjugate = id
#endif

-- | Requires and provides a default definition such that
--
-- @
-- 'conjugate' = 'id'
-- @
class Conjugate a => TrivialConjugate a

instance Conjugate Integer
instance Conjugate Int
instance Conjugate Int64
instance Conjugate Int32
instance Conjugate Int16
instance Conjugate Int8
instance Conjugate Word
instance Conjugate Word64
instance Conjugate Word32
instance Conjugate Word16
instance Conjugate Word8
instance Conjugate Double
instance Conjugate Float
instance Conjugate CFloat
instance Conjugate CDouble

instance (Conjugate a, RealFloat a) => Conjugate (Complex a) where
  {-# SPECIALIZE instance Conjugate (Complex Float) #-}
  {-# SPECIALIZE instance Conjugate (Complex Double) #-}
  conjugate (a :+ b) = conjugate a :+ negate b

instance TrivialConjugate Integer
instance TrivialConjugate Int
instance TrivialConjugate Int64
instance TrivialConjugate Int32
instance TrivialConjugate Int16
instance TrivialConjugate Int8
instance TrivialConjugate Word
instance TrivialConjugate Word64
instance TrivialConjugate Word32
instance TrivialConjugate Word16
instance TrivialConjugate Word8
instance TrivialConjugate Double
instance TrivialConjugate Float
instance TrivialConjugate CFloat
instance TrivialConjugate CDouble

max_stars_repo_name: ekmett/linear

max_issues_repo_name: ekmett/linear

max_forks_repo_name: ekmett/linear


text:
-- {-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE DataKinds, GADTs, TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FunctionalDependencies #-}

module HBLAS.BLAS.Level3Spec(main, spec) where

import Data.Complex

import Numerical.HBLAS.MatrixTypes as Matrix
import Numerical.HBLAS.BLAS.Level3 as BLAS

import Test.Hspec

main :: IO ()
main = hspec spec

spec :: Spec
spec = do
  gemmSpec
  hemmSpec
  herkSpec
  her2kSpec
  symmSpec
  syrkSpec
  syr2kSpec
  trmmSpec
  trsmSpec

gemmSpec :: Spec
gemmSpec =
  context "?GEMM" $ do
    describe "SGEMM" $ do
      it "2x2 all 1's" $ do
        matmatTest1SGEMM
    describe "DGEMM" $ do
      it "2x2 all 1's" $ do
        matmatTest1DGEMM Matrix.SRow
      it "3x2 and 5x3 all 1's" $ do
        matmatTest2DGEMM Matrix.SRow
      it "2x3^T and 5x3 all 1's" $ do
        matmatTest3DGEMM Matrix.SRow
      it "2x3 and 2x3^T all 1's" $ do
        matmatTest3aDGEMM Matrix.SRow
      it "3x2 and 3x5^T all 1's" $ do
        matmatTest4DGEMM Matrix.SRow
      it "2x3^T and 3x5^T all 1's" $ do
        matmatTest5DGEMM Matrix.SRow
      it "3x2^T and 3x2 all 1's" $ do
        matmatTest6DGEMM Matrix.SRow
      it "2x64^T and 2x64 all 1's" $ do
        matmatTest7DGEMM Matrix.SRow
      it "2x9^T and 2x9 all 1's" $ do
        matmatTest8DGEMM Matrix.SRow
      it "2x2 all 1's (column oriented)" $ do
        matmatTest1DGEMM Matrix.SColumn
      it "3x2 and 5x3 all 1's (column oriented)" $ do
        matmatTest2DGEMM Matrix.SColumn
      it "2x3^T and 5x3 all 1's (column oriented)" $ do
        matmatTest3DGEMM Matrix.SColumn
      it "2x3 and 2x3^T all 1's (column oriented)" $ do
        matmatTest3aDGEMM Matrix.SColumn
      it "3x2 and 3x5^T all 1's (column oriented)" $ do
        matmatTest4DGEMM Matrix.SColumn
      it "2x3^T and 3x5^T all 1's (column oriented)" $ do
        matmatTest5DGEMM Matrix.SColumn
      it "3x2^T and 3x2 all 1's (column oriented)" $ do
        matmatTest6DGEMM Matrix.SColumn
      it "2x64^T and 2x64 all 1's (column oriented)" $ do
        matmatTest7DGEMM Matrix.SColumn
      it "2x9^T and 2x9 all 1's (column oriented)" $ do
        matmatTest8DGEMM Matrix.SColumn
    describe "CGEMM" $ do
      it "2x2 all 1's " $ do
        matmatTest1CGEMM
    describe "ZGEMM" $ do
      it "2x2 all 1's " $ do
        matmatTest1ZGEMM

matmatTest1SGEMM:: IO ()
matmatTest1SGEMM = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const (1.0 :: Float))
    right <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const (1.0 :: Float))
    res   <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const (0.0 :: Float))
    BLAS.sgemm Matrix.NoTranspose Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [2,2,2,2]

matmatTest1DGEMM :: Matrix.SOrientation x -> IO ()
matmatTest1DGEMM or = do
    left  <- Matrix.generateMutableDenseMatrix or (2,2) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix or (2,2) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (2,2) (const 0.0)
    BLAS.dgemm Matrix.NoTranspose Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [2.0,2.0,2.0,2.0]

matmatTest2DGEMM :: Matrix.SOrientation x -> IO ()
matmatTest2DGEMM or = do
    left  <- Matrix.generateMutableDenseMatrix or (3,2) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix or (5,3) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (5,2) (const 0.0)
    BLAS.dgemm Matrix.NoTranspose Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` replicate 10 3

matmatTest3DGEMM :: Matrix.SOrientation x -> IO ()
matmatTest3DGEMM or = do
    left  <- Matrix.generateMutableDenseMatrix or (2,3) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix or (5,3) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (5,2) (const 0.0)
    BLAS.dgemm Matrix.Transpose Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` replicate 10 3

matmatTest3aDGEMM :: Matrix.SOrientation x -> IO ()
matmatTest3aDGEMM or = do
    left  <- Matrix.generateMutableDenseMatrix or (2,3) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix or (2,3) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (3,3) (const 0.0)
    BLAS.dgemm Matrix.NoTranspose Matrix.Transpose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` replicate 9 2

matmatTest4DGEMM :: Matrix.SOrientation x -> IO ()
matmatTest4DGEMM or = do
    left  <- Matrix.generateMutableDenseMatrix or (3,2) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix or (3,5) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (5,2) (const 0.0)
    BLAS.dgemm Matrix.NoTranspose Matrix.Transpose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` replicate 10 3

matmatTest5DGEMM :: Matrix.SOrientation x -> IO ()
matmatTest5DGEMM or = do
    left  <- Matrix.generateMutableDenseMatrix or  (2,3) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix or  (3,5) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or  (5,2) (const 0.0)
    BLAS.dgemm Matrix.Transpose Matrix.Transpose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` replicate 10 3

matmatTest6DGEMM :: Matrix.SOrientation x -> IO ()
matmatTest6DGEMM or = do
    left  <- Matrix.generateMutableDenseMatrix or (3,2) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix or (3,2) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (3,3) (const 0.0)
    BLAS.dgemm Matrix.Transpose Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` replicate 9 2

matmatTest7DGEMM :: Matrix.SOrientation x -> IO ()
matmatTest7DGEMM or = do
    left  <- Matrix.generateMutableDenseMatrix or (2,64) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix or (2,64) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (2,2) (const 0.0)
    BLAS.dgemm Matrix.Transpose Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` replicate 4 64

matmatTest8DGEMM :: Matrix.SOrientation x -> IO ()
matmatTest8DGEMM or = do
    left  <- Matrix.generateMutableDenseMatrix or (2,9) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix or (2,9) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (2,2) (const 0.0)
    BLAS.dgemm Matrix.Transpose Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` replicate 4 9

matmatTest1CGEMM:: IO ()
matmatTest1CGEMM = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 0.0)
    BLAS.cgemm Matrix.NoTranspose Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [2.0,2.0,2.0,2.0]

matmatTest1ZGEMM:: IO ()
matmatTest1ZGEMM = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 0.0)
    BLAS.zgemm Matrix.NoTranspose Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [2.0,2.0,2.0,2.0]

hemmSpec :: Spec
hemmSpec =
  context "?HEMM" $ do
    describe "CHEMM" $ do
      it "3x3 and 2x3 with leftside upper (row oriented)" $ do
        matmatTest1CHEMM
    describe "ZHEMM" $ do
      it "3x3 and 3x2 with rightside lower (column oriented)" $ do
        matmatTest1ZHEMM

-- [1:+0    1:+1    1:+1]   [1 1]   [3:+2    3:+2   ]
-- [1:+(-1) 1:+0    2:+2] * [1 1] = [4:+1    4:+1   ]
-- [1:+(-1) 2:+(-2) 1:+0]   [1 1]   [4:+(-3) 4:+(-3)]
matmatTest1CHEMM :: IO ()
matmatTest1CHEMM = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,3) (\(x, y) -> [1:+0, 1:+1, 1:+1, 0:+0, 1:+0, 2:+2, 0:+0, 0:+0, 1:+0] !! (x + y * 3))
    right <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (const 0.0)
    BLAS.chemm Matrix.LeftSide Matrix.MatUpper 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [3:+2, 3:+2, 4:+1, 4:+1, 4:+(-3), 4:+(-3)]

-- [1 1 1]   [1:+0 1:+(-1) 1:+(-1)]    [3:+2 4:+1 4:+(-3)]
-- [1 1 1] * [1:+1 1:+0    2:+(-2)]  = [3:+2 4:+1 4:+(-3)]
--           [1:+1 2:+2    1:+0   ]
matmatTest1ZHEMM :: IO ()
matmatTest1ZHEMM = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,3) (\(x, y) -> [1:+0, 0:+0, 0:+0, 1:+1, 1:+0, 0:+0, 1:+1, 2:+2, 1:+0] !! (x + y * 3))
    right <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (const 0.0)
    BLAS.zhemm Matrix.RightSide Matrix.MatLower 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [3:+2, 3:+2, 4:+1, 4:+1, 4:+(-3), 4:+(-3)]

herkSpec :: Spec
herkSpec = do
  context "?HERK" $ do
    describe "CHERK" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1CHERK
    describe "ZHERK" $ do
      it "3x3 and 3x2 with lower conjtranspose (column oriented)" $ do
        matmatTest1ZHERK

-- [1 2]   [1 3 5]   [1:+0    1:+1    1:+1]   [6:+0 12:+1 18:+1]
-- [3 4] * [2 4 6] + [1:+(-1) 1:+0    2:+2] = [0:+0 26:+0 41:+2]
-- [5 6]             [1:+(-1) 2:+(-2) 1:+0]   [0:+0 0:+0  62:+0]
matmatTest1CHERK :: IO ()
matmatTest1CHERK = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,3) (\(x, y) -> [1:+0, 1:+1, 1:+1, 0:+0, 1:+0, 2:+2, 0:+0, 0:+0, 1:+0] !! (x + y * 3))
    BLAS.cherk Matrix.MatUpper Matrix.NoTranspose 1.0 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [6:+0, 12:+1, 18:+1, 0:+0, 26:+0, 41:+2, 0:+0, 0:+0, 62:+0]

-- [1:-1 2]   [1:+1 3 5]   [1:+0 1:+(-1) 1:+(-1)]   [7:+0  0:+0  0:+0]
-- [3    4] * [2    4 6] + [1:+1 1:+0    2:+(-2)] = [12:+4 26:+0 0:+0]
-- [5    6]                [1:+1 2:+2    1:+0   ]   [18:+6 41:+2 62:+0]
matmatTest1ZHERK :: IO ()
matmatTest1ZHERK = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (\(x, y) -> [1:+1, 3, 5, 2, 4, 6] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,3) (\(x, y) -> [1:+0, 0:+0, 0:+0, 1:+1, 1:+0, 0:+0, 1:+1, 2:+2, 1:+0] !! (x + y * 3))
    BLAS.zherk Matrix.MatLower Matrix.ConjTranspose 1.0 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [7:+0, 12:+4, 18:+6, 0:+0, 26:+0, 41:+2, 0:+0, 0:+0, 62:+0]

her2kSpec :: Spec
her2kSpec =
  context "?HER2K" $ do
    describe "CHER2K" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1CHER2K
    describe "ZHER2k" $ do
      it "3x3 and 3x2 with lower conjtranspose (column oriented)" $ do
        matmatTest1ZHER2K

-- [1 2]   [1 3 5]       [1:+0    1:+1    1:+1]   [11:+0 23:+1  35:+1]
-- [3 4] * [2 4 6] * 2 + [1:+(-1) 1:+0    2:+2] = [0:+0  51:+0  80:+2]
-- [5 6]                 [1:+(-1) 2:+(-2) 1:+0]   [0:+0  0:+0  123:+0]
matmatTest1CHER2K :: IO ()
matmatTest1CHER2K = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    b <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,3) (\(x, y) -> [1:+0, 1:+1, 1:+1, 0:+0, 1:+0, 2:+2, 0:+0, 0:+0, 1:+0] !! (x + y * 3))
    BLAS.cher2k Matrix.MatUpper Matrix.NoTranspose 1.0 1.0 a b c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [11:+0, 23:+1, 35:+1, 0:+0, 51:+0, 80:+2, 0:+0, 0:+0, 123:+0]

-- [1:-1 2]   [1:+1 3 5]       [1:+0 1:+(-1) 1:+(-1)]   [13:+0  0:+0    0:+0]
-- [3    4] * [2    4 6] * 2 + [1:+1 1:+0    2:+(-2)] = [23:+7  51:+0   0:+0]
-- [5    6]                    [1:+1 2:+2    1:+0   ]   [35:+11 80:+2 123:+0]
matmatTest1ZHER2K :: IO ()
matmatTest1ZHER2K = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (\(x, y) -> [1:+1, 3, 5, 2, 4, 6] !! (x + y * 3))
    b <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (\(x, y) -> [1:+1, 3, 5, 2, 4, 6] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,3) (\(x, y) -> [1:+0, 0:+0, 0:+0, 1:+1, 1:+0, 0:+0, 1:+1, 2:+2, 1:+0] !! (x + y * 3))
    BLAS.zher2k Matrix.MatLower Matrix.ConjTranspose 1.0 1.0 a b c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [13:+0, 23:+7, 35:+11, 0:+0, 51:+0, 80:+2, 0:+0, 0:+0, 123:+0]

symmSpec :: Spec
symmSpec =
  context "?SYMM" $ do
    describe "SSYMM" $ do
      it "2x2 upper all 1's" $ do
        matmatTest1SSYMM
    describe "DSYMM" $ do
      it "2x2 upper all 1's" $ do
        matmatTest1DSYMM Matrix.SRow Matrix.MatUpper
      it "2x2 and 3x2 upper all 1's" $ do
        matmatTest2DSYMM Matrix.SRow Matrix.MatUpper
      it "2x5 and 2x2 upper all 1's" $ do
        matmatTest3DSYMM Matrix.SRow Matrix.MatUpper
      it "2x2 lower all 1's" $ do
        matmatTest1DSYMM Matrix.SRow Matrix.MatLower
      it "2x2 and 3x2 lower all 1's" $ do
        matmatTest2DSYMM Matrix.SRow Matrix.MatLower
      it "2x5 and 2x2 lower all 1's" $ do
        matmatTest3DSYMM Matrix.SRow Matrix.MatLower
      it "2x2 upper all 1's (column oriented)" $ do
        matmatTest1DSYMM Matrix.SColumn Matrix.MatUpper
      it "2x2 and 3x2 upper all 1's (column oriented)" $ do
        matmatTest2DSYMM Matrix.SColumn Matrix.MatUpper
      it "2x5 and 2x2 upper all 1's (column oriented)" $ do
        matmatTest3DSYMM Matrix.SColumn Matrix.MatUpper
      it "2x2 lower all 1's (column oriented)" $ do
        matmatTest1DSYMM Matrix.SColumn Matrix.MatLower
      it "2x2 and 3x2 lower all 1's (column oriented)" $ do
        matmatTest2DSYMM Matrix.SColumn Matrix.MatLower
      it "2x5 and 2x2 lower all 1's (column oriented)" $ do
        matmatTest3DSYMM Matrix.SColumn Matrix.MatLower
    describe "CSYMM" $ do
      it "2x2 all 1's" $ do
        matmatTest1CSYMM
    describe "ZSYMM" $ do
      it "2x2 all 1's" $ do
        matmatTest1ZSYMM

matmatTest1SSYMM:: IO ()
matmatTest1SSYMM = do
    left  <- Matrix.generateMutableUpperTriangular (Matrix.SRow)  (2,2) (const (1.0 :: Float))
    right <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const (1.0 :: Float))
    res   <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const (0.0 :: Float))
    BLAS.ssymm Matrix.LeftSide Matrix.MatUpper 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [2,2,2,2]

matmatTest1DSYMM :: Matrix.SOrientation x -> Matrix.MatUpLo -> IO ()
matmatTest1DSYMM or uplo = do
    left  <- (if uplo == Matrix.MatUpper then Matrix.generateMutableUpperTriangular or (2,2) (const 1.0)
                else Matrix.generateMutableLowerTriangular or (2,2) (const 1.0))
    right <- Matrix.generateMutableDenseMatrix or (2,2) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (2,2) (const 0.0)
    BLAS.dsymm Matrix.LeftSide uplo 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [2.0,2.0,2.0,2.0]

matmatTest2DSYMM :: Matrix.SOrientation x -> Matrix.MatUpLo -> IO ()
matmatTest2DSYMM or uplo = do
    left  <- (if uplo == Matrix.MatUpper then Matrix.generateMutableUpperTriangular or (2,2) (const 1.0)
                else Matrix.generateMutableLowerTriangular or (2,2) (const 1.0))
    right <- Matrix.generateMutableDenseMatrix or (3,2) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (3,2) (const 0.0)
    BLAS.dsymm Matrix.LeftSide uplo 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [2.0,2.0,2.0,2.0,2.0,2.0]

matmatTest3DSYMM :: Matrix.SOrientation x -> Matrix.MatUpLo -> IO ()
matmatTest3DSYMM or uplo = do
    left  <- (if uplo == Matrix.MatUpper then Matrix.generateMutableUpperTriangular or (2,2) (const 1.0)
                else Matrix.generateMutableLowerTriangular or (2,2) (const 1.0))
    right <- Matrix.generateMutableDenseMatrix or (2,5) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix or (2,5) (const 0.0)
    BLAS.dsymm Matrix.RightSide uplo 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` replicate 10 2

matmatTest1CSYMM:: IO ()
matmatTest1CSYMM = do
    left  <- Matrix.generateMutableUpperTriangular (Matrix.SRow)  (2,2) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 0.0)
    BLAS.csymm Matrix.LeftSide Matrix.MatUpper 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [2.0,2.0,2.0,2.0]

matmatTest1ZSYMM:: IO ()
matmatTest1ZSYMM = do
    left  <- Matrix.generateMutableUpperTriangular (Matrix.SRow)  (2,2) (const 1.0)
    right <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 1.0)
    res   <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (const 0.0)
    BLAS.zsymm Matrix.LeftSide Matrix.MatUpper 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
    resList `shouldBe` [2.0,2.0,2.0,2.0]

syrkSpec :: Spec
syrkSpec =
  context "?SYRK" $ do
    describe "SSYRK" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1SSYRK
    describe "DSYRK" $ do
      it "3x3 and 3x2 with lower transpose (column oriented)" $ do
        matmatTest1DSYRK
    describe "CSYRK" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1CSYRK
    describe "ZSYRK" $ do
      it "3x3 and 3x2 with lower transpose (column oriented)" $ do
        matmatTest1ZSYRK

-- [1 2]   [1 3 5]   [1 1 1]   [6 12 18]
-- [3 4] * [2 4 6] + [1 1 2] = [0 26 41]
-- [5 6]             [1 2 1]   [0 0  62]
matmatTest1SSYRK :: IO ()
matmatTest1SSYRK = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,3) (\(x, y) -> [1, 1, 1, 0, 1, 2, 0, 0, 1] !! (x + y * 3))
    BLAS.ssyrk Matrix.MatUpper Matrix.NoTranspose 1.0 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [6, 12, 18, 0, 26, 41, 0, 0, 62]

-- [1 2]   [1 3 5]   [1 1 1]   [6  0  0 ]
-- [3 4] * [2 4 6] + [1 1 2] = [12 26 0 ]
-- [5 6]             [1 2 1]   [18 41 62]
matmatTest1DSYRK :: IO ()
matmatTest1DSYRK = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (\(x, y) -> [1, 3, 5, 2, 4, 6] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,3) (\(x, y) -> [1, 0, 0, 1, 1, 0, 1, 2, 1] !! (x + y * 3))
    BLAS.dsyrk Matrix.MatLower Matrix.Transpose 1.0 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [6, 12, 18, 0, 26, 41, 0, 0, 62]

-- [1 2]   [1 3 5]   [1:+0 1:+1 1:+1]   [6:+0 12:+1 18:+1]
-- [3 4] * [2 4 6] + [1:+1 1:+0 2:+2] = [0:+0 26:+0 41:+2]
-- [5 6]             [1:+1 2:+2 1:+0]   [0:+0 0:+0  62:+0]
matmatTest1CSYRK :: IO ()
matmatTest1CSYRK = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,3) (\(x, y) -> [1:+0, 1:+1, 1:+1, 0:+0, 1:+0, 2:+2, 0:+0, 0:+0, 1:+0] !! (x + y * 3))
    BLAS.csyrk Matrix.MatUpper Matrix.NoTranspose 1.0 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [6:+0, 12:+1, 18:+1, 0:+0, 26:+0, 41:+2, 0:+0, 0:+0, 62:+0]

-- [1:+1 2]   [1:+1 3 5]   [1:+0 1:+1 1:+1]   [6:+0  0:+0  0:+0]
-- [3    4] * [2    4 6] + [1:+1 1:+0 2:+2] = [12:+1 26:+0 0:+0]
-- [5    6]                [1:+1 2:+2 1:+0]   [18:+1 41:+2 62:+0]
matmatTest1ZSYRK :: IO ()
matmatTest1ZSYRK = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (\(x, y) -> [1:+1, 3, 5, 2, 4, 6] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,3) (\(x, y) -> [1:+0, 0:+0, 0:+0, 1:+1, 1:+0, 0:+0, 1:+1, 2:+2, 1:+0] !! (x + y * 3))
    BLAS.zsyrk Matrix.MatLower Matrix.Transpose 1.0 1.0 a c -- TODO: Matrix.ConjTranspose is invalid to pass to cblas_zsyrk
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [5:+2, 12:+4, 18:+6, 0:+0, 26:+0, 41:+2, 0:+0, 0:+0, 62:+0]

syr2kSpec :: Spec
syr2kSpec =
  context "?SYR2K" $ do
    describe "SSYR2K" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1SSYR2K
    describe "DSYR2K" $ do
      it "3x3 and 3x2 with lower transpose (column oriented)" $ do
        matmatTest1DSYR2K
    describe "CSYR2K" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1CSYR2K
    describe "ZSYR2K" $ do
      it "3x3 and 3x2 with lower conjtranspose (column oriented)" $ do
        matmatTest1ZSYR2K

-- [1 2]   [1 3 5]       [1 1 1]   [11 23  35]
-- [3 4] * [2 4 6] * 2 + [1 1 2] = [0  51  80]
-- [5 6]                 [1 2 1]   [0  0  123]
matmatTest1SSYR2K :: IO ()
matmatTest1SSYR2K = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    b <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,3) (\(x, y) -> [1, 1, 1, 0, 1, 2, 0, 0, 1] !! (x + y * 3))
    BLAS.ssyr2k Matrix.MatUpper Matrix.NoTranspose 1.0 1.0 a b c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [11, 23, 35, 0, 51, 80, 0, 0, 123]

-- [1 2]   [1 3 5]       [1 1 1]   [11 0  0  ]
-- [3 4] * [2 4 6] * 2 + [1 1 2] = [23 51 0  ]
-- [5 6]                 [1 2 1]   [35 80 123]
matmatTest1DSYR2K :: IO ()
matmatTest1DSYR2K = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (\(x, y) -> [1, 3, 5, 2, 4, 6] !! (x + y * 3))
    b <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (\(x, y) -> [1, 3, 5, 2, 4, 6] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,3) (\(x, y) -> [1, 0, 0, 1, 1, 0, 1, 2, 1] !! (x + y * 3))
    BLAS.dsyr2k Matrix.MatLower Matrix.Transpose 1.0 1.0 a b c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [11, 23, 35, 0, 51, 80, 0, 0, 123]

-- [1 2]   [1 3 5]       [1:+0 1:+1 1:+1]   [11:+0 23:+1 35:+1 ]
-- [3 4] * [2 4 6] * 2 + [1:+1 1:+0 2:+2] = [0:+0  51:+0 80:+2 ]
-- [5 6]                 [1:+1 2:+2 1:+0]   [0:+0  0:+0  123:+0]
matmatTest1CSYR2K :: IO ()
matmatTest1CSYR2K = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    b <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,3) (\(x, y) -> [1:+0, 1:+1, 1:+1, 0:+0, 1:+0, 2:+2, 0:+0, 0:+0, 1:+0] !! (x + y * 3))
    BLAS.csyr2k Matrix.MatUpper Matrix.NoTranspose 1.0 1.0 a b c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [11:+0, 23:+1, 35:+1, 0:+0, 51:+0, 80:+2, 0:+0, 0:+0, 123:+0]

-- [1:+1 2]   [1:+1 3 5]       [1:+0 1:+1 1:+1]   [9:+4   0:+0  0:+0  ]
-- [3    4] * [2    4 6] * 2 + [1:+1 1:+0 2:+2] = [23:+7  51:+0 0:+0  ]
-- [5    6]                    [1:+1 2:+2 1:+0]   [35:+11 80:+2 123:+0]
matmatTest1ZSYR2K :: IO ()
matmatTest1ZSYR2K = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (\(x, y) -> [1:+1, 3, 5, 2, 4, 6] !! (x + y * 3))
    b <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,2) (\(x, y) -> [1:+1, 3, 5, 2, 4, 6] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3,3) (\(x, y) -> [1:+0, 0:+0, 0:+0, 1:+1, 1:+0, 0:+0, 1:+1, 2:+2, 1:+0] !! (x + y * 3))
    BLAS.zsyr2k Matrix.MatLower Matrix.Transpose 1.0 1.0 a b c -- TODO: Matrix.ConjTranspose is invalid to pass to cblas_zsyr2k
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [9:+4, 23:+7, 35:+11, 0:+0, 51:+0, 80:+2, 0:+0, 0:+0, 123:+0]

trmmSpec :: Spec
trmmSpec =
  context "?TRMM" $ do
    describe "STRMM" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1STRMM
    describe "DTRMM" $ do
      it "3x3 and 3x2 with lower transpose (column oriented)" $ do
        matmatTest1DTRMM
    describe "CTRMM" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1CTRMM
    describe "ZTRMM" $ do
      it "3x3 and 3x2 with lower conjtranspose (column oriented)" $ do
        matmatTest1ZTRMM

-- [1 1 1]   [1 4]   [6 15]
-- [0 1 2] * [2 5] = [8 17]
-- [0 0 1]   [3 6]   [3 6 ]
matmatTest1STRMM :: IO ()
matmatTest1STRMM = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3, 3) (\(x, y) -> [1, 1, 1, 0, 1, 2, 0, 0, 1] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2, 3) (\(x, y) -> [1, 4, 2, 5, 3, 6] !! (x + y * 2))
    BLAS.strmm Matrix.LeftSide Matrix.MatUpper Matrix.NoTranspose Matrix.MatUnit 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [6, 15, 8, 17, 3, 6]

-- [1 2 3]   [2 1 1]   [2 5  11]
-- [4 5 6] * [0 2 2] = [8 14 26]
--           [0 0 2]
matmatTest1DTRMM :: IO ()
matmatTest1DTRMM = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 3) (\(x, y) -> [2, 0, 0, 1, 2, 0, 1, 2, 2] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 2) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 3))
    BLAS.dtrmm Matrix.RightSide Matrix.MatLower Matrix.Transpose Matrix.MatNonUnit 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [2, 8, 5, 14, 11, 26]

-- [1:+0 1:+1 1:+1]   [1 2]   [9:+8   12:+10]
-- [0:+0 1:+0 2:+2] * [3 4] = [13:+10 16:+12]
-- [0:+0 0:+0 1:+0]   [5 6]   [5:+0   6:+0  ]
matmatTest1CTRMM :: IO ()
matmatTest1CTRMM = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3, 3) (\(x, y) -> [1:+0, 1:+1, 1:+1, 0:+0, 1:+0, 2:+2, 0:+0, 0:+0, 1:+0] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2, 3) (\(x, y) -> [1, 2, 3, 4, 5, 6] !! (x + y * 2))
    BLAS.ctrmm Matrix.LeftSide Matrix.MatUpper Matrix.NoTranspose Matrix.MatUnit 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [9:+8, 12:+10, 13:+10, 16:+12, 5:+0, 6:+0]

-- [1:+1 3 5]   [1:+0 1:+1 1:+1]   [1:+1 3:+2 11:+8 ]
-- [2    4 6] + [0:+0 1:+0 2:+2] = [2:+0 6:+2 16:+10]
--              [0:+0 0:+0 1:+0]
matmatTest1ZTRMM :: IO ()
matmatTest1ZTRMM = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 3) (\(x, y) -> [1:+0, 0:+0, 0:+0, 1:+(-1), 1:+0, 0:+0, 1:+(-1), 2:+(-2), 1:+0] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 2) (\(x, y) -> [1:+1, 3, 5, 2, 4, 6] !! (x + y * 3))
    BLAS.ztrmm Matrix.RightSide Matrix.MatLower Matrix.ConjTranspose Matrix.MatNonUnit 1.0 a c -- TODO: Matrix.ConjTranspose is invalid to pass to cblas_zsyr2k
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [1:+1, 2:+0, 3:+2, 6:+2, 11:+8, 16:+10]

trsmSpec :: Spec
trsmSpec =
  context "?TRSM" $ do
    describe "STRSM" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1STRSM
    describe "DTRSM" $ do
      it "3x3 and 3x2 with lower transpose (column oriented)" $ do
        matmatTest1DTRSM
    describe "CTRSM" $ do
      it "3x3 and 2x3 with upper no transpose (row oriented)" $ do
        matmatTest1CTRSM
    describe "ZTRSM" $ do
      it "3x3 and 3x2 with lower conjtranspose (column oriented)" $ do
        matmatTest1ZTRSM

-- [1 1 1]   [1 4]   [6 15]
-- [0 1 2] * [2 5] = [8 17]
-- [0 0 1]   [3 6]   [3 6 ]
matmatTest1STRSM :: IO ()
matmatTest1STRSM = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3, 3) (\(x, y) -> [1, 1, 1, 0, 1, 2, 0, 0, 1] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2, 3) (\(x, y) -> [6, 15, 8, 17, 3, 6] !! (x + y * 2))
    BLAS.strsm Matrix.LeftSide Matrix.MatUpper Matrix.NoTranspose Matrix.MatUnit 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [1, 4, 2, 5, 3, 6]

-- [1 2 3]   [2 1 1]   [2 5  11]
-- [4 5 6] * [0 2 2] = [8 14 26]
--           [0 0 2]
matmatTest1DTRSM :: IO ()
matmatTest1DTRSM = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 3) (\(x, y) -> [2, 0, 0, 1, 2, 0, 1, 2, 2] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 2) (\(x, y) -> [2, 5, 11, 8, 14, 26] !! (x + y * 3))
    BLAS.dtrsm Matrix.RightSide Matrix.MatLower Matrix.Transpose Matrix.MatNonUnit 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [1, 4, 2, 5, 3, 6]

-- [1:+0 1:+1 1:+1]   [1 2]   [9:+8   12:+10]
-- [0:+0 1:+0 2:+2] * [3 4] = [13:+10 16:+12]
-- [0:+0 0:+0 1:+0]   [5 6]   [5:+0   6:+0  ]
matmatTest1CTRSM :: IO ()
matmatTest1CTRSM = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3, 3) (\(x, y) -> [1:+0, 1:+1, 1:+1, 0:+0, 1:+0, 2:+2, 0:+0, 0:+0, 1:+0] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2, 3) (\(x, y) -> [9:+8, 12:+10, 13:+10, 16:+12, 5:+0, 6:+0] !! (x + y * 2))
    BLAS.ctrsm Matrix.LeftSide Matrix.MatUpper Matrix.NoTranspose Matrix.MatUnit 1.0 a c
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [1, 2, 3, 4, 5, 6]

-- [1:+1 3 5]   [1:+0 1:+1 1:+1]   [1:+1 3:+2 11:+8 ]
-- [2    4 6] + [0:+0 1:+0 2:+2] = [2:+0 6:+2 16:+10]
--              [0:+0 0:+0 1:+0]
matmatTest1ZTRSM :: IO ()
matmatTest1ZTRSM = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 3) (\(x, y) -> [1:+0, 0:+0, 0:+0, 1:+(-1), 1:+0, 0:+0, 1:+(-1), 2:+(-2), 1:+0] !! (x + y * 3))
    c <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 2) (\(x, y) -> [1:+1, 3:+2, 11:+8, 2:+0, 6:+2, 16:+10] !! (x + y * 3))
    BLAS.ztrsm Matrix.RightSide Matrix.MatLower Matrix.ConjTranspose Matrix.MatNonUnit 1.0 a c -- TODO: Matrix.ConjTranspose is invalid to pass to cblas_zsyr2k
    resList <- Matrix.mutableVectorToList $ _bufferDenMutMat c
    resList `shouldBe` [1:+1, 2, 3, 4, 5, 6]

max_stars_repo_name: schnecki/hblas

max_issues_repo_name: schnecki/hblas

max_forks_repo_name: schnecki/hblas


text:
{-# LANGUAGE BangPatterns         #-}
{-# LANGUAGE CPP                  #-}
{-# LANGUAGE DataKinds            #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE GADTs                #-}
{-# LANGUAGE RankNTypes           #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE TypeOperators        #-}
{-# LANGUAGE UndecidableInstances #-}
{-|
Module      : Grenade.Core.Shape
Description : Dependently typed shapes of data which are passed between layers of a network
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental

-}
module Grenade.Core.Shape (
    S (..)
  , Shape (..)
#if MIN_VERSION_singletons(2,6,0)
  , SShape (..)
#else
  , Sing (..)
#endif

  , randomOfShape
  , fromStorable
  , fromStorableV
  , isVectorShape
  ) where

#if MIN_VERSION_singletons(2,6,0)
import           Data.Kind                    (Type)
import           Data.Singletons.TypeLits     (SNat (..))
#endif

import           Control.DeepSeq              (NFData (..))
import           Data.Maybe                   (fromMaybe)
import           Data.Proxy
import           Data.Serialize
import           Data.Singletons
import           Data.Singletons.TypeLits
import           Data.Vector.Storable         (Vector)
import qualified Data.Vector.Storable         as V
import           GHC.TypeLits                 hiding (natVal)
import qualified Numeric.LinearAlgebra        as NLA
import           Numeric.LinearAlgebra.Static hiding (zipWithVector)
import qualified Numeric.LinearAlgebra.Static as H
import           System.Random.MWC
import           Unsafe.Coerce                (unsafeCoerce)

import           Debug.Trace

import           Grenade.Types
import           Grenade.Utils.Vector

-- | The current shapes we accept.
--   at the moment this is just one, two, and three dimensional
--   Vectors/Matricies.
--
--   These are only used with DataKinds, as Kind `Shape`, with Types 'D1, 'D2, 'D3.
data Shape
  = D1 Nat
  -- ^ One dimensional vector
  | D2 Nat Nat
  -- ^ Two dimensional matrix. Row, Column.
  | D3 Nat Nat Nat
  -- ^ Three dimensional matrix. Row, Column, Channels.

-- | Concrete data structures for a Shape.
--
--   All shapes are held in contiguous memory.
--   3D is held in a matrix (usually row oriented) which has height depth * rows.
data S (n :: Shape) where
  S1D :: ( KnownNat len )
      => R len
      -> S ('D1 len)

  S2D :: ( KnownNat rows, KnownNat columns )
      => L rows columns
      -> S ('D2 rows columns)

  S3D :: ( KnownNat rows
         , KnownNat columns
         , KnownNat depth
         , KnownNat (rows * depth))
      => L (rows * depth) columns
      -> S ('D3 rows columns depth)

  -- HMatrix instances
  S1DV :: (KnownNat len) -- ^ Always a row-vector
       => !(V.Vector RealNum)
       -> S ('D1 len)

  -- HMatrix instances
  S2DV :: (KnownNat rows, KnownNat columns)
       => !(V.Vector RealNum) -- ^ Vector in row or column major accoring to @order@ in @Grenade.Utils.Conversion@. Use Data.Matrix here?
       -> S ('D2 rows columns)

  -- HMatrix instances
  -- S3DH :: (KnownNat rows, KnownNat columns, KnownNat depth, KnownNat (rows * depth))
  --      => HBLAS.MDenseMatrix RealWorld 'HBLAS.Row Double
  --      -> S ('D3 rows columns depth)

isVectorShape :: S x -> Bool
isVectorShape S1DV{} = True
isVectorShape S2DV{} = True
isVectorShape _      = False

instance Show (S n) where
  show (S1D x) = "S1D " ++ show x
  show (S2D x) = "S2D " ++ show x
  show (S3D x) = "S3D " ++ show x
  show (S1DV x) = "S1DV " ++ show x
  show inp@(S2DV x) = "S2DV" ++ show (sz inp) ++ " " ++ show x
    where
      sz ::
           forall rows cols. (KnownNat rows, KnownNat cols)
        => (S ('D2 rows cols))
        -> (Int, Int)
      sz _ =
        let rows = fromIntegral $ natVal (Proxy :: Proxy rows)
            cols = fromIntegral $ natVal (Proxy :: Proxy cols)
         in (rows, cols)

-- Singleton instances.
--
-- These could probably be derived with template haskell, but this seems
-- clear and makes adding the KnownNat constraints simple.
-- We can also keep our code TH free, which is great.
#if MIN_VERSION_singletons(2,6,0)
-- In singletons 2.6 Sing switched from a data family to a type family.
type instance Sing = SShape

data SShape :: Shape -> Type where
  D1Sing :: Sing a -> SShape ('D1 a)
  D2Sing :: Sing a -> Sing b -> SShape ('D2 a b)
  D3Sing :: KnownNat (a * c) => Sing a -> Sing b -> Sing c -> SShape ('D3 a b c)
#else
data instance Sing (n :: Shape) where
  D1Sing :: Sing a -> Sing ('D1 a)
  D2Sing :: Sing a -> Sing b -> Sing ('D2 a b)
  D3Sing :: KnownNat (a * c) => Sing a -> Sing b -> Sing c -> Sing ('D3 a b c)
#endif

instance KnownNat a => SingI ('D1 a) where
  sing = D1Sing sing
instance (KnownNat a, KnownNat b) => SingI ('D2 a b) where
  sing = D2Sing sing sing
instance (KnownNat a, KnownNat b, KnownNat c, KnownNat (a * c)) => SingI ('D3 a b c) where
  sing = D3Sing sing sing sing

instance SingI x => Num (S x) where
  (+) = n2 (+)
  (-) = n2 (-)
  (*) = n2 (*)
  abs = n1 abs
  signum = n1 signum
  fromInteger x = nk (fromInteger x)

instance SingI x => Fractional (S x) where
  (/) = n2 (/)
  recip = n1 recip
  fromRational x = nk (fromRational x)

instance SingI x => Floating (S x) where
  pi = nk pi
  exp = n1 exp
  log = n1 log
  sqrt = n1 sqrt
  (**) = n2 (**)
  logBase = n2 logBase
  sin = n1 sin
  cos = n1 cos
  tan = n1 tan
  asin = n1 asin
  acos = n1 acos
  atan = n1 atan
  sinh = n1 sinh
  cosh = n1 cosh
  tanh = n1 tanh
  asinh = n1 asinh
  acosh = n1 acosh
  atanh = n1 atanh

--
-- I haven't made shapes strict, as sometimes they're not needed
-- (the last input gradient back for instance)
--
instance NFData (S x) where
  rnf (S1D x)   = rnf x
  rnf (S2D x)   = rnf x
  rnf (S3D x)   = rnf x
  rnf (S1DV !v) = rnf v
  rnf (S2DV !v) = rnf v
  -- rnf (S3DH !_) = ()

-- | Generate random data of the desired shape
randomOfShape :: forall x . (SingI x) => IO (S x)
randomOfShape = do
  seed :: Int <- withSystemRandom . asGenST $ \gen -> uniform gen
  return $ case (sing :: Sing x) of
    D1Sing SNat ->
        S1D (randomVector seed Uniform * 2 - 1)

    D2Sing SNat SNat ->
        S2D (uniformSample seed (-1) 1)

    D3Sing SNat SNat SNat ->
        S3D (uniformSample seed (-1) 1)

-- | Generate a shape from a Storable Vector.
--
--   Returns Nothing if the vector is of the wrong size.
fromStorable :: forall x. SingI x => Vector RealNum -> Maybe (S x)
fromStorable xs =
  case (sing :: Sing x) of
    D1Sing SNat           -> S1D <$> H.create xs
    D2Sing SNat SNat      -> S2D <$> mkL xs
    D3Sing SNat SNat SNat -> S3D <$> mkL xs
  where
    mkL ::
         forall rows columns. (KnownNat rows, KnownNat columns)
      => Vector RealNum
      -> Maybe (L rows columns)
    mkL v =
      let rows = fromIntegral $ natVal (Proxy :: Proxy rows)
          columns = fromIntegral $ natVal (Proxy :: Proxy columns)
       in if rows * columns == V.length v
            then H.create $ NLA.reshape columns v
            else Nothing

fromStorableV :: forall x. SingI x => Vector RealNum -> S x
fromStorableV v =
  case (sing :: Sing x) of
    D1Sing SNat           -> S1DV v
    D2Sing SNat SNat      -> S2DV v
    D3Sing SNat SNat SNat -> error "unexpected case in fromStorableS"

instance SingI x => Serialize (S x) where
  put i =
    case i of
      S1D x -> put (1 :: Int) >> (putListOf put . NLA.toList . H.extract $ x)
      S2D x -> put (1 :: Int) >> (putListOf put . NLA.toList . NLA.flatten . H.extract $ x)
      S3D x -> put (1 :: Int) >> (putListOf put . NLA.toList . NLA.flatten . H.extract $ x)
      S1DV x -> put (2 :: Int) >> putListOf put (V.toList x)
      S2DV x -> put (2 :: Int) >> putListOf put (V.toList x)
  get = do
    (nr :: Int) <- get
    case nr of
      1 -> do
        Just i <- fromStorable . V.fromList <$> getListOf get
        return i
      2 -> fromStorableV . V.fromList <$> getListOf get
      _ -> error "unexpected case in get in Serialize instance in Shape.hs"

-- Helper function for creating the number instances
n1 :: ( forall a. Floating a => a -> a ) -> S x -> S x
n1 f (S1D x)  = S1D (f x)
n1 f (S2D x)  = S2D (f x)
n1 f (S3D x)  = S3D (f x)
n1 f (S1DV x) = S1DV (mapVector f x)
n1 f (S2DV x) = S2DV (mapVector f x)

-- helper function for creating the number instances
n2 :: ( forall a. Floating a => a -> a -> a ) -> S x -> S x -> S x
n2 f (S1D x) (S1D y)    = S1D (f x y)
n2 f (S2D x) (S2D y)    = S2D (f x y)
n2 f (S3D x) (S3D y)    = S3D (f x y)
n2 f (S1DV x) (S1DV y)  = S1DV (zipWithVector f x y)
n2 f (S2DV x) (S2DV y)  = S2DV (zipWithVector f x y)
n2 f (S1D x) y@S1DV {}  = n2 f (S1DV $ extract x) y
n2 f x@S1DV {} (S1D y)  = n2 f x (S1DV $ unsafeCoerce y)
n2 f x@S2D{} y@(S2DV _) = n2 f (toS2DV x) y
n2 f x@(S2DV _) y@S2D{} = n2 f x (toS2DV y)

toS2DV :: S ('D2 i j) -> S ('D2 i j)
toS2DV (S2D x)  = S2DV $ V.concat $ map H.extract . H.toColumns $ x
toS2DV (S2DV x) = S2DV x

-- Helper function for creating the number instances
nk :: forall x. SingI x => RealNum -> S x
nk x = case (sing :: Sing x) of
  D1Sing SNat ->
    S1D (konst x)

  D2Sing SNat SNat ->
    S2D (konst x)

  D3Sing SNat SNat SNat ->
    S3D (konst x)

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
-- some tests of the interface for pure
-- computations with inplace updates

import Numeric.LinearAlgebra
import Data.Packed.ST
import Data.Packed.Convert

import Data.Array.Unboxed
import Data.Array.ST
import Control.Monad.ST
import Control.Monad

main = sequence_[
    print test1,
    print test2,
    print test3,
    print test4,
    test5,
    test6,
    print test7,
    test8,
    test0]

-- helper functions
vector l = fromList l :: Vector Double
norm v = pnorm PNorm2 v

-- hmatrix vector and matrix
v = vector [1..10]
m = (5><10) [1..50::Double]

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

-- vector creation by in-place updates on a copy of the argument
test1 = fun v

fun :: Element t => Vector t -> Vector t
fun x = runSTVector $ do
    a <- thawVector x
    mapM_ (flip (modifyVector a) (+57)) [0 .. dim x `div` 2 - 1]
    return a

-- another example: creation of an antidiagonal matrix from a list
test2 = antiDiag 5 8 [1..] :: Matrix Double

antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b
antiDiag r c l = runSTMatrix $ do
    m <- newMatrix 0 r c
    let d = min r c - 1
    sequence_ $ zipWith (\i v -> writeMatrix m i (c-1-i) v) [0..d] l
    return m

-- using vector or matrix functions on mutable objects requires freezing:
test3 = g1 v

g1 x = runST $ do
    a <- thawVector x
    writeVector a (dim x -1) 0
    b <- freezeVector a
    return (norm b)

--  another possibility:
test4 = g2 v

g2 x = runST $ do
    a <- thawVector x
    writeVector a (dim x -1) 0
    t <- liftSTVector norm a
    return t

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

-- haskell arrays
hv = listArray (0,9) [1..10::Double]
hm = listArray ((0,0),(4,9)) [1..50::Double]

-- conversion from standard Haskell arrays
test5 = do
    print $ norm (vectorFromArray hv)
    print $ norm v
    print $ rcond (matrixFromArray hm)
    print $ rcond m

-- conversion to mutable ST arrays
test6 = do
    let y = clearColumn m 1
    print y
    print (matrixFromArray y)

clearColumn x c = runSTUArray $ do
    a <- mArrayFromMatrix x
    forM_ [0..rows x-1] $ \i->
        writeArray a (i,c) (0::Double)
    return a

-- hmatrix functions applied to mutable ST arrays
test7 = unitary (listArray (1,4) [3,5,7,2] :: UArray Int Double)

unitary v = runSTUArray $ do
    a <- thaw v
    n <- norm `fmap` vectorFromMArray a
    b <- mapArray (/n) a
    return b

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

-- (just to check that they are not affected)
test0 = do
    print v
    print m
    --print hv
    --print hm

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

histogram n ds = runSTVector $ do
    h <- newVector (0::Double) n -- number of bins
    let inc k = modifyVector h k (+1)
    mapM_ inc ds
    return h

-- check that newVector is really called with a fresh new array
histoCheck ds = runSTVector $ do
    h <- newVector (0::Double) 15 -- > constant for this test
    let inc k = modifyVector h k (+1)
    mapM_ inc ds
    return h

hc = fromList [1 .. 15::Double]

-- check that thawVector creates a new array
histoCheck2 ds = runSTVector $ do
    h <- thawVector hc
    let inc k = modifyVector h k (+1)
    mapM_ inc ds
    return h

test8 = do
    let ds = [0..14]
    print $ histogram 15 ds
    print $ histogram 15 ds
    print $ histogram 15 ds
    print $ histoCheck ds
    print $ histoCheck ds
    print $ histoCheck ds
    print $ histoCheck2 ds
    print $ histoCheck2 ds
    print $ histoCheck2 ds
    putStrLn "----------------------"

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
module Network.Network.SOM
( SOM(..)

) where

import Network.Neuron
import Network.Layer
import Network.Network

import Numeric.LinearAlgebra

data SOM = SOM { map :: Matrix (Vector Double) }

instance Network (SOM) where
  predict inputs network = inputs

  createNetwork transformation g defs = SOM $ fromLists [[fromList [1.0] :: Vector Double]]

max_stars_repo_name: AkatsukiSirius/LambdaNet

max_issues_repo_name: world-admin/LambdaNet

max_forks_repo_name: world-admin/LambdaNet


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.Normal
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The normal distribution.  This is a continuous probability
-- distribution that describes data that cluster around a mean.

module Statistics.Distribution.Normal
    (
      NormalDistribution
    -- * Constructors
    -- , normalDistr
    --, normalDistrE
    , standard
    ) where

import Data.Data             (Data, Typeable)
import GHC.Generics          (Generic)
import Numeric.MathFunctions.Constants (m_sqrt_2, m_sqrt_2_pi)
import Numeric.SpecFunctions (erfc, invErfc)

import qualified Statistics.Distribution as D
import Statistics.Internal

-- | The normal distribution.
data NormalDistribution = ND {
      mean       :: {-# UNPACK #-} !Double
    , stdDev     :: {-# UNPACK #-} !Double
    , ndPdfDenom :: {-# UNPACK #-} !Double
    , ndCdfDenom :: {-# UNPACK #-} !Double
    } deriving (Eq, Typeable, Data, Generic)

instance Show NormalDistribution where
  showsPrec i (ND m s _ _) = defaultShow2 "normalDistr" m s i
instance Read NormalDistribution where
  readPrec = defaultReadPrecM2 "normalDistr" normalDistrE

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

instance D.ContDistr NormalDistribution where
    logDensity    = logDensity
    quantile      = quantile
    complQuantile = complQuantile

-- | Standard normal distribution with mean equal to 0 and variance equal to 1
standard :: NormalDistribution
standard = ND { mean       = 0.0
              , stdDev     = 1.0
              , ndPdfDenom = log m_sqrt_2_pi
              , ndCdfDenom = m_sqrt_2
              }

-- | Create normal distribution from parameters.
--
-- IMPORTANT: prior to 0.10 release second parameter was variance not
-- standard deviation.
normalDistrE :: Double            -- ^ Mean of distribution
             -> Double            -- ^ Standard deviation of distribution
             -> Maybe NormalDistribution
normalDistrE m sd
  | sd > 0    = Just ND { mean       = m
                        , stdDev     = sd
                        , ndPdfDenom = log $ m_sqrt_2_pi * sd
                        , ndCdfDenom = m_sqrt_2 * sd
                        }
  | otherwise = Nothing

logDensity :: NormalDistribution -> Double -> Double
logDensity d x = (-xm * xm / (2 * sd * sd)) - ndPdfDenom d
    where xm = x - mean d
          sd = stdDev d

cumulative :: NormalDistribution -> Double -> Double
cumulative d x = erfc ((mean d - x) / ndCdfDenom d) / 2

complCumulative :: NormalDistribution -> Double -> Double
complCumulative d x = erfc ((x - mean d) / ndCdfDenom d) / 2

quantile :: NormalDistribution -> Double -> Double
quantile d p
  | p == 0         = -inf
  | p == 1         = inf
  | p == 0.5       = mean d
  | p > 0 && p < 1 = x * ndCdfDenom d + mean d
  | otherwise      =
    error $ "Statistics.Distribution.Normal.quantile: p must be in [0,1] range. Got: "++show p
  where x          = - invErfc (2 * p)
        inf        = 1/0

complQuantile :: NormalDistribution -> Double -> Double
complQuantile d p
  | p == 0         = inf
  | p == 1         = -inf
  | p == 0.5       = mean d
  | p > 0 && p < 1 = x * ndCdfDenom d + mean d
  | otherwise      =
    error $ "Statistics.Distribution.Normal.complQuantile: p must be in [0,1] range. Got: "++show p
  where x          = invErfc (2 * p)
        inf        = 1/0

max_stars_repo_name: runeksvendsen/hs-gauge

max_issues_repo_name: runeksvendsen/hs-gauge

max_forks_repo_name: runeksvendsen/hs-gauge


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

module AI.Layer
( LayerDefinition(..)
, Layer(..)
, Connectivity
, RandomTransform
, Randomization

, showableToLayer

, createLayer
, scaleLayer
, randomizeFully
, randomizeLocally
, connectFully
, connectLocally

, randomList
, boxMuller
, normals
, uniforms
, boundedUniforms
) where

import           Data.Binary           (Binary (..), decode, encode)
import           AI.Neuron
import           Numeric.LinearAlgebra
import           System.Random

-- | The LayerDefinition type is an intermediate type initialized by the
--   library user to define the different layers of the network.
data LayerDefinition g = LayerDefinition { neuronDef   :: Neuron
                                         , neuronCount :: Int
                                         , connect     :: Connectivity
                                         , randomize   :: Randomization g
                                         }

-- | The Layer type, which stores the weight matrix, the bias matrix, and
--   a neuron type.
data Layer = Layer { weightMatrix :: Matrix Double
                   , biasVector   :: Vector Double
                   , neuron       :: Neuron
                   } deriving Show

instance Binary Layer where
  put Layer{..} = do put weightMatrix; put biasVector
  get = do weightMatrix <- get; biasVector <- get; return Layer{..}

-- | Connectivity is the type alias for a function that defines the connective
--   matrix for two layers (fully connected, convolutionally connected, etc.)
--   and takes in the number of output and input neurons
type Connectivity = Int -> Int -> Matrix Double

-- | Randomiation is the type alias for a function that defines
--   the initial random values for the weight matrix and bias vector
--   for two layers and takes in a random transformation on an infinite
--   stream of uniformly generated numbers, a source of entropy,
--   the number of output neurons, and the number of input neurons
type Randomization g = g -> RandomTransform -> Int -> Int -> (Matrix Double, Vector Double)

-- | ConvolutionalSettings is a type alias for the receptive field size,
--   stride, zero-padding, the number of filters, the number of dimensions,
--   the weight and height of the input and output fields
type ConvolutionalSettings = Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int

-- | A random transformation type alias. It is a transformation defined on an
--   infinite list of uniformly distributed random numbers, and returns a list
--   distributed on the transforming distribution.
type RandomTransform = [Double] -> [Double]

-- | The createLayer function takes in a random transformation on an infinite
--   stream of uniformly generated numbers, a source of entropy, and two
--   layer definitions, one for the previous layer and one for the next layer.
--   It returns a layer defined by the Layer type -- a weight matrix, a bias
--   vector, and a neuron type.
createLayer :: (RandomGen g)
  => RandomTransform -> g -> LayerDefinition g -> LayerDefinition g -> Layer
createLayer t g layerDef layerDef' =
  Layer (randomMatrix * connectivity i j)
        (randomVector * bias)
        (neuronDef layerDef)
  where (randomMatrix, randomVector) = randomize layerDef' g t i j
        i = neuronCount layerDef'
        j = neuronCount layerDef
        connectivity = connect layerDef'
        bias = i |> repeat 1 -- bias connectivity (full)

scaleLayer :: Double -> Layer -> Layer
scaleLayer factor l =
  Layer (factor `scale` weightMatrix l) (factor `scale` biasVector l) (neuron l)

-- | The randomizeFully function takes in a source of entropy, the number of output
--   neurons, and the number of input neurons, and returns a tuple of the
--   a fully random matrix, and a fully random vector
randomizeFully :: (RandomGen g) => Randomization g
randomizeFully g t i j = (randomMatrix, randomVector)
  where randomMatrix = (i >< j) (randomList t g')
        randomVector = i |> randomList t g''
        (g', g'') = split g

-- | The randomizeLocally function takes in ConvolutionalSettings a source
--   of entropy, a random transform, the number of output and input neurons.
--   It returns a tuple with a random matrix and a random vector
randomizeLocally :: (RandomGen g) => Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Randomization g
randomizeLocally f s p k d w1 h1 w2 h2 g t i j = (randomMatrix, randomVector)
  where randomMatrix = fromLists (locallyRandomList f s p k d w1 h1 w2 h2 g' t i j) :: Matrix Double
        randomVector = i |> randomList t g''
        (g', g'') = split g

-- | The locallyRandomList function takes in ConvolutionalSettings a source
--   of entropy, a random transform, the number of output and input neurons.
--   It returns a list of lists of locally random numbers
locallyRandomList :: (RandomGen g) => Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> g -> RandomTransform -> Int -> Int -> [[Double]]
locallyRandomList f s p k d w1 h1 w2 h2 g t i j =
  if k == 0 then []
  else filterValues ++ nextFilterValues
    where filterValues = [replicate (rowZeroOffset + colZeroOffset) 0
                          ++ take (j - rowZeroOffset - colZeroOffset) (randomList t g')
                          | n <- [0..div i k-1],
                          let rowSize = w1 + 2 * p,
                          let postsynPerFilter = rem n (div i k),
                          let rowZeroOffset = (1 + s) * quot postsynPerFilter w2 * rowSize,
                          let colZeroOffset = (1 + s) * mod n w2]
          nextFilterValues = locallyRandomList f s p (k - 1) k w1 h1 w2 h2 g'' t (i - div i k) j
          (g', g'') = split g

-- | The connectFully function takes the number of input neurons for a layer, i,
--   and the number of output neurons of a layer, j, and returns an i x j
--   connectivity matrix for a fully connected network.
connectFully :: Connectivity
connectFully i j = (i >< j) (repeat 1)

-- | The connectLocally function takes in ConvolutionalSettings and the number
--   of output and input neurons
connectLocally :: Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Connectivity
connectLocally f s p k d w1 h1 w2 h2 i j =
  repmat (fromLists conn :: Matrix Double) k d
  where conn = [replicate rowZeroOffset 0
                ++ take (f * rowSize) (cycle fieldArea)
                ++ replicate (rowSize * colSize - rowSize * f - rowZeroOffset) 0
                | n <- [0.. quot i k-1],
                  let rowSize = w1 + 2 * p,
                  let colSize = h1 + 2 * p,
                  let rowZeroOffset = (1 + s) * quot n w2 * rowSize,
                  let fieldAreaZeroOffset = (1 + s) * mod n w2,
                  let fieldArea = replicate fieldAreaZeroOffset 0
                                   ++ replicate f 1
                                   ++ replicate (rowSize - f - fieldAreaZeroOffset) 0]

-- | To go from a showable to a layer, we also need a neuron type,
--   which is an unfortunate restriction owed to Haskell's inability to
--   serialize functions.
showableToLayer :: (RandomGen g) => (Layer, LayerDefinition g) -> Layer
showableToLayer (s, d) = Layer (weightMatrix s) (biasVector s) (neuronDef d)

-- | Initialize an infinite random list given a random transform and a source
--   of entroy.
randomList :: RandomGen g => RandomTransform -> g -> [Double]
randomList transform = transform . randoms

-- | Define a transformation on the uniform distribution to generate
--   normally distributed numbers in Haskell (the Box-Muller transform)
boxMuller :: Double -> Double -> (Double, Double)
boxMuller x1 x2 = (z1, z2)
  where z1 = sqrt ((-2) * log x1) * cos (2 * pi * x2)
        z2 = sqrt ((-2) * log x1) * sin (2 * pi * x2)

-- | This is a function of type RandomTransform that transforms a list of
--   uniformly distributed numbers to a list of normally distributed numbers.
normals :: RandomTransform
normals (x1:x2:xs) = z1:z2:normals xs
  where (z1, z2) = boxMuller x1 x2
normals _ = []

-- | A non-transformation to return a list of uniformly distributed numbers
--   from a list of uniformly distributed numbers. It's really a matter of
--   naming consistency. It generates numbers on the range (0, 1]
uniforms :: RandomTransform
uniforms xs = xs

-- | An affine transformation to return a list of uniforms on the range
--   (a, b]
boundedUniforms :: (Double, Double) -> [Double] -> [Double]
boundedUniforms (lower, upper) = map affine
  where affine x = lower + x * (upper - lower)

max_stars_repo_name: jbarrow/LambdaNet

max_issues_repo_name: jbarrow/LambdaNet

max_forks_repo_name: jbarrow/LambdaNet


text:
-- | Example originally inspired by
-- <https://ro-che.info/articles/2015-12-05-testing-fft Roman Cheplyaka>.

{-# LANGUAGE GADTs #-}
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns #-}

module Main where

import Data.Complex
import qualified Language.R as R
import Language.R (R)
import Language.R.QQ

-- Call R's FFT
r_fft :: [Complex Double] -> R s [Complex Double]
r_fft nums = do
    R.dynSEXP <$> [r| fft(nums_hs) |]

main :: IO ()
main = R.withEmbeddedR R.defaultConfig $ do
    result <- R.runRegion $ r_fft [1,2,1]
    print result

max_stars_repo_name: reactormonk/HaskellR

max_issues_repo_name: reactormonk/HaskellR

max_forks_repo_name: reactormonk/HaskellR