text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MultiWayIf #-}
{-# LANGUAGE ParallelListComp #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}
module Language.Grappa.Inference.OldVariationalInference where

import Control.Monad ( forM, replicateM )
import Data.List ( transpose )

import Control.Lens ( (.~) )
import qualified Data.Colour as G
import qualified Data.Colour.Names as G
import qualified Data.Default.Class as G
import qualified Data.Matrix as Dm
import qualified Graphics.Rendering.Chart as G
import qualified Graphics.Rendering.Chart.Backend.Diagrams as G
import qualified Numeric.AD as Ad
import qualified Numeric.LinearAlgebra as La
import Numeric.LinearAlgebra ( Matrix, Vector, (<>), (<.>), (#>) )
import qualified Statistics.Distribution as S
import qualified Statistics.Distribution.StudentT as S
import qualified System.Random.MWC as MWC
import qualified Test.QuickCheck as Qc
import qualified Test.QuickCheck.Monadic as Qc
import qualified Test.QuickCheck.Property as Qc

import Language.Grappa.Distribution
import Language.Grappa.Inference.GradientDescent
import Language.Grappa.Rand.MWCRandM

import Debug.Trace

----------------------------------------------------------------
-- * BBVI algorithms
----------------------------------------------------------------
--
-- $bbviAlgorithms
-- Here is a summary of ways @conathan@ thought of to improve
-- BBVI. Only the first, simplest observation has been implemented, by
-- 'iteratedBbvi'.
--
-- When running BBVI to learn goal values @lambda@ that
-- are very far from the randomly chosen start values for
-- @lambda@, simply rerunning BBVI, starting from the previous
-- @lambda@, works well. For example, in
--
-- > runBbviTestUsingLazyRand 1
-- >   2      {- singleFactorBbvi2 -}
-- >   3      {- bbviTest3 -}
-- >   0.01   {- epsilon -}
-- >   100    {- num trials -}
-- >   [ 1000 {- mu -}
-- >   , 10 ] {- sigma -}
--
-- when lambda is initialized via
--
-- > replicateM 2 $ distSample (Normal 0 10)
--
-- the guesses are likely to be close to zero and far from 1000,
-- the goal value for @mu@.
--
-- The goal is to learn @lambda = [1000,10]@, and starting with
-- random lambda near zero, and then redefining @sampleLambda@ to
-- return the result, and iterating, we get the following sequence
-- of discovered lambda values:
--
-- @
-- sampleLambda = return $
--   -- [266.51484560478195,7.379162331856864]
--   -- [484.67445818993525,7.68279736309357]
--   -- [695.344938596016,8.239851183556492]
--   -- [874.483037696744,9.262297809983037]
--   -- [976.4159168688128,9.889695342218715]
--   -- [999.2564614475032,10.024071543947]
--   -- [999.9992140408281,9.999553691817553]
--   [999.99991729097,10.000438169812853]
-- @
--
-- As the result gets closer to the correct answer @[1000,10]@,
-- the incremental BBVI calls converge much more quickly. So, the
-- above corresponds to 8 runs of BBVI, but only take 3 or 4 times
-- as long as the first run.
--
-- With this improvement via iteration in mind, an obvious idea is
-- to simply implement a meta algorithm that does this iteration
-- for us: run BBVI in a loop, reusing the last @lambda@ to seed
-- the current run, until the @lambda@ computed by this meta
-- algorithm stabilizes (or some number of iterations are reached,
-- since for some functions, e.g. @sin@-like functions, we would
-- never converge).
--
-- And why does this meta algorithm work? Because when we restart
-- BBVI the learning rate / step size is large again, so we can
-- make lots of early progress before the decreasing learning rate
-- artificially slows us down.
--
-- This observation motivates a more general approach/goal: how
-- can we adjust the learning rate in a more intelligent way, so
-- that we don't slow down too early, but do slow down when we
-- reach a local optimum? ("Get while the gettin' is good" and
-- then stop.)
--
-- First, some ideas for detecting that we're not making progress
-- anymore:
--
-- 1. the difference between two consecutive values of lambda is
--    small. This is the approach the BBVI paper suggests. As
--    we've seen, it's misleading in conjunction with the
--    decreasing learning rates suggested by the paper.
--
-- 2. a whole group (tail of the sequence of lambdas) is
--    "clustered". The work here is defining "clustered". For
--    example, if the average distance (or max distance) between a
--    bunch of lambdas was close to the min distance (or average
--    distance), then the points are clustered.
--
-- 3.  the value of the objective function (ELBO) is not increasing
--     fast enough / in proportion to the effort. The work here
--     includes
--
--     - evaluating / estimating the objective function: turns out
--       this is quite similar to estimating the gradient of the
--       objective, and using a control variate should be similar
--       here. Indeed, we can define a control variate by adding a
--       multiple of the score function (as we do for estimating the
--       gradient itself.
--
--     - defining "fast enough / in proportion to the effort". For
--       example, we could define a "velocity of improvement" as the
--       ratio of the change in the objective to the step size.
--
-- Second, some ideas for ensuring we do keep making progress when
-- we can, by not decreasing the step size unnecessarily:
--
-- 1.  don't decrease the step size until the objective function stops
--     improving "fast enough", as determined by (3) above.
--
-- 2.  split the AdaGrad step size into two parts, a scalar learning
--     rate (@eta@ or @rho@ in the papers) and normalized AdaGrad
--     learning rate vector. That is, in place of
--
--     > diag(G^t)^{-1/2}
--
--     we use
--
--     > eta^t (d^t / || d^t ||_2)
--
--     for some @d^t@ based on @diag(G^t)^{-1/2}@.
--
--     While this approach punts a bit, by depending on a new
--     scalar learning rate @eta@, it simpler in that it decouples
--     the per component learning rate / emphasis part of AdaGrad,
--     a vector, from the scalar decrease in the step size. We can
--     use a notion of decreasing progress as above to wind down
--     the scalar learning rate.
--
--     For the @d^t@ we can simply use @diag(G^t)^{-1/2}@ itself,
--     or perhaps a forgetful version, e.g.
--
--     > diag(G^{t-k,t})^{-1/2}
--
--     for
--
--     > G^{t1,t2} := g^{t1} g^{t1}^T + ... + g^{t1} g^{t1}^T
--
--     i.e.
--
--     > G^t = G^{1,t}.
--
-- So, there are a lot of interesting things we could try ...

-- | The type of BBVI algorithms.
--
-- The @m@ is always monadic, but adding a @Monad m@ constraint on the
-- RHS here causes problems with "impredicative polymorphism" later,
-- when we try to make a list of @Bbvi m a@.
type Bbvi m a =
  R {-^ @epsilon@  -} ->
  Int {-^ @numSamples@ -} ->
  (Int -> [R] -> [R]) {-^ @projectLambda@ -} ->
  m [R] {-^ @sampleLambda@ -} ->
  ([R] -> m a) {-^ @sampleQ@  -} ->
  (a -> [R] -> m R) {-^ @evalLogQ@ -} ->
  (a -> [R] -> m [R]) {-^ @evalGradLoqQ@ -} ->
  (a -> m R) {-^ @evalLogP@ -} ->
  m [R]

-- | Black Box Variational Inference Algorithm 1.
--
-- From
-- http://www.cs.columbia.edu/~blei/papers/RanganathGerrishBlei2014.pdf.
--
-- We assume that @evalLogP@ is already (partially) applied to the
-- data @x@, unlike in the paper.
--
-- Here @epsilon@ is @0.01@ in the paper, and @numSamples@ is the
-- unspecified @S@ in the paper.
--
-- The @projectLambda@ function is used to project computed @lambda@
-- values back into the valid range of values for lambda.
--
-- === /Alternative approach/
--
-- As an alternative to defining a non-identity @projectLambda@, you
-- can "just" reparameterize to make @lambda@ range over all of
-- @R^n@. This is complicated and not suggested; see `bbviTest2` for
-- an example of reparameterizing.
--
-- You may be able to use functions from our 'Continuous' class to do
-- the transformations.
--
-- For example, normal distributions only support positive values for
-- sigma; we can work around this by using @fromRealLbInfty 0@ and
-- @toRealLbInfty 0@.
--
-- Alternatively, it might be possible to make the @q@-related
-- functions signal some kind of error when some part of lambda is out
-- of range. Not sure how to best recover here, but coming up with a
-- good way to recover could be much simpler than reparameterizing.
bbvi1 :: Monad m => Bbvi m a
bbvi1 epsilon numSamples projectLambda
  sampleLambda sampleQ evalLogQ evalGradLogQ evalLogP = do
  lambda <- projectLambda 1 <$> sampleLambda
  loop lambda 1
  where
    loop lambda t = do
      zs <- replicateM numSamples (sampleQ lambda)
      elboGradTerms <- forM zs $ \z -> do
        glq <- evalGradLogQ z lambda
        lq <- evalLogQ z lambda
        lp <- evalLogP z
        return $ (lp - lq) `scaleV` glq
      let elboGrad = (1 / fromIntegral numSamples) `scaleV`
                     sumV elboGradTerms
      let lambda' = projectLambda t $ lambda `addV` (rho t `scaleV` elboGrad)
      let delta = ell1 $ lambda `subV` lambda'
      if delta >= epsilon
      then loop lambda' (t + 1)
      else trace ("t = "++show t) $ return lambda'

    -- Paper says any sequence converging in @\ell^2@ and diverging in
    -- @\ell^1@ will do for @rho@ here.
    rho t = 1/fromIntegral t

-- | Single factor Black Box Variational Inference Algorithm 2,.
--
-- The n-factor algorithm can be implemented in terms of the single
-- factor algorithm, assuming the factors can be sampled
-- independently, which seems likely since their pdfs are independent
-- by assumption. The single factor Algorithm 2 is pretty similar to
-- Algorithm 1, except we use control variates to reduce the
-- variation.
--
-- See 'bbvi1' for more docs.
--
-- The use of @projectLambda@ in this implementation only corresponds
-- to the AdaGrad paper when the domain of lambda is a closed (hyper)
-- rectangle, but (partially) open rectangles should be close
-- enough. For non-rectangular domains, the projection function in the
-- AdaGrad paper takes @diagG@ into account via a "Mahalanobis norm"
-- (same as the @\ell^2@ norm for rectangular domains in this
-- context). I expect this doesn't matter much ...
--
-- TODO: change this into a stream algorithm, so that epsilon does not
-- have to be passed in, and instead the caller can just
--
-- > dropWhile (error > epsilon)
--
-- Chad suggested this approach earlier for coordinate descent.
singleFactorBbvi2 :: Monad m => Bbvi m a
singleFactorBbvi2 epsilon numSamples projectLambda
  sampleLambda sampleQ evalLogQ evalGradLogQ evalLogP = do
  lambda <- projectLambda 1 <$> sampleLambda
  let diagG = map (const 0) lambda
  loop diagG lambda 1
  where
    loop diagG lambda t = do
      zs <- replicateM numSamples (sampleQ lambda)
      hs <- forM zs $ \z -> do
        evalGradLogQ z lambda
      fs <- forM (zip zs hs) $ \(z, h) -> do
        lq <- evalLogQ z lambda
        lp <- evalLogP z
        return $ (lp - lq) `scaleV` h

      -- In the paper they say that they used 1000 for @numSamples@
      -- (called @S@ in the paper) and 100 estimates (i.e. 1/10th of
      -- the 1000) for the approximate Cov and Var in the computation
      -- of @\hat{a}^*_d@. So, we're generalizing this to using 1/10th
      -- of the samples.
      let numSamples' = if numSamples <= 10
                        then numSamples
                        else numSamples `div` 10
      let hs' = take numSamples' hs
      let fs' = take numSamples' fs
      let hs't = transpose hs'
      let fs't = transpose fs'
      let as = [ cov fCol hCol / var hCol | fCol <- fs't | hCol <- hs't ]

      let elboGradTerms = [ f `subV` aH
                          | f <- fs
                          | aH <- map (multV as) hs ]
      let elboGrad = (1 / fromIntegral numSamples) `scaleV` sumV elboGradTerms
      let diagG' = diagG `addV` (elboGrad `multV` elboGrad)
      -- Square root of pseudo inverse of @diagG'@.
      let learningRate = eta `scaleV`
                         [ if x == 0 then 0 else x**(-0.5) | x <- diagG' ]
      let dLambda = learningRate `multV` elboGrad
      let lambda' = projectLambda t $ lambda `addV` dLambda

      let delta = ell1 $ lambda `subV` lambda'
      if trace ("delta = "++show delta++"\ndiagG = "++show diagG'++"\ndLambda = "++show dLambda++"\nlambda = "++show lambda') $
         delta >= epsilon
      then loop diagG' lambda' (t + 1)
      else trace ("t = "++show t) $ return lambda'

    -- I think the AdaGrad paper said this parameter doesn't matter
    -- much (e.g. the theorem about the regret going to zero doesn't
    -- depend on an optimal choice for eta). I think this is the value
    -- they used though.
    eta = 2**0.5

    -- Note that @cov@ and @var@ are unnormalized -- i.e. not divided
    -- by the length of the input vectors -- since the normalization
    -- factors would just cancel in the calculation of @as@ as a ratio
    -- of @cov@ to @var@.
    var xs = cov xs xs
    cov xs ys = sum [ (x - meanXs) * (y - meanYs) | x <- xs | y <- ys ]
      where
        meanXs = mean xs
        meanYs = mean ys
    mean xs = sum xs / fromIntegral (length xs)

-- | Iterate BBVI until convergence or round limit is reached.
--
-- The @bbvi@ argument is expected to be a partially applied to all
-- arguments except for @sampleLambda@.
iterateBbvi ::
  Monad m =>
  Int ->
  R ->
  m [R] ->
  (m [R] -> m [R]) ->
  m [R]
iterateBbvi maxRounds epsilon sampleLambda bbvi = do
  result <- bbvi sampleLambda
  go 2 result
  where
    go round result
      | round >= maxRounds = return result
      | otherwise = do
          result' <- bbvi (return result)
          if ell1 (result `subV` result') <= epsilon
            then return result'
            else trace (show result') $ go (round + 1) result'

-- | Create an iterated version of an existing BBVI algorithm.
--
-- A helper to call @iterateBbvi@ for you, producing something of type
-- @Bbvi@' (useful for tests where I want to put BBVI algs in a list).
iteratedBbvi :: Monad m => Int -> Bbvi m a -> Bbvi m a
iteratedBbvi maxRounds bbvi epsilon numSamples projectLambda
  sampleLambda sampleQ evalLogQ evalGradLogQ evalLogP =
  iterateBbvi maxRounds epsilon sampleLambda
  (\sample -> bbvi epsilon numSamples projectLambda
              sample sampleQ evalLogQ evalGradLogQ evalLogP)

-- | A variety of BBVI algs for use in testing.
bbvis :: Monad m => [Bbvi m a]
bbvis =  [ bbvi1
         , singleFactorBbvi2
         , iteratedBbvi 100 bbvi1
         , iteratedBbvi 100 singleFactorBbvi2 ]

----------------------------------------------------------------
-- * Variational families
----------------------------------------------------------------
--
-- $variationalFamilies
-- The code needed to use a specific variational family (the @q \in
-- Q@) with a BBVI algorithm is independent of the target distribution
-- (the @p@) being approximated.

----------------------------------------------------------------
-- ** Normal

sampleNormal :: SampleableIn m Normal => [R] -> m R
sampleNormal [mu, sigma] = distSample (Normal mu sigma)

evalLogNormalDensity :: Monad m => R -> [R] -> m R
evalLogNormalDensity z [mu, sigma] = return $
  probToLogR $ distDensity (Normal mu sigma) z

evalGradLogNormalDensity :: Monad m => R -> [R] -> m [R]
evalGradLogNormalDensity z [mu, sigma] = return
  -- Calculated by hand from @normalDensity@ ...
  [ (z - mu) / sigma
  , (z - mu)**2 / sigma**3 - 1 / sigma ]

projectLambdaNormal :: Int -> [R] -> [R]
projectLambdaNormal = projectRectangle
  [ (OpenBound (-1/0), OpenBound (1/0)),
    (OpenBound 0, OpenBound (1/0)) ]

----------------------------------------------------------------
-- ** Multivariate normal
--
-- $mvNormal
-- Everywhere below @k@ is the dimension of the mv normal.

-- | Convert the flat list of params used by the bbvi algorithm into
-- the vector and matrix params used by the mv normal distribution.
paramsToDistParamsMVNormal :: Int -> [Double] -> (Vector Double, Matrix Double)
paramsToDistParamsMVNormal k muAndC =
  if length muAndC /= expectedLength
  then error $
       "paramsToDistParamsMVNormal: "++
       "the parameter list is the wrong length! Expected: "++
       show expectedLength++", actual: "++show (length muAndC)
  else (mu, c)
  where
    -- An upper triangular @k x k@ matrix has @k + 1@ choose 2
    -- potentially non-zero entries.
    expectedLength = k + ((k * (k+1)) `div` 2)

    (muList, cList) = splitAt k muAndC
    mu = La.fromList muList
    c = (k La.>< k) $ upperTriangular k cList

-- | Build full list of entries in upper triangular matrix from list
-- of its potentially non-zero entries.
--
-- Builds up the full list for the matrix by filling in zeros for
-- undefined entries. The undefined entries are the lower tringular,
-- i.e. where the row index is larger than the column index.
upperTriangular :: (Num a) => Int -> [a] -> [a]
upperTriangular k xs0 = zeroPad rcs xs0
  where
    rcs = [ (r, c) | r <- [0..(k-1)], c <- [0..(k-1)] ]
    zeroPad []           []         = []
    zeroPad ((r,c):rcs') xs@(x:xs') =
      if r > c
      then 0 : zeroPad rcs' xs
      else x : zeroPad rcs' xs'

sampleMVNormal :: SampleableIn m MVNormal => Int -> [R] -> m (Vector R)
sampleMVNormal k muAndC = distSample (MVNormal mu c)
  where
    (mu, c) = paramsToDistParamsMVNormal k muAndC

-- | The @k@ is the dimension of the mv normal.
evalLogMVNormalDensity ::
  Monad m => Int -> Vector R -> [R] -> m R
evalLogMVNormalDensity k z muAndC = return $
  probToLogR $ distDensity (MVNormal mu c) z
  where
    (mu, c) = paramsToDistParamsMVNormal k muAndC

-- | The log pdf is
--
-- > lpdf = -1/2 ( ln |C^T C| + (x-mu)^T(C^T C)^{-1}(x-mu) + k ln (2 pi) ) ,
--
-- where @|M|@ is the absolute value of the determinant of @M@.
--
-- The partials with respect to @mu@'s components are easy. The
-- product rule gives
--
-- > d_{mu_i} lpdf = -1/2 ( m^T(C^T C)^{-1}(x-mu) + (x-mu)^T(C^T C)^{-1}m )
--
-- for @m := - d_{mu_i} mu@, i.e. a vector with @m_j = - (Kronecker) delta i j@.
--
-- The partials with respect to @C@'s components are trickier, but
-- simplified by the fact that the determinant of @C@ is simple.
--
-- Let @d@ be some partial @d/dC_{ij}@. It comes down to calculating
--
-- > d ln |C^T C|
--
-- and
--
-- > d (x-mu)^T(C^T C)^{-1}(x-mu) .
--
--
-- For the first term we have @|C^T C|@ = @|C|^2@ and so
--
-- > ln |C^T C| = 2 ln |C|
--
-- and so
--
-- > d ln |C^T C| = (2/|C|) d |C| .
--
-- Since @C@ is upper triangular with positive diagonal, we have
--
-- > |C| = C_{11} .. C_{kk}
--
-- and so for @d = d/dC_{ij}@ we have
--
-- > d |C| = if i == j then |C|/C_{ij} else 0 .
--
-- and so
--
-- > d ln |C^T C| = if i == j then 2/C_{ij} else 0 .
--
--
-- For the second term we have
--
-- > (C^T C)^{-1} = (C^{-1}) (C^{-1})^T
--
-- and so
--
-- > d (C^T C)^{-1} = (d C^{-1}) (C^{-1})^T + C^{-1} (d C^{-1})^T ,
--
-- and so we just need the derivative of the inverse:
--
-- > d C^{-1} = - C^{-1} (d C) C^{-1}
--
-- see http://planetmath.org/derivativeofinversematrix.
--
-- Putting that all back together, with @B := C^{-1}@, we have
--
-- > d (C^T C)^{-1} = - (B dC B B^T + B B^T dC^T B^T)
-- >                = - B (dC B + (dC B)^T) B^T
--
-- and for @d = d/d(C_{ij})@ we have
--
-- > (dC)_{st} = delta i s * delta j t .
evalGradLogMVNormalDensity ::
  Monad m => Int -> Vector Double -> [R] -> m [R]
evalGradLogMVNormalDensity k x muAndC =
  return $ map (* (-1/2)) $ muGrads ++ cGrads
  where
    (mu, c) = paramsToDistParamsMVNormal k muAndC
    delta i j = if i == j then 1 else 0

    muGrads = [ d <.> cInvFun (x - mu) +
                (x - mu) <.> cInvFun d
              | i <- [ 0 .. k-1 ]
              , d <- [ La.fromList $ map (((-1) *) . delta i) [ 0 .. k-1 ] ] ]
    -- Efficient left multiplication of a vector by @(C^T C)^{-1}@.
    cInvFun = La.flatten . La.cholSolve c . La.asColumn

    cGrads = [ lnGrad i j + prodGrad i j
             | i <- [ 0 .. k-1 ]
             , j <- [ i .. k-1 ] ]
    lnGrad i j = if i == j then 2 / (c `La.atIndex` (i, j)) else 0
    prodGrad i j = (x - mu) <.> middle #> (x - mu)
      where
        middle = (-1) * b <> (dC i j <> b + La.tr (dC i j <> b)) <> La.tr b
    -- The inverse of @C@ as a matrix.
    (b, _) = La.invlndet c
    -- The derivative of @C@ w.r.t. @C_{ij}@.
    dC i j = La.assoc (k, k) 0 [ ((i, j), 1) ]

-- | The diagonal entries of the Cholesky matrix must be positive. The
-- other values are unconstrained.
--
-- The @k@ is the dimension of the mv normal.
projectLambdaMVNormal :: Int -> Int -> [R] -> [R]
projectLambdaMVNormal k = projectRectangle $ muBounds ++ cBounds
  where
    reals = (OpenBound (-1/0), OpenBound (1/0))
    positiveReals = (OpenBound 0, OpenBound (1/0))

    muBounds = replicate k reals
    cBounds = flip concatMap [ k, k - 1 .. 1 ] $ \rowLength ->
      take rowLength $ [ positiveReals ] ++ repeat reals

----------------------------------------------------------------
-- *** Test mv normal implementation.
--
-- $testMvNormal
--
-- The gradient we calculated by hand was quite complicated, so
-- compare it with the gradient we get from automatic differentiation
-- (would have been much faster in human time to just use automatic
-- differentiation in the first place, but now we have very high
-- confidence the gradient is correct).

-- | An overloaded (i.e. @Floating@) implementation of the log pdf.
--
-- This allows us to use automatic differentiation and compare the
-- result with the manual, by hand gradient
-- 'evalGradLogMVNormalDensity' above.
logMVNormalDensityOverloaded ::
  (Floating a, Eq a) => Int -> [a] -> [a] -> a
logMVNormalDensityOverloaded k xList muAndC =
  ((- 0.5) *) $
  2 * log detC +
  fromIntegral k * log (2 * pi) +
  -- The @Dm@ library uses 1-based indexing.
  prod Dm.! (1,1)
  where
    (muList, cList) = splitAt k muAndC
    x = Dm.fromList k 1 xList
    mu = Dm.fromList k 1 muList
    c = Dm.fromList k k $ upperTriangular k cList
    detC = foldr (*) 1 $ Dm.getDiag c
    sigma = Dm.transpose c * c
    Right sigmaInv = Dm.inverse sigma
    prod = (Dm.transpose $ x - mu) * sigmaInv * (x - mu)

-- | Use automatic differentiation to check the manually computed mv
-- normal gradient.
compareGrads :: Monad m => Int -> [R] -> [R] -> m [R]
compareGrads k x muAndC = do
  manual <- evalGradLogMVNormalDensity k (La.fromList x) muAndC
  -- AD needs the function to differentiate to overloaded / take
  -- @Floating@ arguments.
  let x' :: Floating a => [a]
      x' = map (fromRational . toRational) x
  let auto = Ad.grad (logMVNormalDensityOverloaded k x') muAndC
  -- The rescaling of errors is essential here: without it, for very
  -- large gradient entries (e.g. size 10^10), there can be large
  -- absolute differences between gradient entries. We care about
  -- large relative differences, which can't be explained away by
  -- rounding error. Dividing by 1 here avoids division by zero, and
  -- means we avoid spurious differences detected when we divide by
  -- very small numbers.
  return $ [ abs (m - a) / (abs m + abs a + 1) | m <- manual | a <- auto ]

-- | QuickCheck property asserting 'compareGrads' for given mv normal
-- dimension @k@.
prop_evalGradLogMVNormalDensity_correct :: Int -> Qc.Property
prop_evalGradLogMVNormalDensity_correct k =
  -- forAll' "k" Qc.arbitrary $ \(Qc.Positive (k :: Int)) ->
  forAll' "xList" (Qc.vector k) $ \xList ->
  forAll' "muList" (Qc.vector k) $ \muList ->
  forAll' "cList" (genCList k) $ \cList ->
  Qc.monadicIO $ do
    deltas <- Qc.run $ compareGrads k (i xList) (i $ muList ++ cList)
    Qc.monitor (Qc.counterexample $ "(Index, Delta): " ++ show (argMax deltas))
    Qc.assert $ all isSmall deltas
  where
    -- Return index of max entry, and max entry.
    argMax :: Ord a => [a] -> (Int, a)
    argMax (x:xs) = go 1 (0,x) xs
      where
        go _ p [] = p
        go j (i,x) (x':xs) = go (j+1) (if x < x' then (j, x') else (i, x)) xs

    i :: [Integer] -> [R]
    i = map fromInteger
    -- Generate a list representation of a Cholesky factor. The
    -- strategy is to recursively generate row tails from the diagonal
    -- entry onwards, with non-zero diagonal entries and arbitrary
    -- other entries.
    genCList :: Int -> Qc.Gen [Integer]
    genCList k = do
      rows <- forM [k,k-1..1] $ \len -> do
        Qc.Positive diagEntry <- Qc.arbitrary
        rest <- Qc.vectorOf (len-1) Qc.arbitrary
        return $ diagEntry : rest
      return $ concat rows

    isSmall delta = abs delta <= 0.001

    -- A version of @Qc.forAll@ with named quantifiers. This makes the
    -- counterexamples much easier to read.
    forAll' name gen pf =
      -- The library provided @Qc.forAll@ that I adapted to create
      -- this @forAll'@ has a @Qc.again@ here. But I don't understand
      -- why I need it, and more importantly, I get an import error,
      -- so leaving it out ...
      --
      -- Qc.again $
      Qc.MkProperty $
      gen >>= \x ->
        Qc.unProperty $
        Qc.counterexample (name++": "++show x) (pf x)

-- | Check the QuickCheck props @numTrials@ many times.
--
-- For larger @k@ this fails, but I think it's due to numerical
-- instability. For small @k@, e.g. less than or equal to 5, we can
-- pass hundreds of thousands of tests. For large @k@, e.g. @k = 30@,
-- we fail pretty fast.
checkProps :: Int -> Int -> IO ()
checkProps k numTrials = do
  qc $ prop_evalGradLogMVNormalDensity_correct k
  where
    qc = Qc.quickCheckWith (Qc.stdArgs { Qc.maxSuccess = numTrials })

----------------------------------------------------------------
-- * Tests
----------------------------------------------------------------

-- | Even simpler version of the 'bbv1Test2' below.
--
-- Here we approximate a normal by a normal. The standard deviation is
-- fixed, and we fit the mean. The input @mu@ and @sigma@ are the
-- parameters of the normal that we are trying to approximate; our
-- goal is to discover @mu@.
--
-- We make @mu@ the parameter to discover, because it ranges over all
-- reals and so no transformations are needed.
--
-- Returns the error in the discovered @mu@, and the discovered @mu@.
bbviTest1 ::
  (SampleableIn m Normal) =>
  Int -> R -> Int -> [R] -> m (R, [R])
bbviTest1 algNum epsilon numSamples [mu, sigma] = do
  [mu'] <- (bbvis !! (algNum - 1)) epsilon numSamples projectLambda
           sampleLambda sampleQ evalLogQ evalGradLogQ evalLogP
  let error' = ell1 $ [mu'] `subV` [mu]
  return (error', [mu'])
  where
    -- We can see that bbvi1 at least kind of works by seeding @lambda@
    -- with a value very close to the right answer.
    --
    -- sampleLambda = replicateM 1 $ distSample (Normal mu 0.1)
    sampleLambda = defaultSampleLambda 1
    sampleQ [mu'] = distSample (Normal mu' sigma)
    evalLogQ z [mu'] = return $
      probToLogR $ distDensity (Normal mu' sigma) z
    evalGradLogQ z [mu'] = return [ (z - mu') / sigma ]
    evalLogP z = return $
      probToLogR $ distDensity (Normal mu sigma) z
    projectLambda _ lambda = lambda

-- | Test 'bbv1' by approximating a linear combination of normals by a
-- normal ... except it turns out that a linear combination of normals
-- *is* a normal, and the only way to analytically weight its samples
-- is to know this, so I'm really just going to approximate a normal
-- by a normal. At least we can expect this to succeed :D
--
-- Formulas for linear combination of normals are here:
-- https://www.statlect.com/probability-distributions/normal-distribution-linear-combinations.
--
-- Note that we have to transform the standard deviation param sigma,
-- per the comments above in 'bbvi1'.
bbviTest2 ::
  (SampleableIn m Normal) =>
  Int -> R -> Int -> [R] -> m (R, [R])
bbviTest2 algNum epsilon numSamples [mu, sigma] = do
  [mu', sigma'] <- fRLambda <$>
    (bbvis !! (algNum - 1)) epsilon numSamples projectLambda
    sampleLambda sampleQ evalLogQ evalGradLogQ evalLogP
  -- In general, it might only make sense to compute the error in the
  -- real domain.
  let error' = ell1 $ [mu', sigma'] `subV` [mu, sigma]
  return (error', [mu', sigma'])
  where
    projectLambda _ lambda = lambda
    -- Since lambda is always a vector of reals, it's not clear why it
    -- should be a param to 'bbvi1'.
    --
    -- !!!: THE @fromRealLbInfty 0@ UNDERFLOWS TO ZERO FOR LARGE VALUES.
    --
    -- ???: does the delta in the output change when we increase
    -- numSampls???
    sampleLambda = defaultSampleLambda 2

    sampleQOrig = sampleNormal
    sampleQ = sampleQOrig . fRLambda

    evalLogQOrig = evalLogNormalDensity
    evalLogQ z = evalLogQOrig z . fRLambda

    -- Calculated by hand from @normalDensity@ ...
    evalGradLogQOrig = evalGradLogNormalDensity
    evalGradLogQ z lambda =
      -- Chain rule specialized to vars transformed independently
      -- (i.e. the derivative matrix of the transformation @fRLambda@
      -- is diagonal).
      zipWith (*) (fRLambda' lambda) <$>
      evalGradLogQOrig z (fRLambda lambda)

    evalLogP z = evalLogNormalDensity z [mu, sigma]

    -- Transformations from and to reals for sigma param of
    -- normal. Note: the 'Continuous' class transforms the support /
    -- output, not the params / input of a distribution.
    fRSigma = fromRealLbInfty 0
    fRSigma' = probToR . fromRealLbInfty' 0
    -- tRSigma = toRealLbInfty 0

    fRLambda [mu', sigma'R] = [mu', fRSigma sigma'R]
    -- Here @1@ is the derivative of the identity.
    fRLambda' [_mu', sigma'R] = [1, fRSigma' sigma'R]
    -- tRLambda [mu', sigma'] = [mu', tRSigma sigma']

-- | Like `bbviTest2`, but use a non-trivial `projectLambda`, instead
-- of a reparameterization.
bbviTest3 ::
  (SampleableIn m Normal) =>
  Int -> R -> Int -> [R] -> m (R, [R])
bbviTest3 algNum epsilon numSamples [mu, sigma] = do
  [mu', sigma'] <-
    (bbvis !! (algNum - 1)) epsilon numSamples projectLambda
    sampleLambda sampleQOrig evalLogQOrig evalGradLogQOrig evalLogP
  -- In general, it might only make sense to compute the error in the
  -- real domain.
  let error' = ell1 $ [mu', sigma'] `subV` [mu, sigma]
  return (error', [mu', sigma'])
  where
    -- Since lambda is always a vector of reals, it's not clear why it
    -- should be a param to 'bbvi1'.
    --
    -- !!!: THE @fromRealLbInfty 0@ UNDERFLOWS TO ZERO FOR LARGE VALUES.
    sampleLambda = defaultSampleLambda 2
    sampleQOrig = sampleNormal
    evalLogQOrig = evalLogNormalDensity
    evalGradLogQOrig = evalGradLogNormalDensity
    projectLambda = projectLambdaNormal

    evalLogP z = evalLogNormalDensity z [mu, sigma]

-- | Learn a normal that approximates a Student's t-distribution.
--
-- The parameters are @mu@ real and @nu@ positive real.
--
-- We want the mean of the learned normal to be the mean of the
-- t-distribution. I don't know what to expect for the std deviation
-- of the learned normal tho. The variance of the t-distribution is
-- @nu / (nu - 2)@ for @nu > 2@ and infinite or undefined otherwise.
--
-- See https://en.wikipedia.org/wiki/Student's_t-distribution.
bbviTest4 ::
  SampleableIn m Normal =>
  Int -> R -> Int -> [R] -> m (R, [R])
bbviTest4 algNum epsilon numSamples [mu, nu] = do
  [mu', sigma'] <-
    (bbvis !! (algNum - 1)) epsilon numSamples projectLambda
    sampleLambda sampleQ evalLogQ evalGradLogQ evalLogP
  -- In general, it might only make sense to compute the error in the
  -- real domain.
  let error' = ell1 $ [mu'] `subV` [mu]
  return (error', [mu', sigma'])
  where
    sampleQ = sampleNormal
    evalLogQ = evalLogNormalDensity
    evalGradLogQ = evalGradLogNormalDensity
    sampleLambda = defaultSampleLambda 2
    projectLambda = projectLambdaNormal

    evalLogP z = return $
      S.logDensity (S.studentT nu) (z - mu)

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

-- | Learn a multivariate normal that approximates a linear
-- regression.
--
-- The linear regression has parameters slope @m@, @y@-intercept @b@,
-- and error standard deviation @sigma@. The generative model
-- comprises a arbitrary distribution of the parameters @m@, @b@, and
-- @sigma@, say
--
-- > m ~ Uniform 0 10
-- > b ~ Uniform 0 10
-- > sigma ~ Uniform 0 10 ,
--
-- and a noisy linear distribution on the dependent variable @y@ given
-- the independent variable @x@, i.e.
--
-- > e ~ Normal 0 sigma
-- > y := m*x + b + e ,
--
-- or equivalently
--
-- > y ~ Normal (m*x + b) e .
--
-- So, the joint distribution given data
--
-- > xys = [(x1,y1),...,(xn,yn)]
--
-- is
--
-- > P[xys,m,b,sigma] = P[xys|m,b,sigma] * P[m,b,sigma] ,
--
-- where
--
-- > P[m,b,sigma] = 1
--
-- and
--
-- > P[xys|m,b,sigma] = p(x1,y1) * ... * p(xn,yn) ,
--
-- for
--
-- > p(x,y) =
-- > Pr[y = m*x + b + e] =
-- > Pr[e = y - m*x - b] =
-- > distDensity (Normal (m*x + b) sigma) y .
--
-- The @xysFlat@ argument is
--
-- > [x1,y1,x2,y2,...,xn,yn] .
bbviTest5 ::
  (SampleableIn m Normal, SampleableIn m MVNormal) =>
  Int -> R -> Int -> [R] -> m (R, [R])
bbviTest5 algNum epsilon numSamples (m:b:sigma:xysFlat) = do
  muAndC'@(m':b':sigma':_cList') <-
    (bbvis !! (algNum - 1)) epsilon numSamples projectLambda
    sampleLambda sampleQ evalLogQ evalGradLogQ evalLogP
  let error' = ell1 $ [m',b',sigma'] `subV` [m,b,sigma]
  return (error', muAndC')
  where
    -- Three dimensions, for @m@, @b@, and @sigma@.
    k = 3
    muSize = k
    cSize = (k * (k+1)) `div` 2

    sampleQ = sampleMVNormal k
    evalLogQ = evalLogMVNormalDensity k
    evalGradLogQ = evalGradLogMVNormalDensity k
    sampleLambda = defaultSampleLambda (muSize + cSize)
    projectLambda = projectLambdaMVNormal k

    xys = go xysFlat
      where
        go [] = []
        go (x:y:xysFlat') = (x,y) : go xysFlat'

    evalLogP mBSigma = return $ linRegLogDensity xys m' b' sigma'
      where
        [m', b', sigma'] = La.toList mBSigma

-- | A version of 'bbviTest5' that fixes the noise deviation sigma, so
-- that the mv normal only needs two params. The sigma passed in will
-- be used, and should correspond to the supplied data in @xysFlat@.
bbviTest7 ::
  (SampleableIn m Normal, SampleableIn m MVNormal) =>
  Int -> R -> Int -> [R] -> m (R, [R])
bbviTest7 algNum epsilon numSamples (m:b:sigma:xysFlat) = do
  muAndC'@(m':b':_cList') <-
    (bbvis !! (algNum - 1)) epsilon numSamples projectLambda
    sampleLambda sampleQ evalLogQ evalGradLogQ evalLogP
  let error' = ell1 $ [m',b'] `subV` [m,b]
  return (error', muAndC')
  where
    -- Two dimensions, for @m@, @b@.
    k = 2
    muSize = k
    cSize = (k * (k+1)) `div` 2

    sampleQ = sampleMVNormal k
    evalLogQ = evalLogMVNormalDensity k
    evalGradLogQ = evalGradLogMVNormalDensity k
    sampleLambda = defaultSampleLambda (muSize + cSize)
    projectLambda = projectLambdaMVNormal k

    xys = go xysFlat
      where
        go [] = []
        go (x:y:xysFlat') = (x,y) : go xysFlat'

    evalLogP mB = return $ linRegLogDensity xys m' b' sigma
      where
        [m', b'] = La.toList mB

linRegLogDensity ::
  [(R, R)] -> R -> R -> R -> R
linRegLogDensity xys m b sigma =
  -- An alternative approach:
  --
  -- The @abs sigma'@ is here to avoid using negative sigma
  -- values. Not sure what the right approach is here: we could
  -- use a different @projectLambda@ implementation that restricts
  -- to positive third component (i.e. the mean of the sigma
  -- values), but that won't stop them from being
  -- negative.
  {-
  sum [ logFromLogFloat (distDensity (Normal (m'*x + b') (abs sigma')) y)
      | (x, y) <- xys ]
  -}
  -- We penalize negative sigmas, by returning a very small value
  -- here.
  if sigma > 0
  then
    (if logProb < logVerySmallValue
     then trace ("XXX: logProb is very small: "++show logProb)
     else id)
    logProb
  else
    -- Is this value is really small enough? Could try @-1/0@.
    --
    -- Update: using negative infinity, or more generally very large
    -- negative values here, causes various crashes (for finite values
    -- an underflow in @LogFloat@, and for negative infinity a
    -- malformed Cholesky factor).
    logVerySmallValue
    where
      logVerySmallValue = -10e10
      logProb = sum [ probToLogR $ distDensity (Normal (m*x + b) sigma) y
                    | (x, y) <- xys ]

simpleLinRegData :: Int -> R -> R -> [(R, R)]
simpleLinRegData numPoints m b =
  [ (x, m*x+b) | x <- [ 1 .. fromIntegral numPoints ] ]

genLinRegDataIO :: Int -> R -> R -> R -> IO [(R, R)]
genLinRegDataIO numSamples m b sigma = do
  g <- MWC.create
  runRand g $ genLinRegData numSamples m b sigma

genLinRegData ::
  ( SampleableIn m Normal, SampleableIn m Uniform ) =>
  Int -> R -> R -> R -> m [(R, R)]
genLinRegData numSamples m b sigma = do
  replicateM numSamples $ do
    -- Use integral @x@s to make the output a little easier to inspect
    -- visually.
    x <- fromInteger . round <$> distSample (Uniform (-1000) 1000)
    y <- distSample (Normal (m*x + b) sigma)
    return (x,y)

-- scatterPlot :: [(R, R)] -> ???
-- | Based on http://indiana.edu/~ppaml/HakaruTutorial.html and
-- https://github.com/timbod7/haskell-chart/wiki/example%202.
--
-- Writes output to an SVG file, regardless of supplied file extension
-- :P The examples I referred to use the Cairo backend, but I'm using
-- the Diagrams backend, so that might be the issue. So, I changed
-- this function to supply the extension, to avoid confusion.
scatterPlot :: String -> String -> [(R, R)] -> IO (G.PickFn ())
scatterPlot title file xys =
  G.renderableToFile G.def (file++".svg") $ G.toRenderable layout
  where
    plot = G.plot_points_style .~ G.filledCircles 2 (G.opaque G.purple)
           $ G.plot_points_values .~ xys
           $ G.def

    layout = G.layout_title .~ title
             $ G.layout_plots .~ [G.toPlot plot]
             $ G.layout_x_axis . G.laxis_generate .~ G.scaledAxis G.def (-4, 4)
             $ G.layout_y_axis . G.laxis_generate .~ G.scaledAxis G.def (-4, 4)
             $ G.def

-- | https://en.wikipedia.org/wiki/File:MultivariateNormal.png
--
-- I can't figure out how to set the aspect ratio :P
wikipediaExample :: IO ()
wikipediaExample = do
  g <- MWC.create
  xyVectors <- runRand g $ replicateM numSamples (distSample (MVNormal mu c))
  let xys = [ (x, y) | v <- xyVectors, [x,y] <- [ La.toList v ] ]
  _ <- scatterPlot ("Wikipedia Test ("++show numSamples++" samples)") "/tmp/out" xys
  return ()
  where
    mu = La.fromList [0, 0]
    c = La.chol (La.trustSym $ (2 La.>< 2) [1, 0.6, 0.6, 2])
    numSamples = 10000

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

-- | Sanity test: match an mv normal to an mv normal.
--
-- This works pretty well with iterated bbvi 2, e.g. with
--
-- > runBbviTestUsingLazyRand 1 4 6 1 1000 [1,2,3,1,0,0,2,0,3]
--
-- Of course, this doesn't prove much.
bbviTest6 ::
  (SampleableIn m Normal, SampleableIn m MVNormal) =>
  Int -> R -> Int -> [R] -> m (R, [R])
bbviTest6 algNum epsilon numSamples muAndC = do
  muAndC' <-
    (bbvis !! (algNum - 1)) epsilon numSamples projectLambda
    sampleLambda sampleQ evalLogQ evalGradLogQ evalLogP
  let error' = ell1 $ muAndC' `subV` muAndC
  return (error', muAndC')
  where
    k = 3
    muSize = k
    cSize = (k * (k+1)) `div` 2

    sampleQ = sampleMVNormal k
    evalLogQ = evalLogMVNormalDensity k
    evalGradLogQ = evalGradLogMVNormalDensity k
    sampleLambda = defaultSampleLambda (muSize + cSize)
    projectLambda = projectLambdaMVNormal k

    evalLogP z = evalLogMVNormalDensity k z muAndC

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

runBbviTestUsingLazyRand ::
  Int -> Int -> Int -> R -> Int -> [R] -> IO [(R, [R])]
runBbviTestUsingLazyRand numTrials algNum testNum epsilon numSamples testParams = do
  g <- MWC.create
  runRand g $ replicateM numTrials $ test algNum epsilon numSamples testParams
  where
    test = [ bbviTest1
           , bbviTest2
           , bbviTest3
           , bbviTest4
           , bbviTest5 -- Mv normal for linear regression: learn m, b, sigma.
           , bbviTest6 -- Mv normal for mv normal.
           , bbviTest7 -- Mv normal for linear regression: learn m, b.
           ] !! (testNum - 1)

test5 :: IO [(R, [R])]
test5 = runBbviTestUsingLazyRand 1 3 5 0.1 1000
  ([m,b,sigma] ++ concat [ [x, y] | (x, y) <- simpleLinRegData 1000 m b ])
  where
    m = 1
    b = 2
    sigma = 0.001


text:
module Main where
  import Data.Picture
  import System.Environment
  import Data.Either
  import Numeric.LinearAlgebra (rows, cols)

  data Options = Options { file :: FilePath
                         , output :: FilePath
                         , argInvert :: Bool
                         , argGrayscale :: Bool
                         , argRotate :: Double
                         , argFade :: Double
                         , argContrast :: Double
                         , argGamma :: Int
                         , argBrightness :: Double
                         , argCompress :: Int
                         , argScale :: Double
                         }

  opts = Options { file = ""
                 , output = "output.png"
                 , argInvert = False
                 , argGrayscale = False
                 , argRotate = 0
                 , argFade = 100
                 , argContrast = 0
                 , argGamma = 1
                 , argBrightness = 0
                 , argCompress = 0
                 , argScale = 1
                 }

  main :: IO ()
  main = do
    args <- getArgs

    if null args
      then do
        putStrLn "Usage: picedit <input> [OPTIONS]"
        putStrLn "Options:"
        putStrLn "  --contrast <n> - a number between -255 and 255"
        putStrLn "  --brightness <n> - a number between -255 and 255"
        putStrLn "  --gamma <n>"
        putStrLn "  --fade <n> - a number between 0 and 100"
        putStrLn "  --rotate <n> - rotate image by n degrees"
        putStrLn "  --grayscale - turn the image grayscale"
        putStrLn "  --invert - invert (negative) the image"
        putStrLn "  --scale <n> - scale the image using nearest-neighbor interpolation"
        putStrLn "  --compress <n> - approximate the (width - n)-th rank of image using SVD, note: this is not size compression, a number between 0 (no compression) and image width (full compression)"
        putStrLn "  --output <filename> - output name, defaults to output.png"
      else do
        let options = parseArgs args opts

        pic <- readPicture (file options)
        Right other <- readPicture ("output.png")

        case pic of
          Left err -> print err
          Right p -> do
            let edited = rotate (argRotate options) Nothing
                       . fade (argFade options / 100)
                       . contrast (argContrast options)
                       . gamma (argGamma options)
                       . brightness (argBrightness options)
                       . conditionalFn grayscale (argGrayscale options)
                       . conditionalFn invert (argInvert options)
                       . scale (argScale options)
                       . compress (argCompress options) $ p
            writePicturePng (output options) edited

        return ()

    where
      conditionalFn f True = f
      conditionalFn f False = id

  parseArgs :: [String] -> Options -> Options
  parseArgs [] opts = opts
  parseArgs ("--invert":rest) opts = parseArgs rest (opts { argInvert = True })
  parseArgs ("--grayscale":rest) opts = parseArgs rest (opts { argGrayscale = True })
  parseArgs ("--rotate":n:rest) opts = parseArgs rest (opts { argRotate = read n })
  parseArgs ("--fade":n:rest) opts = parseArgs rest (opts { argFade = read n })
  parseArgs ("--contrast":n:rest) opts = parseArgs rest (opts { argContrast = read n })
  parseArgs ("--brightness":n:rest) opts = parseArgs rest (opts { argBrightness = read n })
  parseArgs ("--gamma":n:rest) opts = parseArgs rest (opts { argGamma = read n })
  parseArgs ("--compress":n:rest) opts = parseArgs rest (opts { argCompress = read n })
  parseArgs ("--scale":n:rest) opts = parseArgs rest (opts { argScale = read n })
  parseArgs ("--output":n:rest) opts = parseArgs rest (opts { output = n })
  parseArgs (name:rest) opts = parseArgs rest (opts { file = name })

ext: hs

max_stars_repo_name: mdibaiee/picedit

max_issues_repo_name: mdibaiee/picedit

max_forks_repo_name: mdibaiee/picedit


text:
{-# LANGUAGE FlexibleContexts #-}

module Lib.Practice.P8 (
    f8
) where

import Numeric.LinearAlgebra
import Lib.Utility

f8 = do
    putStrLn "f8"

y = (2><10) [0.1, 0.05, 0.6, 0.0, 0.05, 0.1, 0.0, 0.1, 0.0, 0.0,   0.1, 0.05, 0.1, 0.0, 0.05, 0.1, 0.0, 0.6, 0.5, 0.0] :: Matrix Double
t = asRow $ fromList [0,0,1,0,0,0,0,0,0,0] :: Matrix Double
meanSquaredError y t = map ((/ 2) . sumElements) (toRows $ (y - t) ** 2)
crossEntropyError y t = map (negate . sumElements) $ toRows $ cmap log (y + delta) * t
    where delta = 1e-7

function1 x = (0.01 * x ** 2) + 0.1 * x
function1_xs = fromList [0,0.1..20.0] :: Vector Double
function1_ys = cmap function1 function1_xs
function2 :: Vector Double -> Double
function2 xs = sumElements $ xs ** 2
function2_initXs = fromList [-3.0, 4.0] :: Vector Double


text:
module HaskovSpec where

import Haskov (fromList,imap,hmatrix,walk,walkFrom,steady,steadyState,statesI)

import System.Random
import Test.Hspec
import qualified Data.Set as Set
import qualified Numeric.LinearAlgebra.Data as Dat
import Data.List (intercalate)

import           Control.Monad (unless,when)

testTransitions =
  [ (("A", "B"), 0.3)
  , (("A", "A"), 0.7)
  , (("B", "C"), 1.0)
  , (("C", "A"), 1.0)
  ]

spec :: Spec
spec =

  describe "haskov" $ do

--    describe "steadyState" $ do
--
--      it "should always contain positive numbers" $ do
--        let
--          transitions = [ (("A", "B"), 1.0), (("B", "C"), 1.0), (("C", "A"), 1.0)]
--          haskov = fromList transitions
--          s = steadyState haskov
--        --mapM (\(s, p) -> s `shouldSatisfy` (> 0))
--        print s
--        True `shouldBe` True

    describe "A haskov walk" $ do

    -- could use quickcheck for better tests

      it "should never do invalid transitions when transitions are ordered" $ do
        let
          transitions = [(("A", "B"), 1.0), (("B", "C"), 1.0), (("C", "A"), 1.0)]
          valid = map fst transitions
          markov = fromList transitions
        gen <- getStdGen
        res <- walk 10 markov gen
        expectOnlyValidTransitions transitions res

      it "should never do invalid transitions when transtions are not ordered" $ do
        let
          transitions = [ (("C", "A"), 1.0), (("A", "B"), 1.0), (("B", "C"), 1.0)]
          valid = map fst transitions
          haskov = fromList transitions
        gen <- getStdGen
        res <- walk 3 haskov gen
        expectOnlyValidTransitions transitions res

    describe "A haskov walkFrom" $ do

      it "starts with the head initial state" $ do
        let
          markov = fromList testTransitions
        gen <- getStdGen
        res <- walkFrom "B" 10 markov gen
        head res `shouldBe` "B"

      it "just some test" $ do
        let
          markov = fromList testTransitions
          index = imap markov
          matrix = hmatrix markov
          start = "A"
        gen <- getStdGen
        res <- walkFrom "B" 10 markov gen
        --putStrLn $ "index: " ++ ( show index)
        --putStrLn $ "matrix: " ++ ( show matrix)
        --putStrLn $ "result from " ++ start ++ ": " ++ ( show res)
        return ()

expectTrue :: HasCallStack => String -> Bool -> Expectation
expectTrue msg b = unless b (expectationFailure msg)

expectOnlyValidTransitions transitions actual = do
  let
    valid = map fst transitions
    actualTransitions = zip actual (drop 1 actual)
    actualUnique = Set.fromList actualTransitions
    invalid = Set.filter (\e -> not $ elem e valid) actualUnique
    len = length invalid
    invalidMsg inv = "invalid transitions encountered: " ++ ( intercalate ", " (Set.toList (Set.map (\(a,b) -> a ++ "->" ++ b) inv)))
  expectTrue (invalidMsg invalid) $ len == 0

max_stars_repo_name: francisdb/haskov

max_issues_repo_name: francisdb/haskov

max_forks_repo_name: francisdb/haskov


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveDataTypeable    #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE LambdaCase            #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternSynonyms       #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeInType            #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE ViewPatterns          #-}

module Backprop.Learn.Model.Function (
  -- * Types
    ParamFunc(..)
  , ParamFuncP, pattern PFP, runParamFuncP
  , FixedFunc, pattern FF, runFixedFunc
  , paramMap, learnParam
  , idLearn
  -- ** Combinators
  , dimapPF, lmapPF, rmapPF, compPF, parPF
  -- *** Chain
  , (.-), nilPF, onlyPF
  -- * Activation functions
  -- | See <https://en.wikipedia.org/wiki/Activation_function>
  --
  -- ** Maps
  -- *** Unparameterized
  , step
  , logistic
  , softsign
  , reLU
  , softPlus
  , bentIdentity
  , siLU
  , softExponential
  , sinc
  , gaussian
  , tanh
  , atan
  , sin
  , vmap
  , vmap'
  -- *** Parameterized
  , liftUniform
  , isru
  , preLU
  , sreLU
  , sreLUPFP
  , eLU
  , isrLU
  , apl
  , aplPFP
  -- ** Mixing
  , softMax
  , maxout
  , kSparse
  ) where

import           Backprop.Learn.Model.Class
import           Control.Category
import           Data.Foldable
import           Data.Proxy
import           Data.Type.Length
import           Data.Type.Mayb hiding                        (type (<$>))
import           Data.Type.Tuple hiding                       (T2(..), T3(..))
import           Data.Typeable
import           GHC.TypeNats
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop hiding (tr)
import           Prelude hiding                               ((.), id)
import           Type.Class.Known
import           Type.Family.List
import qualified Data.Type.Tuple                              as T
import qualified Data.Vector.Sized                            as SV
import qualified Data.Vector.Storable.Sized                   as SVS
import qualified Numeric.LinearAlgebra                        as HU
import qualified Numeric.LinearAlgebra.Static                 as H
import qualified Numeric.LinearAlgebra.Static.Vector          as H

-- | An unparameterized function.  Has a 'Category' instance.
--
-- A @'FixedFunc' a b@ essentially the same as a:
--
-- @
-- forall s. 'Reifies' s 'W' => 'BVar' s a -> 'BVar' s b
-- @
--
-- And the 'FF' pattern (and 'runFixedFunc' extractor) witness this.
--
-- It is usually better to just use the functions directly, with
-- combinators like 'Dimap', 'LMap', 'RMap', 'dimapPF', 'lmapPF', 'rmapPF',
-- etc.  This is just provided to let you work nicely with 'ParamFunc'
-- combinators.
type FixedFunc = ParamFunc 'Nothing

-- | 'FF' and 'runFixedFunc' witness the fact that a @'FixedFunc' a b@ is
-- just a function from @'BVar' s a@ to @'BVar' s b@.
pattern FF :: (forall s. Reifies s W => BVar s a -> BVar s b) -> FixedFunc a b
pattern FF { runFixedFunc } <- (getFF->runFixedFunc) where
    FF f = PF (const f)
{-# COMPLETE FF #-}

getFF :: forall a b. FixedFunc a b -> (forall s. Reifies s W => BVar s a -> BVar s b)
getFF ff = runParamFunc ff N_

-- | Identity model, useful for using with other combinators.
idLearn :: FixedFunc a a
idLearn = FF id

-- | A @'ParamFunc' p a b@ is a parameterized function from @a@ to @b@,
-- potentially with trainable parameter @p@.
--
-- A utility wrapper for a deterministic and stateless model.
newtype ParamFunc p a b =
    PF { runParamFunc :: forall s. Reifies s W => Mayb (BVar s) p -> BVar s a -> BVar s b
       }
  deriving (Typeable)

instance ( Num a, Num b, MaybeC Num p, KnownMayb p)
      => Learn a b (ParamFunc p a b) where
    type LParamMaybe (ParamFunc p a b) = p

    runLearn p = stateless . runParamFunc p

instance Category (ParamFunc p) where
    id = PF (const id)
    f . g = PF $ \p -> runParamFunc f p
                     . runParamFunc g p

-- | Convenient type synonym for a 'ParamFunc' with parameters.
--
-- Mostly made to be easy to construct/deconstruct with 'PFP' and
-- 'runParamFuncP'.
type ParamFuncP p = ParamFunc ('Just p)

pattern PFP :: (forall s. Reifies s W => BVar s p -> BVar s a -> BVar s b)
            -> ParamFuncP p a b
pattern PFP { runParamFuncP } <- (getPFP->(getWF->runParamFuncP))
  where
    PFP f = PF $ \(J_ p) -> f p
{-# COMPLETE PFP #-}

newtype WrapFun p a b = WF { getWF :: forall s. Reifies s W => BVar s p -> BVar s a -> BVar s b }

getPFP :: ParamFuncP p a b -> WrapFun p a b
getPFP pf = WF (\case p -> runParamFunc pf (J_ p))

-- | Create a 'ParamFunc' from any instance of 'Learn' that does not have
-- state.
learnParam
    :: forall l a b. (Learn a b l, NoState l)
    => l
    -> ParamFunc (LParamMaybe l) a b
learnParam l = PF (runLearnStateless l)

-- | 'ParamFuncP' taking a singleton list; meant to be used with '.-'
onlyPF
    :: forall p a b. (KnownMayb p, MaybeC Backprop p)
    => ParamFunc p a b
    -> ParamFuncP (T (MaybeToList p)) a b
onlyPF f = PFP $ \ps -> case knownMayb @p of
                          N_   -> runParamFunc f N_
                          J_ _ -> runParamFunc f (J_ (isoVar tOnly onlyT ps))

-- | Compose two 'ParamFuncP's on lists.
(.-)
    :: forall ps qs a b c. (ListC (Backprop <$> ps), ListC (Backprop <$> qs), Known Length ps)
    => ParamFuncP (T ps) b c
    -> ParamFuncP (T qs) a b
    -> ParamFuncP (T (ps ++ qs)) a c
f .- g = PFP $ \ps -> runParamFuncP f (ps ^^. tTake @ps @qs known)
                    . runParamFuncP g (ps ^^. tDrop @ps @qs known)
infixr 9 .-

-- | The identity of '.-'
nilPF :: ParamFuncP (T '[]) a a
nilPF = id

-- | Pre- and post-compose a 'ParamFunc'
dimapPF
    :: (forall s. Reifies s W => BVar s a -> BVar s b)
    -> (forall s. Reifies s W => BVar s c -> BVar s d)
    -> ParamFunc p b c
    -> ParamFunc p a d
dimapPF f g h = PF $ \p -> g . runParamFunc h p . f

-- | Precompose a 'ParamFunc'
lmapPF
    :: (forall s. Reifies s W => BVar s a -> BVar s b)
    -> ParamFunc p b c
    -> ParamFunc p a c
lmapPF f = dimapPF f id

-- | Postcompose a 'ParamFunc'
rmapPF
    :: (forall s. Reifies s W => BVar s b -> BVar s c)
    -> ParamFunc p a b
    -> ParamFunc p a c
rmapPF = dimapPF id

-- | Compose two 'ParamFunc's sequentially
compPF
    :: forall p q a b c. (MaybeC Backprop p, MaybeC Backprop q, KnownMayb p, KnownMayb q)
    => ParamFunc p a b
    -> ParamFunc q b c
    -> ParamFunc (TupMaybe p q) a c
compPF f g = PF $ \pq ->
    let (p, q) = splitTupMaybe @_ @p @q (\(T2B p' q') -> (p', q')) pq
    in  runParamFunc g q . runParamFunc f p

-- | Compose two 'ParamFunc's in parallel
parPF
    :: forall p q a b c d.
     ( MaybeC Backprop p
     , MaybeC Backprop q
     , KnownMayb p
     , KnownMayb q
     , Backprop a
     , Backprop b
     , Backprop c
     , Backprop d
     )
    => ParamFunc p a c
    -> ParamFunc q b d
    -> ParamFunc (TupMaybe p q) (T.T2 a b) (T.T2 c d)
parPF f g = PF $ \pq (T2B x y) ->
    let (p, q) = splitTupMaybe @_ @p @q (\(T2B p' q') -> (p', q')) pq
    in  T2B (runParamFunc f p x)
            (runParamFunc g q y)

-- TODO: replace all of these with manual ops?

-- | Softmax normalizer
softMax
    :: (KnownNat i, Reifies s W)
    => BVar s (R i)
    -> BVar s (R i)
softMax x = expx / konst (norm_1V expx)
  where
    expx = exp x

-- | Logistic function
--
-- \[
-- \sigma(x) = \frac{1}{1 + e^{-x}}
-- \]
logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))

-- | Binary step / heaviside step
--
-- To use with vectors ('R'), use 'vmap''.
--
-- \[
-- \begin{cases}
-- 0 & \text{for } x < 0 \\
-- 1 & \text{for } x \ge 0
-- \end{cases}
-- \]
step :: (Ord a, Num a) => a -> a
step x | x < 0     = 0
       | otherwise = 1

-- | Softsign activation function
--
-- \[
-- \frac{x}{1 + \lvert x \rvert}
-- \]
softsign :: Fractional a => a -> a
softsign x = x / (1 + abs x)

-- | Inverse square root unit
--
-- \[
-- \frac{x}{\sqrt{1 + \alpha x^2}}
-- \]
--
-- See 'liftUniform' to make this compatible with 'PFP'.
--
-- You can also just use this after partially applying it, to fix the
-- parameter (and not have it trained).
isru
    :: Floating a
    => a        -- ^ α (scaling parameter)
    -> a        -- ^ x
    -> a
isru α x = x / sqrt (1 + α * x * x)

-- | Rectified linear unit.
--
-- To use with vectors ('R'), use 'vmap''.
--
--
-- \[
-- \max(0,x)
-- \]
--
-- @
-- 'reLU' = 'preLU' 0
-- @
--
reLU :: (Num a, Ord a) => a -> a
reLU x | x < 0     = 0
       | otherwise = x

-- | Parametric rectified linear unit
--
-- To use with vectors ('R'), use 'vmap''.
--
-- If scaling parameter is a fixed (and not learned) parameter, this is
-- typically called a leaky recitified linear unit (typically with
-- α = 0.01).
--
-- To use as a learned parameter:
--
-- @
-- 'vmap' . 'preLU' :: 'BVar' s Double -> 'BVar' s ('R' n) -> BVar s (R n)
-- @
--
-- This can be give directly to 'PFP'.
--
-- To fix the paramater ("leaky"), just partially apply a parameter:
--
-- @
-- 'preLU' 0.01           :: 'BVar' s ('R' n) -> BVar s (R n)
-- preLU ('realToFrac' α) :: BVar s (R n) -> BVar s (R n)
-- @
--
-- See also 'rreLU'.
--
-- \[
-- \begin{cases}
-- \alpha x & \text{for } x < 0 \\
-- x        & \text{for } x \ge 0
-- \end{cases}
-- \]
preLU
    :: (Num a, Ord a)
    => a        -- ^ α (scaling parameter)
    -> a        -- ^ x
    -> a
preLU α x | x < 0     = α * x
          | otherwise = x

-- | Exponential linear unit
--
-- To use with vectors ('R'), use 'vmap''.
--
-- To use as a learned parameter:
--
-- @
-- 'vmap' . 'eLU' :: 'BVar' s Double -> 'BVar' s ('R' n) -> BVar s (R n)
-- @
--
-- This can be give directly to 'PFP'.
--
-- To fix the paramater, just partially apply a parameter:
--
-- @
-- 'vmap'' ('eLU' 0.01) :: 'BVar' s ('R' n) -> BVar s (R n)
-- @
--
-- \[
-- \begin{cases}
-- \alpha (e^x - 1) & \text{for } x < 0 \\
-- x                & \text{for } x \ge 0
-- \end{cases}
-- \]
eLU :: (Floating a, Ord a)
    => a        -- ^ α (scaling parameter)
    -> a        -- ^ x
    -> a
eLU α x | x < 0     = α * (exp x - 1)
        | otherwise = x

-- | S-shaped rectified linear activiation unit
--
-- See 'sreLUPFP' for an uncurried and uniformly lifted version usable with
-- 'PFP'.
--
-- \[
-- \begin{cases}
-- t_l + a_l (x - t_l) & \text{for } x \le t_l  \\
-- x                   & \text{for } t_l < x < t_r  \\
-- t_r + a_r (x - t_r) & \text{for } x \ge t_r
-- \end{cases}
-- \]
sreLU
    :: (Num a, Ord a)
    => a        -- ^ @t_l@
    -> a        -- ^ @a_l@
    -> a        -- ^ @t_r@
    -> a        -- ^ @a_l@
    -> a        -- ^ @x@
    -> a
sreLU tl al tr ar x
    | x < tl    = tl + al * (x - tl)
    | x > tr    = tr + ar * (x - tr)
    | otherwise = x

-- | An uncurried and uniformly lifted version of 'sreLU' directly usable
-- with 'PFP'.
sreLUPFP
    :: (KnownNat n, Reifies s W)
    => BVar s (T.T2 (T.T2 Double Double) (T.T2 Double Double))
    -> BVar s (R n)
    -> BVar s (R n)
sreLUPFP (T2B (T2B tl al) (T2B tr ar)) = vmap (sreLU tl al tr ar)

-- | Inverse square root linear unit
--
-- To use with vectors ('R'), use 'vmap''.
--
-- \[
-- \begin{cases}
-- \frac{x}{1 + \alpha x^2} & \text{for } x < 0 \\
-- x                        & \text{for } x \ge 0
-- \end{cases}
-- \]
isrLU
    :: (Floating a, Ord a)
    => a        -- ^ α
    -> a        -- ^ x
    -> a
isrLU α x
    | x < 0     = x / sqrt (1 + α * x * x)
    | otherwise = x

-- | Adaptive piecewise linear activation unit
--
-- See 'aplPFP' for an uncurried version usable with 'PFP'.
--
-- \[
-- \max(0, x_i) + \sum_j^M a_i^j \max(0, -x_i + b_i^j)
-- \]
apl :: (KnownNat n, KnownNat m, Reifies s W)
    => BVar s (L n m)     -- ^ a
    -> BVar s (L n m)     -- ^ b
    -> BVar s (R m)       -- ^ x
    -> BVar s (R m)
apl as bs x = vmap' (max 0) x
            + sum (toRows (as * (bs - fromRows (SV.replicate x))))

-- | 'apl' uncurried, to be directly usable with 'PFP'.
aplPFP
    :: (KnownNat n, KnownNat m, Reifies s W)
    => BVar s (T.T2 (L n m) (L n m))
    -> BVar s (R m)
    -> BVar s (R m)
aplPFP (T2B a b) = apl a b

-- | SoftPlus
--
-- \[
-- \ln(1 + e^x)
-- \]
softPlus :: Floating a => a -> a
softPlus x = log (1 + exp x)

-- | Bent identity
--
-- \[
-- \frac{\sqrt{x^2 + 1} - 1}{2} + x
-- \]
bentIdentity :: Floating a => a -> a
bentIdentity x = (sqrt (x * x + 1) - 1) / 2 + x

-- | Sigmoid-weighted linear unit.  Multiply 'logistic' by its input.
--
-- \[
-- x \sigma(x)
-- \]
siLU :: Floating a => a -> a
siLU x = x * logistic x

-- | SoftExponential
--
-- To use with vectors ('R'), use 'vmap''.
--
-- \[
-- \begin{cases}
-- - \frac{\ln(1 - \alpha (x + \alpha))}{\alpha} & \text{for } \alpha < 0  \\
-- x                                             & \text{for } \alpha = 0  \\
-- \frac{e^{\alpha x} - 1}{\alpha} + \alpha      & \text{for } \alpha > 0
-- \end{cases}
-- \]
softExponential
    :: (Floating a, Ord a)
    => a            -- ^ α
    -> a            -- ^ x
    -> a
softExponential α x = case compare α x of
    LT -> - log (1 - α * (x + α)) / α
    EQ -> x
    GT -> (exp (α * x) - 1) / α + α

-- | Sinc
--
-- \[
-- \begin{cases}
-- 1                 & \text{for } x = 0  \\
-- \frac{\sin(x)}{x} & \text{for } x \ne 0
-- \end{cases}
-- \]
sinc :: (Floating a, Eq a) => a -> a
sinc 0 = 1
sinc x = sin x / x

-- | Gaussian
--
-- \[
-- e^{-x^2}
-- \]
gaussian :: Floating a => a -> a
gaussian x = exp (- (x * x))

-- | Maximum of vector.
--
-- Compare to 'norm_InfV', which gives the maximum absolute value.
maxout :: (KnownNat n, Reifies s W) => BVar s (R n) -> BVar s Double
maxout = liftOp1 . op1 $ \x ->
    let n = HU.maxElement . H.extract $ x
    in  ( n
        , \d -> H.vecR . SVS.map (\e -> if e == n then d else 0) . H.rVec $ x
        )

-- | Usable with functions like '*', 'isru', etc. to turn them into a form
-- usable with 'PFP':
--
-- @
-- 'liftUniform' ('*')  :: 'BVar' s 'Double' -> BVar s ('R' n) -> BVar s (R n)
-- liftUniform 'isru' :: BVar s Double -> BVar s (R n) -> BVar s (R n)
-- @
--
-- Basically turns a parmaeterized function on individual elements of
-- into one that shares the same parameter across all elements of the
-- vector.
liftUniform
    :: (Reifies s W, KnownNat n)
    => (BVar s (R n) -> r)
    -> BVar s Double
    -> r
liftUniform f = f . konst

-- | Utility function to make a 'ParamFunc' that maps a parameterized
-- function over every item in the 'R'.  The parameter is shared across the
-- entire map, and trained.
paramMap
    :: KnownNat i
    => (forall s. Reifies s W => Mayb (BVar s) p -> BVar s Double -> BVar s Double)
    -> ParamFunc p (R i) (R i)
paramMap f = PF (vmap . f)

-- -- TODO: vmap can do better.

-- | Keep only the top @k@ values, and zero out all of the rest.
--
-- Useful for postcomposing in between layers (with a logistic function
-- before) to encourage the number of "activated" neurons is kept to be
-- around @k@.  Used in k-Sprase autoencoders (see 'KAutoencoder').
--
-- <http://www.ericlwilkinson.com/blog/2014/11/19/deep-learning-sparse-autoencoders>
kSparse
    :: forall n s. (Reifies s W, KnownNat n)
    => Int                      -- ^ number of items to keep
    -> BVar s (R n)
    -> BVar s (R n)
kSparse k = liftOp1 . op1 $ \xs ->
    let xsSort = HU.sortVector (H.extract xs)
        thresh = xsSort `HU.atIndex` (n - k)
        mask   = H.dvmap (\x -> if x >= thresh then 1 else 0) xs
    in  ( H.dvmap (\x -> if x >= thresh then x else 0) xs
        , (* mask)
        )
  where
    n = fromIntegral $ natVal (Proxy @n)

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE RecordWildCards     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeInType          #-}
{-# LANGUAGE TypeOperators       #-}

#if __GLASGOW_HASKELL__ >= 805
{-# LANGUAGE NoStarIsType #-}
{-# LANGUAGE ExplicitNamespaces #-}
#endif

module Eigen.Matrix.Mutable
  ( -- * Types
    MMatrix(..)
  , MMatrixXf
  , MMatrixXd
  , MMatrixXcf
  , MMatrixXcd
  , IOMatrix
  , STMatrix

    -- * Construction
  , new
  , replicate

    -- * Indexing
  , read
  , write

    -- * Modification
  , set
  , copy

    -- * Modification with Pointers
  , unsafeWith

    -- * Conversion to Vectors
  , vals
  , fromVector
  ) where

import Eigen.Internal
  ( Elem
  , C(..)
  , natToInt
  , Row(..)
  , Col(..)
  )

import Prelude hiding (replicate, read)
import Control.Monad.Primitive (PrimMonad(..))
import Data.Complex (Complex)
import Data.Kind (Type)
import Foreign.C.Types (CInt)
import Foreign.Ptr (Ptr)
import GHC.Exts (RealWorld)
import GHC.TypeLits (Nat, type (*), type (<=), KnownNat)
import qualified Data.Vector.Storable.Mutable as VSM

-- | A mutable matrix. See 'Eigen.Matrix.Matrix' for
--   details about matrix layout.
newtype MMatrix :: Nat -> Nat -> Type -> Type -> Type where
  MMatrix :: Vec (n * m) s a -> MMatrix n m s a

-- | Used internally to track the size and corresponding C type of the matrix.
newtype Vec :: Nat -> Type -> Type -> Type where
  Vec :: VSM.MVector s (C a) -> Vec n s a

-- | Alias for single precision mutable matrix
type MMatrixXf  n m s = MMatrix n m s Float
-- | Alias for double precision mutable matrix
type MMatrixXd  n m s = MMatrix n m s Double
-- | Alias for single precision mutable matrix of complex numbers
type MMatrixXcf n m s = MMatrix n m s (Complex Float)
-- | Alias for double precision mutable matrix of complex numbers
type MMatrixXcd n m s = MMatrix n m s (Complex Double)

-- | A mutable matrix where the state token is specialised to 'RealWorld'.
type IOMatrix n m a   = MMatrix n m RealWorld a
-- | This type does not differ from MSparseMatrix, but might be desirable for readability.
type STMatrix n m s a = MMatrix n m s a

-- | Create a mutable matrix of the given size and fill it with 0 as an initial value.
new :: (PrimMonad p, Elem a, KnownNat n, KnownNat m) => p (MMatrix n m (PrimState p) a)
{-# INLINE new #-}
new = replicate 0

-- | Create a mutable matrix of the given size and fill it with an initial value.
replicate :: forall n m p a. (PrimMonad p, Elem a, KnownNat n, KnownNat m) => a -> p (MMatrix n m (PrimState p) a)
{-# INLINE replicate #-}
replicate !val = do
  let !mm_rows = natToInt @n
      !mm_cols = natToInt @m
      !cval    = toC val
  _vals <- VSM.replicate (mm_rows * mm_cols) cval
  pure (MMatrix $! Vec $! _vals)

-- | Set all elements of the matrix to a given value.
set :: (PrimMonad p, Elem a) => MMatrix n m (PrimState p) a -> a -> p ()
{-# INLINE set #-}
set (MMatrix (Vec !vec)) !val =
  let !cval = toC val
  in VSM.set vec cval

-- | Copy a matrix.
copy :: (PrimMonad p, Elem a) => MMatrix n m (PrimState p) a -> MMatrix n m (PrimState p) a -> p ()
{-# INLINE copy #-}
copy (MMatrix (Vec m1)) (MMatrix (Vec m2)) = VSM.unsafeCopy m1 m2

-- | Yield the element at the given position.
read :: forall n m p a r c. (PrimMonad p, Elem a, KnownNat n, KnownNat r, KnownNat c, r <= n, c <= m)
  => Row r -> Col c -> MMatrix n m (PrimState p) a -> p a
{-# INLINE read #-}
read _ _ (MMatrix (Vec m)) =
  let !row = natToInt @r
      !col = natToInt @c
      !mm_rows = natToInt @n
  in VSM.unsafeRead m (col * mm_rows + row) >>= \ !val -> let !cval = fromC val in pure cval

-- | Replace the element at the given position.
write :: forall n m p a r c. (PrimMonad p, Elem a, KnownNat n, KnownNat r, KnownNat c, r <= n, c <= m)
  => Row r -> Col c -> MMatrix n m (PrimState p) a -> a -> p ()
{-# INLINE write #-}
write _ _ (MMatrix (Vec m)) !val =
  let !row = natToInt @r
      !col = natToInt @c
      !mm_rows = natToInt @n
      !cval = toC val
  in VSM.unsafeWrite m (col * mm_rows + row) cval

-- | Pass a pointer to the matrix's data to the IO action.
--   Modifying dat through the pointer is unsafe if the matrix
--   could have been frozen before the modification.
unsafeWith :: forall n m a b. (KnownNat n, KnownNat m, Elem a) => IOMatrix n m a -> (Ptr (C a) -> CInt -> CInt -> IO b) -> IO b
{-# INLINE unsafeWith #-}
unsafeWith (MMatrix (Vec m)) f =
  let !cmm_rows = toC $! natToInt @n
      !cmm_cols = toC $! natToInt @m
  in VSM.unsafeWith m $ \p -> f p cmm_rows cmm_cols

-- | Return a mutable storable 'VSM.MVector' of the corresponding C types to one's mutable matrix.
vals :: MMatrix n m s a -> VSM.MVector s (C a)
{-# INLINE vals #-}
vals (MMatrix (Vec x)) = x

-- | Create a mutable matrix from a mutable storable 'VSM.MVector'.
fromVector :: Elem a => VSM.MVector s (C a) -> MMatrix n m s a
{-# INLINE fromVector #-}
fromVector x = MMatrix (Vec x)

max_stars_repo_name: nilsalex/eigen

max_issues_repo_name: nilsalex/eigen

max_forks_repo_name: nilsalex/eigen


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.ChiSquared
-- Copyright : (c) 2010 Alexey Khudyakov
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The chi-squared distribution. This is a continuous probability
-- distribution of sum of squares of k independent standard normal
-- distributions. It's commonly used in statistical tests
module Statistics.Distribution.ChiSquared (
          ChiSquared
        , chiSquaredNDF
        -- * Constructors
        , chiSquared
        , chiSquaredE
        ) where

import Control.Applicative
import Data.Data             (Data, Typeable)
import GHC.Generics          (Generic)
import Numeric.SpecFunctions ( incompleteGamma,invIncompleteGamma,logGamma,digamma)
import Numeric.MathFunctions.Constants (m_neg_inf)
import qualified System.Random.MWC.Distributions as MWC

import qualified Statistics.Distribution         as D
import Statistics.Internal

-- | Chi-squared distribution
newtype ChiSquared = ChiSquared
  { chiSquaredNDF :: Int
    -- ^ Get number of degrees of freedom
  }
  deriving (Eq, Typeable, Data, Generic)

instance Show ChiSquared where
  showsPrec i (ChiSquared n) = defaultShow1 "chiSquared" n i
instance Read ChiSquared where
  readPrec = defaultReadPrecM1 "chiSquared" chiSquaredE

-- | Construct chi-squared distribution. Number of degrees of freedom
--   must be positive.
chiSquared :: Int -> ChiSquared
chiSquared n = maybe (error $ errMsg n) id $ chiSquaredE n

-- | Construct chi-squared distribution. Number of degrees of freedom
--   must be positive.
chiSquaredE :: Int -> Maybe ChiSquared
chiSquaredE n
  | n <= 0    = Nothing
  | otherwise = Just (ChiSquared n)

errMsg :: Int -> String
errMsg n = "Statistics.Distribution.ChiSquared.chiSquared: N.D.F. must be positive. Got " ++ show n

instance D.Distribution ChiSquared where
  cumulative = cumulative

instance D.ContDistr ChiSquared where
  density chi x
    | x <= 0    = 0
    | otherwise = exp $ log x * (ndf2 - 1) - x2 - logGamma ndf2 - log 2 * ndf2
    where
      ndf  = fromIntegral $ chiSquaredNDF chi
      ndf2 = ndf/2
      x2   = x/2

  logDensity chi x
    | x <= 0    = m_neg_inf
    | otherwise = log x * (ndf2 - 1) - x2 - logGamma ndf2 - log 2 * ndf2
    where
      ndf  = fromIntegral $ chiSquaredNDF chi
      ndf2 = ndf/2
      x2   = x/2

  quantile = quantile

instance D.Mean ChiSquared where
    mean (ChiSquared ndf) = fromIntegral ndf

instance D.Variance ChiSquared where
    variance (ChiSquared ndf) = fromIntegral (2*ndf)

instance D.MaybeMean ChiSquared where
    maybeMean = Just . D.mean

instance D.MaybeVariance ChiSquared where
    maybeStdDev   = Just . D.stdDev
    maybeVariance = Just . D.variance

instance D.Entropy ChiSquared where
  entropy (ChiSquared ndf) =
    let kHalf = 0.5 * fromIntegral ndf in
    kHalf
    + log 2
    + logGamma kHalf
    + (1-kHalf) * digamma kHalf

instance D.MaybeEntropy ChiSquared where
  maybeEntropy = Just . D.entropy

instance D.ContGen ChiSquared where
    genContVar (ChiSquared n) = MWC.chiSquare n

cumulative :: ChiSquared -> Double -> Double
cumulative chi x
  | x <= 0    = 0
  | otherwise = incompleteGamma (ndf/2) (x/2)
  where
    ndf = fromIntegral $ chiSquaredNDF chi

quantile :: ChiSquared -> Double -> Double
quantile (ChiSquared ndf) p
  | p == 0         = 0
  | p == 1         = 1/0
  | p > 0 && p < 1 = 2 * invIncompleteGamma (fromIntegral ndf / 2) p
  | otherwise      =
    error $ "Statistics.Distribution.ChiSquared.quantile: p must be in [0,1] range. Got: "++show p

max_stars_repo_name: vmchale/statistics

max_issues_repo_name: vmchale/statistics

max_forks_repo_name: vmchale/statistics


text:
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeOperators       #-}

#if MIN_VERSION_base(4,12,0)
{-# LANGUAGE NoStarIsType         #-}
#endif

-- |
-- Module      : Numeric.LinearAlgebra.Static.Vector
-- Copyright   : (c) Justin Le 2018
-- License     : BSD3
--
-- Maintainer  : justin@jle.im
-- Stability   : experimental
-- Portability : non-portable
--
-- Conversions between statically sized types in
-- "Numeric.LinearAlgebra.Static" from /hmatrix/ and /vector-sized/.
--
-- This module is intentionally minimal, exporting only functions that
-- cannot be written without "unsafe" operations.  With these, however, you
-- can easily write other useful combinators by using type-safe operations
-- like 'fmap', 'VS.map', 'liftA2', 'Data.Vector.Generic.Sized.convert',
-- etc.
--

module Numeric.LinearAlgebra.Static.Vector (
  -- * Vector
  -- ** Real
    rVec
  , grVec
  , vecR
  , gvecR
  -- ** Complex
  , cVec
  , gcVec
  , vecC
  , gvecC
  -- * Matrix
  -- ** Real
  , lRows
  , rowsL
  , lCols
  , colsL
  , lVec
  , glVec
  , vecL
  , gvecL
  -- ** Complex
  , mRows
  , rowsM
  , mCols
  , colsM
  , mVec
  , gmVec
  , vecM
  , gvecM
  ) where

import           Data.Foldable
import           Data.Proxy
import           GHC.TypeLits
import           Unsafe.Coerce
import qualified Data.Vector                  as UV
import qualified Data.Vector.Generic          as UVG
import qualified Data.Vector.Generic.Sized    as VG
import qualified Data.Vector.Sized            as V
import qualified Data.Vector.Storable.Sized   as VS
import qualified Numeric.LinearAlgebra        as HU
import qualified Numeric.LinearAlgebra.Static as H

-- | Convert an /hmatrix/ vector (parameterized by its lenth) to
-- a /vector-sized/ storable vector of 'Double's.
--
-- This is normally /O(1)/, but will be /O(n)/ if the 'H.R' was contructed
-- with 'H.konst' or any other replicated-value constructor (like literals
-- and 'fromInteger'/'fromRational').
rVec :: KnownNat n => H.R n -> VS.Vector n H.ℝ
rVec = unsafeCoerce . H.extract

-- | 'rVec', but generalized to work for all types of sized vectors.
--
-- Usually /O(n)/, but if using this with storable vectors, should have the
-- same characteristics as 'rVec' due to rewrite rules.
--
-- @since 0.1.3.0
grVec :: (KnownNat n, UVG.Vector v H.ℝ) => H.R n -> VG.Vector v n H.ℝ
grVec = VG.convert . rVec
{-# NOINLINE[1] grVec #-}
{-# RULES "grVec" grVec = rVec #-}

-- | Convert a /vector-sized/ storable vector to an /hmatrix/ vector
-- (parameterized by its lenth).
--
-- /O(1)/
vecR :: VS.Vector n H.ℝ -> H.R n
vecR = unsafeCoerce

-- | 'vecR', but generalized to work for all types of sized vectors.
--
-- Usually /O(n)/, but if using this with storable vectors, should be /O(1)/
-- due to rewrite rules (but don't rely on this).
--
-- @since 0.1.3.0
gvecR :: UVG.Vector v H.ℝ => VG.Vector v n H.ℝ -> H.R n
gvecR = vecR . VG.convert
{-# NOINLINE[1] gvecR #-}
{-# RULES "gvecR" gvecR = vecR #-}

-- | Convert an /hmatrix/ complex vector (parameterized by its lenth) to
-- a /vector-sized/ storable vector of 'Complex Double's, preserving the
-- length in the type.
--
-- This is normally /O(1)/, but will be /O(n)/ if the 'H.C' was contructed
-- with 'H.konst' or any other replicated-value constructor (like literals
-- and 'fromInteger'/'fromRational').
cVec :: KnownNat n => H.C n -> VS.Vector n H.ℂ
cVec = unsafeCoerce . H.extract

-- | 'cVec', but generalized to work for all types of sized vectors.
--
-- Usually /O(n)/, but if using this with storable vectors, should have the
-- same characteristics as 'cVec' due to rewrite rules.
--
-- @since 0.1.3.0
gcVec :: (KnownNat n, UVG.Vector v H.ℂ) => H.C n -> VG.Vector v n H.ℂ
gcVec = VG.convert . cVec
{-# NOINLINE[1] gcVec #-}
{-# RULES "gcVec" gcVec = cVec #-}

-- | Convert a /vector-sized/ storable vector to an /hmatrix/ complex
-- vector (parameterized by its lenth), preserving the length in the type.
--
-- /O(1)/
vecC :: VS.Vector n H.ℂ -> H.C n
vecC = unsafeCoerce

-- | 'vecC', but generalized to work for all types of sized vectors.
--
-- Usually /O(n)/, but if using this with storable vectors, should be /O(1)/
-- due to rewrite rules (but don't rely on this).
--
-- @since 0.1.3.0
gvecC :: UVG.Vector v H.ℂ => VG.Vector v n H.ℂ -> H.C n
gvecC = vecC . VG.convert
{-# NOINLINE[1] gvecC #-}
{-# RULES "gvecC" gvecC = vecC #-}

-- | Split an /hmatrix/ matrix (parameterized by its dimensions) to
-- a /vector-sized/ boxed vector of its rows (as /hmatrix/ vectors).
--
-- This is normally /O(m*n)/, but can sometimes be /O(m)/ depending on the
-- representation of the 'H.L' being used.
lRows
    :: (KnownNat m, KnownNat n)
    => H.L m n
    -> V.Vector m (H.R n)
lRows = unsafeCoerce
      . UV.fromList
      . HU.toRows
      . H.extract

-- | Join together a /vector-sized/ boxed vector of /hmatrix/ vectors to an
-- /hmatrix/ matrix as its rows.
--
-- /O(m*n)/
rowsL
    :: forall m n. KnownNat n
    => V.Vector m (H.R n)
    -> H.L m n
rowsL = (unsafeCoerce :: HU.Matrix Double -> H.L m n)
      . HU.fromRows
      . map H.extract
      . toList

-- | Split an /hmatrix/ matrix (parameterized by its dimensions) to
-- a /vector-sized/ boxed vector of its columns (as /hmatrix/ vectors).
--
-- This is normally /O(m*n)/, but can sometimes be /O(n)/ depending on the
-- representation of the 'H.L' being used.
lCols
    :: forall m n. (KnownNat m, KnownNat n)
    => H.L m n
    -> V.Vector n (H.R m)
lCols = unsafeCoerce
      . UV.fromList
      . HU.toColumns
      . H.extract

-- | Join together a /vector-sized/ boxed vector of /hmatrix/ vectors to an
-- /hmatrix/ matrix as its columns.
--
-- /O(m*n)/
colsL
    :: forall m n. KnownNat m
    => V.Vector n (H.R m)
    -> H.L m n
colsL = (unsafeCoerce :: HU.Matrix Double -> H.L m n)
      . HU.fromColumns
      . map H.extract
      . toList

-- | Split an /hmatrix/ complex matrix (parameterized by its dimensions) to
-- a /vector-sized/ boxed vector of its rows (as /hmatrix/ complex
-- vectors).
--
-- This is normally /O(m*n)/, but can sometimes be /O(m)/ depending on the
-- representation of the 'H.C' being used.
mRows
    :: forall m n. (KnownNat m, KnownNat n)
    => H.M m n
    -> V.Vector m (H.C n)
mRows = unsafeCoerce
      . UV.fromList
      . HU.toRows
      . H.extract

-- | Join together a /vector-sized/ boxed vector of /hmatrix/ complex
-- vectors to an /hmatrix/ complex matrix as its rows.
--
-- /O(m*n)/
rowsM
    :: forall m n. KnownNat n
    => V.Vector m (H.C n)
    -> H.M m n
rowsM = (unsafeCoerce :: HU.Matrix H.ℂ -> H.M m n)
      . HU.fromRows
      . map H.extract
      . toList

-- | Split an /hmatrix/ complex matrix (parameterized by its dimensions) to
-- a /vector-sized/ boxed vector of its columns (as /hmatrix/ complex
-- vectors).
--
-- This is normally /O(m*n)/, but can sometimes be /O(n)/ depending on the
-- representation of the 'H.C' being used.
mCols
    :: forall m n. (KnownNat m, KnownNat n)
    => H.M m n
    -> V.Vector n (H.C m)
mCols = unsafeCoerce
      . UV.fromList
      . HU.toColumns
      . H.extract

-- | Join together a /vector-sized/ boxed vector of /hmatrix/ complex
-- vectors to an /hmatrix/ complex matrix as its columns.
--
-- /O(m*n)/
colsM
    :: forall m n. KnownNat m
    => V.Vector n (H.C m)
    -> H.M m n
colsM = (unsafeCoerce :: HU.Matrix H.ℂ -> H.M m n)
      . HU.fromColumns
      . map H.extract
      . toList

-- | Shape a /vector-sized/ storable vector of elements into an /hmatrix/
-- matrix.
--
-- /O(1)/
--
-- @since 0.1.1.0
vecL
    :: forall m n. KnownNat n
    => VS.Vector (m * n) H.ℝ
    -> H.L m n
vecL = (unsafeCoerce :: HU.Matrix H.ℝ -> H.L m n)
     . HU.reshape (fromIntegral (natVal (Proxy @n)))
     . unsafeCoerce

-- | 'vecL', but generalized to work for all types of sized vectors.
--
-- Usually /O(n)/, but if using this with storable vectors, should be /O(1)/
-- due to rewrite rules (but don't rely on this).
--
-- @since 0.1.3.0
gvecL
    :: (KnownNat n, UVG.Vector v H.ℝ)
    => VG.Vector v (m * n) H.ℝ
    -> H.L m n
gvecL = vecL . VG.convert
{-# NOINLINE[1] gvecL #-}
{-# RULES "gvecL" gvecL = vecL #-}

-- | Flatten an /hmatrix/ matrix into a /vector-sized/ storable vector of
-- its items.
--
-- This is normally /O(m*n)/, but can sometimes be /O(1)/ depending on the
-- representation of the 'H.L' being used.
--
-- @since 0.1.1.0
lVec
    :: forall m n. (KnownNat m, KnownNat n)
    => H.L m n
    -> VS.Vector (m * n) H.ℝ
lVec = unsafeCoerce
     . HU.flatten
     . H.extract

-- | 'lVec', but generalized to work for all types of sized vectors.
--
-- Usually /O(m*n)/, but if using this with storable vectors, should have the
-- same characteristics as 'lVec' due to rewrite rules.
--
-- @since 0.1.3.0
glVec
    :: (KnownNat m, KnownNat n, UVG.Vector v H.ℝ)
    => H.L m n
    -> VG.Vector v (m * n) H.ℝ
glVec = VG.convert . lVec
{-# NOINLINE[1] glVec #-}
{-# RULES "glVec" glVec = lVec #-}

-- | Shape a /vector-sized/ storable vector of elements into an /hmatrix/
-- complex matrix.
--
-- /O(1)/
--
-- @since 0.1.1.0
vecM
    :: forall m n. KnownNat n
    => VS.Vector (m * n) H.ℂ
    -> H.M m n
vecM = (unsafeCoerce :: HU.Matrix H.ℂ -> H.M m n)
     . HU.reshape (fromIntegral (natVal (Proxy @n)))
     . unsafeCoerce

-- | 'vecM', but generalized to work for all types of sized vectors.
--
-- Usually /O(m*n)/, but if using this with storable vectors, should be /O(1)/
-- due to rewrite rules (but don't rely on this).
--
-- @since 0.1.3.0
gvecM
    :: (KnownNat n, UVG.Vector v H.ℂ)
    => VG.Vector v (m * n) H.ℂ
    -> H.M m n
gvecM = vecM . VG.convert
{-# NOINLINE[1] gvecM #-}
{-# RULES "gvecM" gvecM = vecM #-}

-- | Flatten an /hmatrix/ complex matrix into a /vector-sized/ storable
-- vector of its items.
--
-- This is normally /O(m*n)/, but can sometimes be /O(1)/ depending on the
-- representation of the 'H.M' being used.
--
-- @since 0.1.1.0
mVec
    :: forall m n. (KnownNat m, KnownNat n)
    => H.M m n
    -> VS.Vector (m * n) H.ℂ
mVec = unsafeCoerce
     . HU.flatten
     . H.extract

-- | 'mVec', but generalized to work for all types of sized vectors.
--
-- Usually /O(m*n)/, but if using this with storable vectors, should have the
-- same characteristics as 'mVec' due to rewrite rules.
--
-- @since 0.1.3.0
gmVec
    :: (KnownNat m, KnownNat n, UVG.Vector v H.ℂ)
    => H.M m n
    -> VG.Vector v (m * n) H.ℂ
gmVec = VG.convert . mVec
{-# NOINLINE[1] gmVec #-}
{-# RULES "gmVec" gmVec = mVec #-}

max_stars_repo_name: mstksg/hmatrix-vector-sized

max_issues_repo_name: mstksg/hmatrix-vector-sized

max_forks_repo_name: mstksg/hmatrix-vector-sized


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MagicHash             #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PolyKinds             #-}
{-# LANGUAGE RecordWildCards       #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UnboxedTuples         #-}
{-# LANGUAGE UndecidableInstances  #-}
module Numeric.Matrix.Internal
  ( MatrixTranspose (..)
  , SquareMatrix (..)
  , MatrixDeterminant (..)
  , MatrixInverse (..)
  , Matrix
  , HomTransform4 (..)
  , Mat22f, Mat23f, Mat24f
  , Mat32f, Mat33f, Mat34f
  , Mat42f, Mat43f, Mat44f
  , Mat22d, Mat23d, Mat24d
  , Mat32d, Mat33d, Mat34d
  , Mat42d, Mat43d, Mat44d
  , mat22, mat33, mat44
  , (%*)
  ) where

import GHC.Base
import Numeric.DataFrame.Contraction        ((%*))
import Numeric.DataFrame.Internal.PrimArray
import Numeric.DataFrame.Type
import Numeric.Dimensions
import Numeric.Scalar.Internal
import Numeric.Vector.Internal

-- | Alias for DataFrames of rank 2
type Matrix (t :: l) (n :: k) (m :: k) = DataFrame t '[n,m]

class MatrixTranspose t (n :: k) (m :: k) where
    -- | Transpose Mat
    transpose :: Matrix t n m -> Matrix t m n

class SquareMatrix t (n :: Nat) where
    -- | Mat with 1 on diagonal and 0 elsewhere
    eye :: Matrix t n n
    -- | Put the same value on the Mat diagonal, 0 otherwise
    diag :: Scalar t -> Matrix t n n
    -- | Sum of diagonal elements
    trace :: Matrix t n n -> Scalar t

class MatrixDeterminant t (n :: Nat) where
    -- | Determinant of  Mat
    det :: Matrix t n n -> Scalar t

class MatrixInverse t (n :: Nat) where
    -- | Matrix inverse
    inverse :: Matrix t n n -> Matrix t n n

-- | Operations on 4x4 transformation matrices and vectors in homogeneous coordinates.
--   All angles are specified in radians.
--
--   Note: since version 2 of @easytensor@, DataFrames and matrices are row-major.
--         A good SIMD implementation may drastically improve performance
--         of 4D vector-matrix products of the form @v %* m@, but not so much
--         for products of the form @m %* v@ (due to memory layout).
--         Thus, all operations here assume the former form to benefit more from
--         SIMD in future.
class HomTransform4 t where
    -- | Create a translation matrix from a vector.  The 4th coordinate is ignored.
    --
    --   If @p ! 3 == 1@ and @v ! 3 == 0@, then
    --
    --   > p %* translate4 v == p + v
    --
    translate4  :: Vector t 4 -> Matrix t 4 4
    -- | Create a translation matrix from a vector.
    --
    --   If @p ! 3 == 1@, then
    --
    --   > p %* translate3 v == p + toHomVector v
    --
    translate3  :: Vector t 3 -> Matrix t 4 4
    -- | Rotation matrix for a rotation around the X axis, angle is given in radians.
    --   e.g. @p %* rotateX (pi/2)@ rotates point @p@ around @Ox@ by 90 degrees.
    rotateX     :: t -> Matrix t 4 4
    -- | Rotation matrix for a rotation around the Y axis, angle is given in radians.
    --   e.g. @p %* rotateY (pi/2)@ rotates point @p@ around @Oy@ by 90 degrees.
    rotateY     :: t -> Matrix t 4 4
    -- | Rotation matrix for a rotation around the Z axis, angle is given in radians.
    --   e.g. @p %* rotateZ (pi/2)@ rotates point @p@ around @Oz@ by 90 degrees.
    rotateZ     :: t -> Matrix t 4 4
    -- | Rotation matrix for a rotation around an arbitrary normalized vector
    --   e.g. @p %* rotate (pi/2) v@ rotates point @p@ around @v@ by 90 degrees.
    rotate      :: Vector t 3 -> t -> Matrix t 4 4
    -- | Rotation matrix from the Euler angles roll (axis @Z@), yaw (axis @Y'@), and pitch (axis @X''@).
    --   This order is known as Tait-Bryan angles (@Z-Y'-X''@ intrinsic rotations), or nautical angles, or Cardan angles.
    --
    --   > rotateEuler pitch yaw roll == rotateZ roll %* rotateY yaw %* rotateX pitch
    --
    --   https://en.wikipedia.org/wiki/Euler_angles#Conventions_2
    rotateEuler :: t -- ^ pitch (axis @X''@)
                -> t -- ^ yaw (axis @Y'@)
                -> t -- ^ roll (axis @Z@)
                -> Matrix t 4 4
    -- | Create a transform matrix using up direction, camera position and a point to look at.
    --   Just the same as GluLookAt.
    lookAt      :: Vector t 3 -- ^ The up direction, not necessary unit length or perpendicular to the view vector
                -> Vector t 3 -- ^ The viewers position
                -> Vector t 3 -- ^ The point to look at
                -> Matrix t 4 4
    -- | A perspective symmetric projection matrix. Right-handed coordinate system. (@x@ - right, @y@ - top)
    --   http://en.wikibooks.org/wiki/GLSL_Programming/Vertex_Transformations
    perspective :: t -- ^ Near plane clipping distance (always positive)
                -> t -- ^ Far plane clipping distance (always positive)
                -> t -- ^ Field of view of the y axis, in radians
                -> t -- ^ Aspect ratio, i.e. screen's width\/height
                -> Matrix t 4 4
    -- | An orthogonal symmetric projection matrix. Right-handed coordinate system. (@x@ - right, @y@ - top)
    --   http://en.wikibooks.org/wiki/GLSL_Programming/Vertex_Transformations
    orthogonal  :: t -- ^ Near plane clipping distance
                -> t -- ^ Far plane clipping distance
                -> t -- ^ width
                -> t -- ^ height
                -> Matrix t 4 4
    -- | Add one more dimension and set it to 1.
    toHomPoint  :: Vector t 3 -> Vector t 4
    -- | Add one more dimension and set it to 0.
    toHomVector :: Vector t 3 -> Vector t 4
    -- | Transform a homogenous vector or point into a normal 3D vector.
    --   If the last coordinate is not zero, divide the rest by it.
    fromHom     :: Vector t 4 -> Vector t 3

-- Type abbreviations

type Mat22f = Matrix Float 2 2
type Mat32f = Matrix Float 3 2
type Mat42f = Matrix Float 4 2
type Mat23f = Matrix Float 2 3
type Mat33f = Matrix Float 3 3
type Mat43f = Matrix Float 4 3
type Mat24f = Matrix Float 2 4
type Mat34f = Matrix Float 3 4
type Mat44f = Matrix Float 4 4

type Mat22d = Matrix Double 2 2
type Mat32d = Matrix Double 3 2
type Mat42d = Matrix Double 4 2
type Mat23d = Matrix Double 2 3
type Mat33d = Matrix Double 3 3
type Mat43d = Matrix Double 4 3
type Mat24d = Matrix Double 2 4
type Mat34d = Matrix Double 3 4
type Mat44d = Matrix Double 4 4

-- | Compose a 2x2D matrix
mat22 :: PrimBytes (t :: Type)
      => Vector t 2 -> Vector t 2 -> Matrix t 2 2
mat22 = packDF @_ @2 @'[2]
{-# INLINE mat22 #-}

-- | Compose a 3x3D matrix
mat33 :: PrimBytes (t :: Type)
      => Vector t 3 -> Vector t 3 -> Vector t 3 -> Matrix t 3 3
mat33 = packDF @_ @3 @'[3]
{-# INLINE mat33 #-}

-- | Compose a 4x4D matrix
mat44 :: PrimBytes (t :: Type)
      => Vector t 4
      -> Vector t 4
      -> Vector t 4
      -> Vector t 4
      -> Matrix t 4 4
mat44 = packDF @_ @4 @'[4]
{-# INLINE mat44 #-}

instance ( KnownDim n, KnownDim m
         , PrimArray t (Matrix t n m)
         , PrimArray t (Matrix t m n)
         ) => MatrixTranspose t (n :: Nat) (m :: Nat) where
    transpose df = case uniqueOrCumulDims df of
      Left a -> broadcast a
      Right _
         | wm <- dimVal' @m
         , wn <- dimVal' @n
         , m <- case wm of W# w -> word2Int# w
         , n <- case wn of W# w -> word2Int# w
         -> let f ( I# i,  I# j )
                  | isTrue# (i ==# n) = f ( 0, I# (j +# 1#) )
                  | otherwise         = (# ( I# (i +# 1#), I# j )
                                         , ix# (i *# m +# j) df
                                         #)
            in case gen# (CumulDims [wm*wn, wn, 1]) f (0,0) of (# _, r #) -> r

instance MatrixTranspose (t :: Type) (xn :: XNat) (xm :: XNat) where
    transpose (XFrame (df :: DataFrame t ns))
      | ((D :: Dim n) :* (D :: Dim m) :* U) <- dims @ns
      , Dict <- inferPrimElem df
      = XFrame (transpose df :: Matrix t m n)
    transpose _ = error "MatrixTranspose/transpose: impossible argument"

instance (KnownDim n, PrimArray t (Matrix t n n), Num t)
      => SquareMatrix t n where
    eye
      | n <- dimVal' @n
      = let f 0 = (# n    , 1 #)
            f k = (# k - 1, 0 #)
        in case gen# (CumulDims [n*n, n, 1]) f 0 of
            (# _, r #) -> r
    diag se
      | n <- dimVal' @n
      , e <- unScalar se
      = let f 0 = (# n    , e #)
            f k = (# k - 1, 0 #)
        in case gen# (CumulDims [n*n, n, 1]) f 0 of
            (# _, r #) -> r
    trace df
      | I# n <- fromIntegral $ dimVal' @n
      , n1 <- n +# 1#
      = let f 0# = ix# 0# df
            f k  = ix# k  df + f (k -# n1)
        in scalar $ f (n *# n -# 1#)

max_stars_repo_name: cjay/easytensor

max_issues_repo_name: cjay/easytensor

max_forks_repo_name: cjay/easytensor


text:
import Debug.Trace
import Data.Complex

traceIdShow :: Show a => a -> a
traceIdShow x = trace (show x) x

tracef :: Show b => (a->b) -> a -> a
tracef f x = trace (show $ f x) x

type Coord = Complex Float

data Action = Turn Coord
            | Forward Float
            | Move Coord
  deriving Show

data Ship = Ship Coord Coord
  deriving Show

act :: Ship -> Action -> Ship
act (Ship orientation position) (Turn a)
  = traceIdShow $
    trace ("Turning "++ show a) $
    Ship (orientation * a) position
act (Ship orientation position) (Forward d)
  = traceIdShow $
    trace ("Forward "++ show d) $
    Ship orientation (position + orientation * (d:+0))
act (Ship orientation position) (Move p)
  = traceIdShow $
    trace ("Moving "++ show p) $
    Ship orientation (position + p)

-- part 2 changes the orientation to be a waypoint. Because we
-- represent the orientation as a complex number, the changes we make
-- are minimal: we apply moves to the waypoint rather than the ship's
-- position. Also, when calling act, change the start position of the
-- waypoint.
act2 :: Ship -> Action -> Ship
act2 (Ship waypoint position) (Turn a)
  = traceIdShow $
    trace ("Turning "++ show a) $
    Ship (waypoint * a) position
act2 (Ship waypoint position) (Forward d)
  = traceIdShow $
    trace ("Forward "++ show d) $
    Ship waypoint (position + waypoint * (d:+0))
act2 (Ship waypoint position) (Move p)
  = traceIdShow $
    trace ("Moving "++ show p) $
    Ship (waypoint + p) position

main = interact $ (++ "\n") .  show . length
  . tracef ((\d-> "Part 2: " ++ show d)
            . manhatten
            . foldl act2 (Ship (10:+1) (0:+0)))
  . tracef ((\d-> "Part 1: " ++ show d)
            . manhatten
            . foldl act (Ship (1:+0) (0:+0)))
  . map toAction . lines

toAction (a:vs) | a == 'N' = Move $ 0 :+ v
                | a == 'S' = Move $ 0 :+ (-1 * v)
                | a == 'E' = Move $ v :+ 0
                | a == 'W' = Move $ (-1 * v) :+ 0
                | a == 'L' = Turn $  cis $ v / 180 * pi
                | a == 'R' = Turn $  cis $ -1 * v / 180 * pi
                | a == 'F' = Forward $  v
                | otherwise = error $ "unknown action "++ show a
  where v = read vs :: Float

manhatten :: Ship -> Int
manhatten (Ship _ p) = round $ abs(realPart p) + abs(imagPart p)


text:
{-# language TypeFamilies #-}
{-# language ScopedTypeVariables #-}
{-# language RankNTypes #-}

-- | A <https://en.wikipedia.org/wiki/Monus commutative monoid with monus>
-- is a 'Monoid' equipped with a subtraction operator.
module Data.Monoid.Monus
  ( Monus(..)
  , (-)
  ) where

import Prelude hiding ((-))
import Data.Set (Set)
import Data.Complex (Complex(..))
import Data.Monoid (Any(..),All(..),Sum(..), Endo(..))
import Control.Applicative (liftA2)
import Numeric.Natural (Natural)
import Data.Foldable
import Data.Coerce

import qualified Prelude as P
import qualified Data.Set as S

infixl 6 -

-- | A commutative monoid that supports subtraction. The following
-- laws must hold:
--
-- > x <> (y - x) = y <> (x - y)
-- > (x - y) - z = x - (y <> z)
-- > x - x = mempty
-- > mempty - x = mempty
class Monoid a => Monus a where
  monus :: a -> a -> a

-- | An infix synonym for 'subtraction'.
(-) :: Monus a => a -> a -> a
(-) = monus
{-# INLINE (-) #-}

instance Ord a => Monus (Set a) where
  monus = S.difference
  {-# INLINE monus #-}

-- | Unlike the subtraction provided by the 'Num' instance of
-- 'Natural', this subtraction is total.
instance (a ~ Natural) => Monus (Sum a) where
  monus (Sum x) (Sum y) = Sum (if x > y then x P.- y else 0)
  {-# INLINE monus #-}

-- | Defined as @P - Q = P ∧ ¬Q@
instance Monus Any where
  monus (Any x) (Any y) = case x of
    False -> Any False
    True -> Any (not y)
  {-# INLINE monus #-}

-- | Defined as @P - Q = P ∨ ¬Q@
instance Monus All where
  monus (All x) (All y) = case x of
    False -> All (not y)
    True -> All True
  {-# INLINE monus #-}

instance forall a. Monus a => Monus (Endo a) where
  monus = coerce (liftA2 monus :: (a -> a) -> (a -> a) -> (a -> a))
  {-# INLINE monus #-}

instance Monus () where
  monus _ _ = ()
  {-# INLINE monus #-}

instance (Monus a, Monus b) => Monus (a,b) where
  monus (a1,b1) (a2,b2) = (monus a1 a2,monus b1 b2)
  {-# INLINE monus #-}

instance (Monus a, Monus b, Monus c) => Monus (a,b,c) where
  monus (a1,b1,c1) (a2,b2,c2) = (monus a1 a2,monus b1 b2,monus c1 c2)
  {-# INLINE monus #-}

instance (Monus a, Monus b, Monus c,Monus d) => Monus (a,b,c,d) where
  monus (a1,b1,c1,d1) (a2,b2,c2,d2) =
    (monus a1 a2,monus b1 b2,monus c1 c2,monus d1 d2)
  {-# INLINE monus #-}

instance (Monus a, Monus b, Monus c,Monus d,Monus e) => Monus (a,b,c,d,e) where
  monus (a1,b1,c1,d1,e1) (a2,b2,c2,d2,e2) =
    (monus a1 a2,monus b1 b2,monus c1 c2,monus d1 d2,monus e1 e2)
  {-# INLINE monus #-}

instance Monus b => Monus (a -> b) where
  monus = liftA2 monus
  {-# INLINE monus #-}

instance Monus a => Monus (IO a) where
  monus = liftA2 monus
  {-# INLINE monus #-}

instance Monus a => Monus [a] where
  monus [] _ = []
  monus x [] = x
  monus (x:xs) (y:ys) = monus x y : monus xs ys

instance Monus a => Monus (Maybe a) where
  monus = liftA2 monus
  {-# INLINE monus #-}

max_stars_repo_name: andrewthad/monus

max_issues_repo_name: andrewthad/monus

max_forks_repo_name: andrewthad/monus