Generic programming – writing functions that work for multiple different argument types – is straightforward in Julia. To recap, let’s assume we want to write a function that sums all elements of a vector. We could naively write

```
function vintsum(xs::Vector{Int})
s = 0
for x in xs
s += x
end
s
end
```

Of course, we might also want a function that works for floating-point vectors, so we write

```
function vfloatsum(xs::Vector{Float64})
s = 0.0
for x in xs
s += x
end
s
end
```

This gets boring quickly – the functions are almost identical except for their argument types. Thus the preferred way to implement this function in Julia is

```
function vsum{T<:Number}(xs::Vector{T})
s = zero(T)
for x in xs
s += x
end
s
end
```

This function accepts a vector of any numeric type. We restrict the
argument to numeric vectors so that addition is defined for the type
`T`

. Note that we also use the function `zero(T)`

to initialise the
sum `s`

; this is equivalent to writing `T(0)`

, but is more general,
and is preferred in Julia.

The above is standard in Julia, and could have be taken straight from a Julia tutorial. Nothing new here.

Let’s increase the level of genericity. The function `vsum`

above
works for any numeric vector – but what if the input is not a vector,
but a set? So here we go:

```
function setsum{T<:Number}(xs::Set{T})
s = zero(T)
for x in xs
s += x
end
s
end
```

Again, the function looks identical except for the type of its
argument. In addition to vectors and sets, there are many other
interesting collections – Julia itself features a host of array
types, and the `DataStructures`

package offers many additional
container types. We want to generalize the container type as well.

Our first attempt might be

```
function csum0{C<:TypeConstructor, T<:Number}(xs::C{T})
s = zero(T)
for x in xs
s += x
end
s
end
```

Here, a `TypeConstructor`

is something like `Vector`

, which takes
another type as argument (e.g. `Int`

) to produce a type, as in
`Vector{Int}`

. `Vector`

by itself is not and cannot be the type of any
object in Julia; objects can only have concrete types such as
`Vector{Int}`

, `Vector{Float64}`

, or `Vector{Any}`

. Since `Vector`

is
used to construct a type, we call it here a *type constructor*.
(`Vector`

is in Julia also an abstract type that is a supertype of all
types `Vector{X}`

; this is different, and is not relevant in this
context.)

However, Our `csum0`

function above doesn’t work because Julia’s
pattern matching algorithm for function dispatch does not handle type
constructor expressions such as `C{T}`

. We need to write things
differently:

```
function csum{CT}(xs::CT)
T = eltype(CT)
T<:Number || error("Argument element type $T is not a subtype of Number")
s = zero(T)
for x in xs
s += x
end
s
end
```

Instead of matching a type constructor `C`

and an element type `T`

, we
need to match a single type `CT`

. This means we lose the ability to
restrict `T`

to numeric types, and we need to check this condition
explicitly in the body of the function.

We used the function `eltype`

to obtain the element type `T`

from
`CT`

. As written above, the function `csum`

looks like code that is
commonly found in Julia.

You might notice that we were able to obtain the element `T`

from the
argument type `CT`

, but not the type constructor `C`

. In fact, there
is no standard way in Julia to obtain `C`

– but there is a way to get
(almost) the same effect (see below).

First, let’s look at an example where we actually need to know the
type constructor `C`

. Here we invert each element of a vector:

```
function vinv{T<:Number}(xs::Vector{T})
R = typeof(inv(one(T)))
rs = Vector{R}(length(xs))
for i in eachindex(rs)
rs[i] = inv(xs[i])
end
rs
end
```

`inv(x)`

calculates the reciprocal value `1/x`

. For integer `x`

, the
inverse is a floating-point value, i.e. it has different type. We thus
cannot just use `Vector{T}`

to hold the result – we need to calculate
the type `R`

of the result. We do so by examining a sample
calculation: `one(T)`

yields the value `1`

of type `T`

, `inv`

calculates its reciprocal, and the type of the result is the type `R`

we need to use as result type. Such type calculations are also quite
common in Julia.

On a side node, such type calculations are also commonplace in templated C++ code, but they usually look more complex there:

```
typedef std::decay_t<decltype(inv(declval<T>()))> R;
```

This is unfortunate because it likely drives people away from generic programming in C++. All the more power to a fresh approach!

We now generalize `vinv`

to work with an arbitrary container type `C`

.
Since, as mentioned above, Julia does not provide a way to obtain a
type constructor `C`

from the argument type `CT`

, we need to use a
work-around:

```
function cinv{CT}(xs::CT)
T = eltype(CT)
T<:Number || error("Argument element type $T is not a subtype of Number")
R = typeof(inv(one(T)))
rs = similar(xs, R)
for i in eachindex(rs)
rs[i] = inv(xs[i])
end
rs
end
```

The function `similar`

creates a new container from an existing
container (here `xs`

) and a new type (here `R`

). It also uses the
length of `xs`

to allocate an equal-sized new container. In a way,
`similar`

implements a function that takes a type `C{T}`

and a type
`R`

, and returns a new type `C{R}`

, except that it expects a
*container object* instead of a *container type* as argument, and that
it also returns a new container object instead of just its type.

So far, we used generic objects as arguments and return values. Let’s
go one step further and pass generic *types* as arguments. As example
we’ll be using a “generator” function, i.e. a function that generates
a container as return value from a non-container input. Specifically,
let us write a function that takes a container type constructor and a
single element as input, and returns a container with exactly this
element.

For example, a function that creates a single-element vector would be

```
function vunit{T}(x::T)
rs = Vector{T}()
push!(rs, x)
rs
end
```

A function that works for many containers would be

```
function cunit1{C,T}(::Type{C}, x::T)
rs = C{T}()
push!(rs, x)
rs
end
```

Note that `C`

is not the container type (such as `Vector{Int}`

), but
only the type constructor (e.g. `Vector`

), so that we need to use the
expression `C{T}`

as container type.

It turns out that Julia does not have a uniform internal
representation for type constructors (such as `Vector`

or `Set`

), and
we thus need to use the generic `Type`

as argument type; in Julia, all
type constructors are subtypes of `Type`

.

We can call `cunit1`

e.g. as `cunit1(Vector, 1)`

or `cunit1(Set, 2)`

,
and the result is as expected.

Unfortunately, this function `cunit1`

is not very useful. While it is
straightforward to call `cunit1`

with a manually-provided type
constructor, there is no way to obtain a type constructor from a
container type (or a container object): As mentioned above, there is
no function in standard Julia that takes the type `Vector{Int}`

as
input, and returns the type constructor `Vector`

. We need to look for
an alternative definition of `cunit1`

. (Of course, such a function
could be added to Julia, and I argue that it should.)

In C++, the template pattern matching mechanism is more powerful than in Julia, since one can define template specializations. This makes it possible to determine a type constructor from a type. C++ calls these type constructors “type templates”, and if you use a type template in a type template, you are using “template template parameters” (sic!), a confusing terminology.

It turns out that template template parameters in C++ have a serious disadvantage: the compiler often cannot determine whether two template template parameters are identical, and then needs to generate multiple identical template instantiations. This is highly undesirable. Of course, one can use in C++ the same work-around as in Julia, which I will describe now.

Our work-around is simply to use a type `C{X}`

instead of a type
constructor `C`

, using an arbitrary type `X`

that we will then ignore.
In many cases, it will be convenient to use the types `Void`

or `Any`

for this. (The type `Union{}`

also has its charm, as it by definition
prevents non-empty containers.)

Our actual implementation of `cunit`

looks like this:

```
function cunit{C, T}(::Type{C}, x::T)
rs = similar(C(), T)
push!(rs, x)
rs
end
```

We can call `cunit`

e.g. as `cunit(Vector{Void}, 1)`

or
`cunit(Set{Void}, 2)`

, with result types `Vector{Int}`

and `Set{Int}`

,
respectively. Due to two quirks of Julia – every type constructor is
also an abstract type, and containers usually define fall-back
constructors for abstract types as well – the calls ```
cunit(Vector,
1)
```

and `cunit(Set, 2)`

are also legal, and are interpreted as
`cunit(Vector{Any}, 1)`

and `cunit(Set{Any}, 2)`

, respectively, with
the same results as above.

Note that we needed to call `C()`

in the call to `similar`

,
constructing a throw-away empty container of type `C`

. This is
unfortunate, and could be avoided by extending `similar`

to accept
types as arguments as well, and not just objects.

So why would this higher-order genericity be interesting in practice?
As mentioned above, “regular” generic programming (as used in the
functions `csum`

or `cinv`

) has obvious advantages, and is commonplace
in Julia. But where would a function such as `cunit`

be useful, apart
from type-theoretical discussions on monads, and blog entries
comparing Julia to C++?

In my own work, I came across this issue while experimenting with how vector spaces would best be represented in Julia. In large numerical simulations of e.g. coalescing binary black holes, linear algebra and vectors spaces are everpresent, yet the algebraic concept of a “vector space” is usually not defined anywhere in the source code. Instead, various objects are defined, combined, added, and scaled in an ad-hoc manner that often leads to quite tedious code.

Many types represent vector spaces – arrays obviously do, and so do
certain tuples, most containers, and also many user-defined types. A
“vector space” should thus be represented as *trait* (or *concept* in
C++), i.e. as a property that a type might have. Each type that is a
vector space will need to implement certain functionality: create a
null vector, scale a vector, add two vectors, etc.

In addition to these operations that act on vectors, there exist also
operations on vector spaces: Given two vector spaces `V`

and `U`

, one
can create their sum space, product space, a power space `V^n`

for a
non-negative integer `n`

, and even a fiber vector space `V(U)`

. Those
are vector spaces again, and all these definitions are commonly used
in numerical simulations. (Adaptive mesh refinement introduces sum
spaces where either a coarse or a fine grid is used in a region,
multi-block discretizations are product spaces, power spaces define
grids and meshes, and fiber spaces express grids and meshes holding
vectors at each element.)

Given this, I was looking for generic implementations to create sum, product, power, and fiber spaces. That is, I wanted to implement e.g. a type

```
immutable VProduct0{V1, V2, S}
v1::V1{S}
v2::V2{S}
end
```

to representd a product space of `V1`

and `V2`

(essentially a tuple).
`V1`

and `V2`

represent vector spaces; they should be type
constructors that can be applied to either a scalar type, or to a
another vector space (!): `Vector{Int}`

is a vector space,
`Vector{Vector{Int}}`

is one as well. We thus want to write

```
typealias TwoVectors{T} VProduct{Vector, Vector, T}
```

and then use e.g. `TwoVectors{Float64}`

to hold some data.

Unfortunately, this type definition is not legal in Julia. In a
generic type definition such as this, `V1`

and `V2`

can be datatypes,
but they cannot be type constructors. One thus needs an alternative
way to represent type constructors, and then to define functions
`vnull`

, `vadd`

, `vscale`

, etc. that implement the vector space
operations for `VProduct`

.

To give a maybe more concrete motivation from linear algebra:
Block-structured matrices appear naturally in many algorithms. If one
needs a block-structured matrix where each block is dense, and the
blocks are sparsely distributed (many blocks are zero), then one
should be able to write `TwoLevelMatrix{SparseMatrix, DenseMatrix}`

.
Implementing `TwoLevelMatrix`

then leads to the problems described
here, since both `SparseMatrix`

and `DenseMatrix`

are type
constructors, not types.

Generic programming is ubiquitous in Julia. Higher-order generic programming, i.e. using generic type constructors instead of just generic types, seems alluring, and is able to provide elegant solutions to otherwise tedious coding exercises, at least in other languages. It remains to be seen in how far this is currently possible in Julia. I think most of the bits and pieces are already there, although certain things are still more cumbersome than they should be.

Julia’s competitor in this respect is not C++. While C++ templates are more powerful than Julia’s generic types, and thus better suited in theory, its strange and convoluted template syntax and type semantics make it quite difficult to use in practice. Haskell, of course, supported higher-order generic programming already a decade ago, but is coming from a very different angle. In Haskell, one needs to understand and implement every corner case correctly before the code compiles, and before one can begin to experiment. In Julia, where types are objects that can be manipulated at run time, experimentation is much easier, but one foregoes the safety of a much stronger typing harness. I conjecture that Julia’s approach makes it easier than Haskell’s to work in large, interdisciplinary open-source collaborations.