Блог

Я хочу написать код, который вычисляет объем n-го шара.
С помощью последовательности случайных точек r_i(i = 1, … , N) можно приблизительно определить объем такого шара:
где d — размер, а N — количество точек

У меня возникли проблемы с запуском этого кода

 function MonteCarloHypersphereVolume(radius, dimension, number_of_generations)
    
number_within_sphere = 0;
  for i = 1 : number_of_generations
    randoms = zeros( 1, dimension );
    for j = 1 : dimension 
        randoms(j) = rand(radius * 2) - radius;
    end

    if sum( randoms .^ 2 ) <= radius^2
        number_within_sphere = number_within_sphere   1;
    end
end

approximate_volume = (number_within_sphere/number_of_generations) * (2*radius)^dimension;

end

end
 

Я не уверен, что происходит не так

Итак, есть несколько ошибок, которые бросаются в глаза. У вас на одного больше end , чем вам нужно, и очередь

 randoms(j) = rand(radius * 2) - radius;
 

имеет пару очевидных ошибок, потому что (1) вы индексируете массивы с [] помощью in julia , а не () и (2) rand(n) генерирует n случайные числа, а не случайные числа между 0 и n . Так что, если мы исправим эти проблемы

 function MonteCarloHypersphereVolume(radius, dimension, number_of_generations)
    
    number_within_sphere = 0
    for i = 1 : number_of_generations
        randoms = zeros( 1, dimension )
        for j = 1 : dimension 
            randoms[j] = 2*radius*rand() - radius
        end

        if sum( randoms .^ 2 ) <= radius^2
            number_within_sphere = number_within_sphere   1
        end
    end

    approximate_volume = (number_within_sphere/number_of_generations) * (2*radius)^dimension;

end
 

Теперь у нас есть подход, который, кажется, дает нам правильные цифры, если мы протестируем простой 3d-случай

 julia> MonteCarloHypersphereVolume(1,3,100000)
4.18328

julia> 4/3*pi*1
4.1887902047863905
 

Теперь мы можем провести несколько оптимизаций. В настоящее время это занимает около 13 мс и огромное количество выделений

 julia> using BenchmarkTools

julia> @benchmark MonteCarloHypersphereVolume(1,3,100000)
BechmarkTools.Trial: 375 samples with 1 evaluations.
 Range (min … max):  11.016 ms … 23.764 ms  ┊ GC (min … max): 0.00% … 12.36%
 Time  (median):     12.893 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   13.339 ms ±  1.647 ms  ┊ GC (mean ± σ):  4.24% ±  4.51%

      ▄▃▂█  ▆▇▅▇▁▁ ▂▅ ▁                                        
  ▃▄█▄█████▇███████████▆▇▆▅▇▅█▅▅▄▃▆▃▄▄▄▃▁▃▃▃▁▁▃▄▃▁▃▁▁▁▃▁▁▃▁▃▃ ▄
  11 ms           Histogram: frequency by time        18.9 ms <

 Memory estimate: 21.36 MiB, allocs estimate: 200000.
 

Мы можем добиться гораздо большего успеха, выйдя randoms = zeros( 1, dimension ) за пределы цикла и используя sum функцию таким образом, чтобы ей не нужно было выделять промежуточный массив для квадратов. Мы также можем @inbounds немного ускорить for циклы, учитывая, что здесь не следует ожидать каких-либо выходящих за рамки индексов.

 function montecarlohyperspherevolume(radius, dimension, number_of_generations)
    
    number_within_sphere = 0
    randoms = zeros( 1, dimension )
    @inbounds for i = 1 : number_of_generations
        for j = 1 : dimension 
            randoms[j] = 2*radius*rand() - radius
        end

        if sum(abs2, randoms) <= radius^2 # could also use anonymous function x->x^2 if you didn't wnat to use abs2
            number_within_sphere  = 1
        end
    end

    return approximate_volume = (number_within_sphere/number_of_generations) * (2*radius)^dimension

end
 

Я также внес несколько стилистических изменений. Помимо прочего, в Джулии функции обычно строчные по соглашению (в то время как типы могут быть в верхнем регистре).

В любом случае, это дает нам примерно 6-кратное ускорение и значительно сокращает использование и распределение памяти

 julia> @benchmark montecarlohyperspherevolume(1,3,100000)
BechmarkTools.Trial: 2249 samples with 1 evaluations.
 Range (min … max):  1.846 ms …   4.053 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.080 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.217 ms ± 377.123 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   █▆▅▂▁▂                                                      
  ▅████████▇▆▅▅▅▃▃▃▃▃▃▄▃▃▂▃▂▃▂▂▃▂▂▂▂▂▁▁▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  1.85 ms         Histogram: frequency by time        3.61 ms <

 Memory estimate: 112 bytes, allocs estimate: 1.