haskell - 如何概括抽样框架？

Question

在随机射线跟踪器的情况下，我想将MC集成（路径跟踪，双向路径跟踪）与样本生成（均匀随机，分层，泊松，都会，...）分离。其中大多数已经实现，但是使用起来很繁琐。因此，我放弃了这一点，并尝试通过将采样的计算分为两个阶段来构建更好的东西：在SampleGen中，允许您使用mk1d和mk2d函数请求随机值，然后为它们提供实际的Float由采样算法决定。可以在SampleRun中检查这些值以进行实际计算。这是一些带有分层采样器有趣部分的代码，它的使用是：

{-# LANGUAGE GeneralizedNewtypeDeriving #-}

import Control.Applicative
import Control.Monad.State.Strict
import Control.Monad.Primitive
import System.Random.MWC as MWC

-- allows to construct sampled computations
newtype SampleGen s m a = SampleGen (StateT s m a)
                       deriving ( Functor, Applicative, Monad
                                , MonadState s, MonadTrans )

-- allows to evaluate sampled computations constructed in SampleGen
newtype SampleRun s m a = SampleRun (StateT s m a)
                       deriving ( Functor, Applicative, Monad
                                , MonadState s )

-- a sampled computation, parametrized over the generator's state g,
-- the evaluator's state r,  the underlying monad m and the result
-- type a
type Sampled g r m a = SampleGen g m (SampleRun r m a)

----------------------
-- Stratified Sampling
----------------------

-- | we just count the number of requested 1D samples
type StratGen = Int

-- | the pre-computed values and a RNG for additional ones
type StratRun m = ([Float], Gen (PrimState m))

-- | specialization of Sampled for stratified sampling
type Stratified m a = Sampled StratGen (StratRun m) m a

-- | gives a sampled value in [0..1), this is kind
--   of the "prime" value, upon which all computations
--   are built
mk1d :: PrimMonad m => Stratified m Float
mk1d = do
  n1d <- get
  put $ n1d + 1

  return $ SampleRun $ do
    fs <- gets fst
    if length fs > n1d
      then return (fs !! n1d)
      else gets snd >>= lift . MWC.uniform

-- | gives a pair of stratified values, should really also
--   be a "prime" value, but here we just construct them
--   from two 1D samples for fun
mk2d :: (Functor m, PrimMonad m) => Stratified m (Float, Float)
mk2d = mk1d >>= \f1 -> mk1d >>= \f2 ->
  return $ (,) <$> f1 <*> f2

-- | evaluates a stratified computation
runStratified
  :: (PrimMonad m)
  => Int            -- ^ number of samples
  -> Stratified m a -- ^ computation to evaluate
  -> m [a]          -- ^ the values produced, a list of nsamples values
runStratified nsamples (SampleGen c) = do
  (SampleRun x, n1d) <- runStateT c 0
  -- let's just pretend I'd use n1d to actually
  -- compute stratified samples
  gen <- MWC.create
  replicateM nsamples $ evalStateT x ([{- samples would go here #-}], gen)

-- estimate Pi by Monte Carlo sampling
-- mcPi :: (Functor m, PrimMonad m) => Sampled g r m Float
mcPi :: (Functor m, PrimMonad m) => Stratified m Float
mcPi = do
  v <- mk2d
  return $ v >>= \(x, y) -> return $ if x * x + y * y < 1 then 4 else 0

main :: IO ()
main = do
  vs <- runStratified 10000 mcPi :: IO [Float]
  print $ sum vs / fromIntegral (length vs)

这里缺少的部分是 mcPi函数的当前形式具有类型

mcPi :: (Functor m, PrimMonad m) => Stratified m Float

而实际上应该是这样的

mcPi :: (Functor m, PrimMonad m) => Sampled g r m Float

公认的是， Sampled上的四个类型参数并不是很漂亮，但是至少类似这样会有用。总而言之，我正在寻找可以独立于采样算法来表达 mcPi之类的计算的东西，例如：


统一的随机采样器不需要在 SampleGen阶段保持任何状态，而在 SampleRun阶段仅需要RNG
分层磁盘采样器和泊松磁盘采样器（可能还有其他）都跟踪所需的1D和2D样本的数量，并将它们预先计算为向量，然后将允许它们共享 SampleGen和 SampleRun实现，以便仅在 SampleGen和 SampleRun之间发生的事情不同（向量的实际填充方式）
大都市采样器在其 SampleRun阶段将使用 lazy sample generation技术


我想使用GHC进行编译，因此扩展名 MultiParamTypeClasses和 TypeFamilies对我来说还可以，但是我没有提出任何可以远程使用的东西。

PS：作为动力，有些 pretty pictures。当前形式的代码在 GitHub上

Answer 1

我将从一个截然不同的问题开始，即“代码应该是什么样子？”，然后朝着“采样框架如何组合在一起”这个问题着手。

代码应该是什么样子

mcPi的定义应为

mcPi :: (Num s, Num p) => s -> s -> p
mcPi x y = if x * x + y * y < 1 then 4 else 0

pi的蒙特卡洛估计是，给定两个数字（恰好来自间隔[0..1]），如果pi落在一个圆内，则pi是正方形的面积，否则为0。 pi对计算一无所知。它不知道是否要重复，或有关数字的来源。它确实知道数字应该在正方形上均匀分布，但这是另一个问题的话题。 pi的蒙特卡洛估计只是从样本到估计的函数。

其他随机事物将知道它们是随机过程的一部分。一个简单的随机过程可能是：掷硬币，如果硬币出现“正面”，请再次掷硬币。

simpleRandomProcess :: (Monad m, MonadCoinFlip m) => m Coin
simpleRandomProcess =
    do
        firstFlip <- flipACoin
        case firstFlip of 
            Heads -> flipACoin
            Tails -> firstFlip

这个随机过程希望能够看到类似

data Coin = Heads | Tails

class MonadCoinFlip m where
    flipACoin :: m Coin -- The coin should be fair

随机过程可以根据之前的实验结果更改所需的随机数据量。这表明我们最终将需要提供一个 Monad。

介面

您想“将MC集成（路径跟踪，双向路径跟踪）与样本生成（均匀随机，分层，泊松，都会，...）分离开”。在您的示例中，他们都想采样浮点数。这表明以下课程

class MonadSample m where
    sample :: m Float -- Should be on the interval [0..1)

除了两件事之外，这与 existing MonadRandom class非常相似。 MonadRandom实现本质上需要在自己选择的某个范围内提供统一的随机 Int。您的采样器将在间隔[0..1）上提供未知分布的 Float样本。这足够不同，足以证明拥有自己的新课程。

由于即将出现的 Monad Applicative更改，我将为该类建议一个不同的名称 SampleSource。

class SampleSource f where
    sample :: f Float -- Should be on the interval [0..1)

sample替换代码中的 mk1d。 mk2d也可以被替换，同样也不知道样本的来源。 sample2d的替代品 mk2d可以与任何 Applicative样本源一起使用，不需要成为 Monad。不需要 Monad的原因是，它不会根据采样结果来决定要获取多少个样本，或者要做什么？它的计算结构是提前知道的。

sample2d :: (Applicative f, SampleSource f) => f (Float, Float)
sample2d = (,) <$> sample <*> sample

如果要允许样本源引入维之间的相互作用，例如用于泊松磁盘采样，则需要将其添加到接口中，或者显式枚举维

class SampleSource f where
    sample   :: f Float
    sample2d :: f (Float, Float)
    sample3d :: f (Float, Float, Float)
    sample4d :: f (Float, Float, Float, Float)

或使用某些向量库。

class SampleSource f where
    sample  :: f Float
    samples :: Int -> f (Vector Float)

实施界面

现在，我们需要描述如何将每个样本源用作 SampleSource。例如，我将为最差的示例源之一实现 SampleSource。

newtype ZeroSampleSourceT m a = ZeroSampleSourceT {
    unZeroSampleSourceT :: IdentityT m a
} deriving (MonadTrans, Monad, Functor, MonadPlus, Applicative, Alternative, MonadIO)

instance (Monad m) => SampleSource (ZeroSampleSourceT m a) where
    sample = return 0

runZeroSampleSourceT :: (Monad m) => ZeroSampleSourceT m a -> m a
runZeroSampleSourceT = runIdentityT . unZeroSampleSourceT

当所有 Monad是 Applicative时，我改写

instance (Applicative f) => SampleSource (ZeroSampleSourceT f) where
    sample = pure 0

我还将实现MWC统一 SampleSource。

newtype MWCUniformSampleSourceT m a = MWCUniformSampleSourceT m a {
    unMWCUniformSampleSourceT :: ReaderT (Gen (PrimState m)) m a
} deriving (MonadTrans, Monad, Functor, MonadPlus, Applicative, Alternative, MonadIO)

runMWCUniformSampleSourceT :: MWCUniformSampleSourceT m a -> (Gen (PrimState m)) -> m a
runMWCUniformSampleSourceT = runReaderT . unMWCUniformSampleSourceT

-- MWC's uniform generates floats on the open-closed interval (0,1]
uniformClosedOpen :: PrimMonad m => Gen (PrimState m) -> m Float
uniformClosedOpen = fmap (\x -> x - 2**(-33)) . uniform

instance (PrimMonad m) => SampleSource (MWCUniformSampleSourceT m) where
    sample = MWCUniformSampleSourceT . ReaderT $ uniformClosedOpen

我们不会完全实现 Stratified或 runStratified，因为您的示例代码未包含针对它们的完整实现。

但我想知道会提前使用多少个样本

我不确定您要对“分层”抽样采取什么措施。我不知道分层抽样是预先生成数字，并在数字用尽时使用生成器。如果要为某事物提供单子接口，则将无法提前告知将执行什么，因此您将无法在开始执行计算之前预测一次计算需要多少个样本。如果您只能满足于 Applicative接口，则可以提前测试整个计算将需要多少个样本。

但是泊松磁盘采样需要提前知道要采样多少点

如果单个采样可以同时取决于所需的样本数量和维数，例如在泊松磁盘采样中，则在已知采样后就需要将其传递给采样器。

class SampleSource f where
    sample   :: f Float
    samples  :: Int -> f ([Float])
    sampleN  :: Int -> f (Vector Float)
    samplesN :: Int -> Int -> f ([Vector Float])

您可以将其推广到任意尺寸的任意形状的样本中，这是我们下一步的工作。

带有Monadic解释器的应用查询语言

我们可以非常非常详尽地为样本请求创建 Applicative查询语言。该语言将需要在 Applicative已有功能的基础上添加两个功能。它需要能够重复请求，并且需要将样本请求分组在一起，以识别哪些分组是有意义的。它受以下代码的激励，该代码想要获取6个不同的2d样本，其中 sample2d与我们的第一个定义相同。

take 6 (repeat sample2d)

首先，我们需要能够重复一遍又一遍。最好的方法是如果我们可以写，例如

take 6 (repeat sample) :: SampleSource f => [f Float]

我们需要一种从 [f a]转到 f [a]的方法。这已经存在；是 Data.Traversable的 sequenceA，要求 f为 Applicative。因此，我们已经从 Applicative获得了重复。

sequenceA . take 6 . repeat $ sample2d

要将请求分组在一起，我们将在 mark中添加一个函数，该分组有意义。

sequenceA . take 6 . repeat . mark $ sample2d

还有一个可以标记某些分组的类。如果我们不仅需要分组，还需要更多的含义-例如，如果内部事物是依赖的或独立的，则可以在此处添加。

class Mark f where
    mark :: f a -> f a

如果一切都非常相似，我们可以为可查询的示例源添加一个类

class (Applicative f, Mark f, SampleSource f) => QueryableSampleSouce f where

现在，我们将讨论具有更优化查询语言的monad的想法。在这里，我们将开始使用所有这些GHC特定扩展。特别是 TypeFamilies。

class MonadQuery m where
    type Query m :: * -> *
    interpret :: (Query m a) -> m a

最后是使用 Applicative查询语言的monad示例源的类

class (MonadQuery m, QueryableSampleSource (Query m), SampleSource m, Monad m) => MonadSample m where

在这一点上，我们将要弄清楚应该遵循哪些法律。我建议一些：

interpret sample == sample
interpret (sequenceA a) = sequence (interpret a)

也就是说，没有 mark，样本源就无法对查询做任何特别的事情。这意味着要对泊松磁盘对2d点进行特殊处理和对点集进行特殊处理的查询需要标记两次：

 mark . sequenceA . take 6 . repeat . mark $ sample2d

Applicative查询语言的排序方式与您的 StratGen类型相对应；通过拥有一个很小的 Applicative接口，它可以使您提前查看传入查询的结构。然后 Monad与您的 StratRun类型相对应。