Computing the double Integral using MonteCarlo techniques using Julia

Question

I decided to try and learn Julia for doing scientific computing, and I decided to tackle the problem of finding

$$ \int_{D_{\frac{1}{4}}} x^4 + y^2 dA $$

where \$ D_{\frac{1}{4}} \$ is the part of the unit circle in the first cuadrant.

My code in Julia is the following:

using Distributions
e = 10.0^(-3);
p = 0.85;
variance = 4;

N = floor(Int, variance / ((1-p)*((e/2)^2))) + 1

u = Uniform(0,2);
x = rand(N);
y = rand(N);
z = rand(u, N);

result = sum((x.^2 + y.^2 .<= 1) & (z .<= x.^4 + y.^2))*2.0 / N

which gives the nice result \$ = 0.2945746303294543 \$

I kindly ask for how to improve my implementation, and reduce the footprint of memory (it uses almost 2 to 3gb in RAM).

Answer 1

Try swapping out the three arrays with calling a function iteratively---three Int64 arrays (each element 8 bytes) of ~100M elements is a lot. That is, switch out the body of the sum into a separate function:

function foo(u)
    x = rand()
    y = rand()
    z = rand(u)
    (x^2 + y^2 <= 1) & (z <= x^4 + y^2) ? 1 : 0
end

and replace the sum expression itself with just a loop:

count = 0
for i in 1:N
    count += foo(u)
end
count * 2.0 / N

This runs a bit faster, since the compiler can optimize the internal function, and with a much smaller memory footprint, since we sidestep the arrays entirely.

Answer 2

The following is internally equivalent to @nybblet's answer, but uses nicer syntax instead of a manual loop:

f((x,y,z)) = (x^2 + y^2 <= 1) & (z <= x^4 + y^2)

function computeintegral(e, p, variance)
    N = floor(Int, variance / ((1-p)*((e/2)^2))) + 1
    x = (rand() for _ in 1:N)
    y = (rand() for _ in 1:N)
    z = (2rand() for _ in 1:N)

    sum(f, zip(x, y, z)) * 2.0 / N
end

The trick is to use generators instead of arrays for the random values. Since you are only iterating over them once, consuming one element at a time, no extra space is used:

julia> @btime sum(f, zip(x, y, z)) * 2.0 / N
  11.289 s (11 allocations: 400 bytes)
0.29453870532956655

Of course this depends on sum(f, itr) and zip using iterators in the right (lazy) way.

f makes use of argument destructuring, introduced (I think?) in 0.7.

And I replaced the usage of a Distribution object in z by a simple transformation, but that was just to save adding a package. For more complex domains, it probably makes things more readable.