Structs and parametric constructors

This is a performance benchmark I did on a smaller example

julia> struct g{F<:Function}
       h::F end

julia> struct G
       h::Function end

julia> a = g(x->x^2-2x^7+1)
g{getfield(Main, Symbol("##21#22"))}(getfield(Main, Symbol("##21#22"))())

julia> b = G(x->x^2-2x^7+1)
G(getfield(Main, Symbol("##23#24"))())

julia> gh(x::g{F},a) = x.h(a)
gf (generic function with 1 method)

julia> Gh(x::G,a) = x.h(a)
Gf (generic function with 1 method)

julia> @btime gh(a,17+im)
  32.813 ns (2 allocations: 64 bytes)
-761386735 - 332091598im

julia> @btime Gh(b,17+im)
  33.525 ns (2 allocations: 64 bytes)
-761386735 - 332091598im

I think you missed a where F here. But in any case, cf (note @noinline)

julia> @noinline gh(x::g{F},a) where F = x.h(a)
gh (generic function with 1 method)

julia> @noinline Gh(x::G,a) = x.h(a)
Gh (generic function with 1 method)

julia> @btime gh($a,17+im)
  15.977 ns (0 allocations: 0 bytes)
-761386735 - 332091598im

julia> @btime Gh($b,17+im)
  34.889 ns (2 allocations: 64 bytes)
-761386735 - 332091598im
1 Like

Well, then let me give you the numbers.

With

struct fractal
  xx
  yy
  myfunction::Function
  derivative::Function
  number_of_roots::Int
  roots
end

I get

julia> @code_warntype iterationstep(Fractal,zz)
Body::Any
31 1 ─ %1 = (Base.getfield)(Fractal, :myfunction)::Function                                               │╻ getproperty
   │   %2 = (%1)(z)::Any                                                                                  │ 
   │   %3 = (Base.getfield)(Fractal, :derivative)::Function                                               │╻ getproperty
   │   %4 = (%3)(z)::Any                                                                                  │ 
   │   %5 = (%2 / %4)::Any                                                                                │ 
   │   %6 = (z - %5)::Any                                                                                 │ 
   └──      return %6              
julia> @benchmark iterationstep($Fractal,$zz)
BenchmarkTools.Trial: 
  memory estimate:  224 bytes
  allocs estimate:  7
  --------------
  minimum time:     122.678 ns (0.00% GC)
  median time:      131.516 ns (0.00% GC)
  mean time:        147.694 ns (10.15% GC)
  maximum time:     50.438 μs (99.71% GC)
  --------------
  samples:          10000
  evals/sample:     911

and using

struct fractal{F<:Function,G<:Function}
  xx
  yy
  myfunction::F
  derivative::G
  number_of_roots::Int
  roots
end

the numbers are

julia> @benchmark iterationstep($Fractal,$zz)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     25.975 ns (0.00% GC)
  median time:      26.199 ns (0.00% GC)
  mean time:        26.472 ns (0.00% GC)
  maximum time:     49.658 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     995

not showing the nice but lengthy @code_warntype and noting that “factor two” was a bit optimistic.

However, the original question was about the matching rules for the type parameters. I know realise that my confusion originated from the failing

struct fractal{F<:Function}
  xx
  yy
  myfunction::F
  derivative::F
  number_of_roots::Int
  roots
end

which was actually my first try based on the examples in the manual. In the working case

struct fractal{F<:Function,G<:Function}
  xx
  yy
  myfunction::F
  derivative::G
  number_of_roots::Int
  roots
end

of course the myfunction::F vs. derivative::G makes the match unique. It does not do any positional matching. So the identifier in the struct, i.e. F or G has to be unique which is enforced. Probably one can make sure that the type they reference are the same with an inner constructor ?

Edit: typos …

Yes, indeed, I just pushed an update to Fatou where used an inner constructor, like this

struct fractal{F<:Function,G<:Function}
  xx
  yy
  myfunction::F
  derivative::G
  number_of_roots::Int
  roots
  ## inner constructor here
  fractal(x,y,mf::Function,df::Function,nr,r) = new{typeof(mf),typeof(df)}(x,y,mf,df,nr::Int,r)
end

And what is the real world speed difference for you ?

Let me summarize what I now believe the answer to the original rather confused question is:

The parametric type arguments in the struct definition, i.e. {F<:Function,G<:Function} are matched to the actual outer constructor arguments on a first match basis using the actual argument type and consuming the type parameter afterwards. The match to the fields in the struct is based purely on the uniqueness of the identifiers, i.e. F and G in the example, no positional matching applies. There is no risk of mis-identifying xx or yy which are ranges in the constructor call with these rules.

Edit: Deleted a stupid mistake stemming from translating between actual code and MWE.

I also found out that changing roots to be an array and typing roots as <:AbstractArray was beneficial as the compiler suddenly was able to infer the type in the following function which tries to figure out which root it had converged to:

  ergebnis=0
  for ii=1:Fractal.number_of_roots
    diff=Fractal.roots[ii]-z;
    if ( abs(diff) < Fractal.accuracy*1e2 )
      ergebnis=ii
      break
    end
  end
  return ergebnis
end

Without that I had some Any’s in @code_warntype there.

When I only add the Function parameters to the Define type, I get basically 0% to 10% improvement.

My original timings look like this

julia> newton(:(z^3-1);n=500)|>fatou;
  0.240635 seconds (1.69 M allocations: 58.407 MiB)

With the Define{FT<:Function,QT<:Function,CT<:Function} definition, I don’t get much improvement

julia> newton(:(z^3-1);n=500)|>fatou;
  0.221504 seconds (1.53 M allocations: 50.724 MiB)

However, because in my situation there is additional program logic, I added additional Bool type paramters to the Fatou.Define{FT,QT,CT,M,N} object, which gave me a 10x performance improvement. This time the parameters have additionally Bool pair to help decide on Mandelbrot and Newton fractal mode.

julia> newton(:(z^3-1);n=500)|>fatou;
  0.020936 seconds (542 allocations: 36.784 KiB)

This is because of my iteration function definition, which depends on these true or false values, by including them in the Fatou.Define type parameters, the function can specialize on those values

# define function for computing orbit of a z0 input
function nf(K::Define{FT,QT,CT,M,N},z0::Complex{Float64}) where {FT,QT,CT,M,N}
    M ? (z = K.seed) : (z = z0)
    zn = 0x00
    while (N ? (K.Q(z,z0)::Float64>K.ϵ)::Bool : (K.Q(z,z0)::Float64<K.ϵ))::Bool && K.N>zn
        z = K.F(z,z0)::Complex{Float64}
        zn+=0x01
    end; #end
    # return the normalized argument of z or iteration count
    return (zn::UInt8,z::Complex{Float64})
end

So now it computes much, much faster. However, in my case, having a Function parameter alone did not do the trick, it was also necessary to have other parameters also.

Sprechen sie Deutsch?

OK, from this and roots and related function above (which is still missing the function definition … sorry) the lesson to me is to take the performance tip to heart and to actually check with @code_warntype :slight_smile:

Yes. This is basically the explanation for the mess in the previous post, I had to (or better tried to) translate everything from one language to the other and juggled several versions at the same time … Sorry about that.

Edit: OK, here is the whole untranslated function (to avoid further mess ups), it is probably better than editing above post again which I was contemplating.

function vergleich_zu_nullstellen(Fraktal,z)
  ergebnis=0
  if (Fraktal.index < 5)
    for ii=1:Fraktal.anzahl_nullstellen
      diff=Fraktal.nullstellen[ii]-z;
      if ( abs(diff) < Fraktal.genauigkeit*1e2 )
        ergebnis=ii
        break
      end
    end
  else
    ergebnis=Int(round(real(z)/(pi/2)))
  end 
  return ergebnis
end

This includes an additional if, as you probably know the Newton fractal can also be defined for e.g. sin() which has a single real root and I should simplify that, I know.

1 Like

Well, I am a german speaker also, but probably not everyone on this website is.

What you could do is write a single while loop, which has the if statement inside the condition, as I have done in my code shown above.

Amazingly, with this improvement the difference is enormous for larger n, with old definition

julia> newton(:(z^3-1);n=1500)|>fatou;
 19.215167 seconds (65.07 M allocations: 1.817 GiB, 7.66% gc time)

However, with new definition

julia> newton(:(z^3-1);n=1500)|>fatou;
  0.175065 seconds (541 allocations: 36.706 KiB)

julia> newton(:(z^3-1);n=2500)|>fatou;
  0.490495 seconds (542 allocations: 36.784 KiB)

julia> newton(:(z^3-1);n=3500)|>fatou;
  1.004239 seconds (541 allocations: 36.706 KiB)

Clearly, the memory consumption and allocations are greatly reduced. I believe this is because in my original version, the @threads parallel processing had to assign the memory for Fatou.Define object every time to determine the Bool parameters for Mandelbrot and Newton mode. Having these Bool and Function available in the type parameter for Fatou.Define enables the compiler to leave out the full reference to the original Fatou.Define object, I suspect, and is probably inlining everything now.

So if you include enough type parameters to make the original struct obsolete, I think that’s the trick, so that the bigger object does not need to be re-assigned for parallel computations.

1 Like

I may be restating the conclusion you came to, but for further clarity, the order of the type parameters has no relation to the order of the fields in the struct.

Your struct here would be just as valid with the type parameters rearranged as:

struct fractal{G<:Function,F<:Function}
  xx
  yy
  myfunction::F
  derivative::G
  number_of_roots::Int
  roots
end

When calling the type constructor, that is when the argument order matters; each field is bound to the corresponding argument (e.g. the first arg in the constructor gets bound to xx, the second arg is bound to yy, etc.).

Additionally, in case it wasn’t clear, the reason that

didn’t work is because F represents a single, unique type. Assuming your myfunction and derivative are different,

typeof(myfunction) != typeof(derivative)

and no method matching a type signature with different functions for myfunction and derivative will be found.

1 Like

Hello,

that is nice to hear !

Thank you for spelling this out again. I believe I have somewhat of a problem with the syntax. For me it would currently be much clearer if the parametric types were given in the outer constructor definition and the outer constructor interface had to be explicitely defined once (I realize that means usually more pointless typing), something like

fractal(xx<:Type,yy<:Type,myfunction<:Type,derivative<:Type,3<:Type,roots<:Type)

Probably I just need to work more with the actual syntax, i.e. the ```{…}``, and get accustomed to it :slight_smile:

But may be you or someone else can answer a final question. I now switched everything to using the ranges in Fractal.xx and Fractal.yy which made things slow as they were untyped. It turns out that apparently (including here the unrelated typing of roots for completeness, see above posts)

struct fractal{a<:Any,b<:Any,F<:Function,G<:Function,r<:AbstractArray}
  xx::a
  yy::b
  myfunction::F
  derivative::G
  number_of_roots::Int
  roots::r
end

is sufficient to restore performance, so I assume the compiler is then able to infer type again. Or to give a self contained example

using BenchmarkTools
struct gamma
       a
end
struct alpha{a<:Any}
       a::a
end
function delta(Gamma)
       zz=sqrt(Gamma.a[3])
end
Gamma=gamma(range(1.0,stop=2.0,length=10))
@btime delta($Gamma)
Gamma=alpha(range(1.0,stop=2.0,length=10))
@btime delta($Gamma)

which gives

julia> include("new.jl")
  42.343 ns (2 allocations: 32 bytes)
  16.811 ns (0 allocations: 0 bytes)

@code_warntype indicates that the first case defaults to Any for everything while in the second case the whole range function appears to be inlined.

What would be a better type for the ranges if there is one ? And why do I have to do a seemingly trivial <:Any by hand, I am sure there is a good reason but I would like to know it ? And why should I not type everything as <:Any then (ok, catching errors, but for performance) ?

1 Like

This is because your second example includes the type of a in the type of alpha. Technically, you don’t need to have <:Any in your definition, because you can just leave that out, because Any is assumed

struct alpha{a}
       a::a
end

Should have same behavior. This definition stores the type of a as part of the type.

julia> alpha(range(1.0,stop=2.0,length=10)) |> typeof
alpha{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}

The difference between gamma and alpha is that alpha types will be able to specialize on their type parameter, so that code generated for that type can be inlined. With gamma, the information about the type of a is not present when the code is generated. The type parameter helps the compiler infer field type.

I understand that perfectly :slight_smile: I should reword my question:

Why do I have to do that manually in the first place ?

Adding <:Any is something I would expect as a default behaviour, avoiding the performance pitfalls. This is completely astonishing to me, especially since doing this as a default appears to be trivial, so there must be a very compelling reason not to do it which I would like to know.

This is because Julia compiles specialized methods for specific types. Since your gamma definition has no information about the type of the field a, the compiler will not be able to infer its type. However, because you specifically included the type parameter for alpha, the compiler can infer more about it. You don’t have to use a parameter, because you could specify what the type should be like this too:

struct alpha
  a::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
end

But then, it is not parametric, it is statically typed.

There is no specific reason not to do it, but it also isn’t necessary, since it is automatically assumed.

1 Like

I think you may have misunderstood something: performance suffers when you don’t parametrize the type at all (see the performance tips). Whether you then constrain the type parameter (or leave it unconstrained, either implicitly or with <: Any) is orthogonal to performance.

1 Like

Well, yes :slight_smile: Please let me try again, may be I am a bit slow :slight_smile:

when I do a b=range(...) the type of a is well known

  julia> b=range(1.0,stop=2.0, length=10)
1.0:0.1111111111111111:2.0

julia> typeof(b)
StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}

I pass this to gamma. At the time I call the default constructor with b as an argument the type is known. I do not see a reason why julia could not specialise the type based on the argument by default. Instead I apparently have to explicitely “allow” (?)/ “force” (?) the type specialization by setting <:Any as the parametric type known at construction time (or typing statically with the StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}} which is quite difficult to manually type correctly at first try in by hand :-).

Or to paraphrase:
The way everything else in julia appears to be designed is to specialize at runtime on the concrete types of the arguments given by compiling a specialized version for each combination of types of the actual arguments. Opposed to that the behaviour for constructors is apparently to throw the type information away instead by default, unless not explicitely requested

This appears to be a very strange decision, the behaviour of constructors is diametrically opposed compared to “normal” functions ?

Unfortunately, I can’t give you an explanation of why these design decisions were made, but it seems reasonable to me, since there might only be a few field types that need to be specialized on. For example, in the Fatou.Define object, I only need to specialize on the Function parameters and the Bool parameters, all the other fields do not need type parameters, because they are not referenced in the iteration function, and so it is not necessary to know their types when the iteration function is compiled. Not all of the fields are performance critical, some of them are only used for printing the label of a plot.

As you can see, parametric types can grow to be quite large. If Julia automatically specializes on any ambiguous type, then all of the types would automatically have a huge amount of parameters in their definition, which can be a lot of un-necessary information to pass around to the compiler and in other places. I assume this is why it is that way, for sake of minimalism?

A struct with an untyped field b must be able to hold a b of any type without the type of the struct itself changing. Julia does still specialize the constructor just like any other function, but the struct that it returns must not, by definition, depend on the type of b. If you want the type of the struct to depend on the type of b then you need a type parameter.

A couple of other notes:

  1. There are lots of cases in which untyped fields can be desirable, like when a field might need to store elements of varying types and accessing that field is not critical for performance. You could imagine making parametric structs the default, but aside from breaking all sorts of code, you would still presumably want some way to refer to those parameters by name, in which case you’d be back to exactly the system we have.

  2. You don’t need <:Any when declaring a type parameter. struct F{T} and struct F{T <: Any} are identical.

  3. If the type of a particular field is fixed but simply hard to type, you can just give it a shorter name:

const T = typeof(1:2)

struct F
  b::T
end

# this also works:
struct F
  b::typeof(1:2)
end
4 Likes

Thank you for the explanation. As far as I understand what you write (I do not have a computer science background) the argument might simply boil down to the immutability of the struct ? I see the utility of the current behaviour for general storage containers.

That’s a nice trick.

Mm, not exactly. The same argument applies to both mutable and immutable structs: there are legitimate use cases where you might not want a given field to have a concrete type, so Julia allows this. The fact that it’s the default is ultimately a design choice, but I think it’s a good one: if the language allowe a field to be silently or anonymously parameterized (i.e. didn’t have a corresponding {T} in the struct definition), then we would have no way to refer to the actual type associated with that field, which would make life harder in many cases.

1 Like

Right. I think I see where this argument is going, it starts to make sense. Probably with more experience this will become entirely clear. Thank you !