# Generated (staged) function to avoid initializing an array multiple times for Monte Carlo simulations

**URL:** <https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726>\
**Category:** General Usage\
**Created:** [May 15, 2017, 7:14pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726 "2017-05-15T19:14:07Z")\
**Posts on this page:** 3\
**Page:** 2

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 19, 2017, 1:51pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/21 "2017-05-19T13:51:40Z")

</div>

> [@TsurHerman](#):
>
> I am no expert on type system and type stability and maybe I use the wrong term here, but from experience closures created using generating functions using static variables in the form of Ref’s as in my example , compile to optimal machine code.

Usually. But this is not something you can easily guarantee or reason about because of this bug:

> <https://github.com/JuliaLang/julia/issues/15276>
>
> \`\`\` jl
> using Images: realtype
> 
> function ifi{T\<:Real,K,N}(img::AbstractArray{T,N}…, kern::AbstractArray{K,N}, border::AbstractString, value)
> if border == "circular" && size(img) == size(kern)
> out = real(ifftshift(ifft(fft(img).\*fft(kern))))
> elseif border != "inner"
> prepad = \[div(size(kern,i)-1, 2) for i = 1:N\]
> postpad = \[div(size(kern,i), 2) for i = 1:N\]
> fullpad = \[nextprod(\[2,3\], size(img,i) + prepad\[i\] + postpad\[i\]) - size(img, i) - prepad\[i\] for i = 1:N\]
> A = padarray(img, prepad, fullpad, border, convert(T, value))
> krn = zeros(typeof(one(T)\*one(K)), size(A))
> indexesK = ntuple(d-\>\[size(krn,d)-prepad\[d\]+1:size(krn,d);1:size(kern,d)-prepad\[d\]\], N)::NTuple{N,Vector{Int}}
> AF = ifft(fft(A).\*fft(krn))
> out = Array(realtype(eltype(AF)), size(img))
> end
> out
> end
> \`\`\`
> 
> Test:
> 
> \`\`\` jl
> julia\> @code\_warntype ifi(rand(3,3), rand(3,3), "replicate", 0)
> Variables:
> #self#::#ifi
> img::Array{Float64,2}
> kern::Array{Float64,2}
> border::ASCIIString
> value::Int64
> prepad::Box
> ...
> \`\`\`
> 
> Now comment out the \`indexesK = ...\` line (the output of which is not used at all). Suddenly \`prepad\` is inferred as \`Array{Int, 1}\`.

I didn’t check your exact case, but if @musm says there’s an inference issue and you’re thing “but that should infer?”, you’re probably hitting this.

---

<div class="post-metadata">

**Author:** ![TsurHerman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tsurherman/32/1234_2.png) [@TsurHerman](https://discourse.julialang.org/u/TsurHerman)\
**Post date:** [May 19, 2017, 10:00pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/22 "2017-05-19T22:00:21Z")

</div>

@vmh @ChrisRackauckas @musm

At first I did not see any difference in benchmark (for n=1000) between using `draw_samples!` with a globally pre-allocated vector  
and between a generated `draw_samples` using

> [@TsurHerman](#):
>
> ```julia
> function fsamples_prealoc(g,n,T = Float64)
> samples = Ref{Array{T, 1}}( Array{T, 1}(n) )
> (args...) -&gt; g(samples.x,args...)
> end
> 
> ```

But then I benchmarked it with n=1 and I saw some overhead(~300 ns).After some investigation it turns out that  
it is mostly due to ellipsis operator.  
It also turns out that the Ref trick is unnecessary, I guess because a Vector is already an indirection.  
I used the Ref trick to bypass the built-in boxing operation that produced less optimal code in the past.

After this fix, the code\_native for the closure and for the original method is almost identical.

however I do see in code\_warntype that the return value is ANY , which would suggest type instability, but  
all benchmarks I made show that the generated function is so very slightly **faster** than the original…

my code:

```julia

using Distributions
G = Gamma(1,1)
n=1

samples2 = Array{Float64, 1}(n)

function draw_samples!(samples,dist)
    rand!(dist, samples)
    samples
end

function fsamples_prealoc(g,n)
    samples = Array{Float64, 1}(n)
    (x) -> g(samples,x)
end

draw_samples = fsamples_prealoc(draw_samples!,n)

@code_native draw_samples(G)

@code_native draw_samples!(samples2,G)

function foo2(d,samples)
    x = draw_samples!(samples,d)
    2x
end

function foo(d)
    x = draw_samples(d)
    2x
end

@benchmark foo2(G,samples2)

@benchmark foo(G)

@code_warntype foo2(G,samples2)

@code_warntype foo(G)

```

@ChrisRackauckas @musm maybe you can explain whats going on here and whether it is type stable or not  
and whether it has any implications on performance

Thanks  
Tsur

---

<div class="post-metadata">

**Author:** ![TsurHerman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tsurherman/32/1234_2.png) [@TsurHerman](https://discourse.julialang.org/u/TsurHerman)\
**Post date:** [May 21, 2017, 2:05pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/23 "2017-05-21T14:05:27Z")

</div>

@musm  
After a too long investigation it turns out that it is not type stable  
because draw\_samples2 which I presume is the name you gave the generated function, is not const.

if you declare it using:

`const draw_samples2 = fsamples_prealoc(draw_samples!,n)`

then foo2 becomes type\_stable.

[Previous page](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726.md?page=1)