Best practice to use lots of pre-allocated arrays in a function

General Usage
, ,
1

I am still quite new to Julia, working on numerical modeling and I am starting to write a bigger project than usual.

I have a function that takes, as arguments, lots of pre-allocated arrays and parameters and update them in different ways using different functions. I am wondering what would be the best practice to store them and update them?

In pseudocode:

array1 = zeros(nx)
array2 = zeros(nx)
array3 = zeros(nx)
...

function update!(array1, array2, array3,...., parameters)
    updatearray1!(array1, parameters..)
    updatearray2!(array2, parameters..)
    updatearray3!(array3, parameters..)
end

What bother me here is that I have lots of arguments for my function. I was wondering what is the common practice for storing pre-allocated arrays?

My idea would be to store them all in an array, input that in the function, and use @view at the beginning of the function for each separate array, and inputing that in the internal functions. Is there a better way to do that?

I have started to work with Parameters.jl for storing the constant of my models and it is working great. So I am looking for something similar to Parameters.jl where I can easily update the arrays with their new values.

Tell me if it is confusing. I am just looking for the best practices for dealing with lots of preallocated arrays.

Thanks.

2

Personally, I like creating a “workspace” struct to hold this kind of data. Something like:

struct Myworkspace
  x::Vector{Float64}
  y::Vector{Float64}
end

function update!(workspace::MyWorkspace, parameters)
  # do stuff with `workspace.x` etc.
end

Then you can also write a convenience function if you want:

function update(parameters)
  work = MyWorkspace(...)
  update!(work, parameters)
end

which is less efficient but very convenient for testing or interactive usage.

2 Likes
3

Note that the arrays within a struct are still mutable, so you can use that package for organizing your arrays. Or a named tuple:

julia> myarrays = ( arr1 = [1,2], arr2 = [3,4] )
(arr1 = [1, 2], arr2 = [3, 4])

julia> function f(myarrays)
           myarrays.arr1[1] = 0
       end
f (generic function with 1 method)

julia> f(myarrays)
0

julia> myarrays
(arr1 = [0, 2], arr2 = [3, 4])
1 Like
4

Thx, that’s what I was looking for! I am still quite new to structure, but I am already liking it.

5

Sorry, I replied at the same time. Thx for the tips, that’s actually quite neat!

EDIT

I tried to do what you did with struct, and arrays don’t seem to be updated in that case, at least when I am doing it that way:

julia> using Parameters

julia> @with_kw struct test
           u0::Vector{Float64} = zeros(5)
       end
test

julia> function f(test)
           @unpack u0 = test()
           u0 .= 2
       end
f (generic function with 1 method)

julia> f(test)
5-element Vector{Float64}:
 2.0
 2.0
 2.0
 2.0
 2.0

julia> test().u0
5-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0

Did I miss something?

6

This defines a new type with the name test. By default, its constructor creates an instance which holds a freshly-allocated vector

julia> @with_kw struct test
           u0::Vector{Float64} = zeros(5)
       end

This defines a function that accepts a callable and binds it to the name test within its body

julia> function f(test)
           @unpack u0 = test() # calls the `test` that is passed as argument
           u0 .= 2
       end

The defined function is called on a type constructor:

julia> f(test)

u0 within f is a freshly created array.

test() allocates a fresh array:

julia> test().u0

If you come from Python, you may have the impression that a mutable default argument for a function or a default mutable attribute for a class is created on function/class definition. In Julia, the default argument is an expression which is evaluated each time the function or constructor is called.

What @rdeits and @lmiq advised can be done as:

julia> @with_kw struct TestStruct
           u0::Vector{Float64} = zeros(5)
       end
test

julia> function f(test)
           @unpack u0 = test
           u0 .= 2
       end

julia> test_obj = TestStruct();

julia> f(test_obj)

julia> test_obj.u0
2 Likes
7

Great, it is pretty clear for me now. Thx for deconstructing my code, it helps. It also explains why I had 1 allocation on my f function before.

So my mistake was to use TestStruct() as an argument. Even though I can define default values using Parameters.jl, TestStruct() is still a Type constructor, and I should then build an object from it.

8

As of Julia 1.8.0, working with typed globals is perfectly fine. Despite the disadvantages of using global variables, the introduction of this nice feature can help reduce much of the noise caused by passing the same variables to multiple functions. The following function update!(..) causes 0-allocations:

nx::Int64 = 10
array1::Vector{Float64} = ones(nx)
array2::Vector{Float64} = ones(nx)
array3::Vector{Float64} = ones(nx)

updatearray1!(array,parameters) = array .+= sum(parameters)
updatearray2!(array,parameters) = array .+= minimum(parameters)
updatearray3!(array,parameters) = array .+= maximum(parameters)

function update!(parameters)
   updatearray1!( array1, parameters )
   updatearray2!( array2, parameters )
   updatearray3!( array3, parameters )
end

update!((1,2,3))

 array1 = [7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0] 
 array2 = [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0] 
 array3 = [4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0]

and

@btime update!((1,2,3))
  25.000 ns (0 allocations: 0 bytes)
1 Like