# 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:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![vmh](https://avatars.discourse-cdn.com/v4/letter/v/ecb155/32.png) [@vmh](https://discourse.julialang.org/u/vmh)\
**Post date:** [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/1 "2017-05-15T19:14:07Z")

</div>

Hi everyone,

I have a function that returns an array (vector) of samples from a given distribution. E.g., consider the function below:

```julia
function draw_samples(dist)
    samples = Array{Float64, 1}(1000)
    rand!(dist, samples)
    x
end

```

And I’m calling this function a million times inside a loop for independent Monte Carlo (MC) trials. Really would like to avoid creating a new array (`samples`) in every iteration of my MC loop. Of course, I can always preallocate `samples` in my main program and pass it as an argument to `draw_samples!` to mutate the array (please note the ! at the end). However, I don’t like that style because things get really entangled and that approach is not quite modular. So, I’m considering a generated function, where in its body, I create `samples` and then the return statement is simply a call to the mutating `_draw_samples!` function, with the value of `samples` passed to it. E.g.,

```julia
function _draw_samples!(samples, dist)
    rand!(dist, samples)
    samples
end

@generated function draw_samples_1_allocation(dist)
    samples = Array{Float64, 1}(1000)
    quote
        _draw_samples!($(samples), dist)
    end
end

```

I’d really appreciate your advice on whether this approach makes any sense, what could go wrong with this method, or if there are any alternatives that are recommended over this one. Thank you for your help.

---

<div class="post-metadata">

**Author:** ![Keno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/keno/32/285_2.png) [@Keno](https://discourse.julialang.org/u/Keno)\
**Post date:** [May 15, 2017, 7:33pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/2 "2017-05-15T19:33:35Z")

</div>

Generated functions are the wrong approach to this in general. Try this:

```julia
let samples = Array{Float64, 1000}, inited = false
global draw_samples
function draw_samples(dist)
        inited || rand!(dist, samples)
        inited = true
        samples
end
end

```

Obviously you can adjust this examples as you see fit (e.g. adjust the cache to depend on dist), etc.

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [May 15, 2017, 7:35pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/3 "2017-05-15T19:35:36Z")

</div>

> [@vmh](#):
>
> However, I don’t like that style because things get really entangled and that approach is not quite modular.

There’s nothing non-modular about it. You can just have a context object and each “module” will be interested in a small part of the context.

> [@vmh](#):
>
> what could go wrong with this method,

It’s not reentrant or thread safe.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [May 15, 2017, 7:43pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/4 "2017-05-15T19:43:04Z")

</div>

I don’t see the point with `inited` That will cache the results from the first call which is not what the original code seems to do?

---

<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 15, 2017, 8:07pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/5 "2017-05-15T20:07:07Z")

</div>

Instead of a function use a call overloaded type with the array as a field?

```julia
# Define the type
immutable MyFunction{T}
  samples::T
end

# Define a call
function (p::MyFunction)(dist)
  _draw_samples!(p.samples, dist)
end

# Build the type with the cache
draw_samples_1_allocation = MyFunction(Array{Float64, 1}(1000))

# Now use it like a function
dist = ...
draw_samples_1_allocation(dist)

```

---

<div class="post-metadata">

**Author:** ![vmh](https://avatars.discourse-cdn.com/v4/letter/v/ecb155/32.png) [@vmh](https://discourse.julialang.org/u/vmh)\
**Post date:** [May 15, 2017, 8:47pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/6 "2017-05-15T20:47:41Z")

</div>

Thank you for your response. Would you recommend the preallocation approach even when there are say 20 different functions, each of which needs a different preallocated array (and all of them are being called in one single function)? Currently, I’ve written my code using preallocation in the single function and then mutating the arrays in individual functions but preallocating 20 different arrays in the single function (main program) just didn’t look right to me and I was wondering if I was doing something wrong.

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [May 15, 2017, 8:55pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/7 "2017-05-15T20:55:55Z")

</div>

> [@vmh](#):
>
> Would you recommend the preallocation approach even when there are say 20 different functions, each of which needs a different preallocated array (and all of them are being called in one single function)?

Yes.

---

<div class="post-metadata">

**Author:** ![gasagna](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gasagna/32/1275_2.png) [@gasagna](https://discourse.julialang.org/u/gasagna)\
**Post date:** [May 15, 2017, 9:10pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/8 "2017-05-15T21:10:49Z")

</div>

What are the implications of this being non reentrant?

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [May 15, 2017, 9:26pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/9 "2017-05-15T21:26:34Z")

</div>

As the code gets more complicated you can easily get’s into trouble if it ever recurse back.

---

<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 15, 2017, 9:52pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/10 "2017-05-15T21:52:15Z")

</div>

you could maybe benefit from the functor approach.

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

function fsamples_prealoc(g,n,T = Float64)
    const samples = Array{T, 1}(n)
    (args...) -> g(samples,args...)
end

draw_samples = fsamples_prealoc(draw_samples!,1000)

```

since you have 20 different functions with different requirements for cache size

---

<div class="post-metadata">

**Author:** ![vmh](https://avatars.discourse-cdn.com/v4/letter/v/ecb155/32.png) [@vmh](https://discourse.julialang.org/u/vmh)\
**Post date:** [May 16, 2017, 12:56am UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/11 "2017-05-16T00:56:20Z")

</div>

Thank you for your response. I’d checked the closure method before but unfortunately, it is considerably slower than the direct mutating approach.

> [@TsurHerman](#):
>
> you could maybe benefit from the functor approach.

---

<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 16, 2017, 1:15am UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/12 "2017-05-16T01:15:31Z")

</div>

> [@vmh](#):
>
> I’d checked the closure method before but unfortunately, it is considerably slower than the direct mutating approach.

You mean @TsurHerman’s? Tha because the const isn’t actually doing anything:

> <https://github.com/JuliaLang/julia/issues/5148>
>
> I am wondering why code:
> 
> \`\`\`
> function f()
> const x = "aaa"
> x = 4…
> return x
> end
> println(f())
> \`\`\`
> 
> executes and prints 4 without any warning or error, while
> 
> \`\`\`
> const x = "aaa"
> x = 4
> println(x)
> \`\`\`
> 
> raises:
> \`ERROR: invalid redefinition of constant x\`
> 
> Similarly:
> 
> \`\`\`
> function f()
> const x = "aaa"
> x = "bbb"
> return x
> end
> println(f())
> \`\`\`
> 
> prints bbb without warning, while:
> 
> \`\`\`
> const x = "aaa"
> x = "bbb"
> println(x)
> \`\`\`
> 
> prints bbb, but with warning:
> \`Warning: redefining constant x\`
> 
> I think behavior of constants should be consistent in global and function scopes.
> 
> My Julia version:
> 
> \`\`\`
> Julia Version 0.2.0
> Commit 05c6461 (2013-11-16 23:44 UTC)
> Platform Info:
> System: Windows (x86\_64-w64-mingw32)
> WORD\_SIZE: 64
> BLAS: libopenblas (USE64BITINT DYNAMIC\_ARCH NO\_AFFINITY)
> LAPACK: libopenblas
> LIBM: libopenlibm
> \`\`\`

But a call-overloaded immutable type should be fine.

---

<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 16, 2017, 7:15am UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/13 "2017-05-16T07:15:16Z")

</div>

you can write the closure using Ref’s and then it will be type stable and the compiler will optimise  
it correctly

```julia
function fsamples_prealoc(g,n,T = Float64)
    const samples = Array{T, 1}(n)
    (args...) -> g(samples,args...)
end

```

becomes

```julia
function fsamples_prealoc(g,n,T = Float64)
    samples = Ref{Array{T, 1}}( Array{T, 1}(n) )
    (args...) -> g(samples.x,args...)
end

```

---

<div class="post-metadata">

**Author:** ![vmh](https://avatars.discourse-cdn.com/v4/letter/v/ecb155/32.png) [@vmh](https://discourse.julialang.org/u/vmh)\
**Post date:** [May 18, 2017, 4:23pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/14 "2017-05-18T16:23:51Z")

</div>

@TsurHerman, thank you so much for your response. This is a great and neat suggestion. As someone who knows very little about thread-safe programming, I wonder if your suggestions is thread-safe. Any thoughts on this would be appreciated.

> [@TsurHerman](#):
>
> you can write the closure using Ref’s and then it will be type stable and the compiler will optimise it correctly

---

<div class="post-metadata">

**Author:** ![vmh](https://avatars.discourse-cdn.com/v4/letter/v/ecb155/32.png) [@vmh](https://discourse.julialang.org/u/vmh)\
**Post date:** [May 18, 2017, 4:26pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/15 "2017-05-18T16:26:22Z")

</div>

Hi @ChrisRackauckas, thank you so much for your response. Frankly, I didn’t know enough about parametric types that first I had to read about them to make sure I’d understand your solution. This is a pretty neat solution. Is your suggestion thread-safe? Any thoughts on this would be appreciated. @yuyichao, any thoughts on this approach (out of curiosity)?

> [@ChrisRackauckas](#):
>
> Instead of a function use a call overloaded type with the array as a field?
> 
> Define the type
> 
> ```julia
> immutable MyFunction{T}
> samples::T
> end
> 
> ```
> 
> Define a call
> 
> ```julia
> function (p::MyFunction)(dist)
> _draw_samples!(p.samples, dist)
> end
> 
> ```
> 
> Build the type with the cache  
> `draw_samples_1_allocation = MyFunction(Array{Float64, 1}(1000))`
> 
> Now use it like a function
> 
> ```julia
> dist = ...
> draw_samples_1_allocation(dist)
> 
> ```

---

<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 18, 2017, 4:37pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/16 "2017-05-18T16:37:49Z")

</div>

> [@vmh](#):
>
> Is your suggestion thread-safe?

Depends on how you’re using it. Are you going to be calling the same instance of `MyFunction` on each thread? Then no because they will be sharing the same cache. But in that case you just make one of these per thread, put them all in an array, and access them via `threadid()`. That will work because each function is locally keeping its own cache. Example:

```julia
cached_functions = [MyFunction(Array{Float64, 1}(1000)) for i in Threads.nthreads()]
Threads.@threads for i in 1:1000
  cached_functions[Threads.threadid()](dist)
end

```

Or just make sure you split access to `samples` by the `threadid` in another way.

---

<div class="post-metadata">

**Author:** ![vmh](https://avatars.discourse-cdn.com/v4/letter/v/ecb155/32.png) [@vmh](https://discourse.julialang.org/u/vmh)\
**Post date:** [May 18, 2017, 4:46pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/17 "2017-05-18T16:46:21Z")

</div>

Thank you, Chris, for your prompt response. Your example is very helpful (I’m learning a lot here).

---

<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 18, 2017, 10:29pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/18 "2017-05-18T22:29:51Z")

</div>

Not thread safe as there is only one samples array and all threads access it , so you might ask for certain distribution in one thread and by the time you use it another thread called the function with a different distribution.

to make it thread-safe i assume you would want to change the generating function.  
or better , make another generating function that handles thread safety through replication for example.

```julia
function fthread_safe(g)
    fv = [deepcopy(g) for i=1:Threads.nthreads()]
    G = Ref{typeof(fv)}(fv)
    (args...) -> G.x[Threads.threadid()](args...)
end

```

and since there is only one parameter you can write it very elegantly at the end of another generating function using the |\> operator (I dont know its name in english).

`draw_samples = fsamples_prealoc(draw_samples!,1000) |> fthread_safe`

---

<div class="post-metadata">

**Author:** ![musm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/musm/32/3675_2.png) [@musm](https://discourse.julialang.org/u/musm)\
**Post date:** [May 19, 2017, 4:28am UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/19 "2017-05-19T04:28:18Z")

</div>

[quote=“TsurHerman, post:18, topic:3726”]  
function fsamples\_prealoc(g,n,T = Float64)  
samples = Ref{Array{T, 1}}( Array{T, 1}(n) )  
(args…) → g(samples.x,args…)  
end[/quote]

```julia
function foo2(d)
    x = draw_samples2(d)
    2x
end

```

Looks like `foo2` is not type-stable ?

---

<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, 12:17pm UTC](https://discourse.julialang.org/t/generated-staged-function-to-avoid-initializing-an-array-multiple-times-for-monte-carlo-simulations/3726/20 "2017-05-19T12:17:25Z")

</div>

I didn’t understand what you meant … what is this foo2?

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.

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