How to make a vector of parametric composite types?

[Olivier_Merchiers]

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?

Many thanks in advance.
Olivier

[braydenware]

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}}”.

[Olivier_Merchiers]

Thank you very much for your answer.
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

[mkborregaard]

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)]

[Olivier_Merchiers]

@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?

[mkborregaard]

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)

[Olivier_Merchiers]

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

Thanks again
Olivier

[dpsanders]

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?

[mkborregaard]

@dpsanders

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

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)

[mkborregaard]

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

[Olivier_Merchiers]

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.

[Olivier_Merchiers]

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

@mkborregaard

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

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 
    w_lo    = 2.451e14 # rad/s 
    w_to    = 1.985e14 # rad/s 
    gamma   = 9.934e11 # rad/s 
    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

[dpsanders]

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

[Olivier_Merchiers]

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 
    w_lo    = 2.451e14 # rad/s 
    w_to    = 1.985e14 # rad/s 
    gamma   = 9.934e11 # rad/s 
    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?
Did I misunderstand your suggestion?

[mkborregaard]

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%

[dpsanders]

@mkborregaard

permittivity(Layer(Al(), 0.))

You need to interpolate:

@benchmark permittivity($l)

This then gives the same time.

[dpsanders]

See here for details about interpolation.

[dpsanders]

This seems to be the best solution.

[mkborregaard]

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

[Olivier_Merchiers]

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