How to make a vector of parametric composite types?

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

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

Well maybe your original model is best, I really don’t see what is wrong with it :slight_smile: 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 :slight_smile:

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 :slight_smile:

1 Like

If it was the syntax you liked about the enum type approach you could always define a little macro to define the abstract and immutable types:

macro myenum(args...)
    exs = Expr[]
    push!(exs, :(abstract $(args[1])))
    for arg in args[2:end]
        push!(exs, :(immutable $(arg) <: $(args[1]) end))
    end
    esc(Expr(:block, exs...))
end

The enum example looked 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{Sic}})
    println("pompuous nonsense")
end

With the myenum macro you could write the last, recommended, approach as

@myenum Material Sic Al Vac

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

permittivity(l::Layer) = permittivity(l.material)

function permittivity(x::Al)
    println("pompuous nonsense")
end

Which I actually like even better.

1 Like

Hi @Olivier_Merchiers

If I understand correctly, I think you have settled with solution:

abstract OptProp

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

immutable Layer
          material::OptProp
          thickness::Float64
 end 

I think this is a good idea. I don’t see a speed issue here (I will elaborate below).

I do suggest you change your names, though. I noticed @mkborregaard thought you were talking about aluminum & silicium, but my guess is that you were talking about alumina and silicon carbide (SiC) - otherwise you would have just called your type Si (unless you were worried of confusions with the french “if” :slight_smile: ).

I would also rename OptProp to be Material… since this Al & Vac are not optical properties, but materials.

1 Like

Material{:Alumina} solution

(vs: immutable Alumina < Material; end)

In your particular case, I don’t believe there is anything significantly different (performance-wise) with the Julia implementation the Material{} solution than with the immutable Alumina solution.

I use both solutions quite often myself. I often prefer the immutable Alumina solution simply because it avoids mistakes:

You can accidentally construct objects with typos using parametric types, unless you are careful to raise exceptions when the wrong one is created:

m = Material{:Alumina}() #Great
m = Material{:Aluminq}() #You get a silent error here

This is not an issue using the immutable Alumina solution:

m = Alumina() #Great
m = Aluminq() #Julia will raise an error (very helpful)

Type stability & layers

When you define your Layer object as:

immutable Layer
	material::Material
	thickness::Float64
end 

You can now build material stacks of Vector{Layer} (which are type stable):

interfaceAB = [LayerA ; LayerB]
typeof(interfaceAB)
--> Array{Layer, 1}

Now, as soon as you parameterize Layer, you get type instability:

immutable Layer{T}
	material::Material{T}
	thickness::Float64
end

interfaceAB = [LayerA ; LayerB]
typeof(interfaceAB)
--> Array{Layer, 1}

Note that this time, Layer represents a “not fully concrete type” (almost like an abstract type). A type is only concrete when it is qualified with its parameters (ex: Layer{:Alumina}):

layerAl = Layer(Material{:Alumina}, 1)
interfaceAB = [layerAl; layerAl]
typeof(interfaceAB)
--> Array{Layer{:Alumina}, 1}

This time interfaceAB will generate “type stable” code - because it is a vector of truly concrete types.

There in lies the problem: you will only get “type stable” code if all your types are concrete (ie your materials are all the same).

1 Like

Speed Issues

Okay, now after saying all of this, what about speed?

As I just mentionned, you can only get truly “type stable” code using this methodology if you are building a stackup of the same material (sort of a useless case).

So what you are really tring to do is trying to avoid writing your own if statements - by leveraging Julia’s type system. I do this all the time, and I think it is great (I use this instead of using Dicts myself).

I will try to illustrate what changes between the 2 solutions:

#Parameterized:
immutable MaterialP{T}; end
typealias AluminaP MaterialP{:Alumina}
typealias VacuumP MaterialP{:Vacuum}
immutable LayerP{T}
	material::MaterialP{T}
end
LayerAlP = LayerP(AluminaP())
LayerVacP = LayerP(VacuumP())

#Non-parameterized:
abstract Material
immutable Alumina <: Material; end
immutable Vacuum <: Material; end
immutable Layer
	material::Material
end
LayerAl = Layer(Alumina())
LayerVac = Layer(Vacuum())
function myfun(v::Vector{LayerP})
	return [dosomething(l) for l in v]
end
function myfun(v::Vector{Layer})
	return [dosomething(l) for l in v]
end

Now, the following might be considered type stable, but there is an implicit “if” stmt called just before calling dosomething()… that is unless Julia can statically figure out what types you have (ex: if it inlines myfun()):

#Case 1
myfun([LayerAl, LayerVac])

As far as I know, Julia’s type system would generate very similar code/speed if we used an array of parameterized Layers:

#Case 2
myfun([LayerAlP, LayerVacP])

Where Julia can make its optimization, though, is if it understands that all the elements are of the same type. That way, it does not have to do its implicit if before calling dosomething():

#Case 3
myfun([LayerAlP, LayerAlP])

… Not often a very useful stack, though (you would probably merge these materials yourself into a single thickness).

Now here is where you loose speed with the non-parameterized system: When you create a stack Vector using the the same Layer object twice, you don’t actually get a speedup wrt Case 1:

#Case 4
myfun([LayerAl, LayerAl])

… Again, that’s not really a useful stack, so I don’t see the point in optimizing for this.

Hope this helps.