diff --git a/src/ACME.jl b/src/ACME.jl
index 42d41da3..13715aa6 100644
--- a/src/ACME.jl
+++ b/src/ACME.jl
@@ -14,6 +14,37 @@ using IterTools
 
 import Base.getindex
 
+if VERSION ≥ v"0.6.0"
+    macro pfunction(sig, params, body)
+        esc(Expr(:function, Expr(:where, sig, params.args...), body))
+    end
+else
+    macro pfunction(sig, params, body)
+        ts = copy(params.args)
+        if VERSION < v"0.5.0"
+            for i in eachindex(ts)
+                if isa(ts[i], Expr) && ts[i].head === :comparison && ts[i].args[2] === :<:
+                    ts[i] = Expr(:<:, ts[i].args[1], ts[i].args[3])
+                end
+            end
+        end
+        esc(Expr(:function,
+                 Expr(:call, Expr(:curly, sig.args[1], ts...),
+                      sig.args[2:end]...),
+                 body))
+    end
+end
+
+macro expandafter(mc)
+    args = copy(mc.args)
+    for i in eachindex(args)
+        if isa(args[i], Expr) && args[i].head === :macrocall
+            args[i] = macroexpand(Compat.@__MODULE__, args[i])
+        end
+    end
+    esc(Expr(:macrocall, args...))
+end
+
 include("kdtree.jl")
 include("solvers.jl")
 
@@ -247,7 +278,7 @@ end
 # lines marked with !SV avoid creation of SparseVector by indexing with Ranges
 # instead of Ints; a better way for cross-julia-version compatibilty would be
 # nice; maybe Compat helps in the future...
-function topomat!{T<:Integer}(incidence::SparseMatrixCSC{T})
+@pfunction topomat!(incidence::SparseMatrixCSC{T}) [T<:Integer] begin
     @assert all(x -> abs(x) == 1, nonzeros(incidence))
     @assert all(sum(incidence, 1) .== 0)
 
@@ -289,7 +320,9 @@ function topomat!{T<:Integer}(incidence::SparseMatrixCSC{T})
     tv, ti
 end
 
-topomat{T<:Integer}(incidence::SparseMatrixCSC{T}) = topomat!(copy(incidence))
+@pfunction topomat(incidence::SparseMatrixCSC{T}) [T<:Integer] begin
+    topomat!(copy(incidence))
+ end
 topomat(c::Circuit) = topomat!(incidence(c))
 
 #mutable struct DiscreteModel{Solvers}
@@ -314,13 +347,15 @@ eval(Expr(:type, true, :(DiscreteModel{Solvers}), quote
     solvers::Solvers
     x::Vector{Float64}
 
-    @compat function (::Type{DiscreteModel{Solver}}){Solver}(circ::Circuit, t::Float64)
+    @expandafter @compat @pfunction (::Type{DiscreteModel{Solver}})(circ::Circuit,
+            t::Float64) [Solver] begin
         Base.depwarn("DiscreteModel{Solver}(circ, t) is deprecated, use DiscreteModel(circ, t, Solver) instead.",
                      :DiscreteModel)
         DiscreteModel(circ, t, Solver)
     end
 
-    @compat function (::Type{DiscreteModel{Solvers}}){Solvers}(mats::Dict{Symbol}, nonlinear_eqs::Vector{Expr}, solvers::Solvers)
+    @expandafter @compat @pfunction (::Type{DiscreteModel{Solvers}})(mats::Dict{Symbol},
+            nonlinear_eqs::Vector{Expr}, solvers::Solvers) [Solvers] begin
         model = new{Solvers}()
 
         for mat in (:a, :b, :c, :pexps, :dqs, :eqs, :fqprevs, :fqs, :dy, :ey, :fy, :x0, :q0s, :y0)
@@ -334,8 +369,8 @@ eval(Expr(:type, true, :(DiscreteModel{Solvers}), quote
     end
 end))
 
-function DiscreteModel{Solver}(circ::Circuit, t::Real, ::Type{Solver}=HomotopySolver{CachingSolver{SimpleSolver}};
-                               decompose_nonlinearity=true)
+@pfunction DiscreteModel(circ::Circuit, t::Real, ::Type{Solver}=HomotopySolver{CachingSolver{SimpleSolver}};
+                         decompose_nonlinearity=true) [Solver] begin
     mats = model_matrices(circ, t)
 
     nns = Int[nn(e) for e in circ.elements]
@@ -772,7 +807,8 @@ eval(Expr(:type, false, :(ModelRunner{Model<:DiscreteModel,ShowProgress}), quote
     ycur::Vector{Float64}
     xnew::Vector{Float64}
     z::Vector{Float64}
-    @compat function (::Type{ModelRunner{Model,ShowProgress}}){Model<:DiscreteModel,ShowProgress}(model::Model)
+    @expandafter @compat @pfunction (::Type{ModelRunner{Model,ShowProgress}})(
+            model::Model) [Model<:DiscreteModel,ShowProgress] begin
         ucur = Array{Float64,1}(nu(model))
         ps = Vector{Float64}[Vector{Float64}(np(model, idx)) for idx in 1:length(model.solvers)]
         ycur = Array{Float64,1}(ny(model))
@@ -782,9 +818,12 @@ eval(Expr(:type, false, :(ModelRunner{Model<:DiscreteModel,ShowProgress}), quote
     end
 end))
 
-ModelRunner{Model<:DiscreteModel}(model::Model) = ModelRunner{Model,true}(model)
-ModelRunner{Model<:DiscreteModel,ShowProgress}(model::Model, ::Val{ShowProgress}) =
+@pfunction ModelRunner(model::Model) [Model<:DiscreteModel] begin
+     ModelRunner{Model,true}(model)
+ end
+@pfunction ModelRunner(model::Model, ::Val{ShowProgress}) [Model<:DiscreteModel,ShowProgress] begin
     ModelRunner{Model,ShowProgress}(model)
+end
 
 """
     ModelRunner(model::DiscreteModel, showprogress::Bool = true)
@@ -798,8 +837,9 @@ By default `run!` for the constructed `ModelRunner` will show a progress bar to
 report its progress. This can be disabled by passing `false` as second
 parameter.
 """
-ModelRunner{Model<:DiscreteModel}(model::Model, showprogress::Bool) =
+@pfunction ModelRunner(model::Model, showprogress::Bool) [Model<:DiscreteModel] begin
     ModelRunner{Model,showprogress}(model)
+end
 
 """
     run!(runner::ModelRunner, u::AbstractMatrix{Float64})
@@ -845,18 +885,16 @@ To simulate a circuit without inputs, a matrix with zero rows may be passed.
 The internal state of the  underlying `DiscreteModel` (e.g. capacitor charges)
 is preserved accross calls to `run!`.
 """
-function run!{Model<:DiscreteModel}(runner::ModelRunner{Model,true},
-                                    y::AbstractMatrix{Float64},
-                                    u::AbstractMatrix{Float64})
+@pfunction run!(runner::ModelRunner{Model,true}, y::AbstractMatrix{Float64},
+                u::AbstractMatrix{Float64}) [Model<:DiscreteModel] begin
     checkiosizes(runner, u, y)
     @showprogress "Running model: " for n = 1:size(u, 2)
         step!(runner, y, u, n)
     end
 end
 
-function run!{Model<:DiscreteModel}(runner::ModelRunner{Model,false},
-                                    y::AbstractMatrix{Float64},
-                                    u::AbstractMatrix{Float64})
+@pfunction run!(runner::ModelRunner{Model,false}, y::AbstractMatrix{Float64},
+                u::AbstractMatrix{Float64}) [Model<:DiscreteModel] begin
     checkiosizes(runner, u, y)
     for n = 1:size(u, 2)
         step!(runner, y, u, n)
diff --git a/src/kdtree.jl b/src/kdtree.jl
index 1c792f4b..e46361ac 100644
--- a/src/kdtree.jl
+++ b/src/kdtree.jl
@@ -94,9 +94,11 @@ eval(Expr(:type, true, :(Alts{T}), quote
     number_valid::Int
 end))
 
-Alts{T}(p::Vector{T}) = Alts([AltEntry(1, zeros(p), zero(T))], typemax(T), 0, 1)
+@pfunction Alts(p::Vector{T}) [T] begin
+     Alts([AltEntry(1, zeros(p), zero(T))], typemax(T), 0, 1)
+ end
 
-function init!{T}(alts::Alts{T}, best_dist, best_pidx)
+@pfunction init!(alts::Alts{T}, best_dist, best_pidx) [T] begin
     alts.number_valid = 1
     alts.entries[1].idx = 1
     fill!(alts.entries[1].delta, zero(T))
@@ -161,7 +163,9 @@ function dequeue!(alts::Alts)
     return e
 end
 
-function enqueue!{T}(alts::Alts{T}, new_idx::Int, ref_delta::Vector{T}, delta_update_dim::Int, delta_update_val::T, new_delta_norm::T)
+@pfunction enqueue!(alts::Alts{T}, new_idx::Int, ref_delta::Vector{T},
+                    delta_update_dim::Int, delta_update_val::T,
+                    new_delta_norm::T) [T] begin
     if alts.number_valid == length(alts.entries)
         delta = copy(ref_delta)
         delta[delta_update_dim] = delta_update_val
diff --git a/src/solvers.jl b/src/solvers.jl
index 0c6e9ff1..f849c835 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -13,22 +13,21 @@ eval(Expr(:type, false, :(ParametricNonLinEq{F_eval<:Function,F_setp<:Function,F
     Jp::Matrix{Float64}
     J::Matrix{Float64}
     scratch::Scratch
-    @compat function (::Type{ParametricNonLinEq{F_eval,F_setp,F_calcjp,Scratch}}){
-            F_eval<:Function,F_setp<:Function,F_calcjp<:Function,Scratch
-        }(
+    @expandafter @compat @pfunction (::Type{ParametricNonLinEq{F_eval,F_setp,F_calcjp,Scratch}})(
             func::F_eval, set_p::F_setp, calc_Jp::F_calcjp, scratch::Scratch,
             nn::Integer, np::Integer
-        )
+        ) [F_eval<:Function,F_setp<:Function,F_calcjp<:Function,Scratch] begin
         res = zeros(nn)
         Jp = zeros(nn, np)
         J = zeros(nn, nn)
         return new{F_eval,F_setp,F_calcjp,Scratch}(func, set_p, calc_Jp, res, Jp, J, scratch)
     end
 end))
-ParametricNonLinEq{F_eval<:Function,F_setp<:Function,F_calcjp<:Function,
-                   Scratch}(func::F_eval, set_p::F_setp, calc_Jp::F_calcjp,
-                            scratch::Scratch, nn::Integer, np::Integer) =
+@pfunction ParametricNonLinEq(func::F_eval, set_p::F_setp, calc_Jp::F_calcjp,
+        scratch::Scratch, nn::Integer, np::Integer) [F_eval<:Function,
+        F_setp<:Function,F_calcjp<:Function,Scratch] begin
     ParametricNonLinEq{F_eval,F_setp,F_calcjp,Scratch}(func, set_p, calc_Jp, scratch, nn, np)
+end
 ParametricNonLinEq(func::Function, nn::Integer, np::Integer) =
     ParametricNonLinEq(func, default_set_p, default_calc_Jp,
                        (zeros(np), zeros(nn, np)), nn, np)
@@ -160,8 +159,8 @@ eval(Expr(:type, true, :(SimpleSolver{NLEQ<:ParametricNonLinEq}), quote
     tol::Float64
     tmp_nn::Vector{Float64}
     tmp_np::Vector{Float64}
-    @compat function (::Type{SimpleSolver{NLEQ}}){NLEQ<:ParametricNonLinEq}(
-            nleq::NLEQ, initial_p::Vector{Float64}, initial_z::Vector{Float64})
+    @expandafter @compat @pfunction (::Type{SimpleSolver{NLEQ}})(
+            nleq::NLEQ, initial_p::Vector{Float64}, initial_z::Vector{Float64}) [NLEQ<:ParametricNonLinEq] begin
         z = zeros(nn(nleq))
         linsolver = LinearSolver(nn(nleq))
         last_z = zeros(nn(nleq))
@@ -187,9 +186,10 @@ available Jacobians. Due to the missing global convergence, the `SimpleSolver`
 is rarely useful as such.
 """ SimpleSolver
 
-SimpleSolver{NLEQ<:ParametricNonLinEq}(nleq::NLEQ, initial_p::Vector{Float64},
-                                       initial_z::Vector{Float64}) =
+@pfunction SimpleSolver(nleq::NLEQ, initial_p::Vector{Float64},
+                        initial_z::Vector{Float64}) [NLEQ<:ParametricNonLinEq] begin
     SimpleSolver{NLEQ}(nleq, initial_p, initial_z)
+end
 
 set_resabstol!(solver::SimpleSolver, tol) = solver.tol = tol
 
@@ -254,13 +254,13 @@ eval(Expr(:type, true, :(HomotopySolver{BaseSolver}), quote
     start_p::Vector{Float64}
     pa::Vector{Float64}
     iters::Int
-    @compat function (::Type{HomotopySolver{BaseSolver}}){BaseSolver}(
-            basesolver::BaseSolver, np::Integer)
+    @expandafter @compat @pfunction (::Type{HomotopySolver{BaseSolver}})(
+            basesolver::BaseSolver, np::Integer) [BaseSolver] begin
         return new{BaseSolver}(basesolver, zeros(np), zeros(np), 0)
     end
-    @compat function (::Type{HomotopySolver{BaseSolver}}){BaseSolver}(
+    @expandafter @compat @pfunction (::Type{HomotopySolver{BaseSolver}})(
             nleq::ParametricNonLinEq, initial_p::Vector{Float64},
-            initial_z::Vector{Float64})
+            initial_z::Vector{Float64}) [BaseSolver] begin
         basesolver = BaseSolver(nleq, initial_p, initial_z)
         return HomotopySolver{typeof(basesolver)}(basesolver, np(nleq))
     end
@@ -326,15 +326,15 @@ eval(Expr(:type, true, :(CachingSolver{BaseSolver}), quote
     new_count::Int
     new_count_limit::Int
     alts::Alts{Float64}
-    @compat function (::Type{CachingSolver{BaseSolver}}){BaseSolver}(basesolver::BaseSolver,
-            initial_p::Vector{Float64}, initial_z::Vector{Float64}, nn::Integer)
+    @expandafter @compat @pfunction (::Type{CachingSolver{BaseSolver}})(basesolver::BaseSolver,
+            initial_p::Vector{Float64}, initial_z::Vector{Float64}, nn::Integer) [BaseSolver] begin
          ps_tree = KDTree(hcat(initial_p))
          zs = reshape(copy(initial_z), nn, 1)
          alts = Alts(initial_p)
          return new{BaseSolver}(basesolver, ps_tree, zs, 1, 0, 2, alts)
     end
-    @compat function (::Type{CachingSolver{BaseSolver}}){BaseSolver}(nleq::ParametricNonLinEq,
-            initial_p::Vector{Float64}, initial_z::Vector{Float64})
+    @expandafter @compat @pfunction (::Type{CachingSolver{BaseSolver}})(nleq::ParametricNonLinEq,
+            initial_p::Vector{Float64}, initial_z::Vector{Float64}) [BaseSolver] begin
         basesolver = BaseSolver(nleq, initial_p, initial_z)
         return CachingSolver{typeof(basesolver)}(basesolver, initial_p, initial_z, nn(nleq))
     end