text:
{-# LANGUAGE DataKinds, NamedFieldPuns, GADTs, KindSignatures, TypeFamilies, TypeOperators, TypeApplications, ScopedTypeVariables, MultiParamTypeClasses, OverloadedLabels, InstanceSigs, FlexibleContexts, AllowAmbiguousTypes, FlexibleInstances, RankNTypes, UndecidableInstances, ConstraintKinds, OverlappingInstances, TypeFamilyDependencies #-}

module TraceTypes where

import           Data.Row.Records        hiding ( Disjoint
                                                , sequence
                                                , map
                                                , zip
                                                )
import           Data.Row.Internal              ( Subset
                                                , type (.\\)
                                                , type (.\/)
                                                )
import           Data.Default
import           Numeric.Log                   as Log
import qualified Data.Vector.Sized             as V
import           Data.Vector.Sized              ( Vector )
import           GHC.OverloadedLabels
import           Data.Finite
import           GHC.TypeLits
import           Control.Monad
import           Control.Monad.Bayes.Class     as Bayes
import Statistics.Distribution (logDensity, complCumulative)
import Statistics.Distribution.Gamma
import Statistics.Distribution.Beta
import Statistics.Distribution.Normal

-- The "Constraint" type is intended only
-- for internal use, and is not user-facing.
-- A Constraint of a certain record type is
-- a record where each field type `a` has been
-- replaced with the type `Maybe a`; a value of
-- `Nothing` indicates an unconstrained value,
-- whereas `Just x` indicates an observed (i.e.
-- constrained) value.
type Constraint r = Rec (Map Maybe r)

-- (ConcatTo l r s) means that l, r, and s are all
-- valid record types, and that l .+ r = r .+ l = s.
type ConcatTo l r s
  = ( WellBehaved l
    , WellBehaved r
    , WellBehaved s
    , Subset l s
    , Subset r s
    , (l .+ r) ≈ s
    , (r .+ l) ≈ s
    , s .\\ l ≈ r
    , s .\\ r ≈ l
    )

-- Are the label sets of l and r disjoint?
type Disjoint l r
  = ( ConcatTo l r (l .+ r)
    , ConcatTo (Map Maybe l) (Map Maybe r) (Map Maybe (l .+ r))
    , Forall (Map Maybe l) Default
    , Forall (Map Maybe r) Default
    )

-- Three main types: probability distributions (D), unnormalized density-carrying measures (U),
-- and probabilistic programs (P).
data U m a = U
  { usampler :: m a
  , udensity :: a -> Log Double
  }
data D m a = D
  { dsampler :: m a
  , ddensity :: a -> Log Double
  }
data P m t a = P
  { traced    :: D m (Rec t)
  , retval    :: Rec t -> a
  , constrain :: Constraint t -> U m (Rec t)
  }

class DensityCarrying f m where
  sampler :: f m a -> m a
  density :: f m a -> a -> Log Double

instance DensityCarrying U m where
  sampler = usampler
  density = udensity

instance DensityCarrying D m where
  sampler = dsampler
  density = ddensity

dist sampler density = D { dsampler = sampler, ddensity = density }
udist sampler density = U { usampler = sampler, udensity = density }

-- DEFAULT INFERENCE
observe
  :: forall m a t s
   . (Disjoint t s, Monad m)
  => P m (t .+ s) a
  -> Rec t
  -> U m (Rec s)
observe p obs =
  let
    -- Turn the observation record `obs` of type `t` into a Constraint of
    -- type `t .+ s`, by filling in the `s` piece of the record with `Nothing`.
      constraints :: Constraint (t .+ s)
      constraints = map' Just obs .+ (default' @Default def :: Constraint s)

      -- Obtain a constrained distribution over traces.
      u :: U m (Rec (t .+ s))
      u = constrain p constraints
  in
    -- Modify the distribution u so that it is over only the unobserved values
    -- (Rec s), not full-but-constrained traces (Rec (t .+ s)).
      udist (fmap restrict $ sampler u) (density u . ((.+) obs))

traceSampler = sampler . traced
traceDensity = density . traced

-- UNIT INTERVAL TYPE

newtype Unit = Sigmoid {logit :: Double}

fromUnit :: Unit -> Double
fromUnit (Sigmoid x) = 1.0 / (1.0 + exp (-x))

mkUnit :: Double -> Unit
mkUnit u = if 0 < u && u < 1
  then Sigmoid (log (u / (1.0 - u)))
  else error "Number must be between 0 and 1."

instance Show Unit where
  show x = show (fromUnit x)

------------------------------
-- STPL LANGUAGE CONSTRUCTS --
------------------------------

-- RETURN

pReturn :: Monad m => a -> P m Empty a
pReturn x =
  let sampler = return Data.Row.Records.empty
      density _ = 1
      traced = dist sampler density
      constrain _ = udist sampler density
      retval _ = x
  in  P { traced, retval, constrain }

-- BIND

pBind :: (Disjoint t s, Monad m) => P m t a -> (a -> P m s b) -> P m (t .+ s) b
pBind f g =
  let
    bindSampler = do
      fTrace <- sampler $ traced f
      let x = retval f fTrace
      gTrace <- sampler $ traced (g x)
      return (fTrace .+ gTrace)

    bindDensity tr =
      let fDensity = density (traced f) (restrict tr)
          x        = retval f (restrict tr)
          gDensity = density (traced $ g x) (restrict tr)
      in  fDensity * gDensity

    bindTraced = dist bindSampler bindDensity

    bindRet tr = let x = retval f (restrict tr) in retval (g x) (restrict tr)

    constrainedSampler c = do
      fTrace <- sampler $ constrain f (restrict c)
      let x = retval f fTrace
      gTrace <- sampler $ constrain (g x) (restrict c)
      return (fTrace .+ gTrace)

    constrainedDensity c tr =
      let fDensity = density (constrain f (restrict c)) (restrict tr)
          x        = retval f (restrict tr)
          gDensity = density (constrain (g x) (restrict c)) (restrict tr)
      in  fDensity * gDensity

    bindConstrain c = udist (constrainedSampler c) (constrainedDensity c)
  in
    P { traced = bindTraced, retval = bindRet, constrain = bindConstrain }

-- SAMPLE

singleLabelConstrain
  :: (KnownSymbol l, MonadInfer m)
  => Label l
  -> D m (Rec (l .== a))
  -> Constraint (l .== a)
  -> U m (Rec (l .== a))
singleLabelConstrain lbl traced c = case c .! lbl of
  Just x ->
    let t                  = lbl .== x
        constrainedSampler = do
          score $ density traced t
          return t
        constrainedDensity _ = density traced t
    in  udist constrainedSampler constrainedDensity

  Nothing -> udist (sampler traced) (density traced)

sample :: (KnownSymbol l, MonadInfer m) => Label l -> D m a -> P m (l .== a) a
sample lbl d =
  let sampleTraced = dist
        (do
          x <- sampler d
          return (lbl .== x)
        )
        (\t -> density d (t .! lbl))

      sampleRet t = t .! lbl

      sampleConstrain = singleLabelConstrain lbl sampleTraced
  in  P { traced    = sampleTraced
        , retval    = sampleRet
        , constrain = sampleConstrain
        }

-- WITH PROBABILITY

withProbability
  :: (KnownSymbol l, MonadInfer m)
  => Label l
  -> Unit
  -> P m t a
  -> P m s a
  -> P m (l .== Either (Rec t) (Rec s)) a
withProbability lbl logit f g =
  let p         = fromUnit logit

      wpSampler = do
        r <- random
        if r < p
          then do
            t <- sampler $ traced f
            return (lbl .== Left t)
          else do
            t <- sampler $ traced g
            return (lbl .== Right t)

      wpDensity tr = case tr .! lbl of
        Left  t -> Exp (log p) * density (traced f) t
        Right t -> Exp (log (1.0 - p)) * density (traced g) t

      wpTraced = dist wpSampler wpDensity

      wpRet t = case t .! lbl of
        Left  t -> retval f t
        Right t -> retval g t

      wpConstrain = singleLabelConstrain lbl wpTraced
  in  P { traced = wpTraced, retval = wpRet, constrain = wpConstrain }

-- FOR EACH

forEach
  :: (KnownSymbol l, KnownNat n, MonadInfer m)
  => Label l
  -> Vector n a
  -> (a -> P m t b)
  -> P m (l .== Vector n (Rec t)) (Vector n b)
forEach lbl xs body =
  let
    forEachSampler = do
      traces <- V.mapM (sampler . traced . body) xs
      return (lbl .== traces)

    forEachDensity t = V.product
      (V.map (\(x, tr) -> density (traced $ body x) tr) (V.zip xs (t .! lbl)))

    forEachTraced = dist forEachSampler forEachDensity

    forEachRet t = V.map (\(x, tr) -> retval (body x) tr) (V.zip xs (t .! lbl))

    forEachConstrain = singleLabelConstrain lbl forEachTraced
  in
    P { traced    = forEachTraced
      , retval    = forEachRet
      , constrain = forEachConstrain
      }

-- FOR ... IN RANDOM RANGE ...

forRandomRange
  :: (KnownSymbol l, MonadInfer m)
  => Label l
  -> D m Int
  -> (Int -> P m t a)
  -> P m (l .==  [Rec t]) [a]
forRandomRange lbl d body =
  let
    -- Sample a trace of the for loop.
    forSampler = do
      n      <- sampler d
      traces <- sequence $ map (sampler . traced . body) [0 .. n - 1]
      return (lbl .== traces)

    -- Compute the density of a trace of the for loop.
    forDensity t =
      density d (length $ t .! lbl)
        * product
            (map (\(tr, i) -> density (traced $ body i) tr)
                 (zip (t .! lbl) [0 ..])
            )

    -- The for loop's density-carrying distribution over traces
    forTraced = dist forSampler forDensity

    -- Compute the return value of the for loop, given a trace of an execution.
    -- The return value is just a list of the return values from each iteration.
    forRet t = map (\(tr, i) -> retval (body i) tr) (zip (t .! lbl) [0 ..])

    -- Constrain the random choices made inside the for loop.
    forConstrain = singleLabelConstrain lbl forTraced
  in
    P { traced = forTraced, retval = forRet, constrain = forConstrain }

{- WHILE:

    A stochastic while loop.
    Takes as input:
    - lbl    : The label at which to trace the while loop's random choices.
    - init   : The initial state of the while loop.
    - pi     : A function computing a probability of continuing, based on the current state.
    - pi_max : An upper bound beyond which the probability of continuing is truncated.
    - body   : The body of the while loop, which simulates a new state based on the current state.

    The resulting probabilistic program has type P m (l .== [Rec t]) a:
    - The trace contains, at the specified label, a list of subtraces, one for each iteration.
    - The return value is the final state at the finish of the while loop.

 -}

while
  :: (KnownSymbol l, MonadInfer m)
  => Label l
  -> a
  -> (a -> Unit)
  -> Unit
  -> (a -> P m t a)
  -> P m (l .== [Rec t]) a
while lbl init pi pi_max body =
  let
    -- Sample a trace of the while loop.
    whileSampler state = do
      shouldContinue <- bernoulli $ min (fromUnit (pi state)) (fromUnit pi_max)
      if shouldContinue
        then do
          nextTrace <- sampler (traced $ body state)
          let nextState = retval (body state) nextTrace
          restOfTraces <- whileSampler nextState
          return (lbl .== (nextTrace : (restOfTraces .! lbl)))
        else return (lbl .== [])

    -- Given traces from each iteration, and an initial state,
    -- compute the final state.
    retvalFromTraces init ts = case ts of
      []     -> init
      t : ts -> retvalFromTraces (retval (body init) t) ts

    -- Given a trace of the entire while loop, compute the final
    -- state.
    whileRet t = retvalFromTraces init (t .! lbl)

    -- Given an initial state and traces from each iteration,
    -- compute the density.
    densityFromTraces init ts = case ts of
        [] -> Exp (log (1.0 - (fromUnit (pi init))))
        t : ts ->
          Exp (log (fromUnit (pi init)))
            * (density (traced $ body init) t)
            * densityFromTraces (retval (body init) t) ts

    -- Given a trace of the entire while loop, compute the density.
    whileDensity t = densityFromTraces init (t .! lbl)

    -- The while loop's density-carrying distribution over traces.
    whileTraced    = dist (whileSampler init) whileDensity

    -- Constrain the choices made inside the loop.
    whileConstrain = singleLabelConstrain lbl whileTraced

  in  P { traced = whileTraced, retval = whileRet, constrain = whileConstrain }

-------------------
-- DISTRIBUTIONS --
-------------------

nrm :: MonadInfer m => Double -> Log Double -> D m Double
nrm mu (Exp sigln) = dist (normal mu (exp sigln)) (normalPdf mu (exp sigln))

brn :: MonadSample m => Unit -> D m Bool
brn logit =
  let p = fromUnit logit
  in  dist (bernoulli p) (\b -> Exp . log $ if b then p else (1.0 - p))

gmma :: MonadSample m => Log Double -> Log Double -> D m (Log Double)
gmma shape scale =
  let toDouble (Exp l) = exp l
      toLogDouble d = Exp (log d)
  in
    dist (fmap toLogDouble (gamma (toDouble shape) (toDouble scale)))
         (Exp . (logDensity $ Statistics.Distribution.Gamma.gammaDistr
                                (toDouble shape) (toDouble scale))
              . toDouble)

geom :: MonadSample m => Unit -> D m Int
geom p' =
  let p = fromUnit p'
  in dist (geometric p)
          (\n -> Exp (fromIntegral n * (log p) + (log (1.0 - p))))

pois :: MonadSample m => Int -> D m Int
pois = undefined

bta :: MonadSample m => Log Double -> Log Double -> D m Unit
bta (Exp lna) (Exp lnb) =
  let
    a = exp lna
    b = exp lnb
  in
    dist (fmap mkUnit (beta a b))
         (Exp . (logDensity $ Statistics.Distribution.Beta.betaDistr a b)
              . fromUnit)

unif :: MonadSample m => D m Unit
unif = dist (fmap mkUnit Bayes.random) (\x -> 1.0 :: Log Double)

cat :: (KnownNat n, MonadSample m) => Vector n Unit -> D m (Finite n)
cat probs = undefined

lognormal :: forall m. MonadInfer m => Log Double -> Log Double -> D m (Log Double)
lognormal (Exp logmu) sigma =
  dist (fmap Exp $ sampler $ nrm logmu sigma)
       (\(Exp logx) -> density (nrm logmu sigma :: D m Double) logx / (Exp logx))

truncnrm :: MonadInfer m => Log Double -> Log Double -> D m (Log Double)
truncnrm (Exp logmu) (Exp logsigma) =
  let mu = exp logmu
      sigma = exp logsigma
      unbounded = nrm mu (Exp logsigma)

      truncNormSampler = do
        x <- sampler unbounded
        if x < 0 then truncNormSampler else return (Exp $ log x)

      truncNormPdf (Exp logy) =
        let y = exp logy
            z = complCumulative (normalDistr mu sigma) 0
        in (density unbounded y / (Exp (log z)))
  in dist truncNormSampler truncNormPdf

----------------------------
-- PROGRAMMABLE INFERENCE --
----------------------------

-- IMPORTANCE SAMPLING

importance
  :: (Disjoint t s, MonadInfer m)
  => P m (t .+ s) a
  -> Rec t
  -> P m s b
  -> m (Rec s)
importance p t q =
  let target = observe p t
  in  do
        tr <- sampler $ traced q
        score ((density target tr) / (density (traced q) tr))
        return tr

-- SMC

unroll
  :: (MonadInfer m, Disjoint t s)
  => P m i a
  -> (a -> P m (t .+ s) a)
  -> [Rec t]
  -> m (Rec i, [Rec s])
unroll init next steps = do
  tInit <- sampler $ traced init
  let processStep soFar t = do {
    (s, results) <- soFar;
    tNext <- sampler $ observe (next s) t;
    return (retval (next s) (tNext .+ t), tNext : results)
  }
  (s, l) <- foldl processStep (return (retval init tInit, [])) steps
  return (tInit, reverse l)

particleFilter
  :: (MonadInfer m, Disjoint t s,
      Disjoint Empty i, (Empty .+ i) ≈ i)
  => P m i a
  -> P m i b
  -> (a -> P m (t .+ s) a)
  -> [(Rec t, a -> P m s c)]
  -> m (Rec i, [Rec s])
particleFilter init qInit next steps = do
  tInit <- importance init Data.Row.Records.empty qInit
  let processStep soFar (t, q) = do {
    (s, results) <- soFar;
    tNext <- importance (next s) t (q s);
    return (retval (next s) (tNext .+ t), tNext : results)
  }
  (s, l) <- foldl processStep (return (retval init tInit, [])) steps
  return (tInit, reverse l)

-- MCMC

type K m t (s :: Row *) = U m (Rec t) -> Rec t -> m (Rec t)

mh
  :: forall s t m a
   . (Disjoint t s, MonadInfer m)
  => (Rec (t .+ s) -> P m t a)
  -> K m (t .+ s) t
mh q p old =
  let proposal = do
        new <- sampler (traced $ q old)
        return $ new .+ (restrict old :: Rec s)
      rho (x, y) =
          ((density p y) * (density (traced $ q y) (restrict x)))
            / ((density p x) * (density (traced $ q x) (restrict y)))
  in  do
        proposed <- proposal
        r        <- Bayes.random
        if Exp (log r) < (min 1 (rho (old, proposed)))
          then return proposed
          else return old

seqK :: Monad m => K m t s -> K m t r -> K m t (s .\/ r)
seqK k1 k2 p old = do
  t' <- k1 p old
  k2 p t'

mixK :: MonadSample m => Double -> K m t s -> K m t r -> K m t (s .\/ r)
mixK p k1 k2 model old = do
  r <- Bayes.random
  if r < p then k1 model old else k2 model old

repeatK :: Monad m => Int -> K m t s -> K m t s
repeatK n k p old = if n == 0
  then return old
  else do
    t <- k p old
    repeatK (n - 1) k p t

ifK
  :: (Disjoint t s, Monad m)
  => (Rec t -> Bool)
  -> K m (t .+ s) s
  -> K m (t .+ s) s
ifK b k p old = if (b (restrict old)) then k p old else return old

-- VARIATIONAL (STUBBED)

svi
  :: (Disjoint t s, Monad m)
  => P m (t .+ s) a
  -> Rec t
  -> (theta -> P m s b)
  -> theta
  -> m theta
svi p t q theta = undefined
-- do tr <- traceSampler $ q theta
--    let sco  = density (observe p t) tr
--    let gradient = grad (\th -> traceDensity (q th) tr) theta
--    return theta + stepSize * sco * gradient

trainAmortized
  :: (Disjoint t s, Monad m)
  => P m (t .+ s) a
  -> (theta -> Rec t -> P m s b)
  -> theta
  -> m theta
trainAmortized p q theta = undefined
-- do tr <- traceSampler p
--    let gradient = grad (\th -> traceDensity (q th (restrict tr)) (restrict tr)) theta
--    return theta + step * gradient

ifThenElse c t f = if c then t else f

max_stars_repo_name: probcomp/haskell-trace-types

max_issues_repo_name: probcomp/haskell-trace-types

max_forks_repo_name: probcomp/haskell-trace-types


text:
module Main where

import Test.Hspec
import Test.QuickCheck
import System.Console.Ansigraph
import System.Console.Ansigraph.Examples (wave)
import Data.Complex
import Data.Char (isSpace)
-- for GHC <= 7.8
import Control.Applicative

---- For horizontal graphing tests ----

bars = " ▁▂▃▄▅▆▇█"

eighths :: [Double]
eighths = (/8) <$> [0..8]

{- These two vectors have a maximum element of 8, to maintain a particular normalization,
   while displacing other entries by just under one half. Both vectors should display
   as █▇▆▅▄▃▂▁. Same idea for m1, m_hi, m_low below. -}

v_hi :: [Double]
v_hi = 8 : map (+0.48) (reverse [0..7])

v_low :: [Double]
v_low = 8 : map (\x -> x - 0.48) (reverse [0..7])

{- These functions are for testing the property that:
   simpleRender v == complement <$> simpleRenderR v-}

-- the 'complement' of a particular bar character is the one representing the negation of the
-- corresponding numeric value (when displayed with reversed foreground/background colors).
complement :: Char -> Char
complement c = match c bars (reverse bars)

match :: Eq a => a -> [a] -> [a] -> a
match x (y:ys) (z:zs) = if x == y then z else match x ys zs

rwave = realPart <$> wave

barInvert :: (String,String) -> (String,String)
barInvert (x,y)  = (complement <$> y, complement <$> x)

barInvertC :: (String,String,String,String) -> (String,String,String,String)
barInvertC (x1,y1,x2,y2) = (complement <$> y1,
                            complement <$> x1,
                            complement <$> y2,
                            complement <$> x2)

---- For matrix graphing tests ----

fourths :: [Double]
fourths = (/4) <$> [1..4]

mmap :: (a -> b) -> [[a]] -> [[b]]
mmap = map . map

m1 :: [[Double]]
m1 = [[1,2,3],[3,4,4]]

m_low :: [[Double]]
m_low = [[0.55,1.55,2.55], [2.55,3.55,4]]

m_hi :: [[Double]]
m_hi = [[1.45,2.45,3.45], [3.45,4,4]]

---- The spec ----

main :: IO ()
main = hspec $ do

  describe "renderPV" $ do

    it "maps exact multiples of 1/8 to the correct characters" $
      renderPV eighths `shouldBe` bars

    it "maps scaled multiples of 1/8 to the correct characters" $
      renderPV ((* 1337) <$> eighths) `shouldBe` bars

    it "maps values near tops of rounding regions to the correct characters" $
      renderPV v_hi `shouldBe` reverse bars

    it "maps values near bottoms of rounding regions to the correct characters" $
      renderPV v_low `shouldBe` reverse bars

    it "maps zero vectors to whitespace" $
      renderPV [0,0,0] `shouldBe` "   "

    it "maps non-positive vectors to whitespace" $
      renderPV [-1,-23,-0.02] `shouldBe` "   "

    it "maps null vectors to null strings" $
      renderPV [] `shouldBe` ""

  describe "renderRV" $ do

    it "inverts 'rwave' consistently" $
      let (p,n) = renderRV rwave
      in  renderRV (negate <$> rwave) `shouldBe` (complement <$> n, complement <$> p)

    it "inverts arbitrary vectors consistently" $
      property $ \xs -> renderRV xs == barInvert (renderRV $ negate <$> xs)

    it "maps arbitrary-dimensional zero vectors to whitespace, for positive components" $
      property $ \n -> let (p,m) = renderRV $ replicate n 0
                       in  all isSpace p

    it "maps arbitrary-dimensional zero vectors to solid bocks █, for negative components" $
      property $ \n -> let (p,m) = renderRV $ replicate n 0
                       in  all (== '█') m

  describe "renderCV" $ do

    it "inverts 'wave' consistently" $ do
      let (rp,rn,ip,im) = renderCV wave
      renderCV (negate <$> wave) `shouldBe` (complement <$> rn,
                                             complement <$> rp,
                                             complement <$> im,
                                             complement <$> ip)

    it "inverts arbitrary complex vectors consistently" $
      property $ \zs -> renderCV zs == barInvertC (renderCV $ negate <$> zs)

    it "maps arbitrary dimensional zero vectors to whitespace, for positive components" $
      property $ \n -> let (r1,_,i1,_) = renderCV $ replicate n (0.0 :+ 0.0)
                       in  all isSpace $ r1 ++ i1

    it "maps arbitrary dimensional zero vectors to solid bocks █, for negative components" $
      property $ \n -> let (_,r2,_,i2) = renderCV $ replicate n (0.0 :+ 0.0)
                       in  all (== '█') $ r2 ++ i2

    it "transforms appropriately under exchange of real and complex components" $
      property $ \v -> let v'  = (\(x :+ y) -> (y :+ x)) <$> v
                           (x1,y1,x2,y2)  = renderCV v
                           (x3,y3,x4,y4)  = renderCV v'
                       in  and [x1 == x4,
                                y1 == y4,
                                x2 == x3,
                                y2 == y3]

  describe "matShow" $ do

    it "maps exact multiples of 1/4 to the correct characters" $
      matShow [[1/4,2/4],[3/4,4/4]] `shouldBe` ["░▒", "▓█"]

    it "maps scaled multiples of 1/4 to the correct characters" $ do
      matShow (mmap (* 3097) [[1,2],[3,4]]) `shouldBe` ["░▒", "▓█"]
      matShow m1 `shouldBe` ["░▒▓", "▓██"]

    it "maps values near top of rounding regions to the correct characters" $
      matShow m_hi `shouldBe` ["░▒▓", "▓██"]

    it "maps values near bottom of rounding regions to the correct characters" $
      matShow m_low `shouldBe` ["░▒▓", "▓██"]

ext: hs

max_stars_repo_name: fieldstrength/ansigraph

max_issues_repo_name: fieldstrength/ansigraph

max_forks_repo_name: fieldstrength/ansigraph


text:
{-# LANGUAGE DataKinds #-}
module Grenade.Utils.LinearAlgebra
    ( sumV
    , sumM
    , squareV
    , squareM
    ) where

import           GHC.TypeLits
import           Grenade.Types
import           Numeric.LinearAlgebra.Static

sumV :: (KnownNat n) => R n -> RealNum
sumV v = v <.> 1

sumM :: (KnownNat m, KnownNat n) => L m n -> RealNum
sumM m = (m #> 1) <.> 1

squareV :: (KnownNat n) => R n -> R n
squareV v = dvmap (^ (2 :: Int)) v

squareM :: (KnownNat m, KnownNat n) => L m n ->  L m n
squareM m = dmmap (^ (2 :: Int)) m

-- v :: R 10
-- v = vector [1..10]

-- m :: L 3 3
-- m = matrix [-3..5]

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


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

module Language.Grappa.Test.Inference.GradientDescent where

import           Control.Monad ( forM, forM_, replicateM )
import qualified Numeric.LinearAlgebra as La
import qualified Statistics.Distribution as S
import qualified Statistics.Distribution.StudentT as S
import qualified System.Random.MWC as MWC

import           Test.Hspec
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.Rand.MWCRandM
import           Language.Grappa.Inference.GradientDescent

tests :: Spec
tests = do
  forM_ [1..8] $ \ x ->
    it ("should pass this QC test with k=" ++ show x) $
      prop_evalGradLogMVNormalDensity_correct x

-- | 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] -> [Double]
    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 })

-- | For example:
--
-- > runGdTest (gd2 0.001) (gdTest3 4 [1,42])
runGdTest :: Gd Rand -> GdTest Rand -> IO (Double, [Double])
runGdTest gd GdTest{..} = do
  -- Use 'MWC.create' instead if you want repeatable/predictable
  -- results.
  g <- MWC.createSystemRandom
  lambda <- runRand g $ gd _gdParams
  return (_estimateError lambda, lambda)

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

data GdTest m = GdTest
  { _gdParams      :: GdParams m
  , _estimateError :: [Double] -> Double
  }

-- | Test for objective function @f x = (x - 42)**2@, which is
-- minimized at @x = 42@.
gdTest1 :: SampleableIn m Normal => GdTest m
gdTest1 = GdTest
  { _gdParams = GdParams
    { _projectLambda = const id
    , _sampleLambda = defaultSampleLambda 1
    , _evalObj = \[x] -> return $ (x - center!!0)**2
    , _evalGradObj = \[x] -> return [2 * (x - center!!0)]
    }
  , _estimateError = \lambda -> ell1 $ lambda `subV` center
  }
  where
    center = [42]

-- | Test for objective function @f [x,y] = ([x,y] - [42,10^6])**2@,
-- which is minimized at @[42,10^6]@.
gdTest2 :: SampleableIn m Normal => GdTest m
gdTest2 = GdTest
  { _gdParams = GdParams
    { _projectLambda = const id
    , _sampleLambda = defaultSampleLambda 2
    , _evalObj = \[x,y] -> return $ sqrt ( (x - center!!0)**2
                                         + (y - center!!1)**2 )
    , _evalGradObj = \[x,y] -> return [ 2 * (x - center!!0)
                                      , 2 * (y - center!!1) ]
    }
  , _estimateError = \lambda -> ell1 $ lambda `subV` center
  }
  where
    center = [ 42, 1000000 ]

-- | Test for objective function @f x = (x - center)**k@ where
-- @center@ is an arbitrary vector. The dimension of @center@
-- determines the dimension of @x@. The objective is minimized at
-- @center@ for even powers @k@, and has no minimum for odd powers
-- @k@. You're on your own for fractional powers.
gdTest3 :: SampleableIn m Normal => Double -> [Double] -> GdTest m
gdTest3 k center = GdTest
  { _gdParams = GdParams
    { _projectLambda = const id
    , _sampleLambda = defaultSampleLambda (length center)
    , _evalObj = \lambda -> return $ (ell2 $ v lambda)**k
    , _evalGradObj = \lambda -> return
        [ k * (ell2 (v lambda))**(k-2) * vi
        | vi <- v lambda ]
    }
  , _estimateError = \lambda -> ell1 $ lambda `subV` center
  }
  where
    v lambda = lambda `subV` center

-- | Test for objective function
--
-- > f x = (x - center)**k -
-- >       amplitude * cos (2pi * |x - center| / period)
--
-- I.e. like 'gdTest3', but with an added noise factor of @- cos@ to
-- introduce infinitely many local minima. This is a better test of
-- the line search in 'gd2'. I expect that making the period and
-- amplitude large enough, and the center far enough from the origin,
-- would trick 'gd2' into a local minium, but I can't find parameters
-- that actually do this!
--
-- This objective is minimized at @center@ for even powers @k@, like
-- 'gdTest3'.
gdTest4 :: SampleableIn m Normal =>
  Double -> [Double] -> Double -> Double -> GdTest m
gdTest4 k center amplitude period = GdTest
  { _gdParams = GdParams
    { _projectLambda = const id
    , _sampleLambda = defaultSampleLambda (length center)
    , _evalObj = \lambda ->
        return $ (ell2 $ v lambda)**k -
                 amplitude * cos (2 * pi *
                                  (ell2 $ v lambda) / period)
    , _evalGradObj = \lambda -> return
        [ k * (ell2 (v lambda))**(k-2) * vi +
          sin (ell2 (v lambda)) *
            (2 * pi / period) *
            vi / ell2 (v lambda)
        | vi <- v lambda ]
    }
  , _estimateError = \lambda -> ell1 $ lambda `subV` center
  }
  where
    v lambda = lambda `subV` center

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

max_stars_repo_name: kquick/grappa

max_issues_repo_name: kquick/grappa

max_forks_repo_name: kquick/grappa


text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-|
Module      : Grenade.Layers.Internal.Shrink
Description : Functions to shrink a matrix
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Internal.Shrink (
    shrink_2d,
    shrink_2d_rgba
  ) where

import qualified Data.Vector.Storable        as S

import           Foreign                     (mallocForeignPtrArray, withForeignPtr)
import           Foreign.Ptr                 (Ptr)
import           Numeric.LinearAlgebra       (Matrix, flatten)
import qualified Numeric.LinearAlgebra.Devel as U
import           System.IO.Unsafe            (unsafePerformIO)

import           Data.Word

import           Grenade.Types

shrink_2d :: Int            -- ^ input rows
          -> Int            -- ^ input cols
          -> Int            -- ^ target rows
          -> Int            -- ^ target cols
          -> Matrix RealNum -- ^ input matrix
          -> Matrix RealNum -- ^ shrinked matrix
shrink_2d rows cols rows' cols' m
 = let outMatSize      = rows' * cols'
       vec             = flatten m
   in unsafePerformIO $ do
     --print vec
     outPtr           <- mallocForeignPtrArray outMatSize
     let (inPtr, _, _) = S.unsafeToForeignPtr vec

     withForeignPtr inPtr $ \inPtr' ->
        withForeignPtr outPtr $ \outPtr' ->
          shrink_2d_cpu inPtr' rows cols rows' cols' outPtr'

     let matVec = S.unsafeFromForeignPtr0 outPtr outMatSize
     return (U.matrixFromVector U.RowMajor rows' cols' matVec)
--{-# INLINE shrink_2d #-}

shrink_2d_rgba :: Int            -- ^ input rows
               -> Int            -- ^ input cols
               -> Int            -- ^ target rows
               -> Int            -- ^ target cols
               -> S.Vector Word8 -- ^ vector of word8s in the form rgba, like from a HTML canvas
               -> Matrix RealNum -- shrinked matrix
shrink_2d_rgba rows cols rows' cols' vec
  = let outMatSize = rows' * cols'
    in unsafePerformIO $ do
      outPtr           <- mallocForeignPtrArray outMatSize
      let (inPtr, _, _) = S.unsafeToForeignPtr vec

      -- print vec

      withForeignPtr inPtr $ \inPtr' ->
        withForeignPtr outPtr $ \outPtr' ->
          shrink_2d_rgba_cpu inPtr' rows cols rows' cols' outPtr'

      let matVec = S.unsafeFromForeignPtr0 outPtr outMatSize
      return (U.matrixFromVector U.RowMajor rows' cols' matVec)
--{-# INLINE shrink_2d_rgba #-}

foreign import ccall unsafe
    shrink_2d_cpu :: Ptr RealNum -> Int -> Int -> Int -> Int -> Ptr RealNum -> IO ()

foreign import ccall unsafe
    shrink_2d_rgba_cpu :: Ptr Word8 -> Int -> Int -> Int -> Int -> Ptr RealNum -> IO ()

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE Haskell2010
    , TypeFamilies
    , FlexibleContexts
    , Trustworthy
    , StandaloneDeriving
    , DeriveDataTypeable
    , CPP
 #-}
{-# OPTIONS -Wall -fno-warn-name-shadowing #-}
{-# OPTIONS -cpp  -pgmPcpphs  -optP--cpp #-}

-- | Efficient matrix operations in 100% pure Haskell.
--
-- This package uses miscellaneous implementations,
-- depending on the type of its components. Typically unboxed
-- arrays will perform best, while unboxed arrays give you
-- certain features such as 'Rational' or 'Complex' components.
--
-- The following component types are supported by 'Matrix':
--
-- [@Int@] Uses unboxed arrays internally. 'inv' will always
--      return 'Nothing'.
--
-- [@Integer@] Uses boxed arrays internally. 'inv' will always
--      return 'Nothing'.
--
-- [@Double@ and @Float@] Uses unboxed arrays internally.
--      All matrix operations will work as expected.
--      @Matrix Double@ will probably yield the best peformance.
--
-- [@Rational@] Best choice if precision is what you aim for.
--      Uses boxed arrays internally. All matrix operations will
--      work as expected.
--
-- [@Complex@] Experimental. Uses boxed arrays internally.
--      All matrix operations will work as expected, though
--      finding the inverse of a matrix isa tad less numerically
--      stable than with a @Double@ matrix.
module Numeric.Matrix (

    Matrix,

    MatrixElement (..),

    -- * Matrix property and utility functions.
    (<|>),
    (<->),
    scale,

    -- ** Matrix properties
    isUnit,
    isZero,
    isDiagonal,
    isEmpty,
    isSquare,

    -- ** Conversions
    toDoubleMatrix,
    toComplexMatrix,
    toRationalMatrix

) where

import Control.Applicative ((<$>))
import Control.DeepSeq
import Control.Monad
import Control.Monad.ST

import Data.Function (on)
import Data.Ratio
import Data.Complex
import Data.Maybe
import Data.Int
import Data.Word

import qualified Data.List as L
import Data.Array.IArray
import Data.Array.MArray
import Data.Array.Unboxed
import Data.Array.ST
import Data.STRef
import Data.Binary

import qualified Data.Array.Unsafe as U

import Data.Typeable

import Prelude hiding (any, all, read, map)
import qualified Prelude as P

-- | Matrices are represented by a type which fits best the component type.
-- For example a @Matrix Double@ is represented by unboxed arrays,
-- @Matrix Integer@ by boxed arrays.
--
-- Data instances exist for 'Int', 'Float', 'Double', 'Integer', 'Ratio',
-- and 'Complex'. Certain types do have certain disadvantages, like for
-- example you can not compute the inverse matrix of a @Matrix Int@.
--
-- Every matrix (regardless of the component type) has instances for
-- 'Show', 'Read', 'Num', 'Fractional', 'Eq', 'Typeable', 'Binary',
-- and 'NFData'. This means that you can use arithmetic operations like
-- '+', '*', and '/', as well as functions like 'show', 'read', or 'typeOf'.
--
-- [@Show (Matrix e)@]
-- Note that a Show instance for the component type @e@ must exist.
--
-- [@Read (Matrix e)@]
-- You can read a matrix like so:
--
-- > read "1 0\n0 1\n" :: Matrix Double
--
-- [@Num (Matrix e)@]
-- '+', '-', '*', 'negate', 'abs', 'signum', and 'fromInteger'.
--
-- 'signum' will compute the determinant and return the signum
-- of it.
--
-- 'abs' applies @map abs@ on the matrix (that is, it applies
-- @abs@ on every component in the matrix and returns a new
-- matrix without negative components).
--
-- @fromInteger@ yields a 1-x-1-matrix.
--
-- [@Fractional (Matrix e)@]
-- Only available if there exists an instance @Fractional e@
-- (the component type needs to have a @Fractional@ instance, too).
-- Note that while the 'Num' operations are safe, 'recip' and
-- '/' will fail (with an 'error') if the involved matrix is
-- not invertible or not a square matrix.
--
-- [@NFData (Matrix e)@]
-- Matrices have instances for NFData so that you can use a
-- matrix in parallel computations using the @Control.Monad.Par@
-- monad (see the @monad-par@ package for details).
--
-- [@Typeable (Matrix e)@]
-- Allows you to use matrices as 'Dynamic' values.
--
-- [@Binary (Matrix e)@]
-- Serialize and unserialize matrices using the @binary@ package.
-- See @encode@ and @decode@.
data family Matrix e

#if defined(__GLASGOW_HASKELL__) && (__GLASGOW_HASKELL__ >= 707)
deriving instance Typeable Matrix
#else
deriving instance Typeable1 Matrix
#endif

data instance Matrix Int
    = IntMatrix !Int !Int (Array Int (UArray Int Int))

data instance Matrix Float
    = FloatMatrix !Int !Int (Array Int (UArray Int Float))

data instance Matrix Double
    = DoubleMatrix !Int !Int (Array Int (UArray Int Double))

data instance Matrix Integer
    = IntegerMatrix !Int !Int (Array Int (Array Int Integer))

data instance Matrix (Ratio a)
    = RatioMatrix !Int !Int (Array Int (Array Int (Ratio a)))

data instance Matrix (Complex a)
    = ComplexMatrix !Int !Int (Array Int (Array Int (Complex a)))

instance (MatrixElement e, Show e) => Show (Matrix e) where
    show = unlines . P.map showRow . toList
      where
        showRow = unwords . P.map ((' ':) . show)

instance (Read e, MatrixElement e) => Read (Matrix e) where
    readsPrec _ = (\x -> [(x, "")]) . fromList . P.map (P.map P.read . words) . lines

instance (MatrixElement e) => Num (Matrix e) where
    (+) = plus
    (-) = minus
    (*) = times
    abs         = map abs
    signum      = matrix (1,1) . const . signum . det
    fromInteger = matrix (1,1) . const . fromInteger

instance (MatrixElement e, Fractional e) => Fractional (Matrix e) where
    recip        = fromJust . inv
    fromRational = matrix (1,1) . const . fromRational

instance (MatrixElement e) => Eq (Matrix e) where
    m == n
        | dimensions m == dimensions n
            = allWithIndex (\ix e -> m `at` ix == e) n
        | otherwise = False

instance (MatrixElement e) => NFData (Matrix e) where
    rnf matrix = matrix `deepseq` ()

instance (MatrixElement e, Binary e) => Binary (Matrix e) where

    put m = do
        let (rows, cols) = dimensions m

        put rows >> put cols

        forM_ [1..rows] $ \i -> do
            forM_ [1..cols] $ \j -> do
                put (m `at` (i,j))

    get = do
        rows <- get :: Get Int
        cols <- get :: Get Int

        forM [1..rows] (const (forM [1..cols] (const get)))
            >>= return . fromList

(<|>) :: MatrixElement e => Matrix e -> Matrix e -> Matrix e
-- ^ Joins two matrices horizontally.
--
-- > 1 2 3     1 0 0      1 2 3 1 0 0
-- > 3 4 5 <|> 2 1 0  ->  3 4 5 2 1 0
-- > 5 6 7     3 2 1      5 6 7 3 2 1
m1 <|> m2 = let m = numCols m1
                n1 = numRows m1
                n2 = numRows m2
            in matrix (max n1 n2, m + numCols m2)
              $ \(i,j) -> if j > m
                    then (if i > n2 then 0 else m2 `at` (i,j-m))
                    else (if i > n1 then 0 else m1 `at` (i,j))

(<->) :: MatrixElement e => Matrix e -> Matrix e -> Matrix e
-- ^ Joins two matrices vertically.
--
-- > 1 2 3     1 0 0      1 2 3
-- > 3 4 5 <-> 2 1 0  ->  3 4 5
-- > 5 6 7     3 2 1      5 6 7
-- >                      1 0 0
-- >                      2 1 0
-- >                      3 2 1
m1 <-> m2 = let m = numRows m1
                n1 = numCols m1
                n2 = numCols m2
            in matrix (m + numRows m2, max n1 n2)
              $ \(i,j) -> if i > m
                    then (if j > n2 then 0 else m2 `at` (i-m,j))
                    else (if j > n1 then 0 else m1 `at` (i,j))

scale :: MatrixElement e => Matrix e -> e -> Matrix e
-- ^ Scales a matrix by the given factor.
--
-- > scale s == map (*s)
scale m s = map (*s) m

isUnit, isDiagonal, isZero, isEmpty, isSquare :: MatrixElement e => Matrix e -> Bool

-- | Check whether the matrix consists of all zeros.
--
-- > isZero == all (== 0)
isZero = all (== 0)

-- | Check whether the matrix is an identity matrix.
--
-- > 1 0 0
-- > 0 1 0
-- > 0 0 1 (True)
isUnit m = isSquare m && allWithIndex (uncurry check) m
    where check = \i j e -> if i == j then e == 1 else e == 0

-- | Checks whether the matrix is empty.
--
-- > isEmpty m = numCols == 0 || numRows == 0
isEmpty m = numRows m == 0 || numCols m == 0

-- | Checks whether the matrix is a diagonal matrix.
--
-- > 4 0 0 0
-- > 0 7 0 0
-- > 0 0 3 0
-- > 0 0 0 9 (True)
isDiagonal m = isSquare m && allWithIndex (uncurry check) m
    where check = \i j e -> if i /= j then e == 0 else True

-- | Checks whether the matrix is a square matrix.
--
-- > isSquare == uncurry (==) . dimensions
isSquare m = let (a, b) = dimensions m in a == b

toDoubleMatrix :: (MatrixElement a, Integral a) => Matrix a -> Matrix Double
toDoubleMatrix = map fromIntegral

toRationalMatrix :: (MatrixElement a, Real a) => Matrix a -> Matrix Rational
toRationalMatrix = map toRational

toComplexMatrix :: (MatrixElement a, RealFloat a, Show a) => Matrix a -> Matrix (Complex a)
toComplexMatrix = map (:+ 0)

class Division e where
    divide :: e -> e -> e

instance Division Int     where divide = quot
instance Division Int8    where divide = quot
instance Division Int16   where divide = quot
instance Division Int32   where divide = quot
instance Division Int64   where divide = quot
instance Division Word8   where divide = quot
instance Division Word16  where divide = quot
instance Division Word32  where divide = quot
instance Division Word64  where divide = quot
instance Division Integer where divide = quot
instance Division Float   where divide = (/)
instance Division Double  where divide = (/)
instance Integral a => Division (Ratio a) where divide = (/)
instance RealFloat a => Division (Complex a) where divide = (/)

class (Eq e, Num e) => MatrixElement e where

    -- | Creates a matrix of the given size using a generator
    -- function for the value of each component.
    matrix :: (Int, Int) -> ((Int, Int) -> e) -> Matrix e

    -- | Builds a list from a matrix for the indices for which
    -- the given predicate matches.
    --
    -- > trace == select (uncurry (==))
    select :: ((Int, Int) -> Bool) -> Matrix e -> [e]

    -- | Returns the component at the given position in the matrix.
    -- Note that indices start at one, not at zero.
    at :: Matrix e -> (Int, Int) -> e

    -- | Returns the row at the given index in the matrix.
    -- Note that indices start at one, not at zero.
    row :: Int -> Matrix e -> [e]

    -- | Returns the row at the given index in the matrix.
    -- Note that indices start at one, not at zero.
    col :: Int -> Matrix e -> [e]

    -- | The dimensions of a given matrix.
    dimensions :: Matrix e -> (Int, Int)

    -- | The number of rows in the given matrix.
    numRows :: Matrix e -> Int

    -- | The number of columns in the given matrix.
    numCols :: Matrix e -> Int

    -- | Builds a matrix from a list of lists.
    --
    -- The innermost lists represent the rows. This function will create a m-n-matrix,
    -- where m is the number of rows, which is the minimum length of the row lists
    -- and n is the number of columns, i.e. the length of the outer list.
    --
    -- > fromList [[1,2,3],[2,1,3],[3,2,1]] :: Matrix Rational
    fromList :: [[e]] -> Matrix e

    -- | Turns a matrix into a list of lists.
    --
    -- > (toList . fromList) xs == xs
    --
    -- > (fromList . toList) mat == mat
    toList   :: Matrix e -> [[e]]

    -- | An identity square matrix of the given size.
    --
    -- >>> unit 4
    -- 1 0 0 0
    -- 0 1 0 0
    -- 0 0 1 0
    -- 0 0 0 1
    unit  :: Int -> Matrix e

    -- | A square matrix of the given size consisting of all zeros.
    --
    -- >>> zero 2
    -- 0 0
    -- 0 0
    zero  :: Int -> Matrix e

    -- | A square matrix which trace is the given list, all other components
    -- set to zero.
    --
    -- >>> diag [1,4,7,9]
    -- 1 0 0 0
    -- 0 4 0 0
    -- 0 0 7 0
    -- 0 0 0 9
    diag  :: [e] -> Matrix e

    -- | Check whether the matrix is the empty matrix.
    --
    -- > dimensions empty == (0, 0)
    empty :: Matrix e

    -- | Subtract a matrix from another.
    minus :: Matrix e -> Matrix e -> Matrix e

    -- | Add two matrices.
    --
    -- You may also use the 'Num' instance for matrices,
    -- i.e. @matrix1 + matrix2@ will work, too.
    plus  :: Matrix e -> Matrix e -> Matrix e

    -- | Multiply two matrices /O(n^3)/.
    --
    -- You may also use the 'Num' instance for matrices,
    -- i.e. @matrix1 * matrix2@ will work, too.
    times :: Matrix e -> Matrix e -> Matrix e

    -- | Compute the inverse of a matrix. Returns @Nothing@
    -- if the matrix is not invertible.
    inv   :: Matrix e -> Maybe (Matrix e)

    -- | Applies Bareiss multistep integer-preserving
    -- algorithm for finding the determinant of a matrix.
    -- Returns 0 if the matrix is not a square matrix.
    det       :: Matrix e -> e

    -- | Flip rows and columns.
    --
    -- > 1 8 9                1 2 3
    -- > 2 1 8  --transpose-> 8 1 2
    -- > 3 2 1                9 8 1
    transpose :: Matrix e -> Matrix e

    -- | Compute the rank of a matrix.
    rank      :: Matrix e -> e

    -- | Select the diagonal elements of a matrix as a list.
    --
    -- > 1 8 3
    -- > 3 6 5 --trace-> [1, 6, 2]
    -- > 7 4 2
    trace     :: Matrix e -> [e]

    -- | Select the minor of a matrix, that is the determinant
    -- of the 'minorMatrix'.
    --
    -- > minor = det . minorMatrix
    minor :: MatrixElement e => (Int, Int) -> Matrix e -> e

    -- | Select the minor matrix of a matrix, a matrix that is obtained
    -- by deleting the i-th row and j-th column.
    --
    -- > 10  9 95 45
    -- >  8  7  3 27                        8  3 27
    -- > 13 17 19 23 --minorMatrix (1,2)-> 13 19 23
    -- >  1  2  5  8                        1  5  8
    minorMatrix :: MatrixElement e => (Int, Int) -> Matrix e -> Matrix e

    cofactors :: MatrixElement e => Matrix e -> Matrix e
    adjugate :: MatrixElement e => Matrix e -> Matrix e

    -- | Apply a function on every component in the matrix.
    map :: MatrixElement f => (e -> f) -> Matrix e -> Matrix f

    -- | Apply a predicate on every component in the matrix
    -- and returns True iff all components satisfy it.
    all :: (e -> Bool) -> Matrix e -> Bool

    -- | Apply a predicate on every component in the matrix
    -- and return True if one or more components satisfy it.
    any :: (e -> Bool) -> Matrix e -> Bool

    mapWithIndex :: MatrixElement f => ((Int, Int) -> e -> f) -> Matrix e -> Matrix f
    allWithIndex :: ((Int, Int) -> e -> Bool) -> Matrix e -> Bool
    anyWithIndex :: ((Int, Int) -> e -> Bool) -> Matrix e -> Bool

    unit n  = fromList [[ if i == j then 1 else 0 | j <- [1..n]] | i <- [1..n] ]
    zero n  = matrix (n,n) (const 0)
    empty   = fromList []
    diag xs = matrix (n,n) (\(i,j) -> if i == j then xs !! (i-1) else 0)
      where n = length xs

    select p m = [ at m (i,j) | i <- [1..numRows m]
                              , j <- [1..numCols m]
                              , p (i,j) ]

    at mat (i, j) = ((!! j) . (!! i) . toList) mat

    row i = (!! (i-1)) . toList
    col i = row i . transpose

    numRows = fst . dimensions
    numCols = snd . dimensions
    dimensions m = case toList m of [] -> (0, 0)
                                    (x:xs) -> (length xs + 1, length x)

    adjugate = transpose . cofactors
    transpose mat = matrix (n, m) (\(i,j) -> mat `at` (j,i))
      where (m, n) = dimensions mat
    trace = select (uncurry (==))
    inv _ = Nothing

    minorMatrix (i,j) mat = matrix (numRows mat - 1, numCols mat - 1) $
                \(i',j') -> mat `at` (if i' >= i then i' + 1 else i',
                                      if j' >= j then j' + 1 else j')

    minor ix = det . minorMatrix ix

    cofactors mat = matrix (dimensions mat) $
       \(i,j) -> fromIntegral ((-1 :: Int)^(i+j)) * minor (i,j) mat

    map f = mapWithIndex (const f)
    all f = allWithIndex (const f)
    any f = anyWithIndex (const f)

    mapWithIndex f m = matrix (dimensions m) (\x -> f x (m `at` x))
    allWithIndex f m = P.all id [ f (i, j) (m `at` (i,j))
                                | i <- [1..numRows m], j <- [1..numCols m]]
    anyWithIndex f m = P.any id [ f (i, j) (m `at` (i,j))
                                | i <- [1..numRows m], j <- [1..numCols m]]

    a `plus` b
        | dimensions a /= dimensions b = error "Matrix.plus: dimensions don't match."
        | otherwise = matrix (dimensions a) (\x -> a `at` x + b `at` x)
    a `minus` b
        | dimensions a /= dimensions b = error "Matrix.minus: dimensions don't match."
        | otherwise = matrix (dimensions a) (\x -> a `at` x - b `at` x)
    a `times` b
        | numCols a /= numRows b = error "Matrix.times: `numRows a' and `numCols b' don't match."
        | otherwise = _mult a b

instance MatrixElement Int where
    matrix d g = runST (_matrix IntMatrix arrayST arraySTU d g)
    fromList = _fromList IntMatrix

    at         (IntMatrix _ _ arr) = _at arr
    dimensions (IntMatrix m n _)   = (m, n)
    row i      (IntMatrix _ _ arr) = _row i arr
    col j      (IntMatrix _ _ arr) = _col j arr
    toList     (IntMatrix _ _ arr) = _toList arr
    det        (IntMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)
    rank       (IntMatrix _ _ arr) = runST (_rank thawsBoxed arr)

instance MatrixElement Integer where
    matrix d g = runST (_matrix IntegerMatrix arrayST arrayST d g)
    fromList   = _fromList IntegerMatrix

    at         (IntegerMatrix _ _ arr) = _at arr
    dimensions (IntegerMatrix m n _)   = (m, n)
    row i      (IntegerMatrix _ _ arr) = _row i arr
    col j      (IntegerMatrix _ _ arr) = _col j arr
    toList     (IntegerMatrix _ _ arr) = _toList arr
    det        (IntegerMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr)
    rank       (IntegerMatrix _ _ arr) = runST (_rank thawsBoxed arr)

instance MatrixElement Float where
    matrix d g = runST (_matrix FloatMatrix arrayST arraySTU d g)
    fromList   = _fromList FloatMatrix

    at         (FloatMatrix _ _ arr) = _at arr
    dimensions (FloatMatrix m n _  ) = (m, n)
    row i      (FloatMatrix _ _ arr) = _row i arr
    col j      (FloatMatrix _ _ arr) = _col j arr
    toList     (FloatMatrix _ _ arr) = _toList arr
    det        (FloatMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)
    rank       (FloatMatrix _ _ arr) = runST (_rank thawsBoxed arr)
    inv        (FloatMatrix m n arr) = if m /= n then Nothing else
                                         let x = runST (_inv unboxedST pivotMax arr)
                                         in maybe Nothing (Just . FloatMatrix m n) x

instance MatrixElement Double where
    matrix d g = runST (_matrix DoubleMatrix arrayST arraySTU d g)
    fromList   = _fromList DoubleMatrix

    at         (DoubleMatrix _ _ arr) = _at arr
    dimensions (DoubleMatrix m n _  ) = (m, n)
    row i      (DoubleMatrix _ _ arr) = _row i arr
    col j      (DoubleMatrix _ _ arr) = _col j arr
    toList     (DoubleMatrix _ _ arr) = _toList arr
    det        (DoubleMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)
    rank       (DoubleMatrix _ _ arr) = runST (_rank thawsBoxed arr)
    inv        (DoubleMatrix m n arr) = if m /= n then Nothing else
                                         let x = runST (_inv unboxedST pivotMax arr)
                                         in maybe Nothing (Just . DoubleMatrix m n) x

instance (Show a, Integral a) => MatrixElement (Ratio a) where
    matrix d g = runST (_matrix RatioMatrix arrayST arrayST d g)
    fromList   = _fromList RatioMatrix

    at         (RatioMatrix _ _ arr) = _at arr
    dimensions (RatioMatrix m n _  ) = (m, n)
    row i      (RatioMatrix _ _ arr) = _row i arr
    col j      (RatioMatrix _ _ arr) = _col j arr
    toList     (RatioMatrix _ _ arr) = _toList arr
    det        (RatioMatrix m n arr) = if m /= n then 0 else  runST (_det thawsBoxed arr)
    rank       (RatioMatrix _ _ arr) = runST (_rank thawsBoxed arr)
    inv        (RatioMatrix m n arr) = if m /= n then Nothing else
                                        let x = runST (_inv boxedST pivotMax arr)
                                        in maybe Nothing (Just . RatioMatrix m n) x

instance (Show a, RealFloat a) => MatrixElement (Complex a) where
    matrix d g = runST (_matrix ComplexMatrix arrayST arrayST d g)
    fromList   = _fromList ComplexMatrix

    at         (ComplexMatrix _ _ arr) = _at arr
    dimensions (ComplexMatrix m n _  ) = (m, n)
    row i      (ComplexMatrix _ _ arr) = _row i arr
    col j      (ComplexMatrix _ _ arr) = _col j arr
    toList     (ComplexMatrix _ _ arr) = _toList arr
    det        (ComplexMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr)
    rank       (ComplexMatrix _ _ arr) = runST (_rank thawsBoxed arr)
    inv        (ComplexMatrix m n arr) = if m /= n then Nothing else
                                          let x = runST (_inv boxedST pivotNonZero arr)
                                          in maybe Nothing (Just . ComplexMatrix m n) x

_at :: (IArray a (u Int e), IArray u e)
    => a Int (u Int e) -> (Int, Int) -> e
_at arr (i,j) = arr ! i ! j

_row, _col :: (IArray a (u Int e), IArray u e) => Int -> a Int (u Int e) -> [e]
_row i arr = let row = arr ! i in [ row ! j | j <- [1..(snd (bounds row))] ]
_col j arr = [ arr ! i ! j | i <- [1..(snd (bounds arr))] ]

_toList :: (IArray a e) => Array Int (a Int e) -> [[e]]
_toList = P.map elems . elems

_fromList :: (IArray a (u Int e), IArray u e)
          => (Int -> Int -> a Int (u Int e) -> matrix e) -> [[e]] -> matrix e
_fromList c xs =
    let lengths = P.map length xs
        numCols = if null lengths then 0 else foldl1 min lengths
        numRows = length lengths

    in  c numRows numCols
          $ array (1, numRows)
          $ zip [1..numRows]
          $ P.map (array (1, numCols) . zip [1..numCols]) xs

thawsBoxed :: (IArray a e, MArray (STArray s) e (ST s))
           => Array Int (a Int e)
           -> ST s [STArray s Int e]
thawsBoxed = mapM thaw . elems

thawsUnboxed :: (IArray a e, MArray (STUArray s) e (ST s))
             => Array Int (a Int e)
             -> ST s [STUArray s Int e]
thawsUnboxed = mapM thaw . elems

arrays :: [(u s) Int e]
       -> ST s ((STArray s) Int ((u s) Int e))
arrays list = newListArray (1, length list) list

augment :: (IArray a e, MArray (u s) e (ST s), Num e)
        => ((Int, Int) -> [e] -> ST s ((u s) Int e))
        -> Array Int (a Int e)
        -> ST s (STArray s Int (u s Int e))
augment _ arr = do
    let (_, n) = bounds arr
        row (a,i) = newListArray (1, 2*n)
                                 [ if j > n then (if j == i + n then 1 else 0)
                                            else a ! j
                                 | j <- [1..2*n] ]

    mapM row (zip (elems arr) [1..]) >>= newListArray (1, n)

boxedST :: MArray (STArray s) e (ST s)
        => (Int, Int) -> [e] -> ST s ((STArray s) Int e)
boxedST = newListArray

unboxedST :: MArray (STUArray s) e (ST s)
          => (Int, Int) -> [e] -> ST s ((STUArray s) Int e)
unboxedST = newListArray

arrayST :: MArray (STArray s) e (ST s)
        => (Int, Int) -> e -> ST s ((STArray s) Int e)
arrayST = newArray

arraySTU :: MArray (STUArray s) e (ST s)
         => (Int, Int) -> e -> ST s ((STUArray s) Int e)
arraySTU = newArray

tee :: Monad m => (b -> m a) -> b -> m b
tee f x = f x >> return x

read :: (MArray a1 b m, MArray a (a1 Int b) m) =>
                       a Int (a1 Int b) -> Int -> Int -> m b
read a i j = readArray a i >>= flip readArray j

pivotMax :: Ord v => [(i, v)] -> i
pivotMax = fst . L.maximumBy (compare `on` snd)

pivotNonZero :: (Num v, Eq v) => [(i, v)] -> i
pivotNonZero xs = maybe (fst $ head xs) fst $ L.find ((/= 0) . snd) xs

_inv :: (IArray a e, MArray (u s) e (ST s), Fractional e, Show e, Eq e)
     => ((Int, Int) -> [e] -> ST s ((u s) Int e))
        -- ^ A function for building a new array
     -> ([(Int, e)] -> Int)
        -- ^ A function for selecting pivot elements
     -> Array Int (a Int e)
        -- ^ A matrix as arrays or arrays
     -> ST s (Maybe (Array Int (a Int e)))
_inv mkArrayST selectPivot mat = do
    let m = snd $ bounds mat
        n = 2 * m

        swap a i j = do
            tmp <- readArray a i
            readArray a j >>= writeArray a i
            writeArray a j tmp

    okay <- newSTRef True

    a <- augment mkArrayST mat

    forM_ [1..m] $ \k -> do
        iPivot <- selectPivot <$> zip [k..m]
                              <$> mapM (\i -> abs <$> read a i k) [k..m]

        p <- read a iPivot k
        if p == 0 then writeSTRef okay False else do

            swap a iPivot k

            forM_ [k+1..m] $ \i -> do
                a_i <- readArray a i
                a_k <- readArray a k
                forM_ [k+1..n] $ \j -> do
                    a_ij <- readArray a_i j
                    a_kj <- readArray a_k j
                    a_ik <- readArray a_i k
                    writeArray a_i j (a_ij - a_kj * (a_ik / p))
                writeArray a_i k 0

    invertible <- readSTRef okay

    if invertible then
      do
        forM_ [ m - v | v <- [0..m-1] ] $ \i -> do
            a_i <- readArray a i
            p   <- readArray a_i i
            writeArray a_i i 1
            forM_ [i+1..n] $ \j -> do
                readArray a_i j >>= writeArray a_i j . (/ p)

            unless (i == m) $ do
                forM_ [i+1..m] $ \k -> do
                    a_k <- readArray a k
                    p   <- readArray a_i k

                    forM_ [k..n] $ \j -> do
                        a_ij <- readArray a_i j
                        a_kj <- readArray a_k j
                        writeArray a_i j (a_ij - p * a_kj)

        mapM (\i -> readArray a i >>= getElems
                        >>= return . listArray (1, m) . drop m) [1..m]
            >>= return . Just . listArray (1, m)

      else return Nothing

_rank :: (IArray a e, MArray (u s) e (ST s), Num e, Division e, Eq e)
      => (Array Int (a Int e) -> ST s [(u s) Int e])
      -- ^ A function for thawing a boxed array
      -> Array Int (a Int e)
      -- ^ A matrix given as array of arrays
      -> ST s e
_rank thaws mat = do
    let m = snd $ bounds mat
        n = snd $ bounds (mat ! 1)

        swap a i j = do
            tmp <- readArray a i
            readArray a j >>= writeArray a i
            writeArray a j tmp

    a <- thaws mat >>= arrays

    ixPivot <- newSTRef 1
    prevR   <- newSTRef 1

    forM_ [1..n] $ \k -> do
        pivotRow <- readSTRef ixPivot

        switchRow <- mapM (\i -> read a i k) [pivotRow .. m]
            >>= return . L.findIndex (/= 0)

        when (isJust switchRow) $ do
            let ix = fromJust switchRow + pivotRow
            when (pivotRow /= ix) (swap a pivotRow ix)

            a_p   <- readArray a k
            pivot <- readArray a_p k
            prev  <- readSTRef prevR

            forM_ [pivotRow+1..m] $ \i -> do
                a_i <- readArray a i
                forM_ [k+1..n] $ \j -> do
                    a_ij <- readArray a_i j
                    a_ik <- readArray a_i k
                    a_pj <- readArray a_p j
                    writeArray a_i j ((pivot * a_ij - a_ik * a_pj)
                                        `divide` prev)

            writeSTRef ixPivot (pivotRow + 1)
            writeSTRef prevR pivot

    readSTRef ixPivot >>= return . (+ negate 1) . fromIntegral

_det :: (IArray a e, MArray (u s) e (ST s),
         Num e, Eq e, Division e)
     => (Array Int (a Int e) -> ST s [(u s) Int e])
     -> Array Int (a Int e) -> ST s e
_det thaws mat = do

    let size = snd $ bounds mat

    a <- thaws mat >>= arrays

    signR  <- newSTRef 1
    pivotR <- newSTRef 1

    forM_ [1..size] $ \k -> do
        sign <- readSTRef signR
        unless (sign == 0) $ do

            prev  <- readSTRef pivotR
            pivot <- read a k k >>= tee (writeSTRef pivotR)

            when (pivot == 0) $ do
                s <- forM [(k+1)..size] $ \r -> do
                    a_rk <- read a r k
                    if a_rk == 0 then return 0 else return r
                let sf = filter (>0) s

                when (not $ null sf) $ do
                    let sw = head sf

                    row <- readArray a sw
                    readArray a k >>= writeArray a sw
                    writeArray a k row

                    read a k k >>= writeSTRef pivotR
                    readSTRef signR >>= writeSTRef signR . negate

                when (null sf) (writeSTRef signR 0)

            sign' <- readSTRef signR
            unless (sign' == 0) $ do
                pivot' <- readSTRef pivotR
                forM_ [(k+1)..size] $ \i -> do
                    a_i <- readArray a i
                    forM [(k+1)..size] $ \j -> do
                        a_ij <- readArray a_i j
                        a_ik <- readArray a_i k
                        a_kj <- read a k j
                        writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)

    liftM2 (*) (readSTRef pivotR) (readSTRef signR)

-- TODO: More efficient implementation (decrease the constant factors
-- a little bit by working in the ST monad)
-- [ remark: not sure if this will be faster than lists -> benchmark! ]
_mult :: MatrixElement e => Matrix e -> Matrix e -> Matrix e
_mult a b = let rowsA = numRows a
                rowsB = numRows b
                colsB = numCols b
            in  matrix (rowsA, colsB) (\(i,j) -> L.foldl' (+) 0 [a `at` (i, k) * b `at` (k, j) | k <- [1..rowsB]])

_matrix :: (IArray a1 (u Int e), IArray u e,
            MArray a2 (u Int e) (ST s), MArray a3 e (ST s),
            Num e)
        => (Int -> Int -> a1 Int (u Int e) -> matrix)
        -> ((Int, Int) -> a -> ST s (a2 Int (u Int e)))
        -> ((Int, Int) -> e -> ST s (a3 Int e))
        -> (Int, Int)
        -> ((Int, Int) -> e)
        -> ST s matrix
_matrix c newArray newArrayU (m, n) g = do
    rows <- newArray (1, m) undefined
    forM_ [1..m] $ \i -> do
        cols <- newArrayU (1, n) 0
        forM_ [1..n] $ \j -> do
            writeArray cols j (g (i,j))
        U.unsafeFreeze cols >>= writeArray rows i
    U.unsafeFreeze rows >>= return . c m n

{-# RULES

"det/pow"
    forall a k. det (a ^ k) = (det a) ^ k

 #-}

max_stars_repo_name: phadej/bed-and-breakfast

max_issues_repo_name: phadej/bed-and-breakfast

max_forks_repo_name: phadej/bed-and-breakfast


text:
{-# language DeriveGeneric, DeriveAnyClass, StrictData #-}
{-# language KindSignatures, DataKinds, ScopedTypeVariables #-}
{-|
Module      : GA
Description : Genetic Algorithm for Poly Regression
Copyright   : (c) Fabricio Olivetti de Franca, 2021
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com
Stability   : experimental
Portability : POSIX
|-}

module GA
  ( generational
  , createRndSolution
  , crossover
  , mutate
  , Solution(..)
  , Poly(..)
  , Population
  , Fitness
  ) where

import System.Random
import Numeric.LinearAlgebra (Vector)

import Random
import Control.Monad
import Control.Monad.State.Strict
import GHC.Generics (Generic)
import Control.DeepSeq
import GHC.TypeLits
import Data.Proxy

-- | Data type representing a solution
newtype Poly (n :: Nat) = Poly { _getPoly :: (Bool, Int) } deriving (Show, Generic, NFData)

instance KnownNat n => RandomAlele (Poly n) where
  randomProb p x = Poly <$> biRandomWith p randomBool (randomInt (1, maxK)) (_getPoly x)
    where maxK = fromIntegral $ natVal (Proxy :: Proxy n)
  randomAlele    = Poly <$> randomPair (1, maxK)
    where maxK = fromIntegral $ natVal (Proxy :: Proxy n)

data Solution a = Sol { _chromo    :: [a]                    -- ^ chromossome representation
                      , _coeffs    :: Maybe (Vector Double)  -- ^ coefficients of the regression
                      , _fitness   :: Maybe Double           -- ^ fitness of an individual
                      } deriving (Show, Generic, NFData)

instance Eq (Solution a) where
  (Sol _ _ f1) == (Sol _ _ f2) = f1 == f2
instance Ord (Solution a) where
  (Sol _ _ f1) <= (Sol _ _ f2) = f1 <= f2

type Population a   = [Solution a]
type Fitness a      = Solution a -> Solution a

nVars :: Solution a -> Int
nVars s = length (_chromo s)

-- | Creates a random solution
createRndSolution :: RandomAlele a => Int -> Rnd (Solution a)
createRndSolution nLen =
  do chromo <- replicateM nLen randomAlele
     return $ Sol chromo Nothing Nothing

-- | Tournament between two solutions
tournament :: Population a -> Rnd (Solution a)
tournament pop = do s1 <- chooseRndFrom pop
                    s2 <- chooseRndFrom pop
                    return (min s1 s2)

-- | Creates a new solution from the combination of two solutions
crossover :: Population a -> Rnd (Solution a)
crossover pop = do s1 <- tournament pop
                   s2 <- tournament pop
                   ix <- randomInt (0, nVars s1)
                   return (cx ix s1 s2)

-- | Performs crossover
cx :: Int -> Solution a -> Solution a -> Solution a
cx ix sol1 sol2 = Sol child Nothing Nothing
  where
    recombine xs ys = take ix xs ++ drop ix ys
    child           = recombine (_chromo sol1) (_chromo sol2)

-- | Mutates a single solution
mutate :: RandomAlele a => Double -> Solution a -> Rnd (Solution a)
mutate prob sol = do chromo <- traverse (randomProb prob) (_chromo sol)
                     return $ Sol chromo Nothing Nothing

-- | Select the next generation
generational :: Population a -> Population a -> Rnd (Population a)
generational _ = return

max_stars_repo_name: folivetti/gapoly

max_issues_repo_name: folivetti/gapoly

max_forks_repo_name: folivetti/gapoly


text:
{-# LANGUAGE RankNTypes, BangPatterns, ScopedTypeVariables #-}

module Main where

import Math.Probably.MCMC
import qualified Math.Probably.PDF as PDF
import Math.Probably.RandIO
import Math.Probably.HamMC
import Math.Probably.MALA
import Math.Probably.FoldingStats
import Math.Probably.Sampler
import Control.Applicative
import Control.Monad

import Numeric.LinearAlgebra
import Numeric.AD
import Text.Printf

import Graphics.Gnewplot.Exec
import Graphics.Gnewplot.Types
import Graphics.Gnewplot.Style
import Graphics.Gnewplot.Panels
import Graphics.Gnewplot.Instances
import Graphics.Gnewplot.Histogram
import qualified Data.Vector.Storable as V

import Debug.Trace

import qualified Control.Monad.State.Strict as S

varf i = [10, 0.5, 1, 5, 0.1]!!(i `mod` 5)

corrf i j | i < j = corrf j i
          | otherwise = [-0.5, -0.2, 0, 0.2, 0.5]!!((i+j) `mod` 5)

covf (i,j) | i == j = varf i
           | otherwise = varf i * varf j * corrf i j

d = 3

muMVN = fromList $ replicate d 10
covMVN = buildMatrix d d covf

(invSigmaMVN, (lndetMVN, _)) = invlndet covMVN

pdfMVN = PDF.multiNormalByInv lndetMVN invSigmaMVN muMVN

--pdfI :: [Double] -> Double
pdfI xs = sum $ map f $ zip xs (map varf [0..]) where
  f (x, var) = PDF.normal 10 var x

pdfI2 xs = sum $ map f $ zip xs (map varf [0..]) where
  f (x, var) = PDF.logNormal 0.1 var x

pdf3 [alpha, error_sd] = sum $ map f regrdata where
  f (age, height) = PDF.normalLogPdf (alpha) error_sd  height

pdf2 [alpha, error_sd, beta] = sum $ map f regrdata where
  f (age, height) = PDF.normalLogPdf (alpha+(age*beta)) error_sd  height

regrdata :: (Real a, Floating a) => [(a,a)]
regrdata = [ 35.1**6.3,
 1.5**	2.1,
 35.2**	5.1,
 4.5**	3.1,
 10.1**	3.9,
 80.4**	9.1,
 5.2**	3.2,
 60.1**	8.1,
 70.1**	9.1,
 40.6**	7] where (**) = (,)

-- split HMC paper p 5
nealCov :: Matrix Double
neal_d = 200

nealCov = buildMatrix neal_d neal_d $ \(i,j) -> if i==j then realToFrac (i+1) / 100 else 0
--(2><2) [1, 0.95, 0.95, 2]
(nealCovInv, (nealLnDet,_)) = invlndet nealCov
nealMean = fromList $ replicate neal_d 3
nealPDF v = PDF.multiNormalByInv nealLnDet nealCovInv nealMean v
nealPostGrad v =
  let p = nealPDF v
      tol = 0.001
      f (i,x) = let h = max (tol/10) $ tol * x
                    v' = v V.// [(i,x+h)]
                    v'' = v V.// [(i,x-h)]
                in (nealPDF v' - nealPDF v'')/(2*h)
      grad3 = fromList $ map f $ zip [0..] $ toList v
  in (p,grad3)

covM1 :: Matrix Double
covM1 = (3><3) [1, 0,0,
                0 , 2,0,
                0, 0, 4]

means1, vars2 :: Vector Double
vars2 = fromList [1,2,4]

means1 = fromList [3,3,3]

(covMinv, (covMLnDet,_)) = invlndet covM1

main = runRIO $ do
  --io $ print $ nealCov
  --io $ print $ nealCovInv
{-  vs1 <- sample $ sequence $ replicate 10000 $ multiNormal means1 covM1
  io $ print $ runStat meanSDF $ map (@>0) $ vs1
  io $ print $ runStat meanSDF $ map (@>1) $ vs1
  io $ print $ runStat meanSDF $ map (@>2) $ vs1
  vs1 <- sample $ sequence $ replicate 10000 $ mvnSampler means1 2 vars2
  io $ print $ runStat meanSDF $ map (@>0) $ vs1
  io $ print $ runStat meanSDF $ map (@>1) $ vs1
  io $ print $ runStat meanSDF $ map (@>2) $ vs1
  io $ print $ PDF.multiNormal means1 (scale 2 covM1) (fromList [2,2,2])
--  io $ print $ PDF.multiNormalByInv covMLnDet covMinv means1 (fromList [2,2,2])
--  io $ print $ PDF.multiNormalIndep vars2 means1  (fromList [2,2,2])
  io $ print $ mvnPDF means1 2 vars2  (fromList [2,2,2]) -}
  let vinit = fromList $ replicate neal_d (1)
  --io $ print $ nealPostGrad vinit
  let pinit = fst $ nealPostGrad vinit
  io $ print $ ("pinit", pinit)
  let hmcp0 = HMCPar vinit 100 0 0.15 0 0 True
--  (vs') <- runMalaRioSimple (nealCov, nealCovInv, cholSH nealCov) nealPostGrad 2000 100 10 vinit
  (hmcp1, vs') <- runHMC nealPostGrad 100 hmcp0
  let vs = drop 50 $ reverse vs'
  --io $ forM (vs) $ \v -> print (v, nealPDF v)
  io $ print $ runStat meanSDF $ map (@>20) $ vs
  io $ print $ runStat meanSDF $ map (@>80) $ vs
  let have_ess = calcESSprim $ vs
  io  $ putStrLn $ "ESS=" ++show have_ess
  --io $ print hmcp1
  return ()

thin n [] = []
thin n (x:xs) = x : thin n (drop n xs)

covF i j = pure (-) <*> before meanF (\v -> v@>i * v@>j)
                    <*> fmap (uncurry (*)) (both (before meanF (@>i)) (before meanF (@>j)))

--cache prev, use grad'?

max_stars_repo_name: glutamate/probably

max_issues_repo_name: glutamate/probably

max_forks_repo_name: glutamate/probably


text:
module Main where

import Lib
import Data.Complex
import Graphics.Image

width = 750
height = 500
scaleDivisor = 250

xOffset = 2
yOffset = 1

isInMandelbrotSet :: Complex Double -> Complex Double -> Int -> Int
isInMandelbrotSet c z 0 = 0
isInMandelbrotSet c z itterations | (Data.Complex.magnitude z)  < 4 = isInMandelbrotSet c (z**2 + c) (itterations - 1)
                                  | (Data.Complex.magnitude z) >= 4 = itterations

calculatePixel :: Int -> Int -> Pixel RGB Double
calculatePixel y x = PixelRGB  0 0 value where
    real = (fromIntegral x / scaleDivisor) - xOffset
    imag = (fromIntegral y / scaleDivisor) - yOffset
    value = ((fromIntegral (isInMandelbrotSet (real :+ imag) 0 1000)) / 1000)

main = do
    let image = makeImageR VS (height, width) (\(y, x) -> calculatePixel y x)
    writeImageExact PNG [] "output.png" image

max_stars_repo_name: badnames/haskellbrot

max_issues_repo_name: badnames/haskellbrot

max_forks_repo_name: badnames/haskellbrot


text:
import qualified Numeric.BMML.RVMR             as RVM
import qualified Numeric.LinearAlgebra.HMatrix as H
import           System.Random.MWC             (asGenIO, uniformVector,
                                                withSystemRandom)

f :: Double -> Double
f x = 2.0 * x

generateX :: IO (H.Vector Double)
generateX =
    withSystemRandom . asGenIO $
    \gen ->
         uniformVector gen 12500

generateT :: H.Vector Double -> H.Vector Double
generateT = H.cmap f

main :: IO ()
main = do
  x <- generateX
  noiseX <- H.rand 12500 2
  let trainT = generateT x
      trainX = (H.fromColumns [x]) H.||| noiseX
  print (f 0.5)
  print (RVM.predict (RVM.fit trainX trainT) (H.fromList [0.5, 10.0, 10000.0]))

max_stars_repo_name: DbIHbKA/BMML

max_issues_repo_name: DbIHbKA/BMML

max_forks_repo_name: DbIHbKA/BMML


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
-- FIXME: PropSize and Switch not matched by semantics.
{-# OPTIONS_GHC -Wno-incomplete-patterns #-}

--
-- Copyright (c) 2019, ERICSSON AB
-- All rights reserved.
--
-- Redistribution and use in source and binary forms, with or without
-- modification, are permitted provided that the following conditions are met:
--
--     * Redistributions of source code must retain the above copyright notice,
--       this list of conditions and the following disclaimer.
--     * Redistributions in binary form must reproduce the above copyright
--       notice, this list of conditions and the following disclaimer in the
--       documentation and/or other materials provided with the distribution.
--     * Neither the name of the ERICSSON AB nor the names of its contributors
--       may be used to endorse or promote products derived from this software
--       without specific prior written permission.
--
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
-- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-- SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-- CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-- OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-- OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--

module Feldspar.Core.Semantics (Semantics(..), semantics) where

import Control.Monad
import Control.Monad.Par.Scheds.TraceInternal (yield)
import Data.Array.IArray
import Data.Array.IO
import Data.Array.Unsafe (unsafeFreeze)
import Data.Bits
import Data.Complex
import Data.Function (on)
import Data.IORef
import Data.List
import Feldspar.Core.Representation
import Feldspar.Core.NestedTuples
import Feldspar.Core.Types
import qualified Control.Exception as C
import qualified Control.Monad.Par as CMP
import System.IO.Unsafe

data Semantics a = Sem {semStr :: String, semSem :: a}

semantics :: Op a -> Semantics a
-- Feldspar.Core.Constructs.Array
semantics Append    = Sem "(++)" (++)
semantics GetIx     = Sem "(!)" evalGetIx
semantics GetLength = Sem "getLength" genericLength
semantics SetLength = Sem "setLength"
        (\n as -> genericTake n (as ++ repeat err))
      where
        err = error "reading uninitialized array element"
semantics Parallel = Sem "parallel"
        (\len ixf -> genericTake len $ map ixf [0..])
semantics Sequential = Sem "sequential"
        (\len i step -> genericTake len $
                        snd $ mapAccumL (\a ix -> swap (step ix a)) i [0..])
      where swap t = (nsnd t, nfst t)
semantics SetIx = Sem "setIx" evalSetIx
-- Feldspar.Core.Constructs.Binding
semantics Let = Sem "let" (flip ($))
-- Feldspar.Core.Constructs.Bits
semantics BAnd          = Sem "(.&.)"      (.&.)
semantics BOr           = Sem "(.|.)"      (.|.)
semantics BXor          = Sem "xor"        xor
semantics Complement    = Sem "complement" complement
semantics Bit           = Sem "bit"           (bit . fromIntegral)
semantics SetBit        = Sem "setBit"        (liftIntWord setBit)
semantics ClearBit      = Sem "clearBit"      (liftIntWord clearBit)
semantics ComplementBit = Sem "complementBit" (liftIntWord complementBit)
semantics TestBit       = Sem "testBit"       (liftIntWord testBit)
semantics ShiftLU       = Sem "shiftLU"     (liftIntWord shiftL)
semantics ShiftRU       = Sem "shiftRU"     (liftIntWord shiftR)
semantics ShiftL        = Sem "shiftL"      (liftInt shiftL)
semantics ShiftR        = Sem "shiftR"      (liftInt shiftR)
semantics RotateLU      = Sem "rotateLU"    (liftIntWord rotateL)
semantics RotateRU      = Sem "rotateRU"    (liftIntWord rotateR)
semantics RotateL       = Sem "rotateL"     (liftInt rotateL)
semantics RotateR       = Sem "rotateR"     (liftInt rotateR)
semantics ReverseBits   = Sem "reverseBits" evalReverseBits
semantics BitScan       = Sem "bitScan"  evalBitScan
semantics BitCount      = Sem "bitCount" evalBitCount
-- Feldspar.Core.Constructs.Complex
semantics MkComplex = Sem "complex"   (:+)
semantics RealPart  = Sem "creal"     realPart
semantics ImagPart  = Sem "cimag"     imagPart
semantics Conjugate = Sem "conjugate" conjugate
semantics MkPolar   = Sem "mkPolar"   mkPolar
semantics Magnitude = Sem "magnitude" magnitude
semantics Phase     = Sem "phase"     phase
semantics Cis       = Sem "cis"       cis
-- Feldspar.Core.Constructs.Condition
semantics Condition = Sem "condition" (\ c t e -> if c then t else e)
-- Feldspar.Core.Constructs.ConditionM
semantics ConditionM = Sem "if" ifM
      where
        ifM cond e t = if cond then e else t
-- Feldspar.Core.Constructs.Conversion
semantics F2I     = Sem "f2i"     truncate
semantics I2N     = Sem "i2n"     (fromInteger.toInteger)
semantics B2I     = Sem "b2i"     (\b -> if b then 1 else 0)
semantics Round   = Sem "round"   round
semantics Ceiling = Sem "ceiling" ceiling
semantics Floor   = Sem "floor"   floor
-- Feldspar.Core.Constructs.Elements
semantics EMaterialize    = Sem "materialize" ematerialize
semantics EWrite          = Sem "write" (\ix e -> Elements [(ix, e)])
semantics ESkip           = Sem "skip" (Elements [])
semantics EPar            = Sem "par" (\(Elements l) (Elements r) -> Elements (l ++ r))
semantics EparFor         = Sem "parFor" eparFor
-- Feldspar.Core.Constructs.Eq
semantics Equal    = Sem "(==)" (==)
semantics NotEqual = Sem "(/=)" (/=)
-- Feldspar.Core.Constructs.Error
semantics Undefined    = Sem "undefined" undefined
semantics (Assert msg) = Sem "assert"
        (\cond a -> if cond then a else error ("Assert failed: " ++ msg))
-- Feldspar.Core.Constructs.FFI
semantics (ForeignImport name f) = Sem name $ unEqBox f
-- Feldspar.Core.Constructs.Floating
semantics Pi      = Sem "pi"      Prelude.pi
semantics Exp     = Sem "exp"     Prelude.exp
semantics Sqrt    = Sem "sqrt"    Prelude.sqrt
semantics Log     = Sem "log"     Prelude.log
semantics Pow     = Sem "(**)"    (Prelude.**)
semantics LogBase = Sem "logBase" Prelude.logBase
semantics Sin     = Sem "sin"     Prelude.sin
semantics Tan     = Sem "tan"     Prelude.tan
semantics Cos     = Sem "cos"     Prelude.cos
semantics Asin    = Sem "asin"    Prelude.asin
semantics Atan    = Sem "atan"    Prelude.atan
semantics Acos    = Sem "acos"    Prelude.acos
semantics Sinh    = Sem "sinh"    Prelude.sinh
semantics Tanh    = Sem "tanh"    Prelude.tanh
semantics Cosh    = Sem "cosh"    Prelude.cosh
semantics Asinh   = Sem "asinh"   Prelude.asinh
semantics Atanh   = Sem "atanh"   Prelude.atanh
semantics Acosh   = Sem "acosh"   Prelude.acosh
-- Feldspar.Core.Constructs.Fractional
semantics DivFrac = Sem "(/)" (/)
-- Feldspar.Core.Constructs.Future
semantics MkFuture = Sem "future" FVal
semantics Await    = Sem "await"  unFVal
-- Feldspar.Core.Constructs.Integral
semantics Quot = Sem "quot" quot
semantics Rem  = Sem "rem"  rem
semantics Div  = Sem "div"  div
semantics Mod  = Sem "mod"  mod
semantics IExp = Sem "(^)"  (^)
-- Feldspar.Core.Constructs.Literal
semantics (Literal l) = Sem "lit" l
-- Feldspar.Core.Constructs.Logic
semantics And = Sem "(&&)" (&&)
semantics Or  = Sem "(||)" (||)
semantics Not = Sem "not"  not
-- Feldspar.Core.Constructs.Loop
semantics While = Sem "while" while
      where
        while cond body = do
                            c <- cond
                            when c (body >> while cond body)
semantics For = Sem "for" for
      where
        for 0 _    = return ()
        for l body = forM_ [0..l-1] body
semantics ForLoop = Sem "forLoop" forLoop
      where
        forLoop 0 initial _    = initial
        forLoop l initial body = foldl (flip body) initial [0..l-1]
semantics WhileLoop = Sem "whileLoop" whileLoop
      where
        whileLoop initial cond body = go initial
          where
            go st | cond st   = go $ body st
                  | otherwise = st
-- Feldspar.Core.Constructs.MutableArray
semantics NewArr    = Sem "newMArr"  $ \l -> newArray (mkBounds $ toInteger l)
semantics NewArr_   = Sem "newMArr_" $ \l -> newListArray (mkBounds $ toInteger l)
        [error $ "Undefined element at index " ++ show (i::Integer) | i <- [0..]]
semantics GetArr    = Sem "getMArr"  $ \arr i -> readArray arr (toInteger i)
semantics SetArr    = Sem "setMArr"  $ \arr i -> writeArray arr (toInteger i)
semantics ArrLength = Sem "arrLength" (getBounds >=> \(l,u) -> return $ fromInteger (u-l+1))
-- Feldspar.Core.Constructs.Mutable
semantics Run = Sem "runMutable" unsafePerformIO
semantics Bind = Sem "(>>=)" (>>=)
semantics Then = Sem "(>>)" (>>)
semantics When = Sem "when" when
semantics Return = Sem "return" return
-- Feldspar.Core.Constructs.MutableReference
semantics NewRef = Sem "newRef" newIORef
semantics GetRef = Sem "getRef" readIORef
semantics SetRef = Sem "setRef" writeIORef
semantics ModRef = Sem "modRef" (\r f -> readIORef r >>= writeIORef r . f)
-- Feldspar.Core.Constructs.MutableToPure
semantics RunMutableArray = Sem "runMutableArray" runMutableArrayEval
semantics WithArray       = Sem "withArray"       withArrayEval
-- Nested tuples
semantics Cons  = Sem "cons"  (:*)
semantics Nil   = Sem "nil"   TNil
semantics Car   = Sem "car"   (\ (x :* _) -> x)
semantics Cdr   = Sem "cdr"   (\ (_ :* y) -> y)
semantics Tup   = Sem "tup"   id
-- Feldspar.Core.Constructs.NoInline
semantics NoInline  = Sem "NoInline" id
-- Feldspar.Core.Constructs.Num
semantics Abs  = Sem "abs" abs
semantics Sign = Sem "signum" signum
semantics Add  = Sem "(+)" (+)
semantics Sub  = Sem "(-)" (-)
semantics Mul  = Sem "(*)" (*)
-- Feldspar.Core.Constructs.Ord
semantics LTH = Sem "(<)"  (<)
semantics GTH = Sem "(>)"  (>)
semantics LTE = Sem "(<=)" (<=)
semantics GTE = Sem "(>=)" (>=)
semantics Min = Sem "min"  min
semantics Max = Sem "max"  max
-- Feldspar.Core.Constructs.Par
semantics ParRun    = Sem "runPar" CMP.runPar
semantics ParNew    = Sem "new" CMP.new
semantics ParGet    = Sem "get" CMP.get
semantics ParPut    = Sem "put" CMP.put_
semantics ParFork   = Sem "fork" CMP.fork
semantics ParYield  = Sem "yield" yield
-- Feldspar.Core.Constructs.RealFloat
semantics Atan2   = Sem "atan2" Prelude.atan2
-- Feldspar.Core.Constructs.Save
semantics Save = Sem "save" id

-- | Support for Array
evalGetIx :: (Integral a, Show a) => [p] -> a -> p
evalGetIx as i
    | 0 <= i && i < len = genericIndex as i
    | otherwise = error $ unwords
        [ "getIx: accessing index"
        , show i
        , "outside the bounds of an array of length"
        , show len
        ]
    where
      len = genericLength as

evalSetIx :: (Integral b, Show b) => [a] -> b -> a -> [a]
evalSetIx as i v
    | 0 <= i && i < len = genericTake i as ++ [v] ++ genericDrop (i+1) as
    | otherwise = error $ unwords
        [ "setIx: assigning index"
        , show i
        , "outside the bounds of an array of length"
        , show len
        ]
    where
      len = genericLength as

-- | Support for Bits
liftIntWord :: (a -> Int -> b) -> a -> WordN -> b
liftIntWord f x = f x . fromIntegral

liftInt :: (a -> Int -> b) -> a -> IntN -> b
liftInt f x = f x . fromIntegral

evalReverseBits :: (Num b, FiniteBits b) => b -> b
evalReverseBits b = revLoop b 0 (0 `asTypeOf` b)
  where
    bSz = finiteBitSize b
    revLoop x i n | i >= bSz    = n
                  | testBit x i = revLoop x (i+1) (setBit n (bSz - i - 1))
                  | otherwise   = revLoop x (i+1) n

evalBitScan :: (FiniteBits b) => b -> WordN
evalBitScan b =
   if isSigned b
   then scanLoop b (testBit b (finiteBitSize b - 1)) (finiteBitSize b - 2) 0
   else scanLoop b False (finiteBitSize b - 1) 0
  where
    scanLoop x t i n | i Prelude.< 0            = n
                     | testBit x i Prelude./= t = n
                     | otherwise                = scanLoop x t (i-1) (n+1)

evalBitCount :: (FiniteBits b) => b -> WordN
evalBitCount b = loop b (finiteBitSize b - 1) 0
  where
    loop x i n | i Prelude.< 0 = n
               | testBit x i   = loop x (i-1) (n+1)
               | otherwise     = loop x (i-1) n

-- | Support for Elements
ematerialize :: Length -> Elements a -> [a]
ematerialize l (Elements xs) = map snd xs'
  where xs' = genericTake l $ sortBy (compare `on` fst) xs

eparFor :: Length -> (Index -> Elements a) -> Elements a
eparFor len ixf = Elements $ concatMap (\(Elements vs) -> vs) xs
      where xs = genericTake len $ map ixf [0..]

-- | Support for MultableArray
-- | Calculate array bounds. If the length is zero, flip the arguments to
-- make an empty range
mkBounds :: Integer -> (Integer,Integer)
mkBounds l = (0, pred l)

-- | Support for MutableToPure
runMutableArrayEval :: forall a . Mut (MArr a) -> [a]
runMutableArrayEval m = unsafePerformIO $
                        do marr <- m
                           iarr <- unsafeFreeze marr
                           return (elems (iarr :: Array Integer a))

withArrayEval :: forall a b. MArr a -> ([a] -> Mut b) -> Mut b
withArrayEval ma f
    = do a <- f (elems (unsafePerformIO $ freeze ma :: Array Integer a))
         C.evaluate a

max_stars_repo_name: Feldspar/feldspar-language

max_issues_repo_name: Feldspar/feldspar-language

max_forks_repo_name: Feldspar/feldspar-language


text:
{-# LANGUAGE DeriveDataTypeable    #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards       #-}

-- |
-- Module    : Statistics.Distribution.Normal
-- Copyright : (c) 2021 Alexander E. Zarebski
-- License   : MIT
--
-- Maintainer  : aezarebski@gmail.com
-- Stability   : experimental
-- Portability : portable
--
-- The lognormal distribution.

module Statistics.Distribution.Lognormal
   (
     LognormalDistribution
   -- * constructors
   , lognormalDistr
   , lognormalDistrE
   ) where

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

import qualified Statistics.Distribution         as D

-- | The lognormal distribution.
data LognormalDistribution = LND {
      meanlog :: {-# UNPACK #-} !Double
    , sdlog   :: {-# UNPACK #-} !Double
    } deriving (Eq, Typeable, Data, Generic)

instance D.Distribution LognormalDistribution where
    cumulative      = cumulative

instance D.ContDistr LognormalDistribution where
    logDensity    = logDensity
    quantile      = quantile

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

instance D.Mean LognormalDistribution where
    mean LND {..} = exp (meanlog + 0.5 * sdlog ** 2.0)

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

instance D.Variance LognormalDistribution where
    variance LND {..} = exp (2.0 * meanlog + sdlog ** 2.0) * (exp (sdlog ** 2.0) - 1)

-- | Create lognormal distribution from parameters.
lognormalDistr :: Double            -- ^ Mean of the logarithm
            -> Double            -- ^ Standard deviation of the logarithm
            -> LognormalDistribution
lognormalDistr m sd = maybe (error $ errMsg m sd) id $ lognormalDistrE m sd

-- | Create normal distribution from parameters.
lognormalDistrE :: Double            -- ^ Mean of logarithm
             -> Double            -- ^ Standard deviation of logarithm
             -> Maybe LognormalDistribution
lognormalDistrE mln sdln
  | sdln > 0    = Just LND { meanlog       = mln
                           , sdlog = sdln
                           }
  | otherwise = Nothing

errMsg :: Double -> Double -> String
errMsg _ sd = "Statistics.Distribution.Lognormal.lognormalDistr: standard deviation must be positive. Got " ++ show sd

logDensity :: LognormalDistribution -> Double -> Double
logDensity LND {..} x = - (log x - meanlog) ** 2.0 / (2 * sdlog ** 2.0) - log (x * sdlog * m_sqrt_2_pi)

cumulative :: LognormalDistribution -> Double -> Double
cumulative LND {..} x = 0.5 * (1 + erf ((log x - meanlog) / (sdlog * m_sqrt_2)))

quantile :: LognormalDistribution -> Double -> Double
quantile LND {..} p
  | p == 0 = 0
  | p == 1 = inf
  | p == 0.5 = exp meanlog
  | p > 0 && p < 1 = exp (meanlog + sdlog * m_sqrt_2 * invErf (2 * p - 1))
  | otherwise = error $ "Statistics.Distribution.Lognormal.quantile: p must be in [0,1] range. Got: "++show p
  where inf = 1/0

max_stars_repo_name: aezarebski/timtam

max_issues_repo_name: aezarebski/timtam

max_forks_repo_name: aezarebski/timtam


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TupleSections         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE FlexibleContexts      #-}

import           Control.Applicative
import           Control.Monad
import           Control.Monad.Random

import           Codec.Compression.GZip ( decompress )
import           Data.Serialize ( Get )
import qualified Data.Serialize as Serialize
import qualified Data.ByteString.Lazy as B

import           Data.List ( foldl' )
import           Data.List.Split ( chunksOf )
import           Data.Maybe ( fromMaybe )
#if ! MIN_VERSION_base(4,13,0)
import           Data.Semigroup ( (<>) )
#endif
import           Data.Word ( Word32 , Word8 )
import qualified Data.Vector.Storable as V

import           Numeric.LinearAlgebra ( maxIndex )
import qualified Numeric.LinearAlgebra.Static as SA

import           Options.Applicative
import           System.FilePath ( (</>) )

import           Grenade
import           Grenade.Utils.OneHot

#ifdef COMPLEX
-- It's logistic regression!
--
-- This network is used to show how we can embed a Network as a layer in the larger MNIST
-- type.
--
type FL i o =
  Network
    '[ FullyConnected i o, Logit ]
    '[ 'D1 i, 'D1 o, 'D1 o ]

-- The definition of our convolutional neural network.
-- In the type signature, we have a type level list of shapes which are passed between the layers.
-- One can see that the images we are inputing are two dimensional with 28 * 28 pixels.

-- It's important to keep the type signatures, as there's many layers which can "squeeze" into the gaps
-- between the shapes, so inference can't do it all for us.

-- With the mnist data from Kaggle normalised to doubles between 0 and 1, learning rate of 0.01 and 15 iterations,
-- this network should get down to about a 1.3% error rate.
--
-- /NOTE:/ This model is actually too complex for MNIST, and one should use the type given belowe instead.
--         This one is just here to demonstrate Inception layers in use.
--
type MNIST =
  Network
    '[ Reshape,
       Concat ('D3 28 28 1) Trivial ('D3 28 28 14) (InceptionMini 28 28 1 5 9),
       Pooling 2 2 2 2, Relu,
       Concat ('D3 14 14 3) (Convolution 15 3 1 1 1 1) ('D3 14 14 15) (InceptionMini 14 14 15 5 10), Crop 1 1 1 1, Pooling 3 3 3 3, Relu,
       Reshape, FL 288 80, FL 80 10 ]
    '[ 'D2 28 28, 'D3 28 28 1,
       'D3 28 28 15, 'D3 14 14 15, 'D3 14 14 15, 'D3 14 14 18,
       'D3 12 12 18, 'D3 4 4 18, 'D3 4 4 18,
       'D1 288, 'D1 80, 'D1 10 ]

#else

-- The simpler network can just be dropped in without changing any of the other code
--

type MNIST
  = Network
    '[ Convolution 1 10 5 5 1 1, Pooling 2 2 2 2, Relu
     , Convolution 10 16 5 5 1 1, Pooling 2 2 2 2, Reshape, Relu
     , FullyConnected 256 80, Logit, FullyConnected 80 10, Logit]
    '[ 'D2 28 28, 'D3 24 24 10, 'D3 12 12 10, 'D3 12 12 10
     , 'D3 8 8 16, 'D3 4 4 16, 'D1 256, 'D1 256
     , 'D1 80, 'D1 80, 'D1 10, 'D1 10]

-- ... and this is an even simpler network from  https://crypto.stanford.edu/~blynn/haskell/brain.html
{-
type MNIST
  = Network
      '[Reshape, FullyConnected 784 30, Relu, FullyConnected 30 10, Logit]
      '[ 'D2 28 28, 'D1 784, 'D1 30, 'D1 30, 'D1 10, 'D1 10]
-}

#endif

randomMnist :: MonadRandom m => m MNIST
randomMnist = randomNetwork

convTest :: Int -> FilePath -> Maybe Int -> LearningParameters -> IO ()
convTest iterations dataDir nSamples rate = do
  net0         <- randomMnist
  trainData    <- readMNIST (dataDir </> "train-images-idx3-ubyte.gz")
                            (dataDir </> "train-labels-idx1-ubyte.gz")
  validateData <- readMNIST (dataDir </> "t10k-images-idx3-ubyte.gz")
                            (dataDir </> "t10k-labels-idx1-ubyte.gz")

  foldM_ (runIteration (maybe trainData (`take` trainData) nSamples) validateData) net0 [1..iterations]

    where
  trainEach rate' !network (i, o) = train rate' network i o

  runIteration trainRows validateRows net i = do
    let trained' = foldl' (trainEach ( rate { learningRate = learningRate rate * 0.9 ^ i} )) net trainRows
    print trained'

    putStrLn "Checking..."
    let res      = fmap (\(rowP,rowL) -> (rowL,) $ runNet trained' rowP) validateRows
    let res'     = fmap (\(S1D label, S1D prediction) -> (maxIndex (SA.extract label), maxIndex (SA.extract prediction))) res
    let matched   = length $ filter ((==) <$> fst <*> snd) res'
    let total     = length res'
    let matchedpc = fromIntegral matched / fromIntegral total * 100.0 :: Float
    putStrLn $ "Iteration " ++ show i ++ ": matched " ++ show matched ++ " of " ++ show total ++ " (" ++ show matchedpc ++ "%)"
    return trained'

data MnistOpts = MnistOpts FilePath (Maybe Int) Int LearningParameters

mnist' :: Parser MnistOpts
mnist' = MnistOpts <$> argument str (metavar "DATADIR")
                       -- option to reduce the number of training samples used from 60,000
                       -- to avoid running out of memory
                   <*> option (Just <$> auto) (long "limit_samples_to" <> short 'l' <> value Nothing)
                   <*> option auto (long "iterations" <> short 'i' <> value 15)
                   <*> (LearningParameters
                       <$> option auto (long "train_rate" <> short 'r' <> value 0.01)
                       <*> option auto (long "momentum" <> value 0.9)
                       <*> option auto (long "l2" <> value 0.0005)
                       )

main :: IO ()
main = do
    MnistOpts dataDir nSamples iter rate <- execParser (info (mnist' <**> helper) idm)
    putStr "Training convolutional neural network with "
    putStr $ maybe "all" show nSamples
    putStrLn " samples..."

    convTest iter dataDir nSamples rate

-- Adapted from https://github.com/tensorflow/haskell/blob/master/tensorflow-mnist/src/TensorFlow/Examples/MNIST/Parse.hs
-- Could also have used Data.IDX, although that uses a different Vector variant from that need for fromStorable
readMNIST :: FilePath -> FilePath -> IO [(S ( 'D2 28 28), S ( 'D1 10))]
readMNIST iFP lFP = do
  labels  <- readMNISTLabels lFP
  samples <- readMNISTSamples iFP
  return $ zip
    (fmap (fromMaybe (error "bad samples") . fromStorable) samples)
    (fromMaybe (error "bad labels") . oneHot . fromIntegral <$> labels)

-- | Check's the file's endianess, throwing an error if it's not as expected.
checkEndian :: Get ()
checkEndian = do
  magic <- Serialize.getWord32be
  when (magic `notElem` ([2049, 2051] :: [Word32]))
    $ error "Expected big endian, but image file is little endian."

-- | Reads an MNIST file and returns a list of samples.
readMNISTSamples :: FilePath -> IO [V.Vector Double]
readMNISTSamples path = do
  raw <- decompress <$> B.readFile path
  either fail ( return . fmap (V.map normalize) ) $ Serialize.runGetLazy getMNIST raw
 where
  getMNIST :: Get [V.Vector Word8]
  getMNIST = do
    checkEndian
    -- Parse header data.
    cnt    <- fromIntegral <$> Serialize.getWord32be
    rows   <- fromIntegral <$> Serialize.getWord32be
    cols   <- fromIntegral <$> Serialize.getWord32be
    -- Read all of the data, then split into samples.
    pixels <- Serialize.getLazyByteString $ fromIntegral $ cnt * rows * cols
    return $ V.fromList <$> chunksOf (rows * cols) (B.unpack pixels)

  normalize :: Word8 -> Double
  normalize = (/ 255) . fromIntegral
  -- There are other normalization functions in the literature, such as
  -- normalize = (/ 0.3081) . (`subtract` 0.1307) . (/ 255) . fromIntegral
  -- but we need values in the range [0..1] for the showShape' pretty printer

-- | Reads a list of MNIST labels from a file and returns them.
readMNISTLabels :: FilePath -> IO [Word8]
readMNISTLabels path = do
  raw <- decompress <$> B.readFile path
  either fail return $ Serialize.runGetLazy getLabels raw
 where
  getLabels :: Get [Word8]
  getLabels = do
    checkEndian
    -- Parse header data.
    cnt <- fromIntegral <$> Serialize.getWord32be
    -- Read all of the labels.
    B.unpack <$> Serialize.getLazyByteString cnt