text:
"""
`ParametersDGVIP1`: Parameters for right-hand side function of Discontinuous Galerkin Variational Integrator.

### Parameters

* `DT`: data type
* `TT`: parameter type
* `D`: dimension of the system
* `S`: number of basis nodes
* `R`: number of quadrature nodes

### Fields

* `Θ`:  function of the noncanonical one-form (∂L/∂v)
* `f`:  function of the force (∂L/∂q)
* `g`:  function of the projection ∇ϑ(q)⋅v
* `Δt`: time step
* `b`:  quadrature weights
* `c`:  quadrature nodes
* `m`:  mass matrix
* `a`:  derivative matrix
* `r⁻`: reconstruction coefficients, jump lhs value
* `r⁺`: reconstruction coefficients, jump rhs value
* `t`:  current time
* `q`:  current solution of qₙ
* `q⁻`: current solution of qₙ⁻
* `q⁺`: current solution of qₙ⁺
"""
mutable struct ParametersDGVIP1{DT,TT,D,S,R,ΘT,FT,GT} <: Parameters{DT,TT}
    Θ::ΘT
    f::FT
    g::GT

    Δt::TT

    b::Vector{TT}
    c::Vector{TT}
    m::Matrix{TT}
    a::Matrix{TT}
    r⁻::Vector{TT}
    r⁺::Vector{TT}

    t::TT

    q::Vector{DT}
    q⁻::Vector{DT}
    q⁺::Vector{DT}
end

function ParametersDGVIP1(Θ::ΘT, f::FT, g::GT, Δt::TT,
                b::Vector{TT}, c::Vector{TT}, m::Matrix{TT}, a::Matrix{TT}, r⁻::Vector{TT}, r⁺::Vector{TT},
                q::Vector{DT}, q⁻::Vector{DT}, q⁺::Vector{DT}) where {DT,TT,ΘT,FT,GT}

    @assert length(q)  == length(q⁻)  == length(q⁺)
    @assert length(b)  == length(c)
    @assert length(r⁻) == length(r⁺)

    D = length(q)
    S = length(r⁻)
    R = length(c)

    println()
    println("  Discontinuous Galerkin Variational Integrator")
    println("  =============================================")
    println()
    println("    b = ", b)
    println("    c = ", c)
    println("    m = ", m)
    println("    a = ", a)
    println("    r⁻= ", r⁻)
    println("    r⁺= ", r⁺)
    println()

    ParametersDGVIP1{DT,TT,D,S,R,ΘT,FT,GT}(
                Θ, f, g, Δt, b, c, m, a, r⁻, r⁺, 0, q, q⁻, q⁺)
end

function ParametersDGVIP1(Θ::ΘT, f::FT, g::GT, Δt::TT,
                basis::Basis{TT}, quadrature::Quadrature{TT},
                q::Vector{DT}, q⁻::Vector{DT}, q⁺::Vector{DT}) where {DT,TT,ΘT,FT,GT}

    # compute coefficients
    m = zeros(TT, nnodes(quadrature), nbasis(basis))
    a = zeros(TT, nnodes(quadrature), nbasis(basis))
    r⁻= zeros(TT, nbasis(basis))
    r⁺= zeros(TT, nbasis(basis))

    for i in 1:nbasis(basis)
        for j in 1:nnodes(quadrature)
            m[j,i] = evaluate(basis, i, nodes(quadrature)[j])
            a[j,i] = derivative(basis, i, nodes(quadrature)[j])
        end
        r⁻[i] = evaluate(basis, i, one(TT))
        r⁺[i] = evaluate(basis, i, zero(TT))
    end

    ParametersDGVIP1(Θ, f, g, Δt, weights(quadrature), nodes(quadrature), m, a, r⁻, r⁺, q, q⁻, q⁺)
end

function ParametersDGVIP1(Θ, f, g, Δt, basis, quadrature, q, q⁻)
    q⁺ = zero(q)
    q⁺ .= q
    ParametersDGVIP1(Θ, f, g, Δt, basis, quadrature, q, q⁻, q⁺)
end

"Compute stages of variational partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    cache = IntegratorCacheDGVI{ST,D,S,R}()

    quote
        @assert length(x) == length(b)

        compute_stages!(x, $cache, params)

        compute_rhs!(b, $cache, params)
    end
end

function compute_stages!(x, cache::IntegratorCacheDGVI{ST,D,S}, params::ParametersDGVIP1{DT,TT,D,S}) where {ST,DT,TT,D,S}
    # copy x to X
    for i in 1:S
        for k in 1:D
            cache.X[i][k] = x[D*(i-1)+k]
        end
    end

    # copy x to q̅=qₙ+₁
    for k in 1:D
        cache.q̅[k] = x[D*S+k]
    end

    # copy x to q̅⁺=qₙ+₁⁺
    for k in 1:D
        cache.q̅⁺[k] = x[D*(S+1)+k]
    end

    # compute Q, qₙ⁺ and qₙ₊₁⁻
    compute_stages_q!(cache, params)

    # compute V
    compute_stages_v!(cache, params)

    # compute P and F
    compute_stages_p!(cache, params)

    # compute jump
    compute_stages_λ!(cache, params)
end

"Compute solution at quadrature nodes and across jump."
function compute_stages_q!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local q::ST
    local q⁺::ST
    local q̅⁻::ST

    local X = cache.X
    local Q = cache.Q

    # copy q and q⁻
    cache.q  .= params.q
    cache.q⁻ .= params.q⁻

    # compute Q
    for i in 1:R
        for k in 1:D
            q = 0
            for j in 1:S
                q += params.m[i,j] * X[j][k]
            end
            Q[i][k] = q
        end
    end

    # compute qₙ⁺ and qₙ₊₁⁻
    for k in 1:D
        q⁺ = 0
        q̅⁻ = 0
        for i in 1:S
            q⁺ += params.r⁺[i] * X[i][k]
            q̅⁻ += params.r⁻[i] * X[i][k]
        end
        cache.q⁺[k] = q⁺
        cache.q̅⁻[k] = q̅⁻
    end
end

"Compute velocities at quadrature nodes."
function compute_stages_v!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local v::ST

    for i in 1:R
        for k in 1:D
            v = 0
            for j in 1:S
                v += params.a[i,j] * cache.X[j][k]
            end
            cache.V[i][k] = v / params.Δt
        end
    end
end

"Compute one-form and forces at quadrature nodes."
function compute_stages_p!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local tᵢ::TT

    # compute P=ϑ(Q) and F=f(Q)
    for i in 1:R
        tᵢ = params.t + params.Δt * params.c[i]
        params.Θ(tᵢ, cache.Q[i], cache.V[i], cache.P[i])
        params.f(tᵢ, cache.Q[i], cache.V[i], cache.F[i])
    end
end

function compute_stages_λ!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local t₀::TT = params.t
    local t₁::TT = params.t + params.Δt

    # compute ϕ and ϕ̅
    cache.ϕ  .= 0.5 * (cache.q⁻ .+ cache.q⁺)
    cache.ϕ⁻ .= 0.5 * (cache.q⁻ .+ cache.q )
    cache.ϕ⁺ .= 0.5 * (cache.q  .+ cache.q⁺)

    cache.ϕ̅  .= 0.5 * (cache.q̅⁻ .+ cache.q̅⁺)
    cache.ϕ̅⁻ .= 0.5 * (cache.q̅⁻ .+ cache.q̅ )
    cache.ϕ̅⁺ .= 0.5 * (cache.q̅  .+ cache.q̅⁺)

    # compute λ and λ̅
    cache.λ  .= cache.q⁺ .- cache.q⁻
    cache.λ⁻ .= cache.q  .- cache.q⁻
    cache.λ⁺ .= cache.q⁺ .- cache.q

    cache.λ̅  .= cache.q̅⁺ .- cache.q̅⁻
    cache.λ̅⁻ .= cache.q̅  .- cache.q̅⁻
    cache.λ̅⁺ .= cache.q̅⁺ .- cache.q̅

    # compute ϑ
    params.Θ(t₀, cache.q,  cache.q,  cache.θ)
    params.Θ(t₀, cache.q⁻, cache.q⁻, cache.θ⁻)
    params.Θ(t₀, cache.q⁺, cache.q⁺, cache.θ⁺)

    params.Θ(t₁, cache.q̅,  cache.q̅,  cache.Θ̅)
    params.Θ(t₁, cache.q̅⁻, cache.q̅⁻, cache.Θ̅⁻)
    params.Θ(t₁, cache.q̅⁺, cache.q̅⁺, cache.Θ̅⁺)

    # compute projection
    params.g(t₀, cache.q,  cache.λ,  cache.g)
    params.g(t₀, cache.q⁻, cache.λ⁻, cache.g⁻)
    params.g(t₀, cache.q⁺, cache.λ⁺, cache.g⁺)

    params.g(t₁, cache.q̅,  cache.λ̅,  cache.g̅)
    params.g(t₁, cache.q̅⁻, cache.λ̅⁻, cache.g̅⁻)
    params.g(t₁, cache.q̅⁺, cache.λ̅⁺, cache.g̅⁺)

    # # compute ϑ
    # params.Θ(t₀, cache.ϕ⁻, cache.ϕ⁻, cache.θ⁻)
    # params.Θ(t₀, cache.ϕ⁺, cache.ϕ⁺, cache.θ⁺)
    # params.Θ(t₁, cache.ϕ̅⁻, cache.ϕ̅⁻, cache.Θ̅⁻)
    #
    # # compute projection
    # params.g(t₀, cache.ϕ⁻, cache.λ⁻, cache.g⁻)
    # params.g(t₀, cache.ϕ⁺, cache.λ⁺, cache.g⁺)
    # params.g(t₁, cache.ϕ̅⁻, cache.λ̅⁻, cache.g̅⁻)
end

function compute_rhs!(b::Vector{ST}, cache::IntegratorCacheDGVI{ST,D,S,R},
                params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local z::ST

    # compute b = - [(P-AF)]
    for i in 1:S
        for k in 1:D
            z = 0
            for j in 1:R
                z += params.b[j] * params.m[j,i] * cache.F[j][k] * params.Δt
                z += params.b[j] * params.a[j,i] * cache.P[j][k]
            end

            z += params.r⁺[i] * 0.5 * ( cache.θ[k] + cache.θ⁺[k] )
            z -= params.r⁻[i] * 0.5 * ( cache.Θ̅[k] + cache.Θ̅⁻[k] )

            z += params.r⁺[i] * 0.5 * cache.g⁺[k]
            z += params.r⁻[i] * 0.5 * cache.g̅⁻[k]

            b[D*(i-1)+k] = z
        end
    end

    # compute b = qₙ⁺ - r⁺⋅X
    for k in 1:D
        b[D*S+k] = params.q⁺[k] - cache.q⁺[k]
    end

    for k in 1:D
        b[D*(S+1)+k] = cache.Θ̅⁺[k] - cache.Θ̅⁻[k] - cache.g̅[k]
    end
end

"""
`IntegratorDGVIP1`: Discontinuous Galerkin Variational Integrator.

### Parameters

### Fields

* `equation`: Implicit Ordinary Differential Equation
* `basis`: piecewise polynomial basis
* `quadrature`: numerical quadrature rule
* `Δt`: time step
* `params`: ParametersDGVIP1
* `solver`: nonlinear solver
* `iguess`: initial guess
* `q`: current solution vector for trajectory
* `p`: current solution vector for one-form
* `cache`: temporary variables for nonlinear solver
"""
struct IntegratorDGVIP1{DT,TT,D,S,R,ΘT,FT,GT,VT,FPT,ST,IT,BT<:Basis} <: DeterministicIntegrator{DT,TT}
    equation::IODE{DT,TT,ΘT,FT,GT,VT}

    basis::BT
    quadrature::Quadrature{TT,R}

    Δt::TT

    params::FPT
    solver::ST
    iguess::InitialGuessPODE{DT,TT,VT,FT,IT}

    q::Vector{DT}
    q⁻::Vector{DT}
    q⁺::Vector{DT}

    cache::IntegratorCacheDGVI{DT}
end

function IntegratorDGVIP1(equation::IODE{DT,TT,ΘT,FT,GT,VT}, basis::Basis{TT,P},
                quadrature::Quadrature{TT,R}, Δt::TT;
                interpolation=HermiteInterpolation{DT}) where {DT,TT,ΘT,FT,GT,VT,P,R}

    D = equation.d
    S = nbasis(basis)

    N = D*(S+2)

    # create solution vector for nonlinear solver
    x = zeros(DT,N)

    # create solution vectors
    q  = zeros(DT,D)
    q⁻ = zeros(DT,D)
    q⁺ = zeros(DT,D)

    # create cache for internal stage vectors and update vectors
    cache = IntegratorCacheDGVI{DT,D,S,R}()

    # create params
    params = ParametersDGVIP1(equation.α, equation.f, equation.g,
                Δt, basis, quadrature, q, q⁻, q⁺)

    # create rhs function for nonlinear solver
    function_stages = (x,b) -> function_stages!(x, b, params)

    # create nonlinear solver
    solver = get_config(:nls_solver)(x, function_stages)

    # create initial guess
    iguess = InitialGuessPODE(interpolation, equation, Δt)

    # create integrator
    IntegratorDGVIP1{DT, TT, D, S, R, ΘT, FT, GT, VT, typeof(params), typeof(solver),
                typeof(iguess.int), typeof(basis)}(
                equation, basis, quadrature, Δt, params, solver, iguess,
                q, q⁻, q⁺, cache)
end

function initialize!(int::IntegratorDGVIP1, sol::Union{SolutionPODE, SolutionPDAE}, m::Int)
    @assert m ≥ 1
    @assert m ≤ sol.ni

    # copy initial conditions from solution
    get_initial_conditions!(sol, int.q, int.q⁺, m)
    int.q⁻ .= int.q

    # initialise initial guess
    initialize!(int.iguess, m, sol.t[0], int.q, int.q⁺)
end

function update_solution!(int::IntegratorDGVIP1{DT,TT}, cache::IntegratorCacheDGVI{DT}) where {DT,TT}
    int.q  .= cache.q̅
    int.q⁻ .= cache.q̅⁻
    int.q⁺ .= cache.q̅⁺
end

@generated function initial_guess!(int::IntegratorDGVIP1{DT,TT, D, S, R}, m::Int) where {DT,TT,D,S,R}
    v = zeros(DT,D)
    y = zeros(DT,D)
    z = zeros(DT,D)

    quote
        # compute initial guess
        if nnodes(int.basis) > 0
            for i in 1:S
                evaluate!(int.iguess, m, $y, $z, $v, nodes(int.basis)[i], nodes(int.basis)[i])
                for k in 1:D
                    int.solver.x[D*(i-1)+k] = $y[k]
                end
            end
        else
            for i in 1:S
                for k in 1:D
                    int.solver.x[D*(i-1)+k] = 0
                end
            end
        end

        evaluate!(int.iguess, m, $y, $z, $v, one(TT), one(TT))
        for k in 1:D
            int.solver.x[D*(S+0)+k] = $y[k]
            int.solver.x[D*(S+1)+k] = $y[k]
        end
    end
end

function integrate_step!(int::IntegratorDGVIP1{DT,TT}, sol::Union{SolutionPODE{DT,TT}, SolutionPDAE{DT,TT}}, m::Int, n::Int) where {DT,TT}
    # set time for nonlinear solver
    int.params.t = sol.t[0] + (n-1)*int.Δt

    # compute initial guess
    initial_guess!(int, m)

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute final update
    compute_stages!(int.solver.x, int.cache, int.params)

    # copy solution from cache to integrator
    update_solution!(int, int.cache)

    # copy solution to initial guess for next time step
    update!(int.iguess, m, sol.t[0] + n*int.Δt, int.q, int.q⁺)

    # take care of periodic solutions
    cut_periodic_solution!(int.q,  int.equation.periodicity)
    cut_periodic_solution!(int.q⁻, int.equation.periodicity)
    cut_periodic_solution!(int.q⁺, int.equation.periodicity)

    # copy to solution
    copy_solution!(sol, int.q, int.q⁺, n, m)
end


text:
module SymDiff

include("differentiate.jl")
export differentiate

include("simplify.jl")
export simplify

include("subs.jl")
export subs

end # module

ext: jl


text:
# Paths are an important concept for BIM (and other things)
using IntervalSets
export open_path,
       closed_path,
       open_path_ops,
       closed_path_ops,
       MoveOp,
       MoveToOp,
       LineOp,
       LineToOp,
       ArcOp,
       CloseOp,
       arc_path,
       circular_path,
       rectangular_path,
       centered_rectangular_path,
       open_polygonal_path,
       closed_polygonal_path,
       polygonal_path,
       open_spline_path,
       closed_spline_path,
       path_set,
       open_path_sequence,
       closed_path_sequence,
       path_sequence,
       translate,
       stroke,
       fill,
       path_length,
       location_at_length,
       path_start,
       path_end,
       subpath,
       join_paths,
       subpaths

path_tolerance = Parameter(1e-10)
coincident_path_location(p1::Loc, p2::Loc) = distance(p1, p2) < path_tolerance()

abstract type Path end

import Base.getindex, Base.firstindex, Base.lastindex
getindex(p::Path, i::Real) = location_at_length(p, i)
firstindex(p::Path) = 0
lastindex(p::Path) = path_length(p)
getindex(p::Path, i::ClosedInterval) = subpath(p, i.left, i.right)
path_start(p::Path) = location_at_length(p, 0)
path_end(p::Path) = location_at_length(p, path_length(p))

abstract type OpenPath <: Path end
abstract type ClosedPath <: Path end

struct ArcPath <: OpenPath
    center::Loc
    radius::Real
    start_angle::Real
    amplitude::Real
end
arc_path(center::Loc=u0(), radius::Real=1, start_angle::Real=0, amplitude::Real=pi) =
  false ? #amplitude < 0 ?
    ArcPath(center, radius, start_angle + amplitude, - amplitude) :
    ArcPath(center, radius, start_angle, amplitude)
struct CircularPath <: ClosedPath
    center::Loc
    radius::Real
end
circular_path(Center::Loc=u0(), Radius::Real=1; center::Loc=Center, radius::Real=Radius) = CircularPath(center, radius)
struct RectangularPath <: ClosedPath
    corner::Loc
    dx::Real
    dy::Real
end
rectangular_path(corner::Loc=u0(), dx::Real=1, dy::Real=1) = RectangularPath(corner, dx, dy)
centered_rectangular_path(p, dx, dy) =
  rectangular_path(p-vxy(dx/2, dy/2), dx, dy)

# Sections also use paths, but they are centered at the origin
export rectangular_profile,
       circular_profile,
       top_aligned_rectangular_profile

rectangular_profile(Width::Real=1, Height::Real=1; width::Real=Width, height::Real=Height) =
  centered_rectangular_path(u0(), width, height)

circular_profile(Radius::Real=1; radius::Real=Radius) =
  circular_path(u0(), radius)

top_aligned_rectangular_profile(Width::Real=1, Height::Real=1; width::Real=Width, height::Real=Height) =
  rectangular_path(xy(-width/2,-height), width, height)

struct OpenPolygonalPath <: OpenPath
    vertices::Locs
end
open_polygonal_path(vertices=[u0(), x(), xy(), y()]) = OpenPolygonalPath(vertices)

struct ClosedPolygonalPath <: ClosedPath
    vertices::Locs
end
closed_polygonal_path(vertices=[u0(), x(), xy(), y()]) = ClosedPolygonalPath(ensure_no_repeated_locations(vertices))

PolygonalPath = Union{OpenPolygonalPath, ClosedPolygonalPath}

polygonal_path(vertices) =
  coincident_path_location(vertices[1], vertices[end]) ?
    closed_polygonal_path(vertices[1:end-1]) :
    open_polygonal_path(vertices)

ensure_no_repeated_locations(locs) =
    begin
        @assert (locs[1] != locs[end])
        locs
    end

path_start(path::PolygonalPath) = path.vertices[1]
path_end(path::OpenPolygonalPath) = path.vertices[end]
path_end(path::ClosedPolygonalPath) = path.vertices[1]

join_paths(p1::OpenPolygonalPath, p2::OpenPolygonalPath) =
  coincident_path_location(path_end(p1), path_start(p2)) ?
    polygonal_path([p1.vertices..., p2.vertices...]) :
    error("Paths are non contiguous")

# Splines

struct OpenSplinePath <: OpenPath
    vertices::Locs
    v0::Union{Bool,Vec}
    v1::Union{Bool,Vec}
end
open_spline_path(vertices=[u0(), x(), xy(), y()], v0=false, v1=false) =
    OpenSplinePath(vertices, v0, v1)

struct ClosedSplinePath <: ClosedPath
    vertices::Locs
end
closed_spline_path(vertices=[u0(), x(), xy(), y()]) = ClosedSplinePath(ensure_no_repeated_locations(vertices))

# There is a set of operations over Paths:
# 1. translate a path a given vector
# 2. stroke a path
# 3. fill a (presumably closed) path
# 4. compute a path location given a length from the path beginning
# 5. compute a sub path from a path, a length from the path begining and a length increment
# 6. Produce the meta representation of a path

translate(path::CircularPath, v::Vec) = circular_path(path.center + v, path.radius)
translate(path::ArcPath, v::Vec) = arc_path(path.center + v, path.radius, path.start_angle, path.amplitude)
translate(path::RectangularPath, v::Vec) = rectangular_path(path.corner + v, path.dx, path.dy)
translate(path::OpenPolygonalPath, v::Vec) = open_polygonal_path(translate(path.vertices, v))
translate(path::ClosedPolygonalPath, v::Vec) = closed_polygonal_path(translate(path.vertices, v))
translate(path::OpenSplinePath, v::Vec) = open_polygonal_path(translate(path.vertices, v), path.v0, path.v1)
translate(path::ClosedSplinePath, v::Vec) = closed_polygonal_path(translate(path.vertices, v))
translate(ps::Locs, v::Vec) = map(p->p+v, ps)

path_length(path::CircularPath) = 2*pi*path.radius
path_length(path::ArcPath) = path.radius*abs(path.amplitude)
path_length(path::RectangularPath) = 2*(path.dx + path.dy)
path_length(path::OpenPolygonalPath) = path_length(path.vertices)
path_length(path::ClosedPolygonalPath) = path_length(path.vertices) + distance(path.vertices[end], path.vertices[1])
path_length(ps::Locs) =
  let p = ps[1]
      l = 0.0
    for i in 2:length(ps)
        pp = ps[i]
        l += distance(p, pp)
        p = pp
    end
    l
  end

location_at_length(path::CircularPath, d::Real) =
  loc_from_o_phi(path.center + vpol(path.radius, d/path.radius), d/path.radius + pi/2)
location_at_length(path::ArcPath, d::Real) =
  let Δα = d/path.radius,
      s = sign(path.amplitude)
    Δα <= abs(path.amplitude) + path_tolerance() ?
      loc_from_o_phi(path.center + vpol(path.radius, path.start_angle + Δα*s),
                     path.start_angle + (Δα + pi/2)*s) :
      error("Exceeded path length by ", Δα - path.amplitude)
  end
location_at_length(path::RectangularPath, d::Real) =
    let d = d % (2*(path.dx + path.dy)) # remove multiple periods
        p = path.corner
        for (delta, phi) in zip([path.dx, path.dy, path.dx, path.dy], [0, pi/2, pi, 3pi/2])
            if d - delta < path_tolerance()
                return loc_from_o_phi(add_pol(p, d, phi), phi)
            else
                p = add_pol(p, delta, phi)
                d -= delta
            end
        end
    end
location_at_length(path::OpenPolygonalPath, d::Real) =
    let p = path.vertices[1]
        for i in 2:length(path.vertices)
            pp = path.vertices[i]
            delta = distance(p, pp)
            if d - delta < path_tolerance()
                phi = pol_phi(pp - p)
                return loc_from_o_phi(add_pol(p, d, phi), phi)
            else
                p = pp
                d -= delta
            end
        end
        error("Exceeded path length by ", d)
    end
location_at_length(path::ClosedPolygonalPath, d::Real) =
    let p = path.vertices[1]
        for i in countfrom(1)
            pp = path.vertices[i%length(path.vertices)+1]
            delta = distance(p, pp)
            if d - delta < path_tolerance()
                phi = pol_phi(pp - p)
                return loc_from_o_phi(add_pol(p, d, phi), phi)
            else
                p = pp
                d -= delta
            end
        end
    end

subpath(path::CircularPath, a::Real, b::Real) =
    arc_path(path.center, path.radius, a/path.radius, (b-a)/path.radius)
subpath(path::ArcPath, a::Real, b::Real) =
  let Δα = (b - a)/path.radius
    a/path.radius + Δα <= path.start_angle + path.amplitude + path_tolerance() ?
      arc_path(path.center, path.radius, a/path.radius + path.start_angle, Δα) :
      error("Exceeded path length by ", path.start_angle + path.amplitude - a/path.radius - Δα)
  end
subpath(path::RectangularPath, a::Real, b::Real) =
    subpath(convert(ClosedPolygonalPath, path), a, b)
subpath(path::ClosedPolygonalPath, a::Real, b::Real) =
    subpath(convert(OpenPolygonalPath, path), a, b)
subpath(path::OpenPolygonalPath, a::Real, b::Real) =
    subpath_starting_at(subpath_ending_at(path, b), a)

subpath_starting_at(path::OpenPolygonalPath, d::Real) =
    let pts = path.vertices
        p1 = pts[1]
        for i in 2:length(pts)
            p2 = pts[i]
            delta = distance(p1, p2)
            diff = d - delta
            if diff < 0
                mp = p1 + (p2 - p1)*d/delta
                return open_polygonal_path([mp, pts[i:end]...])
            elseif abs(diff) < 1e-9
                return open_polygonal_path(i == length(pts) ? [pts[i], pts[i]] : pts[i:end])
            else
                p1 = p2
                d -= delta
            end
        end
        abs(d) < path_tolerance() ?
          path :
          error("Exceeded path length by ", d)
    end

subpath_ending_at(path::OpenPolygonalPath, d::Real) =
    let pts = path.vertices
        p1 = pts[1]
        for i in 2:length(pts)
            p2 = pts[i]
            delta = distance(p1, p2)
            diff = d - delta
            if diff < 0
                mp = p1 + (p2 - p1)*d/delta
                return open_polygonal_path([pts[1:i-1]..., mp])
            elseif abs(diff) < path_tolerance()
                return open_polygonal_path(pts[1:i])
            else
              p1 = p2
              d -= delta
            end
        end
        abs(d) < path_tolerance() ?
          path :
          error("Exceeded path length by ", d)
    end

#=
collect_vertices_length(p::Loc, vs::Locs, d::Real) =
    let pp = vs[1]
        dist = distance(p, pp)
        d <= dist ?
        push!(Vector{Loc}(), p + (pp - p)/dist*d) :
        unshift!(collect_vertices_length(pp, vs[2:end], d - dist), pp)
    end

subpath(path::OpenPolygonalPath, a::Real, b::Real) =
    let p = path.vertices[1]
        pts = []
        for i in 2:length(path.vertices)
            pp = path.vertices[i]
            delta = distance(p, pp)
            if a <= delta
                phi = pol_phi(pp - p)
                p0 = add_pol(p, d, phi)
                return open_polygonal_path(
                        unshift!(collect_vertices_length(p0,
                                                         vcat(path.vertices[i:end], path.vertices),
                                                         delta_d),
                                 p0))
            else
                p = pp
                d -= delta
            end
        end
        error("Exceeded path length")
    end

=#

meta_program(p::OpenPolygonalPath) =
    Expr(:call, :open_polygonal_path, meta_program(p.vertices))

# Path can be made of subparts
abstract type PathOp end
#struct MoveToOp <: PathOp loc::Loc end
#struct MoveOp <: PathOp vec::Vec end
#struct LineToOp <: PathOp loc::Loc end
struct LineOp <: PathOp
  vec::Vec
end
#struct CloseOp <: PathOp end
struct ArcOp <: PathOp
    radius::Real
    start_angle::Real
    amplitude::Real
end
#struct LineToXThenToYOp <: PathOp loc::Loc end
#struct LineToYThenToXOp <: PathOp loc::Loc end
struct LineXThenYOp <: PathOp
  vec::Vec
end
struct LineYThenXOp <: PathOp
  vec::Vec
end

struct PathOps <: Path
  start::Loc
  ops::Vector{<:PathOp}
  closed::Bool
end
open_path_ops(start, ops...) = PathOps(start, [ops...], false)
closed_path_ops(start, ops...) = PathOps(start, [ops...], true)

path_length(path::PathOps) =
  let len = mapreduce(path_length, +, path.ops, init=0)
    path.closed ?
      len + distance(path.start, location_at_length(path, len)) :
      len
  end
path_length(op::LineOp) = norm(op.vec)
path_length(op::ArcOp) = op.radius*abs(op.amplitude)

translate(path::PathOps, v::Vec) =
  PathOps(path.start + v, path.ops, path.closed)

#translate_op(op::MoveToOp, v) = MoveToOp(op.loc + v)
#translate_op(op::LineToOp, v) = LineToOp(op.loc + v)
#translate_op(op::LineToXThenToYOp, v) = LineToXThenToYOpp(op.loc + v)
#translate_op(op::LineToYThenToXOp, v) = LineToYThenToXOp(op.loc + v)
#translate_op(op::PathOp, v) = op

location_at_length(path::PathOps, d::Real) =
  let ops = path.ops,
      start = path.start
    for op in ops
      delta = path_length(op)
      if d < delta + path_tolerance()
        return location_at_length(op, start, d)
      else
        start = location_at_length(op, start, delta)
        d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end

location_at_length(op::LineOp, start::Loc, d::Real) =
  start + op.vec*d/path_length(op)
location_at_length(op::ArcOp, start::Loc, d::Real) =
  let center = start - vpol(op.radius, op.start_angle),
      a = d/op.radius*sign(op.amplitude)
    center + vpol(op.radius, op.start_angle + a)
  end

subpath_starting_at(path::PathOps, d::Real) =
  let ops = path.ops,
      start = location_at_length(path, d)
    for i in 1:length(ops)
      op = ops[i]
      delta = path_length(op)
      if d == delta
        return PathOps(start, ops[i+1:end], false)
      elseif d < delta + path_tolerance()
        op = subpath_starting_at(op, d)
        return PathOps(start, [op, ops[i+1:end]...], false)
      else
        d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end
subpath_ending_at(path::PathOps, d::Real) =
  let ops = path.ops
    for i in 1:length(ops)
      op = ops[i]
      delta = path_length(op)
      if d == delta
        return PathOps(path.start, ops[1:i], false)
      elseif d < delta + path_tolerance()
        return PathOps(path.start, [ops[1:i-1]..., subpath_ending_at(op, d)], false)
      else
        d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end

subpath_starting_at(pathOp::LineOp, d::Real) =
  let len = path_length(pathOp.vec)
    LineOp(pathOp.vec*(len-d)/len)
  end
subpath_starting_at(pathOp::ArcOp, d::Real) =
  let a = d/pathOp.radius
    ArcOp(pathOp.radius, pathOp.start_angle + a, pathOp.amplitude - a)
  end

subpath_ending_at(pathOp::LineOp, d::Real) =
  let len = path_length(pathOp.vec)
    LineOp(pathOp.vec*d/len)
  end
subpath_ending_at(pathOp::ArcOp, d::Real) =
  let a = d/pathOp.radius
    ArcOp(pathOp.radius, pathOp.start_angle, a)
  end

subpath(path::PathOps, a::Real, b::Real) =
  subpath_starting_at(subpath_ending_at(path, b), a)

# A path sequence is a sequence of paths where the next element of the sequence
# starts at the same place where the previous element ends.

struct OpenPathSequence <: OpenPath
  paths::Vector{<:Path}
end
open_path_sequence(paths...) =
  OpenPathSequence(ensure_connected_paths([paths...]))
struct ClosedPathSequence <: ClosedPath
  paths::Vector{<:Path}
end
closed_path_sequence(paths...) =
  ClosedPathSequence(ensure_connected_paths([paths...]))

PathSequence = Union{OpenPathSequence, ClosedPathSequence}

path_sequence(paths...) =
  coincident_path_location(path_start(paths[1]), path_end(paths[end])) ?
    closed_path_sequence(paths...) :
    open_path_sequence(paths...)

ensure_connected_paths(paths) = # AML: Finish this
    paths

translate(path::T, v::Vec) where T<:PathSequence =
  T(translate.(path.paths, v))
path_length(path::PathSequence) =
  sum(map(path_length, path.paths))
location_at_length(path::PathSequence, d::Real) =
  let paths = path.paths
    for path in paths
      delta = path_length(path)
      if d <= delta + path_tolerance()
          return location_at_length(path, d)
      else
          d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end

import Base.convert
convert(::Type{PathOps}, path::PathSequence) =
  let start = in_world(location_at_length(path, 0)),
      ops = mapfoldl(convert_to_path_ops, vcat, path.paths, init=PathOp[])
    PathOps(start, ops, coincident_path_location(start, path_end(path)))
  end
convert_to_path_ops(path::ArcPath) =
  [ArcOp(path.radius, path.start_angle, path.amplitude)]
convert_to_path_ops(path::PolygonalPath) =
  [LineOp(in_world(p)-in_world(q)) for (p,q) in zip(path.vertices[[2:end;1]], path.vertices)]
convert_to_path_ops(path::RectangularPath) =
  let p0 = in_world(path.corner),
      p1 = in_world(add_x(path.corner, path.dx)),
      p2 = in_world(add_xy(path.corner, path.dx, path.dy)),
      p3 = in_world(add_y(path.corner, path.dy))
    [LineOp(p1-p0), LineOp(p2-p1), LineOp(p3-p2), LineOp(p0-p3)]
  end

# A path set is a set of independent paths.

struct PathSet <: Path
  paths::Vector{<:Path}
end

# Should we just use tuples instead of arrays?
path_set(paths...) =
  PathSet([paths...])

# Convertions from/to paths
import Base.convert
convert(::Type{OpenPath}, vs::Locs) = open_polygonal_path(vs)
convert(::Type{ClosedPath}, vs::Locs) = closed_polygonal_path(vs)
convert(::Type{Path}, vs::Locs) =
  if vs[1] == vs[end]
    closed_polygonal_path(vs[1:end-1])
  else
    open_polygonal_path(vs)
  end
convert(::Type{ClosedPath}, p::OpenPath) =
  if isa(p.ops[end], CloseOp) || coincident_path_location(path_start(p), path_end(p))
    closed_path(p.ops)
  else
    error("Can't convert to a Closed Path: $p")
  end
convert(::Type{ClosedPath}, ops::Vector{<:PathOp}) =
  if isa(ops[end], CloseOp)
    closed_path(ops)
  else
    error("Can't convert to a Closed Path: $ops")
  end
convert(::Type{OpenPath}, ops::Vector{<:PathOp}) = open_path(ops)
convert(::Type{Path}, ops::Vector{<:PathOp}) =
  if isa(ops[end], CloseOp)
    closed_path(ops)
  else
    open_path(ops)
  end
convert(::Type{ClosedPolygonalPath}, path::RectangularPath) =
  let p = path.corner,
      dx = path.dx,
      dy = path.dy
    closed_polygonal_path([p, add_x(p, dx), add_xy(p, dx, dy), add_y(p, dy)])
  end
# The next one is looses precision. It converts from a circular to a polygonal path
# Note that the level of detail is hardwired. Maybe it will be good to make this a parameter
convert(::Type{ClosedPolygonalPath}, path::CircularPath) =
  let c = path.center,
      r = path.radius
    closed_polygonal_path([c + vpol(r, phi) for phi in division(0, 2pi, 64, false)])
  end
convert(::Type{OpenPolygonalPath}, path::ArcPath) =
  let c = path.center,
      r = path.radius
    open_polygonal_path([c + vpol(r, phi) for phi in division(path.start_angle, path.amplitude, 32, true)])
  end
convert(::Type{OpenPolygonalPath}, path::ClosedPolygonalPath) =
  open_polygonal_path(vcat(path.vertices, [path.vertices[1]]))
convert(::Type{OpenPolygonalPath}, path::RectangularPath) =
  convert(OpenPolygonalPath, convert(ClosedPolygonalPath, path))
convert(::Type{OpenPolygonalPath}, path::OpenSplinePath) = # ERROR: ignores limit vectors
  let interpolator = curve_interpolator(path.vertices),
      fixed_normal(vn, vt) = norm(vn) < path_tolerance() ? SVector{3}(vpol(1, sph_phi(xyz(vt[1],vt[2],vt[3], world_cs))+pi/2).raw[1:3]) : vn
    open_polygonal_path(
      map_division(
        t-> let p = interpolator(t),
                vt = Interpolations.gradient(interpolator, t)[1],
                vn = fixed_normal(Interpolations.hessian(interpolator, t)[1], vt)
                vy = cross(vt, vn)
              loc_from_o_vx_vy(
                xyz(p[1], p[2], p[3], world_cs),
                vxyz(vn[1], vn[2], vn[3], world_cs),
                vxyz(vy[1], vy[2], vy[3], world_cs))
            end,
        0.0, 1.0, 64)) # HACK this must be parameterized!
    end

curve_interpolator(pts::Locs) =
    let pts = map(pts) do p
                let v = in_world(p).raw
                  SVector{3,Float64}(v[1], v[2], v[3])
                end
              end
        Interpolations.scale(
            interpolate(pts, BSpline(Cubic(Natural(OnGrid())))),
            range(0,stop=1,length=size(pts, 1)))
    end

convert(::Type{ClosedPolygonalPath}, path::ClosedPathSequence) =
  let paths = path.paths,
      vertices = []
    for path in paths
      append!(vertices, convert(OpenPolygonalPath, path).vertices[1:end-1])
    end
    closed_polygonal_path(vertices)
  end

convert(::Type{OpenPolygonalPath}, path::OpenPathSequence) =
  let paths = path.paths,
      vertices = []
    for path in paths[1:end-1]
      append!(vertices, convert(OpenPolygonalPath, path).vertices[1:end-1])
    end
    append!(vertices, convert(OpenPolygonalPath, paths[end]).vertices)
    open_polygonal_path(vertices)
  end

#### Utilities
export path_vertices, subpaths, subtract_paths

path_vertices(path::OpenPolygonalPath) = path.vertices
path_vertices(path::ClosedPolygonalPath) = path.vertices
path_vertices(path::OpenSplinePath) = path.vertices
path_vertices(path::ClosedSplinePath) = path.vertices
path_vertices(path::Path) = path_vertices(convert(ClosedPolygonalPath, path))

subpaths(path::OpenPolygonalPath) =
  let ps = path.vertices
    map((p0, p1)->open_polygonal_path([p0, p1]), ps[1:end-1], ps[2:end])
  end
subpaths(path::ClosedPolygonalPath) =
  let ps = path.vertices
      map((p0, p1)->open_polygonal_path([p0, p1]), ps, [ps[2:end]..., ps[1]])
  end
subpaths(path::Path) = subpaths(convert(OpenPolygonalPath, path))

subtract_paths(path1::ClosedPolygonalPath, path2::ClosedPolygonalPath) =
  closed_polygonal_path(
    subtract_polygon_vertices(
      path_vertices(path1),
      path_vertices(path2)))


text:
using Documenter

push!(LOAD_PATH, "../src/")

using CGcoefficient

makedocs(
    sitename="CGcoefficient",
    pages = [
        "Home" => "index.md",
        "Formula" => "formula.md",
        "API" => "api.md"
    ]
)


text:
module ExtendableGrids

using DocStringExtensions
using ElasticArrays
using AbstractTrees

using InteractiveUtils
using SparseArrays
using Printf
using Random
using Dates
using LinearAlgebra

include("adjacency.jl")
export Adjacency,VariableTargetAdjacency,FixedTargetAdjacency
export atranspose,num_targets,num_sources,num_links,append!, max_num_targets_per_source
export asparse,tryfix,makevar

include("vectorofconstants.jl")
export VectorOfConstants

include("typehierarchy.jl")
export AbstractExtendableGridApexType
export typehierarchy

include("elementgeometry.jl")
export AbstractElementGeometry, ElementInfo
export elementgeometries

export AbstractElementGeometry0D
export Vertex0D

export AbstractElementGeometry1D
export Edge1D

export AbstractElementGeometry2D
export Polygon2D,Triangle2D,Quadrilateral2D,Pentagon2D,Hexagon2D,Parallelogram2D,Circle2D

export AbstractElementGeometry3D
export Polyhedron3D,Tetrahedron3D, Hexahedron3D,Parallelepiped3D,Prism3D,TrianglePrism3D,Sphere3D

export AbstractElementGeometry4D
export Polychoron4D,HyperCube4D

export dim_element

include("coordinatesystem.jl")
export coordinatesystems
export AbstractCoordinateSystem
export Cartesian1D,Cartesian2D,Cartesian3D
export Cylindrical2D,Cylindrical3D
export Polar2D,Polar1D ,Spherical3D,Spherical1D

include("extendablegrid.jl")
export ExtendableGrid
export instantiate, veryform
export AbstractGridComponent
export AbstractGridAdjacency,AbstractElementGeometries,AbstractElementRegions
export Coordinates,CellNodes,BFaceNodes,CellGeometries,BFaceGeometries,CellRegions,BFaceRegions
export NumCellRegions,NumBFaceRegions,CoordinateSystem
export AbstractGridFloatArray1D,AbstractGridFloatArray2D
export AbstractGridIntegerArray1D,AbstractGridIntegerArray2D
export index_type, coord_type
export dim_space, dim_grid
export num_nodes, num_cells, num_bfaces, num_cellregions, num_bfaceregions
export gridcomponents
export seemingly_equal

include("subgrid.jl")
export subgrid

include("more.jl")
export EdgeNodes, CellEdges,EdgeCells,BFaceCells,BFaceNormals
export local_celledgenodes,num_edges

include("regionedit.jl")
export cellmask!,bfacemask!

include("simplexgrid.jl")
export simplexgrid, geomspace,glue
export XCoordinates, YCoordinates, ZCoordinates
export writefile

include("tokenstream.jl")
export TokenStream, gettoken, expecttoken,trytoken

end # module


text: abstract type Topology <: Component end


text:
# Tests of this file not covered in other tests
@testset "dual_model_variables.jl" begin

    @testset "push_to_dual_obj_aff_terms!" begin
        primal_model = soc1_test()
        dual_obj_affine_terms = Dict{VI, Float64}()
        ci = CI{VVF, MOI.SecondOrderCone}(2)
        Dualization.push_to_dual_obj_aff_terms!(primal_model, dual_obj_affine_terms, VI(1), ci, 1)
        @test isempty(dual_obj_affine_terms)
    end

    @testset "set_dual_variable_name" begin
        primal_model = soc1_test()
        vi = VI(1)
        Dualization.set_dual_variable_name(primal_model, vi, 1, "con", "")
        @test MOI.get(primal_model, MOI.VariableName(), vi) == "con_1"
        Dualization.set_dual_variable_name(primal_model, vi, 2, "con", "")
        @test MOI.get(primal_model, MOI.VariableName(), vi) == "con_2"
        Dualization.set_dual_variable_name(primal_model, vi, 2, "con", "oi")
        @test MOI.get(primal_model, MOI.VariableName(), vi) == "oicon_2"
    end
end


text:
@testset "board" begin
    brd = board()
    @test startswith(brd.name, "brd_")
    @test brd isa J.Board
    @test brd.xlim == [-10.,10.]
    @test brd.ylim == [-10.,10.]
    @test !brd.axis
    @test !brd.showcopyright
    @test !brd.shownavigation

    brd = board("foo", xlim=[-2,2], ylim=[3,4], axis=true)
    @test brd.xlim == [-2,2]
    @test brd.ylim == [3,4]
    @test brd.axis

    @test_throws AssertionError board(boundingbox=[1,2])

    # --- Get options
    o = J.get_opts(brd)
    @test o.boundingbox == [-2, 4, 2, 3]
    @test o.axis
    @test !o.shownavigation
end

@testset "board-add" begin
    brd = board()
    s1 = slider("a", [[1,-1],[5,-1],[0,1.5,3]])
    s2 = slider("b", [[1,-1],[5,-1],[0,1.5,3]])
    (s1, s2) |> brd
    @test length(brd.objects) == 2
    @test brd.objects[1].name == "a"
    @test brd.objects[2].name == "b"

    brd = board()
    brd ++ s1
    @test length(brd.objects) == 1
    @test brd.objects[1].name == "a"

    brd = board()
    brd ++ (s1, s2)
    @test length(brd.objects) == 2
end

@testset "board-str" begin
    brd = board("brd")
    s1 = slider("a", [[1,-1],[5,-1],[0,1.5,3]])
    brd ++ s1

    s = J.str(brd)
    @test occursin("function val(x)", s)
    @test occursin(".initBoard(\"brd\",", s)
    @test occursin(".create('slider',", s)

    # -- save
    tmp = tempname()
    J.save(brd, tmp)
    s = read(tmp * ".js", String)
    @test s == J.str(brd)

    # -- standalone
    h = J.standalone(brd, full=true)
    @test occursin(s, h)
end


text:
print_bool=false

for experiment=1:size(filenameIBS,1)
    #############
    # Treduce=12
    # Vreduce=Int(Treduce/Ts_IBS)
    IBS[filenameIBS[experiment,1]]=IBS[filenameIBS[experiment,1]][Vreduce:end,:]

    if print_bool
        println(string("Reduce IBS finished, experiment : "),experiment)
    end

end

println("Reduce IBS finished")


text:
"""
$(SIGNATURES)

Helper function to process an individual block when the block is a `HCond` such as `{{ if
showauthor }} {{ fill author }} {{ end }}`.
"""
function convert_html_block(β::HCond, allvars::PageVars)::String
    # check that the relevant bool variable exists
    allconds = (β.init_cond, β.sec_conds...)
    if any(c -> !haskey(allvars, c), allconds)
        throw(HTMLBlockError("At least one of the booleans in a conditional html block" *
                             "could not be found. The block is reproduced below:" *
                             "\n" * "-"^50 * "\n" * β.ss * "\n" * "-"^50 * "\n"))
    end
    # check if there's an "else" clause
    has_else = (length(β.actions) == 1 + length(β.sec_conds) + 1)
    # check the first clause that is verified
    k = findfirst(c -> allvars[c].first, allconds)
    # if none is verified, use the else clause if there is one or do nothing
    if isnothing(k)
        has_else || return ""
        partial = β.actions[end]
    # otherwise run the 1st one which is verified
    else
        partial = β.actions[k]
    end
    # The String(...) is necessary here as to avoid problematic indexing further on
    return convert_html(String(partial), allvars)
end

"""
$(SIGNATURES)

Helper function to process an individual block when the block is a `HIsDef` such as `{{ isdef
author }} {{ fill author }} {{ end }}`. Which checks if a variable exists and if it does, applies
something.
"""
function convert_html_block(β::HCondDef, allvars::PageVars)::String
    hasvar = haskey(allvars, β.vname)
    # check if the corresponding bool is true and if so, act accordingly
    doaction = ifelse(β.checkisdef, hasvar, !hasvar)
    doaction && return convert_html(String(β.action), allvars)
    # default = do nothing
    return ""
end

"""
$(SIGNATURES)

Helper function to process an individual block when the block is a `HIsPage` such as `{{ ispage
path/to/page}} ... {{end}}`. Which checks if the current page is a given one and applies something
if that's the case (useful to handle different layouts on different pages).
"""
function convert_html_block(β::HCondPage, allvars::PageVars)::String
    # current path is relative to /src/ for instance /src/pages/blah.md -> pages/blah.md
    rpath = CUR_PATH[]
    # replace the `pages/` by `pub/`
    rpath = replace(rpath, Regex("^pages$(escape_string(PATH_SEP))") => "pub$(PATH_SEP)")
    # remove the extension
    rel_path = splitext(rpath)[1]
    # compare with β.pnames
    inpage = any(page -> splitext(page)[1] == rel_path, β.pages)
    # check if the corresponding bool is true and if so, act accordingly
    doaction = ifelse(β.checkispage, inpage, !inpage)
    doaction && return convert_html(String(β.action), allvars)
    # default = do nothing
    return ""
end


text:
struct Coordinate{N, T}
    coordinates::NTuple{N, T}
    Coordinate{N}(xs::NTuple{N, T}) where {N, T} = new{N, T}(xs)
    Coordinate(xs::NTuple{N, T}) where {N, T} = new{N, T}(xs)
    Coordinate(xs::T...) where T = new{length(xs), T}(xs)
end

for op in [:+, :-]
    @eval Base.$op(x::Coordinate{N}, y::Coordinate{N}) where N =
        Coordinate(map($op, x.coordinates, y.coordinates))

    @eval Base.$op(x::Coordinate{N}, y::AbstractVector) where N = $op(x, Coordinate(y...))
    @eval Base.$op(x::AbstractVector, y::Coordinate{N}) where N = $op(Coordinate(x...), y)
end

==(a::Coordinate{N, T}, b::Coordinate{N, T}) where {N, T} =
    (a.coordinates === b.coordinates)
==(a::Coordinate, b::Coordinate) = false

Base.:*(x::Coordinate, y::Int) = Coordinate(map(c -> y*c, x.coordinates))
Base.:*(x::Int, y::Coordinate) = y * x

Base.:-(x::Coordinate) = -1 * x

function Base.show(io::IO, x::Coordinate)
    print(io, "Coordinate(", join(x.coordinates, ", "), ")")
end


text:
module ECCodes

import Base: close,keys,getindex
export withmessages

const depfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("libeccodes not properly installed. Please run Pkg.build(\"ECCodes\")")
end

type GribFile
    file::Ptr{Void}
end

function GribFile(fname::String)
    f=ccall((:fopen,"libc"),
            Ptr{Void},(Cstring,Cstring),
            fname,"r")
    if f==C_NULL
        error("Can't open grib file $(fname).")
    end
    return GribFile(f)
end

function close(f::GribFile)
    if f.file!=C_NULL
        errno=ccall((:fclose,"libc"),
                    Cint,(Ptr{Void},),
                    f.file)
        if errno!=0
            error("Error $(errno) when calling fclose.")
        end
        f.file=C_NULL
    end
end

type GribMessage
    handle::Ptr{Void}
end

valid(m::GribMessage)=m.handle!=C_NULL

function nextmessage(f::GribFile)
    errref=Ref{Cint}(0)
    m=ccall((:codes_grib_handle_new_from_file,libeccodes),
            Ptr{Void},(Ptr{Void},Ptr{Void},Ref{Cint}),
            C_NULL,f.file,errref)
    if errref[]!=0
        error("Error $(errref[]) when calling codes_grib_handle_new_from_file.")
    end
    return GribMessage(m)
end

function freemessage(m::GribMessage)
    if m.handle!=C_NULL
        ccall((:codes_handle_delete,libeccodes),
              Cint,(Ptr{Void},),
              m.handle)
        m.handle=C_NULL
    end
end

function keys(m::GribMessage,filter_flags=0)
    BADKEYS=["codedValues","values","bitmap","reservedNeedNotBePresent"]

    r=String[]
    kiter=ccall((:codes_keys_iterator_new,libeccodes),
                Ptr{Void},(Ptr{Void},Culong,Ptr{Void}),
                m.handle,filter_flags,C_NULL)
    while true
        flag=ccall((:codes_keys_iterator_next,libeccodes),
                   Cint,(Ptr{Void},),
                   kiter)
        if flag==0
            break
        end
        kname=ccall((:codes_keys_iterator_get_name,libeccodes),
                    Cstring,(Ptr{Void},),
                    kiter)
        key=unsafe_string(kname)
        if !(key in BADKEYS)
            push!(r,key)
        end
    end
    errno=ccall((:codes_keys_iterator_delete,libeccodes),
                Cint,(Ptr{Void},),
                kiter)
    if errno!=0
        error("Error $(errno) when calling codes_keys_iterator_delete.")
    end
    return r
end

function getstring(m::GribMessage,key::String)
    if !valid(m)
        error("Attempt to read from invalid message.")
    end
    length=Ref{Csize_t}(0)
    ccall((:codes_get_length,libeccodes),
          Cint, (Ptr{Void},Cstring,Ref{Csize_t}),
          m.handle,key,length)
    value=zeros(UInt8,length[]+1) # hopefully this buffer is big enough
    errno=ccall((:codes_get_string,libeccodes),
                Cint, (Ptr{Void},Cstring,Ptr{UInt8},Ref{Csize_t}),
                m.handle,key,value,length)
    if errno==-10
        throw(KeyError(key))
    end
    if errno!=0
        error("Error $(errno) when calling codes_get_string.")
    end
    return unsafe_string(pointer(value))
end

function data(m::GribMessage)
    errref=Ref{Cint}(0)
    iterator=ccall((:codes_grib_iterator_new,libeccodes),
                   Ptr{Void},(Ptr{Void},Culong,Ref{Cint}),
                   m.handle,0,errref)
    lat=Ref{Cdouble}(0)
    lon=Ref{Cdouble}(0)
    val=Ref{Cdouble}(0)
    lats=Float64[]
    lons=Float64[]
    vals=Float64[]
    while true
        ok=ccall((:codes_grib_iterator_next,libeccodes),
                 Cint,(Ptr{Void},Ref{Cdouble},Ref{Cdouble},Ref{Cdouble}),
                 iterator,lat,lon,val)
        if ok==0
            break
        end
        push!(lats,lat[])
        push!(lons,lon[])
        push!(vals,val[])
    end
    ccall((:codes_grib_iterator_delete,libeccodes),
          Cint,(Ptr{Void},),
          iterator)
    return (lats,lons,vals)

end

getindex(m::GribMessage,key::String)=getstring(m,key)

function withmessages(dofunc,fname::String)
    f=GribFile(fname)
    while true
        m=nextmessage(f)
        if !valid(m)
            break
        end
        dofunc(m)
        freemessage(m)
    end
    close(f)
end

end


text:
export data_properties_struct, get_data_properties, tod2map_mpi, baseline2map_mpi, destripe, baselines_covmat

using LinearAlgebra

try
    import MPI
catch
end

@doc raw"""
This structure holds a number of parameters relative to the noise level and
the 1/f baselines measured or simulated for a certain polarimeter.

Field               | Type           | Meaning
:-----------------  |:-------------- |:----------------------------------------------------------------------------------
`pol_number`        | Int            | ID number of the polarimeter
`σ`                 | Float          | σ of white noise
`num_of_samples`    | Int            | total number of samples
`num_of_baselines`  | Int            | total number of 1/f baselines in the measured or simulated TOD for this polarimeter
`baselines_lengths` | Array{Int}     | Array containing the length of each 1/f baseline for this polarimeter

# Example
data_properties(1, 0.002463, 500, 5, [100, 100, 100, 100, 100, 100])

says that we have 5 baselines simulated for polarimeter number 1 (in this case polarimeter STRIP31, with σ = 0.002463), each of length 100 samples.
"""
struct data_properties_struct
    polarimeter::Int
    σ::Float64
    number_of_samples::Int
    number_of_baselines::Int
    baselines_lengths::Array{Int}
end

@doc raw"""
    function get_data_properties(detector_number, σ_k, num_of_baselines, num_of_samples) -> data_properties

This function can be used to build the object `data properties` needed for desriping and calculation of baselines_covmat.

It requires in input 4 arrays containing:
- the ID number of the polarimeters that the current rank will simulate
- the white noise σ for each polarimeter that the current rank will simulate
- the number of 1/f baselines for each polarimeter that the current rank will simulate
- the total number of samples for each polarimeter that the current rank will simulate

N.B. the ID numbers, the number of 1/f baselines and the total number of samples can be obtained using the function `get_chunk_properties`.

It returns an array of `data_properties_struct`, of length equal to the number of polarimeters simulated by current rank.

"""
function get_data_properties(detector_number, σ_k, num_of_baselines, num_of_samples)
    data_properties  = Array{data_properties_struct}(undef, length(detector_number))
    for i in 1:length(detector_number)
        baselines_lengths  = repeat([Int64(num_of_samples[i]/num_of_baselines[i])], num_of_baselines[i])  #we suppose baselines of equal lengths
        data_properties[i]  = data_properties_struct(detector_number[i], σ_k[detector_number[i]], num_of_samples[i], num_of_baselines[i], baselines_lengths)
    end
    data_properties
end

@doc raw"""
    tod2map(pix_idx, tod, num_of_pixels, comm) -> binned_map
    tod2map(pix_idx, tod, num_of_pixels, data_properties, comm) -> binned_map

This function creates a binned map from the time-ordered data kept in
the array `tod`, assuming that each sample is observing the pixel
whose index is in `pix_idx`. The parameter `num_of_pixels` contains
the number of pixels in the Healpix map, and it is used as an upper
bound for the values in `pix_idx`. The parameter `comm` must be a MPI
communicator, or `missing` if you are not using MPI.
This is a MPI based function: each MPI process computes a map from its
available data.  All partial maps are then combined together with MPI.allreduce.
The function returns an array containing the binned map.

If Array of structures `data_properties_struct`is passed to the function, the output binned map will be a weighted binned map.
Each sample will be weighted according to the inverse white noise variance σ^2 of the corrisponding polarimeter.
In this way, the less noisy polarimeters will count more in the estimation of the map.

# Requirements
- The length of the arrays `pix_idx` and `tod` must be the same

"""
tod2map_mpi

function tod2map_mpi(pix_idx, tod, num_of_pixels, comm; unseen = NaN)
    T = eltype(tod)
    N = eltype(pix_idx)

    partial_map = zeros(T, num_of_pixels)
    partial_hits = zeros(N, num_of_pixels)
    binned_map = zeros(T, num_of_pixels)
    hits = zeros(N, num_of_pixels )

    for i in eachindex(pix_idx)
        partial_map[pix_idx[i]] += tod[i]
        partial_hits[pix_idx[i]] += 1
    end

    if(!ismissing(comm))
        binned_map = MPI.allreduce(partial_map, MPI.SUM, comm)
        hits = MPI.allreduce(partial_hits, MPI.SUM, comm)
    else
        binned_map .= partial_map
        hits .= partial_hits
    end

    @inbounds for i in eachindex(binned_map)
        if (hits[i] > 0)
            binned_map[i] = binned_map[i] / hits[i]
        else
            binned_map[i] = unseen
        end
    end

    binned_map
end

function tod2map_mpi(pix_idx, tod, num_of_pixels, data_properties, comm; unseen = NaN)
    T = eltype(tod)

    partial_map = zeros(T, num_of_pixels)
    partial_hits = zeros(T, num_of_pixels)
    binned_map = zeros(T, num_of_pixels)
    hits = zeros(T, num_of_pixels )

    start_idx = 1

    for j in eachindex(data_properties)                 #loop on detectors
        end_idx = start_idx + data_properties[j].number_of_samples-1
        for i in start_idx:end_idx             #loop on samples
            partial_map[pix_idx[i]] += tod[i] * 1/(data_properties[j].σ)^2
            partial_hits[pix_idx[i]] += 1/(data_properties[j].σ)^2
        end
        start_idx += data_properties[j].number_of_samples
    end

    if(!ismissing(comm))
        binned_map = MPI.allreduce(partial_map, MPI.SUM, comm)
        hits = MPI.allreduce(partial_hits, MPI.SUM, comm)
    else
        binned_map .= partial_map
        hits .= partial_hits
    end

    @inbounds for i in eachindex(binned_map)
        if (hits[i] > 0)
            binned_map[i] = binned_map[i] / hits[i]
        else
            binned_map[i] = unseen
        end
    end

    binned_map
end

@doc raw"""
    baseline2map_mpi(pix_idx, baselines, baseline_lengths, num_of_pixels, comm)-> noise_map
    baseline2map_mpi(pix_idx, baselines, num_of_pixels, data_properties, comm) -> noise_map

This is a MPI based function: each MPI process computes a map from its
available data.  All partial maps are then combined together with
MPI.allreduce.
The function returns an array containing the binned map.

If Array of structures `data_properties_struct`is passed to the function (instead of `baseline_lengths`) the output binned map will be a weighted binned map.
Each sample will be weighted according to the inverse white noise variance σ^2 of the corrisponding polarimeter.
In this way, the less noisy polarimeters will count more in the estimation of the map.

# Requirements
- The length of `baselines` and `baseline_lengths` must be the same;
- The value `sum(baseline_lengths)` must be the same as the length of `pix_idx`.
"""
baseline2map_mpi

function baseline2map_mpi(pix_idx, baselines, baseline_lengths, num_of_pixels, comm;
                          unseen=NaN)

    T = eltype(baselines)
    N = eltype(pix_idx)

    partial_map = zeros(T, num_of_pixels)
    partial_hits = zeros(N, num_of_pixels)
    noise_map = zeros(T, num_of_pixels)
    hits = zeros(N, num_of_pixels)

    startidx = 1

    for i in eachindex(baseline_lengths)
        endidx = baseline_lengths[i] + startidx - 1

        for j in startidx:endidx
            partial_map[pix_idx[j]] += baselines[i]
            partial_hits[pix_idx[j]] += 1
        end
        startidx += baseline_lengths[i]
    end

    if(!ismissing(comm))
        noise_map = MPI.allreduce(partial_map, MPI.SUM, comm)
        hits = MPI.allreduce(partial_hits, MPI.SUM, comm)
    else

        noise_map .= partial_map
        hits .= partial_hits
    end

    for i in eachindex(noise_map)
        if (hits[i] > 0)
            noise_map[i] /= hits[i]
        else
            noise_map[i] = unseen
        end
    end

    noise_map
end

function baseline2map_mpi(pix_idx, baselines, num_of_pixels, data_properties, num_of_baselines, comm; unseen=NaN)

    T = eltype(baselines)

    partial_map = zeros(T, num_of_pixels)
    partial_hits = zeros(T, num_of_pixels)
    noise_map = zeros(T, num_of_pixels)
    hits = zeros(T, num_of_pixels)

    startidx = 1
    baseline_idx = 1

    for l in eachindex(data_properties)  #loop on detectors
        for i in eachindex(data_properties[l].baselines_lengths)
            endidx = data_properties[l].baselines_lengths[i] + startidx - 1

            for j in startidx:endidx
                partial_map[pix_idx[j]] += baselines[baseline_idx] * 1/(data_properties[l].σ)^2
                partial_hits[pix_idx[j]] += 1/(data_properties[l].σ)^2
            end
            startidx += data_properties[l].baselines_lengths[i]
            baseline_idx += 1
        end
    end

    if(!ismissing(comm))
        noise_map = MPI.allreduce(partial_map, MPI.SUM, comm)
        hits = MPI.allreduce(partial_hits, MPI.SUM, comm)
    else

        noise_map .= partial_map
        hits .= partial_hits
    end

    for i in eachindex(noise_map)
        if (hits[i] > 0)
            noise_map[i] /= hits[i]
        else
            noise_map[i] = unseen
        end
    end

    noise_map
end

function applyz_and_sum(pix_idx, tod, num_of_pixels, data_properties, num_of_baselines, comm; unseen=NaN)

    @assert length(tod) == length(pix_idx)

    baselines_sum = zeros(eltype(tod), num_of_baselines)

    binned_map = tod2map_mpi(
        pix_idx,
        tod,
        num_of_pixels,
        data_properties,
        comm,
        unseen=unseen)

    baseline_idx = 1
    startidx = 1
    for l in eachindex(data_properties)  #loop on detectors

        for i in eachindex(data_properties[l].baselines_lengths)
            endidx = data_properties[l].baselines_lengths[i] + startidx - 1

            for j in startidx:endidx
                baselines_sum[baseline_idx] += (tod[j] - binned_map[pix_idx[j]]) * 1/(data_properties[l].σ)^2
            end

            startidx += data_properties[l].baselines_lengths[i]
            baseline_idx += 1
        end
    end

    baselines_sum
end

function applya(baselines, pix_idx, num_of_baselines, num_of_pixels, data_properties, comm; unseen=NaN)
    @assert length(baselines) == num_of_baselines

    baselines_sum = zeros(eltype(baselines), num_of_baselines)
    total_sum = zero(eltype(baselines))

    binned_map = baseline2map_mpi(
        pix_idx,
        baselines,
        num_of_pixels,
        data_properties,
        num_of_baselines,
        comm;
        unseen=unseen)

    startidx = 1
    baseline_idx = 1

    for l in eachindex(data_properties)

        for i in eachindex(data_properties[l].baselines_lengths)
            endidx = data_properties[l].baselines_lengths[i] + startidx - 1

            for j in startidx:endidx
                baselines_sum[baseline_idx] += (baselines[baseline_idx] - binned_map[pix_idx[j]]) * 1/(data_properties[l].σ)^2
            end

            startidx += data_properties[l].baselines_lengths[i]
            baseline_idx += 1
        end
    end

    #needed to assure that sum(baselines)==0

    if(!ismissing(comm))
        total_sum = MPI.allreduce([sum(baselines)], MPI.SUM, comm)[1]
    else
        total_sum = sum(baselines)
    end

    baselines_sum .+= total_sum
end

function mpi_dot_prod(x, y, comm)

    @assert eltype(x) == eltype(y)
    result = zero(eltype(x))
    local_sum::eltype(x) = dot(x, y)

    if(!ismissing(comm))
        result = MPI.allreduce([local_sum], MPI.SUM, comm)[1]
    else
        result = local_sum
    end

    return result
end

function conj_grad(baselines_sum, pix_idx, tod,
                   num_of_pixels, data_properties,
                   num_of_baselines, rank, comm;
                   threshold = 1e-9, max_iter=10000)

    T = eltype(tod)
    N = eltype(pix_idx)

    baselines = Array{T}(undef, num_of_baselines)
    r = Array{T}(undef, num_of_baselines)
    r_next = Array{T}(undef, num_of_baselines)
    p = Array{T}(undef, num_of_baselines)
    Ap = Array{T}(undef,num_of_baselines)

    baselines .= 0    #starting baselines
    k = zero(N)       #number of iterations
    convergence_parameter = zero(T)
    rdotr = zero(T)
    rdotr_next = zero(T)

    best_convergence_parameter = zero(T)
    best_baselines = zeros(T, num_of_baselines)
    best_k = zero(N)

    r = baselines_sum - applya(baselines, pix_idx, num_of_baselines, num_of_pixels, data_properties, comm)  #residual
    p .= r
    rdotr = mpi_dot_prod(r, r, comm)
    best_convergence_parameter = sqrt(rdotr)

    if(rdotr == 0)
        return best_baselines
    end

    while true

        Ap =  applya(p, pix_idx,  num_of_baselines, num_of_pixels, data_properties, comm)

        rdotr = mpi_dot_prod(r, r, comm)
        pdotAp = mpi_dot_prod(p, Ap, comm)

        alpha = rdotr / pdotAp
        @. baselines += alpha * p
        @. r_next = r - alpha * Ap

        rdotr_next = mpi_dot_prod(r_next, r_next, comm)

        convergence_parameter = sqrt(rdotr_next)

        if (convergence_parameter < best_convergence_parameter)
            best_convergence_parameter = convergence_parameter
            best_baselines .= baselines
            best_k = k
        end

        if convergence_parameter < threshold  break end
        if k  > max_iter   break end

        beta = rdotr_next/rdotr

        @. p = r_next + beta * p
        r .= r_next
        k += 1
    end

    return best_baselines
end

function destriped_map(baselines, pix_idx, tod, data_properties, num_of_pixels, num_of_baselines, comm; unseen=NaN)
    @assert length(tod) == length(pix_idx)
    tod2map_mpi(pix_idx, tod, num_of_pixels, data_properties, comm) - baseline2map_mpi(pix_idx, baselines, num_of_pixels, data_properties, num_of_baselines, comm)
end

@doc raw"""
    destripe(pix_idx, tod, num_of_pixels, data_properties, rank, comm; threshold = 1e-9, max_iter = 10000) -> (pixels, baselines)

This MPI based function creates a map from a TOD and removes both 1/f
and white noise, using the destriping technique.

The parameters passed to the function have the following meaning:

- `pix_idx`: array containing the indices of the pixels visited by the instrument

- `tod`: the values measured by the polarimeters for each pixel
  (either I, Q, or U)

- `num_of_pixels`: the number of pixels in the map to be
   produced. This is used as an upper limit for the values in `pix_idx`

- `data_properties`: an array of structures `data_properties_struct`,
   holding information on each simulated polarimeter noise level,
   number and length of 1/f baselines and total number of samples.
   It can be obtained by using function `get_data_properties`.

- `rank`: the rank of the current MPI process

- `comm`: the MPI communicator object.

The following arguments are optional:

- `threshold` is used by the conjugated-gradient algorithm. When the
   residual error of the iteration goes below this value, the iteration
   stops. The smaller the value, the more accurate the solution.

- `max_iter` is the maximum number of iterations to be executed in the
   conjugated-gradient algorithm. If the algorithm does not converge
   within this number of iterations, the process will quit without
   having reached the convergence threshold (see the `threshold`
   keyword above).

The function returns a 2-tuple containing the destriped map itself
(`Array{T,1}`, where `T` is the base type used for the TOD in the
argument`tod`) and an array containing the baselines.
Since it is not granted that the sequence of convergence parameters of
the conjugate gradient is monotonically decreasing, the code keeps the
lowest value of them and the corresponding array of baselines.  If the
loop ends because the maximum number of iterations has been reached,
this is the configuration that will be returned to the caller.

# Remarks
- The length of the arrays `pix_idx` and `tod` must be the same;
- If you are not using MPI, pass `missing` to the `comm` parameter.
"""
function destripe(pix_idx, tod, num_of_pixels, data_properties, rank, comm;
                 threshold = 1e-9, max_iter = 10000, unseen=NaN)

    num_of_baselines = 0
    for i in 1:length(data_properties) num_of_baselines  += data_properties[i].number_of_baselines end

    baselines_sum = applyz_and_sum(
            pix_idx,
            tod,
            num_of_pixels,
            data_properties,
            num_of_baselines,
            comm,
            unseen=unseen)

    baselines = conj_grad(
            baselines_sum,
            pix_idx,
            tod,
            num_of_pixels,
            data_properties,
            num_of_baselines,
            rank,
            comm;
            threshold = threshold,
            max_iter = max_iter)

    # once we have an estimate of the baselines, we can build the destriped map
    destr_map = destriped_map(
            baselines,
            pix_idx,
            tod,
            data_properties,
            num_of_pixels,
            num_of_baselines,
            comm,
            unseen=unseen)

    #check that sum(baselines) = 0
    if(!ismissing(comm))
        total_sum = MPI.allreduce([sum(baselines)], MPI.SUM, comm)[1]
    else
        total_sum = sum(baselines)
    end

    if rank==0
        println("The sum of baselines is: $total_sum")
    end

    (destr_map, baselines)
end

@doc raw"""
    baselines_covmat(polarimeters, σ_k, baseline_length_s, fsamp_hz, total_time) -> covariance_matrix

This function produces the covariance matrix of the 1/f baselines computed by the destriper.
As an approximation, we ignore nondiagonal terms (i.e. the correlations between different baselines).
The function output is therefore an array corresponding to the covariance matrix diagonal.
Another assumption is that the white noise variance stays constant over a given baseline.

The baseline error is computed according to equation 29 in https://arxiv.org/abs/0904.3623

The parameters passed to the function have the following meaning:
    -`polarimeters`: the array of polarimeters ID numbers
    -`σ_k`: the array of corresponding white noise σ (in K)
    -`baseline_length_s` = the length (in s) of each 1/f baseline
    -`fsamp_hz` = the sampling frequency (in Hz)
    -`total_time` = the duration (in s) of the observation

The function firstly computes an Array of structures `data_properties_struct`,
which matches the white noise variance σ with the number of 1/f baselines and their lengths, for each polarimeter.
Then it computes the covariance matrix according to these matches.
"""
function baselines_covmat(polarimeters, σ_k, baseline_length_s, fsamp_hz, total_time)

    ALL_data_properties = Array{data_properties_struct}(undef, length(polarimeters))
    baselines_per_pol =  Int64(total_time/baseline_length_s)
    for i in 1:length(polarimeters)
        baselines_lengths  = repeat([baseline_length_s*fsamp_hz], baselines_per_pol)
        ALL_data_properties[i]  = data_properties_struct(polarimeters[i], σ_k[i], sum(baselines_lengths), baselines_per_pol, baselines_lengths)
    end

    covariance_matrix = []
    for i in 1:length(ALL_data_properties)
        partial_covmat = Array{Float64}(undef,length(ALL_data_properties[i].baselines_lengths))

        for j in 1:length(ALL_data_properties[i].baselines_lengths)
            σ = ALL_data_properties[i].σ
            this_baseline_length = ALL_data_properties[i].baselines_lengths[j]

            partial_covmat[j] = σ^2/this_baseline_length
        end
        covariance_matrix = append!(covariance_matrix, partial_covmat)
    end
    return covariance_matrix
end


text:
abstract type MaternLikelihood{T<:Real} <: RegressionLikelihood{T} end

"""
    Matern3_2Likelihood(ρ::Real=1.0)

## Arguments
- `ρ::Real` : lengthscale

---

Matern 3/2 likelihood for regression:
see [wiki page](https://en.wikipedia.org/wiki/Mat%C3%A9rn_covariance_function)
---

For the analytical solution, it is augmented via:
```math
p(y|f,\\omega) = \\mathcal{N}(y|f,\\sigma^2\\omega)
```
Where ``\\omega \\sim \\mathcal{IG}(\\frac{\\nu}{2},\\frac{\\nu}{2})`` where ``\\mathcal{IG}`` is the inverse gamma distribution
See paper [Robust Gaussian Process Regression with a Student-t Likelihood](http://www.jmlr.org/papers/volume12/jylanki11a/jylanki11a.pdf)
"""
struct Matern3_2Likelihood{T<:Real,A<:AbstractVector{T}} <: MaternLikelihood{T}
    ρ::T
    c::A
    θ::A
    function Matern3_2Likelihood{T}(ρ::T) where {T<:Real}
        return new{T,Vector{T}}(ρ)
    end
    function Matern3_2Likelihood{T}(ρ::T, c²::A, θ::A) where {T<:Real,A<:AbstractVector{T}}
        return new{T,A}(ρ, c², θ)
    end
end

function Matern3_2Likelihood(ρ::T=1.0) where {T<:Real}
    return Matern3_2Likelihood{T}(ρ)
end

function implemented(
    ::Matern3_2Likelihood, ::Union{<:AnalyticVI,<:QuadratureVI,<:GibbsSampling}
)
    return true
end

function init_likelihood(
    likelihood::Matern3_2Likelihood{T},
    inference::AbstractInference{T},
    ::Int,
    nSamplesUsed::Int,
) where {T}
    if inference isa AnalyticVI || inference isa GibbsSampling
        Matern3_2Likelihood{T}(
            likelihood.ρ, abs2.(T.(rand(T, nSamplesUsed))), zeros(T, nSamplesUsed)
        )
    else
        Matern3_2Likelihood{T}(likelihood.ρ)
    end
end

function (l::Matern3_2Likelihood)(y::Real, f::Real)
    u = sqrt(3) * abs(y - f) / l.ρ
    return 4 * l.ρ / sqrt(3) * (one(T) + u) * exp(-u)
end

function Distributions.loglikelihood(l::Matern3_2Likelihood, y::Real, f::Real)
    u = sqrt(3) * abs(y - f) / l.ρ
    return log(4 * l.ρ / sqrt(3)) + log(one(T) + u) - u
end

function Base.show(io::IO, l::Matern3_2Likelihood{T}) where {T}
    return print(io, "Matern 3/2 likelihood (ρ = $(l.ρ))")
end

function compute_proba(
    l::Matern3_2Likelihood{T}, μ::AbstractVector{<:Real}, σ²::AbstractVector{<:Real}
) where {T<:Real}
    return μ, max.(σ², zero(T)) .+ 4 * l.ρ^2 / 3
end

## Local Updates ##
function local_updates!(
    l::Matern3_2Likelihood,
    y::AbstractVector,
    μ::AbstractVector{T},
    diag_cov::AbstractVector{T},
) where {T<:Real}
    l.c .= sqrt.(diag_cov + abs2.(μ - y))
    return l.θ .= 3.0 ./ (2.0 .* sqrt.(3) * l.c * l.ρ .+ 2 * l.ρ^2)
end

function sample_local!(l::Matern3_2Likelihood, y::AbstractVector, f::AbstractVector)
    l.c .= rand.(GeneralizedInverseGaussian.(3 / (2 * l.ρ^2), 2.0 .* abs2.(f - y), 1.5))
    set_ω!(l, l.c)
    return nothing
end

@inline ∇E_μ(l::Matern3_2Likelihood, ::AOptimizer, y::AbstractVector) = (2.0 * l.θ .* y,)
@inline ∇E_Σ(l::Matern3_2Likelihood, ::AOptimizer, ::AbstractVector) = (l.θ,)

## ELBO  ##

function expecLogLikelihood(
    l::Matern3_2Likelihood,
    i::AnalyticVI,
    y::AbstractVector,
    μ::AbstractVector,
    diag_cov::AbstractVector,
)
    return throw(error("Not implemented yet"))
end

function InverseGammaKL(l::Matern3_2Likelihood{T}) where {T}
    α_p = l.ν / 2
    β_p = α_p * l.σ^2
    return model.inference.ρ * sum(broadcast(InverseGammaKL, l.α, l.β, α_p, β_p))
end

## PDF and Log PDF Gradients ## (verified gradients)

@inline function ∇loglikehood(l::Matern3_2Likelihood, y::Real, f::Real)
    return 3.0 * (y - f) / (l.ρ * (abs(f - y) * sqrt(3) + l.ρ))
end

@inline function hessloglikelihood(l::Matern3_2Likelihood, y::Real, f::Real)
    return 3.0 / (l.ρ + sqrt(3) * abs(f - y))^2
end


text:
import NaiveNASflux
import NaiveNASflux: AbstractMutableComp, MutableLayer, LazyMutable, weights, bias, select, layer, mutate, hiddenweights, hiddenstate, state, outscale
import Flux: mapchildren

@testset "Mutable computation" begin

    @testset "Dense MutableLayer" begin

        m = MutableLayer(Dense(2,3))

        @test nin(m) == nin(m.layer) == 2
        @test nout(m) == nout(m.layer) == 3
        @test minΔninfactor(m) == 1
        @test minΔnoutfactor(m) == 1
        @test m([1.0, 2.0]) == m.layer([1.0, 2.0])

        m.layer = Dense(3,4)
        bias(m.layer)[1:end] = 1:4
        @test nin(m) == nin(m.layer) == 3
        @test nout(m) == nout(m.layer) == 4
        @test m([1.0, 2.0, 3.0]) == m.layer([1.0, 2.0, 3.0])

        inds = [1,3]
        Wexp, bexp = weights(m.layer)[:, inds], bias(m.layer)
        mutate_inputs(m, inds)
        assertlayer(m.layer, Wexp, bexp)

        inds = [1,2,4]
        Wexp, bexp = weights(m.layer)[inds, :], bias(m.layer)[inds]
        mutate_outputs(m, inds)
        assertlayer(m.layer, Wexp, bexp)

        inds = [1,-1, 2]
        Wexp = hcat(weights(m.layer)[:, 1], zeros(Float32, 3), weights(m.layer)[:, 2])
        mutate_inputs(m, inds)
        assertlayer(m.layer, Wexp, bexp)

        inds = [-1, 1, -1, 3, -1]
        Wexp = permutedims(hcat(zeros(Float32, 3), weights(m.layer)[1, :], zeros(Float32, 3), weights(m.layer)[3, :], zeros(Float32,3)))
        bexp = Float32[0, bias(m.layer)[1], 0, bias(m.layer)[3], 0]
        mutate_outputs(m, inds)
        assertlayer(m.layer, Wexp, bexp)
    end
    @testset "Convolutional layers" begin

        @testset "Conv MutableLayer" begin
            m = MutableLayer(Conv((2,3),(4=>5)))

            @test nin(m) == nin(m.layer) == 4
            @test nout(m) == nout(m.layer) == 5
            @test minΔninfactor(m) == 1
            @test minΔnoutfactor(m) == 1
            input = reshape(collect(Float32, 1:3*4*4), 3, 4, 4, 1)
            @test m(input) == m.layer(input)

            inds = [1,3]
            Wexp, bexp = weights(m.layer)[:,:,inds,:], bias(m.layer)
            mutate_inputs(m, inds)
            assertlayer(m.layer, Wexp, bexp)

            inds = [1,2,4]
            Wexp, bexp = weights(m.layer)[:,:,:,inds], bias(m.layer)[inds]
            mutate_outputs(m, inds)
            assertlayer(m.layer, Wexp, bexp)

            inds = [1,-1, 2]
            wsize = deleteat!(collect(size(weights(m.layer))), 3)

            # Nothing beats working in four dimensions...
            # Stack 3D arrays in a 4:th dimension and then swap dim 3 and 4
            Wexp = permutedims(cat(
            weights(m.layer)[:,:,1,:],
            zeros(Float32, wsize...),
            weights(m.layer)[:,:,2,:], dims=4), [1,2,4,3])

            mutate_inputs(m, inds)
            assertlayer(m.layer, Wexp, bexp)

            inds = [-1, 1, -1, 3, -1]
            wsize = deleteat!(collect(size(weights(m.layer))), 4)

            Wexp = cat(
            zeros(Float32, wsize...),
            weights(m.layer)[:,:,:,1],
            zeros(Float32, wsize...),
            weights(m.layer)[:,:,:,3],
            zeros(Float32,wsize...), dims=4)

            bexp = Float32[0, bias(m.layer)[1], 0, bias(m.layer)[3], 0]
            mutate_outputs(m, inds)
            assertlayer(m.layer, Wexp, bexp)
        end

        @testset "ConvTranspose MutableLayer" begin
            m = MutableLayer(ConvTranspose((2,3),(4=>5)))

            @test nin(m) == nin(m.layer) == 4
            @test nout(m) == nout(m.layer) == 5
            @test minΔninfactor(m) == 1
            @test minΔnoutfactor(m) == 1
            input = reshape(collect(Float32, 1:3*4*4), 3, 4, 4, 1)
            @test m(input) == m.layer(input)

            inputs = [1,3]
            outputs = [1,2,4]
            Wexp, bexp = weights(m.layer)[:,:,outputs,inputs], bias(m.layer)[outputs]
            mutate(m, inputs=inputs, outputs=outputs)
            assertlayer(m.layer, Wexp, bexp)
        end

        @testset "DepthwiseConv MutableLayer" begin
            m = MutableLayer(DepthwiseConv((2,2),(3=>6*3)))

            @test nin(m) == nin(m.layer) == 3
            @test nout(m) == nout(m.layer) == 18

            input = reshape(collect(Float32, 1:3*3*3), 3, 3, 3, 1)
            @test m(input) == m.layer(input)

            ins = [1, 3]
            wouts = [1, 2, 5, 6]
            outs = mapfoldl(i -> 2 * i .+ [-1, -0] ,vcat, wouts)
            Wexp, bexp = weights(m.layer)[:,:,wouts,ins], bias(m.layer)[outs]
            mutate(m, inputs=ins, outputs=outs)
            assertlayer(m.layer, Wexp, bexp)
            @test size(m(ones(Float32, 3,3,2,2)))[3:4] == (8, 2)

            ins = [1, 2, -1]
            outs = [1, 2, -1, -1, -1, -1, 3, 4, -1, -1, -1, -1]
            mutate(m, inputs=ins, outputs=outs)

            @test nin(m) == 3
            @test nout(m) == 12

            @test size(m(ones(Float32, 3,3,3,2)))[3:4] == (12, 2)
        end

        @testset "CrossCor MutableLayer" begin
            m = MutableLayer(CrossCor((2,3),(4=>5)))

            @test nin(m) == nin(m.layer) == 4
            @test nout(m) == nout(m.layer) == 5
            @test minΔninfactor(m) == 1
            @test minΔnoutfactor(m) == 1
            input = reshape(collect(Float32, 1:3*4*4), 3, 4, 4, 1)
            @test m(input) == m.layer(input)

            inputs = [1,3]
            outputs = [1,2,4]
            Wexp, bexp = weights(m.layer)[:,:,inputs, outputs], bias(m.layer)[outputs]
            mutate(m, inputs=inputs, outputs=outputs)
            assertlayer(m.layer, Wexp, bexp)
        end
    end

    @testset "Diagonal MutableLayer" begin
        m = MutableLayer(Flux.Diagonal(4))

        @test nin(m) == nin(m.layer) == nout(m) == nout(m.layer) == 4
        @test minΔninfactor(m) == 1
        @test minΔnoutfactor(m) == 1
        weights(m.layer)[1:end] = 1:4
        bias(m.layer)[1:end] = 1:4

        @test m(Float32[1,2,3,4]) == m.layer(Float32[1,2,3,4])

        inds = [1,3]
        Wexp, bexp = weights(m.layer)[inds], bias(m.layer)[inds]
        mutate_inputs(m, inds)
        assertlayer(m.layer, Wexp, bexp)

        inds = [-1, 2, -1]
        Wexp, bexp = Float32[0, weights(m.layer)[2], 0], Float32[0, bias(m.layer)[2], 0]
        mutate_outputs(m, inds)
        assertlayer(m.layer, Wexp, bexp)
    end

    @testset "Normalization layers" begin

        @testset "LayerNorm MutableLayer" begin
            m = MutableLayer(LayerNorm(3))

            @test nin(m) == nin(m.layer) == nout(m) == nout(m.layer) == 3
            @test minΔninfactor(m) == 1
            @test minΔnoutfactor(m) == 1
            weights(m.layer.diag)[1:end] = 1:3
            bias(m.layer.diag)[1:end] = 1:3

            @test m(Float32[1,2,3]) == m.layer(Float32[1,2,3])

            inds = [1,3]
            Wexp, bexp = weights(m.layer.diag)[inds], bias(m.layer.diag)[inds]
            mutate_inputs(m, inds)
            @test typeof(layer(m)) <: LayerNorm
            assertlayer(m.layer.diag, Wexp, bexp)

            inds = [-1, 2, -1]
            Wexp, bexp = Float32[0, weights(m.layer.diag)[2], 0], Float32[0, bias(m.layer.diag)[2], 0]
            mutate_outputs(m, inds)
            @test typeof(layer(m)) <: LayerNorm
            assertlayer(m.layer.diag, Wexp, bexp)
        end

        array(arr) = arr
        array(arr::TrackedArray) = arr.data
        setpar(::Any, x) = x
        setpar(::TrackedArray, x) = param(x)

        function assertnorm(l, meanexp, varexp)
            @test l.β.data == meanexp
            @test l.γ.data == varexp
            @test l.μ == meanexp
            @test l.σ² == varexp
        end

        @testset "$l MutableLayer" for l in (BatchNorm, InstanceNorm, n -> GroupNorm(n,n))
            m = MutableLayer(l(5))
            l_orig = layer(m)

            @test nin(m) == nin(m.layer) == nout(m) == nout(m.layer) == 5
            @test minΔninfactor(m) == 1
            @test minΔnoutfactor(m) == 1
            m.layer = mapchildren(par -> setpar(par, collect(Float32, 1:5)), layer(m))

            inds = [1,3,4]
            expall = Float32.(inds)
            mutate_inputs(m, inds)
            assertnorm(m.layer, inds, inds)

            mutate_outputs(m, [-1, 2, -1])
            assertnorm(m.layer, [0, 3, 0], [1, 3, 1])
        end

        @testset "GroupNorm MutableLayer with groups" begin

            #Test with groups of 2
            m = MutableLayer(GroupNorm(6,3))
            m.layer = mapchildren(par -> setpar(par, collect(Float32, 1:length(par))), layer(m))
            mutate_inputs(m, [1,2,5,6])
            @test layer(m).μ == [1, 3]
            @test layer(m).σ² == [1, 3]

            # Now when dimensions don't add up: size 8 becomes size 9
            m = MutableLayer(GroupNorm(8,4))
            mutate_inputs(m, [1,3,-1,-1,4,-1,7,-1,8])
            # Current alg for selecting which group to pick in this case is poor, don't wanna test it :)
            @test length(layer(m).μ) == 3
            @test length(layer(m).σ²) == 3

        end
    end

    @testset "Recurrent layers" begin

        function assertrecurrent(l, Wiexp, Whexp, bexp, hexp, sexp)
            assertlayer(l, Wiexp, bexp)
            @test hiddenweights(l) == Whexp
            @test hiddenstate(l) == hexp
            @test state(l) == sexp
        end

        function setparsrnn(l)
            bias(l)[1:end] = collect(Float32, 1:nout(l)*outscale(l))
            hiddenstate(l)[1:end] = collect(Float32, 1:nout(l))
            state(l)[1:end] = collect(Float32, 1:nout(l))
        end

        function setparslstm(l)
            bias(l)[1:end] = collect(Float32, 1:nout(l)*outscale(l))
            foreach(h -> h[1:end] = collect(Float32, 1:nout(l)), hiddenstate(l))
            foreach(h -> h[1:end] = collect(Float32, 1:nout(l)), state(l))
        end

        @testset "RNN MutableLayer" begin
            m = MutableLayer(RNN(3, 4))
            setparsrnn(layer(m))

            @test nin(m) == nin(m.layer) == 3
            @test nout(m) == nout(m.layer) == 4
            @test minΔninfactor(m) == 1
            @test minΔnoutfactor(m) == 1

            inds = [1, 3]
            Wiexp = weights(layer(m))[:, inds]
            Whexp = copy(hiddenweights(layer(m)))
            bexp = copy(bias(layer(m)))
            hexp = copy(hiddenstate(layer(m)))
            sexp = copy(state(layer(m)))
            mutate_inputs(m, inds)
            assertrecurrent(layer(m), Wiexp, Whexp, bexp, hexp, sexp)

            inds = [1,-1, 2]
            Wiexp = permutedims(hcat(weights(layer(m))[1, :], zeros(Float32, 2), weights(layer(m))[2, :]))
            wh = hiddenweights(layer(m))
            Whexp = [wh[1, 1] 0 wh[1, 2]; zeros(Float32, 1, 3); wh[2, 1] 0 wh[2, 2]]
            bexp = Float32[bias(layer(m))[1], 0, bias(layer(m))[2]]
            hexp = Float32[hiddenstate(layer(m))[1], 0, hiddenstate(layer(m))[2]]
            sexp = Float32[state(layer(m))[1], 0, state(layer(m))[2]]
            mutate_outputs(m, inds)
            assertrecurrent(layer(m), Wiexp, Whexp, bexp, hexp, sexp)

            #Sanity check that the layer still seems to work after mutation
            output = m(reshape(collect(Float32, 1:2*10), 2,10))
            @test size(output) == (3, 10)
            @test isnan.(output) == falses(size(output))
        end

        @testset "LSTM MutableLayer" begin
            m = MutableLayer(LSTM(3, 4))
            setparslstm(layer(m))

            @test nin(m) == nin(m.layer) == 3
            @test nout(m) == nout(m.layer) == 4
            @test minΔninfactor(m) == 1
            @test minΔnoutfactor(m) == 1

            inds = [1, 3]
            Wiexp = weights(layer(m))[:, inds]
            Whexp = copy(hiddenweights(layer(m)))
            bexp = copy(bias(layer(m)))
            hexp = copy(hiddenstate(layer(m)))
            sexp = copy(state(layer(m)))
            mutate_inputs(m, inds)
            assertrecurrent(layer(m), Wiexp, Whexp, bexp, hexp, sexp)

            inds = [1,-1, 2]
            wi = weights(layer(m))
            scalerange = (0:outscale(layer(m))-1) .* nout(layer(m))
            Wiexp = permutedims(mapfoldl(offs -> hcat(wi[1+offs, :], zeros(Float32, 2), wi[2+offs, :]), hcat, scalerange))
            wh = hiddenweights(layer(m))
            Whexp = mapfoldl(offs -> [wh[1+offs, 1] 0 wh[1+offs, 2]; zeros(Float32, 1, 3); wh[2+offs, 1] 0 wh[2+offs, 2]], vcat, scalerange)
            bexp = mapfoldl(offs -> Float32[bias(layer(m))[1+offs], 0, bias(layer(m))[2+offs]], vcat, scalerange)
            hexp = map(hs -> Float32[hs[1], 0, hs[2]], hiddenstate(layer(m)))
            sexp = map(hs -> Float32[hs[1], 0, hs[2]], state(layer(m)))
            mutate_outputs(m, inds)
            assertrecurrent(layer(m), Wiexp, Whexp, bexp, hexp, sexp)

            #Sanity check that the layer still seems to work after mutation
            output = m(reshape(collect(Float32, 1:2*10), 2,10))
            @test size(output) == (3, 10)
            @test isnan.(output) == falses(size(output))
        end

        @testset "GRU MutableLayer" begin
            m = MutableLayer(GRU(3, 4))
            setparsrnn(layer(m))

            @test nin(m) == nin(m.layer) == 3
            @test nout(m) == nout(m.layer) == 4
            @test minΔninfactor(m) == 1
            @test minΔnoutfactor(m) == 1

            inds = [1, 3]
            Wiexp = weights(layer(m))[:, inds]
            Whexp = copy(hiddenweights(layer(m)))
            bexp = copy(bias(layer(m)))
            hexp = copy(hiddenstate(layer(m)))
            sexp = copy(state(layer(m)))
            mutate_inputs(m, inds)
            assertrecurrent(layer(m), Wiexp, Whexp, bexp, hexp, sexp)

            inds = [1,-1, 2]
            wi = weights(layer(m))
            scalerange = (0:outscale(layer(m))-1) .* nout(layer(m))
            Wiexp = permutedims(mapfoldl(offs -> hcat(wi[1+offs, :], zeros(Float32, 2), wi[2+offs, :]), hcat, scalerange))
            wh = hiddenweights(layer(m))
            Whexp = mapfoldl(offs -> [wh[1+offs, 1] 0 wh[1+offs, 2]; zeros(Float32, 1, 3); wh[2+offs, 1] 0 wh[2+offs, 2]], vcat, scalerange)
            bexp = mapfoldl(offs -> Float32[bias(layer(m))[1+offs], 0, bias(layer(m))[2+offs]], vcat, scalerange)
            hexp = Float32[hiddenstate(layer(m))[1], 0, hiddenstate(layer(m))[2]]
            sexp = Float32[state(layer(m))[1], 0, state(layer(m))[2]]
            mutate_outputs(m, inds)
            assertrecurrent(layer(m), Wiexp, Whexp, bexp, hexp, sexp)

            #Sanity check that the layer still seems to work after mutation
            output = m(reshape(collect(Float32, 1:2*10), 2,10))
            @test size(output) == (3, 10)
            @test isnan.(output) == falses(size(output))
        end
    end

    @testset "Clone MutableLayer" begin
            m = MutableLayer(Dense(2,3))
            cloned = clone(m)
            @test layer(cloned) !== layer(m)
            @test cloned([1, 2]) == m([1, 2])
    end

    @testset "LazyMutable" begin
        @testset "LazyMutable Dense factory" begin

            struct DenseFactory end
            function NaiveNASflux.dispatch!(m::LazyMutable, ::DenseFactory, x)
                m.mutable = MutableLayer(Dense(nin(m), nout(m)))
                return m(x)
            end
            m = LazyMutable(DenseFactory(), 2, 3)

            @test typeof(m.mutable) == DenseFactory
            expected = m(Float32[2,3])

            @test typeof(m.mutable) == MutableLayer
            @test m(Float32[2,3]) == expected

            #Now mutate before create
            m = LazyMutable(DenseFactory(), 2, 3)

            mutate_inputs(m, 5)
            mutate_outputs(m, 4)

            #No eagerness allowed :)
            @test m.mutable != MutableLayer
            expected = m(Float32[0,1,2,3,4])

            @test typeof(m.mutable) == MutableLayer
            @test m(Float32[0,1,2,3,4]) == expected
        end

        @testset "LazyMutable with Dense MutableLayer" begin
            m = MutableLayer(Dense(3,4))
            mlazy = LazyMutable(m)

            Wexp = weights(layer(m))
            bexp = bias(layer(m))

            @test minΔninfactor(m) == 1
            @test minΔnoutfactor(m) == 1

            mutate_inputs(mlazy, [1, 3])
            assertlayer(layer(m), Wexp, bexp)

            mutate_outputs(mlazy, [2, 4, -1])
            assertlayer(layer(m), Wexp, bexp)

            @test layer(m) == layer(mlazy)
            @test layertype(m) == layertype(mlazy)

            expected = mlazy(Float32[2,3])

            @test nin(mlazy) == nin(m) == 2
            @test nout(mlazy) == nout(m) == 3

            @test expected == m(Float32[2,3])
        end

        @testset "LazyeMutable DepthwiseConv" begin
            m = LazyMutable(MutableLayer(DepthwiseConv((2,2),(3=>6*3))))

            @test nin(m) == nin(layer(m)) == 3
            @test nout(m) == nout(layer(m)) == 18

            input = reshape(collect(Float32, 1:3*3*3), 3, 3, 3, 1)
            @test m(input) == layer(m)(input)

            ins = [1, 3]
            wouts = [1, 2, 5, 6]
            outs = mapfoldl(i -> 2 * i .+ [-1, -0] ,vcat, wouts)
            Wexp, bexp = weights(layer(m))[:,:,wouts,ins], bias(layer(m))[outs]

            mutate_inputs(m, ins)
            mutate_outputs(m, outs)
            @test size(m(ones(Float32, 3,3,2,2)))[3:4] == (8, 2)

            assertlayer(layer(m), Wexp, bexp)
        end

        @testset "LazyMutable reselect" begin
            m = LazyMutable(MutableLayer(Dense(5,5)))

            mutate_inputs(m, [-1, 1, 3, -1, 4])
            @test m.inputs == [-1, 1, 3, -1, 4]

            mutate_inputs(m, [1,2,4,5])
            @test m.inputs == [-1, 1,-1, 4]

            mutate_outputs(m, [2, -1, 3, -1, 4])
            @test m.outputs == [2, -1, 3, -1, 4]

            mutate_outputs(m, [1, 2, 5,-1])
            @test m.outputs == [2, -1, 4, -1]

            @test m(Float32[1,3,5,7]) == layer(m)(Float32[1,3,5,7])
        end

        @testset "Clone" begin
            mlazy = LazyMutable(MutableLayer(Dense(2,3)))
            cloned = clone(mlazy)
            @test layer(cloned) !== layer(mlazy)
            @test cloned([1, 2]) == mlazy([1, 2])
        end

        @testset "Treelike" begin
            m = LazyMutable(MutableLayer(Dense(2,3)))
            visitfun(x) = x
            visitdense = false
            function visitfun(l::TrackedArray)
                visitdense = true
                return l
            end

            mutate_inputs(m, [-1, -1, 1, 2])

            Flux.mapleaves(visitfun, m)
            @test visitdense
        end
    end
end


text:
module packageDevTest1

using ForwardDiff
include("helper.jl")

end


text:
# ------------------- #
# BasisMatrix Type #
# ------------------- #

abstract AbstractBasisMatrixRep
typealias ABSR AbstractBasisMatrixRep

immutable Tensor <: ABSR end
immutable Direct <: ABSR end
immutable Expanded <: ABSR end

type BasisMatrix{BST<:ABSR, TM<:AbstractMatrix}
    order::Matrix{Int}
    vals::Matrix{TM}
end

Base.show{BST}(io::IO, b::BasisMatrix{BST}) =
    print(io, "BasisMatrix{$BST} of order $(b.order)")

Base.ndims(bs::BasisMatrix) = size(bs.order, 2)
# TODO: determine if I want `size(bs, i::Int)` methods to do
# `prod([size(b.vals[1,j], i) for j=1:ndims(b)])`

# not the same if either type parameter is different
function =={BST1<:ABSR,BST2<:ABSR}(::BasisMatrix{BST1}, ::BasisMatrix{BST2})
    false
end

function =={BST<:ABSR,TM1<:AbstractMatrix,TM2<:AbstractMatrix}(::BasisMatrix{BST,TM1},
                                                               ::BasisMatrix{BST,TM2})
    false
end

# if type parameters are the same, then it is the same if all fields are the
# same
function =={BST<:ABSR,TM<:AbstractMatrix}(b1::BasisMatrix{BST,TM},
                                          b2::BasisMatrix{BST,TM})
    b1.order == b2.order && b1.vals == b2.vals
end

# -------------- #
# Internal Tools #
# -------------- #

@inline function _checkx(N, x::AbstractMatrix)
    size(x, 2) != N && error("Basis is $N dimensional, x must have $N columns")
    x
end

@inline function _checkx{T}(N, x::AbstractVector{T})
    # if we have a 1d basis, we can evaluate at each point
    if N == 1
        return x
    end

    # If Basis is > 1d, one evaluation point and reshape to (1,N) if possible...
    if length(x) == N
        return reshape(x, 1, N)
    end

    # ... or throw an error
    error("Basis is $N dimensional, x must have $N elements")
end

@inline function _checkx{T}(N, x::Vector{Vector{T}})
    # for BasisMatrix{Tensor} family. Need one vector per dimension
    if length(x) == N
        return x
    end

    # otherwise throw an error
    error("Basis is $N dimensional, need one Vector per dimension")
end

# common checks for all convert methods
function check_convert(bs::BasisMatrix, order::Matrix)
    d = ndims(bs)
    numbas, d1 = size(order)

    d1 != d && error("ORDER incompatible with basis functions")  # 35-37

    # 39-41
    if any(minimum(order, 1) .< bs.order)
        error("Order of derivative operators exceeds basis")
    end
    return d, numbas, d1
end

"""
Do common transformations to all constructor of `BasisMatrix`

##### Arguments

- `N::Int`: The number of dimensions in the corresponding `Basis`
- `x::AbstractArray`: The points for which the `BasisMatrix` should be
constructed
- `order::Array{Int}`: The order of evaluation for each dimension of the basis

##### Returns

- `m::Int`: the total number of derivative order basis functions to compute.
This will be the number of rows in the matrix form of `order`
- `order::Matrix{Int}`: A `m × N` matrix that, for each of the `m` desired
specifications, gives the derivative order along all `N` dimensions
- `minorder::Matrix{Int}`: A `1 × N` matrix specifying the minimum desired
derivative order along each dimension
- `numbases::Matrix{Int}`: A `1 × N` matrix specifying the total number of
distinct derivative orders along each dimension
- `x::AbstractArray`: The properly transformed points at which to evaluate
the basis

"""
function check_basis_structure(N::Int, x, order)
    order = _check_order(N, order)

    # initialize basis structure (66-74)
    m = size(order, 1)  # by this time order is a matrix
    if m > 1
        minorder = minimum(order, 1)
        numbases = (maximum(order, 1) - minorder) + 1
    else
        minorder = order + zeros(Int, 1, N)
        numbases = fill(1, 1, N)
    end

    x = _checkx(N, x)

    return m, order, minorder, numbases, x
end

# give the type of the `vals` field based on the family type parameter of the
# corresponding basis. `Spline` and `Lin` use sparse, `Cheb` uses dense
# a hybrid must fall back to a generic AbstractMatrix{Float64}
_vals_type(::Type{Spline}) = SparseMatrixCSC{Float64,Int}
_vals_type(::Type{Lin}) = SparseMatrixCSC{Float64,Int}
_vals_type(::Type{Cheb}) = Matrix{Float64}
_vals_type(::Type{BasisFamily}) = AbstractMatrix{Float64}

# conveneince method so we can pass an instance of the type also
_vals_type{TF<:BasisFamily}(::TF) = _vals_type(TF)

# --------------- #
# convert methods #
# --------------- #

# no-op. Don't worry about the order argument.
Base.convert{T<:ABSR}(::Type{T}, bs::BasisMatrix{T}, order=bs.order) = bs

# funbconv from direct to expanded
function Base.convert{TM}(::Type{Expanded}, bs::BasisMatrix{Direct,TM},
                      order=fill(0, 1, size(bs.order, 2)))
    d, numbas, d1 = check_convert(bs, order)

    vals = Array(TM, numbas, 1)

    for i=1:numbas
        vals[i] = bs.vals[order[i, d] - bs.order[d]+1, d]  # 63
        for j=d-1:-1:1
            vals[i] = row_kron(vals[i],
                               bs.vals[order[i, j] - bs.order[j]+1, j])  #65
        end
    end
    BasisMatrix{Expanded,TM}(order, vals)
end

# funbconv from tensor to expanded
function Base.convert{TM}(::Type{Expanded}, bs::BasisMatrix{Tensor,TM},
                      order=fill(0, 1, size(bs.order, 2)))
    d, numbas, d1 = check_convert(bs, order)

    vals = Array(TM, numbas, 1)

    for i=1:numbas  # 54
        vals[i] = bs.vals[order[i, d] - bs.order[d]+1, d]  # 55
        for j=d-1:-1:1  # 56
            # 57
            vals[i] = kron(vals[i], bs.vals[order[i, j] - bs.order[j]+1, j])
        end
    end

    BasisMatrix{Expanded,TM}(order, vals)
end

# funbconv from tensor to direct
# TODO: there is probably a more efficient way to do this, but since I don't
#       plan on doing it much, this will do for now. The basic point is that
#       we need to expand the rows of each element of `vals` so that all of
#       them have prod([size(v, 1) for v in bs.vals])) rows.
function Base.convert{TM}(::Type{Direct}, bs::BasisMatrix{Tensor,TM},
                      order=fill(0, 1, size(bs.order, 2)))
    d, numbas, d1 = check_convert(bs, order)
    vals = Array(TM, numbas, d)
    raw_ind = Array(Vector{Int}, d)

    for j=1:d
        for i=1:size(bs.vals, 1)
            if !(isempty(bs.vals[i, j]))  # 84
                raw_ind[j] = collect(1:size(bs.vals[i, j], 1))  # 85
                break
            end
        end
    end

    ind = gridmake(raw_ind...)  # 90

    for j=1:d, i=1:numbas
        if !(isempty(bs.vals[i, j]))  # 93
            vals[i, j] = bs.vals[i, j][ind[:, j], :]  # 94
        end
    end

    BasisMatrix{Direct,TM}(order, vals)
end

# ------------ #
# Constructors #
# ------------ #

# quick function to take order+vals and return expanded form for 1d problems
function to_expanded(out_order::Matrix{Int}, vals::Array)
    vals = vals[collect(out_order + (1 - minimum(out_order)))]
    vals = reshape(vals, size(out_order))
    BasisMatrix{Expanded,eltype(vals)}(out_order, vals)
end

# method to construct BasisMatrix in direct or expanded form based on
# a matrix of `x` values  -- funbasex
function BasisMatrix{N,BF}(basis::Basis{N,BF}, ::Direct,
                           x::AbstractArray=nodes(basis)[1], order=0)

    m, order, minorder, numbases, x = check_basis_structure(N, x, order)
    # 76-77
    out_order = minorder
    out_format = Direct()
    val_type = _vals_type(BF)
    vals = Array(val_type, maximum(numbases), N)

    # now do direct form, will convert to expanded later if needed
    for j=1:N
        # 126-130
        if (m > 1)
            orderj = unique(order[:, j])
        else
            orderj = order[1, j]
        end

        #131-135
        if length(orderj) == 1
            vals[1, j] = evalbase(basis.params[j], x[:, j], orderj[1])
        else
            vals[orderj-minorder[j]+1, j] =
                evalbase(basis.params[j], x[:, j], orderj)
        end
    end

    # if N == 1, switch to expanded format and return it directly
    # 140-145
    if N == 1  # 1 dimension
        return to_expanded(order, vals)
    end

    # construct Direct Format
    BasisMatrix{Direct,val_type}(out_order, vals)
end

function BasisMatrix(basis::Basis, ::Expanded,
                     x::AbstractArray=nodes(basis)[1], order=0)  # funbasex
    # create direct form, then convert to expanded
    bsd = BasisMatrix(basis, Direct(), x, order)
    convert(Expanded, bsd, bsd.order)
end

function BasisMatrix{N,BF,TV<:AbstractVector}(basis::Basis{N,BF}, ::Tensor,
    x::AbstractVector{TV}=nodes(basis)[2], order=0)

    m, order, minorder, numbases, x = check_basis_structure(N, x, order)
    out_order = minorder
    out_format = Tensor()
    val_type = _vals_type(BF)
    vals = Array(val_type, maximum(numbases), N)

    # construct tensor base
    for j=1:N
        # 113-117
        if (m > 1)
            orderj = unique(order[:, j])
        else
            orderj = order[1, j]
        end

        #118-122
        if length(orderj) == 1
            vals[1, j] = evalbase(basis.params[j], x[j], orderj)
        else
            vals[orderj-minorder[j]+1, j] = evalbase(basis[:params][j], x[j],
                                                     orderj)
        end
    end

    N == 1 ? to_expanded(order, vals) :
             BasisMatrix{Tensor,val_type}(out_order, vals)
end