text:
"""
    R2d(theta)

Two-dimensional rotation matrix by an angle `theta`.
"""
function R2d(theta)
    SMatrix{2, 2}([cos(theta) -sin(theta); sin(theta) cos(theta)])
end

"""
    tblg_commensurate_angle(m, r)

Returns the commensurate angle for a commensurate twisted bilayer graphene structure labelled by co-prime integers `m` and `r`, according to the notation of PRB 86, 155449 (2012) [arXiv:1202.1088].
"""
function tblg_commensurate_angle(m, r)

    if gcd(m, r) != 1
        println("m and r are not coprime")
        nothing
    end

    return acos((3*m^2 + 3*m*r + r^2/2)/(3*m^2 + 3*m*r + r^2))
end

"""
    tblg_commensurate_basis(m, r, a0)

Returns the Bravais basis for a commensurate (m, r) twisted bilayer graphene structure. `a0` is the lattice parameter of an isolated graphene layer.
"""
function tblg_commensurate_basis(m, r, a0)

    a1 = a0*SVector(0.5, sqrt(3)/2)
    a2 = a0*SVector(-0.5, sqrt(3)/2)

    if gcd(r, 3) == 1
        t1 = (m)*a1 + (m+r)*a2
        t2 = -(m+r)*a1 + (2*m+r)*a2

    elseif gcd(r, 3) == 3
        t1 = (m+div(r,3))*a1 + div(r, 3)*a2
        t2 = -div(r, 3)*a1 + (m+2*div(r, 3))*a2
    end

    return t1, t2

end

"""
    tblg_commensurate_number_atoms(m, r)

Number of total carbon atoms in a unit cell of a commensurate (m, r) twisted bilayer graphene structure.
"""
function tblg_commensurate_number_atoms(m, r)

    if gcd(r, 3) == 1
        return 4*(3*m^2 + 3*m*r + r^2)
    elseif gcd(r, 3) == 3
        return 4*(m^2 + m*r + div(r^2, 3))
    end

end

"""
    tblg_commensurate_carbon_positions(m, r, a0)

Determines the position of the carbon atoms inside a unit cell of a commensurate (m, r) twisted bilayer graphene structure. `a0` is the lattice parameter of a single layer graphene. The functions returns as a tuple:
θ, Natls, WZ, [A1, A2], [a1bot, a2bot], [a1top, a2top], [listAbot, listBbot, listAtop, listBtop]
where:
- θ is the commensurate twist angle,
- Natls is number of carbon atoms per layer, per sublattice in the unit cell,
- WZ is a vector of positions, deliminting the Wigner-Seitz unit cell of the tBLG structure,
- [A1, A2] is the Bravais basis of the commensurate structure
- [a1bot, a2bot] is the Bravais basis of the isolated bottom layer
- [a1top, a2top] is the Bravais basis of the isolated top layer
- [listAbot, listBbot, listAtop, listBtop] constains the position of the atoms in each layer and sublattice.
"""
function tblg_commensurate_carbon_positions(m::Int, r::Int, a0)

    # @assert gcd(m, r) == 1 "m and r must be coprimes"

    θ = tblg_commensurate_angle(m, r)
    Natls =div(tblg_commensurate_number_atoms(m, r), 4)

    a1 = a0*SVector(0.5, sqrt(3)/2)
    a2 = a0*SVector(-0.5, sqrt(3)/2)

    A1, A2 = tblg_commensurate_basis(m, r, a0)

    α = atan(A1[2]+A2[2], A1[1]+A2[1])

    #rotate basis vectors
    A1 = R2d(pi/2-α)*A1
    A2 = R2d(pi/2-α)*A2

    # WZ unit cell
    C1 = (A1 + A2)/3
    WZ = [R2d(2*pi*n/6)*C1 for n=0:5]

    # bottom layer basis vectors
    a1bot = R2d(pi/2-α)*a1
    a2bot = R2d(pi/2-α)*a2
    τAbot = R2d(pi/2-α)*SVector(0.0, 0.0)
    τBbot = R2d(pi/2-α)*SVector(0.0, -a0/sqrt(3))

    # top layer basis vectors
    a1top = R2d(θ + pi/2 - α)*a1
    a2top = R2d(θ + pi/2 - α)*a2
    τAtop = R2d(θ + pi/2-α)*SVector(0.0, 0.0)
    τBtop = R2d(θ + pi/2-α)*SVector(0.0, a0/sqrt(3))

    N = ceil(Int, norm(A1)/a0)
    listAbot = SVector{2, Float64}[]
    listBbot = SVector{2, Float64}[]
    listAtop = SVector{2, Float64}[]
    listBtop = SVector{2, Float64}[]

    for n=-N:N, m=-N:N
        pbot = n*a1bot + m*a2bot
        if isinsideconvex(pbot, WZ)
            push!(listAbot, pbot + τAbot)
            push!(listBbot, pbot + τBbot)
        end

        ptop = n*a1top + m*a2top
        if isinsideconvex(ptop, WZ)
            push!(listAtop, ptop + τAtop)
            push!(listBtop, ptop + τBtop)
        end

    end

    # test we captured all the points
    if length(listAbot) != Natls
        println("Warning: Not all A atoms from bottom layer were found")
    end
    if length(listBbot) != Natls
        println("Warning: Not all B atoms from bottom layer were found")
    end
    if length(listAtop) != Natls
        println("Warning: Not all A atoms from top layer were found")
    end
    if length(listBtop) != Natls
        println("Warning: Not all B atoms from top layer were found")
    end

    return θ, Natls, WZ, [A1, A2], [a1bot, a2bot], [a1top, a2top], [listAbot, listBbot, listAtop, listBtop]

end

"""
    Vinter_tblg_sk(x; Vppp0 = -2.7, Vpps0 = 0.48, a0 = 2.46, d = 3.35, r0 = 0.45264)

Slater-Koster parametrization of the interlayer hopping in twisted bilayer graphene as a function of the in-plane separation `x` (a vector).
"""
function Vinter_tblg_sk(x; Vppp0 = -2.7, Vpps0 = 0.48, a0 = 2.46, d = 3.35, r0 = 0.45264)

    r = norm(x)

    return Vpps0*d^2/(r^2+d^2)*exp(-(sqrt(r^2+d^2)-d)/r0) + Vppp0*r^2/(r^2+d^2)*exp(-(sqrt(r^2+d^2)-a0/sqrt(3))/r0)

end

"""
    tblg_commensurate_build(m, r, a0, d, t, hoprange, Vinter::Function=Vinter_tblg_sk)

Builds a TBHamiltonian system for a commensurate (m, r) twisted bilayer graphene structure. The remaining arguments are:
- `a0` is the lattice parameter of single layer graphene;
- `d` is the interlayer seperation;
- `t` is the nearest-neighbour in-plane hopping;
- `hoprange` is a distance cutoff for the interlayer hoppings;
- `Vinter` is a function parametrizing the interlayer hopping as a function of the in-plane seperation.

It resturns: tblg, theta, Napl, WZ, abot, atop; where:
- `tblg`, is the TBHamiltonian struct;
- `theta`, is the twist angle (in radians);
- `Napl`, is the number of atoms per layer and per sublattice in the unit cell;
- `WZ`, are the vertices of the Wigner-Seitz unit cell;
- `abot`, is the Bravais basis of the botton layer;
- `atop`, is the Bravais basis of the top layer.
"""
function tblg_commensurate_build(m, r, a0, d, t, hoprange, Vinter::Function=Vinter_tblg_sk)

    theta, Napl, WZ, A, abot, atop, carbon = tblg_commensurate_carbon_positions(m, r, a0)

    tblg = TBHamiltonian(A[1], A[2])

    # Add orbtials
    ## A top atoms
    for atom in carbon[3]
        add_orbital!(tblg, SVector(atom..., d), :pz)
    end

    ## B top atoms
    for atom in carbon[4]
        add_orbital!(tblg, SVector(atom..., d), :pz)
    end

    ## A bot atoms
    for atom in carbon[1]
        add_orbital!(tblg, SVector(atom..., 0.0), :pz)
    end

    ## B bot atoms
    for atom in carbon[2]
        add_orbital!(tblg, SVector(atom..., 0.0), :pz)
    end

    # Add intralayer hoppings
    ## nearest-neighbours unit cells
    NN = [(0, 0), (1, 0), (0, 1), (-1, 1), (-1, 0), (0, -1), (1, -1)]

    ## create trees for NN atom search
    topR = zeros(2, 2*Napl)
    for i=1:Napl
        topR[:, i] = carbon[3][i]
        topR[:, i + Napl] = carbon[4][i]
    end
    toptree = BallTree(topR)

    botR = zeros(2, 2*Napl)
    for i=1:Napl
        botR[:, i] = carbon[1][i]
        botR[:, i + Napl] = carbon[2][i]
    end
    bottree = BallTree(botR)

    ## add intralayer hopping in top layer
    for (n, m) in NN
        for i=1:2*Napl
            idxs = inrange(toptree, n*A[1] + m*A[2] + topR[:, i], 1.1*a0/sqrt(3), true)
            for j in idxs
                if norm(n*A[1] + m*A[2] + topR[:, i] - topR[:, j])>0.1*a0/sqrt(3)
                    add_hopping!(tblg, j, i, [n, m], -t)
                end
            end
        end
    end

    ## add intralayer hopping in bottom layer
    for (n, m) in NN
        for i=1:2*Napl
            idxs = inrange(bottree, n*A[1] + m*A[2] + botR[:, i], 1.1*a0/sqrt(3), true)
            for j in idxs
                if norm(n*A[1] + m*A[2] + botR[:, i] - botR[:, j])>0.1*a0/sqrt(3)
                    add_hopping!(tblg, 2*Napl + j, 2*Napl + i, [n, m], -t)
                end
            end
        end
    end

    ## add interlayer hoppings
    for (n, m) in NN
        ### from bottom to top
        for i=1:2*Napl
            idxs = inrange(bottree, n*A[1] + m*A[2] + topR[:, i], hoprange, true)
            for j in idxs
                add_hopping!(tblg, 2*Napl + j, i, [n, m], Vinter(n*A[1] + m*A[2] + topR[:, i] - botR[:, j]))
            end
        end

        ### from top to bottom
        for i=1:2*Napl
            idxs = inrange(toptree, n*A[1] + m*A[2] + botR[:, i], hoprange, true)
            for j in idxs
                add_hopping!(tblg, j, 2*Napl + i, [n, m], Vinter(n*A[1] + m*A[2] + botR[:, i] - topR[:, j]))
            end
        end
    end

    sort_hoppings!(tblg)

    return tblg, theta, Napl, WZ, abot, atop

end


text:
import .Makie
import GeometryBasics
export visualize, visualize!

"""
    visualize(pcloud::PointCloud; kwargs...)

Visualize PointCloud `pcloud` at `index`.

Dimension of points in PointCloud `pcloud` must be 3.

### Optional Arguments:
- color (Symbol)       - Color of the marker, default `:blue`
- markersize (Number)  - Size of the marker, default `npoints(pcloud)/5000`

"""
function visualize(p::PointCloud, index::Number = 1; kwargs...)
    points = cpu(p[index])
    size(points, 1) == 3 || error("dimension of points in PointCloud must be 3.")

    kwargs = convert(Dict{Symbol,Any}, kwargs)
    get!(kwargs, :color, :lightgreen)
    get!(kwargs, :markersize, npoints(p) / 5000)

    Makie.meshscatter(points[1, :], points[2, :], points[3, :]; kwargs...)
end

function visualize!(axis3::Makie.Axis3, p::PointCloud, index::Number = 1; kwargs...)
    points = cpu(p[index])
    size(points, 1) == 3 || error("dimension of points in PointCloud must be 3.")

    kwargs = convert(Dict{Symbol,Any}, kwargs)
    get!(kwargs, :color, :lightgreen)
    get!(kwargs, :markersize, npoints(p) / 5000)

    Makie.meshscatter!(axis3, points[1, :], points[2, :], points[3, :]; kwargs...)
end

"""
    visualize(m::TriMesh, index::Int=1; kwargs...)

Visualize mesh at `index` in TriMesh `m`.

### Optional Arguments:
- color (Symbol)       - Color of the marker, default `:red`

"""
function visualize(m::GeometryBasics.Mesh; kwargs...) where {T,R}
    kwargs = convert(Dict{Symbol,Any}, kwargs)
    get!(kwargs, :color, :orange)

    Makie.mesh(GeometryBasics.normal_mesh(m); kwargs...)
end

function visualize!(axis3::Makie.Axis3, m::GeometryBasics.Mesh; kwargs...) where {T,R}
    plt = visualize(m; kwargs...)
    Makie.mesh!(axis3, plt.plot.input_args[1].val; plt.plot.attributes...)
end

visualize(m::TriMesh, index::Int = 1; kwargs...) = visualize(GBMesh(m, index); kwargs...)

"""
     visualize!(axis3::Makie.Axis3, object, args...; kwargs...)
Similar to `Flux3D.visualize` except it accepts `axis3::Makie.Axis3` and update it by rendering `object`.
See also [`visualize`](@ref).
"""
function visualize!(axis3::Makie.Axis3, m::TriMesh, args...; kwargs...)
    plt = visualize(GBMesh(m, args...); kwargs...)
    Makie.mesh!(axis3, plt.plot.input_args[1].val; plt.plot.attributes...)
end

"""
    visualize(v::VoxelGrid, index::Int=1; kwargs...)

Visualize voxel at `index` in VoxelGrid `v`.

### Optional Arguments:
- color (Symbol)       - Color of the marker, default `:red`

"""
function visualize(
    v::VoxelGrid,
    index::Int = 1,
    thresh::Number = 0.49f0;
    algo = :Exact,
    kwargs...,
)
    algo in [:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets] ||
        error("given algo: $(algo) is not supported. Accepted algo are
              {:Exact,:MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets}.")
    kwargs = convert(Dict{Symbol,Any}, kwargs)
    get!(kwargs, :color, :violet)
    method = algo == :Exact ? _voxel_exact : _voxel_algo
    v, f = method(cpu(v[index]), Float32(thresh), algo)

    m = GBMesh(v, f)
    Makie.mesh(GeometryBasics.normal_mesh(m); kwargs...)
end

function visualize!(
    axis3::Makie.Axis3,
    v::VoxelGrid,
    args...;
    kwargs...,
)
    plt = visualize(v, args..., ; kwargs...)
    Makie.mesh!(axis3, plt.plot.input_args[1].val; plt.plot.attributes...)
end

visualize(v::Dataset.AbstractDataPoint; kwargs...) = visualize(v.data; kwargs...)

visualize(v::AbstractCustomObject; kwargs...) =
    error("Define visualize function for custom type: $(typeof(v)).
            Use `import Flux3D.visualize` and define function
            `visualize(v::$(typeof(v)); kwargs...)`")

ext: jl


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.DAE (ModiaMath/DAE/_module.jl)
#

# Utility functions used below
defaultVariableName(model::Any, vcat::VariableCategory, vindex::Int) = vcat == Category_X ? "x[" * string(vindex) * "]" :
                                                                      (vcat == Category_W ? "w[" * string(vindex) * "]" :
                                                                                "der(x[" * string(vindex) * "])" )

#nameVector(name, nNames::Int) = [Symbol(name, "[", i, "]") for i = 1:nNames]
#fcNameVector(fc, names::AbstractVector) = [ Symbol(fc, "(", names[i], ")") for i in eachindex(names)]

nameVector(name, nNames::Int) = [string(name, "[", i, "]") for i = 1:nNames]
fcNameVector(fc, names::AbstractVector) = [ string(fc, "(", names[i], ")") for i in eachindex(names)]

function getVariableName(model::Any, vcat::VariableCategory, vindex::Int, nx=0;
                         xNames::Vector{String}=nameVector("x", nx),
                         derxNames::Vector{String}=fcNameVector("der", xNames),
                         wNames::Vector{String}=String[])
    if vcat == Category_X
        return xNames[vindex]
    elseif vcat == Category_DERX
        return derxNames[vindex]
    elseif vcat == Category_W
        return wNames[vindex]
    else
        error("Wrong argument vcat = ", vcat)
    end
end

function getResultNames(model::Any)
    sim::SimulationState = model.simulationState
    resultNames = ["time"; String[sim.getVariableName(model, Category_X, i)  for i = 1:sim.nx];
                          String[sim.getVariableName(sim, Category_DERX, i) for i = 1:sim.nx];
                          String[sim.getVariableName(sim, Category_W, i)    for i = 1:sim.nw] ]
    return resultNames
end

mutable struct SimulationState
    name::Symbol                      # Name of model
    model::Any                        # Model

    getModelResidues!::Function
    getVariableName::Function
    getResultNames::Function
    storeResult!::Function
    getResult::Function
    eventHandler::EventHandler

    nx::Int   # Length of x-vector
    nd::Int   # fd(der(x),x,t) = 0; nd = length(fd)
    nc::Int   #        fc(x,t) = 0; nc = length(fc)
    nw::Int   # Number of auxiliary variables

    # Other model information
    nz::Int                           # Number of event indicators
    zDir::Vector{Cint}                # zDir[i] =  0: Root is reported for both crossing directions
                                     #         =  1: Root is reported when crossing from negative to positive direction
                                     #         = -1: Root is reported when crossing from positive to negative direction
    sparse::Bool                      # = true, if (sparse) jac present and shall be used
    # jac::Union{SparseMatrixCSC{Float64,Cint},Void}
                                     # Optional sparse Jacobian of DAE: der(f!,y) + cr*der(f!, yp)
                                     # (cr is a constant provided by the integrator)
    @static if VERSION >= v"0.7.0-DEV.2005"
        jac::Nothing
        #cg::Union{ModiaMath.SparseJacobian.ColumnGroups,Void} # if sparse, column groups of Jacobian jac
        cg::Nothing # if sparse, column groups of Jacobian jac
    else
        jac::Void
        #cg::Union{ModiaMath.SparseJacobian.ColumnGroups,Void} # if sparse, column groups of Jacobian jac
        cg::Void # if sparse, column groups of Jacobian jac

    end
    # Model specific information not used by ModiaMath (can be used as default setting for ModiaMath)
    defaultTolerance::Float64    # default relative integration tolerance
    defaultStartTime::Float64    # default start time
    defaultStopTime::Float64     # default stop time
    defaultInterval::Float64     # default integration interval

    # Run-time information
    time::ModiaMath.Time                        # Actual simulation time
    logger::ModiaMath.Logger                    # Logger object
    statistics::ModiaMath.SimulationStatistics  # Statistics object (complete after a full simulation run)
    tolerance::Float64                          # Telative integration tolerance
    startTime::Float64                          # Start time of the simulation
    stopTime::Float64                           # Stop time of the simulation
    interval::Float64                           # Interval of the simulation

    x_start::Vector{Float64}
    x_fixed::Vector{Bool}
    x_nominal::Vector{Float64}
    x_errorControl::Vector{Bool}

    w::Vector{Float64}  # length(w) = nw

    # Auxiliary storage needed during initialization and at events
    eqInfo::ModiaMath.NonlinearEquationsInfo
    nonlinearEquationsMode::Int
    xev_old::Vector{Float64}
    xev_beforeEvent::Vector{Float64}
    xev::Vector{Float64}
    yNonlinearSolver::Vector{Float64}
    derxev::Vector{Float64}
    residues::Vector{Float64}
    rScale::Vector{Float64}
    tev::Float64
    hev::Float64
    FTOL::Float64
    scaleConstraintsAtEvents::Bool
    initialization::Bool
    use_x_fixed::Bool

    # Information updated during simulation
    storeResult::Bool             # = true, if DAE is called to store results of variables
    rawResult::ModiaMath.RawResult
    result::Any

    function SimulationState(name,
                            getModelResidues!::Function,
                            x_start::Vector{Float64},
                            getVariableName::Function=defaultVariableName;
                            nc=0,
                            nz=0,
                            nw=0,
                            zDir::Vector{Int}=fill(0, nz),
                            x_fixed::Vector{Bool}=fill(false, length(x_start)),
                            x_errorControl::Vector{Bool}=fill(true, length(x_start)),
                            hev=1e-8,
                            scaleConstraintsAtEvents::Bool=true,
                            jac=nothing,
                            maxSparsity::Float64=0.1,
                            getResultNames::Function=ModiaMath.getResultNames,
                            storeResult!::Function=ModiaMath.storeRawResult!,
                            getResult::Function=ModiaMath.getStringDictResult,
                            defaultTolerance=1e-4,
                            defaultStartTime=0.0,
                            defaultStopTime=1.0,
                            defaultInterval=(defaultStopTime - defaultStartTime) / 500.0)

        # Check input arguments
        @assert(nc >= 0)
        @assert(nc <= length(x_start))
        @assert(nz >= 0)
        @assert(nw >= 0)
        @assert(length(x_fixed) == length(x_start))
        @assert(length(zDir) == nz)
        @assert(hev > 0.0)
        @assert(length(x_errorControl) == length(x_start))
        nx = length(x_start)

        @static if VERSION >= v"0.7.0-DEV.2005"
            if typeof(jac) != Nothing
                @assert(size(jac, 1) == nx)
                @assert(size(jac, 2) == nx)
            end
        else
            if typeof(jac) != Void
                @assert(size(jac, 1) == nx)
                @assert(size(jac, 2) == nx)
            end
        end

        @assert(0.0 <= maxSparsity <= 1.0)
        @assert(defaultTolerance > 0.0)
        @assert(defaultStopTime >= defaultStartTime)
        @assert(defaultInterval > 0.0)

        # Compute utility elements
        nd = nx - nc
        eventHandler = EventHandler(nz=nz)

        eqInfo = ModiaMath.NonlinearEquationsInfo("???", nx, getEventResidues!)

        # If jac is provided and sparse enough, use sparse methods
        #sparse = typeof(jac) == Void ? false : nnz(jac)/(nx*nx) < maxSparsity
        #jac    = sparse ? convert(SparseMatrixCSC{Float64,Cint}, jac) : nothing
        #cg     = sparse ? ColumnGroups(jac) : nothing
        #if sparse
        #   @assert( nnz(jac) >= nx )   # otherwise jac is singular
        #end
        sparse = false
        jac = nothing
        cg = nothing

        new(Symbol(name), nothing, getModelResidues!, getVariableName, getResultNames,
            storeResult!, getResult, eventHandler, nx, nd, nc, nw, nz, zDir,
            sparse, jac, cg, defaultTolerance, defaultStartTime, defaultStopTime,
            defaultInterval, NaN, ModiaMath.Logger(),
            ModiaMath.SimulationStatistics(nx, sparse, sparse ? cg.ngroups : 0),
            NaN, NaN, NaN, NaN,
            x_start, x_fixed, ones(nx), x_errorControl, zeros(nw),
            eqInfo, 0, zeros(nx), zeros(nx), zeros(nx),
            zeros(nx), zeros(nx), zeros(nx), ones(nx), 0.0, hev, 0.0,
            scaleConstraintsAtEvents, false, false, false)
    end
end

"""
    result = isNearlyEqual(x1,x2)

The function returns true if x1 and x2 are nearly equal
"""
isNearlyEqual(x1::Float64, x2::Float64) = abs(x1) > 1e-8 ? abs(x1 - x2) / max(abs(x1)) <= 1e-8 : abs(x1 - x2) <= 1e-8

function getEventResidues!(eqInfo::ModiaMath.NonlinearEquationsInfo, y::Vector{Float64}, r::Vector{Float64})
    model = eqInfo.extraInfo
    sim::SimulationState = model.simulationState
    @assert(length(y) == eqInfo.ny)
    @assert(length(r) == eqInfo.ny)
    @assert(length(y) == sim.nx)

    residues = sim.residues
    hev      = sim.hev

   # Copy unknowns (y) to xev and derxev
    for i in eachindex(y)
        if sim.use_x_fixed && sim.x_fixed[i]
         # Unknowns are at the left limit
            sim.xev[i]    = sim.xev_old[i]
            sim.derxev[i] = (sim.xev_old[i] - y[i]) / hev
        else
         # Unknowns are at the right limit
            sim.xev[i]    = y[i]
            sim.derxev[i] = (y[i] - sim.xev_old[i]) / hev
        end
    end

   # Compute residues
    sim.time = sim.tev
    Base.invokelatest(sim.getModelResidues!, model, sim.tev, sim.xev, sim.derxev, residues, sim.w)

   # Copy to eq-residues and scale with h
    if sim.scaleConstraintsAtEvents
        for i = 1:sim.nd
            r[i] = residues[i]
        end

        for i = sim.nd + 1:sim.nx
            r[i] = residues[i] / hev
        end
    else
        for i = 1:sim.nx
            r[i] = residues[i]
        end
    end
    eqInfo.lastNorm_r = norm(r, Inf)
    eqInfo.lastrScaledNorm_r = norm(sim.rScale .* r, Inf)

    # println("maxabs(r) = ", maxabs(r))
    # println("residue = "),display(r)
    # println("hev = ", hev, ", x_old = ", sim.xev_old[1], ", x = ", y[1], ", der(x) = ",sim.derxev[1], ", r = ", r[1])

    return nothing
end

function reinitialize!(model, sim::SimulationState, tev::Float64, xevIsConsistent::Bool)
    nx = sim.nx
    nd = sim.nd
    nc = sim.nc

    # Determine consistent initial conditions
    for i in eachindex(sim.xev)
        sim.xev_old[i] = sim.xev[i]
    end

    if !xevIsConsistent
        # xev is potentially not consistent; determine consistent xev
        if ModiaMath.isLogEvents(sim)
            println("        determine consistent DAE variables x (with implicit Euler step; step size = ", sim.hev, ")")
            if sim.initialization
                println("            (fixed attribute defines whether unknowns are left or right limit values)")
            end
        end
        sim.tev = tev
        sim.eqInfo.lastNorm_r        = 1.0
        sim.eqInfo.lastrScaledNorm_r = 1.0

        for i in eachindex(sim.yNonlinearSolver)
            sim.yNonlinearSolver[i] = sim.xev[i]
        end

        if sim.initialization
            sim.use_x_fixed = true
        end

        sim.nonlinearEquationsMode = 1
        ModiaMath.solveNonlinearEquations!(sim.eqInfo, sim.yNonlinearSolver; rScale=sim.rScale, FTOL=sim.FTOL)
        sim.use_x_fixed = false

        # Copy result to sim.xev
        for i = 1:nx
            sim.xev[i] = sim.yNonlinearSolver[i]
        end

        # Only during initialization: Copy fixed = true values
        if sim.initialization
            for i = 1:nx
                if sim.x_fixed[i]
                    sim.xev[i] = sim.xev_old[i]
                end
            end
        end

        # Copy xev to xev_old
        for i = 1:nx
            sim.xev_old[i] = sim.xev[i]
        end
    end

    # Print log message
    if ModiaMath.isLogInfos(sim)
        for i = 1:nx
            if !isNearlyEqual(sim.xev_beforeEvent[i], sim.xev[i])
                xname = sim.getVariableName(model, Category_X, i)
                println("            ", xname, " = ", sim.xev_beforeEvent[i], " changed to ", sim.xev[i])
            end
        end
    end

   # Determine consistent derxev
    sim.tev = tev
    sim.eqInfo.lastNorm_r        = 1.0
    sim.eqInfo.lastrScaledNorm_r = 1.0
    for i in eachindex(sim.yNonlinearSolver)
        sim.yNonlinearSolver[i] = sim.xev[i]
    end

    sim.use_x_fixed = false
    if ModiaMath.isLogEvents(sim)
        println("        determine consistent DAE variables der(x) (with implicit Euler step; step size = ", sim.hev, ")")
    end

    sim.nonlinearEquationsMode = 2
    ModiaMath.solveNonlinearEquations!(sim.eqInfo, sim.yNonlinearSolver; rScale=sim.rScale, FTOL=sim.FTOL)
    sim.nonlinearEquationsMode = 0
end

function eventIteration!(model, sim::SimulationState, tev::Float64)
    eh::EventHandler = sim.eventHandler
    xevIsConsistent::Bool   = false

    # Initialize event iteration
    initEventIteration!(eh, tev)

    # Perform event iteration
    while true
        # Determine event branches
        for i = 1:sim.nx
            sim.xev_beforeEvent[i] = sim.xev[i]
        end
        eh.event = true
        Base.invokelatest(sim.getModelResidues!, model, tev, sim.xev, sim.derxev, sim.residues, sim.w)

        if terminateEventIteration!(eh)
            eh.event = false
            break
        end

        if eh.initial
            eh.initial = false
            xevIsConsistent = sim.nc == 0

            if sim.nc > 0
                if norm(sim.residues[sim.nd + 1:end], Inf) < 1e-3
                    xevIsConsistent = true
                end
            end
        else
            xevIsConsistent = true
            for i = 1:sim.nx
                if sim.xev[i] != sim.xev_beforeEvent[i]
                    xevIsConsistent = false
                    break
                end
            end
        end

        # Fix event branches and determine new consistent sim.xev_start, sim.derxev_start
        eh.event = false
        reinitialize!(model, sim, tev, xevIsConsistent)
    end
end

function initialize!(model, sim::SimulationState, t0::Float64, nt::Int, tolerance::Float64)
    eh::EventHandler = sim.eventHandler
    eh.zDir    = sim.zDir
    eh.logger  = sim.logger
    eh.initial = true
    eh.afterSimulationStart = false
    sim.eqInfo.extraInfo = model
    sim.rawResult = ModiaMath.RawResult(nt, sim.getResultNames(model))
    sim.model     = model
    sim.FTOL      =  eps(Float64)^(1 / 3)   # residue tolerance for nonlinear solver
    # println("... FTOL = ", sim.FTOL)
    sim.initialization = true

    # Initialize auxiliary arrays for event iteration
    # println("... initialize!: scaleConstraintsAtEvents = ", sim.scaleConstraintsAtEvents, ", nc = ", sim.nc, ", nd = ", sim.nd)
    for i in eachindex(sim.xev),
      sim.xev[i]    = sim.x_start[i]
        sim.derxev[i] = 0.0
    end

    if sim.scaleConstraintsAtEvents
        for i = 1:sim.nd
            sim.rScale[i] = 1.0
        end
        for i = sim.nd + 1:sim.nx
            sim.rScale[i] = 1.0 / sim.hev
        end
    else
        for i = 1:sim.nx
            sim.rScale[i] = 1.0
        end
    end
    # println("... initialize!: rScale = ", sim.rScale)

    # Perform initial event iteration
    eventIteration!(model, sim, t0)
    sim.initialization = false

    # Compute auxiliary variables w and store all results
    computeAndStoreResult!(model, sim, t0, sim.xev, sim.derxev)
    eh.afterSimulationStart = true

    return InitInfo(sim.xev, sim.derxev;
                   y_nominal=sim.x_nominal,
                   y_errorControl=sim.x_errorControl,
                   maxTime=eh.maxTime,
                   nextEventTime=eh.nextEventTime,
                   integrateToEvent=eh.integrateToEvent,
                   terminate=eh.restart == Terminate)
end

getResidues!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64}, residues::Vector{Float64}, hcur) =
   Base.invokelatest(sim.getModelResidues!, model, t, y, yp, residues, sim.w)

function computeAndStoreResult!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64})
    sim.storeResult = true
    sim.time        = t
    Base.invokelatest(sim.getModelResidues!, model, t, y, yp, sim.residues, sim.w)
    sim.storeResult = false
    sim.storeResult!(sim.rawResult, model, t, y, yp, sim.w)
    return nothing
end

function terminate!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64})
    eh::EventHandler = sim.eventHandler
    eh.terminal = true
    sim.time    = t
    Base.invokelatest(sim.getModelResidues!, model, t, y, yp, sim.residues, sim.w)
    eh.terminal = false
    #sim.result,nt = Result.getDictResult(sim.rawResult)
    sim.result, nt  = sim.getResult(model, sim.rawResult)
    ModiaMath.set_nResultsForSimulationStatistics!(sim.statistics, nt)
    return sim.result
end

function processEvent!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64}, eventInfo::EventInfo)
    eh::EventHandler = sim.eventHandler

    # Event iteration
    for i in eachindex(sim.xev)
        sim.xev[i]    = y[i]
        sim.derxev[i] = yp[i]
    end
    eventIteration!(model, sim, t)
    for i in eachindex(sim.xev)
        y[i]  = sim.xev[i]
        yp[i] = sim.derxev[i]
    end

    # Compute auxiliary variables w and store all results
    computeAndStoreResult!(model, sim, t, y, yp)

    eventInfo.restart          = eh.restart
    eventInfo.maxTime          = eh.maxTime
    eventInfo.nextEventTime    = eh.nextEventTime
    eventInfo.integrateToEvent = eh.integrateToEvent

    return nothing
end

function getEventIndicators!(model, sim::SimulationState, t::Float64, y::Vector{Float64}, yp::Vector{Float64}, z::Vector{Float64})
    eh::EventHandler = sim.eventHandler
    eh.crossing = true
    sim.time    = t
    Base.invokelatest(sim.getModelResidues!, model, t, y, yp, sim.residues, sim.w)
    eh.crossing = false

    for i = 1:eh.nz
        z[i] = eh.z[i]
    end
    return nothing
end


text:
# function schecklt(a::Int, b1::Int, b2::Int, s::Int)
#     a < b && cld(b1, s) == cld(b2, s)
# end

cilt(a::S, b::S) where {S <: CartesianIndex} = all(a.I .<= b.I)

function cartstrides(A::AbstractArray)
    s = size(A)
    dims = length(s)
    inds = []
    for i = 1:dims
        z = zeros(Int, dims)
        z[i] = 1
        push!(inds, CartesianIndex(z...))
    end
    [inds...]
end

function set_neighbours!(xn::Vector{S}, x::S, tags::Array{UInt8, N}, cs::Vector{S}, I1::S, Iend::S) where {N, S <: CartesianIndex{N}}
    for (i,s) in enumerate(cs)
        if cilt(I1, x-s) && @inbounds tags[x-s] < 0x3
            @inbounds xn[2*i-1] = x-s
        else
            @inbounds xn[2*i-1] = x
        end
        if cilt(x+s, Iend) && @inbounds tags[x+s] < 0x3
            @inbounds xn[2*i] = x+s
        else
            @inbounds xn[2*i] = x
        end
    end
end

@noinline function check1bwd(τ1::Array{R, N}, τ0::Array{R, N}, tags::Array{UInt8, N}, x::S, s::S, I1::S) where {N, R <: AbstractFloat, S <: CartesianIndex{N}}
    if cilt(I1, x-s) && @inbounds tags[x-s] == 0x3
        return true
    else
        return false
    end
end

@noinline function check1fwd(τ1::Array{R, N}, τ0::Array{R, N}, tags::Array{UInt8, N}, x::S, s::S, Iend::S) where {N, R <: AbstractFloat, S <: CartesianIndex{N}}
    if cilt(x+s, Iend) && @inbounds tags[x+s] == 0x3
        return true
    else
        return false
    end
end

@noinline function check2bwd(τ1::Array{R, N}, τ0::Array{R, N}, tags::Array{UInt8, N}, x::S, s::S, I1::S) where {N, R <: AbstractFloat, S <:CartesianIndex{N}}
    if cilt(I1, x-2*s) && @inbounds tags[x-2*s] == 0x3 && @inbounds τ1[x-s]*τ0[x-s] >= τ1[x-2*s]*τ0[x-2*s]
        return true
    else
        return false
    end
end

@noinline function check2fwd(τ1::Array{R, N}, τ0::Array{R, N}, tags::Array{UInt8, N}, x::S, s::S, Iend::S) where {N, R <: AbstractFloat, S <:CartesianIndex{N}}
    if cilt(x+2*s, Iend) && @inbounds tags[x+2*s] == 0x3 && @inbounds τ1[x+s]*τ0[x+s] > τ1[x+2*s]*τ0[x+2*s]
        return true
    else
        return false
    end
end


text:
const CONT = 0
const DISC = 1

struct CDConformity{D} <: Conformity
  cont::NTuple{D,Int}
end

function Conformity(reffe::GenericLagrangianRefFE{<:CDConformity},sym::Symbol)
  if sym == :L2
    L2Conformity()
  else
    @unreachable """\n
    It is not possible to use conformity = $sym on a LagrangianRefFE with CD conformity.

    Either use conformity = :L2 or pass an instance of CDConformity.
    """
  end
end

function get_face_own_nodes(reffe::GenericLagrangianRefFE{<:CDConformity},conf::CDConformity)
  _cd_get_face_own_nodes(reffe,conf)
end

function get_face_own_nodes(reffe::GenericLagrangianRefFE{GradConformity},conf::CDConformity)
  _cd_get_face_own_nodes(reffe,conf)
end

# Constructors

function _CDLagrangianRefFE(::Type{T},p::ExtrusionPolytope{D},order::Int,cont) where {T,D}
  orders = tfill(order,Val{D}())
  _CDLagrangianRefFE(T,p,orders,cont)
end

function _CDLagrangianRefFE(::Type{T},p::ExtrusionPolytope{D},orders,cont) where {T,D}
  cond(c,o) = ( o > 0 || c == DISC )
  @assert all((cond(cont[k],orders[k]) for k in 1:length(orders)))
  _cd_lagrangian_ref_fe(T,p,orders,cont)
end

function _cd_lagrangian_ref_fe(::Type{T},p::ExtrusionPolytope{D},orders,cont) where {T,D}

  @assert isa(p,ExtrusionPolytope)

  prebasis = compute_monomial_basis(T,p,orders)

  nodes, face_own_nodes = cd_compute_nodes(p,orders)
  dofs = LagrangianDofBasis(T,nodes)

  nnodes = length(dofs.nodes)
  ndofs = length(dofs.dof_to_node)

  face_own_nodes = _compute_cd_face_own_nodes(p,orders,cont)
  face_nodes = _compute_face_nodes(p,face_own_nodes)

  face_own_dofs = _generate_face_own_dofs(face_own_nodes, dofs.node_and_comp_to_dof)
  face_dofs = _compute_face_nodes(p,face_own_dofs)

  data = nothing

  conf = CDConformity(Tuple(cont))

  reffe = GenericRefFE{typeof(conf)}(
      ndofs,
      p,
      prebasis,
      dofs,
      conf,
      data,
      face_dofs)

  GenericLagrangianRefFE(reffe,face_nodes)
end

function _cd_get_face_own_nodes(reffe,conf::CDConformity)
  p = get_polytope(reffe)
  orders = get_orders(get_prebasis(reffe))
  cont = conf.cont
  cond(c,o) = ( o > 0 || c == DISC )
  @assert all((cond(cont[k],orders[k]) for k in 1:length(orders)))
  p = get_polytope(reffe)
  dofs = get_dof_basis(reffe)
  @assert is_n_cube(p)
  face_own_nodes = _compute_cd_face_own_nodes(p,orders,cont)
  face_own_nodes
end

function _cd_get_face_own_dofs(reffe,conf::CDConformity)
  p = get_polytope(reffe)
  orders = get_orders(get_prebasis(reffe))
  cont = conf.cont
  cond(c,o) = ( o > 0 || c == DISC )
  @assert all((cond(cont[k],orders[k]) for k in 1:length(orders)))
  p = get_polytope(reffe)
  dofs = get_dof_basis(reffe)
  @assert is_n_cube(p)
  face_own_nodes = _compute_cd_face_own_nodes(p,orders,cont)
  face_own_dofs = _generate_face_own_dofs(face_own_nodes, dofs.node_and_comp_to_dof)
  face_own_dofs
end

function cd_compute_nodes(p::Polytope{D},orders) where D
  active_faces = _active_faces(p,orders)
  nodes = Point{D,Float64}[]
  facenodes = [Int[] for i in 1:num_faces(p)]
  _cd_compute_high_order_nodes_dim_0!(nodes,facenodes,p,active_faces)
  for d in 1:(num_dims(p)-1)
    _cd_compute_high_order_nodes_dim_d!(nodes,facenodes,p,orders,Val{d}(),active_faces)
  end
  if active_faces[end]
    _cd_compute_high_order_nodes_dim_D!(nodes,facenodes,p,orders)
  end

  nodes, facenodes = _move_nodes(p,orders,active_faces,nodes,facenodes)

  (nodes, facenodes)
end

function _cd_compute_high_order_nodes_dim_0!(nodes,facenodes,p,active_faces)
  x = get_vertex_coordinates(p)
  k = 1
  for vertex in 1:num_vertices(p)
    if active_faces[vertex]
      push!(nodes,x[vertex])
      push!(facenodes[vertex],k)
      k += 1
    end
  end
end

@noinline function _cd_compute_high_order_nodes_dim_d!(nodes,facenodes,p,orders,::Val{d},active_faces) where d
  x = get_vertex_coordinates(p)
  offset = get_offset(p,d)
  k = length(nodes)+1
  for iface in 1:num_faces(p,d)
    if active_faces[offset+iface]
      face = Polytope{d}(p,iface)
      face_ref_x = get_vertex_coordinates(face)
      face_prebasis = MonomialBasis(Float64,face,1)
      change = inv(evaluate(face_prebasis,face_ref_x))
      face_shapefuns = linear_combination(change,face_prebasis)
      face_vertex_ids = get_faces(p,d,0)[iface]
      face_x = x[face_vertex_ids]
      face_orders = compute_face_orders(p,face,iface,orders)
      face_interior_nodes = compute_own_nodes(face,face_orders)
      face_high_x = evaluate(face_shapefuns,face_interior_nodes)*face_x
      for xi in 1:length(face_high_x)
        push!(nodes,face_high_x[xi])
        push!(facenodes[iface+offset],k)
        k += 1
      end
    end
  end
end

function _cd_compute_high_order_nodes_dim_D!(nodes,facenodes,p,orders)
  k = length(nodes)+1
  p_high_x = compute_own_nodes(p,orders)
  for xi in 1:length(p_high_x)
    push!(nodes,p_high_x[xi])
    push!(facenodes[end],k)
    k += 1
  end
end

function _compute_cd_face_own_nodes(p::ExtrusionPolytope{D},orders::NTuple{D,<:Int},cont::NTuple{D,<:Int}) where D
  nodes, _ = cd_compute_nodes(p,orders)
  _ord = convert(NTuple{D,Float64},orders)
  _nodes = map(i -> i.data.*_ord,nodes)
  _nodes = map(i -> convert(NTuple{D,Int},i),_nodes)
  perm = Dict(_nodes[i] => i for i = 1:length(_nodes))
  face_owned_nodes = Vector{Int}[]
  for nf in p.dface.nfaces
    anc = nf.anchor
    ext = nf.extrusion
    fns = UnitRange{Int}[]
    for i in 1:length(orders)
      e,a,c,o = ext[i],anc[i],cont[i],orders[i]
      if e==0 && c == DISC
        push!(fns,0:-1)
        break
      elseif e==0 && c == CONT
        push!(fns,a*o:a*o)
      elseif e==1 && c == DISC
        push!(fns,0:o)
      elseif e==1 && c == CONT
        push!(fns,1:o-1)
      end
    end
    n_ijk = map(i->i.I,collect(CartesianIndices(Tuple(fns))))
    if length(n_ijk) == 0
      n_i = Int[]
    else
      n_i = sort(reshape([perm[i] for i in n_ijk],length(n_ijk)))
    end
    push!(face_owned_nodes,n_i)
  end
  face_owned_nodes
end

function _compute_face_nodes(p::ExtrusionPolytope,f_o)
  nfs = p.dface.nf_nfs
  f_n = [Int[] for i in 1:num_faces(p)]
  for i in 1:length(nfs)
    for nf in nfs[i]
      push!(f_n[i],f_o[nf]...)
    end
  end
  f_n
end

function _active_faces(p::Polytope,orders)
  nfs = p.dface.nfaces
  is_valid = Bool[]
  cond = (o,e,a) -> (o > 0 || (e == 0 && a == 0))
  for nf in nfs
    isv = all((cond(orders[k],nf.extrusion[k],nf.anchor[k]) for k in 1:length(orders)))
    push!(is_valid,isv)
  end
  is_valid
end

function _move_nodes(p::ExtrusionPolytope,orders,act_f,nodes,f_own)
  if any(map(i->i==0,orders))
    tr = map(i -> i == 0 ? 1 : 0,orders)
    nfs = p.dface.nfaces
    con(ext,anc) = f -> ( f.extrusion == ext && f.anchor == anc )
    perm = zeros(Int,length(nfs))
    for i in findall(act_f)
      nf  =nfs[i]
      ext = VectorValue(map(+,nf.extrusion,tr))
      anc = nf.anchor
      r = findall(con(ext,anc),nfs)
      @assert (length(r)<=1)
      perm[i] = r[1]
    end

    facenodes = [Int[] for i in 1:num_faces(p)]
    for f in findall(act_f)
      facenodes[perm[f]] = f_own[f]
    end

    tr = VectorValue(map(i -> i == 0 ? 0.5 : 0.0,orders))
    nodes = nodes.+tr

    return nodes, facenodes

  else
    return nodes, f_own
  end
end


text:
export LPSolve

function LPSolve()

 nothing
end


text:
using DistributionsAD

using Combinatorics
using Distributions
using FiniteDifferences
using PDMats
using Requires

# Figure out which AD backend to test
const AD = get(ENV, "AD", "All")
if AD == "All" || AD == "ForwardDiff"
    @eval using ForwardDiff
end
if AD == "All" || AD == "Zygote"
    @eval using Zygote
end
if AD == "All" || AD == "ReverseDiff"
    @eval using ReverseDiff
end
if AD == "All" || AD == "Tracker"
    @eval using Tracker
end

using Random, LinearAlgebra, Test

using Distributions: meanlogdet
using DistributionsAD: TuringUniform, TuringMvNormal, TuringMvLogNormal,
                       TuringPoissonBinomial
using StatsBase: entropy
using StatsFuns: binomlogpdf, logsumexp

const FDM = FiniteDifferences
const GROUP = get(ENV, "GROUP", "All")

# Create positive definite matrix
to_posdef(A::AbstractMatrix) = A * A' + I
to_posdef_diagonal(a::AbstractVector) = Diagonal(a.^2 .+ 1)

@require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
    # Define adjoints for Tracker
    to_posdef(A::TrackedMatrix) = Tracker.track(to_posdef, A)
    Tracker.@grad function to_posdef(A::TrackedMatrix)
        data_A = Tracker.data(A)
        S = data_A * data_A' + I
        function pullback(∇)
            return ((∇ + ∇') * data_A,)
        end
        return S, pullback
    end
end

if GROUP == "All" || GROUP == "Others"
    include("others.jl")
end

if GROUP == "All" || GROUP == "AD"
    include("ad/utils.jl")
    include("ad/distributions.jl")
end


text:
module RTRRT

include("rrt_functions.jl")

end


text:
if lowercase(get(ENV, "CI", "false")) == "true"
    include("install_dependencies.jl")
end

using TensND
using Test
using TimerOutputs
using LinearAlgebra, SymPy, Tensors, OMEinsum, Rotations

macro testsection(str, block)
    return quote
        @timeit "$($(esc(str)))" begin
            @testset "$($(esc(str)))" begin
                $(esc(block))
            end
        end
    end
end

reset_timer!()

opequal(x,y) = x == y || x ≈ y

include("test_bases.jl")
include("test_tens.jl")
include("test_tens_isotropic.jl")
include("test_special_tens.jl")
include("test_coorsystems.jl")

print_timer()
println()


text:
const RECOMPILE_BY_DEFAULT = true

Base.summary(prob::AbstractDiffEqFunction) = string(TYPE_COLOR, nameof(typeof(prob)),
                                                   NO_COLOR, ". In-place: ",
                                                   TYPE_COLOR, isinplace(prob),
                                                   NO_COLOR)
TreeViews.hastreeview(x::DiffEqBase.AbstractDiffEqFunction) = true
function TreeViews.treelabel(io::IO,x::DiffEqBase.AbstractDiffEqFunction,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Base.Text(Base.summary(x)))
end

"""
$(TYPEDEF)
"""
abstract type AbstractODEFunction{iip} <: AbstractDiffEqFunction{iip} end

"""
$(TYPEDEF)
"""
struct ODEFunction{iip,F,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} <: AbstractODEFunction{iip}
  f::F
  mass_matrix::TMM
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  syms::S
  colorvec::TCV
end

"""
$(TYPEDEF)
"""
struct SplitFunction{iip,F1,F2,TMM,C,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} <: AbstractODEFunction{iip}
  f1::F1
  f2::F2
  mass_matrix::TMM
  cache::C
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  syms::S
  colorvec::TCV
end

"""
$(TYPEDEF)
"""
struct DynamicalODEFunction{iip,F1,F2,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} <: AbstractODEFunction{iip}
  f1::F1
  f2::F2
  mass_matrix::TMM
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  syms::S
  colorvec::TCV
end

"""
$(TYPEDEF)
"""
abstract type AbstractDDEFunction{iip} <: AbstractDiffEqFunction{iip} end

"""
$(TYPEDEF)
"""
struct DDEFunction{iip,F,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} <: AbstractDDEFunction{iip}
  f::F
  mass_matrix::TMM
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  syms::S
  colorvec::TCV
end

"""
$(TYPEDEF)
"""
abstract type AbstractDiscreteFunction{iip} <: AbstractDiffEqFunction{iip} end

"""
$(TYPEDEF)
"""
struct DiscreteFunction{iip,F,Ta,S} <: AbstractDiscreteFunction{iip}
  f::F
  analytic::Ta
  syms::S
end

"""
$(TYPEDEF)
"""
abstract type AbstractSDEFunction{iip} <: AbstractDiffEqFunction{iip} end

"""
$(TYPEDEF)
"""
struct SDEFunction{iip,F,G,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,GG,TCV} <: AbstractSDEFunction{iip}
  f::F
  g::G
  mass_matrix::TMM
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  ggprime::GG
  syms::S
  colorvec::TCV
end

"""
$(TYPEDEF)
"""
struct SplitSDEFunction{iip,F1,F2,G,TMM,C,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} <: AbstractSDEFunction{iip}
  f1::F1
  f2::F2
  g::G
  mass_matrix::TMM
  cache::C
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  syms::S
  colorvec::TCV
end

"""
$(TYPEDEF)
"""
abstract type AbstractRODEFunction{iip} <: AbstractDiffEqFunction{iip} end

"""
$(TYPEDEF)
"""
struct RODEFunction{iip,F,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} <: AbstractRODEFunction{iip}
  f::F
  mass_matrix::TMM
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  syms::S
  colorvec::TCV
end

"""
$(TYPEDEF)
"""
abstract type AbstractDAEFunction{iip} <: AbstractDiffEqFunction{iip} end

"""
$(TYPEDEF)
"""
struct DAEFunction{iip,F,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} <: AbstractDAEFunction{iip}
  f::F
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  syms::S
  colorvec::TCV
end

"""
$(TYPEDEF)
"""
abstract type AbstractSDDEFunction{iip} <: AbstractDiffEqFunction{iip} end

"""
$(TYPEDEF)
"""
struct SDDEFunction{iip,F,G,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,GG,TCV} <: AbstractSDDEFunction{iip}
  f::F
  g::G
  mass_matrix::TMM
  analytic::Ta
  tgrad::Tt
  jac::TJ
  jvp::JVP
  vjp::VJP
  jac_prototype::JP
  sparsity::SP
  Wfact::TW
  Wfact_t::TWt
  paramjac::TPJ
  ggprime::GG
  syms::S
  colorvec::TCV
end
######### Backwards Compatibility Overloads

(f::ODEFunction)(args...) = f.f(args...)
(f::ODEFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)
(f::ODEFunction)(::Type{Val{:tgrad}},args...) = f.tgrad(args...)
(f::ODEFunction)(::Type{Val{:jac}},args...) = f.jac(args...)
(f::ODEFunction)(::Type{Val{:Wfact}},args...) = f.Wfact(args...)
(f::ODEFunction)(::Type{Val{:Wfact_t}},args...) = f.Wfact_t(args...)
(f::ODEFunction)(::Type{Val{:paramjac}},args...) = f.paramjac(args...)

function (f::DynamicalODEFunction)(u,p,t)
  ArrayPartition(f.f1(u.x[1],u.x[2],p,t),f.f2(u.x[1],u.x[2],p,t))
end
function (f::DynamicalODEFunction)(du,u,p,t)
  f.f1(du.x[1],u.x[1],u.x[2],p,t)
  f.f2(du.x[2],u.x[1],u.x[2],p,t)
end

(f::SplitFunction)(u,p,t) = f.f1(u,p,t) + f.f2(u,p,t)
(f::SplitFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)
function (f::SplitFunction)(du,u,p,t)
  f.f1(f.cache,u,p,t)
  f.f2(du,u,p,t)
  du .+= f.cache
end

(f::DiscreteFunction)(args...) = f.f(args...)
(f::DiscreteFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)

(f::DAEFunction)(args...) = f.f(args...)
(f::DAEFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)
(f::DAEFunction)(::Type{Val{:tgrad}},args...) = f.tgrad(args...)
(f::DAEFunction)(::Type{Val{:jac}},args...) = f.jac(args...)
(f::DAEFunction)(::Type{Val{:Wfact}},args...) = f.Wfact(args...)
(f::DAEFunction)(::Type{Val{:Wfact_t}},args...) = f.Wfact_t(args...)
(f::DAEFunction)(::Type{Val{:paramjac}},args...) = f.paramjac(args...)

(f::DDEFunction)(args...) = f.f(args...)
(f::DDEFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)

(f::SDEFunction)(args...) = f.f(args...)
(f::SDEFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)
(f::SDEFunction)(::Type{Val{:tgrad}},args...) = f.tgrad(args...)
(f::SDEFunction)(::Type{Val{:jac}},args...) = f.jac(args...)
(f::SDEFunction)(::Type{Val{:Wfact}},args...) = f.Wfact(args...)
(f::SDEFunction)(::Type{Val{:Wfact_t}},args...) = f.Wfact_t(args...)
(f::SDEFunction)(::Type{Val{:paramjac}},args...) = f.paramjac(args...)

(f::SDDEFunction)(args...) = f.f(args...)
(f::SDDEFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)
(f::SDDEFunction)(::Type{Val{:tgrad}},args...) = f.tgrad(args...)
(f::SDDEFunction)(::Type{Val{:jac}},args...) = f.jac(args...)
(f::SDDEFunction)(::Type{Val{:Wfact}},args...) = f.Wfact(args...)
(f::SDDEFunction)(::Type{Val{:Wfact_t}},args...) = f.Wfact_t(args...)
(f::SDDEFunction)(::Type{Val{:paramjac}},args...) = f.paramjac(args...)

(f::SplitSDEFunction)(u,p,t) = f.f1(u,p,t) + f.f2(u,p,t)
(f::SplitSDEFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)
function (f::SplitSDEFunction)(du,u,p,t)
  f.f1(f.cache,u,p,t)
  f.f2(du,u,p,t)
  du .+= f.cache
end

(f::RODEFunction)(args...) = f.f(args...)

######### Basic Constructor

function ODEFunction{iip,true}(f;
                 mass_matrix=I,
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if mass_matrix == I && typeof(f) <: Tuple
                  mass_matrix = ((I for i in 1:length(f))...,)
                 end

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 ODEFunction{iip,typeof(f),typeof(mass_matrix),typeof(analytic),typeof(tgrad),
                 typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),
                 typeof(Wfact),typeof(Wfact_t),
                 typeof(paramjac),typeof(syms),typeof(_colorvec)}(
                 f,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,syms,_colorvec)
end
function ODEFunction{iip,false}(f;
                 mass_matrix=I,
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 ODEFunction{iip,Any,Any,Any,Any,
                 Any,Any,Any,Any,Any,
                 Any,typeof(syms),typeof(_colorvec)}(
                 f,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,syms,_colorvec)
end
ODEFunction{iip}(f; kwargs...) where iip = ODEFunction{iip,RECOMPILE_BY_DEFAULT}(f; kwargs...)
ODEFunction{iip}(f::ODEFunction; kwargs...) where iip = f
ODEFunction(f; kwargs...) = ODEFunction{isinplace(f, 4),RECOMPILE_BY_DEFAULT}(f; kwargs...)
ODEFunction(f::ODEFunction; kwargs...) = f

@add_kwonly function SplitFunction(f1,f2,mass_matrix,cache,analytic,tgrad,jac,jvp,vjp,
                                   jac_prototype,sparsity,Wfact,Wfact_t,paramjac,
                                   syms,colorvec)
  f1 = typeof(f1) <: AbstractDiffEqOperator ? f1 : ODEFunction(f1)
  f2 = ODEFunction(f2)
  SplitFunction{isinplace(f2),typeof(f1),typeof(f2),typeof(mass_matrix),
              typeof(cache),typeof(analytic),typeof(tgrad),typeof(jac),typeof(jvp),typeof(vjp),
              typeof(jac_prototype),typeof(sparsity),
              typeof(Wfact),typeof(Wfact_t),typeof(paramjac),typeof(syms),
              typeof(colorvec)}(f1,f2,mass_matrix,cache,analytic,tgrad,jac,jvp,vjp,
              jac_prototype,sparsity,Wfact,Wfact_t,paramjac,syms,colorvec)
end
function SplitFunction{iip,true}(f1,f2;
                                 mass_matrix=I,_func_cache=nothing,
                                 analytic=nothing,
                                 tgrad = nothing,
                                 jac = nothing,
                                 jvp=nothing,
                                 vjp=nothing,
                                 jac_prototype = nothing,
                                 sparsity=jac_prototype,
                                 Wfact = nothing,
                                 Wfact_t = nothing,
                                 paramjac = nothing,
                                 syms = nothing,
                                 colorvec = nothing) where iip
  SplitFunction{iip,typeof(f1),typeof(f2),typeof(mass_matrix),
                typeof(_func_cache),typeof(analytic),
                typeof(tgrad),typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),
                typeof(Wfact),typeof(Wfact_t),typeof(paramjac),typeof(syms),
                typeof(colorvec)}(
                f1,f2,mass_matrix,_func_cache,analytic,tgrad,jac,jvp,vjp,jac_prototype,
                sparsity,Wfact,Wfact_t,paramjac,syms,colorvec)
end
function SplitFunction{iip,false}(f1,f2; mass_matrix=I,
                                  _func_cache=nothing,analytic=nothing,
                                  tgrad = nothing,
                                  jac = nothing,
                                  jvp=nothing,
                                  vjp=nothing,
                                  jac_prototype = nothing,
                                  sparsity=jac_prototype,
                                  Wfact = nothing,
                                  Wfact_t = nothing,
                                  paramjac = nothing,
                                  syms = nothing,
                                  colorvec = nothing) where iip
  SplitFunction{iip,Any,Any,Any,Any,Any,Any,Any,Any,
                Any,Any,Any,Any,Any}(
                f1,f2,mass_matrix,_func_cache,analytic,tgrad,jac,jvp,vjp,jac_prototype,
                sparsity,Wfact,Wfact_t,paramjac,syms,colorvec)
end
SplitFunction(f1,f2; kwargs...) = SplitFunction{isinplace(f2, 4)}(f1, f2; kwargs...)
SplitFunction{iip}(f1,f2; kwargs...) where iip =
SplitFunction{iip,RECOMPILE_BY_DEFAULT}(ODEFunction(f1),ODEFunction{iip}(f2); kwargs...)
SplitFunction(f::SplitFunction; kwargs...) = f

@add_kwonly function DynamicalODEFunction{iip}(f1,f2,mass_matrix,analytic,tgrad,jac,jvp,vjp,
                                   jac_prototype,sparsity,Wfact,Wfact_t,paramjac,
                                   syms,colorvec) where iip
  f1 = typeof(f1) <: AbstractDiffEqOperator ? f1 : ODEFunction(f1)
  f2 = ODEFunction(f2)
  DynamicalODEFunction{isinplace(f2),typeof(f1),typeof(f2),typeof(mass_matrix),
              typeof(analytic),typeof(tgrad),typeof(jac),typeof(jvp),typeof(vjp),
              typeof(jac_prototype),
              typeof(Wfact),typeof(Wfact_t),typeof(paramjac),typeof(syms),
              typeof(colorvec)}(f1,f2,mass_matrix,analytic,tgrad,jac,jvp,vjp,
              jac_prototype,sparsity,Wfact,Wfact_t,paramjac,syms,colorvec)
end

function DynamicalODEFunction{iip,true}(f1,f2;mass_matrix=I,
                                        analytic=nothing,
                                        tgrad=nothing,
                                        jac=nothing,
                                        jvp=nothing,
                                        vjp=nothing,
                                        jac_prototype=nothing,
                                        sparsity=jac_prototype,
                                        Wfact=nothing,
                                        Wfact_t=nothing,
                                        paramjac = nothing,
                                        syms = nothing,
                                        colorvec = nothing) where iip
  DynamicalODEFunction{iip,typeof(f1),typeof(f2),typeof(mass_matrix),
                typeof(analytic),
                typeof(tgrad),typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),
                typeof(Wfact),typeof(Wfact_t),typeof(paramjac),typeof(syms),
                typeof(colorvec)}(
                f1,f2,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,
                Wfact,Wfact_t,paramjac,syms,colorvec)
end

function DynamicalODEFunction{iip,false}(f1,f2;mass_matrix=I,
                                         analytic=nothing,
                                         tgrad=nothing,
                                         jac=nothing,
                                         jvp=nothing,
                                         vjp=nothing,
                                         jac_prototype=nothing,
                                         sparsity=jac_prototype,
                                         Wfact=nothing,
                                         Wfact_t=nothing,
                                         paramjac = nothing,
                                         syms = nothing,
                                         colorvec = nothing) where iip
       DynamicalODEFunction{iip,Any,Any,Any,Any,Any,Any,Any,
                            Any,Any,Any,Any,Any}(
                            f1,f2,mass_matrix,analytic,tgrad,
                            jac,jvp,vjp,jac_prototype,sparsity,
                            Wfact,Wfact_t,paramjac,syms,colorvec)
end

DynamicalODEFunction(f1,f2=nothing; kwargs...) = DynamicalODEFunction{isinplace(f1, 5)}(f1, f2; kwargs...)
DynamicalODEFunction{iip}(f1,f2; kwargs...) where iip =
DynamicalODEFunction{iip,RECOMPILE_BY_DEFAULT}(ODEFunction{iip}(f1), ODEFunction{iip}(f2); kwargs...)
DynamicalODEFunction(f::DynamicalODEFunction; kwargs...) = f

function DiscreteFunction{iip,true}(f;
                 analytic=nothing, syms=nothing) where iip
                 DiscreteFunction{iip,typeof(f),typeof(analytic),typeof(syms)}(
                 f,analytic,syms)
end
function DiscreteFunction{iip,false}(f;
                 analytic=nothing, syms=nothing) where iip
                 DiscreteFunction{iip,Any,Any,Any}(
                 f,analytic,syms)
end
DiscreteFunction{iip}(f; kwargs...) where iip = DiscreteFunction{iip,RECOMPILE_BY_DEFAULT}(f; kwargs...)
DiscreteFunction{iip}(f::DiscreteFunction; kwargs...) where iip = f
DiscreteFunction(f; kwargs...) = DiscreteFunction{isinplace(f, 4),RECOMPILE_BY_DEFAULT}(f; kwargs...)
DiscreteFunction(f::DiscreteFunction; kwargs...) = f

function SDEFunction{iip,true}(f,g;
                 mass_matrix=I,
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 ggprime = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 SDEFunction{iip,typeof(f),typeof(g),
                 typeof(mass_matrix),typeof(analytic),typeof(tgrad),
                 typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),typeof(Wfact),typeof(Wfact_t),
                 typeof(paramjac),typeof(syms),
                 typeof(ggprime),typeof(_colorvec)}(
                 f,g,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,ggprime,syms,_colorvec)
end
function SDEFunction{iip,false}(f,g;
                 mass_matrix=I,
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 ggprime = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 SDEFunction{iip,Any,Any,Any,Any,Any,
                 Any,Any,Any,Any,Any,
                 Any,typeof(syms),Any,typeof(_colorvec)}(
                 f,g,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,ggprime,syms,_colorvec)
end
SDEFunction{iip}(f,g; kwargs...) where iip = SDEFunction{iip,RECOMPILE_BY_DEFAULT}(f,g; kwargs...)
SDEFunction{iip}(f::SDEFunction,g; kwargs...) where iip = f
SDEFunction(f,g; kwargs...) = SDEFunction{isinplace(f, 4),RECOMPILE_BY_DEFAULT}(f,g; kwargs...)
SDEFunction(f::SDEFunction; kwargs...) = f

@add_kwonly function SplitSDEFunction(f1,f2,g,mass_matrix,cache,analytic,tgrad,jac,jvp,vjp,
                                   jac_prototype,Wfact,Wfact_t,paramjac,
                                   syms,colorvec)
  f1 = typeof(f1) <: AbstractDiffEqOperator ? f1 : SDEFunction(f1)
  f2 = SDEFunction(f2)
  SplitFunction{isinplace(f2),typeof(f1),typeof(f2),typeof(g),typeof(mass_matrix),
              typeof(cache),typeof(analytic),typeof(tgrad),typeof(jac),typeof(jvp),typeof(vjp),
              typeof(Wfact),typeof(Wfact_t),typeof(paramjac),typeof(syms),
              typeof(colorvec)}(f1,f2,mass_matrix,cache,analytic,tgrad,jac,
              jac_prototype,Wfact,Wfact_t,paramjac,syms,colorvec)
end

function SplitSDEFunction{iip,true}(f1,f2,g; mass_matrix=I,
                           _func_cache=nothing,analytic=nothing,
                           tgrad = nothing,
                           jac = nothing,
                           jac_prototype=nothing,
                           sparsity=jac_prototype,
                           jvp=nothing,
                           vjp=nothing,
                           Wfact = nothing,
                           Wfact_t = nothing,
                           paramjac = nothing,
                           syms = nothing,
                           colorvec = nothing) where iip
  SplitSDEFunction{iip,typeof(f1),typeof(f2),typeof(g),
              typeof(mass_matrix),typeof(_func_cache),
              typeof(analytic),
              typeof(tgrad),typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),
              typeof(Wfact),typeof(Wfact_t),typeof(paramjac),typeof(syms),
              typeof(colorvec)}(f1,f2,g,mass_matrix,_func_cache,analytic,
              tgrad,jac,jvp,vjp,jac_prototype,sparsity,
              Wfact,Wfact_t,paramjac,syms,colorvec)
end
function SplitSDEFunction{iip,false}(f1,f2,g; mass_matrix=I,
                            _func_cache=nothing,analytic=nothing,
                            tgrad = nothing,
                            jac = nothing,
                            jvp=nothing,
                            vjp=nothing,
                            jac_prototype=nothing,
                            sparsity=jac_prototype,
                            Wfact = nothing,
                            Wfact_t = nothing,
                            paramjac = nothing,
                            syms = nothing,
                            colorvec = nothing) where iip
  SplitSDEFunction{iip,Any,Any,Any,Any,Any,
                   Any,Any,Any,Any,
                   Any,Any,Any,Any,Any}(
                   f1,f2,g,mass_matrix,_func_cache,analytic,
                   tgrad,jac,jvp,vjp,jac_prototype,sparsity,
                   Wfact,Wfact_t,paramjac,syms,colorvec)
end
SplitSDEFunction(f1,f2,g; kwargs...) = SplitSDEFunction{isinplace(f2, 4)}(f1, f2, g; kwargs...)
SplitSDEFunction{iip}(f1,f2, g; kwargs...) where iip =
SplitSDEFunction{iip,RECOMPILE_BY_DEFAULT}(SDEFunction(f1,g), SDEFunction{iip}(f2,g), g; kwargs...)
SplitSDEFunction(f::SplitSDEFunction; kwargs...) = f

function RODEFunction{iip,true}(f;
                 mass_matrix=I,
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 RODEFunction{iip,typeof(f),typeof(mass_matrix),
                 typeof(analytic),typeof(tgrad),
                 typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),typeof(Wfact),typeof(Wfact_t),
                 typeof(paramjac),typeof(syms),typeof(_colorvec)}(
                 f,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,syms,_colorvec)
end
function RODEFunction{iip,false}(f;
                 mass_matrix=I,
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 RODEFunction{iip,Any,Any,Any,Any,
                 Any,Any,Any,Any,Any,
                 Any,typeof(syms),typeof(_colorvec)}(
                 f,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,syms,_colorvec)
end
RODEFunction{iip}(f; kwargs...) where iip = RODEFunction{iip,RECOMPILE_BY_DEFAULT}(f; kwargs...)
RODEFunction{iip}(f::RODEFunction; kwargs...) where iip = f
RODEFunction(f; kwargs...) = RODEFunction{isinplace(f, 5),RECOMPILE_BY_DEFAULT}(f; kwargs...)
RODEFunction(f::RODEFunction; kwargs...) = f

function DAEFunction{iip,true}(f;
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 DAEFunction{iip,typeof(f),typeof(analytic),typeof(tgrad),
                 typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),typeof(Wfact),typeof(Wfact_t),
                 typeof(paramjac),typeof(syms),typeof(_colorvec)}(
                 f,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,syms,_colorvec)
end
function DAEFunction{iip,false}(f;
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 DAEFunction{iip,Any,Any,Any,
                 Any,Any,Any,Any,Any,
                 Any,typeof(syms),typeof(_colorvec)}(
                 f,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,syms,_colorvec)
end
DAEFunction{iip}(f; kwargs...) where iip = DAEFunction{iip,RECOMPILE_BY_DEFAULT}(f; kwargs...)
DAEFunction{iip}(f::DAEFunction; kwargs...) where iip = f
DAEFunction(f; kwargs...) = DAEFunction{isinplace(f, 5),RECOMPILE_BY_DEFAULT}(f; kwargs...)
DAEFunction(f::DAEFunction; kwargs...) = f

function DDEFunction{iip,true}(f;
                 mass_matrix=I,
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 DDEFunction{iip,typeof(f),typeof(mass_matrix),typeof(analytic),typeof(tgrad),
                 typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),typeof(Wfact),typeof(Wfact_t),
                 typeof(paramjac),typeof(syms),typeof(_colorvec)}(
                 f,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,syms,_colorvec)
end
function DDEFunction{iip,false}(f;
                 mass_matrix=I,
                 analytic=nothing,
                 tgrad=nothing,
                 jac=nothing,
                 jvp=nothing,
                 vjp=nothing,
                 jac_prototype=nothing,
                 sparsity=jac_prototype,
                 Wfact=nothing,
                 Wfact_t=nothing,
                 paramjac = nothing,
                 syms = nothing,
                 colorvec = nothing) where iip

                 if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
                  if iip
                    jac = update_coefficients! #(J,u,p,t)
                  else
                    jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
                  end
                 end

                 if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
                   _colorvec = ArrayInterface.matrix_colors(jac_prototype)
                 else
                   _colorvec = colorvec
                 end

                 DDEFunction{iip,Any,Any,Any,Any,
                 Any,Any,Any,Any,Any,
                 Any,typeof(syms),typeof(_colorvec)}(
                 f,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
                 paramjac,syms,_colorvec)
end
DDEFunction{iip}(f; kwargs...) where iip = DDEFunction{iip,RECOMPILE_BY_DEFAULT}(f; kwargs...)
DDEFunction{iip}(f::DDEFunction; kwargs...) where iip = f
DDEFunction(f; kwargs...) = DDEFunction{isinplace(f, 5),RECOMPILE_BY_DEFAULT}(f; kwargs...)
DDEFunction(f::DDEFunction; kwargs...) = f

function SDDEFunction{iip,true}(f,g;
                                mass_matrix=I,
                                analytic=nothing,
                                tgrad=nothing,
                                jac=nothing,
                                jvp=nothing,
                                vjp=nothing,
                                jac_prototype=nothing,
                                sparsity=jac_prototype,
                                Wfact=nothing,
                                Wfact_t=nothing,
                                paramjac = nothing,
                                ggprime = nothing,
                                syms = nothing,
                                colorvec = nothing)  where iip
  if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
   if iip
     jac = update_coefficients! #(J,u,p,t)
   else
     jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
   end
  end

  if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
    _colorvec = ArrayInterface.matrix_colors(jac_prototype)
  else
    _colorvec = colorvec
  end

  SDDEFunction{iip,typeof(f),typeof(g),
  typeof(mass_matrix),typeof(analytic),typeof(tgrad),
  typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),typeof(Wfact),typeof(Wfact_t),
  typeof(paramjac),typeof(syms),
  typeof(ggprime),typeof(_colorvec)}(
  f,g,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
  paramjac,ggprime,syms,_colorvec)
end

function SDDEFunction{iip,false}(f,g;
                                 mass_matrix=I,
                                 analytic=nothing,
                                 tgrad=nothing,
                                 jac=nothing,
                                 jvp=nothing,
                                 vjp=nothing,
                                 jac_prototype=nothing,
                                 sparsity=jac_prototype,
                                 Wfact=nothing,
                                 Wfact_t=nothing,
                                 paramjac = nothing,
                                 ggprime = nothing,
                                 syms = nothing,
                                 colorvec = nothing) where iip

  if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
   if iip
     jac = update_coefficients! #(J,u,p,t)
   else
     jac = (u,p,t) -> update_coefficients!(deepcopy(jac_prototype),u,p,t)
   end
  end

  if jac_prototype !== nothing && colorvec === nothing && ArrayInterface.fast_matrix_colors(jac_prototype)
    _colorvec = ArrayInterface.matrix_colors(jac_prototype)
  else
    _colorvec = colorvec
  end

  SDDEFunction{iip,Any,Any,Any,Any,Any,
  Any,Any,Any,Any,Any,
  Any,typeof(syms),Any,typeof(_colorvec)}(
  f,g,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,Wfact,Wfact_t,
  paramjac,ggprime,syms,_colorvec)
end
SDDEFunction{iip}(f,g; kwargs...) where iip = SDDEFunction{iip,RECOMPILE_BY_DEFAULT}(f,g; kwargs...)
SDDEFunction{iip}(f::SDDEFunction,g; kwargs...) where iip = f
SDDEFunction(f,g; kwargs...) = SDDEFunction{isinplace(f, 5),RECOMPILE_BY_DEFAULT}(f,g; kwargs...)
SDDEFunction(f::SDDEFunction; kwargs...) = f

########## Existance Functions

# Check that field/property exists (may be nothing)
__has_jac(f) = isdefined(f, :jac)
__has_jvp(f) = isdefined(f, :jvp)
__has_vjp(f) = isdefined(f, :vjp)
__has_tgrad(f) = isdefined(f, :tgrad)
__has_Wfact(f) = isdefined(f, :Wfact)
__has_Wfact_t(f) = isdefined(f, :Wfact_t)
__has_paramjac(f) = isdefined(f, :paramjac)
__has_syms(f) = isdefined(f, :syms)
__has_analytic(f) = isdefined(f, :analytic)
__has_colorvec(f) = isdefined(f, :colorvec)

# compatibility
has_invW(f::AbstractDiffEqFunction) = false
has_analytic(f::AbstractDiffEqFunction) = __has_analytic(f) && f.analytic !== nothing
has_jac(f::AbstractDiffEqFunction) = __has_jac(f) && f.jac !== nothing
has_jvp(f::AbstractDiffEqFunction) = __has_jvp(f) && f.jvp !== nothing
has_vjp(f::AbstractDiffEqFunction) = __has_vjp(f) && f.vjp !== nothing
has_tgrad(f::AbstractDiffEqFunction) = __has_tgrad(f) && f.tgrad !== nothing
has_Wfact(f::AbstractDiffEqFunction) = __has_Wfact(f) && f.Wfact !== nothing
has_Wfact_t(f::AbstractDiffEqFunction) = __has_Wfact_t(f) && f.Wfact_t !== nothing
has_paramjac(f::AbstractDiffEqFunction) = __has_paramjac(f) && f.paramjac !== nothing
has_syms(f::AbstractDiffEqFunction) = __has_syms(f) && f.syms !== nothing
has_colorvec(f::AbstractDiffEqFunction) = __has_colorvec(f) && f.colorvec !== nothing

# TODO: find an appropriate way to check `has_*`
has_jac(f::Union{SplitFunction,SplitSDEFunction}) = has_jac(f.f1)
has_jvp(f::Union{SplitFunction,SplitSDEFunction}) = has_jvp(f.f1)
has_vjp(f::Union{SplitFunction,SplitSDEFunction}) = has_vjp(f.f1)
has_tgrad(f::Union{SplitFunction,SplitSDEFunction}) = has_tgrad(f.f1)
has_Wfact(f::Union{SplitFunction,SplitSDEFunction}) = has_Wfact(f.f1)
has_Wfact_t(f::Union{SplitFunction,SplitSDEFunction}) = has_Wfact_t(f.f1)
has_paramjac(f::Union{SplitFunction,SplitSDEFunction}) = has_paramjac(f.f1)
has_colorvec(f::Union{SplitFunction,SplitSDEFunction}) = has_colorvec(f.f1)

has_jac(f::DynamicalODEFunction) = has_jac(f.f1)
has_jvp(f::DynamicalODEFunction) = has_jvp(f.f1)
has_vjp(f::DynamicalODEFunction) = has_vjp(f.f1)
has_tgrad(f::DynamicalODEFunction) = has_tgrad(f.f1)
has_Wfact(f::DynamicalODEFunction) = has_Wfact(f.f1)
has_Wfact_t(f::DynamicalODEFunction) = has_Wfact_t(f.f1)
has_paramjac(f::DynamicalODEFunction) = has_paramjac(f.f1)
has_colorvec(f::DynamicalODEFunction) = has_colorvec(f.f1)

has_jac(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_jac(f.f)
has_jvp(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_jvp(f.f)
has_vjp(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_vjp(f.f)
has_tgrad(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_tgrad(f.f)
has_Wfact(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_Wfact(f.f)
has_Wfact_t(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_Wfact_t(f.f)
has_paramjac(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_paramjac(f.f)
has_colorvec(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_colorvec(f.f)

######### Additional traits

islinear(f) = false # fallback
islinear(::AbstractDiffEqFunction) = false
islinear(f::ODEFunction) = islinear(f.f)
islinear(f::SplitFunction) = islinear(f.f1)

######### Compatibility Constructor from Tratis

function Base.convert(::Type{ODEFunction}, f)
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  ODEFunction(f;analytic=analytic,tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end
function Base.convert(::Type{ODEFunction{iip}},f) where iip
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  ODEFunction{iip,RECOMPILE_BY_DEFAULT}(f;analytic=analytic,tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end

function Base.convert(::Type{DiscreteFunction},f)
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  DiscreteFunction(f;analytic=analytic,syms=syms)
end
function Base.convert(::Type{DiscreteFunction{iip}},f) where iip
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  DiscreteFunction{iip,RECOMPILE_BY_DEFAULT}(f;analytic=analytic,syms=syms)
end

function Base.convert(::Type{DAEFunction},f)
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  DAEFunction(f;analytic=analytic,tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end
function Base.convert(::Type{DAEFunction{iip}},f) where iip
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  DAEFunction{iip,RECOMPILE_BY_DEFAULT}(f;analytic=analytic,tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end

function Base.convert(::Type{DDEFunction},f)
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  DDEFunction(f;analytic=analytic,syms=syms,colorvec=colorvec)
end
function Base.convert(::Type{DDEFunction{iip}},f) where iip
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  DDEFunction{iip,RECOMPILE_BY_DEFAULT}(f;analytic=analytic,syms=syms,colorvec=colorvec)
end

function Base.convert(::Type{SDEFunction},f,g)
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  SDEFunction(f,g;analytic=analytic,tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end
function Base.convert(::Type{SDEFunction{iip}},f,g) where iip
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  SDEFunction{iip,RECOMPILE_BY_DEFAULT}(f,g;analytic=analytic,
              tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end

function Base.convert(::Type{RODEFunction},f)
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  RODEFunction(f;analytic=analytic,tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end
function Base.convert(::Type{RODEFunction{iip}},f) where iip
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  RODEFunction{iip,RECOMPILE_BY_DEFAULT}(f;analytic=analytic,
              tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end

function Base.convert(::Type{SDDEFunction},f,g)
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  SDDEFunction(f,g;analytic=analytic,tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end

function Base.convert(::Type{SDDEFunction{iip}},f,g) where iip
  if __has_analytic(f)
    analytic = f.analytic
  else
    analytic = nothing
  end
  if __has_jac(f)
    jac = f.jac
  else
    jac = nothing
  end
  if __has_jvp(f)
    jvp = f.jvp
  else
    jvp = nothing
  end
  if __has_vjp(f)
    vjp = f.vjp
  else
    vjp = nothing
  end
  if __has_tgrad(f)
    tgrad = f.tgrad
  else
    tgrad = nothing
  end
  if __has_Wfact(f)
    Wfact = f.Wfact
  else
    Wfact = nothing
  end
  if __has_Wfact_t(f)
    Wfact_t = f.Wfact_t
  else
    Wfact_t = nothing
  end
  if __has_paramjac(f)
    paramjac = f.paramjac
  else
    paramjac = nothing
  end
  if __has_syms(f)
    syms = f.syms
  else
    syms = nothing
  end
  if __has_colorvec(f)
    colorvec = f.colorvec
  else
    colorvec = nothing
  end
  SDDEFunction{iip,RECOMPILE_BY_DEFAULT}(f,g;analytic=analytic,
              tgrad=tgrad,jac=jac,jvp=jvp,vjp=vjp,Wfact=Wfact,
              Wfact_t=Wfact_t,paramjac=paramjac,syms=syms,colorvec=colorvec)
end


text:
function genFcode(v,av,n::Array{ADnode,1},fname=[])
# generate source code for the Autodiff forward pass
# (c) David Barber, University College London 2015
   if isempty(fname)
       iostream=STDOUT
       str="function ("
   else
       iostream=open(fname*"!.jl","w")
       str="function "*fname*"!("
   end
    N=length(n)
    str=str*"v,av)"
    println(iostream,str)

    # forward pass:
    for i=1:N
        if !(n[i].input)
            str="v[$i],av[$i]="*string(n[i].f)*"("
            for j=n[i].parents
                tstr="::"*string(typeof(v[j]))
                str=str*"v[$j]"*tstr*","
            end
            str=str[1:end-1]*")"
            println(iostream,str)
        end
    end
    println(iostream,"end")
    if !isempty(fname)
        close(iostream)
    end
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using Bayesutil
using Test

@testset "Bayesutil.jl" begin
    # Write your own tests here.
    @test Bayesutil.f(2,1) == 7
    @test Bayesutil.f(0,0) == 0
end


text:
function exact_expectation_inside_log_of_lower_bound(; y=y, idx=idx, S=S, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    expectationterm1(;y=y, μα=μα, Σα=Σα, μβ=μβ, Σβ=Σβ, μf=μf, Σf=Σf, S=S, idx=idx) +
    expectationterm2(;μf=μf, Σf=Σf, K=K)     +
    expectationterm3(;μα=μα, Σα=Σα, v²ₐ=v²ₐ) +
    expectationterm4(;μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ)

end

##########
# TERM 1 #
##########

function expectationterm1(;y=y, μα=μα, Σα=Σα, μβ=μβ, Σβ=Σβ, μf=μf, Σf=Σf, S=S, idx=idx)

    N = length(y)

    ## look out for the diagonal and the diff between Σα and Σa
    Ā = Diagonal(μα[idx]) ; Σa = Diagonal(diag(Σα)[idx])

    ## look out for the diagonal and the diff between Σβ and Σb
    b̄ = μβ[idx] ; Σb = Diagonal(diag(Σβ)[idx])

    subterm0 = -0.5*logdet(S) - 0.5*N*log(2π)

    subterm1 = -0.5*dot(y, S\y)

    subterm2 = -0.5*dot(y, -S\Ā*μf) * 2

    subterm3 = -0.5*dot(y, -S\b̄) * 2

    # subterm4 = -0.5*sum((μf.^2 .+ diag(Σf)) .* (Ā.diag.^2 .+ Σa.diag) ./ S.diag)
    subterm4 = -0.5*(μf'*(S\(Ā.^2 + Σa))*μf)  -0.5*tr(Σf*(S\(Ā.^2 + Σa)))

    subterm5 = -0.5*dot(b̄, S\Ā*μf) * 2

    subterm6 = -0.5*dot(b̄, S\b̄) -0.5*tr(S\Σb)

    (subterm0 + subterm1 + subterm2 + subterm3 + subterm4 + subterm5 + subterm6)

end

##########
# TERM 2 #
##########

function expectationterm2(;μf=μf, Σf=Σf, K=K)
    N = length(μf)
    -0.5*N*log(2π) -0.5*logdet(K) - 0.5*μf'*(K\μf) - 0.5*tr(K\Σf)
end

##########
# TERM 3 #
##########

function expectationterm3(;μα=μα, Σα=Σα, v²ₐ=v²ₐ)
    F = length(μα)
    - 0.5*F*log(2π) - 0.5*F*log(v²ₐ) - 0.5*dot(μα, μα)/v²ₐ - 0.5*tr(Σα)/v²ₐ
end

##########
# TERM 4 #
##########

function expectationterm4(;μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ)
    F = length(μβ)
    - 0.5*F*log(2π) - 0.5*F*log(v²ᵦ) - 0.5*dot(μβ,μβ)/v²ᵦ - 0.5*tr(Σβ)/v²ᵦ
end