Haskell: Эффективная итерационная карта с добавленным шумом
Мне интересно, как улучшить временные характеристики добавления белого шума к логистической карте? Шум допускается добавлять только после расчета значений (так как это итерационная карта).
module Generating
import System.Random (Random,randomR,random,mkStdGen,StdGen)
import Data.Random.Normal (mkNormals,normal)
import qualified Data.Vector.Unboxed as U
import Control.Monad.State
genR :: (Random a, Fractional a) => Int -> (a, StdGen)
genR x = randomR (0,1.0) (mkStdGen x)
new ::Double-> Double ->Int -> (Double,Int) -> U.Vector (Double,Int)
new skal r n = U.unfoldrN n go
where
go (x0,g0) = let !eins= (1.0-x0)
!x=x0 `seq` eins `seq` r*x0*eins
!g=g0+1
!noise= skal*(fst $ genR g)
in Just ((x0+noise,g0),(x,g))
fs :: (t, t1) -> t
fs (x,y)=x
first :: U.Vector (Double,Int)->U.Vector Double
first =U.map (\(x,y)->x)
Как видите, я на самом деле хочу только первое значение кортежа, но генератор должен быть обновлен.
Какие-либо предложения? Может быть, государственные монады?
1 ответ
tl;dr: не пытайтесь оптимизировать программу на Haskell без использования профилирования и бенчмаркинга. Добавление случайных восклицательных знаков и seqС почти никогда не будет работать. Большая проблема здесь в том, что StdGen это невероятно медленный генератор случайных чисел, который полностью доминирует во времени выполнения вашей программы. Вы должны заменить его, чтобы добиться значительного прогресса.
Вот более длинный ответ: хорошим первым шагом является установка библиотеки для бенчмаркинга, например criterionи напишите контрольный пример:
import Criterion.Main
...your program above...
vect1 :: (Double, Int) -> U.Vector Double
vect1 = first . new 0.5 1 10000
main = defaultMain [
bench "vect1" $ nf vect1 (0,1)
]
В моем случае результаты выглядят так:
benchmarking vect1
time 8.097 ms (8.071 ms .. 8.125 ms)
1.000 R² (1.000 R² .. 1.000 R²)
mean 8.140 ms (8.124 ms .. 8.162 ms)
std dev 52.90 μs (36.32 μs .. 91.72 μs)
поэтому у нас есть около 8 миллисекунд на цикл для генерации вектора из 10000 элементов.
Теперь давайте избавимся от всех ударов, seqs и промежуточные вычисления, которые вы добавили, чтобы попытаться ускорить процесс:
new :: Double-> Double -> Int -> (Double, Int) -> U.Vector (Double,Int)
new skal r n = U.unfoldrN n go
where
go (x0,g0) = let x = r * x0 * (1-x0)
g = g0 + 1
noise = skal * (fst $ genR g)
in Just ((x0+noise, g0), (x,g))
Повторяю, вот результаты:
time 8.195 ms (8.168 ms .. 8.235 ms)
Ах, они вообще не имели никакого эффекта. Рад, что мы от них избавились.
Теперь стоит отметить, что unfoldrN несет с собой аккумулятор для вас, который держится за ваш g, Вам также не нужно включать g в результате, если вы собираетесь все равно выбросить, мы можем упростить new чтобы:
new :: Double-> Double -> Int -> (Double, Int) -> U.Vector Double
new skal r n = U.unfoldrN n go
where
go (x0,g0) = let x = r * x0 * (1-x0)
g = g0 + 1
noise = skal * (fst $ genR g)
in Just (x0+noise, (x,g))
и бросить first вызов из определения vect1:
vect1 :: (Double, Int) -> U.Vector Double
vect1 = new 0.5 1 10000
Это дает:
time 8.289 ms (8.238 ms .. 8.373 ms)
так что это не имело значения. Несомненно, компилятор смог оптимизировать ненужные лишние Doubleв любом случае, поэтому изменение кода не дало никаких результатов.
Более серьезная проблема с алгоритмом состоит в том, что он использует генераторы очень странным образом. StdGen предназначен для посева и последующего повторного использования для генерации нескольких случайных чисел, а не для генерации свежей из начального числа на основе счетчика. Мы действительно должны переписать new правильно использовать генератор:
new :: Double-> Double -> Int -> (Double, Int) -> U.Vector Double
new skal r n (x0, seed) = U.unfoldrN n go (x0, g0)
where
g0 = mkStdGen seed -- create initial generator from seed
go (x0,g0) = let (eps, g) = randomR (0, 1.0) g0 -- use generator properly
x = r * x0 * (1-x0)
noise = skal * eps
in Just (x0 + noise, (x, g))
хотя опять же, это практически не имеет значения для нашего времени тестирования. Я признаю, что это удивило меня. Я думал, что это будет иметь значительный эффект. Хорошо, что я сравнивал эти изменения, поэтому у меня были реальные объективные доказательства эффекта (или отсутствия эффекта) этого изменения!
Теперь, похоже, пришло время описать нашу программу и посмотреть, на что она тратит свое время.
$ stack ghc -- -prof -fprof-auto -O2 Generating.hs
$ ./Generating -n 100 +RTS -p # run 100 iterations
Если вы посмотрите на Generating.prof файл, который выводится, вы увидите, что большая часть времени проводится в System.Random, вот так:
COST CENTRE MODULE SRC %time %alloc
randomR System.Random System/Random.hs:409:3-27 21.7 24.0
stdNext System.Random System/Random.hs:(518,1)-(528,64) 15.4 16.6
randomIvalInteger System.Random System/Random.hs:(468,1)-(489,76) 12.2 12.0
randomIvalInteger.f System.Random System/Random.hs:(486,8)-(489,76) 11.0 4.8
randomIvalInteger.f.v' System.Random System/Random.hs:489:25-76 7.0 8.6
Оказывается, что стандартный генератор случайных чисел на Haskell работает ужасно медленно, и нам нужно заменить его чем-то более быстрым, чтобы добиться большего прогресса.
mersenne-random-pure64 Пакет обеспечивает быструю реализацию Mersenne Twister, которая производит случайные числа высокого качества, и мы можем переписать new использовать это. Обратите внимание, что randomDouble возвращает равномерное случайное число в интервале [0,1):
import System.Random.Mersenne.Pure64
new :: Double-> Double -> Int -> (Double, Int) -> U.Vector Double
new skal r n (x0, seed) = U.unfoldrN n go (x0, g0)
where
g0 = pureMT (fromIntegral seed)
go (x0,g0) = let (eps, g) = randomDouble g0
x = r * x0 * (1-x0)
noise = skal * eps
in Just (x0 + noise, (x, g))
Повторный бенчмаркинг (перекомпилированный без профилирования) дает:
benchmarking vect1
time 106.7 μs (106.4 μs .. 107.0 μs)
1.000 R² (1.000 R² .. 1.000 R²)
mean 107.1 μs (106.7 μs .. 107.7 μs)
std dev 1.415 μs (842.3 ns .. 2.377 μs)
Обратите внимание, что это 107 микросекунд, так что это примерно в 75 раз быстрее.
На этом я бы остановился, но если вы решите продолжить оптимизацию, убедитесь, что вы часто обращаетесь к результатам профилирования и бенчмаркинга, чтобы убедиться, что ваши изменения дают о себе знать.
Я настоятельно рекомендую Googling для "профилирования программ на haskell" и для "критериальной" библиотеки, а также уделю некоторое время изучению того, как использовать эти инструменты.
Для справки, окончательная программа:
import Criterion.Main
import qualified Data.Vector.Unboxed as U
import System.Random.Mersenne.Pure64
new :: Double-> Double -> Int -> (Double, Int) -> U.Vector Double
new skal r n (x0, seed) = U.unfoldrN n go (x0, g0)
where
g0 = pureMT (fromIntegral seed)
go (x0,g0) = let (eps, g) = randomDouble g0
x = r * x0 * (1-x0)
noise = skal * eps
in Just (x0 + noise, (x, g))
vect1 :: (Double, Int) -> U.Vector Double
vect1 = new 0.5 1 10000
main = defaultMain [
bench "vect1" $ nf vect1 (0,1)
]