我目前正在使用Haskell制作光线跟踪程序.由于我是Haskell的初学者,我不明白IO monad的评估策略.
问题是一长串"IO a"的内存使用情况,在我的代码中是"IO Vec".
列表的每个元素由递归函数计算IO Vec,该函数计算表示像素的颜色.因此,列表的长度等于width x height.
另外,我为像素拍摄了多个样本.总的来说,radiance计算像素值的函数称为width x height x samples时间.
首先,我只是通过使用列表理解来实现这个程序.代码就像,
main = do
...
let ray = (compute ray for every pair of [0..w-1], [0..h-1]
pixels <- sequence [ (sumOfRadiance scene ray samples) | ray <- rays]
在我的理解中,由于在将像素写入文件之前未使用像素,因此Haskell存储一些用于函数调用的数据,pixels其中是一个数组IO Vec.最后,通过调用递归函数radiance来计算像素值,可以增加内存消耗.
如果我改变程序来逐个评估像素值,unsafePerformIO可以防止这种奇怪的使用内存空间.
main = do
...
let ray = (compute ray for every pair of [0..w-1], [0..h-1]
let pixels = [ (unsafePerformIO (sumOfRadiance scene ray samples)) | ray <- rays]
我知道unsafePerformIO是一个糟糕的解决方案,所以我想知道Haskell是否提供了另一种立即评估IO monad内部的方法.以下是我的整个代码(抱歉,它有点长......)
谢谢您的帮助.
-- Small path tracing with Haskell
import System.Environment
import System.Random.Mersenne
import System.IO.Unsafe
import Control.Monad
import Codec.Picture
import Data.Time
import qualified Data.Word as W
import qualified Data.Vector.Storable as V
-- Parameters
eps :: Double
eps = 1.0e-4
inf :: Double
inf = 1.0e20
nc :: Double
nc = 1.0
nt :: Double
nt = 1.5
-- Vec
data Vec = Vec (Double, Double, Double) deriving (Show)
instance (Num Vec) where
(Vec (x, y, z)) + (Vec (a, b, c)) = Vec (x + a, y + b, z + c)
(Vec (x, y, z)) - (Vec (a, b, c)) = Vec (x - a, y - b, z - c)
(Vec (x, y, z)) * (Vec (a, b, c)) = Vec (x * a, y * b, z * c)
abs = undefined
signum = undefined
fromInteger x = Vec (dx, dx, dx) where dx = fromIntegral x
x :: Vec -> Double
x (Vec (x, _, _)) = x
y :: Vec -> Double
y (Vec (_, y, _)) = y
z :: Vec -> Double
z (Vec (_, _, z)) = z
mul :: Vec -> Double -> Vec
mul (Vec (x, y, z)) s = Vec (x * s, y * s, z * s)
dot :: Vec -> Vec -> Double
dot (Vec (x, y, z)) (Vec (a, b, c)) = x * a + y * b + z * c
norm :: Vec -> Vec
norm (Vec (x, y, z)) = Vec (x * invnrm, y * invnrm, z * invnrm)
where invnrm = 1 / sqrt (x * x + y * y + z * z)
cross :: Vec -> Vec -> Vec
cross (Vec (x, y, z)) (Vec (a, b, c)) = Vec (y * c - b * z, z * a - c * x, x * b - a * y)
-- Ray
data Ray = Ray (Vec, Vec) deriving (Show)
org :: Ray -> Vec
org (Ray (org, _)) = org
dir :: Ray -> Vec
dir (Ray (_, dir)) = dir
-- Material
data Refl = Diff
| Spec
| Refr
deriving Show
-- Sphere
data Sphere = Sphere (Double, Vec, Vec, Vec, Refl) deriving (Show)
rad :: Sphere -> Double
rad (Sphere (rad, _, _, _, _ )) = rad
pos :: Sphere -> Vec
pos (Sphere (_ , p, _, _, _ )) = p
emit :: Sphere -> Vec
emit (Sphere (_ , _, e, _, _ )) = e
col :: Sphere -> Vec
col (Sphere (_ , _, _, c, _ )) = c
refl :: Sphere -> Refl
refl (Sphere (_ , _, _, _, refl)) = refl
intersect :: Sphere -> Ray -> Double
intersect sp ray =
let op = (pos sp) - (org ray)
b = op `dot` (dir ray)
det = b * b - (op `dot` op) + ((rad sp) ** 2)
in
if det < 0.0
then inf
else
let sqdet = sqrt det
t1 = b - sqdet
t2 = b + sqdet
in ansCheck t1 t2
where ansCheck t1 t2
| t1 > eps = t1
| t2 > eps = t2
| otherwise = inf
-- Scene
type Scene = [Sphere]
sph :: Scene
sph = [ Sphere (1e5, Vec ( 1e5+1, 40.8, 81.6), Vec (0.0, 0.0, 0.0), Vec (0.75, 0.25, 0.25), Diff) -- Left
, Sphere (1e5, Vec (-1e5+99, 40.8, 81.6), Vec (0.0, 0.0, 0.0), Vec (0.25, 0.25, 0.75), Diff) -- Right
, Sphere (1e5, Vec (50.0, 40.8, 1e5), Vec (0.0, 0.0, 0.0), Vec (0.75, 0.75, 0.75), Diff) -- Back
, Sphere (1e5, Vec (50.0, 40.8, -1e5+170), Vec (0.0, 0.0, 0.0), Vec (0.0, 0.0, 0.0), Diff) -- Front
, Sphere (1e5, Vec (50, 1e5, 81.6), Vec (0.0, 0.0, 0.0), Vec (0.75, 0.75, 0.75), Diff) -- Bottom
, Sphere (1e5, Vec (50,-1e5+81.6,81.6), Vec (0.0, 0.0, 0.0), Vec (0.75, 0.75, 0.75), Diff) -- Top
, Sphere (16.5, Vec (27, 16.5, 47), Vec (0.0, 0.0, 0.0), Vec (1,1,1) `mul` 0.999, Spec) -- Mirror
, Sphere (16.5, Vec (73, 16.5, 78), Vec (0.0, 0.0, 0.0), Vec (1,1,1) `mul` 0.999, Refr) -- Glass
, Sphere (600, Vec (50, 681.6 - 0.27, 81.6), Vec (12, 12, 12), Vec (0, 0, 0), Diff) ] -- Light
-- Utility functions
clamp :: Double -> Double
clamp = (max 0.0) . (min 1.0)
isectWithScene :: Scene -> Ray -> (Double, Int)
isectWithScene scene ray = foldr1 (min) $ zip [ intersect sph ray | sph <- scene ] [0..]
nextDouble :: IO Double
nextDouble = randomIO
lambert :: Vec -> Double -> Double -> (Vec, Double)
lambert n r1 r2 =
let th = 2.0 * pi * r1
r2s = sqrt r2
w = n
u = norm $ (if (abs (x w)) > eps then Vec (0, 1, 0) else Vec (1, 0, 0)) `cross` w
v = w `cross` u
uu = u `mul` ((cos th) * r2s)
vv = v `mul` ((sin th) * r2s)
ww = w `mul` (sqrt (1.0 - r2))
rdir = norm (uu + vv + ww)
in (rdir, 1)
reflect :: Vec -> Vec -> (Vec, Double)
reflect v n =
let rdir = v - (n `mul` (2.0 * n `dot` v))
in (rdir, 1)
refract :: Vec -> Vec -> Vec -> Double -> (Vec, Double)
refract v n orn rr =
let (rdir, _) = reflect v orn
into = (n `dot` orn) > 0
nnt = if into then (nc / nt) else (nt / nc)
ddn = v `dot` orn
cos2t = 1.0 - nnt * nnt * (1.0 - ddn * ddn)
in
if cos2t < 0.0
then (rdir, 1.0)
else
let tdir = norm $ ((v `mul` nnt) -) $ n `mul` ((if into then 1 else -1) * (ddn * nnt + (sqrt cos2t)))
a = nt - nc
b = nt + nc
r0 = (a * a) / (b * b)
c = 1.0 - (if into then -ddn else (tdir `dot` n))
re = r0 + (1 - r0) * (c ** 5)
tr = 1.0 - re
pp = 0.25 + 0.5 * re
in
if rr < pp
then (rdir, (pp / re))
else (tdir, ((1.0 - pp) / tr))
radiance :: Scene -> Ray -> Int -> IO Vec
radiance scene ray depth = do
let (t, i) = (isectWithScene scene ray)
if inf <= t
then return (Vec (0, 0, 0))
else do
r0 <- nextDouble
r1 <- nextDouble
r2 <- nextDouble
let obj = (scene !! i)
let c = col obj
let prob = (max (x c) (max (y c) (z c)))
if depth >= 5 && r0 >= prob
then return (emit obj)
else do
let rlt = if depth < 5 then 1 else prob
let f = (col obj)
let d = (dir ray)
let x = (org ray) + (d `mul` t)
let n = norm $ x - (pos obj)
let orn = if (d `dot` n) < 0.0 then n else (-n)
let (ndir, pdf) = case (refl obj) of
Diff -> (lambert orn r1 r2)
Spec -> (reflect d orn)
Refr -> (refract d n orn r1)
nextRad <- (radiance scene (Ray (x, ndir)) (succ depth))
return $ ((emit obj) + ((f * nextRad) `mul` (1.0 / (rlt * pdf))))
toByte :: Double -> W.Word8
toByte x = truncate (((clamp x) ** (1.0 / 2.2)) * 255.0) :: W.Word8
accumulateRadiance :: Scene -> Ray -> Int -> Int -> IO Vec
accumulateRadiance scene ray d m = do
let rays = take m $ repeat ray
pixels <- sequence [radiance scene r 0 | r <- rays]
return $ (foldr1 (+) pixels) `mul` (1 / fromIntegral m)
main :: IO ()
main = do
args <- getArgs
let argc = length args
let w = if argc >= 1 then (read (args !! 0)) else 400 :: Int
let h = if argc >= 2 then (read (args !! 1)) else 300 :: Int
let spp = if argc >= 3 then (read (args !! 2)) else 4 :: Int
startTime <- getCurrentTime
putStrLn "-- Smallpt.hs --"
putStrLn $ " width = " ++ (show w)
putStrLn $ " height = " ++ (show h)
putStrLn $ " spp = " ++ (show spp)
let dw = fromIntegral w :: Double
let dh = fromIntegral h :: Double
let cam = Ray (Vec (50, 52, 295.6), (norm $ Vec (0, -0.042612, -1)));
let cx = Vec (dw * 0.5135 / dh, 0.0, 0.0)
let cy = (norm $ cx `cross` (dir cam)) `mul` 0.5135
let dirs = [ norm $ (dir cam) + (cy `mul` (y / dh - 0.5)) + (cx `mul` (x / dw - 0.5)) | y <- [dh-1,dh-2..0], x <- [0..dw-1] ]
let rays = [ Ray ((org cam) + (d `mul` 140.0), (norm d)) | d <- dirs ]
let pixels = [ (unsafePerformIO (accumulateRadiance sph r 0 spp)) | r <- rays ]
let pixelData = map toByte $! pixels `seq` (foldr (\col lst -> [(x col), (y col), (z col)] ++ lst) [] pixels)
let pixelBytes = V.fromList pixelData :: V.Vector W.Word8
let img = Image { imageHeight = h, imageWidth = w, imageData = pixelBytes } :: Image PixelRGB8
writePng "image.png" img
endTime <- getCurrentTime
print $ diffUTCTime endTime startTime
luq*_*qui 10
首先,我认为有一个错误.当你谈到来自
pixels <- sequence [ (sumOfRadiance scene ray samples) | ray <- rays]
至
pixels <- sequence [ (unsafePerformIO (sumOfRadiance scene ray samples)) | ray <- rays]
这没有意义.类型不应该匹配 - sequence只有在组合一堆类型的东西时才有意义m a.这样做是正确的
let pixels = [ unsafePerformIO (sumOfRadiance scene ray samples) | ray <- rays ]
我会有点骑士认为这就是你所做的,而你在输入问题时却犯了一个错误.
如果是这种情况,那么您实际需要的是一种更懒惰而不是立即执行IO操作的方法.该调用会强制所有操作立即运行,而版本只是创建一个未运行的操作列表(实际上列表本身是懒惰生成的,因此它不会同时存在),并且操作单独运行因为他们的结果是必要的.sequenceunsafePerformIO
看来您需要的原因IO是生成随机数.随机性可能是一种痛苦 - 通常MonadRandom是工作,但它仍然会产生行为之间的顺序依赖,可能仍然不够懒惰(我试试看 - 如果你使用它,你会得到再现性 -即使在遵守monad法则的重构之后,相同的种子也会产生相同的结果.
如果MonadRandom不起作用并且您需要以更加按需的方式生成随机数,那么就可以创建自己的随机性monad,它与您的unsafePerformIO解决方案做同样的事情,但是以适当的方式封装.我打算告诉你我认为是Haskell Way To Cheat的方式.首先,一个可爱的纯实现草图:
-- A seed tells you how to generate random numbers
data Seed = ...
splitSeed :: Seed -> (Seed, Seed)
random :: Seed -> Double
-- A Cloud is a probability distribution of a's, or an a which
-- depends on a random seed. This monad is just as lazy as a
-- pure computation.
newtype Cloud a = Cloud { runCloud :: Seed -> a }
deriving (Functor)
instance Monad Cloud where
return = Cloud . const
m >>= f = Cloud $ \seed ->
let (seed1, seed2) = splitSeed seed in
runCloud (f (runCloud m seed1)) seed2
(我认为我做对了.重点是,在每次绑定时,你将种子分成两部分,然后将一部分传递到左边,另一部分传递到右边.)
现在这是一个完全纯粹的随机性实现......有几个捕获.(1)不存在非平凡的splitSeed,将严格尊重法律的单子,和(2)即使我们要被打破的规律,依据分裂随机数生成器可以是相当缓慢的.但是,如果我们放弃决定论,如果我们关心的是我们从分布中得到一个很好的抽样而不是完全相同的结果,那么我们就不需要严格遵守monad定律.在那一点上,我们作弊并假装有一个合适的Seed类型:
data Seed = Seed
splitSeed Seed = (Seed, Seed)
-- Always NOINLINE functions with unsafePerformIO to keep the
-- optimizer from messing with you.
{-# NOINLINE random #-}
random Seed = unsafePerformIO randomIO
我们应该在模块中隐藏它以保持抽象屏障. Cloud并且runCloud不应该被暴露,因为它们允许我们违反纯度; 只暴露
runCloudIO :: Cloud a -> IO a
runCloudIO = return . runCloud
这在技术上并不需要IO,但传达这不是确定性的.然后你可以在Cloudmonad中构建你需要的任何值,并在主程序中运行一次.
Seed如果它没有任何信息,您可能会问我们为什么会有类型.好吧,我认为splitSeed这只是对纯度的点头,实际上并没有做任何事情 - 你可以删除它 - 但我们需要Cloud成为一个函数类型,以便隐含的懒惰缓存不会破坏我们的语义.除此以外
let foo = random in liftM2 (,) foo foo
总会返回一对有两个相同组件的对,因为随机值确实与之相关foo.我不确定这些事情,因为此时我们正在与优化器进行战争,需要进行一些实验.
快乐的作弊.:-)