3

在随机光线追踪器的上下文中,我想将 MC 集成(路径追踪、双向路径追踪)与样本生成(均匀随机、分层、泊松、大都会......)分离。其中大部分已经实现,但是使用起来很乏味。所以我放弃了它并尝试通过将采样计算分成两个阶段来构建更好的东西:SampleGen你可以使用mk1dand函数请求一个随机值,然后由采样算法mk2d提供实际的 s 。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中只需要一个RNGSampleRun
  • 分层和泊松圆盘采样器(可能还有其他)都跟踪所需的一维和二维样本的数量并将它们预先计算为一个向量,并且它们将被允许共享一个SampleGenSampleRun实现,只是在两者之间发生的事情有所不同SampleGenSampleRun(如何实际填充向量)
  • Metropolis 采样器将在其阶段使用惰性样本生成技术SampleRun

我想使用 GHC 来编译它,所以扩展喜欢MultiParamTypeClasses并且TypeFamilies对我来说没问题,但我没有想出任何远程可用的东西。

PS:作为动机,一些漂亮的图片。当前形式的代码在GitHub 上

4

1 回答 1

3

我将从一个完全不同的问题开始,“代码应该是什么样子”?然后转向“采样框架是如何组合在一起的”这个问题?

代码应该是什么样子

的定义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 是正方形的面积,如果它们落在一个圆内,否则为 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)

这与现有的 MonadRandom 类非常相似,除了两件事。一个MonadRandom实现本质上需要Int在它自己选择的某个范围内提供一个均匀随机的。您的采样器将Float在区间 [0..1) 上提供未知分布的样本。这足以证明拥有自己的新课程是合理的。

由于即将发生的Monad Applicative变化,我将为这个类建议一个不同的名称,SampleSource.

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

samplemk1d在您的代码中替换。mk2d也可以更换,同样不知道样品的来源是什么。sample2d, 的替换mk2d, 将适用于任何Applicative示例源,它不需要是Monad. 它不需要 a 的原因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

当所有Monads 都是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

我们不会完全实现Stratifiedor runStratified,因为您的示例代码不包含它们的完整实现。

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

我不确定你到底想用“分层”抽样做什么。预先生成数字,并在这些数字用完时使用生成器并不是我所理解的分层抽样。如果您要为某些东西提供一元接口,您将无法提前知道将执行什么,因此您将无法在开始执行之前预测计算需要多少样本。如果您只能满足于一个Applicative接口,那么您可以提前测试整个计算需要多少样本。

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

如果单次采样可以同时依赖于所需样本的数量和维度的数量,例如在 Poisson Disk 采样中,则需要在已知时将它们传递给采样器。

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]。这已经存在;it's Data.TraversablesequenceA这要求fbe 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查询语言的单子样本源类

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类型相对应。

于 2014-07-16T23:33:19.120 回答