2
votes

I have been writing stochastic PDE simulations in Julia, and as my problems have become more complicated, the number of independent parameters has increased. So what starts with,

myfun(N,M,dt,dx,a,b)

eventually becomes

myfun(N,M,dt,dx,a,b,c,d,e,f,g,h)

and it results in (1) messy code, (2) increased chance of error due to misplaced function arguments, (3) inability to generalise for use in other functions.

(3) is important, because I have made simple parallelisation of my code to evaluate many different runs of the PDEs. So I would like to convert my functions into a form:

myfun(args)

where args contains all the relevant arguments. The problem I am finding with Julia, is that creating a struct containing all my relevant parameters as attributes slows things down considerably. I think this is due to the continual accessing of the struct attributes. As a simple (ODE) working example,

function example_fun(N,dt,a,b)
  V = zeros(N+1)
  U = 0
  z = randn(N+1)
  for i=2:N+1
    V[i] = V[i-1]*(1-dt)+U*dt
    U = U*(1-dt/a)+b*sqrt(2*dt/a)*z[i]
  end
  return V
end

If I try to rewrite this as,

function example_fun2(args)
  V = zeros(args.N+1)
  U = 0
  z = randn(args.N+1)
  for i=2:args.N+1
    V[i] = V[i-1]*(1-args.dt)+U*args.dt
    U = U*(1-args.dt/args.a)+args.b*sqrt(2*args.dt/args.a)*z[i]
  end
  return V
end

Then while the function call looks elegant, it is cumbersome to rework accessing every attribute from the class and also this continual accessing of attributes slows the simulation down. What is a better solution? Is there a way to simply 'unpack' the attributes of a struct so they do not have to be continually accessed? And if so, how would this be generalised?

edit: I am defining the struct I use as follows:

struct Args
  N::Int64
  dt::Float64
  a::Float64
  b::Float64
end

edit2: I have realised that structs with Array{} attributes can give rise to a performance difference if you do not specify the dimensions of the array in the struct definition. For example, if c is a one-dimensional array of parameters,

struct Args_1
  N::Int64
  c::Array{Float64}
end

will give far worse performance in f(args) than f(N,c). However, if we specify that c is a one-dimensional array in the struct definition,

struct Args_1
  N::Int64
  c::Array{Float64,1}
end

then the performance penalty disappears. This issue and the type instability shown in my function definitions seem to account for the performance difference I encountered when using a struct as the function argument.

2
It absolutely should not be slower just because you are accessing the fields of a struct. As the answerer has suggested, it is possible you have introduced type instability in your struct definition. Could you post how you define your struct? Also, one other alternative to your own struct would be splatting. There might be a performance hit to using splatting though - I'm not certain, you'd have to try it. Personally, I think the struct solution is the right one.Colin T Bowers
I agree, my initial thoughts were that it should not be slower. The example function I give only has the args method being very slightly slower on average, but when the complexity of the simulation increases (ie PDEs), the time difference can be at least an order of magnitude. I'm currently using @benchmark from BenchmarksTools to more rigorously see how exactly the computation time changes between the simulations I use. I'm not familiar with splatting, so I'll look into that, thanks.Arsibalt
Sorry just saw this and may not get a chance to properly respond over weekend (am on phone now). Given how you are defining the struct there should be zero time difference. As in, the compiled code will be the same since the types of the parameters are known beforehand. I strongly suspect the problem is being caused by something else. Hopefully you can spot something using btimeColin T Bowers
In the background Julia is building a tuple equivalent to that struct, so its exactly the same time. You can use Parameters.jl to make this easier as well. If the timing is slower, it's because of type-instability or other user issues, not the design itself.Chris Rackauckas
Parameters.jl solves the tidiness aspect quite well, thanks. The difference in speed does seem to be due to type-instabilities that I had not accounted for.Arsibalt

2 Answers

2
votes

In your code there is a type instability, related to U which is initialized as an 0 (integer), but if you replace it with 0. (floating point number), the type-instability disapears.

For the original versions (with "U=0"), function example_fun takes 801.933 ns (for the parameters 10,0.1,2.,3.) and example_fun2 925.323 ns (for similar values).

In the type-stable version (U=0.), both take 273 ns (+/5 ns). Thus this a substantial speed-up and there is no more a penalty of combining the arguments in the type args.

Here is the complete function:

function example_fun2(args)
     V = zeros(args.N+1)
     U = 0.
     z = randn(args.N+1)
     for i=2:args.N+1
       V[i] = V[i-1]*(1-args.dt)+U*args.dt
       U = U*(1-args.dt/args.a)+args.b*sqrt(2*args.dt/args.a)*z[i]
     end
     return V
end
2
votes

Maybe you did not declare the types of the parameters of the type declaration of args?

Consider this small example:

struct argstype
    N
    dt
end
myfun(args) = args.N * args.dt

myfun is not type-stable can the type of the return type cannot be inferred:

@code_warntype myfun(argstype(10,0.1))
Variables:
  #self# <optimized out>
  args::argstype

Body:
  begin 
      return ((Core.getfield)(args::argstype, :N)::Any * (Core.getfield)(args::argstype, :dt)::Any)::Any
  end::Any

However, if you declare the types, then code becomes type-stable:

struct argstype2
    N::Int
    dt::Float64
end

@code_warntype myfun(argstype2(10,0.1))
Variables:
  #self# <optimized out>
  args::argstype2

Body:
  begin 
      return (Base.mul_float)((Base.sitofp)(Float64, (Core.getfield)(args::argstype2, :N)::Int64)::Float64, (Core.getfield)(args::argstype2, :dt)::Float64)::Float64
  end::Float64

You see that the inferred return type of Float64. With parametric types (https://docs.julialang.org/en/v0.6.3/manual/types/#Parametric-Types-1), your code still remains generic and type-stable at the same time:

struct argstype3{T1,T2}
    N::T1
    dt::T2
end

@code_warntype myfun(argstype3(10,0.1))
Variables:
  #self# <optimized out>
  args::argstype3{Int64,Float64}

Body:
  begin 
      return (Base.mul_float)((Base.sitofp)(Float64, (Core.getfield)(args::argstype3{Int64,Float64}, :N)::Int64)::Float64, (Core.getfield)(args::argstype3{Int64,Float64}, :dt)::Float64)::Float64
  end::Float64