# How to make a vector of parametric composite types?

Hello everyone,
I 've been using Julia for a few months now, and I really love it, especially the type system.

I’m working on a code to compute the optical reflectivity of a multilayer. The multilayer is a stack of homogeneous plane parallel media (layers). Each layer has a thickness and optical properties.

For now I’m using singleton types for the optical properties, but it seems this is not recommended and I would like to use value types instead. But I don’t manage to implement this.
I’m using Julia 0.5.0.

This is how I’m handling the simulations now.
I first define an abstract type OptProp

``````abstract  OptProp
``````

I then define singleton types for different materials

``````immutable Sic <:  OptProp end
immutable Al <:  OptProp end
immutable Vac <:  OptProp end

``````

This allows me to dispatch on the permittivity function to obtain the optical properties for each of those materials.
For instance for the optical properties of Aluminium I use:

``````function permittivity(material :: Al(), frequency)
some code
end
``````

and similarly for other materials.

I then define the layer type

``````immutable Layer
material :: OptProp
thickness :: Float64
end
``````

A multilayer is then a vector of the Layer type. As an example, an interface between vacuum and alumium would be defined as

``````interface = [Layer(Vac(),0.0); Layer(Al(),0.0)]
``````

More complex multilayers can be constructed adding layers with other materials and thickness.
This works fine, but I’m not sure this is the recommended way.

I tried to implement the “value type” -approach, but I didn’t succeed.
Here is what I tried:

``````immutable Mater{T}<: OptProp
end

immutable Layer{T}
material :: Mater{T}
thickness :: Float64
end
``````

This new version of the Layer type does not work.
If I do for a semi-infinite medium of aluminium:

``````Layer(Mater{:al},0.0)
``````

I get a MethodError : no method matching Layer{T}(::Type{Mater{:al}},:: Float64)

I have three questions:

1. Is there an overhead defining vectors of composite types containing singleton types (first approach)?
2. Is it recommended/more efficient to use value types for the material properties, instead of the singleton types?
3. If the answer to question 2 is yes, how would you actually implement it?

Olivier

1 Like

The error arises from a mismatch between the singleton type and the instance of the singleton type. You can avoid it by either modifying the constructor call Layer(Mater{:al}, 0.0) to look like “Layer(Mater{:al}(),0.0)” (i.e. creating an instance of the Mater{:al}) OR by changing the field in the type definition of Layer to “material :: Type{Mater{T}}”.

1 Like

That works.

It seems, however, that the value type approach generates much more allocation than the singleton type one.

Value type :

``````@time ml = [Layer(Mater{:vac},0.0) ; Layer(Mater{:al},0.0)]
0.000294 seconds (52 allocations: 2.125 KB)
``````

Where the layer type was defined as

``````immutable Layer{T}
material :: Type{Mater{T}}
thickness :: Float64
end
``````

with

``````immutable Mater{T} <: OptProp end
``````

For the singleton type approach I get:

``````@time ml = [Layer(Vac(),0.0) ; Layer(Al(),0.0)]
0.000012 seconds (7 allocations: 320 bytes)
``````

I used @code_warntype on both approaches and it seems the ‘value type’ method is type unstable while the singleton one is not.

Is this expected behaviour?
Thanks again

That is interesting, I get the same. I wonder why this is.
BTW is there a reason for not using the slightly (imho) simpler syntax

``````abstract OptProp

immutable Al <: OptProp end
immutable Vac <: OptProp end

immutable Layer{T <: OptProp}
material::Type{T}
thickness::Float64
end

@time ml = [Layer(Al,0.0) ; Layer(Vac,0.0)]
``````

@mkborregaard Thanks a lot for pointing out the simpler syntax.

With this third approach I get the same timing and allocations as for the ‘value type’ approach

``````@time ml = [Layer(Al,0.0) ; Layer(Vac,0.0)]
0.000267 seconds (52 allocations: 2.125 KB)
``````

Doing `@code_warntype`, indicates a type instability on `Type{T}`

For completeness I have to add that the results in my answer to @braydenware were obtained on a windows machine.
The results in this post were obtained on a different machine with UBUNTU 14.04.
But apparently the results are the same in both cases.

From this, it seems I should stick with the singleton type approach with the syntax

``````abstract  OptProp

immutable Sic <:  OptProp end
immutable Al <:  OptProp end
immutable Vac <:  OptProp end

immutable Layer
material :: OptProp
thickness :: Float64
end

``````

Which I don’t find very satisfying since for every new material, I have to introduce a new immutable type.

Furthermore, if I understand the documentation correctly, for better performance, fields of composite types should be concrete types or parametric which is in contradiction with the above example.
Am I missing some crucial point here?

So I am not sure exactly what you want. Maybe instead you are looking for something like this?

``````@enum Material Sic Al Vac

immutable Layer
material::Material
thickness::Float64
end

permittivity(l::Layer) = permittivity(Val{l.material})

function permittivity(::Type{Val{Al}})
println("The permittivity of Aluminum is grand indeed")
end

function permittivity(::Type{Val{Sic}})
println("But even grander is Silicium")
end

#Then
l = Layer(Al, 4)
permittivity(l)
l2 = Layer(Sic, 4)
permittivity(l2)
``````
``````@time ml = [Layer(Al,0.0) ; Layer(Vac,0.0)]
0.000008 seconds (7 allocations: 336 bytes)
``````
3 Likes

Woow, thanks a lot!
This is really nice.
It will make the code much more user friendly.

Thanks again
Olivier

1 Like

I think this will cause calculating the `permitivity` to be slow, since it will be doing dynamic (runtime) dispatch.

Would it be possible to use a dictionary for the permitivities instead?

1 Like

You’re right.

``````@enum Material Sic Al Vac

immutable Layer
material::Material
thickness::Float64
end

permittivity(l::Layer) = permittivity(Val{l.material})
permittivity (generic function with 1 method)

function permittivity(::Type{Val{Sic}})
2
end

using BenchmarkTools

@benchmark permittivity(l2)
BenchmarkTools.Trial:
memory estimate:  16 bytes
allocs estimate:  1
--------------
minimum time:     331.807 ns (0.00% GC)
median time:      335.516 ns (0.00% GC)
mean time:        346.447 ns (0.45% GC)
maximum time:     8.241 μs (93.95% GC)
--------------
samples:          10000
evals/sample:     223
time tolerance:   5.00%
memory tolerance: 1.00%
``````

vs

``````abstract OptProp

immutable Al <: OptProp end

immutable Layer{T <: OptProp}
material::Type{T}
thickness::Float64
end

function permittivity{T<:Al}(x::Layer{T})
2
end
l1 = Layer(Al, 0.)

using BenchmarkTools
@benchmark permittivity(l1)
BenchmarkTools.Trial:
memory estimate:  0 bytes
allocs estimate:  0
--------------
minimum time:     12.587 ns (0.00% GC)
median time:      12.725 ns (0.00% GC)
mean time:        12.890 ns (0.00% GC)
maximum time:     46.025 ns (0.00% GC)
--------------
samples:          10000
evals/sample:     999
time tolerance:   5.00%
memory tolerance: 1.00%
``````

Thanks! (and sorry @Olivier_Merchiers)

1 Like

@Olivier_Merchiers can you give a little more info - e.g. why do you use dispatch on the permittivity function?

Thanks for this correction.

But there remains the issue with allocation and time for the vector of layers with this implementation:

``````abstract OptProp

immutable Al <: OptProp end
immutable Vac <: OptProp end

immutable Layer{T <: OptProp}
material::Type{T}
thickness::Float64
end

@benchmark ml = [Layer(Al,0.0); Layer(Vac,0.0)]
BenchmarkTools.Trial:
memory estimate:  1.97 kb
allocs estimate:  48
--------------
minimum time:     10.937 μs (0.00% GC)
median time:      11.656 μs (0.00% GC)
mean time:        11.802 μs (0.00% GC)
maximum time:     72.133 μs (0.00% GC)
--------------
samples:          10000
evals/sample:     1
time tolerance:   5.00%
memory tolerance: 1.00%

``````

While with

``````abstract OptProp

immutable Al <: OptProp end
immutable Vac <: OptProp end

immutable Layer
material::OptProp
thickness::Float64
end

@benchmark ml = [Layer(Al(),0.0); Layer(Vac(),0.0)]
BenchmarkTools.Trial:
memory estimate:  160.00 bytes
allocs estimate:  3
--------------
minimum time:     43.866 ns (0.00% GC)
median time:      47.852 ns (0.00% GC)
mean time:        61.659 ns (21.21% GC)
maximum time:     2.199 μs (96.82% GC)
--------------
samples:          10000
evals/sample:     990
time tolerance:   5.00%
memory tolerance: 1.00%

``````

This is quite important for me, since I need to pass such multilayer structures to other functions that will compute reflectivity and other physical quantities.

Finally, doing the test with the permittivity function

``````function permittivity(x::Al)
2
end

``````

I get

``````@benchmark permittivity(Al())
BenchmarkTools.Trial:
memory estimate:  0.00 bytes
allocs estimate:  0
--------------
minimum time:     1.409 ns (0.00% GC)
median time:      1.680 ns (0.00% GC)
mean time:        1.711 ns (0.00% GC)
maximum time:     90.449 ns (0.00% GC)
--------------
samples:          10000
evals/sample:     1000
time tolerance:   5.00%
memory tolerance: 1.00%

``````

Which is almost 10 times faster.

@mkborregaard Sorry, I missed this last question before submitting the other reply.

Yes of course.
I need to make reflectivity calculations for many different materials. These materials have optical properties that are defined by their permittivity function.
With every material corresponds a unique permittivity function which returns a complex number.
This number is used to compute the reflectivity.

To keep my code generic, and not using if-then-else in a big permittivity function containing all possible definitions for each material, I wanted to use multiple dispatch to get the right permittivity method corresponding with each material.

To give some examples:
For vacuum this would be

``````permittivity(material::Vacuum) = 1.0 +im*0.0
``````

For frequency dependent (hence w in the parameter list) optical properties you would have something like

``````function permittivity(material::Cbn,w)
eps_fin = 4.46 + 0.0*im
return eps_fin*(w^2-w_lo^2 + im*gamma*w)/(w^2-w_to^2 + im*gamma*w)
end
``````

Let me know if this explanation is not clear

1 Like

That’s why I suggested to use a dictionary instead of the big `if-else`.

How would you implement such a Dict?

If the permittivity for each material was a constant value, I think I could build such a dict.
But this is not the case here, since the permittivity for each material is a function of the frequency as shown in this example:

``````function permittivity(material::Cbn,w)
eps_fin = 4.46 + 0.0*im
return eps_fin*(w^2-w_lo^2 + im*gamma*w)/(w^2-w_to^2 + im*gamma*w)
end
``````

where w is the frequency in rad/s and can take any real value between 0 and infinity.

That would imply a dict of functions?

Well maybe your original model is best, I really don’t see what is wrong with it But the reason your last permittivity function is fast is because there is no dispatch at all.
You could have concrete parametric functions in your first example like

``````    abstract OptProp

immutable Al <: OptProp end
immutable Vac <: OptProp end

immutable Layer{T <: OptProp}
material::T
thickness::Float64
end

# and then
permittivity{T <: Al}(l::Layer{T}) = 2

# or the clearer
permittivity(l::Layer) = permittivity(l.material)
permittivity(x::Al) = 2
``````

Note that in benchmark it is important how to set it up for comparison

``````l = Layer(Al(), 0.)
@benchmark permittivity(l)
BenchmarkTools.Trial:
memory estimate:  0 bytes
allocs estimate:  0
--------------
minimum time:     15.057 ns (0.00% GC)
median time:      15.652 ns (0.00% GC)
mean time:        15.877 ns (0.00% GC)
maximum time:     46.872 ns (0.00% GC)
--------------
samples:          10000
evals/sample:     998
time tolerance:   5.00%
memory tolerance: 1.00%

@benchmark permittivity(Layer(Al(), 0.))
BenchmarkTools.Trial:
memory estimate:  0 bytes
allocs estimate:  0
--------------
minimum time:     1.500 ns (0.00% GC)
median time:      1.502 ns (0.00% GC)
mean time:        1.506 ns (0.00% GC)
maximum time:     3.348 ns (0.00% GC)
--------------
samples:          10000
evals/sample:     1000
time tolerance:   5.00%
memory tolerance: 1.00%
``````
2 Likes

You need to interpolate:

``````@benchmark permittivity(\$l)
``````

This then gives the same time.

1 Like

See here for details about interpolation.

1 Like

This seems to be the best solution.

Thanks, this is great - in trying to help I also (re-)learn some useful things myself. I only hope this hasn’t created too much noise for anyone reading this thread

1 Like

Thanks a lot for your help @mkborregaard, @dpsanders.

This was really instructive. Things are much clearer now to me.

With the interpolation, the benchmark results for the vector of Layers become also identical for the different approaches even with the enum type approach. Too bad the enum types are slow for the permittivity calculation. I really liked them

1 Like