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).


2 Answers 2


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

and replace the sum expression itself with just a loop:

count = 0
for i in 1:N
    count += foo(u)
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.


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

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)

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.