text:
{-# LANGUAGE CPP, DataKinds, GADTs, KindSignatures, RankNTypes,
             StandaloneDeriving, TypeOperators #-}
#if !MIN_VERSION_base(4,16,0)
{-# LANGUAGE IncoherentInstances #-}
#endif
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}
-- | The second component of dual numbers, @Delta@, with its semantics.
-- Neel Krishnaswami calls it \"sparse vector expressions\",
-- and indeed even in the simplest case of an objective function
-- defined on scalars only, the codomain of the function that computes
-- gradients from such delta expressions is a set of vectors, because
-- the gradient of an @R^n@ to @R@ function is an @R^n@ vector.
--
-- The \'sparsity\' is less obvious when the domain of the function consists
-- of multiple vectors, matrices and tensors and when the expressions themselves
-- contain vectors, matrices and tensors. However, a single tiny delta
-- expression (e.g., a sum of two variables) may denote a vector of matrices.
-- Even a delta expression containing a big matrix denotes something much
-- bigger: a whole vector of such matrices and more.
--
-- The algebraic structure here is an extension of vector space.
-- The crucial extra constructor of a variable is used both to represent
-- sharing in order to avoid exponential blowup and to replace the one-hot
-- access to parameters with something cheaper and more uniform.
-- A lot of the remaining additional structure is for introducing
-- and reducing dimensions (ranks).
--
-- This is an internal API now, superseded by "HordeAd.Core.DualClass"
-- that permits other kinds of second component of dual numbers,
-- e.g., the same as primal component, for fast computation
-- of forward derivatives (because @derivativeFromDelta@ below,
-- computing derivatives from delta-expressions, is slow once the expressions
-- grow large enough to affect cache misses).
module HordeAd.Internal.Delta
  ( -- * Abstract syntax trees of the delta expressions
    Delta0 (..), Delta1 (..), Delta2 (..), DeltaX (..), DeltaS (..)
  , -- * Delta expression identifiers
    DeltaId, toDeltaId, convertDeltaId
  , -- * Evaluation of the delta expressions
    DeltaBinding
  , DeltaState (..)
  , Domain0, Domain1, Domain2, DomainX, Domains
  , gradientFromDelta, derivativeFromDelta, ppBindings
  , bindInState0, bindInState1, bindInState2, bindInStateX
  , isTensorDummy
  , CodeOut(..)
  ) where

import Prelude

import           Control.Exception (assert)
import           Control.Monad (unless, zipWithM_)
import           Control.Monad.ST.Strict (ST, runST)
import qualified Data.Array.Convert
import qualified Data.Array.Dynamic as OTB
import qualified Data.Array.DynamicS as OT
import qualified Data.Array.Internal
import qualified Data.Array.Internal.DynamicG
import qualified Data.Array.Internal.DynamicS
import qualified Data.Array.Shaped as OSB
import qualified Data.Array.ShapedS as OS
import           Data.Kind (Type)
import           Data.List (foldl')
import           Data.Proxy (Proxy)
import qualified Data.Strict.Vector as Data.Vector
import qualified Data.Strict.Vector.Autogen.Mutable as Data.Vector.Mutable
import qualified Data.Vector.Generic as V
import qualified Data.Vector.Generic.Mutable as VM
import qualified Data.Vector.Storable.Mutable
import           GHC.TypeLits (KnownNat, Nat, natVal, type (+))
import           Numeric.LinearAlgebra (Matrix, Numeric, Vector, (#>), (<.>))
import qualified Numeric.LinearAlgebra as HM
import           Text.Show.Pretty (ppShow)

import qualified HordeAd.Internal.MatrixOuter as MO
import           HordeAd.Internal.OrthotopeOrphanInstances (liftVS2, liftVT2)

-- * Abstract syntax trees of the delta expressions

-- | This is the grammar of delta-expressions at tensor rank 0, that is,
-- at scalar level. The first few operations have analogues
-- at the level of vectors, matrices and arbitrary tensors.
--
-- For each choice of the underlying scalar type @r@,
-- we have several primitive differentiable types based on the scalar:
-- the scalar type @r@ itself, @Vector r@, @Matrix r@ and tensors.
-- Many operations span the ranks and so span the datatypes, which makes
-- the datatypes mutually recursive.
--
-- The @Outline@ constructors represent primitive numeric function applications
-- for which we delay computing and forgo inlining of the derivative.
data Delta0 r =
    Zero0
  | Scale0 r (Delta0 r)
  | Add0 (Delta0 r) (Delta0 r)
  | Var0 (DeltaId r)

  | SumElements0 (Delta1 r) Int  -- ^ see Note [SumElements0]
  | Index0 (Delta1 r) Int Int  -- ^ second integer is the length of the vector

  | Dot0 (Vector r) (Delta1 r)  -- ^ Dot0 v vd == SumElements0 (Scale1 v vd) n

  | FromX0 (DeltaX r)  -- ^ one of many conversions
  | FromS0 (DeltaS '[] r)

  | Outline0 CodeOut [r] [Delta0 r]
  | Delay0 ~(Delta0 r)

deriving instance (Show r, Numeric r) => Show (Delta0 r)

-- | This is the grammar of delta-expressions at tensor rank 1, that is,
-- at vector level.
data Delta1 r =
    Zero1
  | Scale1 (Vector r) (Delta1 r)
  | Add1 (Delta1 r) (Delta1 r)
  | Var1 (DeltaId (Vector r))

  | Seq1 (Data.Vector.Vector (Delta0 r))  -- ^ "unboxing" conversion
  | Konst1 (Delta0 r) Int  -- ^ length; needed only for forward derivative
  | Append1 (Delta1 r) Int (Delta1 r)  -- ^ the length of the first argument
  | Slice1 Int Int (Delta1 r) Int  -- ^ last integer is the length of argument
  | SumRows1 (Delta2 r) Int  -- ^ the integer is the number of columns
  | SumColumns1 (Delta2 r) Int  -- ^ the integer is the number of rows

  | M_VD1 (Matrix r)
          (Delta1 r)  -- ^ M_VD1 m vd == SumRows1 (M_MD2 m (AsRow2 vd))
  | MD_V1 (Delta2 r)
          (Vector r)  -- ^ MD_V1 md v == SumRows1 (MD_M2 md (asRow v))

  | FromX1 (DeltaX r)
  | forall len. KnownNat len
    => FromS1 (DeltaS '[len] r)

    -- unsorted and undocumented yet
  | Reverse1 (Delta1 r)
  | Flatten1 Int Int (Delta2 r)
  | FlattenX1 OT.ShapeL (DeltaX r)
  | forall sh. OS.Shape sh
    => FlattenS1 (DeltaS sh r)

  | Outline1 CodeOut [Vector r] [Delta1 r]
  | Delay1 ~(Delta1 r)

deriving instance (Show r, Numeric r) => Show (Delta1 r)

-- | This is the grammar of delta-expressions at tensor rank 2, that is,
-- at matrix level.
data Delta2 r =
    Zero2
  | Scale2 (Matrix r) (Delta2 r)
  | Add2 (Delta2 r) (Delta2 r)
  | Var2 (DeltaId (Matrix r))

  | FromRows2 (Data.Vector.Vector (Delta1 r))  -- ^ "unboxing" conversion again
  | FromColumns2 (Data.Vector.Vector (Delta1 r))
  | Konst2 (Delta0 r) (Int, Int)  -- ^ size; needed only for forward derivative
  | Transpose2 (Delta2 r)
  | M_MD2 (Matrix r) (Delta2 r)  -- ^ matrix-(matrix-expression) multiplication
  | MD_M2 (Delta2 r) (Matrix r)  -- ^ (matrix-expression)-matrix multiplication
  | RowAppend2 (Delta2 r) Int (Delta2 r)  -- ^ row-length of first argument
  | ColumnAppend2 (Delta2 r) Int (Delta2 r)  -- ^ col-length of first argument
  | RowSlice2 Int Int (Delta2 r) Int  -- ^ last arg is row-length of the matrix
  | ColumnSlice2 Int Int (Delta2 r) Int  -- ^ column-length of the matrix

  | AsRow2 (Delta1 r)  -- ^ AsRow2 vd == FromRows2 (V.replicate n vd)
  | AsColumn2 (Delta1 r)  -- ^ AsColumn2 vd == FromColumns2 (V.replicate n vd)

  | FromX2 (DeltaX r)
  | forall rows cols. (KnownNat rows, KnownNat cols)
    => FromS2 (DeltaS '[rows, cols] r)

    -- unsorted and undocumented yet
  | Flipud2 (Delta2 r)
  | Fliprl2 (Delta2 r)
  | Reshape2 Int (Delta1 r)
  | Conv2 (Matrix r) (Delta2 r)

  | Outline2 CodeOut [Matrix r] [Delta2 r]
  | Delay2 ~(Delta2 r)

deriving instance (Show r, Numeric r) => Show (Delta2 r)

-- | This is the grammar of delta-expressions at arbitrary tensor rank.
--
-- Warning: not tested enough nor benchmarked.
data DeltaX r =
    ZeroX
  | ScaleX (OT.Array r) (DeltaX r)
  | AddX (DeltaX r) (DeltaX r)
  | VarX (DeltaId (OT.Array r))

  | KonstX (Delta0 r) OT.ShapeL  -- ^ size; needed only for forward derivative
  | AppendX (DeltaX r) Int (DeltaX r)
      -- ^ Append two arrays along the outermost dimension.
      -- All dimensions, except the outermost, must be the same.
      -- The integer argument is the outermost size of the first array.
  | SliceX Int Int (DeltaX r) Int
      -- ^ Extract a slice of an array along the outermost dimension.
      -- The extracted slice must fall within the dimension.
      -- The last argument is the outermost size of the argument array.
  | IndexX (DeltaX r) Int Int
      -- ^ The sub-tensors at the given index of the outermost dimension.
      -- The second integer is the length of the dimension.
  | RavelFromListX [DeltaX r]
      -- ^ Create a tensor from a list treated as the outermost dimension.
  | ReshapeX OT.ShapeL OT.ShapeL (DeltaX r)
      -- ^ Change the shape of the tensor from the first to the second.

  | From0X (Delta0 r)
  | From1X (Delta1 r)
  | From2X (Delta2 r) Int
  | forall sh. OS.Shape sh
    => FromSX (DeltaS sh r)

  | OutlineX CodeOut [OT.Array r] [DeltaX r]
  | DelayX ~(DeltaX r)

deriving instance (Show r, Numeric r) => Show (DeltaX r)

-- | This is the grammar of delta-expressions at arbitrary tensor rank,
-- the fully typed Shaped version.
--
-- Warning: not tested enough nor benchmarked.
data DeltaS :: [Nat] -> Type -> Type where
  ZeroS :: DeltaS sh r
  ScaleS :: OS.Array sh r -> DeltaS sh r -> DeltaS sh r
  AddS :: DeltaS sh r -> DeltaS sh r -> DeltaS sh r
  VarS :: DeltaId (OS.Array sh r) -> DeltaS sh r

  KonstS :: Delta0 r -> DeltaS sh r
  AppendS :: (OS.Shape sh, KnownNat m, KnownNat n)
          => DeltaS (m ': sh) r -> DeltaS (n ': sh) r
          -> DeltaS ((m + n) ': sh) r
    -- ^ Append two arrays along the outermost dimension.
  SliceS :: (KnownNat i, KnownNat n, KnownNat k, OS.Shape rest)
         => Proxy i -> Proxy n -> DeltaS (i + n + k ': rest) r
         -> DeltaS (n ': rest) r
    -- ^ Extract a slice of an array along the outermost dimension.
  IndexS :: (KnownNat ix, KnownNat k, OS.Shape rest)
         => DeltaS (ix + 1 + k ': rest) r -> Proxy ix -> DeltaS rest r
    -- ^ The sub-tensors at the given index of the outermost dimension.
  RavelFromListS :: (KnownNat k, OS.Shape rest)
                 => [DeltaS rest r] -> DeltaS (k : rest) r
    -- ^ Create a tensor from a list treated as the outermost dimension.
  ReshapeS :: (OS.Shape sh, OS.Shape sh', OS.Size sh ~ OS.Size sh')
           => DeltaS sh r -> DeltaS sh' r
    -- ^ Change the shape of the tensor.

  From0S :: Delta0 r -> DeltaS '[] r
  From1S :: Delta1 r -> DeltaS '[n] r
  From2S :: KnownNat cols
         => Proxy cols -> Delta2 r -> DeltaS '[rows, cols] r
  FromXS :: DeltaX r -> DeltaS sh r

  OutlineS :: CodeOut -> [OS.Array sh r] -> [DeltaS sh r] -> DeltaS sh r
  DelayS :: ~(DeltaS sh r) -> DeltaS sh r

instance Show (DeltaS sh r) where
  show _ = "a DeltaS delta expression"

-- * Delta expression identifiers

newtype DeltaId a = DeltaId Int
  deriving Show

toDeltaId :: Int -> DeltaId a
toDeltaId = DeltaId

convertDeltaId :: DeltaId (OT.Array r) -> DeltaId (OS.Array sh r)
convertDeltaId (DeltaId i) = DeltaId i

-- The key property is that it preserves the phantom type.
succDeltaId :: DeltaId a -> DeltaId a
succDeltaId (DeltaId i) = DeltaId (succ i)

-- * Evaluation of the delta expressions

-- | Binding at one of the ranks, with a given underlying scalar.
--
-- The 'DeltaId' components could be re-computed on the fly in 'buildFinMaps',
-- but it costs more (they are boxed, so re-allocation is expensive)
-- than storing them here at the time of binding creation and accessing
-- in `buildFinMaps`.
data DeltaBinding r =
    DeltaBinding0 (DeltaId r) (Delta0 r)
  | DeltaBinding1 (DeltaId (Vector r)) (Delta1 r)
  | DeltaBinding2 (DeltaId (Matrix r)) (Delta2 r)
  | DeltaBindingX (DeltaId (OT.Array r)) (DeltaX r)

data DeltaState r = DeltaState
  { deltaCounter0 :: DeltaId r
  , deltaCounter1 :: DeltaId (Vector r)
  , deltaCounter2 :: DeltaId (Matrix r)
  , deltaCounterX :: DeltaId (OT.Array r)
  , deltaBindings :: [DeltaBinding r]
  }

-- | Helper definitions to shorten type signatures. Note that these
-- differ from their counterparts in all other modules, because the type
-- argument here is the underlying scalar (e.g., @Double),
-- while elsewhere it's the dual component of dual numbers from
-- rank 0 (scalar) level (e.g., @Delta0 Double@).
-- By chance, these definitions and definitions from other modules
-- coincide in case of "forward derivatives computed on the spot"
-- where @r@ is @Double@ and @Double@ is also the dual component.
--
-- More generally, @r@ in this module tends to refer to the underlying
-- scalar type, while in all other modules it refers to the rank 0 dual
-- component type.
type Domain0 r = Vector r

type Domain1 r = Data.Vector.Vector (Vector r)

type Domain2 r = Data.Vector.Vector (Matrix r)

type DomainX r = Data.Vector.Vector (OT.Array r)

type Domains r = (Domain0 r, Domain1 r, Domain2 r, DomainX r)

-- | Delta expressions naturally denote forward derivatives,
-- as encoded in function 'derivativeFromDelta'. However, we are more
-- interested in computing gradients, which is what @gradientFromDelta@ does.
-- The two functions are bound by the equation from Lemma 5 from the paper
-- "Provably correct, asymptotically efficient, higher-order reverse-mode
-- automatic differentiation":
--
-- > dt <.> derivativeFromDelta st d ds = gradientFromDelta st d dt <.> ds
--
-- where @\<.\>@ denotes generalized dot product (multiplying
-- all tensors element-wise and summing the results),
-- @st@ contains bindings of delta variables and @d@ is the top level
-- delta expression from translation of the objective function @f@ to dual
-- numbers, @ds@ belongs to the domain of @f@ and @dt@ to the codomain.
-- We omitted for clarity the @dim0@, @dim1@, @dim2@ and @dimX@ arguments
-- that are the lengths of vectors of the tensors in the domain of @f@.
--
-- Intuitively, @ds@ is a tiny perturbation of the arguments of @f@,
-- for which we compute the derivative, that is, the induced change
-- in the result of @f@. Similarly, @dt@ is a tiny perturbation of the
-- result of @f@, for which we compute the gradient, that is, the change
-- of arguments of @f@ sufficient to cause the perturbation.
-- Note that the scaling factor @r@ in functions @eval*@ in @gradientFromDelta@
-- locally plays the role of @dt@, just as the argument @parameters@
-- in @eval*@ in @derivativeFromDelta@ corresponds to @ds@.
--
-- Let's first discuss in detail the semantics of delta-expressions
-- in terms of forward derivatives, since it's more straightforward.
-- Let @r@ be the type of underlying scalars. Let @f@ be a mathematical
-- differentiable function that takes a collection of type @C@
-- of arguments and produces a single result of type @r@.
-- Let a dual number counterpart of @f@ applied to a collection
-- of parameters @P@ of type @C@ be represented as a Haskell value @b@.
-- Let @d :: Delta0 r@ be the closed delta expression that is the second
-- component of @b@, let @ds@ belong to @C@. The semantics of @d@ is a linear
-- function from @C@ to @r@ that is the derivative of @f@ at point @P@
-- with respect to the perturbation @ds@. The mathematical formula
-- for the derivative follows straightforwardly the syntactic form
-- of the delta expression @d@ (see 'derivativeFromDelta').
--
-- Let's now describe the semantics of closed delta expression @d@
-- as the gradient of @f@ at point @P@ with respect to a @dt@ that belongs
-- to @r@. Here the semantics of @d@ is a collection of four finite maps
-- (vectors) @v0@, @v1@, @v2@, @vX@, corresponding to @C@,
-- each map @vi@ taking indexes of type @DeltaId ai@ to values of type @ai@,
-- where @a0@ is @r@, @a1@ is @Vector r@, @a2@ is @Matrix r@
-- and @aX@ is the type of tensors of @r@.
-- The value of @vi@ at index @DeltaId k@ is the partial derivative
-- of function @f@ at @P@ with respect to its parameter of type @ai@.
-- The position of the @ai@ parameter is represented by @DeltaId k@
-- (in other words, the partial derivative is with respect to a variable
-- quantity tagged with @DeltaId k@) and its value comes from @dt@.
--
-- The semantics of a delta expression that is not closed but contains
-- occurrences of variables that do not correspond to parameters of @f@ is only
-- defined in the context of four vectors that contain values associated
-- to its free variables or, alternatively, of bindings from which the values
-- can be computed, or of a mixture of both. This semantics does not change
-- if a bound expression is substituted for a variable instead of being used
-- to compute a value. (Note however that a computed value can't be
-- substituted for all occurrences of the variable in an expression,
-- because the "computing backwards" trick, needed to get gradients
-- from derivative expressions, computes a value for each occurrence
-- of a variable separately and sums over all occurrences instead
-- of substituting a single value into each occurrence.)
--
-- Function @gradientFromDelta@ computes the four vectors described above.
-- Requested lengths of the vectors are given in the first few arguments.
-- The delta state contains a list of mutually-referencing delta bindings
-- that are to be evaluated, in the given order, starting with the top-level
-- binding of a scalar type provided in the next argument and with respect
-- to perturbation @dt@ (usually set to @1@) in the last argument.
gradientFromDelta :: (Eq r, Numeric r, Num (Vector r))
                  => (CodeOut -> [r] -> [Delta0 r] -> Delta0 r)
                  -> (CodeOut -> [Vector r] -> [Delta1 r] -> Delta1 r)
                  -> (CodeOut -> [Matrix r] -> [Delta2 r] -> Delta2 r)
                  -> (CodeOut -> [OT.Array r] -> [DeltaX r] -> DeltaX r)
                  -> (forall sh. OS.Shape sh
                      => CodeOut -> [OS.Array sh r] -> [DeltaS sh r]
                      -> DeltaS sh r)
                  -> Int -> Int -> Int -> Int -> DeltaState r -> Delta0 r -> r
                  -> Domains r
gradientFromDelta inlineDerivative0 inlineDerivative1 inlineDerivative2
                  inlineDerivativeX inlineDerivativeS
                  dim0 dim1 dim2 dimX st deltaTopLevel dt =
  -- This is morally @V.create@ and so totally safe,
  -- but we can't just call @V.create@ thrice, because it would run
  -- the @ST@ action thrice, so we inline and extend @V.create@ here.
  runST $ do
    (finMap0, finMap1, finMap2, finMapX) <-
      buildFinMaps inlineDerivative0 inlineDerivative1 inlineDerivative2
                   inlineDerivativeX inlineDerivativeS
                   st deltaTopLevel dt
    v0 <- V.unsafeFreeze $ VM.take dim0 finMap0
    v1 <- V.unsafeFreeze $ VM.take dim1 finMap1
    v2 <- V.unsafeFreeze $ VM.take dim2 finMap2
    vX <- V.unsafeFreeze $ VM.take dimX finMapX
    -- Convert to normal matrices, but only the portion of vector
    -- that is not discarded.
    return (v0, v1, V.map MO.convertMatrixOuterOrNull v2, vX)

-- | Create vectors (representing finite maps) that hold delta-variable
-- values. They are initialized with dummy values so that it's cheap to check
-- if any update has already been performed to a cell (allocating big matrices
-- filled with zeros is too costly, especially if never used in an iteration,
-- and adding to such matrices and especially using them as scaling factors
-- is wasteful). The vectors are longer than those representing objective
-- function parameters (e.g., @deltaCounter0@ vs @dim0@), because variables
-- represent not only parameters, but also the bindings that prevent blowup
-- via delta-expression duplication.
initializeFinMaps :: forall s r. Numeric r
                  => DeltaState r
                  -> ST s ( Data.Vector.Storable.Mutable.MVector s r
                          , Data.Vector.Mutable.MVector s (Vector r)
                          , Data.Vector.Mutable.MVector s (MO.MatrixOuter r)
                          , Data.Vector.Mutable.MVector s (OT.Array r) )
initializeFinMaps st = do
  let DeltaId counter0 = deltaCounter0 st
      DeltaId counter1 = deltaCounter1 st
      DeltaId counter2 = deltaCounter2 st
      DeltaId counterX = deltaCounterX st
  finMap0 <- VM.replicate counter0 0  -- correct value
  finMap1 <- VM.replicate counter1 (V.empty :: Vector r)  -- dummy value
  finMap2 <- VM.replicate counter2 MO.emptyMatrixOuter  -- dummy value
  finMapX <- VM.replicate counterX dummyTensor
  return (finMap0, finMap1, finMap2, finMapX)

buildFinMaps :: forall s r. (Eq r, Numeric r, Num (Vector r))
             => (CodeOut -> [r] -> [Delta0 r] -> Delta0 r)
             -> (CodeOut -> [Vector r] -> [Delta1 r] -> Delta1 r)
             -> (CodeOut -> [Matrix r] -> [Delta2 r] -> Delta2 r)
             -> (CodeOut -> [OT.Array r] -> [DeltaX r] -> DeltaX r)
             -> (forall sh. OS.Shape sh
                 => CodeOut -> [OS.Array sh r] -> [DeltaS sh r]
                 -> DeltaS sh r)
             -> DeltaState r -> Delta0 r -> r
             -> ST s ( Data.Vector.Storable.Mutable.MVector s r
                     , Data.Vector.Mutable.MVector s (Vector r)
                     , Data.Vector.Mutable.MVector s (MO.MatrixOuter r)
                     , Data.Vector.Mutable.MVector s (OT.Array r) )
buildFinMaps inlineDerivative0 inlineDerivative1 inlineDerivative2
             inlineDerivativeX inlineDerivativeS
             st deltaTopLevel dt = do
  (finMap0, finMap1, finMap2, finMapX) <- initializeFinMaps st
  let addToVector :: Vector r -> Vector r -> Vector r
      addToVector r = \v -> if V.null v then r else v + r
      addToMatrix :: MO.MatrixOuter r -> MO.MatrixOuter r -> MO.MatrixOuter r
      addToMatrix r = \v -> if MO.nullMatrixOuter v then r else MO.plus v r
      addToArray :: OT.Array r -> OT.Array r -> OT.Array r
      addToArray r = \v -> if isTensorDummy v then r else liftVT2 (+) v r
      addToArrayS :: OS.Shape sh => OS.Array sh r -> OT.Array r -> OT.Array r
      addToArrayS r = \v -> let rs = Data.Array.Convert.convert r
                            in if isTensorDummy v
                               then rs
                               else liftVT2 (+) v rs
      eval0 :: r -> Delta0 r -> ST s ()
      eval0 !r = \case
        Zero0 -> return ()
        Scale0 k d -> eval0 (k * r) d
        Add0 d e -> eval0 r d >> eval0 r e
        Var0 (DeltaId i) -> VM.modify finMap0 (+ r) i

        SumElements0 vd n -> eval1 (HM.konst r n) vd
        Index0 (Var1 (DeltaId i)) ix k -> do
          let f v = if V.null v
                    then HM.konst 0 k V.// [(ix, r)]
                    else v V.// [(ix, v V.! ix + r)]
          VM.modify finMap1 f i
            -- this would be an asymptotic optimization compared to
            -- the general case below, if not for the non-mutable update,
            -- which involves copying the whole vector, so it's just
            -- several times faster (same allocation, but not adding vectors)
        Index0 d ix k -> eval1 (HM.konst 0 k V.// [(ix, r)]) d

        Dot0 v vd -> eval1 (HM.scale r v) vd

        FromX0 d -> evalX (OT.scalar r) d
        FromS0 d -> evalS (OS.scalar r) d

        Outline0 codeOut primalArgs dualArgs ->
          eval0 r $ inlineDerivative0 codeOut primalArgs dualArgs
        Delay0 d -> eval0 r d
      eval1 :: Vector r -> Delta1 r -> ST s ()
      eval1 !r = \case
        Zero1 -> return ()
        Scale1 k d -> eval1 (k * r) d
        Add1 d e -> eval1 r d >> eval1 r e
        Var1 (DeltaId i) -> VM.modify finMap1 (addToVector r) i

        Seq1 lsd -> V.imapM_ (\i d -> eval0 (r V.! i) d) lsd
        Konst1 d _n -> V.mapM_ (`eval0` d) r
        Append1 d k e -> eval1 (V.take k r) d >> eval1 (V.drop k r) e
        Slice1 i n d len ->
          eval1 (HM.konst 0 i V.++ r V.++ HM.konst 0 (len - i - n)) d
        SumRows1 dm cols -> eval2 (MO.asColumn r cols) dm
        SumColumns1 dm rows -> eval2 (MO.asRow r rows) dm

        M_VD1 m dRow ->
          mapM_ (`eval1` dRow)
                (MO.toRows (MO.MatrixOuter (Just m) (Just r) Nothing))
        MD_V1 md row -> eval2 (MO.MatrixOuter Nothing (Just r) (Just row)) md

        FromX1 d -> evalX (OT.fromVector [V.length r] r) d
        FromS1 d -> evalS (OS.fromVector r) d

        Reverse1 d -> eval1 (V.reverse r) d
        Flatten1 rows cols d ->
          eval2 (MO.MatrixOuter (Just $ rows HM.>< cols $ V.toList r)
                                Nothing Nothing)
                d
        FlattenX1 sh d -> evalX (OT.fromVector sh r) d
        FlattenS1 d -> evalS (OS.fromVector r) d

        Outline1 codeOut primalArgs dualArgs ->
          eval1 r $ inlineDerivative1 codeOut primalArgs dualArgs
        Delay1 d -> eval1 r d
      eval2 :: MO.MatrixOuter r -> Delta2 r -> ST s ()
      eval2 !r = \case
        Zero2 -> return ()
        Scale2 k d -> eval2 (MO.multiplyWithOuter k r) d
        Add2 d e -> eval2 r d >> eval2 r e
        Var2 (DeltaId i) -> VM.modify finMap2 (addToMatrix r) i

        FromRows2 lvd -> zipWithM_ eval1 (MO.toRows r) (V.toList lvd)
        FromColumns2 lvd -> zipWithM_ eval1 (MO.toColumns r) (V.toList lvd)
        Konst2 d _sz -> mapM_ (V.mapM_ (`eval0` d)) $ MO.toRows r
        Transpose2 md -> eval2 (MO.transpose r) md  -- TODO: test!
        M_MD2 m md ->
          let mo = MO.MatrixOuter (Just $ HM.tr' m) Nothing Nothing
          in eval2 (MO.matMul mo r) md
        MD_M2 md m ->
          let mo = MO.MatrixOuter (Just $ HM.tr' m) Nothing Nothing
          in eval2 (MO.matMul r mo) md
        RowAppend2 d k e -> eval2 (MO.takeRows k r) d
                            >> eval2 (MO.dropRows k r) e
        ColumnAppend2 d k e -> eval2 (MO.takeColumns k r) d
                               >> eval2 (MO.dropColumns k r) e
        RowSlice2 i n d rows ->
          assert (MO.rows r == n) $
          let cols = MO.cols r
          in eval2 (MO.konst 0 i cols
                    `MO.rowAppend` r
                    `MO.rowAppend` MO.konst 0 (rows - i - n) cols)
                   d
        ColumnSlice2 i n d cols ->
          assert (MO.cols r == n) $
          let rows = MO.rows r
          in eval2 (MO.konst 0 rows i
                    `MO.columnAppend` r
                    `MO.columnAppend` MO.konst 0 rows (cols - i - n))
                   d

        AsRow2 dRow -> mapM_ (`eval1` dRow) (MO.toRows r)
        AsColumn2 dCol -> mapM_ (`eval1` dCol) (MO.toColumns r)

        FromX2 d -> evalX (OT.fromVector [MO.rows r, MO.cols r]
                                         (V.concat $ MO.toRows r)) d
        FromS2 d -> evalS (OS.fromVector $ V.concat $ MO.toRows r) d

        Flipud2 d -> eval2 (MO.flipud r) d
        Fliprl2 d -> eval2 (MO.fliprl r) d
        Reshape2 _cols d -> eval1 (V.concat $ MO.toRows r) d
        Conv2 m md ->
          let mor = MO.convertMatrixOuter r
              convolved = HM.corr2 m mor
              moc = MO.MatrixOuter (Just convolved) Nothing Nothing
          in eval2 moc md

        Outline2 codeOut primalArgs dualArgs ->
          eval2 r $ inlineDerivative2 codeOut primalArgs dualArgs
        Delay2 d -> eval2 r d
      evalX :: OT.Array r -> DeltaX r -> ST s ()
      evalX !r = \case
        ZeroX -> return ()
        ScaleX k d -> evalX (liftVT2 (*) k r) d
        AddX d e -> evalX r d >> evalX r e
        VarX (DeltaId i) -> VM.modify finMapX (addToArray r) i

        KonstX d _sz -> mapM_ (`eval0` d) $ OT.toList r
        AppendX d k e -> case OT.shapeL r of
          n : _ -> evalX (OT.slice [(0, k)] r) d
                   >> evalX (OT.slice [(k, n - k)] r) e
          [] -> error "evalX: appending a 0-dimensional tensor"
        SliceX i n d len -> case OT.shapeL r of
          n' : rest ->
            assert (n' == n) $
            evalX (OT.concatOuter [ OT.constant (i : rest) 0
                                  , r
                                  , OT.constant (len - i - n : rest) 0 ])
                  d
          [] -> error "evalX: slicing a 0-dimensional tensor"
        IndexX d ix len ->
          let rest = OT.shapeL r
          in evalX (OT.concatOuter [ OT.constant (ix : rest) 0
                                   , OT.reshape (1 : rest) r
                                   , OT.constant (len - ix - 1 : rest) 0 ])
                   d  -- TODO: optimize for Var case
        RavelFromListX ld -> do
          let lr = OTB.toList $ OT.unravel r
          mapM_ (uncurry evalX) (zip lr ld)
        ReshapeX sh _sh' d -> evalX (OT.reshape sh r) d

        From0X d -> eval0 (OT.unScalar r) d
        From1X d -> eval1 (OT.toVector r) d
        From2X d cols ->
          eval2 (MO.MatrixOuter (Just $ HM.reshape cols $ OT.toVector r)
                                Nothing Nothing)
                d
        FromSX d -> evalS (Data.Array.Convert.convert r) d

        OutlineX codeOut primalArgs dualArgs ->
          evalX r $ inlineDerivativeX codeOut primalArgs dualArgs
        DelayX d -> evalX r d
      evalS :: OS.Shape sh
            => OS.Array sh r -> DeltaS sh r -> ST s ()
      evalS !r = \case
        ZeroS -> return ()
        ScaleS k d -> evalS (liftVS2 (*) k r) d
        AddS d e -> evalS r d >> evalS r e
        VarS (DeltaId i) -> VM.modify finMapX (addToArrayS r) i

        KonstS d -> mapM_ (`eval0` d) $ OS.toList r
        AppendS (d :: DeltaS (k ': rest) r) (e :: DeltaS (l ': rest) r) ->
          evalS (OS.slice @'[ '(0, k) ] r) d
          >> evalS (OS.slice @'[ '(k, l) ] r) e
        SliceS (_ :: Proxy i) _ (d :: DeltaS (i_plus_n_plus_k ': rest) r) ->
          evalS (OS.constant @(i ': rest) 0
                 `OS.append` r
                 `OS.append` OS.constant 0)
                d
        IndexS (d :: DeltaS (ix_plus_1_plus_k ': rest) r) (_ :: Proxy ix) ->
          evalS (OS.constant @(ix ': rest) 0
                 `OS.append` OS.reshape r
                 `OS.append` OS.constant 0)
                d  -- TODO: optimize for Var case
        RavelFromListS ld -> do
          let lr = OSB.toList $ OS.unravel r
          mapM_ (uncurry evalS) (zip lr ld)
        ReshapeS d -> evalS (OS.reshape r) d

        From0S d -> eval0 (OS.unScalar r) d
        From1S d -> eval1 (OS.toVector r) d
        From2S proxyCols d ->
          eval2 (MO.MatrixOuter
                   (Just $ HM.reshape (fromInteger $ natVal proxyCols)
                         $ OS.toVector r)
                   Nothing Nothing)
                d
        FromXS d -> evalX (Data.Array.Convert.convert r) d

        OutlineS codeOut primalArgs dualArgs ->
          evalS r $ inlineDerivativeS codeOut primalArgs dualArgs
        DelayS d -> evalS r d

  eval0 dt deltaTopLevel

  let evalUnlessZero :: DeltaBinding r -> ST s ()
      evalUnlessZero (DeltaBinding0 (DeltaId i) d) = do
        r <- finMap0 `VM.read` i
        unless (r == 0) $  -- we init with exactly 0.0 so the comparison works
          eval0 r d
      evalUnlessZero (DeltaBinding1 (DeltaId i) d) = do
        r <- finMap1 `VM.read` i
        unless (V.null r) $
          eval1 r d
      evalUnlessZero (DeltaBinding2 (DeltaId i) d) = do
        r <- finMap2 `VM.read` i
        unless (MO.nullMatrixOuter r) $
          eval2 r d
      evalUnlessZero (DeltaBindingX (DeltaId i) d) = do
        r <- finMapX `VM.read` i
        unless (isTensorDummy r) $
          evalX r d
  mapM_ evalUnlessZero (deltaBindings st)
  return (finMap0, finMap1, finMap2, finMapX)

-- | Forward derivative computation via forward-evaluation of delta-expressions
-- (which is surprisingly competitive to the direct forward method,
-- until the allocation of deltas gets large enough to affect cache hits).
-- This is the directional derivative, calculated for the point,
-- at which the delta expression was computed (which is the point
-- represented by the parameters of the objective function and used
-- to compute it's dual number result) and along the direction vector(s)
-- given in the last parameter called @ds@.
derivativeFromDelta
  :: forall r. (Numeric r, Num (Vector r))
  => (CodeOut -> [r] -> [Delta0 r] -> Delta0 r)
  -> (CodeOut -> [Vector r] -> [Delta1 r] -> Delta1 r)
  -> (CodeOut -> [Matrix r] -> [Delta2 r] -> Delta2 r)
  -> (CodeOut -> [OT.Array r] -> [DeltaX r] -> DeltaX r)
  -> (forall sh. OS.Shape sh
      => CodeOut -> [OS.Array sh r] -> [DeltaS sh r]
      -> DeltaS sh r)
  -> DeltaState r -> Delta0 r -> Domains r -> r
derivativeFromDelta inlineDerivative0 inlineDerivative1 inlineDerivative2
                    inlineDerivativeX inlineDerivativeS
                    st deltaTopLevel
                    _ds@(params0Init, params1Init, params2Init, paramsXInit) =
  let eval0 :: Domains r -> Delta0 r -> r
      eval0 parameters@(params0, _, _, _) = \case
        Zero0 -> 0
        Scale0 k d -> k * eval0 parameters d
        Add0 d e -> eval0 parameters d + eval0 parameters e
        Var0 (DeltaId i) -> params0 V.! i

        SumElements0 vd _n -> HM.sumElements $ eval1 parameters vd
        Index0 d ix _k -> eval1 parameters d V.! ix

        Dot0 vr vd -> vr <.> eval1 parameters vd

        FromX0 d -> OT.unScalar $ evalX parameters d
        FromS0 d -> OS.unScalar $ evalS parameters d

        Outline0 codeOut primalArgs dualArgs ->
          eval0 parameters $ inlineDerivative0 codeOut primalArgs dualArgs
        Delay0 d -> eval0 parameters d
      eval1 :: Domains r -> Delta1 r -> Vector r
      eval1 parameters@(_, params1, _, _) = \case
        Zero1 -> 0
        Scale1 k d -> k * eval1 parameters d
        Add1 d e -> eval1 parameters d + eval1 parameters e
        Var1 (DeltaId i) -> params1 V.! i

        Seq1 lsd -> V.convert $ V.map (eval0 parameters) lsd
        Konst1 d n -> HM.konst (eval0 parameters d) n
        Append1 d _k e -> eval1 parameters d V.++ eval1 parameters e
        Slice1 i n d _len -> V.slice i n $ eval1 parameters d
        SumRows1 dm _cols ->
          V.fromList $ map HM.sumElements $ HM.toRows $ eval2 parameters dm
        SumColumns1 dm _rows ->
          V.fromList $ map HM.sumElements $ HM.toColumns $ eval2 parameters dm

        M_VD1 m dRow -> m #> eval1 parameters dRow
        MD_V1 md row -> eval2 parameters md #> row

        FromX1 d -> OT.toVector $ evalX parameters d
        FromS1 d -> OS.toVector $ evalS parameters d

        Reverse1 d -> V.reverse $ eval1 parameters d
        Flatten1 _rows _cols d -> HM.flatten $ eval2 parameters d
        FlattenX1 _sh d -> OT.toVector $ evalX parameters d
        FlattenS1 d -> OS.toVector $ evalS parameters d

        Outline1 codeOut primalArgs dualArgs ->
          eval1 parameters $ inlineDerivative1 codeOut primalArgs dualArgs
        Delay1 d -> eval1 parameters d
      eval2 :: Domains r -> Delta2 r -> Matrix r
      eval2 parameters@( _, _, params2, _) = \case
        Zero2 -> 0
        Scale2 k d -> k * eval2 parameters d
        Add2 d e -> eval2 parameters d + eval2 parameters e
        Var2 (DeltaId i) -> params2 V.! i

        FromRows2 lvd ->
          HM.fromRows $ map (eval1 parameters) $ V.toList lvd
        FromColumns2 lvd ->
          HM.fromColumns $ map (eval1 parameters) $ V.toList lvd
        Konst2 d sz -> HM.konst (eval0 parameters d) sz
        Transpose2 md -> HM.tr' $ eval2 parameters md
        M_MD2 m md -> m HM.<> eval2 parameters md
        MD_M2 md m -> eval2 parameters md HM.<> m
        RowAppend2 d _k e -> eval2 parameters d HM.=== eval2 parameters e
        ColumnAppend2 d _k e -> eval2 parameters d HM.||| eval2 parameters e
        RowSlice2 i n d _rows ->
          HM.takeRows n $ HM.dropRows i $ eval2 parameters d
        ColumnSlice2 i n d _cols ->
          HM.takeColumns n $ HM.dropColumns i $ eval2 parameters d

        AsRow2 dRow -> HM.asRow $ eval1 parameters dRow  -- TODO: risky
        AsColumn2 dCol -> HM.asColumn $ eval1 parameters dCol  -- TODO: risky

        FromX2 d ->
          let t = evalX parameters d
          in case OT.shapeL t of
            [_rows, cols] -> HM.reshape cols $ OT.toVector t
            _ -> error "eval2: wrong tensor dimensions"
        FromS2 d ->
          let t = evalS parameters d
          in case OS.shapeL t of
            [_rows, cols] -> HM.reshape cols $ OS.toVector t
            _ -> error "eval2: wrong tensor dimensions"

        Flipud2 d -> HM.flipud $ eval2 parameters d
        Fliprl2 d -> HM.fliprl $ eval2 parameters d
        Reshape2 cols d -> HM.reshape cols $ eval1 parameters d
        Conv2 m md -> HM.conv2 m $ eval2 parameters md

        Outline2 codeOut primalArgs dualArgs ->
          eval2 parameters $ inlineDerivative2 codeOut primalArgs dualArgs
        Delay2 d -> eval2 parameters d
      evalX :: Domains r -> DeltaX r -> OT.Array r
      evalX parameters@( _, _, _, paramsX) = \case
        ZeroX -> 0
        ScaleX k d -> k * evalX parameters d
        AddX d e -> evalX parameters d + evalX parameters e
        VarX (DeltaId i) -> paramsX V.! i

        KonstX d sz -> OT.constant sz $ eval0 parameters d
        AppendX d _k e -> evalX parameters d `OT.append` evalX parameters e
        SliceX i n d _len -> OT.slice [(i, n)] $ evalX parameters d
        IndexX d ix _len -> OT.index (evalX parameters d) ix
        RavelFromListX ld ->
          let la = map (evalX parameters) ld
              sh = case la of
                a : _ -> length la : OT.shapeL a
                [] -> []
          in OT.ravel $ OTB.fromList sh la
        ReshapeX _sh sh' d -> OT.reshape sh' $ evalX parameters d

        From0X d -> OT.scalar $ eval0 parameters d
        From1X d -> let v = eval1 parameters d
                    in OT.fromVector [V.length v] v
        From2X d cols -> let l = eval2 parameters d
                         in OT.fromVector [HM.rows l, cols] $ HM.flatten l
        FromSX d -> Data.Array.Convert.convert $ evalS parameters d

        OutlineX codeOut primalArgs dualArgs ->
          evalX parameters $ inlineDerivativeX codeOut primalArgs dualArgs
        DelayX d -> evalX parameters d
      evalS :: OS.Shape sh => Domains r -> DeltaS sh r -> OS.Array sh r
      evalS parameters@( _, _, _, paramsX) = \case
        ZeroS -> 0
        ScaleS k d -> k * evalS parameters d
        AddS d e -> evalS parameters d + evalS parameters e
        VarS (DeltaId i) -> Data.Array.Convert.convert $ paramsX V.! i

        KonstS d -> OS.constant $ eval0 parameters d
        AppendS d e -> evalS parameters d `OS.append` evalS parameters e
        SliceS (_ :: Proxy i) (_ :: Proxy n) d ->
          OS.slice @'[ '(i, n) ] $ evalS parameters d
        IndexS d proxyIx ->
          OS.index (evalS parameters d) (fromInteger $ natVal proxyIx)
        RavelFromListS ld ->
          let la = map (evalS parameters) ld
          in OS.ravel $ OSB.fromList la
        ReshapeS d -> OS.reshape $ evalS parameters d

        From0S d -> OS.scalar $ eval0 parameters d
        From1S d -> OS.fromVector $ eval1 parameters d
        From2S _ d -> OS.fromVector $ HM.flatten $ eval2 parameters d
        FromXS d -> Data.Array.Convert.convert $ evalX parameters d

        OutlineS codeOut primalArgs dualArgs ->
          evalS parameters $ inlineDerivativeS codeOut primalArgs dualArgs
        DelayS d -> evalS parameters d
      evalUnlessZero :: Domains r -> DeltaBinding r -> Domains r
      evalUnlessZero parameters@(!params0, !params1, !params2, !paramsX) = \case
        DeltaBinding0 (DeltaId i) d ->
          let v = eval0 parameters d
          in (params0 V.// [(i, v)], params1, params2, paramsX)
        DeltaBinding1 (DeltaId i) d ->
          let v = eval1 parameters d
          in (params0, params1 V.// [(i, v)], params2, paramsX)
        DeltaBinding2 (DeltaId i) d ->
          let v = eval2 parameters d
          in (params0, params1, params2 V.// [(i, v)], paramsX)
        DeltaBindingX (DeltaId i) d ->
          let v = evalX parameters d
          in (params0, params1, params2, paramsX V.// [(i, v)])
      parameters1 = runST $ do
        (finMap0, finMap1, outerFinMap2, finMapX) <- initializeFinMaps st
        -- We use normal hmatrix matrices rather than the sparse replacement.
        finMap2 <- VM.replicate (VM.length outerFinMap2) (HM.fromRows [])
        -- TODO: the following coredumps without the @VM.take@; it's a shame
        -- there's no copying of a smaller vector into a larger one in the API.
        -- Perhaps use https://hackage.haskell.org/package/base-4.16.0.0/docs/Foreign-Marshal-Array.html#v:copyArray?
        V.unsafeCopy (VM.take (V.length params0Init) finMap0) params0Init
        V.unsafeCopy (VM.take (V.length params1Init) finMap1) params1Init
        V.unsafeCopy (VM.take (V.length params2Init) finMap2) params2Init
        V.unsafeCopy (VM.take (V.length paramsXInit) finMapX) paramsXInit
        v0 <- V.unsafeFreeze finMap0
        v1 <- V.unsafeFreeze finMap1
        v2 <- V.unsafeFreeze finMap2
        vX <- V.unsafeFreeze finMapX
        return (v0, v1, v2, vX)
      parametersB = foldl' evalUnlessZero parameters1
                           (reverse $ deltaBindings st)
  in eval0 parametersB deltaTopLevel

-- | This is yet another semantics of delta-expressions and their
-- bindings --- by pretty-printing as texts.
ppBindings :: (Show r, Numeric r) => Bool -> DeltaState r -> Delta0 r -> String
ppBindings reversed st deltaTopLevel =
  let pp = if reversed
           then foldl' (\ !l b -> l ++ ppBinding "where" b)
                       ["COMPUTE " ++ ppShow deltaTopLevel ++ "\n"]
           else foldl' (\ !l b -> ppBinding "let" b ++ l)
                       ["in " ++ ppShow deltaTopLevel ++ "\n"]
  in concat $ pp $ deltaBindings st

ppBinding :: (Show r, Numeric r) => String -> DeltaBinding r -> [String]
ppBinding prefix = \case
  DeltaBinding0 (DeltaId i) d ->
    [prefix ++ "0 DeltaId_", show i, " = ", ppShow d, "\n"]
  DeltaBinding1 (DeltaId i) d ->
    [prefix ++ "1 DeltaId_", show i, " = ", ppShow d, "\n"]
  DeltaBinding2 (DeltaId i) d ->
    [prefix ++ "2 DeltaId_", show i, " = ", ppShow d, "\n"]
  DeltaBindingX (DeltaId i) d ->
    [prefix ++ "X DeltaId_", show i, " = ", ppShow d, "\n"]

bindInState0 :: Delta0 r -> DeltaState r -> (DeltaState r, DeltaId r)
{-# INLINE bindInState0 #-}
bindInState0 u' st =
  let dId = deltaCounter0 st
      !binding = DeltaBinding0 dId u'
  in ( st { deltaCounter0 = succDeltaId dId
          , deltaBindings = binding : deltaBindings st
          }
     , dId )

bindInState1 :: Delta1 r -> DeltaState r -> (DeltaState r, DeltaId (Vector r))
{-# INLINE bindInState1 #-}
bindInState1 u' st =
  let dId = deltaCounter1 st
      !binding = DeltaBinding1 dId u'
  in ( st { deltaCounter1 = succDeltaId dId
          , deltaBindings = binding : deltaBindings st
          }
     , dId )

bindInState2 :: Delta2 r -> DeltaState r -> (DeltaState r, DeltaId (Matrix r))
{-# INLINE bindInState2 #-}
bindInState2 u' st =
  let dId = deltaCounter2 st
      !binding = DeltaBinding2 dId u'
  in ( st { deltaCounter2 = succDeltaId dId
          , deltaBindings = binding : deltaBindings st
          }
     , dId )

bindInStateX :: DeltaX r -> DeltaState r -> (DeltaState r, DeltaId (OT.Array r))
{-# INLINE bindInStateX #-}
bindInStateX u' st =
  let dId = deltaCounterX st
      !binding = DeltaBindingX dId u'
  in ( st { deltaCounterX = succDeltaId dId
          , deltaBindings = binding : deltaBindings st
          }
     , dId )

{- Note [SumElements0]
~~~~~~~~~~~~~~~~~~~~~~

The second argument of SumElements0 is the length of the vector
to be summed. Given that we sum a delta-expression representing
a vector, we can't call Vector.length on it, so the length needs
to be recorded in the constructor. Alternatively, it could be
recorded in the Delta1 argument to SumElements0. This is what
shaped tensors do at the type level, so for DeltaS the argument
would not be needed.

Sum of vector elements can be implemented using a delta-expression
primitive SumElements0 as well as without this primitive, referring
only to the primitive Index0:

https://github.com/Mikolaj/horde-ad/blob/d069a45773ed849913b5ebd0345153072f304fd9/src/HordeAd/Core/DualNumber.hs#L125-L143

which is confirmed by tests to be equivalent in three different
implementations:

https://github.com/Mikolaj/horde-ad/blob/d069a45773ed849913b5ebd0345153072f304fd9/test/TestSingleGradient.hs#L116-L128

and proved to be prohibitively slow in the two implementations
that don't use the SumElements0 primitive in benchmarks (despite
an ingenious optimization of the common case of Index0 applied to a variable):

https://github.com/Mikolaj/horde-ad/blob/d069a45773ed849913b5ebd0345153072f304fd9/bench/BenchProdTools.hs#L178-L193
-}

dummyTensor :: Numeric r => OT.Array r
dummyTensor =  -- an inconsistent tensor array
  Data.Array.Internal.DynamicS.A
  $ Data.Array.Internal.DynamicG.A []
  $ Data.Array.Internal.T [] (-1) V.empty

isTensorDummy :: OT.Array r -> Bool
isTensorDummy (Data.Array.Internal.DynamicS.A
                 (Data.Array.Internal.DynamicG.A _
                    (Data.Array.Internal.T _ (-1) _))) = True
isTensorDummy _ = False

data CodeOut =
    PlusOut | MinusOut | TimesOut | NegateOut | AbsOut | SignumOut
  | DivideOut | RecipOut
  | ExpOut | LogOut | SqrtOut | PowerOut | LogBaseOut
  | SinOut | CosOut | TanOut | AsinOut | AcosOut | AtanOut
  | SinhOut | CoshOut | TanhOut | AsinhOut | AcoshOut | AtanhOut
  | Atan2Out
  deriving Show

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
-- Slope estimates

module Model.Slope
  ( estimateSlope
  , median
  , getSurroundingSlopes
  ) where

import           Data.Vector.Algorithms.Intro (partialSort)
import qualified Data.Vector.Unboxed          as V
import           Statistics.Test.Types
import           Statistics.Test.WilcoxonT    (wilcoxonMatchedPairTest)
import           Statistics.Types

import           Types.Indices

estimateSlope
  :: IVector Index1 Double
  -> IIntSet Index1
  -> Int
  -> Index1
  -> Double
estimateSlope ds badSegments radius i
  | V.length nonZeroSlopes < 5 = 0 -- p < 0.10 impossible with l.t. 5 slopes
  | pValue' < 0.10             = median slopes
  | otherwise                  = 0
  where
    slopes = V.fromList $ getSurroundingSlopes ds badSegments radius i
    nonZeroSlopes = V.filter (/= 0) slopes
    pairs = V.zip nonZeroSlopes (V.replicate (V.length nonZeroSlopes) 0)
    pValue' = pValue $ testSignificance
            $ wilcoxonMatchedPairTest SamplesDiffer pairs

-- Fails if vector is empty
median :: V.Vector Double -> Double
median v
  | odd n = v' V.! middleIndex
  | otherwise = mid (v' V.! pred middleIndex) (v' V.! middleIndex)
  where
    n = V.length v
    middle = n `quot` 2 + 1
    middleIndex = pred middle
    v' = V.modify (`partialSort` middle) v
    mid x y = (x+y)/2

getSurroundingSlopes
  :: IVector Index1 Double
  -> IIntSet Index1
  -> Int
  -> Index1
  -> [Double]
getSurroundingSlopes ds badSegments radius i =
  map (ivIndex ds) $ filter (not . flip iisMember badSegments) [i0..i1]
  where
    (i0, i1) = runIndexInterval $ iiBoundByIVector ds
             $ IndexInterval (iTranslate (-radius) i, iTranslate radius i)

ext: hs

max_stars_repo_name: awjchen/SonoCleaner

max_issues_repo_name: awjchen/SonoCleaner

max_forks_repo_name: awjchen/SonoCleaner


text:
module Jumps (findJumps,getJumps,flt2dbl) where

import Prelude hiding (filter,length,replicate,take,(++),splitAt)

import Data.Function (on)
import Data.Ord (comparing)
import Data.Foldable (foldMap')
import Data.Bifunctor (bimap)
import Data.Vector.Unboxed (filter,slice,length,convert,splitAt)
import Data.Vector (replicate,take,force,modify,(++),maxIndexBy)
import Data.Vector.Mutable (write)

import qualified Data.Vector.Unboxed as U
import qualified Data.Vector as V

import Numeric.Sum (sumVector,kbn)
import Statistics.Sample (mean)
import GHC.Float (float2Double)

data Partition = Partition
    { start :: Int
    , num :: Int
    , level :: Double
    , _loss :: Double
    , _split :: Int
    , diff :: Double
    , _levels :: (Double,Double)
    , _losses :: (Double,Double)
    } deriving Show

data Split = Split
    { _num :: Int
    , _lvls :: (Double,Double)
    , _lsses :: (Double,Double)
    , loss :: Double
    }

instance Eq Split where
    (==) = (==) `on` loss

instance Ord Split where
    (<=) = (<=) `on` loss

instance Semigroup Split where
    (<>) = min

instance Monoid Split where
    mempty = Split 0 (0,0) (0,0) (1 / 0) -- 1 / 0 = Infinity

mkPart :: U.Vector Double -> Int -> Double -> Double -> Partition
mkPart xs a lvl lss
    | m <= 1 = Partition a m lvl lss a 0 (0,0) (0,0)
    | otherwise =
        case flip foldMap' [1..m-1] $ \i ->
            let filt = filter (not . isNaN)
                (x1,x2) = bimap filt filt $ splitAt i xs
                m1 = mean x1
                m2 = mean x2
                s1 = sumVector kbn $ U.map (\x -> (x - m1) ^ (2 :: Int)) x1
                s2 = sumVector kbn $ U.map (\x -> (x - m2) ^ (2 :: Int)) x2
            in Split i (m1,m2) (s1,s2) (s1 + s2)
        of Split n lvls lses lss' -> Partition a m lvl lss n (lss - lss') lvls lses
    where m = length xs

splitPart :: U.Vector Double -> Partition -> (Partition,Partition)
splitPart inData (Partition a m _ _ n _ (lvl1,lvl2) (lss1,lss2)) =
    ( mkPart (slice a n inData) a lvl1 lss1
    , mkPart (slice (a+n) (m-n) inData) (a+n) lvl2 lss2
    )

findJumps :: U.Vector Double -> V.Vector Partition
findJumps inData =
    let noNaN = filter (not . isNaN) inData
        lev = mean noNaN
        los = sumVector kbn $ U.map (\x -> (x - lev) ^ (2 :: Int)) noNaN
        vec = replicate 50 $ mkPart inData 0 lev los
    in loop vec 1 0 (diff $ vec V.! 0)
    where
        loop :: V.Vector Partition -> Int -> Int -> Double -> V.Vector Partition
        loop parts k splitInd lossDiff
            | lossDiff <= 120 = force $ take k parts
            | k `mod` 50 == 0 = loop (parts ++ take 50 parts) k splitInd lossDiff
            | otherwise =
                let (x,y) = splitPart inData $ parts V.! splitInd
                    parts' = modify (\v -> write v splitInd x >> write v k y) parts
                    splitInd' = maxIndexBy (comparing diff) parts'
                    lossDiff' = diff $ parts' V.! splitInd'
                in loop parts' (k + 1) splitInd' lossDiff'

getJumps :: V.Vector Partition -> U.Vector (Double,Int,Int)
getJumps = convert . (V.map $ (,,) <$> level <*> start <*> num)

flt2dbl :: U.Vector Float -> U.Vector Double
flt2dbl = U.map float2Double


text:
--Hasquil registers
--Copyright Laurence Emms 2018

module Register (Quantum(..),
                 Classical(..)) where

import Data.Complex

data Quantum = QubitRegister Int |
               MetaQubitRegister String

instance Show Quantum where
    show (QubitRegister i) = show i
    show (MetaQubitRegister s) = s

data Classical a = Register Int |
                   Range Int Int |
                   RealConstant a |
                   ComplexConstant (Complex a) |
                   MetaRegister String

instance (Floating a, Show a, Ord a) => Show (Classical a) where
    show (Register i) = "[" ++ (show i) ++ "]"
    show (Range i j) = "[" ++ (show i) ++ "-" ++ (show j) ++ "]"
    show (RealConstant r) = show r
    show (ComplexConstant (p :+ q))
        | q >= 0 = (show p) ++ "+" ++ (show q) ++ "i"
        | otherwise = (show p) ++ (show q) ++ "i"
    show (MetaRegister s) = s

max_stars_repo_name: WhatTheFunctional/Hasquil

max_issues_repo_name: WhatTheFunctional/Hasquil

max_forks_repo_name: WhatTheFunctional/Hasquil


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell #-}

module Day15 where

import Control.Lens
import Control.Lens.TH

import qualified Data.ByteString as BS
import Data.Attoparsec.ByteString (Parser)
import qualified Data.Attoparsec.ByteString as DAB
import qualified Data.Attoparsec.ByteString.Char8 as DABC8
import Data.Attoparsec.ByteString.Char8 (decimal, inClass, digit, endOfLine, char, anyChar, letter_ascii, notChar)

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

import qualified Data.List as DL
import qualified Data.List.Split as DLS

import Debug.Trace

import Control.Arrow

import Data.Array.IArray (Array)
import qualified Data.Array.IArray as IA

import Data.Set (Set)
import qualified Data.Set as Set
import Data.Either (isRight)
import Data.Maybe (mapMaybe, catMaybes)
import Data.Bool (bool)

import Data.Vector (Vector, (!))
import qualified Data.Vector as Vec

import Data.Function (on)

import Control.Monad.State.Strict (State)
import qualified Control.Monad.State.Strict as State

import qualified Control.Monad as CM

import qualified Control.Foldl as L

import Data.Bifunctor (bimap)

import Data.Complex
import Data.Functor

import Data.Word
import Data.Bits
import Data.Semigroup
import Data.List.NonEmpty (NonEmpty(..))
import qualified Data.List.NonEmpty as DLNE

type Input = Integer

solver input fTurn = State.evalState (code 6 (fromIntegral $ length input + 1)) $ -- turns start at length input + 1
	Map.fromList $ zip input $ fmap (,Nothing) [1..] --Insert the seed input with their order
	where
	-- Monadic loop. Maybe avoid custom recursion
	-- It returns the fTurn-th value
	code :: Integer -> Integer -> State (Map Integer (Integer, Maybe Integer)) Integer
	code prev turn
		| turn > fTurn = pure prev
		| otherwise =
			-- Check is the previous value is repeated
			State.get >>= \m -> case Map.lookup prev m of
				Just (_, Nothing) -> do --not repeated
					State.modify (Map.insertWith inserter 0 (turn, Nothing))
					code 0 (turn + 1)
				Just (l, Just ll) -> do -- repeated
					State.modify (Map.insertWith inserter (l-ll) (turn, Nothing))
					code (l - ll) (turn + 1)

	-- Update the map keeping the previous position
	-- DONE: Use a tupple (Integer, Maybe Integer)
	-- DONE: avoid error (removed safety check)
	inserter :: (Integer, Maybe Integer) -> (Integer, Maybe Integer) -> (Integer, Maybe Integer)
	inserter (new, Nothing) (prev, _) = (new, Just prev)
		-- | new > prev =  (new, Just prev)
		-- | otherwise = error $ show (prev, new)

runSolution :: FilePath -> IO ()
runSolution _ = do
	putStrLn "**Day 15**"
	print $ solver input 30000000
	where
	input :: [Input]
	--input = [0,3,6]
	input = [0,13,1,16,6,17]


text:
import Test.Framework (defaultMain)

import qualified Statistics.Iteratee.Tests as T

main :: IO ()
main = defaultMain T.tests

max_stars_repo_name: JohnLato/iter-stats

max_issues_repo_name: JohnLato/iter-stats

max_forks_repo_name: JohnLato/iter-stats


text:
-- example 3.22

import Statistics.Distribution
import Statistics.Distribution.Poisson

-- three accidents per month means six accidents per two-months:
pd :: PoissonDistribution
pd = poisson 6.0

p :: Int -> Double
p = probability pd

main = do
  putStrLn $ "the mean =   " ++ show (mean pd)
  putStrLn $ "the stdDev = " ++ show (stdDev pd)
  putStrLn $ "the probability of ten accidents in two months: " ++ show (p 10)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
{-# LANGUAGE DataKinds       #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE ViewPatterns    #-}

import           Data.Maybe
import           Numeric.Hamilton
import           Numeric.LinearAlgebra.Static
import qualified Data.Vector.Sized            as V

doublePendulum :: System 4 2
doublePendulum = mkSystem' masses coordinates potential
  where
    masses :: R 4
    masses = vec4 1 1 2 2
    coordinates
        :: Floating a
        => V.Vector 2 a
        -> V.Vector 4 a
    coordinates (V2 θ1 θ2) = V4 (sin θ1)            (-cos θ1)
                                (sin θ1 + sin θ2/2) (-cos θ1 - cos θ2/2)
    potential
        :: Num a
        => V.Vector 4 a
        -> a
    potential (V4 _ y1 _ y2) = (y1 + 2 * y2) * 5    -- assuming g = 5

pattern V2 :: a -> a -> V.Vector 2 a
pattern V2 x y <- (V.toList->[x,y])
  where
    V2 x y = fromJust (V.fromList [x,y])

pattern V4 :: a -> a -> a -> a -> V.Vector 4 a
pattern V4 x y z a <- (V.toList->[x,y,z,a])
  where
    V4 x y z a = fromJust (V.fromList [x,y,z,a])

config0 :: Config 2
config0 = Cfg (vec2 1 0  )  -- initial positions
              (vec2 0 0.5)  -- initial velocities

phase0 :: Phase 2
phase0 = toPhase doublePendulum config0

evolution :: [Phase 2]
evolution = evolveHam' doublePendulum phase0 [0,0.1 .. 1]

evolution' :: [Phase 2]
evolution' = iterate (stepHam 0.1 doublePendulum) phase0

positions :: [R 2]
positions = phsPositions <$> evolution'

positions' :: [R 4]
positions' = underlyingPos doublePendulum <$> positions

main :: IO ()
main = withRows (take 25 positions) (disp 4)

max_stars_repo_name: mstksg/inCode

max_issues_repo_name: mstksg/inCode

max_forks_repo_name: mstksg/inCode


text:
module European (
    European(European),
    PutCall(PUT, CALL),
    black76Call,
    black76Put) where

import Dates
import Market
import Instrument
import Statistics.Distribution
import Statistics.Distribution.Normal

-- |Enum defining a put (option to sell) or call (option to buy)
data PutCall = PUT | CALL deriving (Eq, Show, Read)

-- |The 1976 formulation of the Black-Scholes formula for pricing a European
-- |call option. Parameters discount d, forward f, variance v, strike k,
black76Call :: Double -> Double -> Double -> Double -> Double
black76Call d f v k = d * ((n d1) * f - (n d2) * k) where
    d1 = (log (f/k) + 0.5 * v) / sqrt v
    d2 = d1 - v
    n = cumulative standard

-- |Put call parity says that C - P = S - K, or that here
-- |P = C + d (k - f)
black76Put :: Double -> Double -> Double -> Double -> Double
black76Put d f v k = (black76Call d f v k) + d * (k - f)

-- |European option
data European = European {
    putOrCall :: PutCall,
    expiry :: Date,
    creditEntity :: String,      -- e.g. Currency for a risk-free option
    underlying :: String,
    underlyingEntity :: String,  -- e.g. Currency of the underlying asset
    strike :: Double } deriving (Eq, Show, Read)

-- |Price a European given a discount, forward and vol surface
priceGivenCurves :: European -> Date -> Discount -> Forward -> VolSurface -> Double
priceGivenCurves e valueDate d f v =
    let
        k = (strike e)
        date = (expiry e)
        dates = [date]   -- only one date of interest (ignoring settlement)
        discount = head (d dates)
        forward = head (f dates)
        vol = v f date k
        t = act365 valueDate date
        variance = vol^2 * t
    in
        case putOrCall e of
            PUT -> black76Put discount forward variance k
            CALL -> black76Call discount forward variance k

-- |Generic pricer for European, where it finds the data it needs in the market
instance Priceable European where
    price e m = do
        let ul = (underlying e)
        let valDate = valueDate m
        discount <- findDiscount m (creditEntity e)
        forward <- findEquityForward m ul (underlyingEntity e)
        volSurface <- findVol m ul
        return $ priceGivenCurves e valDate discount forward volSurface

max_stars_repo_name: MarcusRainbow/Ffinar

max_issues_repo_name: MarcusRainbow/Ffinar

max_forks_repo_name: MarcusRainbow/Ffinar


text:
module Widget.Percentile
  ( Percentile()
  , mkPercentile
  ) where

import Import
import qualified Format.Percentile as P
import qualified Format.FastTime   as D

import Control.Monad (when)
import Statistics.Distribution.Empirical

data Percentile = Percentile (Zipper EmpiricalDistribution) (Zipper NominalDiffTime) (Zipper NominalDiffTime) (Zipper Bool) Window

instance Widget Percentile where
  update                        = updatePercentile
  handle                        = handlePercentile
  redraw                        = redrawPercentile
  window (Percentile _ _ _ _ w) = w

mkPercentile :: [(EmpiricalDistribution, NominalDiffTime)] -> Window -> Percentile
mkPercentile p = Percentile (listToZipper $ map (\(ed, _) -> ed) p)
                            (listToZipper $ map (\(_, c) -> c) p)
                            (listToZipper [toEnum 0])
                            (listToZipper [True])

updatePercentile :: NominalDiffTime -> Percentile -> Curses Percentile
updatePercentile deltaT (Percentile z bs ts v w) =
  let ts' = replace (curs ts + deltaT) ts
      p   = Percentile z bs ts' v w
  in redrawPercentile p >> return p

handlePercentile :: TimerAction -> Percentile -> Curses Percentile
handlePercentile Advance (Percentile z bs ts v w) =
  return $ Percentile (next z) (next bs) (push 0 ts) (push True v) w

handlePercentile Reverse p@(Percentile z bs ts v w) =
  return $
    if null $ left z
      then p
      else
        let c   = curs ts
            ts' = trunc $ prev ts
            v'  = if (not $ curs $ prev $ prev v) && (not $ null $ left $ prev v)
                    then trunc $ prev v
                   else replace True $ trunc $ prev v
        in Percentile (prev z) (prev bs) (replace (curs ts' + c) ts') v' w

handlePercentile Skip    (Percentile z bs ts v w) =
  return $ Percentile (next z) (next bs) (push (toEnum 0) ts) (push False $ replace False v) w

handlePercentile Reset   (Percentile z bs ts v w) =
  return $ mkPercentile (zipWith (\ed c -> (ed, c)) (zipperToList z) (zipperToList bs)) w

handlePercentile ValidInvalid (Percentile z bs ts v w) =
  return $ Percentile z bs ts (replace (not $ curs v) v) w

handlePercentile _ t = return t

redrawPercentile :: Percentile -> Curses ()
redrawPercentile p@(Percentile z bs t v w) = do
  cid <- newColorID ColorBlack ColorDefault 1

  updateWindow w $ do
    (rows, columns) <- windowSize

    erase

    moveCursor 0 0
    drawLineV Nothing rows

    let (lt, ct, rt) = sizeZipperToScreen rows t
        (lb, cb, rb) = sizeZipperToScreen rows bs
        (lv, cv, rv) = sizeZipperToScreen rows v
        (lz, cz, rz) = sizeZipperToScreen rows z

        t'  = (lt ++ ct:rt)
        bs' = (lb ++ cb:rb)
        v'  = (lv ++ cv:rv)
        distrs = (lz ++ cz:rz)

        percentileF distr best current =
          if current == 0
            then "----"
            else if current >= best && best > 0
                   then show $ P.percentile distr current
                   else show $ D.fastTime $ current - best
        percentileFs = zipWith3 percentileF distrs bs' t'

    mapM_ (go cid)
      $ zip3 [0..] v' percentileFs

  where go cid (a, valid, p) = do
          when valid $ do
            moveCursor a 1
            drawString p

max_stars_repo_name: argiopetech/timer

max_issues_repo_name: argiopetech/timer

max_forks_repo_name: argiopetech/timer


text:
module Main where

import qualified Control.Monad.Random as RndM
import qualified Numeric.LinearAlgebra as LA
import MachineLearning.Clustering
import qualified Data.Vector as V

main = do
  putStrLn "\n=== K-Means Sample Application ===\n"
  -- Step 1. Data loading.
  -- Step 1.1 Training Data loading.
  x <- LA.loadMatrix "linear_regression/data.txt"
  -- Step 2. Learning
  let numIters = 25
      numClusters = 3
  clusterList <- RndM.evalRandIO $ kmeans numIters x numClusters

  -- Step 3. Display results
  print $ "Successfully groupped into " ++ show (V.length clusterList) ++ " clusters"

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Quadratic_Formula (csqrt, quadratic) where

import Data.Complex

csqrt z@(a :+ b) = let (mag, ph) = polar z in mkPolar (sqrt mag) (ph/2)

quadratic a b c = (
    quad (conjugate a) (conjugate b) (conjugate c) (+),
    quad (conjugate a) (conjugate b) (conjugate c) (-)
    ) where quad a b c op = (-b `op` csqrt ( b^2 - 4 * a * c))/(2 * a)


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

module Foreign.BLAS.Level3 (
    BLAS3(..),
    ) where

import Data.Complex
import Foreign( Ptr, Storable, with )

import Foreign.BLAS.Types
import Foreign.BLAS.Level2
import Foreign.BLAS.Double
import Foreign.BLAS.Zomplex

-- | Types with matrix-matrix operations.
class (BLAS2 a) => BLAS3 a where
    gemm  :: Trans -> Trans -> Int -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    symm  :: Side -> Uplo -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    hemm  :: Side -> Uplo -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    trmm  :: Side -> Uplo -> Trans -> Diag -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    trsm  :: Side -> Uplo -> Trans -> Diag -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> IO ()
    syrk  :: Uplo -> Trans -> Int -> Int -> a -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    syr2k :: Uplo -> Trans -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    herk  :: Uplo -> Trans -> Int -> Int -> a -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
    her2k :: Uplo -> Trans -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()

withEnum :: (Enum a, Storable a) => Int -> (Ptr a -> IO b) -> IO b
withEnum = with . toEnum
{-# INLINE withEnum #-}

instance BLAS3 Double where
    gemm transa transb m n k alpha pa lda pb ldb beta pc ldc =
        withTrans transa $ \ptransa ->
        withTrans transb $ \ptransb ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            dgemm ptransa ptransb pm pn pk palpha pa plda pb pldb pbeta pc pldc
    {-# INLINE gemm #-}

    symm side uplo m n alpha pa lda pb ldb beta pc ldc =
        withSide side $ \pside ->
        withUplo uplo $ \puplo ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            dsymm pside puplo pm pn palpha pa plda pb pldb pbeta pc pldc
    {-# INLINE symm #-}

    hemm = symm
    {-# INLINE hemm #-}

    trmm side uplo transa diag m n alpha pa lda pb ldb =
        withSide side $ \pside ->
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withDiag diag $ \pdiag ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
            dtrmm pside puplo ptransa pdiag pm pn palpha pa plda pb pldb
    {-# INLINE trmm #-}

    trsm side uplo transa diag m n alpha pa lda pb ldb =
        withSide side $ \pside ->
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withDiag diag $ \pdiag ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
            dtrsm pside puplo ptransa pdiag pm pn palpha pa plda pb pldb
    {-# INLINE trsm #-}

    syrk uplo transa n k alpha pa lda beta pc ldc =
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            dsyrk puplo ptransa pn pk palpha pa plda pbeta pc pldc
    {-# INLINE syrk #-}

    syr2k uplo transa n k alpha pa lda pb ldb beta pc ldc =
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            dsyr2k puplo ptransa pn pk palpha pa plda pb pldb pbeta pc pldc
    {-# INLINE syr2k #-}

    herk = syrk
    {-# INLINE herk #-}

    her2k = syr2k
    {-# INLINE her2k #-}

instance BLAS3 (Complex Double) where
    gemm transa transb m n k alpha pa lda pb ldb beta pc ldc =
        withTrans transa $ \ptransa ->
        withTrans transb $ \ptransb ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            zgemm ptransa ptransb pm pn pk palpha pa plda pb pldb pbeta pc pldc
    {-# INLINE gemm #-}

    symm side uplo m n alpha pa lda pb ldb beta pc ldc =
        withSide side $ \pside ->
        withUplo uplo $ \puplo ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            zsymm pside puplo pm pn palpha pa plda pb pldb pbeta pc pldc
    {-# INLINE symm #-}

    hemm side uplo m n alpha pa lda pb ldb beta pc ldc =
        withSide side $ \pside ->
        withUplo uplo $ \puplo ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            zhemm pside puplo pm pn palpha pa plda pb pldb pbeta pc pldc
    {-# INLINE hemm #-}

    trmm side uplo transa diag m n alpha pa lda pb ldb =
        withSide side $ \pside ->
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withDiag diag $ \pdiag ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
            ztrmm pside puplo ptransa pdiag pm pn palpha pa plda pb pldb
    {-# INLINE trmm #-}

    trsm side uplo transa diag m n alpha pa lda pb ldb =
        withSide side $ \pside ->
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withDiag diag $ \pdiag ->
        withEnum m $ \pm ->
        withEnum n $ \pn ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
            ztrsm pside puplo ptransa pdiag pm pn palpha pa plda pb pldb
    {-# INLINE trsm #-}

    syrk uplo transa n k alpha pa lda beta pc ldc =
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            zsyrk puplo ptransa pn pk palpha pa plda pbeta pc pldc
    {-# INLINE syrk #-}

    syr2k uplo transa n k alpha pa lda pb ldb beta pc ldc =
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            zsyr2k puplo ptransa pn pk palpha pa plda pb pldb pbeta pc pldc
    {-# INLINE syr2k #-}

    herk uplo transa n k alpha pa lda beta pc ldc =
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            zherk puplo ptransa pn pk palpha pa plda pbeta pc pldc
    {-# INLINE herk #-}

    her2k uplo transa n k alpha pa lda pb ldb beta pc ldc =
        withUplo uplo $ \puplo ->
        withTrans transa $ \ptransa ->
        withEnum n $ \pn ->
        withEnum k $ \pk ->
        with alpha $ \palpha ->
        withEnum lda $ \plda ->
        withEnum ldb $ \pldb ->
        with beta $ \pbeta ->
        withEnum ldc $ \pldc ->
            zher2k puplo ptransa pn pk palpha pa plda pb pldb pbeta pc pldc
    {-# INLINE her2k #-}

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# OPTIONS_GHC -Wno-incomplete-patterns #-}

module UnboxedSpec where

import Data.Complex
import Data.Proxy

import Test.QuickCheck

import Category
import Functor
import Hask()
import Unboxed

type UA = Int *#* Float
type UB = Double
type UC = Complex Double
type UD = Int

prop_Unboxed_Discretization_inv :: Fun UA UB -> FnProp UA
prop_Unboxed_Discretization_inv f =
    checkFnEqual (law_Discretization_inv (UQCFun f))

prop_Unboxed_Discretization_approx :: Fun UA UB -> FnProp UA
prop_Unboxed_Discretization_approx (Fn f) =
    checkFnEqual (law_Discretization_approx (Proxy @(-#>)) f)

prop_Unboxed_MId_id :: UA -> Property
prop_Unboxed_MId_id = law_MId_id (Proxy @Unboxed)

prop_Unboxed_MCompose_compose :: Fun UB UC -> Fun UA UB -> UA -> Property
prop_Unboxed_MCompose_compose f g = law_MCompose_compose (UQCFun f) (UQCFun g)

prop_Unboxed_MCompose_left_id :: Fun UA UB -> UA -> Property
prop_Unboxed_MCompose_left_id f = law_MCompose_left_id (UQCFun f)

prop_Unboxed_MCompose_right_id :: Fun UA UB -> UA -> Property
prop_Unboxed_MCompose_right_id f = law_MCompose_right_id (UQCFun f)

prop_Unboxed_MCompose_assoc ::
    Fun UC UD -> Fun UB UC -> Fun UA UB -> UA -> Property
prop_Unboxed_MCompose_assoc f g h =
    law_MCompose_assoc (UQCFun f) (UQCFun g) (UQCFun h)

prop_Unboxed_Functor_id :: FnProp (Float *#* UA)
prop_Unboxed_Functor_id = checkFnEqual law_Functor_id

prop_Unboxed_Functor_comp ::
    (UB -#> UC) -> (UA -#> UB) -> FnProp (Float *#* UA)
prop_Unboxed_Functor_comp f g = checkFnEqual (law_Functor_comp f g)


text:
{-# LANGUAGE DataKinds        #-}
{-# LANGUAGE TypeApplications #-}

module Main ( main ) where

import           Criterion.Main
import           Data.List.NonEmpty (NonEmpty(..))
import           Data.Massiv.Array as A hiding (zipWith)
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U
import           Geography.MapAlgebra
import           GHC.TypeLits
import           Graphics.ColorSpace
import qualified Numeric.LinearAlgebra as LA
import           Prelude hiding (zipWith)

---

main :: IO ()
main = do
  img   <- fileY "data/gray512.tif"
  imgF  <- fileY "data/gray512.tif"
  rgba  <- fileRGB "data/512x512.tif"
  a512  <- fileRGB "data/512x512.tif"
  -- a1024 <- fileRGB "data/1024x1024.tif"
  -- a2048 <- fileRGB "data/2048x2048.tif"
  -- a4096 <- fileRGB "data/4096x4096.tif"
  -- let aHuge :: Raster S u 46500 46500 Word8
  --     aHuge = fromFunction S Par (\(r :. c) -> fromIntegral $ r * c)
  defaultMain
    [
      creation
    , io
    , colouring img
    , massivOps img
    , localOps rgba a512 img
    , hmatrix
    , conversions img
    , focalOps img imgF
    , compositeOps a512 -- a1024 a2048 a4096
    , histograms img
    ]

creation :: Benchmark
creation = bgroup "Raster Creation"
           [ bench "constant 256x256"     $ nf (_array . constantB) 5
           , bench "constant 512x512"     $ nf (_array . constantB') 5
           , bench "fromFunction 256x256" $ nf (_array . functionB) (\(r :. c) -> r * c)
           , bench "fromFunction 512x512" $ nf (_array . functionB') (\(r :. c) -> r * c)
           , bench "fromVector Unboxed Int - 256x256" $ nf (_array . vectorB) uv
           , bench "fromVector Boxed Int - 256x256"   $ nf (_array . vectorB') bv ]
  where uv = U.fromList ([1..65536] :: [Int])
        bv = V.fromList ([1..65536] :: [Int])

io :: Benchmark
io = bgroup "IO"
     [ bench "fromRGBA 512x512" $ nfIO (_array . _red <$> fileRGB @Word8 @512 @512 "data/512x512.tif")
     , bench "fromRGBA (Word8 -> Double) 512x512" $ nfIO (_array . _red <$> fileRGB @Double @512 @512 "data/512x512.tif")
     , bench "fromGray Multiband 512x512"  $ nfIO (_array <$> fileY @Word8 @512 @512 "data/512x512.tif")
     , bench "fromGray Singleband 512x512" $ nfIO (_array <$> gray "data/gray512.tif") ]

colouring :: Raster S p 512 512 Word8 -> Benchmark
colouring img = bgroup "Colouring"
                [ bench "classify 512x512"  $ nf (_array . strict S . classify invisible cr . lazy) img
                , bench "grayscale 512x512" $ nf (_array . strict S . grayscale . lazy) img ]
  where cr = greenRed [25, 50, 75, 100, 125, 150, 175, 200, 225, 255]

massivOps :: Raster S p 512 512 Word8 -> Benchmark
massivOps img = bgroup "Massiv Operations"
                [ bench "strict S . lazy" $ nf (_array . strict S . lazy) img
                , bench "strict U . lazy" $ nf (_array . strict U . lazy) img
                , bench "strict P . lazy" $ nf (_array . strict P . lazy) img ]

localOps :: RGBARaster p 512 512 Word8 -> RGBARaster p 512 512 Double -> Raster S p 512 512 Word8 -> Benchmark
localOps (RGBARaster r g b _) (RGBARaster rF gF bF _) img = bgroup "Local Operations"
  [ bench "fmap (+ 17) . lazy" $ nf (_array . strict S . fmap (+ 17) . lazy) img
  , bench "zipWith (+)"    $ nf (_array . strict S . zipWith (+) r) g
  , bench "zipWith (/)"    $ nf (_array . strict S . zipWith (/) rF) gF
  , bench "(+)"            $ nf (_array . strict S . (+ lazy r)) (lazy g)
  , bench "(/)"            $ nf (_array . strict S . (/ lazy rF)) (lazy gF)
  -- , bench "(.+)"           $ nf (\g' -> computeAs S $ _array r .+ g') (_array g)
  -- , bench "(./)"           $ nf (\g' -> computeAs S $ _array rF ./ g') (_array gF)
  , bench "lmax"           $ nf (_array . strict S . lmax img) img
  , bench "lmin"           $ nf (_array . strict S . lmin img) img
  , bench "lmean (Word8)"  $ nf (_array . strict S . lmean @Word8 @Double) rs
  , bench "lmean (Double)" $ nf (_array . strict S . lmean @Double @Double) rsF
  , bench "lvariety"       $ nf (_array . strict S . lvariety) rs
  , bench "lmajority"      $ nf (_array . strict S . lmajority) rs
  , bench "lminority"      $ nf (_array . strict S . lminority) rs
  , bench "lvariance (Word8)"  $ nf (fmap (_array . strict S) . lvariance) rs
  , bench "lvariance (Double)" $ nf (fmap (_array . strict S) . lvariance) rsF ]
  where rs  = lazy r :| [lazy g, lazy b]
        rsF = lazy rF :| [lazy gF, lazy bF]

hmatrix :: Benchmark
hmatrix = bgroup "HMatrix"
          [ bench "linearSolveLS" $ nf (LA.linearSolveLS zing) (LA.col [8,8,8,8,8,8,8,8,8])
          , bench "manual - MxM"  $ nf (leftPseudo <>) (LA.col [8,8,8,8,8,8,8,8,8])
          , bench "manual - MxV"  $ nf (leftPseudo LA.#>) (LA.vector [8,8,8,8,8,8,8,8,8]) ]

conversions :: Raster S p 512 512 Word8 -> Benchmark
conversions img = bgroup "Numeric Conversion"
                  [ bench "Double -> Double via id"         $ nf id tau
                  , bench "Double -> Double via realToFrac" $ nf (realToFrac @Double @Double) tau
                  , bench "Word -> Double via realToFrac"   $ nf (realToFrac @Word8 @Double) 5
                  , bench "realToFrac on Raster"            $ nf (_array . strict S . fmap (realToFrac @Word8 @Double) . lazy) img ]

focalOps :: Raster S p 512 512 Word8 -> Raster S p 512 512 Double -> Benchmark
focalOps img imgF = bgroup "Focal Operations"
               [ bgroup "fsum"
                 [ bench "512"    $ nf (_array . strict S . fsum) img
                 -- , bench "46500"  $ nf (_array . strict S . fsum) huge
                 -- , bench "maybing" $ nf (_array . strict B . fmap Just . lazy) img
                 -- , bench "nodata" $ nf (_array . strict S . nodatafsum) img
                 ]
               , bgroup "fmean"
                 [ bench "Word8"  $ nf (_array . strict S . fmean @Word8 @Double) img
                 , bench "Double" $ nf (_array . strict S . fmean @Double @Double) imgF
                 ]
               , bench "fmax"        $ nf (_array . strict S . fmax) img
               , bench "fmin"        $ nf (_array . strict S . fmin) img
               , bench "fmajority"   $ nf (_array . strict S . fmajority) img
               , bench "fminority"   $ nf (_array . strict S . fminority) img
               , bench "fvariety"    $ nf (_array . strict S . fvariety) img
               , bench "fpercentage" $ nf (_array . strict S . fpercentage) img
               , bench "fpercentile" $ nf (_array . strict S . fpercentile) img
               , bench "flinkage"    $ nf (_array . strict S . flinkage) img
               , bench "flength"     $ nf (_array . strict S . flength . flinkage) img
               , bench "fpartition"  $ nf (_array . strict B . fpartition) img
               , bench "fshape"      $ nf (_array . strict B . fshape) img
               , bench "ffrontage"   $ nf (_array . strict S . ffrontage . fshape) img
               , bench "farea"       $ nf (_array . strict S . farea . fshape) img
               , bgroup "fvolume"
                 [ bench "Word8 -> Double" $ nf (_array . strict S . fvolume . strict S . fmap (realToFrac @Word8 @Double) . lazy) img
                 , bench "Double" $ nf (_array . strict S . fvolume) imgF
                 ]
               , bgroup "fgradient"
                 [ bench "Word8"  $ nf (_array . strict S . fgradient . strict S . fmap wtod . lazy) img
                 , bench "Double" $ nf (_array . strict S . fgradient) imgF
                 ]
               , bgroup "faspect"
                 [ bench "Unsafe (Word8)"  $ nf (_array . strict S . faspect' . strict S . fmap wtod . lazy) img
                 , bench "Unsafe (Double)" $ nf (_array .strict S . faspect') imgF
                 , bench "Safe (Word8)"    $ nf (_array . strict B . faspect . strict S . fmap wtod . lazy) img
                 , bench "Safe (Double)"   $ nf (_array . strict B . faspect) imgF
                 ]
               , bgroup "fdownstream"
                 [ bench "Word8"  $ nf (_array . strict S . fdownstream . strict S . fmap wtod . lazy) img
                 , bench "Double" $ nf (_array. strict S . fdownstream) imgF
                 ]
               , bgroup "fupstream"
                 [ bench "Word8"  $ nf (_array . strict S . fupstream . strict S . fdownstream . strict S . fmap wtod . lazy) img
                 , bench "Double" $ nf (_array . strict S . fupstream . strict S . fdownstream) imgF
                 ]
               ]

wtod :: Word8 -> Double
wtod = realToFrac

compositeOps :: RGBARaster p 512 512 Double
  -- -> RGBARaster p 1024 1024 Double
  -- -> RGBARaster p 2048 2048 Double
  -- -> RGBARaster p 4096 4096 Double
  -> Benchmark
compositeOps i@(RGBARaster r g _ _) =
  bgroup "Composite Operations"
  [ bench "NDVI" $ nf (_array . strict S . ndvi (lazy g)) (lazy r)
  , bench "EVI (512)"    $ nf (_array . strict S . evi) i
  -- , bench "EVI (1024)"   $ nf (_array . strict S . evi) a1024
  -- , bench "EVI (2048)"   $ nf (_array . strict S . evi) a2048
  -- , bench "EVI (4096)"   $ nf (_array . strict S . evi) a4096
  , bench "EVI + Colour" $ nf (_array . strict S . classify invisible cr . evi) i
  , bench "EVI + Colour + PNG (D)" $ nf (png . classify invisible cr . evi) i
  , bench "EVI + Colour + PNG (S)" $ nf (png . strict S . classify invisible cr . evi) i ]
  where cr = greenRed $ fmap (10 ^) ([1..10] :: [Int])

fromRight :: Either a b -> b
fromRight (Right b) = b
fromRight _         = error "Was Left"

constantB :: Int -> Raster S p 256 256 Int
constantB = constant S Par

constantB' :: Int -> Raster S p 512 512 Int
constantB' = constant S Par

functionB :: (Ix2 -> Int) -> Raster S p 256 256 Int
functionB = fromFunction S Par

functionB' :: (Ix2 -> Int) -> Raster S p 512 512 Int
functionB' = fromFunction S Par

vectorB :: U.Vector Int -> Raster U p 256 256 Int
vectorB = fromRight . fromVector Par

vectorB' :: V.Vector Int -> Raster B p 256 256 Int
vectorB' = fromRight . fromVector Par

-- doubles :: Raster U p 512 512 Double
-- doubles = fromRight . fromVector Par $ U.fromList ([1..262144] :: [Double])

zing :: LA.Matrix Double
zing = LA.matrix 3 [ -0.5, -0.5, 1
                   , -0.5, 0, 1
                   , -0.5, 0.5, 1
                   , 0, -0.5, 1
                   , 0, 0, 1
                   , 0, 0.5, 1
                   , 0.5, -0.5, 1
                   , 0.5, 0, 1
                   , 0.5, 0.5, 1 ]

gray :: FilePath -> IO (Raster S p 512 512 Word8)
gray fp = do
  i <- fromGray fp
  case i of
    Left err  -> error err
    Right img -> pure img

fileY :: (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (Raster S p r c a)
fileY fp = do
  i <- fromGray fp
  case i of
    Left err  -> error err
    Right img -> pure img
{-# INLINE fileY #-}

fileRGB :: (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (RGBARaster p r c a)
fileRGB fp = do
  i <- fromRGBA fp
  case i of
    Left err  -> error err
    Right img -> pure img
{-# INLINE fileRGB #-}

-- | See: https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index
ndvi :: (KnownNat r, KnownNat c) => Raster D p r c Double -> Raster D p r c Double -> Raster D p r c Double
ndvi nir red = (nir - red) / (nir + red)
{-# INLINE ndvi #-}

-- | See: https://en.wikipedia.org/wiki/Enhanced_vegetation_index
evi :: (Fractional a, Storable a, KnownNat r, KnownNat c) => RGBARaster p r c a -> Raster D p r c a
evi (RGBARaster r g b _) = 2.5 * (numer / denom)
  where nir   = lazy g  -- fudging it.
        numer = nir - lazy r
        denom = nir + (6 * lazy r) - (7.5 * lazy b) + 1
{-# INLINE evi #-}

-- nodatafsum :: Raster S p 1753 1760 Double -> Raster DW p 1753 1760 Double
-- nodatafsum = fmap (maybe 0 ((/9) . getSum)) . fmonoid . strict B . fmap check . lazy
--   where check 0 = Nothing
--         check n = Just $ Sum n
-- {-# INLINE nodatafsum #-}

histograms :: Raster S p 512 512 Word8 -> Benchmark
histograms r = bgroup "Histogram Algorithms"
  [ bench "Storable Super Mutable" $ nf (_histogram . histogram) r
  ]

max_stars_repo_name: fosskers/mapalgebra

max_issues_repo_name: fosskers/mapalgebra

max_forks_repo_name: fosskers/mapalgebra


text:
{-# LANGUAGE FlexibleContexts #-}

module Numeric.HHT.Internal.FFT (
    fft
  , ifft
  ) where

import           Data.Complex
import qualified Data.Array.CArray         as CA
import qualified Data.Array.IArray         as IA
import qualified Data.Ix                   as Ix
import qualified Data.Vector.Generic       as VG
import qualified Data.Vector.Generic.Sized as SVG
import qualified Foreign.Storable          as FS
import qualified Math.FFT                  as FFT
import qualified Math.FFT.Base             as FFT

fft :: (FFT.FFTWReal a, VG.Vector v (Complex a))
    => SVG.Vector v n (Complex a)
    -> SVG.Vector v n (Complex a)
fft = SVG.withVectorUnsafe $
        fromCA
      . FFT.dft
      . toCA

ifft
    :: (FFT.FFTWReal a, VG.Vector v (Complex a))
    => SVG.Vector v n (Complex a)
    -> SVG.Vector v n (Complex a)
ifft = SVG.withVectorUnsafe $
        fromCA
      . FFT.idft
      . toCA

fromCA
    :: (FS.Storable a, VG.Vector v (Complex a))
    => CA.CArray Int (Complex a)
    -> v (Complex a)
fromCA v = VG.generate (Ix.rangeSize (IA.bounds v)) (v IA.!)

toCA
    :: (FS.Storable a, VG.Vector v (Complex a))
    => v (Complex a)
    -> CA.CArray Int (Complex a)
toCA v = IA.listArray (0, VG.length v - 1) (VG.toList v)

max_stars_repo_name: mstksg/emd

max_issues_repo_name: mstksg/emd

max_forks_repo_name: mstksg/emd


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Pendulum.Model (
  Pendulum(..),
  mkPendlumSystem,
  stepPendulum
  ) where

import qualified Data.Vector.Sized as V
import           GHC.TypeLits (KnownNat)
import           Numeric.Hamilton
import           Numeric.LinearAlgebra.Static

data Pendulum n = Pendulum
  { _lengthes :: V.Vector n Double
  , _masses   :: V.Vector n Double
  , _config   :: Config n
  }

gravity :: Floating a => a
gravity = 9.8

toFloating :: Floating a => Double -> a
toFloating = fromRational . toRational

duplicate :: [a] -> [a]
duplicate = concatMap (\x -> [x, x])

merge :: [a] -> [a] -> [a]
merge []     ys     = ys
merge xs     []     = xs
merge (x:xs) (y:ys) = x : y : merge xs ys

evenElems :: [a] -> [a]
evenElems []       = []
evenElems (_:[])   = []
evenElems (_:x:xs) = x : evenElems xs

toVector :: forall a m. KnownNat m => [a] -> V.Vector m a
toVector = (\(Just x) -> x) . V.fromList

mkPendlumSystem :: forall m n. (KnownNat m, KnownNat n) => Pendulum n -> System m n
mkPendlumSystem p =
  let lengthes = V.toList (_lengthes p)
      masses   = V.toList (_masses p)
      masses' :: R m
      masses' = vector (duplicate masses)
      coordinates :: Floating a => V.Vector n a -> V.Vector m a
      coordinates v =
        let ths = V.toList v
            ls = map toFloating lengthes
            xs =                scanl1 (+) (zipWith (*) ls (map sin ths))
            ys = (* (-1.0)) <$> scanl1 (+) (zipWith (*) ls (map cos ths))
         in toVector (merge xs ys)
      potential :: Floating a => V.Vector m a -> a
      potential v =
        let ys = evenElems (V.toList v)
            ms = map toFloating masses
         in gravity * foldl1 (+) (zipWith (*) ms ys)
   in mkSystem' masses' coordinates potential

stepPendulum :: forall m n. (KnownNat m, KnownNat n) => Double -> System m n -> Pendulum n -> Pendulum n
stepPendulum dt system pd =
  let phase  = toPhase system (_config pd)
      phase' = stepHam dt system phase
   in pd {_config = fromPhase system phase'}


text:
-- | Internal utilities.

module Math.Probably.Utils where

import Data.Map (Map)
import qualified Data.Map as Map
import Data.Maybe
import qualified Data.Vector.Storable as V
import Math.Probably.Types
import Numeric.LinearAlgebra
import Statistics.Distribution
import Statistics.Distribution.Normal

lookupDefault  :: Ord k => a -> k -> Map k a -> a
lookupDefault d k m = fromMaybe d (Map.lookup k m)

-- | A spherical Gaussian distribution.
sphereGauss :: ContinuousParams -> ContinuousParams -> Double -> Double
sphereGauss xs m sd = product $ zipWith density normalDists xsAsList where
  xsAsList    = toList xs
  meanAsList  = toList m
  normalDists = map (`normalDistr` sd) meanAsList

-- | Scalar-vector multiplication.
(.*) :: Double -> ContinuousParams -> ContinuousParams
z .* xs = mapVector (* z) xs

-- | Scalar-vector subtraction.
(.-) :: ContinuousParams -> ContinuousParams -> ContinuousParams
xs .- ys = zipVectorWith (-) xs ys

-- | Scalar-vector addition.
(.+) :: ContinuousParams -> ContinuousParams -> ContinuousParams
xs .+ ys = zipVectorWith (+) xs ys

-- | The leapfrog integrator.
leapfrog :: Gradient -> Particle -> Double -> Particle
leapfrog glTarget (q, r) e = (qf, rf) where
  rm = adjustMomentum glTarget e (q, r)
  qf = adjustPosition e (rm, q)
  rf = adjustMomentum glTarget e (qf, rm)

-- | Adjust momentum according to a half-leapfrog step.
adjustMomentum :: Gradient -> Double -> Particle -> ContinuousParams
adjustMomentum glTarget e (t, r) = r .+ ((e / 2) .* glTarget t)

-- | Adjust position according to a half-leapfrog step.
adjustPosition :: Double -> Particle -> ContinuousParams
adjustPosition e (r, t) = t .+ (e .* r)

-- | A target augmented by momentum auxilliary variables.
auxilliaryTarget :: (ContinuousParams -> Double) -> Particle -> Double
auxilliaryTarget lTarget (t, r) =
  lTarget t - 0.5 * innerProduct r r

innerProduct :: ContinuousParams -> ContinuousParams -> Double
innerProduct xs ys = V.sum $ V.zipWith (*) xs ys

indicate :: Integral a => Bool -> a
indicate True  = 1
indicate False = 0

fi :: (Integral a, Num b) => a -> b
fi = fromIntegral

max_stars_repo_name: glutamate/probably-baysig

max_issues_repo_name: silky/probably-baysig

max_forks_repo_name: silky/probably-baysig


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

-- | Marshaling between Feldspar and C99 types
--
module Feldspar.Compiler.Marshal
  ( SA(..)
  , allocSA
  , Marshal(..)
  )
  where

import System.Plugins.MultiStage
import Feldspar.Core.Types (IntN(..), WordN(..), Tuple(..), (:*), TNil)
import Feldspar.Core.NestedTuples (npair)

import Data.Int (Int32)
import Data.Word (Word32)
import Data.Complex (Complex(..))
import Data.Default
import Control.Applicative

import Foreign.Ptr
import Foreign.Marshal (free, new, newArray, peekArray,mallocBytes)
import Foreign.Storable (Storable(..))
import Foreign.Storable.Tuple ()
import qualified Foreign.Storable.Record as Store

instance Reference IntN        where type Ref IntN        = IntN
instance Reference WordN       where type Ref WordN       = WordN
instance Reference (Complex a) where type Ref (Complex a) = Complex a

instance Marshal IntN        where type Rep IntN        = IntN
instance Marshal WordN       where type Rep WordN       = WordN
instance Marshal (Complex a) where type Rep (Complex a) = Complex a

instance Default (Ptr a) where def = nullPtr

instance (Storable (Rep a), Marshal a) => Marshal [a]
  where
    type Rep [a] = SA (Rep a)
    to xs = do
        let len  = fromIntegral $ length xs
        let size = fromIntegral $ sizeOf (undefined :: Rep a)
        ys <- mapM to xs
        buffer <- newArray ys
        return $ SA buffer len
    from p | elems p == 0 = return []
    from p = go p
      where
        go SA{..} = do
          res <- mapM from =<< peekArray (fromIntegral elems) buf
          free buf
          return res

-- | Buffer descriptor for Feldspar arrays
data SA a = SA { buf   :: Ptr a
               , elems :: Int32
               }
  deriving (Eq, Show)

instance Default (SA a) where
    def = SA nullPtr def

allocSA :: forall a. Storable a => Int -> IO (Ptr (SA a))
allocSA len = do
    let size  = fromIntegral $ sizeOf (undefined :: a)
    let bytes = len * size
    buffer <- mallocBytes bytes
    new $ SA buffer (fromIntegral len)

storeSA :: Storable a => Store.Dictionary (SA a)
storeSA = Store.run $ SA
    <$> Store.element buf
    <*> Store.element elems

instance Storable a => Storable (SA a)
  where
    sizeOf    = Store.sizeOf    storeSA
    alignment = Store.alignment storeSA
    peek      = Store.peek      storeSA
    poke      = Store.poke      storeSA

instance Reference (Ptr a)
  where
    type Ref (Ptr a) = Ptr a
    ref   = return
    deref = return

instance (Storable a) => Reference (SA a)
  where
    type Ref (SA a) = Ptr (SA a)
    ref   = new
    deref = peek

instance (Storable (a,b)) => Reference (a,b)
  where
    type Ref (a,b) = Ptr (a,b)
    ref   = new
    deref = peek

instance (Storable (a,b,c)) => Reference (a,b,c)
  where
    type Ref (a,b,c) = Ptr (a,b,c)
    ref   = new
    deref = peek

instance (Storable (a,b,c,d)) => Reference (a,b,c,d)
  where
    type Ref (a,b,c,d) = Ptr (a,b,c,d)
    ref   = new
    deref = peek

instance (Storable (a,b,c,d,e)) => Reference (a,b,c,d,e)
  where
    type Ref (a,b,c,d,e) = Ptr (a,b,c,d,e)
    ref   = new
    deref = peek

instance (Storable (a,b,c,d,e,f)) => Reference (a,b,c,d,e,f)
  where
    type Ref (a,b,c,d,e,f) = Ptr (a,b,c,d,e,f)
    ref   = new
    deref = peek

instance (Storable (a,b,c,d,e,f,g)) => Reference (a,b,c,d,e,f,g)
  where
    type Ref (a,b,c,d,e,f,g) = Ptr (a,b,c,d,e,f,g)
    ref   = new
    deref = peek

instance ( Marshal a
         , Marshal b
         ) => Marshal (a,b)
  where
    type Rep (a,b) = (Rep a,Rep b)
    to (a,b)   = (,) <$> to a <*> to b
    from (a,b) = (,) <$> from a <*> from b

instance ( Marshal a
         , Marshal b
         , Marshal c
         ) => Marshal (a,b,c)
  where
    type Rep (a,b,c) = (Rep a,Rep b,Rep c)
    to (a,b,c)   = (,,) <$> to a <*> to b <*> to c
    from (a,b,c) = (,,) <$> from a <*> from b <*> from c

instance ( Marshal a
         , Marshal b
         , Marshal c
         , Marshal d
         ) => Marshal (a,b,c,d)
  where
    type Rep (a,b,c,d) = (Rep a,Rep b,Rep c,Rep d)
    to (a,b,c,d) =
      (,,,) <$> to a <*> to b <*> to c <*> to d
    from (a,b,c,d) =
      (,,,) <$> from a <*> from b <*> from c <*> from d

instance ( Marshal a
         , Marshal b
         , Marshal c
         , Marshal d
         , Marshal e
         ) => Marshal (a,b,c,d,e)
  where
    type Rep (a,b,c,d,e) = (Rep a,Rep b,Rep c,Rep d,Rep e)
    to (a,b,c,d,e) =
      (,,,,) <$> to a <*> to b <*> to c <*> to d <*> to e
    from (a,b,c,d,e) =
      (,,,,) <$> from a <*> from b <*> from c <*> from d <*> from e

instance ( Marshal a
         , Marshal b
         , Marshal c
         , Marshal d
         , Marshal e
         , Marshal f
         ) => Marshal (a,b,c,d,e,f)
  where
    type Rep (a,b,c,d,e,f) = (Rep a,Rep b,Rep c,Rep d,Rep e,Rep f)
    to (a,b,c,d,e,f) =
      (,,,,,) <$> to a <*> to b <*> to c <*> to d <*> to e <*> to f
    from (a,b,c,d,e,f) =
      (,,,,,) <$> from a <*> from b <*> from c <*> from d <*> from e <*> from f

instance ( Marshal a
         , Marshal b
         , Marshal c
         , Marshal d
         , Marshal e
         , Marshal f
         , Marshal g
         ) => Marshal (a,b,c,d,e,f,g)
  where
    type Rep (a,b,c,d,e,f,g) = (Rep a,Rep b,Rep c,Rep d,Rep e,Rep f,Rep g)
    to (a,b,c,d,e,f,g) =
      (,,,,,,) <$> to a <*> to b <*> to c <*> to d <*> to e <*> to f <*> to g
    from (a,b,c,d,e,f,g) =
      (,,,,,,) <$> from a <*> from b <*> from c <*> from d <*> from e <*> from f <*> from g

instance (Marshal a, Marshal b) => Marshal (Tuple (a :* b :* TNil)) where
  type Rep (Tuple (a :* b :* TNil)) = (Rep a, Rep b)
  to (a :* b :* TNil) = (,) <$> to a <*> to b
  from (a, b) = npair <$> from a <*> from b

max_stars_repo_name: Feldspar/feldspar-compiler

max_issues_repo_name: Feldspar/feldspar-compiler

max_forks_repo_name: Feldspar/feldspar-compiler


text:
{- |
Description :  Convenience functions to simulate Gamma rate heterogeneity
Copyright   :  (c) Dominik Schrempf 2017
License     :  GPLv3

Maintainer  :  dominik.schrempf@gmail.com
Stability   :  unstable
Portability :  non-portable (not tested)

Sites evolve at different speeds. One way to model this phenomenon is to use
Gamma rate heterogeneity. This module provides the necessary definitions.

* Changelog

-}

module GammaRate where

import           Numeric.Integration.TanhSinh
import           Statistics.Distribution
import           Statistics.Distribution.Gamma

-- | For a given number of rate categories 'n' and a shape parameter 'alpha'
-- (the rate or scale is set such that the mean is 1.0), return a list of rates
-- that represent the respective categories. Use the mean rate for each
-- category.
getMeans :: Int -> Double -> [Double]
getMeans n alpha = means ++ lastMean
  where gamma = gammaDistr alpha (1.0/alpha)
        quantiles = [ quantile gamma (fromIntegral i / fromIntegral n) | i <- [0..n] ]
        -- Calculate the mean rate. Multiplication with the number of rate
        -- categories 'n' is necessary because in each n-quantile the
        -- probability mass is 1/n.
        meanFunc x = fromIntegral n * x * density gamma x
        -- Only calculate the first (n-1) categories with normal integration.
        means = [ intAToB meanFunc (quantiles !! i) (quantiles !! (i+1)) | i <- [0..n-2] ]
        -- The last category has to be calculated with an improper integration.
        lastMean = [intAToInf meanFunc (quantiles !! (n-1))]

-- | The error of integration.
eps :: Double
eps = 1e-6

-- | The integration method to use
method :: (Double -> Double ) -> Double -> Double -> [Result]
method = parSimpson

-- | Helper function for a normal integral from 'a' to 'b'.
intAToB :: (Double -> Double) -> Double -> Double -> Double
intAToB f a b = result . absolute eps $ method f a b

-- | Helper function for an improper integral from 'a' to infinity.
intAToInf :: (Double -> Double) -> Double -> Double
intAToInf f a = (result . absolute eps $ nonNegative method f) - intAToB f 0 a

max_stars_repo_name: dschrempf/bmm-simulate

max_issues_repo_name: pomo-dev/bmm-simulate

max_forks_repo_name: pomo-dev/bmm-simulate