How to make a vector of parametric composite types?

[mkborregaard]

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.

[MA_Laforge]

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

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

[MA_Laforge]

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)

[MA_Laforge]

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

[MA_Laforge]

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.

[Olivier_Merchiers]

Thanks a lot @mkborregaard for showing me the macro for the nice syntax. This is what makes me like julia so much

Thanks a lot @MA_Laforge for this very detailed explanation. I stumbled upon this issue of type stability for the vectors, but I could not formulate the question properly. You kind of read my mind!
Just to be sure I understand correctly the last part:

Non-parametrized : generates type a stable vector since all elements are of type Layer.
Parametrized : not type stable for a vector in which the materials are different, type stable when all materials are identical

Case 1: type stable, but if-statement is called.

myfun([LayerAl, LayerVac])

Case 2: not type stable but same speed as case 1, because of the if-statement in case 1?

myfun([LayerAlP, LayerVacP])

Case 3: All layers are identical, hence type stable but not useful.

myfun([LayerAlP, LayerAlP])

Case 4 : No speedup compared to case 1 as you get going from case 2 to case 3

myfun([LayerAl, LayerAl])

To finish :

@MA_Laforge

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

You are right about Sic, but with Al I did mean Aluminum.
And I definitely will change the names as you suggest.

[MA_Laforge]

Good to hear I did not add too much confusion: Since you appear to have understood the broad strokes, I will go a bit more into detail (sorry, I had to over-simplify to get my point across):

Re: type stable

Technically, I believe “type stable” simply means Julia can figure out the exact concrete type of a function’s return value knowing only the types of its inputs parameters.

So… technically, you are not supposed to say “generates a type stable vector”… I think the more correct statement is that a particular Vector constructor (a function) was type stable. Again: this is true when Julia can determine the type of the result simply by looking at the types of the constructor’s input parameters.

In reality: type stability is only a subset of the problem (albeit one that people have difficulty with). What you really should be asking is: “is Julia able to figure out which (type-sensitive) method to call given the types of your function’s input parameters”. This is especially true when you start executing long, complex loops.

Now for your questions:

Non-parametrized : generates type a stable vector since all elements are of type Layer.

Technically, that is what I said in order to keep close to your wording. Indeed: it is very likely that Julia “calls” a type-stable constructor when you create the non-parameterized Vector{Layer}.

However, what is more important is how efficiently Julia will then use this vector as opposed to the vector made up of LayerP{:MATERIAL} (the parameterized version). What I was saying here is that both solutions are pretty much of the same efficiency… that is unless you have a vector where all the layers are of type LayerP{:MATERIALX}. In that case, Julia does not have to generate an implicit “if” statement to dispatch on the LayerP object type - because it knows the exact type of all the elements (ie, they are all layers that use MATERIALX).

I will refrain to answer your question about non-parameterized types because it will not help here. Feel free to re-ask once you have read this post.

Case 1: type stable, but if-statement is called.
myfun([LayerAl, LayerVac])

Again, that is what I said in an attempt to keep close to your wording.

Technically, since I don’t know what dosomething() returns… so I cannot say if this myfun() method is type stable… (Remember, type stability is only a subset of the issue).

What I can say is that in this non-parameterized version of Layers, the vector:

lstack = [LayerAl, LayerVac]

Is a fully concrete Vector{Layer}. That means that whenever we call a function on any of those Layer objects themselves, Julia knows exactly which function to call (does not have to do an implicit “if” statement yet…). For example, if you just need to sum up the thicknesses, you can call:

thickness(v::Vector{Layer}) = sum([l.thickness for l in v])

Julia will be very efficient because it does not need an implicit “if” statement anywhere.

However, when you finally run an operation that depends on the the material itself, that’s when Julia starts to add in implicit “if” statements. For example:

effective_permittivity(v::Vector{Layer}) = #Something that depends on .material

Sadly, very little can be done to improve Julia’s calculations on Vectors when you are structuring the solution the way you are. There might be alternatives out there, but your solution would likely be significantly different. I personally find your type hierarchy to be relatively good at the moment.

The main reason why you will not see an improvement by parameterizing your Layer object is the following:

• In most cases (not all), efficiencies gained by parameterizing a type are observed when you have to execute the same code over a large number of that exact type (ie a large Vector of that type). In these cases, Julia will be compiling type-specific code that will often be faster.
• You could in theory parameterize Layer for different Materials in order to gain computational efficiency. But this added efficiency will mostly be observed when you execute code on a Vector of Layer{MaterialX} (all the parameter types are the same). Since there is no real interest in solving problems where a stack is made up of all of the same parameter (or material, in this case), I don’t see a reason for this “optimization” (added complexity).

[MA_Laforge]

@Olivier_Merchiers: You are right about Sic, but with Al I did mean Aluminum.

I withdraw my statement: If you are actually using the chemical composition of the material, I personally think this is a valid shorthand (definitely a widely accepted standard) - as long as the names do not conflict with anything.

In that case, I would use SiC instead of Sic - because that’s how you typically write its chemical formula.

…But I still have a problem with Mater to mean Material. I thought it was the name of a layer that was used to “mate” two other layers - maybe some sort of epoxy, or something.

Also

I am a bit sorry if my first set of posts was misleading a bit. It can be hard to communicate a new concept when you “don’t really speak the same language”. I was trying to target your intent more than addressing the actual words you were saying.

Feel free to ask your questions a second time now that I tried to clear up some of my over-simplifications.

[ChrisRackauckas]

@MA_Laforge

Technically, I believe “type stable” simply means Julia can figure out the exact concrete type of a function’s return value knowing only the types of its inputs parameters.

I would expand the definition a bit. A good working definition is something like this:

A function in Julia is type-stable if the type of every variable inside of the function is deducible and static.

Your definition only requires that the return type is deducible. In many cases that can be the easy part. For example, if it’s returning a Vector{Int}, sometimes that’s easy to deduce, but the fact that you have type instabilities within the function can make things slow (since it will have to do a bunch of unboxing before converting things to ints). The instabilities just won’t “leak out” because it’s still able to deduce the return type, but it can still be the root of performance issues.

Essentially this means that, for type-stable functions:

1. Make each variable have a single type within a function. (If you want to do a conversion, have a x and an _x version, or something like that. Or use dispatches to convert before an inner function).

2. Make each variable type deducible at compile-time. This means that things which invoke the runtime like @eval can cause issues and should have type-assertions. Another thing that can cause issues here are global variables, or calling other functions for which the return type cannot be deduced.

If you follow those 2 rules, then the compiler can understand the entire function as though all of the variables and function calls are static, and it will optimize it like a C or Fortran function. That is what is typically meant by type-stable.

[MA_Laforge]

Thanks for that @ChrisRackauckas!

Your definition for “type-stable” does sound much better. I must have mis-interpreted it to mean only a subset of what it is.

@Olivier_Merchiers:

I maintain that my suggestions still have alot of validity, but I think you should use the definition of “type stability” given by @ChrisRackauckas.

So my major errata is that all my suggestions do fall under “making your code type stable”… and I would currently no longer claim that “type stablility is only a subset of the problem” (since I have no examples outside its scope).

[MA_Laforge]

A Cleaner Solution?

Though I forsee no real issue with the performance of the code you presented at the beginning, I too would like to propose my own alternative. This alternative is probably not significantly faster, but it does seem a bit more natural, and might simplify your overall solution:

abstract PermittivityModel

type Debye <: PermittivityModel
	coeff1::Real
	#...
end
type DjordjevicSarkar <: PermittivityModel
	coeff1::Real
	#...
end

type Layer
	model::PermittivityModel
	thickness::Float64
end

#Now, your materials are simply constants - instead of separate types:
const SiC_Debye = Debye(...)
const SiC_DS = DjordjevicSarkar(...)
const SiC = SiC_DS #Default SiC model
#...

This way, you don’t actually have to write separate permittivity() functions for each material. You just need one per PermittivityModel:

function permittivity(model::DjordjevicSarkar, frequency)
     #Some code
end

You might also want to add in a property/type hierarchy to more naturally represent isotropic/anisotropic materials as well… But I don’t remember how it relates to the models, etc - nor do I know if this is relevant to your solution.

And building your layer stackup is pretty much the same:

mystack = [Layer(SiC_DS, 1e-6), Layer(Vacuum_DS, 5e-6)]

Now here is the neat thing: There is now a potentially useful reason to parameterize Layer:

type Layer{T<:PermittivityModel}
	model::T
	thickness::Float64
end

Unlike before where the parameter was the material itself, there are actual useful cases where you might want a solution for a multi-layer stack where all materials use the same permittivity model. This way, Julia might remove a couple of implicit if statements along the way.

That being said, I am not convinced your particular solution will really notice this reduction in clock cycles. I would instead pick this type hierarchy because I feel it is a bit cleaner and removes a few redundancies in modelling the materials.

[Olivier_Merchiers]

Thanks a lot again @MA_Laforge and @ChrisRackauckas for taking the time to explain in such a detailed way the subtleties of Julia for a newbie
I will keep these posts very carefully, because they add understanding I could not extract from the docs. Thanks again!

@MA_Laforge

In most cases (not all), efficiencies gained by parameterizing a type are observed when you have to execute the same code over a large number of that exact type (ie a large Vector of that type). In these cases, Julia will be compiling type-specific code that will often be faster.

My stacks are quite simple (about 3 layers), so as you point out, I guess the effect will not be really visible.

@MA_Laforge

In that case, I would use SiC instead of Sic - because that’s how you typically write its chemical formula.

I totally agree. [quote=“MA_Laforge, post:28, topic:2439”]
I am a bit sorry if my first set of posts was misleading a bit. It can be hard to communicate a new concept when you “don’t really speak the same language”. I was trying to target your intent more than addressing the actual words you were saying.
[/quote]

Please, don’t apologize. I think you put things in a very pedagogical way. So thanks again for that.

Thanks also for proposing the cleaner solution. I was going to ask about it. You read again my mind!
I like it much more that way.
I actually half started implementing something similar, but your way is much cleaner.

[mkborregaard]

Maybe this section of the manual is relevant to the discussion http://docs.julialang.org/en/latest/manual/performance-tips.html#The-dangers-of-abusing-multiple-dispatch-(aka,-more-on-types-with-values-as-parameters)-1

[andyferris]

If your stacks are “quite simple” (you said 3 layers?), then I would strongly recommend using a tuple to describe a collection of layers, not a vector.

The reason is that the compiler tracks the type of each element of a tuple, but not for vectors. It can then areange dispatch on each element at compile time. You can achieve optimal speed this way - as if you wrote a program specifically for each stack (Julia’s great, this way). Some functions like map and ntuple can help here, but I imagine you’ll need to write at least one “Lisp style” recursive function to iterate through the layers.

If this seems of interest, let me know, and I’ll be more detailed.

[mkborregaard]

I certainly would like to hear more about this approach, if only out of interest.

[andyferris]

OK, I’ll try be quick. First, “lispy recursion” is demonstrated by the pattern here: 0.6: how do you broadcast `call` over a vector of functions? - #12 by andyferris. (Note the timing differences with/without type stability in that thread.)

It uses an @inline function and splatting ... to cleverly build up, piece by piece, each answer to f(x) into a tuple.

[andyferris]

So the idea is to make a stack as a being a tuple of Layers. However, to keep everything nice and conceptualized, we can make a Stack type. There are lots of ways of organizing things, but this is one suggestion. Here I’ll perform a “reduction” over the layers - I’ll give something akin to the optical thickness here, but you can do more complex things (reflectance would presumably be a bit more complex).

abstract AbstractMaterial

abstract Sic <: AbstractMaterial
abstract Al <: AbstractMaterial
...

const eps_0 = ...

function permittivity(::Type{Sic}, frequency)
    ...
end
...

immutable Layer{Material <: AbstractMaterial}
   thickness::Float64
end

immutable Stack{Layers <: Tuple{Vararg{Layer}}}
    layers::Layers
end
Stack(x...) = Stack{typeof(x)}(x) # convenience constructor

s = Stack(Layer{Sic}(1e-5), Layer{Al}(2e-5))

function optical_thickness(s::Stack, frequency)
    _sum_optical_thickness(frequency, 0.0, s.layers...)
end

_sum_optical_thickness(frequency, thickness) = thickness
@inline function _sum_optical_thickness{Material}(frequency, thickness, layer::Layer{Material}, other_layers...) 
    _sum_optical_thickness(frequency, thickness + layer.thickness*permittivity(Material, frequency)/eps_0, other_layers...)
end

optical_thickness(s)

I haven’t tested any of this, and I probably got the physics wrong, but I’d like to see a performance comparison.