text:
{-# LANGUAGE CPP              #-}
{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns -fno-warn-missing-signatures #-}
{-# LANGUAGE DataKinds        #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE RankNTypes       #-}
{-# LANGUAGE TypeFamilies     #-}
{-# LANGUAGE TypeOperators    #-}
{-# LANGUAGE ViewPatterns     #-}

-----------------------------------------------------------------------------
{- |
Module      :  Numeric.LinearAlgebra.Tests
Copyright   :  (c) Alberto Ruiz 2007-14
License     :  BSD3
Maintainer  :  Alberto Ruiz
Stability   :  provisional

Some tests.

-}

module Test.Numeric.LinearAlgebra.Tests(
--  module Numeric.LinearAlgebra.Tests.Instances,
--  module Numeric.LinearAlgebra.Tests.Properties,
   qCheck,
   utest,
   runTests,
   runBenchmarks
 , binaryTests
-- , findNaN
--, runBigTests
) where

import           Data.List                                   (foldl1')
import           Numeric.LinearAlgebra
import           Numeric.LinearAlgebra.Devel
import           Numeric.LinearAlgebra.Static                (L)
import           System.Info
import           Test.HUnit                                  hiding (State, Testable, test,
                                                              (~:))
import           Test.Numeric.LinearAlgebra.Tests.Instances
import           Test.Numeric.LinearAlgebra.Tests.Properties
#if MIN_VERSION_base(4,11,0)
import           Prelude                                     hiding ((<>), (^))
#else
import           Prelude                                     hiding ((^))
#endif
import           Control.Applicative
import           Control.Arrow                               ((***))
import           Control.DeepSeq                             (NFData (..))
import           Control.Monad                               (when)
import           Control.Monad                               (ap)
import           Debug.Trace
import           Numeric.LinearAlgebra.Devel                 (unsafeFromForeignPtr,
                                                              unsafeToForeignPtr)
import qualified Prelude
import           System.CPUTime
import           System.Exit
import           Text.Printf

import           Test.QuickCheck                             (Arbitrary, Property,
                                                              Testable, arbitrary, choose,
                                                              classify, coarbitrary,
                                                              maxSize,
                                                              quickCheckWithResult, shrink,
                                                              sized, stdArgs, vector)
import qualified Test.QuickCheck                             as T

import           Test.QuickCheck.Test                        (isSuccess)

--eps = peps :: Float
--i = 0:+1 :: Complex Float

qCheck n x = do
    r <- quickCheckWithResult stdArgs {maxSize = n} x
    when (not $ isSuccess r) (exitFailure)

a ^ b = a Prelude.^ (b :: Int)

utest str b = TestCase $ assertBool str b

feye n = flipud (ident n) :: Matrix Float

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

detTest1 = det m == 26
        && det mc == 38 :+ (-3)
        && det (feye 2) == -1
    where
        m = (3><3)
            [ 1, 2, 3
            , 4, 5, 7
            , 2, 8, 4 :: Float
            ]
        mc = (3><3)
            [ 1, 2, 3
            , 4, 5, 7
            , 2, 8, iC
            ]

detTest2 = inv1 |~| inv2 && [det1] ~~ [det2]
  where
    m = complex (feye 6)
    inv1 = inv m
    det1 = det m
    (inv2,(lda,sa)) = invlndet m
    det2 = sa * exp lda

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

nd1 = (3><3) [ 1/2, 1/4, 1/4
             , 0/1, 1/2, 1/4
             , 1/2, 1/4, 1/2 :: Float]

nd2 = (2><2) [1, 0, 1, 1:: Complex Float]

expmTest1 = expm nd1 :~14~: (3><3)
 [ 1.762110887278176
 , 0.478085470590435
 , 0.478085470590435
 , 0.104719410945666
 , 1.709751181805343
 , 0.425725765117601
 , 0.851451530235203
 , 0.530445176063267
 , 1.814470592751009 ]

expmTest2 = expm nd2 :~15~: (2><2)
 [ 2.718281828459045
 , 0.000000000000000
 , 2.718281828459045
 , 2.718281828459045 ]

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

mbCholTest = utest "mbCholTest" (ok1 && ok2) where
    m1 = (2><2) [2,5,5,8 :: Float]
    m2 = (2><2) [3,5,5,9 :: Complex Float]
    ok1 = mbChol (trustSym m1) == Nothing
    ok2 = mbChol (trustSym m2) == Just (chol $ trustSym m2)

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

triTest = utest "triTest" ok1 where

  a :: Matrix R
  a = (4><4)
    [
       4.30,  0.00,  0.00, 0.00,
      -3.96, -4.87,  0.00, 0.00,
       0.40,  0.31, -8.02, 0.00,
      -0.27,  0.07, -5.95, 0.12
    ]

  w :: Matrix R
  w = (4><2)
    [
      -12.90, -21.50,
       16.75,  14.93,
      -17.55,   6.33,
      -11.04,   8.09
    ]

  v :: Matrix R
  v = triSolve Lower a w

  e :: Matrix R
  e = (4><2)
    [
      -3.0000, -5.0000,
      -1.0000,  1.0000,
       2.0000, -1.0000,
       1.0000,  6.0000
    ]

  ok1 = (norm_Inf . flatten $ e - v) <= 1e-13

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

triDiagTest = utest "triDiagTest" (ok1 && ok2) where

  dL, d, dU :: Vector Float
  dL =  fromList [3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0]
  d  =  fromList [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  dU =  fromList [4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0]

  b :: Matrix R
  b = (9><3)
    [
      1.0,   1.0,   1.0,
      1.0,  -1.0,   2.0,
      1.0,   1.0,   3.0,
      1.0,  -1.0,   4.0,
      1.0,   1.0,   5.0,
      1.0,  -1.0,   6.0,
      1.0,   1.0,   7.0,
      1.0,  -1.0,   8.0,
      1.0,   1.0,   9.0
    ]

  y :: Matrix R
  y = (9><9)
    [
      1.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      3.0, 1.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      0.0, 3.0, 1.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 3.0, 1.0, 4.0, 0.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 3.0, 1.0, 4.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 3.0, 1.0, 4.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 1.0, 4.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 1.0, 4.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 1.0
    ]

  x :: Matrix R
  x = triDiagSolve dL d dU b

  z :: Matrix C
  z = (4><4)
    [
      1.0 :+ 1.0, 4.0 :+ 4.0, 0.0 :+ 0.0, 0.0 :+ 0.0,
      3.0 :+ 3.0, 1.0 :+ 1.0, 4.0 :+ 4.0, 0.0 :+ 0.0,
      0.0 :+ 0.0, 3.0 :+ 3.0, 1.0 :+ 1.0, 4.0 :+ 4.0,
      0.0 :+ 0.0, 0.0 :+ 0.0, 3.0 :+ 3.0, 1.0 :+ 1.0
    ]

  zDL, zD, zDu :: Vector C
  zDL = fromList [3.0 :+ 3.0, 3.0 :+ 3.0, 3.0 :+ 3.0]
  zD  = fromList [1.0 :+ 1.0, 1.0 :+ 1.0, 1.0 :+ 1.0, 1.0 :+ 1.0]
  zDu = fromList [4.0 :+ 4.0, 4.0 :+ 4.0, 4.0 :+ 4.0]

  zB :: Matrix C
  zB = (4><3)
    [
      1.0 :+ 1.0,   1.0  :+   1.0,  1.0 :+ (-1.0),
      1.0 :+ 1.0, (-1.0) :+ (-1.0), 1.0 :+ (-1.0),
      1.0 :+ 1.0,   1.0  :+   1.0,  1.0 :+ (-1.0),
      1.0 :+ 1.0, (-1.0) :+ (-1.0), 1.0 :+ (-1.0)
    ]

  u :: Matrix C
  u = triDiagSolve zDL zD zDu zB

  ok1 = (maximum $ map abs $ concat $ toLists $ b - (y <> x)) <= 1e-15
  ok2 = (maximum $ map magnitude $ concat $ toLists $ zB - (z <> u)) <= 1e-15

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

triDiagRegression = utest "triDiagRegression" ok where
  minusOnes, twos :: Vector R
  minusOnes = fromList [-1, -1]
  twos      = fromList [2, 2, 2]
  k :: Matrix R
  k = (3><3)
    [  2, -1,  0
    , -1,  2, -1
    ,  0, -1,  2
    ]

  b :: Matrix R
  b = (3><1) [10, 10, 10]

  tridiag = triDiagSolve minusOnes twos minusOnes b
  simple = linearSolve k b

  ok = case simple of
    Just m  -> tridiag |~| m
    Nothing -> False

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

randomTestGaussian = (unSym c) :~3~: unSym (snd (meanCov dat))
  where
    a = (3><3) [1,2,3,
                2,4,0,
               -2,2,1]
    m = 3 |> [1,2,3]
    c = mTm a
    dat = gaussianSample 7 (10^6) m c

randomTestUniform = c :~2~: unSym (snd (meanCov dat))
  where
    c = diag $ 3 |> map ((/12).(^2)) [1,2,3]
    dat = uniformSample 7 (10^6) [(0,1),(1,3),(3,6)]

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

rot :: Float -> Matrix Float
rot a = (3><3) [ c,0,s
               , 0,1,0
               ,-s,0,c ]
    where c = cos a
          s = sin a

rotTest = fun (10^5) :~11~: rot 5E4
    where fun n = foldl1' (<>) (map rot angles)
              where angles = toList $ linspace n (0,1)

---------------------------------------------------------------------
-- vector <= 0.6.0.2 bug discovered by Patrick Perry
-- http://trac.haskell.org/vector/ticket/31

offsetTest = y == y' where
    x = fromList [0..3 :: Float]
    y = subVector 1 3 x
    (f,o,n) = unsafeToForeignPtr y
    y' = unsafeFromForeignPtr f o n

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

normsVTest = TestList [
    utest "normv2CD" $ norm2PropC v
--  , utest "normv2CF" $ norm2PropC (single v)
#ifndef NONORMVTEST
  , utest "normv2D"  $ norm2PropR x
--  , utest "normv2F"  $ norm2PropR (single x)
#endif
  , utest "normv1CD" $ norm_1 v          == 8
--  , utest "normv1CF" $ norm_1 (single v) == 8
  , utest "normv1D"  $ norm_1 x          == 6
--  , utest "normv1F"  $ norm_1 (single x) == 6

  , utest "normvInfCD" $ norm_Inf v          == 5
--  , utest "normvInfCF" $ norm_Inf (single v) == 5
  , utest "normvInfD"  $ norm_Inf x          == 3
--  , utest "normvInfF"  $ norm_Inf (single x) == 3

 ] where v = fromList [1,-2,3:+4] :: Vector (Complex Float)
         x = fromList [1,2,-3] :: Vector Float
#ifndef NONORMVTEST
         norm2PropR a = norm_2 a =~= sqrt (udot a a)
#endif
         norm2PropC a = norm_2 a =~= realPart (sqrt (a `dot` a))
         a =~= b = fromList [a] |~| fromList [b]

normsMTest = TestList [
    utest "norm2mCD" $ norm_2 v          =~= 8.86164970498005
--  , utest "norm2mCF" $ norm_2 (single v) =~= 8.86164970498005
  , utest "norm2mD"  $ norm_2 x          =~= 5.96667765076216
--  , utest "norm2mF"  $ norm_2 (single x) =~= 5.96667765076216

  , utest "norm1mCD" $ norm_1 v          == 9
--  , utest "norm1mCF" $ norm_1 (single v) == 9
  , utest "norm1mD"  $ norm_1 x          == 7
--  , utest "norm1mF"  $ norm_1 (single x) == 7

  , utest "normmInfCD" $ norm_Inf v          == 12
--  , utest "normmInfCF" $ norm_Inf (single v) == 12
  , utest "normmInfD"  $ norm_Inf x          == 8
--  , utest "normmInfF"  $ norm_Inf (single x) == 8

  , utest "normmFroCD" $ norm_Frob v          =~= 8.88819441731559
--  , utest "normmFroCF" $ norm_Frob (single v) =~~= 8.88819441731559
  , utest "normmFroD"  $ norm_Frob x          =~= 6.24499799839840
--  , utest "normmFroF"  $ norm_Frob (single x) =~~= 6.24499799839840

 ] where v = (2><2) [1,-2*iC,3:+4,7] :: Matrix (Complex Float)
         x = (2><2) [1,2,-3,5] :: Matrix Float
         a =~= b = fromList [a] :~10~: fromList [b]
--       a =~~= b = fromList [a] :~5~: fromList [b]

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

sumprodTest = TestList [
    utest "sumCD" $ sumElements z            == 6
  , utest "sumCF" $ sumElements (single z)   == 6
  , utest "sumD"  $ sumElements v            == 6
  , utest "sumF"  $ sumElements (single v)   == 6

  , utest "prodCD" $ prodProp z
  , utest "prodCF" $ prodProp (single z)
  , utest "prodD"  $ prodProp v
  , utest "prodF"  $ prodProp (single v)
 ] where v = fromList [1,2,3] :: Vector Float
         z = fromList [1,2-iC,3+iC]
         prodProp x = prodElements x == product (toList x)

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

chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where
    ms = [ diag (fromList [1,2,3 :: Float])
         , konst 3 (3,5)
         , (5><10) [1 .. ]
         , konst 5 (10,2)
         ]

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

conjuTest m = cmap conjugate (flatten (conj (tr m))) == flatten (tr m)

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

newtype State s a = State { runState :: s -> (a,s) }

instance Functor (State s)
  where
    fmap f x = pure f <*> x

instance Applicative (State s)
  where
    pure = return
    (<*>) = ap

instance Monad (State s) where
    return a = State $ \s -> (a,s)
    m >>= f = State $ \s -> let (a,s') = runState m s
                            in runState (f a) s'

state_get :: State s s
state_get = State $ \s -> (s,s)

state_put :: s -> State s ()
state_put s = State $ \_ -> ((),s)

evalState :: State s a -> s -> a
evalState m s = let (a,s') = runState m s
                in seq s' a

newtype MaybeT m a = MaybeT { runMaybeT :: m (Maybe a) }

instance Monad m => Functor (MaybeT m)
  where
    fmap f x = pure f <*> x

instance Monad m => Applicative (MaybeT m)
  where
    pure = return
    (<*>) = ap

instance Monad m => Monad (MaybeT m) where
    return a = MaybeT $ return $ Just a
    m >>= f  = MaybeT $ do
                        res <- runMaybeT m
                        case res of
                                 Nothing -> return Nothing
                                 Just r  -> runMaybeT (f r)
    fail _   = MaybeT $ return Nothing

lift_maybe m = MaybeT $ do
                        res <- m
                        return $ Just res

-- apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs
--successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool
successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (size v - 1) v))) (v ! 0)
   where stp e  = do
                  ep <- lift_maybe $ state_get
                  if t e ep
                     then lift_maybe $ state_put e
                     else (fail "successive_ test failed")

-- operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input
--successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b
successive f v = evalState (mapVectorM stp (subVector 1 (size v - 1) v)) (v ! 0)
   where stp  e = do
                  ep <- state_get
                  state_put e
                  return $ f ep e

succTest = utest "successive" $
       successive_ (>) (fromList [1 :: Float,2,3,4]) == True
    && successive_ (>) (fromList [1 :: Float,3,2,4]) == False
    && successive (+) (fromList [1..10 :: Float]) == 9 |> [3,5,7,9,11,13,15,17,19]

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

findAssocTest = utest "findAssoc" ok
  where
    ok = m1 == m2
    m1 = assoc (6,6) 7 $ zip (find (>0) (ident 5 :: Matrix Float)) [10 ..] :: Matrix Float
    m2 = diagRect 7 (fromList[10..14]) 6 6

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

condTest = utest "cond" ok
  where
    ok = step v * v == cond v 0 0 0 v
    v = fromList [-7 .. 7 ] :: Vector Float

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

conformTest = utest "conform" ok
  where
    ok = 1 + row [1,2,3] + col [10,20,30,40] + (4><3) [1..]
         == (4><3) [13,15,17
                   ,26,28,30
                   ,39,41,43
                   ,52,54,56]

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

accumTest = utest "accum" ok
  where
    x = ident 3 :: Matrix Float
    ok = accum x (+) [((1,2),7), ((2,2),3)]
         == (3><3) [1,0,0
                   ,0,1,7
                   ,0,0,4]
         &&
         toList (flatten x) == [1,0,0,0,1,0,0,0,1]

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

convolutionTest = utest "convolution" ok
  where
--    a = fromList [1..10]               :: Vector Float
    b = fromList [1..3]                :: Vector Float
    c = (5><7) [1..]                   :: Matrix Float
--    d = (3><3) [0,-1,0,-1,4,-1,0,-1,0] :: Matrix Float
    ok =  separable (corr b) c == corr2 (outer b b) c
       && separable (conv b) c == conv2 (outer b b) c

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

sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix))

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

staticTest = utest "static" (fst $ checkT (undefined :: L 3 5))

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

intTest = utest "int ops" (fst $ checkT (undefined :: Matrix I))

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

modularTest = utest "modular ops" (fst $ checkT (undefined :: Matrix (Mod 13 I)))

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

indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3
  where
    l = map g (toList (f x))
    a1 = maximum l
    b1 = minimum l
    a2 = x `atIndex` maxIndex x
    b2 = x `atIndex` minIndex x
    a3 = maxElement x
    b3 = minElement x

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

_sliceTest = TestList
    [ testSlice (chol . trustSym)  (gen 5 :: Matrix R)
    , testSlice (chol . trustSym)  (gen 5 :: Matrix C)
    , testSlice qr    (rec :: Matrix R)
    , testSlice qr    (rec :: Matrix C)
    , testSlice hess  (agen 5 :: Matrix R)
    , testSlice hess  (agen 5 :: Matrix C)
    , testSlice schur (agen 5 :: Matrix R)
    , testSlice schur (agen 5 :: Matrix C)
    , testSlice lu    (agen 5 :: Matrix R)
    , testSlice lu    (agen 5 :: Matrix C)
    , testSlice (luSolve (luPacked (agen 5 :: Matrix R))) (agen 5)
    , testSlice (luSolve (luPacked (agen 5 :: Matrix C))) (agen 5)
    , test_lus (agen 5 :: Matrix R)
    , test_lus (agen 5 :: Matrix C)

    , testSlice eig   (agen 5 :: Matrix R)
    , testSlice eig   (agen 5 :: Matrix C)
    , testSlice (eigSH . trustSym) (gen 5 :: Matrix R)
    , testSlice (eigSH . trustSym) (gen 5 :: Matrix C)
    , testSlice eigenvalues   (agen 5 :: Matrix R)
    , testSlice eigenvalues   (agen 5 :: Matrix C)
    , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix R)
    , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix C)

    , testSlice svd           (rec :: Matrix R)
    , testSlice thinSVD       (rec :: Matrix R)
    , testSlice compactSVD     (rec :: Matrix R)
    , testSlice leftSV        (rec :: Matrix R)
    , testSlice rightSV       (rec :: Matrix R)
    , testSlice singularValues (rec :: Matrix R)

    , testSlice svd           (rec :: Matrix C)
    , testSlice thinSVD       (rec :: Matrix C)
    , testSlice compactSVD     (rec :: Matrix C)
    , testSlice leftSV        (rec :: Matrix C)
    , testSlice rightSV       (rec :: Matrix C)
    , testSlice singularValues (rec :: Matrix C)

    , testSlice (linearSolve (agen 5:: Matrix R)) (agen 5)
    , testSlice (flip linearSolve (agen 5:: Matrix R)) (agen 5)

    , testSlice (linearSolve (agen 5:: Matrix C)) (agen 5)
    , testSlice (flip linearSolve (agen 5:: Matrix C)) (agen 5)

    , testSlice (linearSolveLS (ogen 5:: Matrix R)) (ogen 5)
    , testSlice (flip linearSolveLS (ogen 5:: Matrix R)) (ogen 5)

    , testSlice (linearSolveLS (ogen 5:: Matrix C)) (ogen 5)
    , testSlice (flip linearSolveLS (ogen 5:: Matrix C)) (ogen 5)

    , testSlice (linearSolveSVD (ogen 5:: Matrix R)) (ogen 5)
    , testSlice (flip linearSolveSVD (ogen 5:: Matrix R)) (ogen 5)

    , testSlice (linearSolveSVD (ogen 5:: Matrix C)) (ogen 5)
    , testSlice (flip linearSolveSVD (ogen 5:: Matrix C)) (ogen 5)

    , testSlice (linearSolveLS (ugen 5:: Matrix R)) (ugen 5)
    , testSlice (flip linearSolveLS (ugen 5:: Matrix R)) (ugen 5)

    , testSlice (linearSolveLS (ugen 5:: Matrix C)) (ugen 5)
    , testSlice (flip linearSolveLS (ugen 5:: Matrix C)) (ugen 5)

    , testSlice (linearSolveSVD (ugen 5:: Matrix R)) (ugen 5)
    , testSlice (flip linearSolveSVD (ugen 5:: Matrix R)) (ugen 5)

    , testSlice (linearSolveSVD (ugen 5:: Matrix C)) (ugen 5)
    , testSlice (flip linearSolveSVD (ugen 5:: Matrix C)) (ugen 5)

    , testSlice ((<>) (ogen 5:: Matrix R)) (gen 5)
    , testSlice (flip (<>) (gen 5:: Matrix R)) (ogen 5)
    , testSlice ((<>) (ogen 5:: Matrix C)) (gen 5)
    , testSlice (flip (<>) (gen 5:: Matrix C)) (ogen 5)
    , testSlice ((<>) (ogen 5:: Matrix Float)) (gen 5)
    , testSlice (flip (<>) (gen 5:: Matrix Float)) (ogen 5)
    , testSlice ((<>) (ogen 5:: Matrix (Complex Float))) (gen 5)
    , testSlice (flip (<>) (gen 5:: Matrix (Complex Float))) (ogen 5)
    , testSlice ((<>) (ogen 5:: Matrix I)) (gen 5)
    , testSlice (flip (<>) (gen 5:: Matrix I)) (ogen 5)
    , testSlice ((<>) (ogen 5:: Matrix Z)) (gen 5)
    , testSlice (flip (<>) (gen 5:: Matrix Z)) (ogen 5)

    , testSlice ((<>) (ogen 5:: Matrix (I ./. 7))) (gen 5)
    , testSlice (flip (<>) (gen 5:: Matrix (I ./. 7))) (ogen 5)
    , testSlice ((<>) (ogen 5:: Matrix (Z ./. 7))) (gen 5)
    , testSlice (flip (<>) (gen 5:: Matrix (Z ./. 7))) (ogen 5)

    , testSlice (flip cholSolve (agen 5:: Matrix R)) (chol $ trustSym $ gen 5)
    , testSlice (flip cholSolve (agen 5:: Matrix C)) (chol $ trustSym $ gen 5)
    , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix R)) (agen 5)
    , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix C)) (agen 5)

    , ok_qrgr        (rec :: Matrix R)
    , ok_qrgr        (rec :: Matrix C)
    , testSlice (test_qrgr 4 tau1) qrr1
    , testSlice (test_qrgr 4 tau2) qrr2
    ]
  where
    QR qrr1 tau1 = qrRaw (rec :: Matrix R)
    QR qrr2 tau2 = qrRaw (rec :: Matrix C)

    test_qrgr n t x = qrgr n (QR x t)

    ok_qrgr x = TestCase . assertBool "ok_qrgr" $ simeq 1E-15 q q'
      where
        (q,_) = qr x
        atau = qrRaw x
        q' = qrgr (rows q) atau

    simeq eps a b =  not $ magnit eps (norm_1 $ flatten (a-b))

    test_lus m = testSlice f lup
      where
        f x = luSolve (LU x p) m
        (LU lup p) = luPacked m

    gen :: Numeric t => Int -> Matrix t
    gen n = diagRect 1 (konst 5 n) n n

    agen :: (Numeric t, Num (Vector t))=> Int -> Matrix t
    agen n = gen n + fromInt ((n><n)[0..])

    ogen :: (Numeric t, Num (Vector t))=> Int -> Matrix t
    ogen n = gen n === gen n

    ugen :: (Numeric t, Num (Vector t))=> Int -> Matrix t
    ugen n = takeRows 3 (gen n)

    rec :: Numeric t => Matrix t
    rec = subMatrix (0,0) (4,5) (gen 5)

    testSlice f x@(size->sz@(r,c)) =
      TestList . map (TestCase . assertEqual "" (f x)) $ (map f (g y1 ++ g y2))
      where
        subm = subMatrix
        g y = [ subm (a*r,b*c) sz y | a <-[0..2], b <- [0..2]]
        h z = fromBlocks (replicate 3 (replicate 3 z))
        y1  = h x
        y2  = (tr . h . tr) x

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

-- | All tests must pass with a maximum dimension of about 20
--  (some tests may fail with bigger sizes due to precision loss).
runTests :: Int  -- ^ maximum dimension
         -> IO ()
runTests n = do
    let test :: forall t . T.Testable t => t -> IO ()
        test p = qCheck n p
    putStrLn "------ index"
    test( \m -> indexProp id flatten (single (m :: RM)) )
    test( \v -> indexProp id id (single (v :: Vector Float)) )
    test( \m -> indexProp id flatten (m :: RM) )
    test( \v -> indexProp id id (v :: Vector Float) )
    test( \m -> indexProp magnitude flatten (single (m :: CM)) )
    test( \v -> indexProp magnitude id (single (v :: Vector (Complex Float))) )
    test( \m -> indexProp magnitude flatten (m :: CM) )
    test( \v -> indexProp magnitude id (v :: Vector (Complex Float)) )
    putStrLn "------ mult Float"
    test (multProp1 10 . rConsist)
    test (multProp1 10 . cConsist)
    test (multProp2 10 . rConsist)
    test (multProp2 10 . cConsist)
--    putStrLn "------ mult Float"
--    test (multProp1  6 . (single *** single) . rConsist)
--    test (multProp1  6 . (single *** single) . cConsist)
--    test (multProp2  6 . (single *** single) . rConsist)
--    test (multProp2  6 . (single *** single) . cConsist)
    putStrLn "------ sub-trans"
    test (subProp . rM)
    test (subProp . cM)
    putStrLn "------ ctrans"
    test (conjuTest . cM)
    test (conjuTest . zM)
    putStrLn "------ lu"
    test (luProp    . rM)
    test (luProp    . cM)
    putStrLn "------ inv (linearSolve)"
    test (invProp   . rSqWC)
    test (invProp   . cSqWC)
    putStrLn "------ luSolve"
    test (linearSolveProp (luSolve.luPacked) . rSqWC)
    test (linearSolveProp (luSolve.luPacked) . cSqWC)
    putStrLn "------ ldlSolve"
    test (linearSolvePropH (ldlSolve.ldlPacked) . rSymWC)
    test (linearSolvePropH (ldlSolve.ldlPacked) . cSymWC)
    putStrLn "------ cholSolve"
    test (linearSolveProp (cholSolve.chol.trustSym) . rPosDef)
    test (linearSolveProp (cholSolve.chol.trustSym) . cPosDef)
    putStrLn "------ luSolveLS"
    test (linearSolveProp linearSolveLS . rSqWC)
    test (linearSolveProp linearSolveLS . cSqWC)
    test (linearSolveProp2 linearSolveLS . rConsist)
    test (linearSolveProp2 linearSolveLS . cConsist)
    putStrLn "------ pinv (linearSolveSVD)"
    test (pinvProp  . rM)
    test (pinvProp  . cM)
    putStrLn "------ det"
    test (detProp   . rSqWC)
    -- test (detProp   . cSqWC)
    putStrLn "------ svd"
    test (svdProp1  . rM)
    test (svdProp1  . cM)
    test (svdProp1a svd . rM)
    test (svdProp1a svd . cM)
--    test (svdProp1a svdRd)
    test (svdProp1b svd . rM)
    test (svdProp1b svd . cM)
--    test (svdProp1b svdRd)
    test (svdProp2 thinSVD . rM)
    test (svdProp2 thinSVD . cM)
--    test (svdProp2 thinSVDRd)
--    test (svdProp2 thinSVDCd)
    test (svdProp3  . rM)
    test (svdProp3  . cM)
    test (svdProp4  . rM)
    test (svdProp4  . cM)
    test (svdProp5a)
    test (svdProp5b)
    test (svdProp6a)
    test (svdProp6b)
    test (svdProp7  . rM)
    test (svdProp7  . cM)
--    putStrLn "------ svdCd"
#ifdef NOZGESDD
--    putStrLn "Omitted"
#else
--    test (svdProp1a svdCd)
--    test (svdProp1b svdCd)
#endif
    putStrLn "------ eig"
    test (eigSHProp . rHer)
    test (eigSHProp . cHer)
    test (eigProp   . rSq)
    test (eigProp   . cSq)
    test (eigSHProp2 . rHer)
    test (eigSHProp2 . cHer)
    test (eigProp2   . rSq)
    test (eigProp2   . cSq)
    putStrLn "------ geig"
    test (uncurry geigProp . rSq2WC)
    test (uncurry geigProp . cSq2WC)
    putStrLn "------ nullSpace"
    test (nullspaceProp . rM)
    test (nullspaceProp . cM)
    putStrLn "------ qr"
    test (qrProp     . rM)
    test (qrProp     . cM)
    test (rqProp     . rM)
--    test (rqProp     . cM)
    test (rqProp1     . cM)
    test (rqProp2     . cM)
--    test (rqProp3     . cM)
    putStrLn "------ hess"
    test (hessProp   . rSq)
    test (hessProp   . cSq)
    putStrLn "------ schur"
    test (schurProp2 . rSq)
    test (schurProp1 . cSq)
    putStrLn "------ chol"
    test (cholProp   . rPosDef)
    test (cholProp   . cPosDef)
--    test (exactProp  . rPosDef)
--    test (exactProp  . cPosDef)
    putStrLn "------ expm"
    test (expmDiagProp . complex. rSqWC)
    test (expmDiagProp . cSqWC)
    putStrLn "------ vector operations - Float"
    test (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM))
    test $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::CM)) . liftMatrix makeUnitary
    test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))
    test (\u -> cos u * tan u |~| sin (u::RM))
    test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary
--    putStrLn "------ vector operations - Float"
--    test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM))
--    test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary
--    test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM))
--    test (\u -> cos u * tan u |~~| sin (u::FM))
--    test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary
    putStrLn "------ read . show"
    test (\m -> (m::RM) == read (show m))
    test (\m -> (m::CM) == read (show m))
    test (\m -> toRows (m::RM) == read (show (toRows m)))
    test (\m -> toRows (m::CM) == read (show (toRows m)))
    test (\m -> (m::FM) == read (show m))
    test (\m -> (m::ZM) == read (show m))
    test (\m -> toRows (m::FM) == read (show (toRows m)))
    test (\m -> toRows (m::ZM) == read (show (toRows m)))
    putStrLn "------ some unit tests"
    c <- runTestTT $ TestList
        [ utest "1E5 rots" rotTest
        , utest "det1" detTest1
        , utest "invlndet" detTest2
        , utest "expm1" (expmTest1)
        , utest "expm2" (expmTest2)
        , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)
        , utest "arith2" $ ((scalar (1+iC) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*iC-51) :: CM)
        , utest "arith3" $ exp (scalar iC * ones(10,10)*pi) + 1 |~| 0
        , utest "<\\>"   $ (3><2) [2,0,0,3,1,1::Float] <\> 3|>[4,9,5] |~| 2|>[2,3]
--        , utest "gamma" (gamma 5 == 24.0)
--        , besselTest
--        , exponentialTest
        , utest "randomGaussian" randomTestGaussian
        , utest "randomUniform" randomTestUniform
        , utest "buildVector/Matrix" $
                        complex (10 |> [0::Float ..]) == build 10 id
                     && ident 5 == build (5,5) (\r c -> if r==c then 1::Float else 0)
        , utest "rank" $  rank ((2><3)[1,0,0,1,5*peps,0::Float]) == 1
                       && rank ((2><3)[1,0,0,1,7*peps,0::Float]) == 2
        , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
        , mbCholTest
        , triTest
        , triDiagTest
        , triDiagRegression
        , utest "offset" offsetTest
        , normsVTest
        , normsMTest
        , sumprodTest
        , chainTest
        , succTest
        , findAssocTest
        , condTest
        , conformTest
        , accumTest
        , convolutionTest
        , sparseTest
        , staticTest
        , intTest
        , modularTest
        -- , sliceTest
        ]
    when (errors c + failures c > 0) exitFailure
    return ()

-- single precision approximate equality
-- infixl 4 |~~|
-- a |~~| b = a :~6~: b

makeUnitary v | realPart n > 1    = v / scalar n
              | otherwise = v
    where n = sqrt (v `dot` v)

binaryTests :: IO ()
binaryTests = do
  let test :: forall t . T.Testable t => t -> IO ()
      test = qCheck 100
  test vectorBinaryRoundtripProp
  test staticVectorBinaryRoundtripProp
  qCheck 30 matrixBinaryRoundtripProp
  qCheck 30 staticMatrixBinaryRoundtripProp

-- -- | Some additional tests on big matrices. They take a few minutes.
-- runBigTests :: IO ()
-- runBigTests = undefined

{-
-- | testcase for nonempty fpu stack
findNaN :: Int -> Bool
findNaN n = all (bugProp . eye) (take n $ cycle [1..20])
  where eye m = ident m :: Matrix ( Float)
-}

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

-- | Performance measurements.
runBenchmarks :: IO ()
runBenchmarks = do
    solveBench
    subBench
    mkVecBench
    multBench
    cholBench
    luBench
    luBench_2
    svdBench
    eigBench
    putStrLn ""

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

time msg act = do
    putStr (msg++" ")
    t0 <- getCPUTime
    act `seq` putStr " "
    t1 <- getCPUTime
    printf "%6.2f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Float)) :: IO ()
    return ()

timeR msg act = do
    putStr (msg++" ")
    t0 <- getCPUTime
    putStr (show act)
    t1 <- getCPUTime
    printf "%6.2f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Float)) :: IO ()
    return ()

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

manymult n = foldl1' (<>) (map rot2 angles) where
    angles = toList $ linspace n (0,1)
    rot2 :: Float -> Matrix Float
    rot2 a = (3><3) [ c,0,s
                    , 0,1,0
                    ,-s,0,c ]
        where c = cos a
              s = sin a

multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Float))

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

manyvec0 xs = sum $ map (\x -> x + x**2 + x**3) xs
manyvec1 xs = sumElements $ fromRows $ map (\x -> fromList [x,x**2,x**3]) xs
manyvec5 xs = sumElements $ fromRows $ map (\x -> vec3 x (x**2) (x**3)) xs

manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs
manyvec3 xs = sum $ map (norm_2 . (\x -> fromList [x,x**2,x**3])) xs

manyvec4 xs = sum $ map (norm_2 . (\x -> vec3 x (x**2) (x**3))) xs

vec3 :: Float -> Float -> Float -> Vector Float
vec3 a b c = runSTVector $ do
    v <- newUndefinedVector 3
    writeVector v 0 a
    writeVector v 1 b
    writeVector v 2 c
    return v

mkVecBench = do
    let n = 1000000
        xs = toList $ linspace n (0,1::Float)
    putStr "\neval data... "; print (sum xs)
    timeR "listproc        " $ manyvec0 xs
    timeR "fromList matrix " $ manyvec1 xs
    timeR "vec3 matrix     " $ manyvec5 xs
    timeR "listproc norm   " $ manyvec2 xs
    timeR "norm fromList   " $ manyvec3 xs
    timeR "norm vec3       " $ manyvec4 xs

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

subBench = do
    putStrLn ""
    let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (size v -1) v))
    time "0.1M subVector   " (g (konst 1 (1+10^5) :: Vector Float) ! 0)
    let f = foldl1' (.) (replicate (10^5) (fromRows.toRows))
    time "subVector-join  3" (f (ident  3 :: Matrix Float) `atIndex` (0,0))
    time "subVector-join 10" (f (ident 10 :: Matrix Float) `atIndex` (0,0))

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

multBench = do
    let a = ident 1000 :: Matrix Float
    let b = ident 2000 :: Matrix Float
    a `seq` b `seq` putStrLn ""
    time "product of 1M different 3x3 matrices" (manymult (10^6))
    putStrLn ""
    time "product of 1M constant  1x1 matrices" (multb 1)
    time "product of 1M constant  3x3 matrices" (multb 3)
    --time "product of 1M constant  5x5 matrices" (multb 5)
    time "product of 1M const.  10x10 matrices" (multb 10)
    --time "product of 1M const.  15x15 matrices" (multb 15)
    time "product of 1M const.  20x20 matrices" (multb 20)
    --time "product of 1M const.  25x25 matrices" (multb 25)
    putStrLn ""
    time "product (1000 x 1000)<>(1000 x 1000)" (a<>a)
    time "product (2000 x 2000)<>(2000 x 2000)" (b<>b)

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

eigBench = do
    let m = reshape 1000 (randomVector 777 Uniform (1000*1000))
        s = m + tr m
    m `seq` s `seq` putStrLn ""
    time "eigenvalues  symmetric 1000x1000" (eigenvaluesSH (trustSym m))
    time "eigenvectors symmetric 1000x1000" (snd $ eigSH (trustSym m))
    time "eigenvalues  general   1000x1000" (eigenvalues m)
    time "eigenvectors general   1000x1000" (snd $ eig m)

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

svdBench = do
    let a = reshape 500  (randomVector 777 Uniform (3000*500))
        b = reshape 1000 (randomVector 777 Uniform (1000*1000))
        fv (_,_,v) = v `atIndex` (0,0)
    a `seq` b `seq` putStrLn ""
    time "singular values  3000x500" (singularValues a)
    time "thin svd         3000x500" (fv $ thinSVD a)
    time "full svd         3000x500" (fv $ svd a)
    time "singular values 1000x1000" (singularValues b)
    time "full svd        1000x1000" (fv $ svd b)

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

solveBenchN n = do
    let x = uniformSample 777 (2*n) (replicate n (-1,1))
        a = tr x <> x
        b = asColumn $ randomVector 666 Uniform n
    a `seq` b `seq` putStrLn ""
    time ("svd solve " ++ show n) (linearSolveSVD a b)
    time (" ls solve " ++ show n) (linearSolveLS a b)
    time ("    solve " ++ show n) (linearSolve a b)
--    time (" LU solve " ++ show n) (luSolve (luPacked a) b)
    time ("LDL solve " ++ show n) (ldlSolve (ldlPacked (trustSym a)) b)
    time ("cholSolve " ++ show n) (cholSolve (chol $ trustSym a) b)

solveBench = do
    solveBenchN 500
    solveBenchN 1000
    solveBenchN 1500

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

cholBenchN n = do
    let x = uniformSample 777 (2*n) (replicate n (-1,1))
        a = tr x <> x
    a `seq` putStr ""
    time ("chol " ++ show n) (chol $ trustSym a)

cholBench = do
    putStrLn ""
    cholBenchN 1200
    cholBenchN 600
    cholBenchN 300
--    cholBenchN 150
--    cholBenchN 50

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

luBenchN f n x msg = do
    let m = diagRect 1 (fromList (replicate n x)) n n
    m `seq` putStr ""
    time (msg ++ " "++ show n) (rnf $ f m)

luBench = do
    putStrLn ""
    luBenchN luPacked  1000 (5::R)          "luPacked  Float    "
    luBenchN luPacked' 1000 (5::R)          "luPacked' Float    "
    luBenchN luPacked' 1000 (5::Mod 9973 I) "luPacked' I mod 9973"
    luBenchN luPacked' 1000 (5::Mod 9973 Z) "luPacked' Z mod 9973"

luBenchN_2 f g n x msg = do
    let m = diagRect 1 (fromList (replicate n x)) n n
        b = flipud m
    m `seq` b `seq` putStr ""
    time (msg ++ " "++ show n) (f (g m) b)

luBench_2 = do
    putStrLn ""
    luBenchN_2 luSolve  luPacked  500 (5::R)          "luSolve .luPacked  Float    "
    luBenchN_2 luSolve' luPacked' 500 (5::R)          "luSolve'.luPacked' Float    "
    luBenchN_2 luSolve' luPacked' 500 (5::Mod 9973 I) "luSolve'.luPacked' I mod 9973"
    luBenchN_2 luSolve' luPacked' 500 (5::Mod 9973 Z) "luSolve'.luPacked' Z mod 9973"

max_stars_repo_name: schnecki/hmatrix-float

max_issues_repo_name: schnecki/hmatrix-float

max_forks_repo_name: schnecki/hmatrix-float


text:
{-# LANGUAGE DataKinds #-}

module SCG.Graph () where
{--
import GHC.TypeLits
import Algebra.Graph
import Backprop
import Numeric.LinearAlgebra.Static

{--
Optimizing loss functions over random variables is often intractable
due to the loss functions and their gradients being either a sum over
an exponential number of latent variable configurations or non-analytic,
high-dimensional integrals.

Usually this is resolved by using problem-specific MC Grad Estimators.

There is a general technique for this!

Stochastic Computation Graphs:
- Allows derivation of unbiased gradient estimators for general expected losses.
- Estimator can be computed as the grad of a differentiable "surrogate loss" through backprop
- Variance reduction techniques can be applied to this problem formulation
- Hessian-free methods and majorization-minimization algorithms can be generalized to the SCG framework

!! The main modification to backprop is to introduce extra gradient signals at the stochastic nodes.
--}

-- Gradient Estimators for a Single Random Variable

type RV = Vector Double

type Theta = Matrix

parameterizedProbDist :: Theta -> RV

costfn :: RV -> Double
costfn rv = 2.0

-- scoreFunctionEstimator: By Log-Derivative Trick (valid if p_x_theta is continuous function of theta, though not necessarily for x)
-- also known as REINFORCE or likelihood ratio estimator
ddTheta_expect_over_x_f_x x f p theta = expectation $ x $ f x $ ddTheta $ log $ p x theta

-- if x is deterministic differentiable (perhaps representable as a type constraint? find one)
-- function of theta and another rv z, (i.e, x(z, theta)) we can use
-- !!Pathwise Derivative
-- only valid if f(x(z, theta)) is continous function of theta for all z (another type constraint?)
ddTheta_expect_over_x_f_ztheta f x z theta = expectation $ z $ ddTheta $ f (x z theta)

-- Theta might appear inside the expectation and the prob dist!
-- ddtheta Expectation_z~p(., theta)[ f( x(z, theta) ) ]
-- then two terms in the gradient estimator:

ddtheta_expect_over_z_from_pOfTheta_f_of_x_of_z_theta p f x z theta =
  expectation map (\z -> pointwiseEstimate) sample (p z theta)
  where pointwiseEstimate = ddTheta $ f (x z theta) + (ddTheta $ log (p z theta)) * (f $ x z theta)

{--
1. SF can be used even if f is discontinuous or x is a discrete rv
2. SF only requires trajectories, PD requires f'(x)
3. SF has higher variance than PD, unless f is rough as in time-series problems with exploding gradients
4. PD has a deterministic limit, SF doesn't.
--}

{--
  STOCHASTIC COMPUTATION GRAPHS:
    Directed, Acyclic Graph with Three Kinds of nodes:
    - Input Nodes, set externally including the parameters we differentiate with respect to
    - Deterministic Nodes: pure functions of their parents
    - Stochastic Nodes: Distributed Conditionally on their parents.

    Each parent v of a non-input node w is connected to it by a directed edge (v, w)

  THE STRUCTURE OF THE GRAPH FULLY SPECIFIES:
    - What estimator we will use, SF or PD or a combination thereof
    - Nodes arranged in series are multiplicative terms only
    - Nodes arranged in parallel lead to sums over mulplicative terms
    -
--}
{--
data Node a = InputNode a | DeterministicNode a | StochasticNode a

newtype InputNode a = IO a

newtype DeterministicNode v e = ([v] -> e -> b)

newtype StochasticNode v e = ConditionalDistribution [v] e

data DirectedEdge a = DirectedEdge
  { _v :: Node,
    _w :: DeterministicNode | StochasticNode
  }
--}

-- TODO: Create a directed acyclic graph DS using alga
-- TODO: Define Node types which typecheck (initially for a scalar parameterization of the relevant types)
-- TODO: Add backprop for folds over the deterministic nodes? derive an instance of Backprop for edges?
-- TODO: Add the right gradient estimation procedure for combinations of deterministic and stochastic nodes
-- TODO: Figure out how to do dataflow programming over the Graph
--}

ext: hs

max_stars_repo_name: faezs/opt-expect

max_issues_repo_name: faezs/opt-expect

max_forks_repo_name: faezs/opt-expect


text:
module Ch01.Broadcast where

import Numeric.LinearAlgebra

as = matrix 2 [1,2,3,4]
bs = matrix 2 [10,20]
r = as * bs

xs :: Matrix Int
xs = (3><2) [51,55,14,19,0,4]


text:
-- * This module defines how to run and extract time/space data
--   from benchmarks.
--
module Test.BigOh.Benchmark
  ( -- * Running benchmarks
    runInputs
  , getTimes
  , getStdDevs
  , defaultBenchmarkConfig
  ) where

import           Control.Monad.IO.Class
import           Criterion.Internal
import           Criterion.Main
import           Criterion.Measurement
import           Criterion.Monad
import           Criterion.Types
import           Statistics.Resampling.Bootstrap

import           Test.BigOh.Generate

defaultBenchmarkConfig = defaultConfig { timeLimit = 1, resamples = 100 }

runOne :: Config -> Benchmarkable -> IO Report
runOne cfg x
  = withConfig cfg
  $ do liftIO initializeTime
       runAndAnalyseOne 0 "" x

runInputs :: Config -> [(Benchmarkable, Input a)] -> IO [(Input a, Report)]
runInputs cfg xs
  = zip (map snd xs) <$> mapM (runOne cfg . fst) xs

getTimes :: [(Input a, Report)] -> [(Int, Double)]
getTimes = map go
  where
    go (i, report)
      = (inputSize i, estPoint $ anMean $ reportAnalysis report)

getStdDevs :: [(Input a, Report)] -> [Double]
getStdDevs = map (estPoint . anStdDev . reportAnalysis . snd)

max_stars_repo_name: tranma/shitty-complexity

max_issues_repo_name: tranma/shitty-complexity

max_forks_repo_name: tranma/shitty-complexity


text:
{-# LANGUAGE AllowAmbiguousTypes   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE TypeOperators         #-}
module Frames.Visualization.VegaLite.ParameterPlots
  (
    -- * Types for named parameters
    ParameterEstimate(..)
  , NamedParameterEstimate(..)
    -- * Plots
  , parameterPlot
  , parameterPlotMany
  , parameterPlotFlex

    -- * Configuration Re-exports
  , ViewConfig(..)
  , AxisBounds(..)
  , TimeEncoding(..)

    -- * Re-exports
  , DateTime(..)
  )
where

import qualified Frames.Visualization.VegaLite.Data
                                               as D

import           Graphics.Vega.VegaLite.Configuration
                                                ( ViewConfig(..)
                                                , configuredVegaLite
                                                , AxisBounds(..)
                                                , TimeEncoding(..)
                                                )
import qualified Graphics.Vega.VegaLite.Compat as VegaCompat

import           Data.Text                      ( Text )
import qualified Data.Text                     as T
import qualified Graphics.Vega.VegaLite        as GV
import           Graphics.Vega.VegaLite         ( DateTime(..) )

import           Text.Printf                    ( printf )

import qualified Frames                        as F
import qualified Frames.Melt                   as F
import qualified Data.Vinyl.TypeLevel          as V
import qualified Data.Vinyl                    as V

import qualified Frames.Transform              as FT

import qualified Statistics.Types              as S

-- | A type to represent the details of a parameter estimate
data ParameterEstimate = ParameterEstimate  { value :: Double, confidence :: (Double, Double) } deriving (Show)
data NamedParameterEstimate = NamedParameterEstimate { name :: Text, pEstimate :: ParameterEstimate } deriving (Show)

type Name = "Name" F.:-> T.Text
type Estimate = "Estimate" F.:-> Double
type Lo = "Lo" F.:-> Double
type Hi = "Hi" F.:-> Double

unpackParameterEstimate
  :: forall pe rs
   . (F.ElemOf rs pe, V.KnownField pe, V.Snd pe ~ NamedParameterEstimate)
  => F.Record rs
  -> F.Record (rs V.++ '[Name, Estimate, Lo, Hi])
unpackParameterEstimate row = FT.mutate (f . F.rgetField @pe) row
 where
  f :: NamedParameterEstimate -> F.Record '[Name, Estimate, Lo, Hi]
  f (NamedParameterEstimate n (ParameterEstimate v (lo, hi))) =
    n F.&: v F.&: lo F.&: hi F.&: V.RNil

-- | Plot parameters with error bars
--  Flex version handles a foldable of results so we can, e.g.,
--  1. Compare across time or other variable in the data
--  2. Compare different fitting methods for the same data

-- TODO: Fix replicated y-axis ticks.  There since we need the difference between them (trailing "'"s) to provide y-offset
-- to the different results.  Otherwise they would overlap in y.  Maybe fix by switching to layers and only having Y-axis labels
-- on 0th?

type NullText = "NullText" F.:-> T.Text

parameterPlot
  :: forall pe rs f
   . ( Foldable f
     , Functor f
     , V.KnownField pe
     , V.KnownField NullText
     , F.ElemOf (rs V.++ '[NullText]) pe
     , F.ElemOf (rs V.++ '[NullText]) NullText
     , F.ElemOf (rs V.++ '[NullText] V.++ '[Name, Estimate, Lo, Hi]) NullText
     , F.ElemOf (rs V.++ '[NullText] V.++ '[Name, Estimate, Lo, Hi]) Estimate
     , F.ElemOf (rs V.++ '[NullText] V.++ '[Name, Estimate, Lo, Hi]) Lo
     , F.ElemOf (rs V.++ '[NullText] V.++ '[Name, Estimate, Lo, Hi]) Hi
     , F.ElemOf (rs V.++ '[NullText] V.++ '[Name, Estimate, Lo, Hi]) Name
     , V.Snd pe ~ NamedParameterEstimate
     )
  => T.Text
  -> S.CL Double
  -> ViewConfig
  -> f (D.Row rs)
  -> GV.VegaLite
parameterPlot title cl vc rows = parameterPlotFlex @NullText @pe
  False
  title
  cl
  vc
  (labeledRows rows)
 where
  nullText = ("" F.&: V.RNil) :: F.Record '[NullText]
  labeledRows x = fmap (FT.mutate (const nullText)) x

parameterPlotMany
  :: forall l pe rs f
   . ( Foldable f
     , Functor f
     , D.DataFieldOf rs l
     , F.ElemOf rs pe
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) l
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) Estimate
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) Lo
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) Hi
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) Name
     , V.KnownField pe
     , V.Snd l ~ T.Text
     , V.Snd pe ~ NamedParameterEstimate
     )
  => T.Text
  -> S.CL Double
  -> ViewConfig
  -> f (D.Row rs)
  -> GV.VegaLite
parameterPlotMany = parameterPlotFlex @l @pe True

parameterPlotFlex
  :: forall l pe rs f
   . ( Foldable f
     , Functor f
     , D.DataFieldOf rs l
     , F.ElemOf rs pe
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) l
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) Estimate
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) Lo
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) Hi
     , F.ElemOf (rs V.++ '[Name, Estimate, Lo, Hi]) Name
     , V.KnownField pe
     , V.Snd l ~ T.Text
     , V.Snd pe ~ NamedParameterEstimate
     )
  => Bool -- ^ legend?
  -> T.Text -- ^ title
  -> S.CL Double -- ^ confidence level
  -> ViewConfig -- ^ size, etc.
  -> f (D.Row rs) -- ^ data
  -> GV.VegaLite
parameterPlotFlex haveLegend title cl vc rows
  = let
      unpackedRows = fmap (unpackParameterEstimate @pe @rs) rows
      dat          = D.recordsToVLData (F.rcast @'[Name, l, Lo, Hi, Estimate])
                                       D.defaultParse
                                       unpackedRows
      xLabel =
        "Estimate (with "
          <> (T.pack $ printf "%2.0f" (100 * S.confidenceLevel cl))
          <> "% confidence error bars)"
      estimateXEnc = GV.position
        GV.X
        [ D.pName @Estimate
        , GV.PmType GV.Quantitative
        , GV.PAxis [GV.AxTitle xLabel]
        ]
      estimateYEnc = GV.position GV.Y [D.pName @Name, GV.PmType GV.Ordinal]
      handleLegend l = if haveLegend then l else (GV.MLegend []) : l
      estimateColorEnc =
        GV.color $ handleLegend $ [D.mName @l, GV.MmType GV.Nominal]
      estimateEnc  = estimateXEnc . estimateYEnc . estimateColorEnc
      estConfLoEnc = GV.position
        GV.X
        [D.pName @Lo, GV.PmType GV.Quantitative, GV.PAxis [GV.AxTitle ""]]
      estConfHiEnc = GV.position
        GV.X2
        [D.pName @Hi, GV.PmType GV.Quantitative, GV.PAxis [GV.AxTitle ""]]
      estConfEnc =
        estConfLoEnc . estConfHiEnc . estimateYEnc . estimateColorEnc
      estSpec  = GV.asSpec [(GV.encoding . estimateEnc) [], GV.mark GV.Point []]
      confSpec = GV.asSpec [(GV.encoding . estConfEnc) [], GV.mark GV.Rule []]
    in
      configuredVegaLite
        vc
        [VegaCompat.title title, GV.layer [estSpec, confSpec], dat]

max_stars_repo_name: adamConnerSax/Frames-hvega

max_issues_repo_name: adamConnerSax/Frames-hvega

max_forks_repo_name: adamConnerSax/Frames-hvega


text:
module Internal
  where
import Data.Complex

type Cplx = Complex Double

toCplx :: Double -> Cplx
toCplx x = x :+ 0.0

getPhiK :: Cplx -> (Cplx, Int)
getPhiK phi
  | realPart phi > pi/2 =
    until (\(x,_) -> realPart x <= pi/2) (\(x,k) -> (x-pi,k+1)) (phi,0)
  | realPart phi < -pi/2 =
    until (\(x,_) -> realPart x >= -pi/2) (\(x,k) -> (x+pi,k-1)) (phi,0)
  | otherwise = (phi,0)

atanC :: Cplx -> Cplx
atanC z = i * (log(1-i*z) - log(1+i*z)) / 2
  where
    i = 0.0 :+ 1.0

max_stars_repo_name: stla/elliptic

max_issues_repo_name: stla/elliptic

max_forks_repo_name: stla/elliptic


text:
module STFT (stft, istft) where

import Data.Foldable (toList)
import Data.Complex
import Numeric.FFT
import DSP.Window (hanning)

-- Sublists having heads at positions [1,1+noverlap..] of length framesz
split :: Int -> Int -> [a] -> [[a]]
split framesz hop = reverse . split_aux [] where
  split_aux res l
    | length newSeg < framesz = res
    | otherwise = split_aux (newSeg:res) (drop hop l)
    where  newSeg = take framesz l

hannWin :: [Complex Double] ->  [Complex Double]
hannWin l = zipWith (*) hann l
  where hann = map (:+0) $ toList $ hanning $ length l + 1

-- STFT function, calculate X(u,f)
stft :: Int -> Int -> [Complex Double] -> [[Complex Double]]
stft framesz hop = map dft . map hannWin . split framesz hop

-- Overlap function, the list must be no empty
overlap :: Num a => Int -> [[a]] -> [a]
overlap hop = foldr1 overlap_aux where
  overlap_aux seg signal = seg1 ++ (zipWith (+) seg2 signal1) ++ signal2
    where (seg1, seg2) = splitAt hop seg
          noverlap = length seg2
          (signal1, signal2) = splitAt noverlap signal

-- ISTFT function, calculate x(t)
istft :: Int -> [[Complex Double]] -> [Complex Double]
istft hop = overlap hop . map idft


text:
{-# Language DerivingVia                #-}
{-# Language FlexibleInstances          #-}
{-# Language GeneralizedNewtypeDeriving #-}
{-# Language MultiParamTypeClasses      #-}
{-# Language UndecidableInstances       #-}

module Geometry.Testing (
    DInt (..)
  , Parameter (..)
  , Position (..)
  , AnyPoint (..)
  , AnyCircle (..)
  , AnyAngle (..)
  , AnyLine (..)
  , AnyRay (..)
  , AnySegment (..)
  , Motion (..), appMotion
  , Nontrivial (..)
  , NonDegenerate (..)
  , (<==>)
) where

import Data.Complex
import Data.Foldable
import Data.Monoid
import Test.QuickCheck hiding (scale)
import Test.QuickCheck.Modifiers
import Data.Tuple.Extra (second)

import Geometry
import Geometry.Circle (mkCircle)
import Geometry.Base (roundUp)

infix 0 <==>
a <==> b = (a ==> b) .&&. (b ==> a)

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

{- |  A type for an interger double number.

>>>  sample (arbitrary :: Gen DInt)
DInt 0.0
DInt 2.0
DInt 0.0
DInt (-1.0)
DInt 7.0
...
-}
newtype DInt = DInt Double deriving Show

instance Arbitrary DInt where
  arbitrary = DInt . fromIntegral <$> (arbitrary :: Gen Int)
  shrink (DInt l) = DInt <$> shrink l

{- |  A type for a double number parameterizing the `Manifold`.

>>>  sample (arbitrary :: Gen Parameter)
Parameter 0.6977012699848539
Parameter 0.5321234247634057
Parameter 9.832672492232597e-2
Parameter 0.9510866292762408
...
-}
newtype Parameter = Parameter { getParameter :: Double} deriving Show

instance Arbitrary Parameter where
  arbitrary = Parameter <$> choose (0, 1)
  shrink (Parameter l) = Parameter <$> shrink l

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

{- |  A type for an affine position (`XY`, `Cmp`, `Point` etc).

>>>  sample (arbitrary :: Gen (Position XY))
Position (0.0,0.0)
Position (-1.0,2.0)
Position (-3.0,3.0)
Position (6.0,4.0)
Position (-1.0,-1.0
...
-}
newtype Position a = Position a
  deriving (Show, Affine, Trans, Metric, Eq)

instance (Eq a, Pnt a, Arbitrary a) =>
         Arbitrary (Position a) where
  arbitrary = Position . roundUp 1 <$> arbitrary
  shrink = shrinkPos 0.5

shrinkPos :: (Eq a, Affine a, Metric a) => Double -> a -> [a]
shrinkPos d x = res
  where res = filter (/= x) $
              map (roundUp d) $
              takeWhile (\p -> distance x p >= d/2) $
              map (\s -> asCmp $ (1 - s) * cmp x) $
              iterate (/2) 1

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

instance Arbitrary Direction where
  arbitrary = oneof [asDeg  <$> arbitrary, asCmp <$> arbitrary]
  shrink = shrinkPos 1

------------------------------------------------------------
instance Arbitrary Point where
  arbitrary = Point <$> arbitrary
  shrink = shrinkPos 1

{- |  A wrapper for a `Point`. -}
newtype AnyPoint = AnyPoint Point
  deriving (Show, Arbitrary)

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

instance Arbitrary Circle where
  arbitrary = do (Position c) <- arbitrary
                 (Position r) <- arbitrary :: Gen (Position Cmp)
                 return (mkCircle (norm r) c # rotateAt' c (angle r))

  shrink cir =
    do Position c <- shrink (Position (center cir))
       r <- shrink (radius cir)
       return $ mkCircle r c

{- |  A wrapper for nondegenative `Circle`. -}
newtype AnyCircle = AnyCircle Circle
  deriving Show
  deriving Arbitrary via Nontrivial Circle

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

instance Arbitrary Angle where
  arbitrary = do (Position p) <- arbitrary
                 Angle p <$> arbitrary <*> arbitrary

  shrink an =
    do Position p <- shrink (Position (refPoint an))
       s <- shrink (angleStart an)
       e <- shrink (angleEnd an)
       return $ Angle p s e

{- |  A wrapper for nontrivial `Angle`. -}
newtype AnyAngle = AnyAngle Angle
  deriving Show
  deriving Arbitrary via Nontrivial Angle

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

instance Arbitrary Segment where
  arbitrary =
    do Position p1 <- arbitrary
       Position p2 <- arbitrary
       return $ Segment (p1, p2)

  shrink l = let (p1, p2) = refPoints l
    in do Position p1' <- shrink (Position p1)
          Position p2' <- shrink (Position p2)
          return $ Segment (p1', p2')

instance Arbitrary Line where
  arbitrary = asLine <$> (arbitrary :: Gen Segment)
  shrink l = asLine <$> shrink (asSegment l)

instance Arbitrary Ray where
  arbitrary = asRay <$> (arbitrary :: Gen Segment)
  shrink l = asRay <$> shrink (asSegment l)

{- |  A wrapper for nontrivial `Line`. -}
newtype AnyLine = AnyLine Line
  deriving Show
  deriving Arbitrary via Nontrivial Line

{- |  A wrapper for nontrivial `Ray`. -}
newtype AnyRay = AnyRay Ray
  deriving Show
  deriving Arbitrary via Nontrivial Ray

{- |  A wrapper for nontrivial `Segment`. -}
newtype AnySegment = AnySegment Segment
  deriving Show
  deriving Arbitrary via Nontrivial Segment

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

{- |  A wrapper for a motion `translate`, `rotate` or `reflect`. -}
newtype Motion a = Motion (String, Endo a)

appMotion :: Motion a -> a -> a
appMotion (Motion (_, m)) = appEndo m

instance Show (Motion a) where
  show (Motion (s, _)) = s

instance Trans a => Arbitrary (Motion a) where
  arbitrary = Motion . fold <$> listOf (oneof motions)
    where
      motions = [ label "Tr " translate <$> (arbitrary :: Gen XY)
                , label "Rot " rotate <$> arbitrary
                , label "Ref " reflect <$> arbitrary ]
      label l t x = (l <> show x <> " ", Endo $ t x)

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

{- |  A wrapper for a nontrivial `Figure`. -}
newtype Nontrivial a = Nontrivial a deriving
  ( Eq
  , Show
  , Figure
  , Affine
  , Trans
  , Curve
  , ClosedCurve
  , Manifold
  , PiecewiseLinear
  , Polygonal
  , APoint
  , Linear
  , Circular
  , Angular
  )

instance (Arbitrary a, Figure a) => Arbitrary (Nontrivial a) where
  arbitrary = Nontrivial <$> arbitrary `suchThat` isNontrivial
  shrink (Nontrivial l) = Nontrivial <$> filter isNontrivial (shrink l)

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

{- |  A wrapper for a nondegenerate `Figure`. -}
newtype NonDegenerate a = NonDegenerate a deriving
  ( Eq
  , Show
  , Figure
  , Affine
  , Trans
  , Curve
  , ClosedCurve
  , Manifold
  , PiecewiseLinear
  , Polygonal
  , APoint
  , Linear
  , Circular
  , Angular
  )

instance (Arbitrary a, PiecewiseLinear a) => Arbitrary (NonDegenerate a) where
  arbitrary = NonDegenerate <$> arbitrary `suchThat` isNondegenerate
  shrink (NonDegenerate l) = NonDegenerate <$> filter isNondegenerate (shrink l)

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

instance Arbitrary Triangle where
  arbitrary = (`suchThat` isNontrivial) $ do
    Position p1 <- arbitrary
    Position p2 <- arbitrary
    Position p3 <- arbitrary
    pure $ Triangle [p1, p2, p3]

  shrink (Triangle [p1, p2, p3]) = filter isNontrivial $ do
    p1' <- shrink p1
    p2' <- shrink p2
    p3' <- shrink p3
    pure $ Triangle [p1', p2', p3']

instance Arbitrary RightTriangle where
  arbitrary = (`suchThat` isNondegenerate) $ do
    Positive a <- arbitrary
    Positive b <- arbitrary
    m <- arbitrary
    pure $ aRightTriangle # scaleX a # scaleY b # appMotion m
  shrink (RightTriangle t) = RightTriangle <$> shrink t

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
module Vec (module Vec, module L) where

import Foreign
import qualified Data.Vector as V
import qualified Data.Vector.Storable as S
import qualified Data.Vector.Unboxed as U
import Numeric.LinearAlgebra as L hiding ((<>))

type UVec = U.Vector
type SVec = S.Vector
type BVec = V.Vector

type Vec2 = L.Vector R
type Vec3 = L.Vector R
type Vec4 = L.Vector R

type Lin2 = (Vec2, Vec2)
type Tgl2 = (Vec2, Vec2, Vec2)

type Lin3 = (Vec3, Vec3)
type Tgl3 = (Vec3, Vec3, Vec3)

type Mat2 = L.Matrix R
type Mat3 = L.Matrix R
type Mat4 = L.Matrix R

{-# INLINE vec #-}
vec :: Storable a => [a] -> L.Vector a
vec = L.fromList

{-# INLINE lis #-}
lis :: Storable a => L.Vector a -> [a]
lis = L.toList

rx :: (Storable a, Floating a) => a -> Matrix a
rx a = let
  s = sin a
  c = cos a
  in (3><3) [ 1, 0 , 0,
              0, c, -s,
              0, s,  c ]

ry :: (Storable a, Floating a) => a -> Matrix a
ry a = let
  s = sin a
  c = cos a
  in (3><3) [  c, 0, s,
               0, 1, 0,
              -s, 0, c ]

rz :: (Storable a, Floating a) => a -> Matrix a
rz a = let
  s = sin a
  c = cos a
  in (3><3) [ c, -s, 0,
              s,  c, 0,
              0,  0, 1 ]


text:
module Fungine.Graphics.Geometry where

import Fungine.Prelude hiding (rotate)
import Numeric.LinearAlgebra

data Point2d = Point2d Float Float deriving Show
data Point3d = Point3d Float Float Float deriving Show

cross3 :: Point3d -> Point3d -> Point3d
cross3 (Point3d x y z) (Point3d x' y' z') =
  Point3d (y * z' - z * y') (z * x' - x * z') (x * y' - y * x')

magnitude2 :: Point2d -> Float
magnitude2 (Point2d x y) = sqrt (x ** 2 + y ** 2)

magnitude3 :: Point3d -> Float
magnitude3 (Point3d x y z) = sqrt (x ** 2 + y ** 2 + z ** 2)

unit2 :: Point2d -> Point2d
unit2 p@(Point2d x y) = let m = magnitude2 p in Point2d (x / m) (y / m)

unit3 :: Point3d -> Point3d
unit3 p@(Point3d x y z) = let m = magnitude3 p in Point3d (x / m) (y / m) (z / m)

data Rectangle = Rectangle Point2d Point2d deriving Show

rectangleDim :: Rectangle -> (Float, Float)
rectangleDim (Rectangle (Point2d tx ty) (Point2d bx by)) = (bx - tx, by - ty)

rectangleOffset :: Rectangle -> (Float, Float) -> Rectangle
rectangleOffset (Rectangle (Point2d tx ty) (Point2d bx by)) (xoff, yoff) =
  Rectangle (Point2d (tx + xoff) (ty + yoff)) (Point2d (bx + xoff) (by + yoff))

mkRectangle :: (Float, Float) -> Float -> Float -> Rectangle
mkRectangle (tx, ty) w h = Rectangle (Point2d tx ty) (Point2d (tx + w) (ty + h))

rectangleTopLeft :: Rectangle -> (Float, Float)
rectangleTopLeft (Rectangle (Point2d tx ty) _) = (tx, ty)

data Box = Box Point3d Point3d

translate :: Point3d -> Matrix Float
translate (Point3d x y z) = (4 >< 4) [1, 0, 0, x, 0, 1, 0, y, 0, 0, 1, z, 0, 0, 0, 1]

rotate :: Float -> Point3d -> Matrix Float
rotate ang vec =
  let
    t               = 1 - cos ang
    s               = sin ang
    c               = cos ang
    (Point3d x y z) = unit3 vec
  in (4 >< 4)
    [ t * (x ** 2) + c
    , t * x * y - s * z
    , t * x * z + s * y
    , 0
    , t * x * y + s * z
    , t * (y ** 2) + c
    , t * y * z - s * x
    , 0
    , t * x * z - s * y
    , t * y * z + s * x
    , t * (z ** 2) + c
    , 0
    , 0
    , 0
    , 0
    , 1
    ]

perspective :: Float -> Float -> Float -> Float -> Matrix Float
perspective fov aspect nearDist farDist
  | fov <= 0
  = ident 4
  | aspect <= 0
  = ident 4
  | otherwise
  = let
      uh            = 1 / tan (fov / 2)
      frustrumDepth = farDist - nearDist
    in (4 >< 4)
      [ uh / aspect
      , 0
      , 0
      , 0
      , 0
      , uh
      , 0
      , 0
      , 0
      , 0
      , farDist * (1 / frustrumDepth)
      , (negate farDist * nearDist) * (1 / frustrumDepth)
      , 0
      , 0
      , 1
      , 0
      ]

scale :: Float -> Matrix Float
scale f = (4 >< 4) [f, 0, 0, 0, 0, f, 0, 0, 0, 0, f, 0, 0, 0, 0, 1]

scale3 :: Float -> Float -> Float -> Matrix Float
scale3 x y z = (4 >< 4) [x, 0, 0, 0, 0, y, 0, 0, 0, 0, z, 0, 0, 0, 0, 1]

-- give distance, look at, pitch, yaw, roll
-- in our game, camera is on a sphere that can be zoomed in or out
-- no roll
camera :: Float -> Point3d -> Float -> Float -> Matrix Float
camera zoom (Point3d fx fy fz) pitch yaw =
  let
    y = zoom * sin pitch
    z = zoom * cos pitch
  in translate (Point3d fx (fy - y) (fz - z)) * rotate yaw (Point3d 0 1 0)

max_stars_repo_name: jmdurr/fungine

max_issues_repo_name: jmdurr/fungine

max_forks_repo_name: jmdurr/fungine


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Data.Pass.Util.Beta
-- Copyright   :  (C) 2012-2013 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  experimental
-- Portability :  non-portable (GADTs, Rank2Types)
--
----------------------------------------------------------------------------
module Data.Pass.Util.Beta
  ( Beta(..)
  , cumulative
  , density
  , quantile
  ) where

import Numeric.SpecFunctions
import Numeric.MathFunctions.Constants (m_NaN)

data Beta = Beta {-# UNPACK #-} !Double {-# UNPACK #-} !Double
  deriving (Eq,Read,Show)

cumulative :: Beta -> Double -> Double
cumulative (Beta a b) x
  | x <= 0 = 0
  | otherwise = incompleteBeta a b x
{-# INLINE cumulative #-}

-- invert a monotone function
invertMono :: (Double -> Double) -> Double -> Double -> Double -> Double
invertMono f l0 h0 b = go l0 h0 where
  go l h
    | h - l < epsilon = m
    | otherwise = case compare (f m) b of
      LT -> go m h
      EQ -> m
      GT -> go l m
    where m = l + (h-l)/2
          epsilon = 1e-12
{-# INLINE invertMono #-}

density :: Beta -> Double -> Double
density (Beta a b) x
  | a <= 0 || b <= 0 = m_NaN
  | x <= 0 = 0
  | x >= 1 = 0
  | otherwise = exp $ (a-1)*log x + (b-1)*log (1-x) - logBeta a b
{-# INLINE density #-}

quantile :: Beta -> Double -> Double
quantile d p
  | p < 0 = error $ "probability must be positive. Got " ++ show p
  | p > 1 = error $ "probability must be less than 1. Got " ++ show p
  | otherwise = invertMono (cumulative d) 0 1 p
{-# INLINE quantile #-}

max_stars_repo_name: ekmett/multipass

max_issues_repo_name: ekmett/multipass

max_forks_repo_name: ekmett/multipass