text:
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeOperators       #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}

module Test.Grenade.Layers.Convolution where

import           Data.Constraint
import           Data.Proxy
import           Data.Singletons
import           Data.Singletons.Prelude.Num            ((%*))
import           GHC.TypeLits
import           Unsafe.Coerce

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

#if MIN_VERSION_base(4,9,0)
import           Data.Kind                              (Type)
#endif

import           Hedgehog
import qualified Hedgehog.Gen                           as Gen

import           Test.Hedgehog.Compat
import           Test.Hedgehog.Hmatrix
import           Test.Hedgehog.TypeLits

import qualified Numeric.LinearAlgebra                  as LA
import qualified Numeric.LinearAlgebra.Static           as H

import           Grenade.Core
import           Grenade.Layers.Convolution
import           Grenade.Types
import           Grenade.Utils.ListStore

import qualified Test.Grenade.Layers.Internal.Reference as Reference

data OpaqueConvolution :: Type where
     OpaqueConvolution :: Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns -> OpaqueConvolution

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

data OpaqueBiasConvolution :: Type where
     OpaqueBiasConvolution :: Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns -> OpaqueBiasConvolution

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

genConvolution :: ( KnownNat channels
                  , KnownNat filters
                  , KnownNat kernelRows
                  , KnownNat kernelColumns
                  , KnownNat strideRows
                  , KnownNat strideColumns
                  , KnownNat kernelFlattened
                  , kernelFlattened ~ (kernelRows * kernelColumns * channels)
                  ) => Gen (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns)
genConvolution = Convolution <$> uniformSample <*> pure mkListStore

genBiasConvolution :: ( KnownNat channels
                      , KnownNat filters
                      , KnownNat kernelRows
                      , KnownNat kernelColumns
                      , KnownNat strideRows
                      , KnownNat strideColumns
                      , KnownNat kernelFlattened
                      , kernelFlattened ~ (kernelRows * kernelColumns * channels)
                      ) => Gen (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns)
genBiasConvolution = BiasConvolution <$> uniformSample <*> randomVector <*> pure mkListStore

genOpaqueConvolution :: Gen OpaqueConvolution
genOpaqueConvolution = do
    channels <- genNat
    filters  <- genNat
    kernel_h <- genNat
    kernel_w <- genNat
    stride_h <- genNat
    stride_w <- genNat
    case (channels, filters, kernel_h, kernel_w, stride_h, stride_w) of
       ( SomeNat (pch :: Proxy ch), SomeNat  (_   :: Proxy fl),
         SomeNat (pkr :: Proxy kr), SomeNat  (pkc :: Proxy kc),
         SomeNat (_   :: Proxy sr), SomeNat  (_   :: Proxy sc)) ->
          let p1 = singByProxy pkr
              p2 = singByProxy pkc
              p3 = singByProxy pch
          in  case p1 %* p2 %* p3 of
            SNat -> OpaqueConvolution <$> (genConvolution :: Gen (Convolution 'WithoutBias padding ch fl kr kc sr sc))

genOpaqueBiasConvolution :: Gen OpaqueBiasConvolution
genOpaqueBiasConvolution = do
    channels <- genNat
    filters  <- genNat
    kernel_h <- genNat
    kernel_w <- genNat
    stride_h <- genNat
    stride_w <- genNat
    case (channels, filters, kernel_h, kernel_w, stride_h, stride_w) of
       ( SomeNat (pch :: Proxy ch), SomeNat  (_   :: Proxy fl),
         SomeNat (pkr :: Proxy kr), SomeNat  (pkc :: Proxy kc),
         SomeNat (_   :: Proxy sr), SomeNat  (_   :: Proxy sc)) ->
          let p1 = singByProxy pkr
              p2 = singByProxy pkc
              p3 = singByProxy pch
          in  case p1 %* p2 %* p3 of
            SNat -> OpaqueBiasConvolution <$> (genBiasConvolution :: Gen (Convolution 'WithBias padding ch fl kr kc sr sc))

prop_conv_no_padding_forward_and_backward_correct :: Property
prop_conv_no_padding_forward_and_backward_correct = withTests 20 $ property $ do
  OpaqueConvolution (convLayer'@(Convolution kernels _) :: Convolution 'WithoutBias padding channels filters kernelRows kernelCols strideRows strideCols) <- forAll genOpaqueConvolution

  let convLayer = unsafeCoerce convLayer' :: Convolution 'WithoutBias 'NoPadding channels filters kernelRows kernelCols strideRows strideCols
      ok stride kernel = [extent | extent <- [(kernel + 1) .. 30 ], (extent - kernel) `mod` stride == 0]
      kr = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
      kc = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
      sr = fromIntegral $ natVal (Proxy :: Proxy strideRows)
      sc = fromIntegral $ natVal (Proxy :: Proxy strideCols)
      fs = fromIntegral $ natVal (Proxy :: Proxy filters   ) :: Int
      cs = fromIntegral $ natVal (Proxy :: Proxy channels  ) :: Int

  -- generate input size
  er <- forAll (Gen.element (ok sr kr))
  ec <- forAll (Gen.element (ok sc kc))

  -- calculate output size
  let rr = ((er - kr) `div` sr) + 1
      rc = ((ec - kc) `div` sc) + 1

  case (someNatVal er, someNatVal ec, someNatVal rr, someNatVal rc) of
      ( Just (SomeNat (pinr :: Proxy inRows)), Just (SomeNat (_  :: Proxy inCols)), Just (SomeNat (pour :: Proxy outRows)), Just (SomeNat (_ :: Proxy outCols))) ->
        case ( singByProxy pinr %* singByProxy (Proxy :: Proxy channels)
             , singByProxy pour %* singByProxy (Proxy :: Proxy filters)
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strideRows * (outRows - 1) <= (inRows - kernelRows + 0) ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inRows - kernelRows + 0) <= (outRows * strideRows ) - 1 ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strideCols * (outCols - 1) <= (inCols - kernelCols + 0) ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inCols - kernelCols + 0) <= (outCols * strideCols ) - 1 ) ) ) of
                    (SNat, SNat, Dict, Dict, Dict, Dict) -> do
                      input :: S ('D3 inRows inCols channels) <- forAll (S3D <$> uniformSample)

                      let (tape, out@(S3D out') :: S ('D3 outRows outCols filters)) = runForwards convLayer input
                          input'   = (\(S3D x) -> H.extract x) input
                          refOut   = Reference.convForwards input' cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) (fromIntegral rr) (fromIntegral rc) (fromIntegral sr) (fromIntegral sc)

                      H.extract out' `isSimilarMatrixTo` refOut

                      let (Convolution' dW, S3D dx :: S ('D3 inRows inCols channels))
                            = runBackwards convLayer tape out
                          (refdx, refdw) = Reference.convBackProp input' cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) refOut (fromIntegral rr) (fromIntegral rc) (fromIntegral sr) (fromIntegral sc)

                      H.extract dx `isSimilarMatrixTo` refdx
                      H.extract dW `isSimilarMatrixTo` refdw

prop_conv_same_upper_padding_forward_and_backward_correct :: Property
prop_conv_same_upper_padding_forward_and_backward_correct = withTests 20 $ property $ do
  OpaqueConvolution (convLayer'@(Convolution kernels _) :: Convolution 'WithoutBias padding channels filters kernelRows kernelCols strideRows strideCols) <- forAll genOpaqueConvolution

  let convLayer = unsafeCoerce convLayer' :: Convolution 'WithoutBias 'SameUpper channels filters kernelRows kernelCols strideRows strideCols
      ok stride kernel = [extent | extent <- [(kernel + 1) .. 30 ], (extent - kernel) `mod` stride == 0]
      kr = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
      kc = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
      sr = fromIntegral $ natVal (Proxy :: Proxy strideRows)
      sc = fromIntegral $ natVal (Proxy :: Proxy strideCols)
      fs = fromIntegral $ natVal (Proxy :: Proxy filters   ) :: Int
      cs = fromIntegral $ natVal (Proxy :: Proxy channels  ) :: Int

  -- generate input size
  er <- forAll (Gen.element (ok sr kr))
  ec <- forAll (Gen.element (ok sc kc))

  annotate $ show (kr, kc, sr, sc, fs, cs)

  -- calculate padding size
  let pw = fromIntegral $ (ec - 1) * sc + kc - ec
      ph = fromIntegral $ (er - 1) * sr + kr - er
      padt = div ph 2
      padl = div pw 2
      padb = ph - padt
      padr = pw - padl

  case (someNatVal er, someNatVal ec) of
      ( Just (SomeNat (pinr :: Proxy inRows)), Just (SomeNat (_  :: Proxy inCols)) ) ->
        case ( singByProxy pinr %* singByProxy (Proxy :: Proxy channels)
             , singByProxy pinr %* singByProxy (Proxy :: Proxy filters) ) of
                    (SNat, SNat) -> do
                      input@(S3D inp') :: S ('D3 inRows inCols channels) <- forAll (S3D <$> uniformSample)

                      let (tape, out@(S3D out') :: S ('D3 inRows inCols filters)) = runForwards convLayer input
                          input'   = (\(S3D x) -> H.extract x) input
                          refOut   = Reference.convForwardsWithPadding input' cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) (fromIntegral er) (fromIntegral ec) (fromIntegral sr) (fromIntegral sc) padl padt padr padb

                      H.extract out' `isSimilarMatrixTo` refOut

                      let (Convolution' dW, S3D dx :: S ('D3 inRows inCols channels))
                            = runBackwards convLayer tape out
                          (refdx, refdw)   = Reference.convBackPropWithPadding (H.extract inp') cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) refOut (fromIntegral er) (fromIntegral ec) (fromIntegral sr) (fromIntegral sc) padl padt padr padb

                      H.extract dx `isSimilarMatrixTo` refdx
                      H.extract dW `isSimilarMatrixTo` refdw

prop_conv_same_lower_padding_forward_and_backward_correct :: Property
prop_conv_same_lower_padding_forward_and_backward_correct = withTests 20 $ property $ do
  OpaqueConvolution (convLayer'@(Convolution kernels _) :: Convolution 'WithoutBias padding channels filters kernelRows kernelCols strideRows strideCols) <- forAll genOpaqueConvolution

  let convLayer = unsafeCoerce convLayer' :: Convolution 'WithoutBias 'SameLower channels filters kernelRows kernelCols strideRows strideCols
      ok stride kernel = [extent | extent <- [(kernel + 1) .. 30 ], (extent - kernel) `mod` stride == 0]
      kr = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
      kc = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
      sr = fromIntegral $ natVal (Proxy :: Proxy strideRows)
      sc = fromIntegral $ natVal (Proxy :: Proxy strideCols)
      fs = fromIntegral $ natVal (Proxy :: Proxy filters   ) :: Int
      cs = fromIntegral $ natVal (Proxy :: Proxy channels  ) :: Int

  -- generate input size
  er <- forAll (Gen.element (ok sr kr))
  ec <- forAll (Gen.element (ok sc kc))

  annotate $ show (kr, kc, sr, sc, fs, cs)

  -- calculate padding size
  let pw = fromIntegral $ (ec - 1) * sc + kc - ec
      ph = fromIntegral $ (er - 1) * sr + kr - er
      padb = div ph 2
      padr = div pw 2
      padt = ph - padb
      padl = pw - padr

  case (someNatVal er, someNatVal ec) of
      ( Just (SomeNat (pinr :: Proxy inRows)), Just (SomeNat (_  :: Proxy inCols)) ) ->
        case ( singByProxy pinr %* singByProxy (Proxy :: Proxy channels)
             , singByProxy pinr %* singByProxy (Proxy :: Proxy filters) ) of
                    (SNat, SNat) -> do
                      input@(S3D inp') :: S ('D3 inRows inCols channels) <- forAll (S3D <$> uniformSample)

                      let (tape, out@(S3D out') :: S ('D3 inRows inCols filters)) = runForwards convLayer input
                          input'   = (\(S3D x) -> H.extract x) input
                          refOut   = Reference.convForwardsWithPadding input' cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) (fromIntegral er) (fromIntegral ec) (fromIntegral sr) (fromIntegral sc) padl padt padr padb

                      H.extract out' `isSimilarMatrixTo` refOut

                      let (Convolution' dW, S3D dx :: S ('D3 inRows inCols channels))
                            = runBackwards convLayer tape out
                          (refdx, refdw)   = Reference.convBackPropWithPadding (H.extract inp') cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) refOut (fromIntegral er) (fromIntegral ec) (fromIntegral sr) (fromIntegral sc) padl padt padr padb

                      H.extract dx `isSimilarMatrixTo` refdx
                      H.extract dW `isSimilarMatrixTo` refdw

prop_conv_explicit_padding_forward_and_backward_correct :: Property
prop_conv_explicit_padding_forward_and_backward_correct = withTests 20 $ property $ do
  OpaqueConvolution (convLayer'@(Convolution kernels _) :: Convolution 'WithoutBias padding channels filters kernelRows kernelCols strideRows strideCols) <- forAll genOpaqueConvolution

  let ok stride kernel = [extent | extent <- [(kernel + 1) .. 30 ], (extent - kernel) `mod` stride == 0]
      kr = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
      kc = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
      sr = fromIntegral $ natVal (Proxy :: Proxy strideRows)
      sc = fromIntegral $ natVal (Proxy :: Proxy strideCols)
      fs = fromIntegral $ natVal (Proxy :: Proxy filters   ) :: Int
      cs = fromIntegral $ natVal (Proxy :: Proxy channels  ) :: Int

  padl <- forAll $ choose 0 kc
  padt <- forAll $ choose 0 kr
  padr <- forAll $ choose 0 kc
  padb <- forAll $ choose 0 kr

  -- generate input size
  er <- forAll (Gen.element (ok sr kr))
  ec <- forAll (Gen.element (ok sc kc))

  let rr = ((padt + er + padb - kr) `div` sr) + 1
      rc = ((padl + ec + padr - kc) `div` sc) + 1

  case (someNatVal er, someNatVal ec, someNatVal rr, someNatVal rc, someNatVal padl, someNatVal padt, someNatVal padr, someNatVal padb) of
      ( Just (SomeNat (pinr :: Proxy inRows)), Just (SomeNat (_  :: Proxy inCols)), Just (SomeNat (poutr :: Proxy outRows)), Just (SomeNat (_  :: Proxy outCols)), Just (SomeNat (_  :: Proxy padl)), Just (SomeNat (_  :: Proxy padt)), Just (SomeNat (_  :: Proxy padr)), Just (SomeNat (_  :: Proxy padb)) ) ->
        case ( singByProxy pinr %* singByProxy (Proxy :: Proxy channels)
             , singByProxy poutr %* singByProxy (Proxy :: Proxy filters)
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strideRows * (outRows - 1) <= (inRows - kernelRows + (padt + padb)) ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inRows - kernelRows + (padt + padb)) <= (outRows * strideRows ) - 1 ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strideCols * (outCols - 1) <= (inCols - kernelCols + (padl + padr)) ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inCols - kernelCols + (padl + padr)) <= (outCols * strideCols ) - 1 ) ) ) of
                    (SNat, SNat, Dict, Dict, Dict, Dict) -> do
                      input@(S3D inp') :: S ('D3 inRows inCols channels) <- forAll (S3D <$> uniformSample)

                      let convLayer = unsafeCoerce convLayer' :: Convolution 'WithoutBias ('Padding padl padt padr padb) channels filters kernelRows kernelCols strideRows strideCols

                      let (tape, out@(S3D out') :: S ('D3 outRows outCols filters)) = runForwards convLayer input
                          input'   = (\(S3D x) -> H.extract x) input
                          refOut   = Reference.convForwardsWithPadding input' cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) (fromIntegral rr) (fromIntegral rc) (fromIntegral sr) (fromIntegral sc) (fromIntegral padl) (fromIntegral padt) (fromIntegral padr) (fromIntegral padb)

                      H.extract out' `isSimilarMatrixTo` refOut

                      let (Convolution' dW, S3D dx :: S ('D3 inRows inCols channels))
                            = runBackwards convLayer tape out
                          (refdx, refdw)   = Reference.convBackPropWithPadding (H.extract inp') cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) refOut (fromIntegral rr) (fromIntegral rc) (fromIntegral sr) (fromIntegral sc) (fromIntegral padl) (fromIntegral padt) (fromIntegral padr) (fromIntegral padb)

                      H.extract dx `isSimilarMatrixTo` refdx
                      H.extract dW `isSimilarMatrixTo` refdw

prop_bias_conv_no_padding_forward_and_backward_correct :: Property
prop_bias_conv_no_padding_forward_and_backward_correct = withTests 20 $ property $ do
  OpaqueBiasConvolution (convLayer'@(BiasConvolution kernels bias _) :: Convolution 'WithBias padding channels filters kernelRows kernelCols strideRows strideCols) <- forAll genOpaqueBiasConvolution

  let convLayer = unsafeCoerce convLayer' :: Convolution 'WithBias 'NoPadding channels filters kernelRows kernelCols strideRows strideCols
      ok stride kernel = [extent | extent <- [(kernel + 1) .. 30 ], (extent - kernel) `mod` stride == 0]
      kr = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
      kc = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
      sr = fromIntegral $ natVal (Proxy :: Proxy strideRows)
      sc = fromIntegral $ natVal (Proxy :: Proxy strideCols)
      fs = fromIntegral $ natVal (Proxy :: Proxy filters   ) :: Int
      cs = fromIntegral $ natVal (Proxy :: Proxy channels  ) :: Int

  -- generate input size
  er <- forAll (Gen.element (ok sr kr))
  ec <- forAll (Gen.element (ok sc kc))

  -- calculate output size
  let rr = ((er - kr) `div` sr) + 1
      rc = ((ec - kc) `div` sc) + 1

  case (someNatVal er, someNatVal ec, someNatVal rr, someNatVal rc) of
      ( Just (SomeNat (pinr :: Proxy inRows)), Just (SomeNat (_  :: Proxy inCols)), Just (SomeNat (pour :: Proxy outRows)), Just (SomeNat (_ :: Proxy outCols))) ->
        case ( singByProxy pinr %* singByProxy (Proxy :: Proxy channels)
             , singByProxy pour %* singByProxy (Proxy :: Proxy filters)
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strideRows * (outRows - 1) <= (inRows - kernelRows + 0) ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inRows - kernelRows + 0) <= (outRows * strideRows ) - 1 ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strideCols * (outCols - 1) <= (inCols - kernelCols + 0) ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inCols - kernelCols + 0) <= (outCols * strideCols ) - 1 ) ) ) of
                    (SNat, SNat, Dict, Dict, Dict, Dict) -> do
                      input@(S3D inp') :: S ('D3 inRows inCols channels) <- forAll (S3D <$> uniformSample)

                      let (tape, out@(S3D out') :: S ('D3 outRows outCols filters)) = runForwards convLayer input
                          input'   = (\(S3D x) -> H.extract x) input
                          refOut   = Reference.biasConvForwards input' cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) (H.extract bias) (fromIntegral rr) (fromIntegral rc) (fromIntegral sr) (fromIntegral sc)

                      H.extract out' `isSimilarMatrixTo` refOut

                      let (BiasConvolution' dW db, S3D dx :: S ('D3 inRows inCols channels))
                            = runBackwards convLayer tape out
                          (refdx, refdw, refdb) = Reference.biasConvBackProp (H.extract inp') cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) refOut (fromIntegral rr) (fromIntegral rc) (fromIntegral sr) (fromIntegral sc)

                      H.extract dx `isSimilarMatrixTo` refdx
                      H.extract dW `isSimilarMatrixTo` refdw
                      H.extract db `isSimilarVectorTo` refdb

prop_bias_conv_same_upper_forward_and_backward_correct :: Property
prop_bias_conv_same_upper_forward_and_backward_correct = withTests 20 $ property $ do
  OpaqueBiasConvolution (convLayer'@(BiasConvolution kernels bias _) :: Convolution 'WithBias padding channels filters kernelRows kernelCols strideRows strideCols) <- forAll genOpaqueBiasConvolution

  let convLayer = unsafeCoerce convLayer' :: Convolution 'WithBias 'SameUpper channels filters kernelRows kernelCols strideRows strideCols
      ok stride kernel = [extent | extent <- [(kernel + 1) .. 30 ], (extent - kernel) `mod` stride == 0]
      kr = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
      kc = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
      sr = fromIntegral $ natVal (Proxy :: Proxy strideRows)
      sc = fromIntegral $ natVal (Proxy :: Proxy strideCols)
      fs = fromIntegral $ natVal (Proxy :: Proxy filters   ) :: Int
      cs = fromIntegral $ natVal (Proxy :: Proxy channels  ) :: Int

  -- generate input size
  er <- forAll (Gen.element (ok sr kr))
  ec <- forAll (Gen.element (ok sc kc))

  -- calculate padding size
  let pw = fromIntegral $ (ec - 1) * sc + kc - ec
      ph = fromIntegral $ (er - 1) * sr + kr - er
      padt = div ph 2
      padl = div pw 2
      padb = ph - padt
      padr = pw - padl

  case (someNatVal er, someNatVal ec) of
      ( Just (SomeNat (pinr :: Proxy inRows)), Just (SomeNat (_  :: Proxy inCols)) ) ->
        case ( singByProxy pinr %* singByProxy (Proxy :: Proxy channels)
             , singByProxy pinr %* singByProxy (Proxy :: Proxy filters) ) of
                    (SNat, SNat) -> do
                      input@(S3D inp') :: S ('D3 inRows inCols channels) <- forAll (S3D <$> uniformSample)

                      let (tape, out@(S3D out') :: S ('D3 inRows inCols filters)) = runForwards convLayer input
                          input'   = (\(S3D x) -> H.extract x) input
                          refOut   = Reference.biasConvForwardsWithPadding input' cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) (H.extract bias) (fromIntegral er) (fromIntegral ec) (fromIntegral sr) (fromIntegral sc) padl padt padr padb

                      H.extract out' `isSimilarMatrixTo` refOut

                      let (BiasConvolution' dW db, S3D dx :: S ('D3 inRows inCols channels))
                            = runBackwards convLayer tape out
                          (refdx, refdw, refdb) = Reference.biasConvBackPropWithPadding (H.extract inp') cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) refOut (fromIntegral er) (fromIntegral ec) (fromIntegral sr) (fromIntegral sc) padl padt padr padb

                      H.extract dx `isSimilarMatrixTo` refdx
                      H.extract dW `isSimilarMatrixTo` refdw
                      H.extract db `isSimilarVectorTo` refdb

prop_bias_conv_same_lower_forward_and_backward_correct :: Property
prop_bias_conv_same_lower_forward_and_backward_correct = withTests 20 $ property $ do
  OpaqueBiasConvolution (convLayer'@(BiasConvolution kernels bias _) :: Convolution 'WithBias padding channels filters kernelRows kernelCols strideRows strideCols) <- forAll genOpaqueBiasConvolution

  let convLayer = unsafeCoerce convLayer' :: Convolution 'WithBias 'SameLower channels filters kernelRows kernelCols strideRows strideCols
      ok stride kernel = [extent | extent <- [(kernel + 1) .. 30 ], (extent - kernel) `mod` stride == 0]
      kr = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
      kc = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
      sr = fromIntegral $ natVal (Proxy :: Proxy strideRows)
      sc = fromIntegral $ natVal (Proxy :: Proxy strideCols)
      fs = fromIntegral $ natVal (Proxy :: Proxy filters   ) :: Int
      cs = fromIntegral $ natVal (Proxy :: Proxy channels  ) :: Int

  -- generate input size
  er <- forAll (Gen.element (ok sr kr))
  ec <- forAll (Gen.element (ok sc kc))

  -- calculate padding size
  let pw = fromIntegral $ (ec - 1) * sc + kc - ec
      ph = fromIntegral $ (er - 1) * sr + kr - er
      padb = div ph 2
      padr = div pw 2
      padt = ph - padb
      padl = pw - padr

  case (someNatVal er, someNatVal ec) of
      ( Just (SomeNat (pinr :: Proxy inRows)), Just (SomeNat (_  :: Proxy inCols)) ) ->
        case ( singByProxy pinr %* singByProxy (Proxy :: Proxy channels)
             , singByProxy pinr %* singByProxy (Proxy :: Proxy filters) ) of
                    (SNat, SNat) -> do
                      input@(S3D inp') :: S ('D3 inRows inCols channels) <- forAll (S3D <$> uniformSample)

                      let (tape, out@(S3D out') :: S ('D3 inRows inCols filters)) = runForwards convLayer input
                          input'   = (\(S3D x) -> H.extract x) input
                          refOut   = Reference.biasConvForwardsWithPadding input' cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) (H.extract bias) (fromIntegral er) (fromIntegral ec) (fromIntegral sr) (fromIntegral sc) padl padt padr padb

                      H.extract out' `isSimilarMatrixTo` refOut

                      let (BiasConvolution' dW db, S3D dx :: S ('D3 inRows inCols channels))
                            = runBackwards convLayer tape out
                          (refdx, refdw, refdb) = Reference.biasConvBackPropWithPadding (H.extract inp') cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) refOut (fromIntegral er) (fromIntegral ec) (fromIntegral sr) (fromIntegral sc) padl padt padr padb

                      H.extract dx `isSimilarMatrixTo` refdx
                      H.extract dW `isSimilarMatrixTo` refdw
                      H.extract db `isSimilarVectorTo` refdb

prop_bias_conv_explicit_padding_forward_and_backward_correct :: Property
prop_bias_conv_explicit_padding_forward_and_backward_correct = withTests 20 $ property $ do
  OpaqueBiasConvolution (convLayer'@(BiasConvolution kernels bias _) :: Convolution 'WithBias padding channels filters kernelRows kernelCols strideRows strideCols) <- forAll genOpaqueBiasConvolution

  let ok stride kernel = [extent | extent <- [(kernel + 1) .. 30 ], (extent - kernel) `mod` stride == 0]
      kr = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
      kc = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
      sr = fromIntegral $ natVal (Proxy :: Proxy strideRows)
      sc = fromIntegral $ natVal (Proxy :: Proxy strideCols)
      fs = fromIntegral $ natVal (Proxy :: Proxy filters   ) :: Int
      cs = fromIntegral $ natVal (Proxy :: Proxy channels  ) :: Int

  padl <- forAll $ choose 0 kc
  padt <- forAll $ choose 0 kr
  padr <- forAll $ choose 0 kc
  padb <- forAll $ choose 0 kr

  -- generate input size
  er <- forAll (Gen.element (ok sr kr))
  ec <- forAll (Gen.element (ok sc kc))

  let rr = ((padt + er + padb - kr) `div` sr) + 1
      rc = ((padl + ec + padr - kc) `div` sc) + 1

  case (someNatVal er, someNatVal ec, someNatVal rr, someNatVal rc, someNatVal padl, someNatVal padt, someNatVal padr, someNatVal padb) of
      ( Just (SomeNat (pinr :: Proxy inRows)), Just (SomeNat (_  :: Proxy inCols)), Just (SomeNat (poutr :: Proxy outRows)), Just (SomeNat (_  :: Proxy outCols)), Just (SomeNat (_  :: Proxy padl)), Just (SomeNat (_  :: Proxy padt)), Just (SomeNat (_  :: Proxy padr)), Just (SomeNat (_  :: Proxy padb)) ) ->
        case ( singByProxy pinr %* singByProxy (Proxy :: Proxy channels)
             , singByProxy poutr %* singByProxy (Proxy :: Proxy filters)
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strideRows * (outRows - 1) <= (inRows - kernelRows + (padt + padb)) ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inRows - kernelRows + (padt + padb)) <= (outRows * strideRows ) - 1 ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strideCols * (outCols - 1) <= (inCols - kernelCols + (padl + padr)) ) )
             , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inCols - kernelCols + (padl + padr)) <= (outCols * strideCols ) - 1 ) ) ) of
                    (SNat, SNat, Dict, Dict, Dict, Dict) -> do
                      input@(S3D inp') :: S ('D3 inRows inCols channels) <- forAll (S3D <$> uniformSample)

                      let convLayer = unsafeCoerce convLayer' :: Convolution 'WithBias ('Padding padl padt padr padb) channels filters kernelRows kernelCols strideRows strideCols

                      let (tape, out@(S3D out') :: S ('D3 outRows outCols filters)) = runForwards convLayer input
                          input'   = (\(S3D x) -> H.extract x) input
                          refOut   = Reference.biasConvForwardsWithPadding input' cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) (H.extract bias) (fromIntegral rr) (fromIntegral rc) (fromIntegral sr) (fromIntegral sc) (fromIntegral padl) (fromIntegral padt) (fromIntegral padr) (fromIntegral padb)

                      H.extract out' `isSimilarMatrixTo` refOut

                      let (BiasConvolution' dW db, S3D dx :: S ('D3 inRows inCols channels))
                            = runBackwards convLayer tape out
                          (refdx, refdw, refdb)   = Reference.biasConvBackPropWithPadding (H.extract inp') cs (fromIntegral er) (fromIntegral ec) (H.extract kernels) fs (fromIntegral kr) (fromIntegral kc) refOut (fromIntegral rr) (fromIntegral rc) (fromIntegral sr) (fromIntegral sc) (fromIntegral padl) (fromIntegral padt) (fromIntegral padr) (fromIntegral padb)

                      H.extract dx `isSimilarMatrixTo` refdx
                      H.extract dW `isSimilarMatrixTo` refdw
                      H.extract db `isSimilarVectorTo` refdb

prop_bias_convolution_same_as_torch :: Property
prop_bias_convolution_same_as_torch = withTests 1 $ property $ H.extract output `isSimilarMatrixTo` refOutput
  where
    refOutput  = LA.fromLists . concat $ outputList :: LA.Matrix RealNum
    S3D output = snd $ runForwards layer input :: S ('D3 4 4 3)

    layer   = BiasConvolution filters bias mkListStore :: Convolution 'WithBias 'NoPadding 1 3 4 4 2 2
    filters = H.tr . H.fromList . concat . concat $ filtersList
    bias    = H.fromList biasList
    input   = S2D . H.fromList . concat $ inputList :: S ('D2 10 10)

    inputList :: [[RealNum]]
    inputList = [[-1.8798149824142456,  0.1028530076146126,  1.8676880598068237,
                  -0.0617060922086239, -0.8356814980506897,  0.6487482190132141,
                   0.8218213915824890, -0.7828793525695801, -0.8224498033523560,
                  -0.9689628481864929
                 ]
                ,[-0.1040988788008690,  1.4789071083068848,  1.6646575927734375,
                   0.0237357039004564, -0.5236989855766296,  0.6900013089179993,
                  -2.3320152759552002, -1.6244263648986816, -1.0726952552795410,
                  -0.4659616351127625
                 ]
                ,[ 0.0606177598237991,  1.2164361476898193, -1.0005422830581665,
                   2.0863819122314453, -1.5312120914459229,  0.4628386795520782,
                  -1.5035923719406128,  0.5512505769729614, -0.6526298522949219,
                   0.6975844502449036
                 ]
                ,[ 1.5345604419708252,  0.2354260832071304, -0.2829766273498535,
                  -1.1268494129180908,  0.1847444474697113, -0.6780525445938110,
                   0.9897040128707886,  0.1145756468176842,  0.5400218367576599,
                  -0.6138708591461182
                 ]
                ,[-1.0893152952194214, -0.0248561538755894,  0.0916035994887352,
                   0.4344462454319000, -2.2688574790954590,  0.5042927861213684,
                  -1.2440685033798218,  0.4743961393833160,  0.2804549634456635,
                   0.1743697375059128
                 ]
                ,[-0.2276446968317032, -1.5715813636779785, -0.2782900631427765,
                   0.2772432267665863,  0.5334178209304810, -0.3267174661159515,
                  -0.6893687248229980, -0.8983678221702576,  1.7771846055984497,
                   0.6515364050865173
                 ]
                ,[-0.8124530911445618,  1.0634243488311768,  0.0320023111999035,
                  -1.3274892568588257,  0.3162937462329865,  0.6602431535720825,
                   0.7445837259292603,  0.1550943553447723,  0.1773122400045395,
                   0.5073298215866089
                 ]
                ,[ 0.8331560492515564,  1.1417356729507446,  1.3021006584167480,
                  -1.2554157972335815, -0.4929245114326477,  0.7357274293899536,
                   0.4516215324401855,  0.9766013026237488,  0.4637960195541382,
                  -1.2788373231887817
                 ]
                ,[-0.1299438923597336,  0.7541506290435791, -0.6450557708740234,
                  -0.0040270495228469,  0.7661917805671692,  0.0933827608823776,
                   0.3463796675205231,  0.1497271358966827,  0.6318596005439758,
                  -1.9423550367355347
                 ]
                ,[ 0.3258689343929291,  0.5653873682022095, -0.1593302637338638,
                   0.5531405806541443, -0.4313830733299255,  0.9392206668853760,
                   1.1189296245574951,  0.9076719880104065, -0.3691771030426025,
                   2.3205711841583252
                 ]
                ]

    filtersList :: [[[RealNum]]]
    filtersList = [[[ 0.0762235820293427, -0.0693399608135223,  0.1232089996337891, -0.1839968860149384]
                   ,[-0.2333435118198395,  0.0558746755123138,  0.1272304952144623,  0.1712160110473633]
                   ,[ 0.2303934693336487,  0.0637164115905762,  0.0874992907047272, -0.2228358089923859]
                   ,[-0.2300853431224823, -0.1301925778388977,  0.1110948324203491, -0.1208976805210114]
                   ]
                  ,[[-0.1482009291648865, -0.1966925859451294, -0.0440416038036346, -0.0799964368343353]
                   ,[ 0.0987019836902618,  0.0791423320770264,  0.2103689610958099, -0.2410672008991241]
                   ,[ 0.0028940439224243,  0.0728487372398376,  0.0047851204872131,  0.0087434351444244]
                   ,[ 0.2424267232418060, -0.1152470409870148,  0.2046474218368530, -0.1154336035251617]
                   ]
                  ,[[-0.1728554666042328,  0.2191164791584015,  0.0494877994060516, -0.1877064704895020]
                   ,[-0.0716654360294342,  0.1852594912052155,  0.1268731951713562,  0.1372097134590149]
                   ,[ 0.2338261902332306,  0.0943823754787445, -0.2169064879417419, -0.1956105828285217]
                   ,[-0.0422836244106293,  0.1209872066974640,  0.1442286074161530, -0.2474469542503357]
                   ]
                  ]

    biasList :: [RealNum]
    biasList = [-0.1065079271793365,  0.1141500771045685, -0.0434862375259399]

    outputList :: [[[RealNum]]]
    outputList = [[[-0.4326016604900360, -0.5381838679313660, -0.8197864294052124, -0.3077610433101654]
                  ,[-1.4160603284835815, -0.8261992335319519, -1.2307931184768677, -0.6055138707160950]
                  ,[-0.1438199430704117, -0.8861438035964966, -0.8439695835113525,  0.2755392491817474]
                  ,[-0.4444667994976044, -0.7047327160835266,  0.4867877364158630, -0.4625016748905182]
                  ]
                 ,[[ 1.2657278776168823,  0.0472350381314754,  0.3799144029617310,  0.2348928004503250]
                  ,[ 0.1615301817655563,  0.0406027995049953,  0.5752753019332886,  0.8285622000694275]
                  ,[ 0.5149132013320923,  0.4453703761100769,  0.2957187891006470,  0.5147784948348999]
                  ,[ 0.8481521010398865, -0.2448751479387283, -0.3152585625648499,  0.2738023698329926]
                  ]
                 ,[[ 1.0434051752090454, -0.3550303280353546,  0.0058417855761945, -0.6423963904380798]
                  ,[-1.1366062164306641,  0.8248907923698425, -0.0979268252849579, -0.2169524729251862]
                  ,[ 0.5547238588333130, -0.8525719046592712, -0.1231328472495079,  0.9879407286643982]
                  ,[ 0.8017534613609314, -1.2922899723052979,  0.5693746805191040, -0.4196736514568329]
                  ]
                 ]

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

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
-----------------------------------------------------------
-- |
-- module:                      Math.Complex
-- copyright:                   (c) 2017 HE, Tao
-- license:                     MIT
-- maintainer:                  sighingnow@gmail.com
--
-- Coordinate `Data.Complex` with foundation's numeric and primitive type support.
--
{-# OPTIONS_GHC -Wno-orphans #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UnboxedTuples #-}

module Math.Complex
  ( Complex (..)
  , conjugate
  ) where

import Foundation
import Foundation.Class.Storable
import Foundation.Primitive

import Data.Complex (Complex(..), conjugate)
import GHC.Exts

instance Additive a => Additive (Complex a) where
    azero = azero :+ azero
    (a :+ b) + (c :+ d) = (a + c) :+ (b + d)

instance (Subtractive a, Difference a ~ a) => Subtractive (Complex a) where
    type Difference (Complex a) = (Complex a)
    (a :+ b) - (c :+ d) = (a - c) :+ (b - d)

instance (Additive a, Subtractive a, Difference a ~ a, Multiplicative a) => Multiplicative (Complex a) where
    midentity = midentity :+ azero
    (a :+ b) * (c :+ d) = (a * c - b * d) :+ (b * c + a * d)

instance (Additive a, Subtractive a, Difference a ~ a, Divisible a) => Divisible (Complex a) where
    (a :+ b) / (c :+ d) = ((a * c + b * d) / (c * c + d * d)) :+ ((b * c - a * d) / (c * c + d * d))

offsetComplex :: Offset (Complex a) -> (# Int#, Int# #)
offsetComplex !(Offset (I# i)) = (# n, n +# 1# #)
    where !n = uncheckedIShiftL# i 1#

{-# INLINE offsetComplex #-}

instance PrimType (Complex Float) where
    primSizeInBytes _ = primSizeInBytes (Proxy :: Proxy Float) + primSizeInBytes (Proxy :: Proxy Float)
    {-# INLINE primSizeInBytes #-}

    primShiftToBytes _ = 3 -- TODO may be wrong
    {-# INLINE primShiftToBytes #-}

    primBaUIndex ba n = F# (indexFloatArray# ba n1) :+ F# (indexFloatArray# ba n2)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primBaUIndex #-}

    primMbaURead mba n = primitive $ \s1 -> let !(# s2, r1 #) = readFloatArray# mba n1 s1
                                                !(# s3, r2 #) = readFloatArray# mba n2 s2
                                             in (# s3, F# r1 :+ F# r2 #)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primMbaURead #-}

    primMbaUWrite mba n ((F# w1) :+ (F# w2)) = primitive $ \s1 -> let !s2 = writeFloatArray# mba n1 w1 s1
                                                                   in (# writeFloatArray# mba n2 w2 s2, () #)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primMbaUWrite #-}

    primAddrIndex addr n = F# (indexFloatOffAddr# addr n1) :+ F# (indexFloatOffAddr# addr n2)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primAddrIndex #-}

    primAddrRead addr n = primitive $ \s1 -> let !(# s2, r1 #) = readFloatOffAddr# addr n1 s1
                                                 !(# s3, r2 #) = readFloatOffAddr# addr n2 s2
                                              in (# s3, F# r1 :+ F# r2 #)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primAddrRead #-}

    primAddrWrite addr n ((F# w1) :+ (F# w2)) = primitive $ \s1 -> let !s2 = writeFloatOffAddr# addr n1 w1 s1
                                                                    in (# writeFloatOffAddr# addr n2 w2 s2, () #)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primAddrWrite #-}

instance PrimType (Complex Double) where
    primSizeInBytes _ = primSizeInBytes (Proxy :: Proxy Double) + primSizeInBytes (Proxy :: Proxy Double)
    {-# INLINE primSizeInBytes #-}

    primShiftToBytes _ = 5 -- TODO may be wrong
    {-# INLINE primShiftToBytes #-}

    primBaUIndex ba n = D# (indexDoubleArray# ba n1) :+ D# (indexDoubleArray# ba n2)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primBaUIndex #-}

    primMbaURead mba n = primitive $ \s1 -> let !(# s2, r1 #) = readDoubleArray# mba n1 s1
                                                !(# s3, r2 #) = readDoubleArray# mba n2 s2
                                             in (# s3, D# r1 :+ D# r2 #)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primMbaURead #-}

    primMbaUWrite mba n ((D# w1) :+ (D# w2)) = primitive $ \s1 -> let !s2 = writeDoubleArray# mba n1 w1 s1
                                                                   in (# writeDoubleArray# mba n2 w2 s2, () #)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primMbaUWrite #-}

    primAddrIndex addr n = D# (indexDoubleOffAddr# addr n1) :+ D# (indexDoubleOffAddr# addr n2)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primAddrIndex #-}

    primAddrRead addr n = primitive $ \s1 -> let !(# s2, r1 #) = readDoubleOffAddr# addr n1 s1
                                                 !(# s3, r2 #) = readDoubleOffAddr# addr n2 s2
                                              in (# s3, D# r1 :+ D# r2 #)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primAddrRead #-}

    primAddrWrite addr n ((D# w1) :+ (D# w2)) = primitive $ \s1 -> let !s2 = writeDoubleOffAddr# addr n1 w1 s1
                                                                    in (# writeDoubleOffAddr# addr n2 w2 s2, () #)
        where !(# n1, n2 #) = offsetComplex n
    {-# INLINE primAddrWrite #-}

instance Storable (Complex Float) where
    peek (Ptr addr) = primAddrRead addr (Offset 0)
    poke (Ptr addr) = primAddrWrite addr (Offset 0)

instance Storable (Complex Double) where
    peek (Ptr addr) = primAddrRead addr (Offset 0)
    poke (Ptr addr) = primAddrWrite addr (Offset 0)

instance StorableFixed (Complex Float) where
    size _ = size (Proxy :: Proxy Float) + size (Proxy :: Proxy Float)
    alignment _ = alignment (Proxy :: Proxy Float)

instance StorableFixed (Complex Double) where
    size _ = size (Proxy :: Proxy Double) + size (Proxy :: Proxy Double)
    alignment _ = alignment (Proxy :: Proxy Double)

ext: hs

max_stars_repo_name: sighingnow/computations

max_issues_repo_name: sighingnow/computations

max_forks_repo_name: sighingnow/computations


text:
module Parser (
  readExpr
              )
  where

import Text.ParserCombinators.Parsec hiding (spaces)
import System.Environment
import Control.Monad
import Numeric (readOct, readHex, readFloat)
import Data.Char (digitToInt)
import Data.Ratio ((%), Rational)
import Data.Complex (Complex (..))

data LispVal = Atom String
             | List [LispVal]
             | DottedList [LispVal] LispVal
             | Integer Integer
             | String String
             | Bool Bool
             | Ratio Rational
             | Complex (Complex Float)
             | Float Float
             | Character Char
  deriving (Eq, Show)

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

spaces :: Parser ()
spaces = skipMany1 space

myEscapes :: Parser Char
myEscapes = do
  char '\\'
  char_escaped <- oneOf ['"', 'n', 't', 'r', '\\']
  return $ case char_escaped of
    '\\' -> '\\'
    '"'  -> '"'
    'n'  -> '\n'
    'r'  -> '\r'
    't'  -> '\t'

-- test escaping quotes with: "\"\\\"\"" hahaha
parseString :: Parser LispVal
parseString = do
  char '"'
  x <- many (myEscapes <|> noneOf ['\\', '"'])
  char '"'
  return $ String x

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

parseDecimal :: Parser LispVal
parseDecimal = liftM (Integer . read) $ many1 digit

readBinRec :: [Int] -> Integer
readBinRec (lsb : rest) = (toInteger lsb) + 2 * readBinRec rest
readBinRec [] = 0

readBinReal :: String -> Integer
readBinReal = readBinRec . reverse . map digitToInt

parseBinary :: Parser LispVal
parseBinary = liftM (Integer. readBinReal) $ many1 (oneOf "01")

readOctReal :: String -> Integer
readOctReal s = fst ((readOct s) !! 0)
parseOctal :: Parser LispVal
parseOctal = liftM (Integer . readOctReal) $ many1 (oneOf "01234567")

readHexReal :: String -> Integer
readHexReal s = fst ((readHex s) !! 0)
parseHex :: Parser LispVal
parseHex = liftM (Integer . readHexReal) $ many1 (oneOf "0123456789ABCDEFabcdef")

parsePrefixedNumber :: Parser LispVal
parsePrefixedNumber = do
  char '#'
  radix <- oneOf ['b', 'o', 'd', 'x']
  case radix of
    'b' -> parseBinary
    'o' -> parseOctal
    'd' -> parseDecimal
    'x' -> parseHex

parseFloat :: Parser LispVal
parseFloat = do
  whole_num <- many1 digit
  char '.'
  fractional_bit <- many1 digit
  return . Float . fst $ (readFloat (whole_num ++ "." ++ fractional_bit) !! 0 )

parseRational :: Parser LispVal
parseRational = do
  top <- many1 digit
  char '/'
  bottom <- many1 digit
  return $ Ratio (read top % read bottom)

floatCast :: LispVal -> Float
floatCast (Float f) = f
floatCast (Integer i) = fromIntegral i

parseComplex :: Parser LispVal
parseComplex = do
  first_bit <- try parseFloat <|> parseDecimal
  char '+'
  second_bit <- try parseFloat <|> parseDecimal
  char 'i'
  (return . Complex) ((floatCast first_bit) :+ (floatCast second_bit))

-- #b (binary), #o (octal), #d (decimal), and #x (hexadecimal)
parseNumber :: Parser LispVal
parseNumber =
  try parseRational
  <|> try parseComplex
  <|> try parseFloat
  <|> try parsePrefixedNumber
  <|> parseDecimal

{-
parseNumber' :: Parser LispVal
parseNumber' = do
  digits <- many1 digit
  return ((Number . read) digits)

parseNumber'' :: Parser LispVal
parseNumber'' = many1 digit >>= (return . Number . read)
-}

{-Characters are objects that represent printed characters such as letters and digits. Characters are written using the notation #\<character> or #\<character name>. For example:

#\a	; lower case letter
#\A	; upper case letter
#\(	; left parenthesis
#\	; the space character
#\space	; the preferred way to write a space
#\newline	; the newline character
Case is significant in #\<character>, but not in #\<character name>. If <character> in #\<character> is alphabetic, then the character following <character> must be a delimiter character such as a space or parenthesis. This rule resolves the ambiguous case where, for example, the sequence of characters ``#\space'' could be taken to be either a representation of the space character or a representation of the character ``#\s'' followed by a representation of the symbol ``pace.''

Characters written in the #\ notation are self-evaluating. That is, they do not have to be quoted in programs. Some of the procedures that operate on characters ignore the difference between upper case and lower case. The procedures that ignore case have ``-ci'' (for ``case insensitive'') embedded in their names.
-}
parseChar :: Parser LispVal
parseChar = do
  char '#' >> char '\\'
  rest <- (string "space") <|> (string "newline") <|> (anyChar >>= (\x -> return  [x]))
  return (Character $ case rest of
    "space" -> ' '
    "newline" -> '\n'
    c -> c !! 0)

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

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

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

{-parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
    char '`'
    x <- parseExpr
    return $ List [Atom "quote", x]-}
parseExpr :: Parser LispVal
parseExpr = try parseAtom
         <|> try parseChar
         <|> try parseNumber
         <|> try parseString
         <|> try parseQuoted
         <|> do char '('
                x <- try parseList <|> parseDottedList
                char ')'
                return x

readExpr :: String -> String
readExpr input = case parse parseExpr "lisp" input of
    Left err -> "No match: " ++ show err ++ ", input: " ++ input
    Right val -> "Found value: " ++ show val

max_stars_repo_name: laudiacay/write-you-a-scheme

max_issues_repo_name: laudiacay/write-you-a-scheme

max_forks_repo_name: laudiacay/write-you-a-scheme


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE LambdaCase          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE TypeApplications    #-}

{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

module Neldermead (
  Simplex(..)
  , haftkaGurdalSimplex
  ) where

import Data.Ord (comparing)
import Data.Function (on)
import Data.Bool (bool)
import Data.Maybe (fromMaybe)
import Data.Monoid (Sum)
import qualified Data.Foldable as FLD
import qualified Data.Proxy as PXY
import Data.Coerce

import Control.Arrow ((&&&))
import Control.Monad.State.Strict

import qualified Numeric.Interval.NonEmpty as IVL

import qualified Data.Finite as F
import Data.Singletons (Sing, sing, fromSing)
import Data.Singletons.TypeLits

import qualified Data.Vector.Sized as V
import qualified Data.Vector.Storable.Sized as SV
--import qualified Data.Vector.Mutable.Sized as MV
import Numeric.LinearAlgebra.Static
import Numeric.LinearAlgebra.Static.Vector
import GHC.TypeNats

import qualified Control.Lens as L
import qualified Control.Lens.Iso as LI
import qualified Control.Lens.Traversal as LT

--TODO in future - mark Vertex with the objective function
--in any constructor, must take objective and then mark it

--TODO in future, allow more flexibility in objective?
--TODO better field names or allow overloaded
--TODO lenses
data NMVertex n = NMVertex { getX :: !(R n)
                           , getFx :: ℝ
                           }

newtype NMSimplex n = NMSimplex { fOrderedVertices :: V.Vector (n + 1) (NMVertex n) }

newtype Simplex n = Simplex { xA :: L n (n + 1) }
    deriving (Show)

newtype BoxConstraints n = BoxConstraints { unBoxConstraints :: V.Vector n (IVL.Interval ℝ) }

data NMScalars = NMScalars {
  reflectWeight :: ℝ
  , expandWeight :: ℝ
  , contractWeight :: ℝ
  , shrinkWeight :: ℝ
}

--TODO convergeTest should be monadic to allow writing details
--
data NMEnv n = MEnv {
  scalars :: NMScalars
  , simplexConvergeTest :: NMSimplex n -> Bool
}

newtype OnBounds = OnBounds { unOnBounds :: Bool }

makeNMVertex :: KnownNat n => (R n -> ℝ) -> R n -> NMVertex n
makeNMVertex f x = NMVertex { getX = x, getFx = f x }

--TODO actually need to change case 2 to keep if better than second worst (not best)
nextNMSimplex :: forall n. (KnownNat n, 1 <= n)
  => NMScalars -> (R n -> ℝ) -> NMSimplex n -> NMSimplex n
nextNMSimplex env f (NMSimplex vtxs)
  | comparing getFx vtxReflect vtxl == LT = undefined --expand case
  | comparing getFx vtxReflect vtxl2 == LT = undefined -- just keep swap out reflect and re-order
  | otherwise = undefined --shrink case
  where
    vtxl = V.head vtxs
    vtxl2 = V.index vtxs $ F.natToFinite (PXY.Proxy @1)
    vtxh = V.last vtxs
    xc = centroid $ getX <$> FLD.toList vtxs
    vtxReflect = makeNMVertex f $ affineCombWith (getX vtxh) (reflectWeight env) xc

--TODO may need MonadReader and MonadState, and allow objective to be in any monad
rVecIso :: KnownNat n => LI.Iso' (R n) (SV.Vector n ℝ)
rVecIso = LI.iso rVec vecR

projectToBox :: KnownNat n => BoxConstraints n -> R n -> R n
projectToBox (BoxConstraints bc) =
  L.over rVecIso $ SV.imap $ IVL.clamp . V.index bc

-- Return bool is weather any dimension was clamped
projectToBoxState :: KnownNat n
  => BoxConstraints n -> R n -> State OnBounds (R n)
projectToBoxState (BoxConstraints bc) = L.mapMOf rVecIso
  $ SV.imapM $ \i x -> do
    let bci = V.index bc i
    modify $ coerce (|| x `IVL.notMember` bci)
    return $ IVL.clamp bci x

--awkward because cant deconstruct interval into inf sup
--and cant zip storable and not storable
rightAngledSimplexInBox :: forall n. KnownNat n
  => BoxConstraints n -> R n -> R n -> Simplex n
rightAngledSimplexInBox (BoxConstraints bc) step x0 = Simplex
  $ ((x0 `outer` 1) + diag step') ||| col x0
  where
    (step' :: R n) = vecR
      $ SV.generate
      $ \i -> let stepBounds = V.index bc i - IVL.singleton (SV.index vx0 i)
                  si = abs $ SV.index vstep i
              in clampStep ((IVL.inf &&& IVL.sup) stepBounds) si
    clampStep (stepLb, stepUb) stepi
      | stepi <= stepUb = stepi
      | -stepi >= stepLb = -stepi
      | -stepLb < stepUb = max 0 $ 0.5 * stepUb
      | otherwise = min 0 $ 0.5 * stepLb
    vstep = rVec step
    vx0 = rVec x0

haftkaGurdalSimplex :: forall n. KnownNat n => ℝ -> R n -> Simplex n
haftkaGurdalSimplex size x0 = Simplex . colsL
  $ V.generate $ \j -> x0 L.& rVecIso L.%~ (SV.imap $ \i x -> x + step i j)
  where
    step i j
      | toInteger j == toInteger n' = 0
      | toInteger j == toInteger i = q + (c * l)
      | otherwise = q
    q = c * (sqrt (l + 1) - 1)
    c = size / (l * sqrt 2)
    l = fromIntegral n'
    n' = fromSing $ (sing :: Sing n)

centroid :: forall n t. (KnownNat n, Foldable t) => t (R n) -> R n
centroid = (/l) . FLD.foldl' (+) 0
  where
    l = fromIntegral . fromSing $ (sing :: Sing n)

centroidL :: forall m n. (KnownNat m, KnownNat n) => L m n -> R m
centroidL = (#> (1.0 / l))
  where
    l = fromIntegral . fromSing $ (sing :: Sing m)

realWithinRelTol :: ℝ -> ℝ -> ℝ -> Bool
realWithinRelTol relTol x y = (<= 0.5 * (x `absAdd` y)) . abs $ x - y
  where absAdd = (+) `on` abs

realWithinTols :: Maybe ℝ -> Maybe ℝ -> ℝ -> ℝ -> Bool
realWithinTols relTol absTol x y = relTolCheck && absTolCheck
  where
    relTolCheck = fromMaybe True
      $ (\t -> realWithinRelTol t x y) <$> relTol
    absTolCheck = fromMaybe True $ (<= abs (x - y)) <$> absTol

simplexRadiusTo :: (KnownNat n, Foldable t) => t (R n) -> R n -> R n
simplexRadiusTo xs c = vecR
  $ SV.generate
  $ \i -> FLD.foldl' (\acc -> (max `on` abs) acc . distci i . rVec) 0 xs
  where
    distci i v = SV.index v i - SV.index vc i
    vc = rVec c

--center + max step to given vertex for each dimension
simplexRadiusToL :: KnownNat n => Simplex n -> R n -> R n
simplexRadiusToL (Simplex xs) c = vecR $ SV.generate (V.index v)
  where
    v = (SV.maximumBy (comparing abs) . rVec) <$> lRows xsDiff -- Vector, need to convert to (R n)
    xsDiff = xs - (c `outer` 1) -- vector not promoted to matrix

affineCombWith :: KnownNat n => R n -> ℝ -> R n -> R n
affineCombWith x cy y = x + rnScale cy (y - x)

--put in helper class, maybe use coerce or something
rnScale :: KnownNat n => ℝ -> R n -> R n
rnScale c = L.over rVecIso $ SV.map (c *)

--TODO write version using mean
--toRows, map mean, list to vec
centroidL2 :: forall m n. Sing m -> Sing n -> L m n -> R m
centroidL2 (sm@SNat) (sn@SNat) = (#> (1.0 / l))
  where
    l = fromIntegral . fromSing $ (sing :: Sing m)

max_stars_repo_name: brnzhg/howod

max_issues_repo_name: brnzhg/howod

max_forks_repo_name: brnzhg/howod


text:
{-# LANGUAGE BangPatterns #-}
module NLP.Albemarle.GloVe (
  ContextVector,
  TrainingApproach,
  adagrad,
  deemph,
  gradient,
  randomModel,
  train
) where
import Numeric.LinearAlgebra (Matrix, Vector, dot)
import qualified Data.Vector.Storable as SVec
import qualified Numeric.LinearAlgebra as HMatrix
import qualified Numeric.LinearAlgebra.Data as HMatrix
import qualified Data.IntMap as IntMap
import Lens.Micro
import Lens.Micro.TH

-- Constants: these are known probably not worth parameterizing because their
-- current settings seem to work the best in all situations.
deemph_alpha = 0.75 -- Cooccurance frequencies are raised to this power
deemph_cap = 100 -- Cooccurance frequencies are clamped to [1..deemph_cap]
learning_rate = 0.05 -- Adagrad learning rate: it's not super sensitive

-- | One input or output context vector, as well as the sum of squared gradients
--   for each dimension, as needed for adagrad
data ContextVector = ContextVector {
  _bias :: !Double, -- ^ Bias associated with the vector
  _embedding :: !Vector Double, -- ^ The input or output vector
  _biasHist :: !Double, -- ^ The bias's SSE (calling it history) - for adagrad
  _embeddingHist :: !Vector Double -- ^ All other SSE's - for adagrad
}
-- | A function computing a training step. This is merely a synonym for clarity.
type TrainingApproach = ContextVector -> ContextVector -> Double
  -> (ContextVector, ContextVector)
makeLenses ''ContextVector

-- | Measure the gradient of error between two vectors. It's used as part of a
--   numerical optimization. (e.g. Adagrad)
gradient :: ContextVector -- ^ The source context vector
  -> ContextVector -- ^ The target context vector
  -> Double -- ^ Cooccurance count of the two vectors
  -> Double -- ^ The gradient of error between vectors
gradient source target edgefrequency = let
  commonness = dot (source^.embedding) (target^.embedding) -- context vector
  logfreq_error = (source^.bias) * (target^.bias) - log edgefrequency -- bias
  in commonness + logfreq_error -- both sources

-- | Deemphasize high frequencies, to avoid overweighting stopwords
deemph :: Double -> Double
deemph f = min 1 ((f/deemph_cap) ** deemph_alpha)

-- | An adaptation of Adagrad to GloVe. The difference from a typical optimizer
--   is that it is intended to operate on large streaming sources of vector
--   pairs.
adagrad :: ContextVector -- ^ Source vector
  -> ContextVector -- ^ Target vector
  -> Double -- ^ Cooccurance frequency
  -> (ContextVector, ContextVector) -- ^ Modified source and target vectors
adagrad source target edgefrequency history = let
  grad = learning_rate * gradient source target edgefrequency
  source_grad = grad * (source^.embedding)
  target_grad = grad * (target^.embedding)
  repair s t sgrad tgrad = ContextVector {
    _embedding = (s^.embedding) - (tgrad / sqrt (s^.embeddingHist))
    _bias = (s^.bias) - (grad / sqrt (s^.biasHist))
    _embeddingHist = (s^.embeddingHist) + (sgrad * sgrad)
    _biasHist = (s^.biasHist) + (grad * grad)
  }
  in (repair source target source_grad target_grad,
      repair target source target_grad source_grad)

-- | Train a model using a specific training approach, an existing model, and
--   some training cooccurances.
train :: IntMap ContextVector -- ^ The initial parameters
  -> TrainingApproach -- ^ The approach to training (like adagrad)
  -> [(Int, Int, Double)] -- ^ (Probably lazy) list of cooccurances
  -> IntMap ContextVector -- ^ Trained model
train trainer params [] = params
train trainer !params (s,t,v):insts = let
  get = (IntMap.!)
  (reps, rept) = trainer (get s) (get t) v
  in train trainer (IntMap.insert s reps $ IntMap.insert t rept params) insts

-- | Make a model out of random vectors (for bootstrapping)
randomModel :: Int -- | Number of words
  -> Int -- | Vector length
  -> IO (IntMap ContextVector) -- | Resulting map
randomModel height width = do
  mat <- HMatrix.rand height width
  return $! IntMap.fromList $ zip [0..] $ HMatrix.toRows mat

max_stars_repo_name: SeanTater/albemarle

max_issues_repo_name: SeanTater/albemarle

max_forks_repo_name: SeanTater/albemarle


text:
module Fungine.Render.Component
  ( render
  )
where

import Fungine.Prelude
import Fungine.Component
import Fungine.Graphics.Geometry as G
import Graphics.Rendering.OpenGL.GL.CoordTrans as GLC hiding (Matrix)
import Graphics.Rendering.OpenGL hiding (Matrix)
import qualified Data.Map as M
import qualified Data.Vector.Storable as VS
import Codec.Picture
import Fungine.Error
import Numeric.LinearAlgebra
import Foreign.Storable (sizeOf)
import Foreign.Ptr
import Data.Text (pack)
import Paths_fungine

-- TODO error handling
loadImage :: UIInfo -> FilePath -> IO (CanError UIInfo)
loadImage info fp = do
  p <- getDataDir
  i <- readImage (p ++ "/textures/" ++ fp)
  case (uiFreeBuffers info, i) of
    (_       , Left e    ) -> pure (Fungine.Error.Error (pack e))
    ((b : bs), Right dyni) -> do
      let (Image w h v) = convertRGBA8 dyni
      textureBinding Texture2D $= Just (TextureObject (fromIntegral b))
      VS.unsafeWith v $ \ptr -> texImage2D
        Texture2D
        NoProxy
        0
        RGBA8
        (TextureSize2D (fromIntegral w) (fromIntegral h))
        0
        (PixelData RGBA UnsignedByte ptr)
      pure $ Success $ info
        { uiLoads       = M.insert
          (UIImage fp)
          (UILoadedImage ((w, h), TextureObject (fromIntegral b)))
          (uiLoads info)
        , uiFreeBuffers = bs
        }
    ([], _) -> pure $ Fungine.Error.Error "out of free buffers for ui texture binding"

loadUITextures :: UIInfo -> UIComponent e -> IO (CanError UIInfo)
loadUITextures info c = case (M.lookup (uiRenderable c) (uiLoads info), uiRenderable c) of
  (Nothing, UINoRender  ) -> go (Success info) (uiChildren c)
  (Nothing, UIText _    ) -> go (Success info) (uiChildren c)
  (Nothing, (UIImage fp)) -> do
    img <- loadImage info fp
    go img (uiChildren c)
  (Just r, re) -> go (Success info { uiLoads = M.insert re r (uiLoads info) }) (uiChildren c)
 where
  go (Fungine.Error.Error t    ) _          = pure $ Fungine.Error.Error t
  go (Success             info') []         = pure $ Success info'
  go (Success             info') (cc : ccs) = do
    info'' <- loadUITextures info' cc
    go info'' ccs

returnUnusedBuffers :: UIInfo -> CanError UIInfo -> IO (CanError UIInfo)
returnUnusedBuffers old (Success new) =
  let mp = M.difference (uiLoads old) (uiLoads new)
  in
    do
      let objs = catMaybes $ M.foldr (\r m -> getObjectName r : m) [] mp
      pure $ Success $ new { uiFreeBuffers = objs ++ uiFreeBuffers new }

 where
  getObjectName (UILoadedImage (_, TextureObject i)) = Just i
  getObjectName _ = Nothing

returnUnusedBuffers _ (Fungine.Error.Error t) = pure $ Fungine.Error.Error t

collectEvents :: CanError UIInfo -> UIComponent e -> [e]
collectEvents _ _ = [] -- TODO actually collect events

render :: UIInfo -> UIComponent e -> IO (CanError (UIInfo, [e]))
render info c = do
  blend $= Enabled
  blendFunc $= (SrcAlpha, OneMinusSrcAlpha)
  -- generate a scale matrix based on center of screen
  let (x, y) = rectangleTopLeft $ uiBounds info
  let (w, h) = rectangleDim $ uiBounds info
  let
    view =
      tr
        $                        (G.scale3 (2 / w) (2 / h) 1)
        Numeric.LinearAlgebra.<> (G.translate (Point3d (x - w / 2) (y + h / 2) 0))
  info'  <- loadUITextures info c
  info'' <- returnUnusedBuffers info info'
  ce     <- renderComponent info'' view c
  whenError ce putStrLn -- TODO better error handling
  pure $ fmap (\i -> (i, collectEvents info'' c)) ce

-- TODO clean this function up, too long
renderComponent :: CanError UIInfo -> Matrix Float -> UIComponent e -> IO (CanError UIInfo)
renderComponent (Fungine.Error.Error t) _ _ = pure $ Fungine.Error.Error t
renderComponent (Success info) view c = -- TODO use view matrix
  let
    rot    = G.rotate (uiRotation c) (Point3d 0 0 1) -- TODO use rot
    (w, h) = rectangleDim (uiRectangle c info)
    (x, y) = rectangleTopLeft (uiRectangle c info)
    verts =
      [ x + w
      , negate (y + h)
      , 1
      , 0
      , 0
      , 0
      , 0
      , 1
      , 1
      , 0
      , 0
      , 1 -- x,y,z,r,g,b,a,uv1,uv2,norm
      , x + w
      , negate y
      , 1
      , 0
      , 0
      , 0
      , 0
      , 1
      , 0
      , 0
      , 0
      , 1
      , x
      , negate y
      , 1
      , 0
      , 0
      , 0
      , 0
      , 0
      , 0
      , 0
      , 0
      , 1
      , x
      , negate (y + h)
      , 1
      , 0
      , 0
      , 0
      , 0
      , 0
      , 1
      , 0
      , 0
      , 1
      ] :: [Float]
    shapes = [0, 1, 3, 1, 2, 3] :: [Word32]
    tx     = M.lookup (uiRenderable c) (uiLoads info)
  in case tx of
    Just (UILoadedImage (_, tobj)) -> do
      currentProgram $= Just (uiShader info)
      vao <- genObjectName
      vbo <- genObjectName
      ebo <- genObjectName
      bindVertexArrayObject $= Just vao
      bindBuffer ArrayBuffer $= Just vbo
      VS.unsafeWith (VS.fromList verts) $ \ptr ->
        bufferData ArrayBuffer
          $= (fromIntegral (length verts * sizeOf (0 :: Float)), ptr, DynamicDraw)
      bindBuffer ElementArrayBuffer $= Just ebo
      VS.unsafeWith (VS.fromList shapes) $ \ptr ->
        bufferData ElementArrayBuffer
          $= (fromIntegral (length shapes * sizeOf (0 :: Word32)), ptr, DynamicDraw)
      vertexAttribPointer (AttribLocation 0)
        $= ( ToFloat
           , VertexArrayDescriptor 3 Float (fromIntegral (12 * sizeOf (0 :: Float))) nullPtr
           )
      vertexAttribArray (AttribLocation 0) $= Enabled
      vertexAttribPointer (AttribLocation 1)
        $= ( ToFloat
           , VertexArrayDescriptor
             4
             Float
             (fromIntegral (12 * sizeOf (0 :: Float)))
             (plusPtr nullPtr (3 * sizeOf (0 :: Float)))
           )
      vertexAttribArray (AttribLocation 1) $= Enabled
      vertexAttribPointer (AttribLocation 2)
        $= ( ToFloat
           , VertexArrayDescriptor
             2
             Float
             (fromIntegral (12 * sizeOf (0 :: Float)))
             (plusPtr nullPtr (7 * sizeOf (0 :: Float)))
           )
      vertexAttribArray (AttribLocation 2) $= Enabled
      vertexAttribPointer (AttribLocation 3)
        $= ( ToFloat
           , VertexArrayDescriptor
             3
             Float
             (fromIntegral (12 * sizeOf (0 :: Float)))
             (plusPtr nullPtr (9 * sizeOf (0 :: Float)))
           )
      vertexAttribArray (AttribLocation 3) $= Enabled
      activeTexture $= TextureUnit 0
      textureBinding Texture2D $= Just tobj
      textureWrapMode Texture2D S $= (Repeated, Repeat)
      textureWrapMode Texture2D T $= (Repeated, Repeat)
      textureFilter Texture2D $= ((Linear', Nothing), Linear')

      mx  <- newMatrix ColumnMajor (VS.toList $ flatten view) :: IO (GLmatrix GLfloat)
      uni <- uniformLocation (uiShader info) "transform"
      uniform uni $= mx
      drawElements Triangles 6 UnsignedInt nullPtr
      foldM
        (\info' ch -> renderComponent info' view ch)
        (Success info { uiBounds = uiRectangle c info })
        (uiChildren c)
    _ -> case uiRenderable c of
      UINoRender -> foldM
        (\info' ch -> renderComponent info' view ch)
        (Success info { uiBounds = uiRectangle c info })
        (uiChildren c)
      _ -> return $ Fungine.Error.Error "could not find texture"

max_stars_repo_name: jmdurr/fungine

max_issues_repo_name: jmdurr/fungine

max_forks_repo_name: jmdurr/fungine


text:
{-|
Module      : Univariate
Description : Univariate polynomials (notably cyclotomics)
Copyright   : (c) Matthew Amy, 2020
Maintainer  : matt.e.amy@gmail.com
Stability   : experimental
Portability : portable
-}

module Feynman.Algebra.Polynomial.Univariate(
  Univariate,
  Cyclotomic,
  var,
  constant,
  (*|),
  primitiveRoot,
  evaluate
  )
where

import Data.List
import Data.Map(Map)
import qualified Data.Map as Map
import Data.Complex
import Data.Maybe(maybe)

import qualified Feynman.Util.Unicode as Unicode
import Feynman.Algebra.Base
import Feynman.Algebra.Polynomial

{-------------------------------
 Univariate polynomials
 -------------------------------}

-- | Univariate polynomials over the ring 'r'
data Univariate r = Univariate { getCoeffs :: !(Map Integer r) } deriving (Eq, Ord)

instance (Eq r, Num r, Show r) => Show (Univariate r) where
  show = showWithName "x"

instance Degree (Univariate r) where
  degree = maybe 0 (fromIntegral . fst) . Map.lookupMax . getCoeffs

-- | Print a univariate polynomial with a particular variable name
showWithName :: (Eq r, Num r, Show r) => String -> Univariate r -> String
showWithName x p
    | p == 0    = "0"
    | otherwise = intercalate " + " $ map showTerm (Map.assocs $ getCoeffs p)
    where showTerm (expt, a)
            | expt == 0  = show a
            | a    == 1  = showExpt expt
            | a    == -1 = "-" ++ showExpt expt
            | otherwise  = show a ++ showExpt expt
          showExpt expt
            | expt == 1  = x
            | otherwise  = Unicode.sup x expt

instance (Eq r, Num r) => Num (Univariate r) where
  (+)           = add
  (*)           = mult
  negate        = Univariate . Map.map negate . getCoeffs
  abs           = id
  signum        = id
  fromInteger 0 = Univariate Map.empty
  fromInteger i = Univariate $ Map.singleton 0 (fromInteger i)

-- | Normalize a univariate polynomial
normalize :: (Eq r, Num r) => Univariate r -> Univariate r
normalize = Univariate . Map.filter (/=0) . getCoeffs

-- | The unique univariate variable
var :: Num r => Univariate r
var = Univariate $ Map.singleton 1 1

-- | Constant polynomial
constant :: (Eq r, Num r) => r -> Univariate r
constant a
  | a == 0    = Univariate $ Map.empty
  | otherwise = Univariate $ Map.singleton 0 a

-- | Multiply by a scalar
(*|) :: (Eq r, Num r) => r -> Univariate r -> Univariate r
(*|) 0 = \_p -> zero
(*|) a = Univariate . Map.map (a*) . getCoeffs

-- | Add two univariate polynomials
add :: (Eq r, Num r) => Univariate r -> Univariate r -> Univariate r
add p = normalize . Univariate . Map.unionWith (+) (getCoeffs p) . getCoeffs

-- | Multiply two univariate polynomials
mult :: (Eq r, Num r) => Univariate r -> Univariate r -> Univariate r
mult p = normalize . Map.foldrWithKey (\expt a -> add (mulTerm expt a)) 0 . getCoeffs
  where mulTerm expt a = Univariate . Map.mapKeysMonotonic (+ expt) . Map.map (* a) $ getCoeffs p

{-------------------------------
 Cyclotomics
 -------------------------------}
-- | Cyclotomic polynomials over the ring 'r'
data Cyclotomic r = Cyc { getOrder :: !Integer, getPoly :: !(Univariate r) } deriving (Eq, Ord)

instance (Eq r, Num r, Show r) => Show (Cyclotomic r) where
  show p = showWithName (Unicode.sub Unicode.zeta (getOrder p)) $ getPoly p

instance Degree (Cyclotomic r) where
  degree = degree . getPoly

instance (Eq r, Num r) => Num (Cyclotomic r) where
  p + q = reduceOrder $ Cyc m (p' + q')
    where (Cyc m p', Cyc _ q') = unifyOrder p q
  p * q = reduceOrder $ Cyc m (p' * q')
    where (Cyc m p', Cyc _ q') = unifyOrder p q
  negate (Cyc m p) = Cyc m $ negate p
  abs    p = p
  signum p = p
  fromInteger i = Cyc 1 (fromInteger i)

-- | Unify the order of two cyclotomics
unifyOrder :: Cyclotomic r -> Cyclotomic r -> (Cyclotomic r, Cyclotomic r)
unifyOrder (Cyc n p) (Cyc m q)
  | n == m    = (Cyc n p, Cyc m q)
  | otherwise = (Cyc r p', Cyc r q')
  where r  = lcm n m
        p' = Univariate . Map.mapKeysMonotonic ((r `div` n) *) $ getCoeffs p
        q' = Univariate . Map.mapKeysMonotonic ((r `div` m) *) $ getCoeffs q

-- | Rewrite the cyclotomic in lowest order
reduceOrder :: (Eq r, Num r) => Cyclotomic r -> Cyclotomic r
reduceOrder (Cyc m c) = Cyc m' c'
  where m' = m `div` d
        c' = Univariate . Map.mapKeysMonotonic (`div` d) $ getCoeffs c
        d  = foldr gcd m . Map.keys $ getCoeffs c

-- | Construct the cyclotomic polynomial \(\zeta_i\)
primitiveRoot :: Num r => Integer -> Cyclotomic r
primitiveRoot i = Cyc i var

-- | Convert to a complex number
evaluate :: (Real r, RealFloat f) => Cyclotomic r -> Complex f
evaluate (Cyc m p) = Map.foldrWithKey f (0.0 :+ 0.0) $ getCoeffs p
  where f root coeff = (mkPolar (realToFrac coeff) (expnt root) +)
        expnt root   = 2.0*pi*(fromInteger root)/(fromInteger m)

max_stars_repo_name: PariaNaghavi/feynman

max_issues_repo_name: PariaNaghavi/feynman

max_forks_repo_name: PariaNaghavi/feynman


text:
{-# LANGUAGE TypeApplications #-}
module ArrayFire.SparseSpec where

import qualified ArrayFire       as A
import           Data.Int
import           Data.Word
import           Data.Complex
import           Data.Proxy
import           Foreign.C.Types
import           Test.Hspec

spec :: Spec
spec =
  describe "Sparse spec" $ do
    it "Should create a sparse array" $ do
      (1+1) `shouldBe` 2
      -- A.createSparseArrayFromDense (A.matrix @Double (10,10) [1..]) A.CSR
      --   `shouldBe`
      --      A.vector @Double 10 [0..]

max_stars_repo_name: danielkroeni/arrayfire-haskell

max_issues_repo_name: danielkroeni/arrayfire-haskell

max_forks_repo_name: danielkroeni/arrayfire-haskell


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE BangPatterns #-}
-- |
-- Module      : Statistics.Correlation.Pearson
--
module Statistics.Correlation
    ( -- * Pearson correlation
      pearson
    , pearsonMatByRow
      -- * Spearman correlation
    , spearman
    , spearmanMatByRow
    ) where

import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import Statistics.Matrix
import Statistics.Sample
import Statistics.Test.Internal (rankUnsorted)

----------------------------------------------------------------
-- Pearson
----------------------------------------------------------------

-- | Pearson correlation for sample of pairs. Exactly same as
-- 'Statistics.Sample.correlation'
pearson :: (G.Vector v (Double, Double), G.Vector v Double)
        => v (Double, Double) -> Double
pearson = correlation
{-# INLINE pearson #-}

-- | Compute pairwise pearson correlation between rows of a matrix
pearsonMatByRow :: Matrix -> Matrix
pearsonMatByRow m
  = generateSym (rows m)
      (\i j -> pearson $ row m i `U.zip` row m j)
{-# INLINE pearsonMatByRow #-}

----------------------------------------------------------------
-- Spearman
----------------------------------------------------------------

-- | compute spearman correlation between two samples
spearman :: ( Ord a
            , Ord b
            , G.Vector v a
            , G.Vector v b
            , G.Vector v (a, b)
            , G.Vector v Int
            , G.Vector v Double
            , G.Vector v (Double, Double)
            , G.Vector v (Int, a)
            , G.Vector v (Int, b)
            )
         => v (a, b)
         -> Double
spearman xy
  = pearson
  $ G.zip (rankUnsorted x) (rankUnsorted y)
  where
    (x, y) = G.unzip xy
{-# INLINE spearman #-}

-- | compute pairwise spearman correlation between rows of a matrix
spearmanMatByRow :: Matrix -> Matrix
spearmanMatByRow
  = pearsonMatByRow . fromRows . fmap rankUnsorted . toRows
{-# INLINE spearmanMatByRow #-}


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}

-- |
-- Module      :  Test.Convolution
-- Copyright   :  (c) 2020 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Test.Convolution (
    convolutionTests,
    linearConvolutionTests,
    cyclicConvolutionTests
  ) where

import Data.Complex
import Test.HUnit ((@?=))
import Test.Hspec

import Spiral.Convolution
import Spiral.Exp
import Spiral.SPL

convolutionTests :: Spec
convolutionTests = describe "Convolution" $ do
    linearConvolutionTests
    cyclicConvolutionTests

linearConvolutionTests :: Spec
linearConvolutionTests = describe "Linear Convolution" $ do
  sequence_ [linearConvolutionTest (Standard i) i | i <- [1..16::Int]]
  sequence_ [linearConvolutionTest (ToomCook i) i | i <- [1..16::Int]]
  linearConvolutionTest (Tensor [2,2] [Standard 2, ToomCook 2]) 4
  linearConvolutionTest (Tensor [3,5] [Standard 3, ToomCook 5]) 15
  linearConvolutionTest (Tensor [2,2,2,2] [ToomCook 2, ToomCook 2, ToomCook 2, ToomCook 2]) 16
  sequence_ [linearConvolutionTest (Lift 5 i (Standard i)) 5 | i <- [6..10::Int]]
  sequence_ [linearConvolutionTest (FromCyclic i (ConvolutionTheorem (2*i - 2))) i | i <- [2..8::Int]]

linearConvolutionTest :: LinearConvolution (Exp (Complex Double)) -> Int -> Spec
linearConvolutionTest lin n =
      it (show lin) $ toMatrix (convolve lin vector) @?= expected_result
  where
      l :: Int
      l = 2 * n - 1

      vector :: [Exp (Complex Double)]
      vector = [fromIntegral i | i <- [1..n]]

      -- expected_result :: Matrix M (Exp (Complex Double))
      expected_result = toMatrix $ transpose $ fromLists
          [replicate i 0 ++ vector ++ replicate (l - n - i) 0 | i <- [0..n-1]]

cyclicConvolutionTests :: Spec
cyclicConvolutionTests = describe "Cyclic Convolution" $ do
  sequence_ [cyclicConvolutionTest (ConvolutionTheorem i) | i <- [3..8::Int]]
  cyclicConvolutionTest (Winograd 2 [Standard 1, Standard 1])
  cyclicConvolutionTest (Winograd 3 [Standard 1, Standard 2])
  cyclicConvolutionTest (Winograd 3 [Standard 1, Lift 2 4 (Standard 4)])
  cyclicConvolutionTest (Winograd 3 [Standard 1, Lift 2 3 (Standard 3)])
  cyclicConvolutionTest (AgarwalCooley 2 3 (ConvolutionTheorem 2) (Winograd 3 [Standard 1, Standard 2]))
  cyclicConvolutionTest (AgarwalCooley 4 3 (ConvolutionTheorem 4) (Winograd 3 [Standard 1, Standard 2]))
  cyclicConvolutionTest (SplitNesting [2, 3] [(Winograd 2 [Standard 1, Standard 1]), (Winograd 3 [Standard 1, Standard 2])])
  cyclicConvolutionTest (SplitNesting [3, 2] [(Winograd 3 [Standard 1, Standard 2]), (Winograd 2 [Standard 1, Standard 1])])

cyclicConvolutionTest :: CyclicConvolution (Exp (Complex Double)) -> Spec
cyclicConvolutionTest cyc = it (show cyc) $
    toMatrix (convolve cyc vector) @?= toMatrix (circ vector)
  where
    n :: Int
    n = getSize cyc

    vector :: [Exp (Complex Double)]
    vector = [fromIntegral i | i <- [1..n]]

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-# LANGUAGE FlexibleContexts, TypeOperators, RankNTypes, GeneralizedNewtypeDeriving, DeriveFunctor, PolyKinds, GADTs #-}

module Qubit where

import Linear.V2
import Linear.V1
import Linear.Vector
import Vec
import Data.Complex
import Data.Functor.Compose
import Data.Functor.Product
import Data.Functor.Identity
import Cat
import Data.Coerce

type Kron = Compose
type Spin = V2
type QuBit = V2
type C = Complex Double

type DSum = Product

newtype FKron f g a = FKron [(f a, g a)]

type Fock a =  Q [a] -- a "fock space" over a. No assumption of anti symmettric or symmettric.
newtype Anti a = Anti (Fock a)
-- for antisymmettric space, keep list sorted (canonical form).
-- liftAnti1 :: (a -> Q a) -> ([a] -> Q [a])
-- liftAnti2 :: (a,a) -> Q (a,a) -> ([a] -> Q [a])
--

-- this is fock system as a sparse indictaor matrix. We can also do the form of occupation number  Q (V Bool)

vac :: Fock a
vac = pure []

type SHO = Q Int

shovac :: SHO
shovac = pure 0
adag :: Int -> SHO
adag n = W [(n + 1, sqrt (fromInteger (toInteger n) + 1))]
agad 0 = mempty
agad n = W [(n - 1, sqrt n)]

raise = LinOp adag
lower = LinOp agad

-- boson system can be seen as collection of oscillators. Q (V Int).  Positive integers really.

type Boson v = Q (v Int) -- accepts a functor describing all the oscillators. A Vec for example.
-- bvac = pure zero -- if v is a Vector type, it probably has a zero.
--

-- traversable ordering. The vector is put into canocnical form with adag appearing first in the traversal coming first.
type Fermion v = Q (v Bool) -- fermion needs to prefix count to figure out signs. This means that v has to be traversable? Foldable?
fvac :: Applicative v => Fermion v
fvac = pure (pure False)
{-
-- what about systems with different types of particles.
particle-hole  Q ([a], [a])

AntiOp a = AntiOp [a] -> Q [a]
-- indexed linear operators. a -> AntiOp a
psi :: Ord a => a -> (a -> Anti a)
psidag :: Ord a => a -> (a -> Anti a)

phi :: Ord a => a -> (a -> Sym a)
phidag :: Ord a => a -> (a -> Sym a)

-- hmm. Something kind of funny about these.
lift1 :: (i -> Double) -> (i -> Fock i)
lift2 :: (i -> i -> Double) -> (i -> Fock i)
symlift -- mention that there is implcit symmettrization
2 natural lifts. cnnected by logarithms. UxUxUxU vs HIII + IHII + IIHI + IIIH

sum1 :: (i -> i -> Q i) -> (i -> Q i) -- takes an indexed operator and turns it into an operator
sum1 :: (i -> AntiOp i) -> (AntiOp i)
sum2
sum :: Addable b => (i -> b) -> b

turn all of this into a typeclass functions. Then do this form, but also full initial form for manipulation.

dyson

-}
-- instance (Additive f, Additive g) => Additive (Compose f g) where
-- instance Num f (g a) => Num (Compose f g a)

-- (Kron (Kron V2 V2) V2) Complex Double  -- a 3 qubit vector space
swap :: (Traversable f, Applicative g) => Kron f g a -> Kron g f a
swap = Compose . sequenceA . getCompose

-- kron :: (Num a, Functor f, Functor g) => f a -> g a -> Kron f g a
-- kron f g = Compose $ fmap (\amp1 -> fmap (\amp2 -> amp1 * amp2) g) f
kron' :: (Num a, Functor f, Functor g) => f a -> g a -> Kron f g a
kron' x y = Compose (outer x y)

densify ::  (Additive (Kron f g), Num a, Functor f, Functor g) => FKron f g a -> Kron f g a
densify (FKron xs) = sumV $ map (uncurry kron') xs

kron''' :: (Applicative f', Applicative g', Traversable f, Traversable g') =>
    (f a -> f' a) -> (g a -> g' a) -> (Kron f g a -> Kron f' g' a)
kron''' lf lg (Compose fga) = Compose $ sequenceA $ (fmap lf) $ sequenceA  $ (fmap lg fga)

type LinOp' f g a = f a -> g a
-- SVD let's you re free
-- but also [indicatro vectro, g a]
-- and perhaps randomize sampling.
    {-
sigx :: Spin C -> Spin C
sigx (V2 up down) = V2 down up

sigz :: Spin C -> Spin C
sigz (V2 (a :+ b) (c :+ d)) = V2 ()

sigz :: Spin C -> Spin C
sigz (V2 up down) = V2 up (-down)
-}
{-
class Add a where
    (+++) :: a -> a -> a
    negate' :: a -> a

instance Num a => Add a where
    +++ = +
instance    => Additive f
-}

flip' :: Spin a -> Spin a
flip' (V2 up down) = V2 down up

sigx = flip'

sigz :: Num a => Spin a -> Spin a
sigz (V2 up down) = V2 up (negate down)

-- phase :: Complex a => Spin a -> Spin a
type f ~> g = forall a. (Functor f, Functor g) => f a -> g a
type f *** g = Compose f g
-- (Functor f, Functor g, Functor h, Functor k)
-- we only need g to be a functor
-- see below. If we remove the type signature, this implementation generalizes to another type
par :: Functor g => (forall a. f a -> g a) -> (forall b. h b -> k b) -> (forall c. (Compose f h c) -> (Compose g k c))
par f g = Compose . (fmap g) . f . getCompose -- insufficient though. We can't write linear operators this way

-- I should used V1, not Identity
-- These are all also coerce.
rightUnitor :: Functor f => Compose f Identity a -> f a
rightUnitor (Compose f) = fmap runIdentity f

rightUnitor' :: f ~> Compose f Identity
rightUnitor' = Compose . fmap Identity

leftUnitor' :: f ~> Compose Identity f
leftUnitor' = Compose . Identity

leftUnitor :: Identity *** f ~> f
leftUnitor = runIdentity . getCompose

-- should be coerce. Confusing role problems
assoc :: Functor f => (Kron (Kron f g) h) ~> (Kron f (Kron g h))
assoc = Compose . (fmap Compose) . getCompose . getCompose

assoc' :: Functor f =>  (Kron f (Kron g h)) ~> (Kron (Kron f g) h)
assoc' (Compose x)  = Compose $ Compose $ (fmap getCompose x)

-- if the functor in question happens to be representable
    -- then we can use our previous defitions of monoidal catgoeries.
    -- Reader Functors
-- Num c => prod :: Reader a c -> Reader b c -> Reader (a,b) c
-- Num
-- newtype ReadCat c a b = ReadCat (Reader a c -> Reader b c)
-- newtype LinOp c a b = LinOp ((a -> c) -> (b -> c))
-- par (LinOp f) (LinOp g) = LinOp \h -> h (\(a,b) ->   )
-- FMapN
-- fmapN @n
{-Compose f g

assoc :: Compose (Compose f g) h a -> Compose f (Compose g h) a

assoc = Compose . (fmap Compose) . getCompose . getCompose

Compose f Id a -> f a

rightunit = (fmap runIdentity) . getCompose

The morphisms of the functor category are natural transformations

par :: (f ~> g) -> (h~> k) -> ((Compose f h) ~> (Compose g k))

-}

newtype Kron' k l f a = Kron' ((k (l f)) a) deriving Functor

-- deriving via compose
{-
newtype V2' f a = V2' (Compose V2 f a)  deriving Functor
newtype V1' f a = V1' (Compose V1 f a) deriving Functor
-}
{-
newtype V2' f a = V2' (V2 (f a))  deriving Functor
newtype V1' f a = V1' (f a) deriving Functor
-}

type V2' f a = Kron V2 f a
type V1' f a = Kron V1 f a

lowerv1 (Compose (V1 f)) = f -- It's a coerce. V1 is assuredly a newtype
{-
exvsa :: (Additive f, Num a) => V2' f a
exvsa = zero
-}

-- constant vecotrs a morphisms from unit
-- Using f ~ one = (V1 1)
exconst :: (a ~ Double, Functor f) => V1' f a -> V2' f a
exconst x = let x' = lowerv1 x in Compose $ V2 (1.4 *^ x') (7.3 *^ x')

one = V1 1
attemp = exconst (Compose (V1 one))

-- Applicative composes, so Additive probably should too.
-- Wait is this in the library?
-- Why is applicative g necessary?
{-
Ok. There doesn't appear to be a compose
But V2 etc have a Num instance
and liftU2 has a default implementation that uses Applicative
I probably don't need Applcative g if I define a liftU2 using g's liftU2

-}
-- I should check to make sure the lift makes sense
-- Why isn't this in the library? There must be a good reason.
-- Is there an additive instance for Product?

instance (Additive g, Additive f) => Additive (Compose f g) where
    (Compose v) ^+^ (Compose w) = Compose (liftU2 (^+^) v w) -- I  we might be able to used
    zero = zero -- what?
    liftU2 f (Compose x) (Compose y) = Compose $ liftU2 (liftU2 f) x y
    liftI2 f (Compose x) (Compose y) = Compose $ liftI2 (liftI2 f) x y

type LinOp1 f g a = forall k. Additive k => Kron f k a -> Kron g k a

-- This feels less right to me. But why?
type LinOp1' f g a = forall k. Functor k => Kron k f a -> Kron k g a
{-
It feels less right perhaps because it will change how we do things so little.
We just have to fmap.
No, but saying we always have direct access to the scalar isn't right.
-}

-- reassoc :: Kron f (Kron g k) a -> Kron (Kron f g) k a
-- reassoc = coerce

--reassoc' ::  Kron (Kron f g) k a -> Kron f (Kron g k) a
--reassoc' = -- coerce

-- yikes. coerce will help a bit. Quite a bit.
-- Also Fixing the unnecessary applicative instances.
-- Additive k' is expected. Applicative is coming from Additve (Kron )  constraint.

{-kron'' :: (Applicative g, Applicative k', Additive k', Functor f', Functor f, Functor g') =>
  (forall k. Additive k => Kron f k a -> Kron f' k a) -> (forall k. Additive k => Kron g k a -> Kron g' k a) -> (Kron (Kron f g) k') a -> (Kron (Kron f' g') k') a
-}

type Qubit = V2

kron'' :: (Additive f, Additive g, Additive k, Additive f', Additive g') =>
  LinOp1 f f' a -> LinOp1 g g' a -> Kron (Kron f g) k a -> Kron (Kron f' g') k a
kron'' lf lg fgk = let v = (Qubit.assoc fgk) in Qubit.assoc' (Compose $ fmap lg $ getCompose (lf v))

sigx' :: LinOp1 Qubit Qubit C
sigx' (Compose (V2 up down)) = Compose $ V2 down up

sigz' :: LinOp1 Qubit Qubit C
sigz' (Compose (V2 up down)) = Compose  $ V2 up ((-1) *^ down)

sigy' :: LinOp1 Qubit Qubit C
sigy' (Compose (V2 up down)) = Compose $ V2 ((-i) *^ down) (i *^ up) where i = 0 :+ 1

swap' :: (Traversable f, Applicative g) => LinOp1 (Kron f g) (Kron g f) a
swap' (Compose (Compose fgk)) = Compose $ Compose $ sequenceA fgk

cnot :: LinOp1 (Kron Qubit Qubit) (Kron Qubit Qubit) a
cnot (Compose (Compose (V2 (V2 up1 down1) v))) = Compose $ Compose $ V2 (V2 down1 up1) v

phase :: Double -> LinOp1 Qubit Qubit C
phase phi (Compose (V2 up down)) = Compose $ V2 up ((cis phi) *^ down)

mtest :: LinOp1 (Kron Qubit Qubit) (Kron Qubit Qubit) C
mtest = (left sigx') . (left sigy') . (right sigz') . swap'

left :: (Additive f, Additive k, Additive g) => LinOp1 f g a -> LinOp1 (Kron f k) (Kron g k) a
left l = kron'' l id -- Qubit.assoc' . l . Qubit.assoc

right :: (Additive k, Additive f, Additive g) => LinOp1 f g a -> LinOp1 (Kron k f) (Kron k g) a
right l = kron'' id l -- (Compose (Compose fkk)) = Compose $ Compose $ (fmap (getCompose . l . Compose) fkk)

leftUnit :: Kron V1 (Kron f k) a -> Kron f k a
leftUnit (Compose (V1 x)) = x

{-
(f a -> g a) is a natural transfromation
(Fractional a => f a -> g a) is a numerical transfromation

Horizontal composition in Functor land. Sure.
Even in applicative
The issue is that Good linear ops are going to hold numbers. So ti feels wrong to quanitfy over them
One option is to add a type class constraint
(Num a) or (Fractional a)
This will work because linear lifts these

-- this is the same as par written above.

-- Functor f => (t -> f a) -> (a -> b) -> t -> f b
The oppotite ordering
--  Functor f' => (f' b -> t') -> (a -> b) -> f' a -> t'
The more general type as inferred is
f b ~ t'

(f' a -> f a) -> (a -> b) -> (f' a -> f b)
vs
(f' b -> f b) -> (a -> b) -> (f' a -> f b)
vs.
forall

kron l l' x = l' (fmap l x)

(forall f. Additive f => Kron V2 f a -> Kron V3 f a) -> () -> (forall p. Kron f (Kron k p) -> Kron f' (Kron k' p) )
We want to allow the field a to unify across them all.
The scalar is generic in a certain sense, but not in another.

This is the lifted Kron'
We don't want Additive1, something went screwy there.
Additive (Kron' )

The problem with Num based Kron: If we define LinOp
(LinOp f g = forall Num a => f a -> g a
Why Num, Why not Fractional. This is a problem.
And really we want to be threading the scalar a throughout the entire system.

kron

And what if we unified the two orderings

Hmm. These are pretty similar to van Laarhoven lenses. This observation is porbably poison
Swtich the arguments
(a -> b) -> Lens
Functor f => (t -> f a) -> (t -> f b)
type Lens s t a b = forall f. Functor f => (a -> f b) -> s -> f t
Lens t a t b

kron :: (Fractional a => f a -> g a) -> (Fractional k) -> LinOp
kron l l' (Compose fk) =  Compose (fmap l' (l fk))

Constrained natural transformation
kron  forall const. (const a => f a -> g a) -> (forall b. const (k b) => k b -> l b)

(f Double -> g Double) ->
(forall f. V2' f Double -> V2' f Double)

kron :: LinOp f g -> LinOp k l -> LinOp (Kron f k) (Kron g l)

-}

{-

-- can't build. Well we can build vzero
exv2' :: (Additive f, Num a) => V2' f a
exv2' = V2' (Compose ( ?  ))

-}
-- constants are given as morphisms
-- exv2' :: (Additive f, Num a) => V1' f a -> V2' f a
-- exv2' (V1' f) = (V2' (V2   )

{-
There is no additive instance for Compose?
Huh.

class (forall f. Additive f => Additive k f) => Additive1 k where

I'm doing something insane.
Yea, The 1 pattern. I've been here before and decided it was insane.
Something in purescript. Surely it wasn't Additive1 was it?

-}
{-
class Additive1 k where
    vadd :: (Additive f, Num a) => k f a -> k f a -> k f a
    vzero :: (Additive f, Num a) => k f a
instance Additive1 V2' where -- (Functor (V2' f)) =>

    vadd (V2' (V2 v w)) (V2' (V2 x y)) = V2' (V2 (v ^+^ x)  (w ^+^ y)   )
    vzero = V2' (V2 zero zero)

    --vadd (V2' v) (V2' w) = V2' (v ^+^ w)
    --vzero = V2' zero

instance Additive1 V1' where -- (Functor (V2' f)) =>
    vadd (V1' v) (V1' w) = V1' (v ^+^ w)
    vzero = V1' zero

instance (Additive1 k, Additive1 l) => Additive1 (Kron' k l) where
   vadd (Kron' x) (Kron' y) = Kron' (vadd x y)
   vzero = Kron' vzero
   -}
{-
Deriv functor for V2'
-}


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLabels      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-|
Module      : Grenade.Layers.LeakyRelu
Description : Rectifying linear unit layer
Copyright   : (c) Manuel Schneckenreither, 2020
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.LeakyRelu (
    LeakyRelu (..)
  ) where

import           Control.DeepSeq                (NFData (..))
import           Data.Proxy                     (Proxy(Proxy))
import           Data.Maybe                     (fromJust)
import           Data.Serialize
import           GHC.Generics                   (Generic)
import           GHC.TypeLits
import qualified Numeric.LinearAlgebra.Static   as H

import           Grenade.Core
import           Grenade.Onnx
import           Grenade.Types
import           Grenade.Layers.Internal.Activations

-- | A rectifying linear unit.
--   A layer which can act between any shape of the same dimension, acting as a
--   diode on every neuron individually.
newtype LeakyRelu = LeakyRelu RealNum
  deriving (Generic, NFData, Show)

instance UpdateLayer LeakyRelu where
  type Gradient LeakyRelu = ()
  runUpdate _ x _ = x
  reduceGradient _ = ()

instance RandomLayer LeakyRelu where
  createRandomWith _ _ = return $ LeakyRelu 0.01

instance Serialize LeakyRelu where
  put (LeakyRelu alpha) = put alpha
  get = LeakyRelu <$> get

instance (KnownNat i) => Layer LeakyRelu ('D1 i) ('D1 i) where
  type Tape LeakyRelu ('D1 i) ('D1 i) = S ('D1 i)

  runForwards (LeakyRelu alpha) (S1D y) = (S1D y, (S1D . fromJust . H.create) relu)
    where
      w    = fromIntegral $ natVal (Proxy :: Proxy i)
      relu = leakyRelu1d w alpha (H.extract y)

  runBackwards (LeakyRelu alpha) (S1D y) (S1D dEdy) = ((), S1D (relu' y * dEdy))
    where
      relu' = H.dvmap (\a -> if a < 0 then alpha else 1)

instance (KnownNat i, KnownNat j) => Layer LeakyRelu ('D2 i j) ('D2 i j) where
  type Tape LeakyRelu ('D2 i j) ('D2 i j) = S ('D2 i j)

  runForwards (LeakyRelu alpha) (S2D y) = (S2D y, (S2D . fromJust . H.create) relu)
    where
      h    = fromIntegral $ natVal (Proxy :: Proxy i)
      w    = fromIntegral $ natVal (Proxy :: Proxy j)
      y'   = H.extract y
      relu = leakyRelu 1 h w alpha y'

  runBackwards (LeakyRelu alpha) (S2D y) (S2D dEdy) = ((), S2D (relu' y * dEdy))
    where
      relu' = H.dmmap (\a -> if a < 0 then alpha else 1)

instance (KnownNat i, KnownNat j, KnownNat k) => Layer LeakyRelu ('D3 i j k) ('D3 i j k) where

  type Tape LeakyRelu ('D3 i j k) ('D3 i j k) = S ('D3 i j k)

  runForwards (LeakyRelu alpha) (S3D y) = (S3D y, (S3D . fromJust . H.create) relu)
    where
      c    = fromIntegral $ natVal (Proxy :: Proxy k)
      h    = fromIntegral $ natVal (Proxy :: Proxy i)
      w    = fromIntegral $ natVal (Proxy :: Proxy j)
      y'   = H.extract y
      relu = leakyRelu c h w alpha y'

  runBackwards (LeakyRelu alpha) (S3D y) (S3D dEdy) = ((), S3D (relu' y * dEdy))
    where
      relu' = H.dmmap (\a -> if a < 0 then alpha else 1)

instance OnnxOperator LeakyRelu where
  onnxOpTypeNames _ = ["LeakyRelu"]

instance OnnxLoadable LeakyRelu where
  loadOnnxNode _ node = do
    alpha <- readFloatAttributeToRealNum "alpha" node

    return $ LeakyRelu alpha