text:
# copied & modified from brutus

"""
    code_ircode_by_signature([pass, ]sig; world=get_world_counter(), interp=NativeInterpreter(world))

Get `IRCode` by given signature, one can use the first argument to transform the `IRCode` during interpretation.
"""
function code_ircode_by_signature(pass, @nospecialize(sig); world=get_world_counter(), interp=NativeInterpreter(world))
    return map(Base._methods_by_ftype(sig, -1, world)) do data
        mi = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance}, (Any, Any, Any), data[3], data[1], data[2])
        code_ircode_by_mi(pass, mi; world, interp)
    end
end

function code_ircode_by_signature(@nospecialize(sig); world=get_world_counter(), interp=NativeInterpreter(world))
    return code_ircode_by_signature(default_julia_pass, sig; world, interp)
end

"""
    code_ircode([pass, ]f, types; world=get_world_counter(), interp=NativeInterpreter(world))

Get `IRCode` by given function `f` and its argument types `types`. An option argument `pass`
can be specified as a transform function on IRCode during type inference.
"""
function code_ircode(@nospecialize(f), @nospecialize(types); world=get_world_counter(), interp=NativeInterpreter(world))
    return code_ircode(default_julia_pass, f, types; world, interp)
end

function code_ircode(@nospecialize(pass), @nospecialize(f), @nospecialize(types); world=get_world_counter(), interp=NativeInterpreter(world))
    return [code_ircode_by_mi(pass, mi; world, interp) for mi in method_instances(f, types, world)]
end

"""
    code_ircode_by_mi(f, mi::MethodInstance; world=get_world_counter(), interp=NativeInterpreter(world))

Return the `IRCode` object along with inferred return type.

# Arguments

- `f(ir::IRCode, sv::OptimizationState) -> IRCode`: optimization passes to run.
- `mi::MethodInstance`: method instance.

# Kwargs

- `world::Int`: world number, default is calling `Core.Compiler.get_world_counter`.
- `interp::AbstractInterpreter`: the interpreter to use for inference.
"""
function code_ircode_by_mi(f, mi::MethodInstance; world=get_world_counter(), interp=NativeInterpreter(world))
    return typeinf_lock() do
        result = Core.Compiler.InferenceResult(mi)
        @static if VERSION < v"1.8-"
            frame = Core.Compiler.InferenceState(result, false, interp)
        else
            frame = Core.Compiler.InferenceState(result, :local, interp)
        end
        frame === nothing && return nothing

        if typeinf(interp, frame)
            opt_params = OptimizationParams(interp)
            sv = OptimizationState(frame, opt_params, interp)
            preserve_coverage = coverage_enabled(sv.mod)
            @static if VERSION < v"1.8-"
                ci = sv.src; nargs = sv.nargs - 1;
                ir = convert_to_ircode(ci, copy_exprargs(ci.code), preserve_coverage, nargs, sv)
                ir = slot2reg(ir, ci, nargs, sv)
            else
                ci = sv.src
                ir = convert_to_ircode(ci, sv)
                ir = slot2reg(ir, ci, sv)
            end
            # passes
            ir = f(ir, sv)
            sv.src.inferred = true
        end

        frame.inferred || return nothing
        @static if VERSION < v"1.8-"
            resize!(ir.argtypes, sv.nargs)
        else
            svdef = sv.linfo.def
            nargs = isa(svdef, Method) ? Int(svdef.nargs) : 0
            resize!(ir.argtypes, nargs)
        end
        return ir => widenconst(result.result)
    end
end

"""
    code_ircode_by_mi(mi::MethodInstance; world=get_world_counter(), interp=NativeInterpreter(world))

The default `code_ircode_by_mi` that uses the default Julia compiler optimization passes.
See also [`code_ircode_by_mi`](@ref).
"""
function code_ircode_by_mi(mi::MethodInstance; world=get_world_counter(), interp=NativeInterpreter(world))
    return code_ircode_by_mi(default_julia_pass, mi; world, interp)
end


text:
@testset "Triangle.jl Interface" begin
    @testset "Basic Triangulation" begin

        @testset "basic_triangulation_vertices" begin
            points = Array{Float64,2}([0. 0.; 1. 0.; 0. 1.])
            @test Triangle.basic_triangulation_vertices(points)[1] == points
        end

        @testset "basic_triangulation with point map" begin
            points = Array{Float64,2}([0. 0.; 1. 0.; 0. 1.])
            points_map = [1, 2, 3]
            @test Triangle.basic_triangulation(points,points_map)[1] == points_map
        end

        @testset "basic_triangulation_vertices with point map" begin
            points = Array{Float64,2}([0. 0.; 1. 0.; 0. 1.])
            points_map = [1, 2, 3]
            @test Triangle.basic_triangulation_vertices(points,points_map)[1] == points
        end

        # Triforce
        @testset "basic_triangulation with unordered point map" begin
            points = Array{Float64,2}([0. 0.; 4. 0.; 2. 3.; 8. 0.; 6. 3.; 4. 6.])
            points_map = [1, 2, 4, 3, 5, 6]
            result_tri = [[1, 2, 4],[4, 2, 5],[5, 6, 4],[3, 5, 2]]
            @test Triangle.basic_triangulation(points,points_map) == result_tri
        end

        # Tetris L (no boundary)
        @testset "basic_triangulation on L tetris piece" begin
            points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
            points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
            triangles = Triangle.basic_triangulation(points,points_map)
            @test length(triangles) == 5
        end

        @testset "basic_triangulation_vertices with point map on L tetris piece" begin
            points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
            points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
            triangles = Triangle.basic_triangulation_vertices(points,points_map)
            @test length(triangles) == 5
        end

         # Tetris L (all boundary)
        @testset "constrained_triangulation" begin
            points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
            points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
            edges_list = Array{Int64,2}([1 2; 2 3; 3 4; 4 5; 5 6; 6 1])
            triangles = Triangle.constrained_triangulation(points,points_map,edges_list)
            @test length(triangles) == 4
        end

        @testset "constrained_triangulation_vertices" begin
            points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
            points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
            edges_list = Array{Int64,2}([1 2; 2 3; 3 4; 4 5; 5 6; 6 1])
            triangles = Triangle.constrained_triangulation_vertices(points,points_map,edges_list)
            @test length(triangles) == 4
        end

        # Tetris L (custom boundary)
        @testset "constrained_triangulation with boundary" begin
            points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
            points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
            edges_list = Array{Int64,2}([1 2; 2 3; 3 4; 4 5; 5 6; 6 1])
            edge_boundary = [false, false, true, true, false, false]
            triangles = Triangle.constrained_triangulation(points,points_map,edges_list,edge_boundary)
            @test length(triangles) == 4
        end

        @testset "constrained_triangulation_vertices with boundary" begin
            points = [0. 0.; 0. 3.; 1. 3.; 1. 1.; 2. 1.; 2. 0.]
            points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
            edges_list = Array{Int64,2}([1 2; 2 3; 3 4; 4 5; 5 6; 6 1])
            edge_boundary = [false, false, true, true, false, false]
            triangles = Triangle.constrained_triangulation_vertices(points,points_map,edges_list,edge_boundary)
            @test length(triangles) == 4
        end

        # Triforce (without center triangle)
        @testset "constrained_triangulation with boundary and holes" begin
            points = [0. 0.; 4. 0.; 2. 3.; 8. 0.; 6. 3.; 4. 6.]
            points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
            edges_list = Array{Int64,2}([1 2; 2 3; 3 1; 2 4; 4 5; 5 2; 3 5; 5 6; 6 3])
            edge_boundary = [false,true,false,false,false,true,true,false,false]
            holes_list = [4. 2.]
            triangles = Triangle.constrained_triangulation(points,points_map,edges_list,edge_boundary,holes_list)
            @test length(triangles) == 3
        end

        @testset "constrained_triangulation_vertices with boundary and holes" begin
            points = [0. 0.; 4. 0.; 2. 3.; 8. 0.; 6. 3.; 4. 6.]
            points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
            edges_list = Array{Int64,2}([1 2; 2 3; 3 1; 2 4; 4 5; 5 2; 3 5; 5 6; 6 3])
            edge_boundary = [false,true,false,false,false,true,true,false,false]
            holes_list = [4. 2.]
            triangles = Triangle.constrained_triangulation_vertices(points,points_map,edges_list,edge_boundary,holes_list)
            @test length(triangles) == 3
        end
    end
end

ext: jl


text:
import CommonRLInterface

"""
Represents an abstract environment for reinforcement learning agents. Has several functions that need to be implemented to work.
All interfaces expect an abstract environment!
"""
abstract type AbstractEnvironment end

"""
    start!(env::AbstractEnvironment, args...)

Function to start the passed environment `env`. There are three variants. Two which start the environment from a random start state (as implemented with reset!) and another which starts the environment from a provided start state. These three variants call the `reset!` functions of the same call signiture.

Returns the starting state of the environment.
"""
function start!(env::AbstractEnvironment, args...)
    reset!(env, args...)
    return get_state(env)
end

function start!(env::CommonRLInterface.AbstractEnv)
    reset!(env)
    return CommonRLInterface.observe(env)
end

"""
    step!(env::AbstractEnvironment, action, args...)

Update the state of the passed environment `env` based on the underlying dynamics and the action.
"""
function step!(env::AbstractEnvironment, action, args...) # -> env, state, reward, terminal
    environment_step!(env, action, args...)
    return get_state(env), get_reward(env), is_terminal(env)
end

function step!(env::CommonRLInterface.AbstractEnv, action) # -> env, state, reward, terminal
    rew = CommonRLInterface.act!(env, action)
    return CommonRLInterface.observe(env), rew, CommonRLInterface.terminated(env)
end

#---------------------------#
#   Interface functions
#---------------------------#
"""
    reset!(env::AbstractEnvironment, args...)

Reset the environment to initial conditions based on the random number generator.

You can implement with or without a personally defined RNG. If you choose to not implement with a personally maintained RNG remember this is not a thread safe function.
"""
function reset!(env::AbstractEnvironment, args...)
    # reset!(env::AbstractEnvironment, Random.GLOBAL_RNG; kwargs...)
    @error "Implement reset! for environment $(typeof(env))"
end

"""
    environment_step!(env::AbstractEnvironment, action, args...)

Update the state of the environment based on the underlying dynamics and the action. This is not used directly, but through the step function.

You can implement with or without a personally defined RNG. If you choose to not implement with a personally maintained RNG remember this is not a thread safe function.
"""
function environment_step!(env::AbstractEnvironment, action, args...)
    @error "Implement environment_step for environment $(typeof(env))"
end

"""
    get_reward(env::AbstractEnvironment)

Retrieve reward for the current state of the environment.
"""
function get_reward(env::AbstractEnvironment) # -> determines if the agent_state is terminal
    @error "Implement get_reward for environment $(typeof(env))"
end

"""
    is_terminal(env::AbstractEnvironment)

Check if the environment is in a terminal state
"""
function is_terminal(env::AbstractEnvironment) # -> determines if the agent_state is terminal
    @error "Implement is_terminal for environment $(typeof(env))"
end

"""
    get_state(env::AbstractEnvironment)

Retrieve the current state of the environment
"""
function get_state(env::AbstractEnvironment) # -> determines if the agent_state is terminal
    @error "Implement get_state for environment $(typeof(env))"
end

"""
    get_actions(env::AbstractEnvironment)

Returns the set of actions available to take.
"""
# get_actions
function get_actions(env::AbstractEnvironment)
    return nothing
end


text:
using PyCall
using PyPlot
using FEniCS

function fenics_plot_scalar(f)
    py"""
    from fenics import *
    def FF(f):
        import matplotlib.pyplot as plt
        plt.figure()
        c = plot(f, colobar=True)
        plt.colorbar(c)
        plt.show()
    """
    py"FF"(f.pyobject)
end

function fenics_plot_vec2d(f)
    py"""
    from fenics import *
    def FF(f):
        import matplotlib.pyplot as plt
        fx = dot(f,Constant((1.0, 0.0)))
        fy = dot(f,Constant((0.0, 1.0)))
        plt.figure()
        plt.subplot(1,2,1)
        c = plot(fx, colobar=True)
        plt.colorbar(c)
        plt.subplot(1,2,2)
        c = plot(fy, colobar=True)
        plt.colorbar(c)
        plt.tight_layout()
        plt.show()
    """
    py"FF"(f.pyobject)
end

# Calculate strain
function fenics_strain(u)
    0.5 * (nabla_grad(u) + FEniCS.Transpose(nabla_grad(u)))
end

# Calculate
function fenics_stress(u, lambda, mu)
     lambda * nabla_div(u) * Identity(2) + 2 * mu * fenics_strain(u)
end

function fenics_solve(V, f)
    bc = DirichletBC(V, Constant((0, 0)), "on_boundary")
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(fenics_stress(u, lambda, mu), fenics_strain(v)) * dx
    L = FEniCS.dot(f, v) * dx
    u = FeFunction(V)
    lvsolve(a, L, u, bc)
    u
end

function fenics_interior_eval(f, x, y)
    py"""
    def G(f, x, y):
        return f(x, y)
    """
    py"G"(f.pyobject, x, y)
end

function fenics_tensor_components(V, f)
    fxx_fxy = project(FEniCS.dot(f, Constant((1,0))), V)
    fxy_fyy = project(FEniCS.dot(f, Constant((0,1))), V)
    fxx_fxy, fxy_fyy
end

function fenics_eval_u(fenics_u, idx, x, y)
    n = length(idx)
    u = zeros(2*n)
    for idx_i in 1:n
        i = idx[idx_i]
        u[(2*idx_i - 1):(2*idx_i)] = fenics_interior_eval(fenics_u, x[i], y[i])
    end
    u
end

function fenics_eval_s(sxx_sxy, sxy_syy, idx, x, y)
    n = length(idx)
    s = zeros(3*n)
    for idx_i in 1:n
        i = idx[idx_i]
        sxx, sxy = fenics_interior_eval(sxx_sxy, x[i], y[i])
        sxy2, syy = fenics_interior_eval(sxy_syy, x[i], y[i])
        s[3*idx_i - 2] = sxx
        s[3*idx_i - 1] = syy
        s[3*idx_i - 0] = sxy
    end
    s
end

function stress_to_trac(S, idx, xn, yn)
    n = length(idx)
    t = zeros(2*n)
    for idx_i in 1:n
        i = idx[idx_i]
        sxx, syy, sxy = S[3*idx_i-2:3*idx_i]
        tx = sxx * xn[i] + sxy * yn[i]
        ty = sxy * xn[i] + syy * yn[i]
        t[2*idx_i - 1] = tx
        t[2*idx_i] = ty
    end
    t
end

function twopanel(xgrid, ygrid, npts, U, S, idx, els; figsize=(12,3), ylim=[-20000, 2000])
    # Start of nice visualization
    figure(figsize=figsize)
    subplot(1, 2, 1)
    field = sqrt.(U[:, 1].^2 + U[:, 2].^2)
    ncontours = 10
    lowfield = 0
    highfield = 100
    ncontours = LinRange(lowfield, highfield, 11)

    xlim = [-20000 20000]
    scale = 1.0
    fieldmax = maximum(@.abs(field))
    contourf(reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
             reshape(field, npts, npts), ncontours,
             vmin=lowfield, vmax=highfield,
             cmap = get_cmap("magma"))
    clim(lowfield, highfield)
    colorbar(fraction=0.020, pad=0.05, extend = "both", label = L"$||u||$ (m)")
    contour(reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
            reshape(field, npts, npts), ncontours,
            vmin=lowfield, vmax=highfield,
            linewidths=0.25, colors="w")
    plotelements(els)
    gca().set_aspect("equal")
    gca().set_xlim([xlim[1], xlim[2]])
    gca().set_ylim([ylim[1], ylim[2]])
    gca().set_xticks([-20000, -10000, 0, 10000, 20000])
    gca().set_yticks([-20000, -10000, 0])
    xlabel(L"$x$ (m)")
    ylabel(L"$y$ (m)")
    xv = [els.x1[idx["T"]] ; els.x2[idx["T"]][end] ; 20000 ; -20000]
    yv = [els.y1[idx["T"]] ; els.y2[idx["T"]][end]; 2000 ; 2000]
    fill(xv, yv, "lightgray", zorder=30)

    # Find principle stress orientations
    stressdiff = zeros(length(U[:, 1]))
    for i in 1:length(U[:,1])
        mat = [S[i, 1] S[i, 3]; S[i, 3] S[i, 2]]
        ev = eigvals(mat)
        stressdiff[i] = log10(abs(ev[2]-ev[1]))
    end

    subplot(1, 2, 2)
    field = stressdiff
    xlim = [-20000 20000]
    lowfield = 5
    highfield = 10
    ncontours = LinRange(lowfield,highfield, 11)
    contourf(reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
             reshape(field, npts, npts), ncontours,
             vmin=lowfield, vmax=highfield,
             cmap=get_cmap("viridis"))
    clim(lowfield, highfield)
    colorbar(fraction=0.020, pad=0.05, extend="both",
             label=L"$\Delta \sigma$ (Pa)")
    contour(reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
            reshape(field, npts, npts), ncontours,
            vmin=lowfield, vmax=highfield,
            linewidths=0.25, colors="w")
    plotelements(els)
    gca().set_aspect("equal")
    gca().set_xlim([xlim[1], xlim[2]])
    gca().set_ylim([ylim[1], ylim[2]])
    gca().set_xticks([-20000, -10000, 0, 10000, 20000])
    gca().set_yticks([-20000, -10000, 0])
    xlabel(L"$x$ (m)")
    ylabel(L"$y$ (m)")
    xv = [els.x1[idx["T"]] ; els.x2[idx["T"]][end] ; 20000 ; -20000]
    yv = [els.y1[idx["T"]] ; els.y2[idx["T"]][end]; 2000 ; 2000]
    fill(xv, yv, "lightgray", zorder=30)
    tight_layout()
    return nothing
end


text:
@testset "CategoricalTree is a collection" begin
    xs = collect(1:8)
    ct = CategoricalTree(xs)

    @test isempty(CategoricalTree([]))
    @test !isempty(ct)

    @test length(ct) == length(xs)

end

@testset "CategoricalTree is an indexable collection" begin
    xs = collect(1:8)
    ct = CategoricalTree(xs)

    @test firstindex(ct) == firstindex(xs)
    @test lastindex(ct) == lastindex(xs)
end


text:
abstract type AbstractLinregSolver end

"""
    SolveQR()

Pass as `method` to [`linregress`](@ref) to use QR factorization to solve the
linear regression system. Slower but more accurate for ill-conditioned systems
than [`SolveCholesky`](@ref).
"""
struct SolveQR <: AbstractLinregSolver end

"""
    SolveCholesky()

Pass as `method` to [`linregress`](@ref) to use Cholesky factorization to solve
the linear regression system. Faster but less accurate for ill-conditioned
systems than [`SolveQR`](@ref).
"""
struct SolveCholesky <: AbstractLinregSolver end

struct LinearRegressor{T<:Union{Matrix{<:Real},Vector{<:Real}}}
    intercept::Bool
    coeffs::T
end

"""
    coef(lr::LinearRegressor)

Returns the linear regression coefficients of a [`linregress`](@ref) result. If
`lr.intercept` is `true`, this includes the bias (intercept) in the last
position.
"""
coef(lr::LinearRegressor) = lr.coeffs
islinear(lr::LinearRegressor) = true

slope(lr::LinearRegressor) = lr.intercept ? _slope(lr.coeffs) : lr.coeffs
bias(lr::LinearRegressor) = lr.intercept ? _bias(lr.coeffs) : _zero_bias(lr.coeffs)

_slope(coef::AbstractMatrix) = @view coef[1:end-1, :]
_slope(coef::AbstractVector) = @view coef[1:end-1]
_bias(coef::AbstractMatrix) = @view coef[end:end, :]
_bias(coef::AbstractVector) = coef[end]
_zero_bias(coef::AbstractMatrix) = zero(eltype(lr.coeffs), 1, size(lr.coeffs, 2))
_zero_bias(coef::AbstractVector) = zero(eltype(lr.coeffs))

function (lr::LinearRegressor)(X::AbstractMatrix)
    if lr.intercept
        return X * _slope(lr.coeffs) .+ _bias(lr.coeffs)
    else
        return X * lr.coeffs
    end
end

function (lr::LinearRegressor{<:Vector})(x::AbstractVector)
    if lr.intercept
        return x'_slope(lr.coeffs) + _bias(lr.coeffs)
    else
        return dot(x, lr.coeffs)
    end
end

if VERSION < v"1.5"
    _vec(x) = collect(vec(x))
else
    _vec(x) = vec(x)
end

function (lr::LinearRegressor{<:Matrix})(x::AbstractVector)
    if lr.intercept
        return _vec(x'_slope(lr.coeffs) + _bias(lr.coeffs))
    else
        return _vec(x'lr.coeffs)
    end
end

@doc raw"""
    linregress(X, y; intercept=true, method=SolveQR())
    linregress(X, y, weights; intercept=true, method=SolveQR())

Do (possibly weighted) linear regression to obtain coefficients β such that `ŷ
= X * β` minimizes `‖ŷ - y‖²`. In the default case, this corresponds to solving
`X \ y`.

Returns a `LinearRegressor`, which can be passed to [`coef`](@ref) to extract
the coefficients β or called with a vector or matrix `X` to predict at a single
point or set of points.

Currently assumes that, if `X` is a Matrix, that `size(X) == (length(y),
num_coefs)` (i.e., each row describes the features for one observation).

## Keyword Arguments:

- If `intercept` is `true`, implicitly adds a column of ones to `X` to model a
  bias term.

- `method` determines how to solve the linear regression system (see
  [`SolveQR`](@ref) and [`SolveCholesky`](@ref)).
"""
function linregress(X, y, weights=nothing; intercept=true, method=SolveQR())
    size(X, 1) == size(y, 1) || throw(DimensionMismatch("size of X and y does not match in first dimension"))
    if intercept
        X = _append_bias_column(method, X)
    else
        X = _maybe_convert(method, X)
    end
    coeffs = if weights === nothing
        _lin_solve(method, X, y)
    else
        length(weights) == size(y, 1) || throw(DimensionMismatch("size of y and weights does not match in first dimension"))
        W = Diagonal(weights)
        _lin_solve(method, X, y, W)
    end
    return LinearRegressor(intercept, coeffs)
end

function _lin_solve(solver::AbstractLinregSolver, X, y, W)
    # √W X \ √W y
    Wsqrt = sqrt(W)
    return _lin_solve(solver, Wsqrt * X, Wsqrt * y)
end

function _lin_solve(::SolveQR, X, y)
    return X \ y
end

function _lin_solve(::SolveCholesky, X, y, W)
    # X' W X \ (X' W y)
    return ldiv!(cholesky!(Hermitian(X' * W * X)), X' * (W * y))
end

function _lin_solve(::SolveCholesky, X, y)
    # X'X \ X'y
    return ldiv!(cholesky!(Hermitian(X'X)), X'y)
end

function _append_bias_column(::AbstractLinregSolver, X)
    # avoids any conversions at this point
    ones_column = ones(eltype(X), size(X, 1))
    return [X ones_column]
end

function _append_bias_column(::SolveCholesky, X)
    # already promotes to the correct type
    ones_column = ones(LinearAlgebra.choltype(X), size(X, 1))
    return [X ones_column]
end

_maybe_convert(::AbstractLinregSolver, X) = X
_maybe_convert(::SolveCholesky, X) = convert(Array{LinearAlgebra.choltype(X)}, X)


text:
struct Scaler{T} <: Number
    s::T
end

function Scaler(x::Vector{T}; robust::Bool=false) where T
    s = robust ? 1.482602median(abs.(x .- median(x))) : std(x)
    Scaler(s)
end

(s::Scaler)(x::Vector) = x ./ s.s

mutable struct StepScaler{S,F} <: AbstractStep
    selections::S
    params::Union{Nothing, Dict{Symbol,Scaler}}
    robust::Bool
    trained::Bool
    skip::Bool
    prehook::F
end

function transform!(s::StepScaler, df)
    for (k,f) in s.params
        df[k] = f(df[k])
    end
    df
end
function transform(s::StepScaler, df)
    transform!(s, copy(df))
end
function fit!(s::StepScaler, df)
    s.params = Dict(col => Scaler(s.prehook(df[col])) for col in getselectionkeys(df, s.selections))
    s.trained = true
end
function step_scale!(r::Recipe, s...; robust::Bool=false, abs::Bool=false, skip::Bool=false, prehook=identity)
    push!(r.steps,
          StepScaler([s...],
                     nothing,
                     robust,
                     false,
                     skip,
                     prehook))
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# trigonometric functions, remap to keep consistent API with Base
@_remap broadcasted(::typeof(sin),  x::NDArray) sin(x)
@_remap broadcasted(::typeof(cos),  x::NDArray) cos(x)
@_remap broadcasted(::typeof(tan),  x::NDArray) tan(x)
@_remap broadcasted(::typeof(asin), x::NDArray) arcsin(x)
@_remap broadcasted(::typeof(acos), x::NDArray) arccos(x)
@_remap broadcasted(::typeof(atan), x::NDArray) arctan(x)

# hyperbolic functions, remap to keep consistent API with Base
@_remap broadcasted(::typeof(sinh),  x::NDArray) sinh(x)
@_remap broadcasted(::typeof(cosh),  x::NDArray) cosh(x)
@_remap broadcasted(::typeof(tanh),  x::NDArray) tanh(x)
@_remap broadcasted(::typeof(asinh), x::NDArray) arcsinh(x)
@_remap broadcasted(::typeof(acosh), x::NDArray) arccosh(x)
@_remap broadcasted(::typeof(atanh), x::NDArray) arctanh(x)


text:
"""
    DataFrames.getmaxwidths(df::AbstractDataFrame,
                            io::IO,
                            rowindices1::AbstractVector{Int},
                            rowindices2::AbstractVector{Int},
                            rowlabel::Symbol,
                            rowid::Union{Integer, Nothing},
                            show_eltype::Bool,
                            buffer::IOBuffer)

Calculate, for each column of an AbstractDataFrame, the maximum
string width used to render the name of that column, its type, and the
longest entry in that column -- among the rows of the data frame
will be rendered to IO. The widths for all columns are returned as a
vector.

Return a `Vector{Int}` giving the maximum string widths required to render
each column, including that column's name and type.

NOTE: The last entry of the result vector is the string width of the
implicit row ID column contained in every `AbstractDataFrame`.

# Arguments
- `df::AbstractDataFrame`: The data frame whose columns will be printed.
- `io::IO`: The `IO` to which `df` is to be printed
- `rowindices1::AbstractVector{Int}: A set of indices of the first
  chunk of the AbstractDataFrame that would be rendered to IO.
- `rowindices2::AbstractVector{Int}: A set of indices of the second
  chunk of the AbstractDataFrame that would be rendered to IO. Can
  be empty if the AbstractDataFrame would be printed without any
  ellipses.
- `rowlabel::AbstractString`: The label that will be used when rendered the
  numeric ID's of each row. Typically, this will be set to "Row".
- `rowid`: Used to handle showing `DataFrameRow`.
- `show_eltype`: Whether to print the column type
   under the column name in the heading.
- `buffer`: buffer passed around to avoid reallocations in `ourstrwidth`
"""
function getmaxwidths(df::AbstractDataFrame,
                      io::IO,
                      rowindices1::AbstractVector{Int},
                      rowindices2::AbstractVector{Int},
                      rowlabel::Symbol,
                      rowid::Union{Integer, Nothing},
                      show_eltype::Bool,
                      buffer::IOBuffer,
                      truncstring::Int)
    maxwidths = Vector{Int}(undef, size(df, 2) + 1)

    undefstrwidth = ourstrwidth(io, "#undef", buffer, truncstring)

    ct = show_eltype ? batch_compacttype(Any[eltype(c) for c in eachcol(df)]) : String[]
    j = 1
    for (col_idx, (name, col)) in enumerate(pairs(eachcol(df)))
        # (1) Consider length of column name
        # do not truncate column name
        maxwidth = ourstrwidth(io, name, buffer, 0)

        # (2) Consider length of longest entry in that column
        for indices in (rowindices1, rowindices2), i in indices
            if isassigned(col, i)
                maxwidth = max(maxwidth, ourstrwidth(io, col[i], buffer, truncstring))
            else
                maxwidth = max(maxwidth, undefstrwidth)
            end
        end
        if show_eltype
            # do not truncate eltype name
            maxwidths[j] = max(maxwidth, ourstrwidth(io, ct[col_idx], buffer, 0))
        else
            maxwidths[j] = maxwidth
        end
        j += 1
    end

    # do not truncate rowlabel
    if rowid isa Nothing
        rowmaxwidth1 = isempty(rowindices1) ? 0 : ndigits(maximum(rowindices1))
        rowmaxwidth2 = isempty(rowindices2) ? 0 : ndigits(maximum(rowindices2))
        maxwidths[j] = max(max(rowmaxwidth1, rowmaxwidth2),
                           ourstrwidth(io, rowlabel, buffer, 0))
    else
        maxwidths[j] = max(ndigits(rowid), ourstrwidth(io, rowlabel, buffer, 0))
    end

    return maxwidths
end

"""
    show(io::IO, mime::MIME, df::AbstractDataFrame)

Render a data frame to an I/O stream in MIME type `mime`.

# Arguments
- `io::IO`: The I/O stream to which `df` will be printed.
- `mime::MIME`: supported MIME types are: `"text/plain"`, `"text/html"`, `"text/latex"`,
  `"text/csv"`, `"text/tab-separated-values"` (the last two MIME types do not support
   showing `#undef` values)
- `df::AbstractDataFrame`: The data frame to print.

Additionally selected MIME types support passing the following keyword arguments:
- MIME type `"text/plain"` accepts all listed keyword arguments and therir behavior
  is identical as for `show(::IO, ::AbstractDataFrame)`
- MIME type `"text/html"` accepts `summary` keyword argument which
  allows to choose whether to print a brief string summary of the data frame.

# Examples
```jldoctest
julia> show(stdout, MIME("text/latex"), DataFrame(A = 1:3, B = ["x", "y", "z"]))
\\begin{tabular}{r|cc}
\t& A & B\\\\
\t\\hline
\t& Int64 & String\\\\
\t\\hline
\t1 & 1 & x \\\\
\t2 & 2 & y \\\\
\t3 & 3 & z \\\\
\\end{tabular}
14

julia> show(stdout, MIME("text/csv"), DataFrame(A = 1:3, B = ["x", "y", "z"]))
"A","B"
1,"x"
2,"y"
3,"z"
```
"""
Base.show(io::IO, mime::MIME, df::AbstractDataFrame)
Base.show(io::IO, mime::MIME"text/html", df::AbstractDataFrame;
          summary::Bool=true, eltypes::Bool=true) =
    _show(io, mime, df, summary=summary, eltypes=eltypes)
Base.show(io::IO, mime::MIME"text/latex", df::AbstractDataFrame; eltypes::Bool=true) =
    _show(io, mime, df, eltypes=eltypes)
Base.show(io::IO, mime::MIME"text/csv", df::AbstractDataFrame) =
    printtable(io, df, header = true, separator = ',')
Base.show(io::IO, mime::MIME"text/tab-separated-values", df::AbstractDataFrame) =
    printtable(io, df, header = true, separator = '\t')
Base.show(io::IO, mime::MIME"text/plain", df::AbstractDataFrame; kwargs...) =
    show(io, df; kwargs...)

##############################################################################
#
# HTML output
#
##############################################################################

function digitsep(value::Integer)
    # Adapted from https://github.com/IainNZ/Humanize.jl
    value = string(abs(value))
    group_ends = reverse(collect(length(value):-3:1))
    groups = [value[max(end_index - 2, 1):end_index]
              for end_index in group_ends]
    return join(groups, ',')
end

function html_escape(cell::AbstractString)
    cell = replace(cell, "&"=>"&")
    cell = replace(cell, "<"=>"<")
    cell = replace(cell, ">"=>">")
    # Replace quotes so that the resulting string could also be used in the attributes of
    # HTML tags
    cell = replace(cell, "\""=>""")
    cell = replace(cell, "'"=>"'")
    return cell
end

function _show(io::IO, ::MIME"text/html", df::AbstractDataFrame;
               summary::Bool=true, eltypes::Bool=true, rowid::Union{Int, Nothing}=nothing)
    _check_consistency(df)

    # we will pass around this buffer to avoid its reallocation in ourstrwidth
    buffer = IOBuffer(Vector{UInt8}(undef, 80), read=true, write=true)

    if rowid !== nothing
        if size(df, 2) == 0
            rowid = nothing
        elseif size(df, 1) != 1
            throw(ArgumentError("rowid may be passed only with a single row data frame"))
        end
    end

    mxrow, mxcol = size(df)
    if get(io, :limit, false)
        tty_rows, tty_cols = displaysize(io)
        mxrow = min(mxrow, tty_rows)
        maxwidths = getmaxwidths(df, io, 1:mxrow, 0:-1, :X, nothing, true, buffer, 0) .+ 2
        mxcol = min(mxcol, searchsortedfirst(cumsum(maxwidths), tty_cols))
    end

    cnames = _names(df)[1:mxcol]
    write(io, "
")
    if summary
        write(io, "
$(digitsep(nrow(df))) rows × $(digitsep(ncol(df))) columns")
        if mxcol < size(df, 2)
            write(io, " (omitted printing of $(size(df, 2)-mxcol) columns)")
        end
        write(io, "
")
    end
    write(io, "
")
    write(io, "
")
    write(io, "
")
    write(io, "	")
    for column_name in cnames
        write(io, "$(html_escape(String(column_name)))	")
    end
    write(io, "
")
    if eltypes
        write(io, "
")
        write(io, "	")
        # We put a longer string for the type into the title argument of the  element,
        # which the users can hover over. The limit of 256 characters is arbitrary, but
        # we want some maximum limit, since the types can sometimes get really-really long.
        types = Any[eltype(df[!, idx]) for idx in 1:mxcol]
        ct, ct_title = batch_compacttype(types), batch_compacttype(types, 256)
        for j in 1:mxcol
            s = html_escape(ct[j])
            title = html_escape(ct_title[j])
            write(io, "$s	")
        end
        write(io, "
")
    end
    write(io, "
")
    write(io, "
")
    for row in 1:mxrow
        write(io, "
")
        if rowid === nothing
            write(io, "$row	")
        else
            write(io, "$rowid	")
        end
        for column_name in cnames
            if isassigned(df[!, column_name], row)
                cell_val = df[row, column_name]
                if ismissing(cell_val)
                    write(io, "missing	")
                elseif cell_val isa Markdown.MD
                    write(io, "")
                    show(io, "text/html", cell_val)
                    write(io, "	")
                elseif cell_val isa SHOW_TABULAR_TYPES
                    write(io, "")
                    cell = sprint(ourshow, cell_val, 0)
                    write(io, html_escape(cell))
                    write(io, "	")
                else
                    cell = sprint(ourshow, cell_val, 0)
                    write(io, "$(html_escape(cell))	")
                end
            else
                write(io, "#undef	")
            end
        end
        write(io, "
")
    end
    if size(df, 1) > mxrow
        write(io, "
")
        write(io, "⋮	")
        for column_name in cnames
            write(io, "⋮	")
        end
        write(io, "
")
    end
    write(io, "
")
    write(io, "
")
    write(io, "
")
end

function Base.show(io::IO, mime::MIME"text/html", dfr::DataFrameRow;
                   summary::Bool=true, eltypes::Bool=true)
    r, c = parentindices(dfr)
    summary && write(io, "
DataFrameRow ($(length(dfr)) columns)
")
    _show(io, mime, view(parent(dfr), [r], c), summary=false, eltypes=eltypes, rowid=r)
end

function Base.show(io::IO, mime::MIME"text/html", dfrs::DataFrameRows;
                   summary::Bool=true, eltypes::Bool=true)
    df = parent(dfrs)
    summary && write(io, "
$(nrow(df))×$(ncol(df)) DataFrameRows
")
    _show(io, mime, df, summary=false, eltypes=eltypes)
end

function Base.show(io::IO, mime::MIME"text/html", dfcs::DataFrameColumns;
                   summary::Bool=true, eltypes::Bool=true)
    df = parent(dfcs)
    if summary
        write(io, "
$(nrow(df))×$(ncol(df)) DataFrameColumns
")
    end
    _show(io, mime, df, summary=false, eltypes=eltypes)
end

function Base.show(io::IO, mime::MIME"text/html", gd::GroupedDataFrame)
    N = length(gd)
    keys = html_escape(join(string.(groupcols(gd)), ", "))
    keystr = length(gd.cols) > 1 ? "keys" : "key"
    groupstr = N > 1 ? "groups" : "group"
    write(io, "
$(nameof(typeof(gd))) with $N $groupstr based on $keystr: $keys
")
    if N > 0
        nrows = size(gd[1], 1)
        rows = nrows > 1 ? "rows" : "row"

        identified_groups = [html_escape(string(col, " = ",
                                                repr(first(gd[1][!, col]))))
                             for col in gd.cols]

        write(io, "
First Group ($nrows $rows): ")
        join(io, identified_groups, ", ")
        write(io, "
")
        show(io, mime, gd[1], summary=false)
    end
    if N > 1
        nrows = size(gd[N], 1)
        rows = nrows > 1 ? "rows" : "row"

        identified_groups = [html_escape(string(col, " = ",
                                                repr(first(gd[N][!, col]))))
                             for col in gd.cols]

        write(io, "
⋮
")
        write(io, "
Last Group ($nrows $rows): ")
        join(io, identified_groups, ", ")
        write(io, "
")
        show(io, mime, gd[N], summary=false)
    end
end

##############################################################################
#
# LaTeX output
#
##############################################################################

function latex_char_escape(char::Char)
    if char == '\\'
        return "\\textbackslash{}"
    elseif char == '~'
        return "\\textasciitilde{}"
    else
        return string('\\', char)
    end
end

function latex_escape(cell::AbstractString)
    replace(cell, ['\\','~', '#', '$', '%', '&', '_', '^', '{', '}']=>latex_char_escape)
end

function _show(io::IO, ::MIME"text/latex", df::AbstractDataFrame;
               eltypes::Bool=true, rowid=nothing)
    _check_consistency(df)

    # we will pass around this buffer to avoid its reallocation in ourstrwidth
    buffer = IOBuffer(Vector{UInt8}(undef, 80), read=true, write=true)

    if rowid !== nothing
        if size(df, 2) == 0
            rowid = nothing
        elseif size(df, 1) != 1
            throw(ArgumentError("rowid may be passed only with a single row data frame"))
        end
    end

    mxrow, mxcol = size(df)
    if get(io, :limit, false)
        tty_rows, tty_cols = get(io, :displaysize, displaysize(io))
        mxrow = min(mxrow, tty_rows)
        maxwidths = getmaxwidths(df, io, 1:mxrow, 0:-1, :X, nothing, true, buffer, 0) .+ 2
        mxcol = min(mxcol, searchsortedfirst(cumsum(maxwidths), tty_cols))
    end

    cnames = _names(df)[1:mxcol]
    alignment = repeat("c", mxcol)
    write(io, "\\begin{tabular}{r|")
    write(io, alignment)
    mxcol < size(df, 2) && write(io, "c")
    write(io, "}\n")
    write(io, "\t& ")
    header = join(map(c -> latex_escape(string(c)), cnames), " & ")
    write(io, header)
    mxcol < size(df, 2) && write(io, " & ")
    write(io, "\\\\\n")
    write(io, "\t\\hline\n")
    if eltypes
        write(io, "\t& ")
        ct = batch_compacttype(Any[eltype(df[!, idx]) for idx in 1:mxcol])
        header = join(latex_escape.(ct), " & ")
        write(io, header)
        mxcol < size(df, 2) && write(io, " & ")
        write(io, "\\\\\n")
        write(io, "\t\\hline\n")
    end
    for row in 1:mxrow
        write(io, "\t")
        write(io, @sprintf("%d", rowid === nothing ? row : rowid))
        for col in 1:mxcol
            write(io, " & ")
            if !isassigned(df[!, col], row)
                print(io, "\\emph{\\#undef}")
            else
                cell = df[row, col]
                if ismissing(cell)
                    print(io, "\\emph{missing}")
                elseif cell isa Markdown.MD
                    print(io, strip(repr(MIME("text/latex"), cell)))
                elseif cell isa SHOW_TABULAR_TYPES
                    print(io, "\\emph{")
                    print(io, latex_escape(sprint(ourshow, cell, 0, context=io)))
                    print(io, "}")
                else
                    if showable(MIME("text/latex"), cell)
                        show(io, MIME("text/latex"), cell)
                    else
                        print(io, latex_escape(sprint(ourshow, cell, 0, context=io)))
                    end
                end
            end
        end
        mxcol < size(df, 2) && write(io, " & \$\\dots\$")
        write(io, " \\\\\n")
    end
    if size(df, 1) > mxrow
        write(io, "\t\$\\dots\$")
        for col in 1:mxcol
            write(io, " & \$\\dots\$")
        end
        mxcol < size(df, 2) && write(io, " & ")
        write(io, " \\\\\n")
    end
    write(io, "\\end{tabular}\n")
end

function Base.show(io::IO, mime::MIME"text/latex", dfr::DataFrameRow; eltypes::Bool=true)
    r, c = parentindices(dfr)
    _show(io, mime, view(parent(dfr), [r], c), eltypes=eltypes, rowid=r)
end

Base.show(io::IO, mime::MIME"text/latex", dfrs::DataFrameRows; eltypes::Bool=true) =
	_show(io, mime, parent(dfrs), eltypes=eltypes)
Base.show(io::IO, mime::MIME"text/latex", dfcs::DataFrameColumns; eltypes::Bool=true) =
	_show(io, mime, parent(dfcs), eltypes=eltypes)

function Base.show(io::IO, mime::MIME"text/latex", gd::GroupedDataFrame)
    N = length(gd)
    keys = join(latex_escape.(string.(groupcols(gd))), ", ")
    keystr = length(gd.cols) > 1 ? "keys" : "key"
    groupstr = N > 1 ? "groups" : "group"
    write(io, "$(nameof(typeof(gd))) with $N $groupstr based on $keystr: $keys\n\n")
    if N > 0
        nrows = size(gd[1], 1)
        rows = nrows > 1 ? "rows" : "row"

        identified_groups = [latex_escape(string(col, " = ",
                                                 repr(first(gd[1][!, col]))))
                             for col in gd.cols]

        write(io, "First Group ($nrows $rows): ")
        join(io, identified_groups, ", ")
        write(io, "\n\n")
        show(io, mime, gd[1])
    end
    if N > 1
        nrows = size(gd[N], 1)
        rows = nrows > 1 ? "rows" : "row"

        identified_groups = [latex_escape(string(col, " = ",
                                                 repr(first(gd[N][!, col]))))
                             for col in gd.cols]

        write(io, "\n\$\\dots\$\n\n")
        write(io, "Last Group ($nrows $rows): ")
        join(io, identified_groups, ", ")
        write(io, "\n\n")
        show(io, mime, gd[N])
    end
end

##############################################################################
#
# MIME: text/csv and text/tab-separated-values
#
##############################################################################

escapedprint(io::IO, x::SHOW_TABULAR_TYPES, escapes::AbstractString) =
    escapedprint(io, summary(x), escapes)
escapedprint(io::IO, x::Any, escapes::AbstractString) =
    escapedprint(io, sprint(print, x), escapes)
escapedprint(io::IO, x::AbstractString, escapes::AbstractString) =
    escape_string(io, x, escapes)

function printtable(io::IO,
                    df::AbstractDataFrame;
                    header::Bool = true,
                    separator::Char = ',',
                    quotemark::Char = '"',
                    missingstring::AbstractString = "missing",
                    nothingstring::AbstractString = "nothing")
    _check_consistency(df)
    n, p = size(df)
    etypes = eltype.(eachcol(df))
    if header
        cnames = _names(df)
        for j in 1:p
            print(io, quotemark)
            print(io, cnames[j])
            print(io, quotemark)
            if j < p
                print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    quotestr = string(quotemark)
    for i in 1:n
        for j in 1:p
            cell = df[i, j]
            if ismissing(cell)
                print(io, missingstring)
            elseif isnothing(cell)
                print(io, nothingstring)
            else
                if cell isa Markdown.MD
                    print(io, quotemark)
                    r = repr(cell)
                    escapedprint(io, chomp(r), quotestr)
                    print(io, quotemark)
                elseif !(etypes[j] <: Real)
                    print(io, quotemark)
                    escapedprint(io, cell, quotestr)
                    print(io, quotemark)
                else
                    print(io, cell)
                end
            end
            if j < p
                print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    nothing
end

function Base.show(io::IO, mime::MIME"text/csv", dfr::DataFrameRow)
    r, c = parentindices(dfr)
    show(io, mime, view(parent(dfr), [r], c))
end

function Base.show(io::IO, mime::MIME"text/tab-separated-values", dfr::DataFrameRow)
    r, c = parentindices(dfr)
    show(io, mime, view(parent(dfr), [r], c))
end

Base.show(io::IO, mime::MIME"text/csv",
          dfs::Union{DataFrameRows, DataFrameColumns}) =
    show(io, mime, parent(dfs))
Base.show(io::IO, mime::MIME"text/tab-separated-values",
          dfs::Union{DataFrameRows, DataFrameColumns}) =
    show(io, mime, parent(dfs))

function Base.show(io::IO, mime::MIME"text/csv", gd::GroupedDataFrame)
    isfirst = true
    for sdf in gd
        printtable(io, sdf, header = isfirst, separator = ',')
        isfirst && (isfirst = false)
    end
end

function Base.show(io::IO, mime::MIME"text/tab-separated-values", gd::GroupedDataFrame)
    isfirst = true
    for sdf in gd
        printtable(io, sdf, header = isfirst, separator = '\t')
        isfirst && (isfirst = false)
    end
end


text:
export SimpleSparseBasis, SparseBasis, MaxBasis

# prototypes for space transforms and cutoffs
function transform end
function transform_d end
function transform_dd end
function inv_transform end

# ------------------------------------------------------------
#  Abstract polynomial degree business

"""
interface functions for `OneParticleBasis`
"""
function add_into_A! end

"""
interface functions for `OneParticleBasis`
"""
function add_into_A_dA! end

"""
`function degree` : compute some notion of degree of
the `arg` argument.
"""
function degree end

"""
`function scaling(b, p)`:

a scaling factor for a basis functions ϕ, which gives a rought estimate on
   the magnitude of ∇ᵖϕ e.g.,
```
ϕ = r^n Ylm
```
has scaling factor `n^p + l^p`, though sharper estimates are also possible.
"""
function scaling end

using LinearAlgebra: Diagonal

diagonal_regulariser(basis; diff = 0) = Diagonal(scaling(basis, diff))

"""
every scalar basis must implement this
"""
function rand_radial end

@doc raw"""
```math
\begin{aligned}
a &= A(x)
b &= B(a) = B(A(x))
∂b/∂x &= ∂b/∂a ⋅ ∂a/∂x
\end{aligned}
```
In code this becomes
```julia
a = evaluate(A, x)
b = evaluate(B, a)
∂a∂x = evaluate_d(A, x)
∂b∂x = frule_evaluate(B, a, ∂a∂x)
```
"""
function frule_evaluate! end

@doc raw"""
```math
\begin{aligned}
a &= A(x)
b &= B(a) = B(A(x))
∂b/∂x &= ∂b/∂a ⋅ ∂a/∂x
\end{aligned}
```
In code this becomes
```julia
a = evaluate(A, x)
b = evaluate(B, a)
∂b∂a = evaluate_d(B, a)
∂b∂x = rrule_evaluate(∂b∂a, A, x)
```
"""
function rrule_evaluate! end


text:
using CreditApproval
using Test

@testset "CreditApproval.jl" begin
end


text:
# broadcasting

using Base.Broadcast: BroadcastStyle, Broadcasted

struct oneArrayStyle{N} <: AbstractGPUArrayStyle{N} end
oneArrayStyle(::Val{N}) where N = oneArrayStyle{N}()
oneArrayStyle{M}(::Val{N}) where {N,M} = oneArrayStyle{N}()

BroadcastStyle(::Type{<:oneArray{T,N}}) where {T,N} = oneArrayStyle{N}()

Base.similar(bc::Broadcasted{oneArrayStyle{N}}, ::Type{T}) where {N,T} =
    similar(oneArray{T}, axes(bc))

Base.similar(bc::Broadcasted{oneArrayStyle{N}}, ::Type{T}, dims...) where {N,T} =
    oneArray{T}(undef, dims...)


text:
# lets test long approxConv chains

using Test
using IncrementalInference
using Statistics
using TensorCast

##

@testset "zdim size of factor is correct" begin

##

fg = initfg()
addVariable!(fg, :x0, ContinuousEuclid{3})
addVariable!(fg, :x1, ContinuousEuclid{3})

lr3 =  LinearRelative(MvNormal([1;0;0.0],diagm(ones(3))))
addFactor!(fg, [:x0; :x1], lr3, graphinit=false)

ccw = IIF._getCCW(fg, :x0x1f1)

@test IIF._getZDim(ccw) == 3

res = calcFactorResidualTemporary(lr3,
                                  (ContinuousEuclid{3},ContinuousEuclid{3}),
                                  [0;0;0.5],
                                  (zeros(3), [0;0;1.0])  )
#

@test (sum(abs.(res)) - 0.5) < 1e-10

##

end

@testset "approxConv basic and chains" begin

##

fg = generateCanonicalFG_Kaess()

# from a prior to neighbor
pts_ = approxConv(fg, :x1f1, :x1)

# lazy
@cast pts[i,j] := pts_[j][i]

@test Statistics.mean(pts) |> abs  < 0.4
@test 0.5 < Statistics.std(pts) < 1.5

# set a value in graph to start things off
initManual!(fg, :x1, pts_)

# legacy case where relative to neighbor
pts_ = approxConv(fg, :x1x2f1, :x2)
@cast pts[i,j] := pts_[j][i]

@test Statistics.mean(pts) |> abs  < 0.7
@test 0.7 < Statistics.std(pts) < 2

# along a chain of variables
pts_ = approxConv(fg, :x1, :x3)
@cast pts[i,j] := pts_[j][i]

@test Statistics.mean(pts) |> abs  < 1.5
@test 1.3 < Statistics.std(pts) < 3.0

# from a prior down the chain of variables
pts_ = approxConv(fg, :x1f1, :l2)
@cast pts[i,j] := pts_[j][i]

@test Statistics.mean(pts) |> abs  < 1.5
@test 1.6 < Statistics.std(pts) < 4.0

##
end

@testset "test approxConvBelief with partial prior" begin
##

fg = initfg()

addVariable!(fg, :x0, ContinuousEuclid{2})
pp = PartialPrior(Normal(),(2,))
addFactor!(fg, [:x0], pp, graphinit=false)

approxConvBelief(fg, :x0f1, :x0)

@info "second test with more complicated manifolds in testSpecialEuclidean2Mani.jl"

##
end

#


text:
@auto_hash_equals struct Generator
    coords::LatLon #latitude and longitude
    volt::Vector{Real} #voltages (may be multiplie values,
                                # our source gives no breakdown) # kV
    tech::AbstractString #technology type
    cap::Real #nominal capacity # MW
    pfactor::Real #nominal power factor
    minload::Real #minimum load # MW
    scap::Real #summer nominal capacity # MW
    wcap::Real #winter nominal capacity # MW
    shut2loadtime::AbstractString #time from cold shutdown to full load
    status::AbstractString #"OP" - Operational; "SB" - Standby;
                           # "OA" and "OS" - Out of service
end

function Generator{T <: Real}(
        coords::Tuple{T,T},
        volt,
        tech,
        cap,
        pfactor,
        minload,
        scap,
        wcap,
        shut2loadtime,
        status
    )
    return Generator(
        LatLon(coords[1], coords[2]),
        volt,
        tech,
        cap,
        pfactor,
        minload,
        scap,
        wcap,
        shut2loadtime,
        status
    )
end

coords(g::Generator) = g.coords
voltage(g::Generator) = g.volt
technology(g::Generator) = g.tech
capacity(g::Generator) = g.cap
power_factor(g::Generator) = g.pfactor
minimum_load(g::Generator) = g.minload
summer_capacity(g::Generator) = g.scap
winter_capacity(g::Generator) = g.wcap
status(g::Generator) = g.status


text:
__precompile__()
module COFF

# This package implements the ObjFileBase interface
import ObjFileBase
import ObjFileBase: sectionsize, sectionoffset, readheader, debugsections, deref,
    endianness, strtab_lookup, readmeta, isrelocatable, sectionname, isundef,
    symbolvalue, handle, symname, symbolnum, intptr, deref

# Reexports from ObjFileBase
export sectionsize, sectionoffset, readheader, debugsections

using StructIO

import Base: show, ==, *

export Relocations

include("constants.jl")

########## COFF.jl - An implementation of the PE/COFF File format ###############

#
# Represents the actual PE/COFF file
#
abstract COFFObjectHandle <: ObjFileBase.ObjectHandle
type COFFHandle{T<:IO} <: COFFObjectHandle
    # The IO object. This field is speciallized on to avoid dispatch performance
    # hits, especially when operating on an IOBuffer, which is an important
    # usecase for in-memory files
    io::T
    # position(io) of the start of the file in the io stream.
    start::Int
    # position(io) of the COFF header. header == start iff the file is a COFF object
    # file (as opposed to a PE image file)
    header::Int
    # A uniqued strtab will be filled in on demand
    strtab
    COFFHandle(io,start,header) = new(io,start,header)
end
ObjFileBase.handle(h::COFFHandle) = h
COFFHandle{T<:IO}(io::T,start,header) = COFFHandle{T}(io,start,header)
__init__() = push!(ObjFileBase.ObjHandles, COFFHandle)

endianness(x::COFFHandle) = :NativeEndian

show(io::IO, x::COFFHandle) = print(io, "COFF Object Handle")

import Base: read, readuntil, readbytes, write, seek, seekstart, position, eof

eof(handle::COFFHandle) = eof(handle.io)
seek{T<:IO}(io::COFFHandle{T},pos::Integer) = seek(io.io,io.start+pos)
seekstart(io::COFFHandle) = seek(io.io,io.start)
position{T<:IO}(io::COFFHandle{T}) = position(io.io)-io.start

import StructIO: unpack

unpack{T,ioT<:IO}(h::COFFHandle{ioT},::Type{T}) = unpack(h.io,T,:NativeEndian)
pack{T,ioT<:IO}(h::COFFHandle{ioT},::Type{T}) = pack(h.io,T,:NativeEndian)

#
# COFF Header
#
@struct immutable COFFHeader
    Machine::UInt16
    NumberOfSections::UInt16
    TimeDateStamp::UInt32
    PointerToSymbolTable::UInt32
    NumberOfSymbols::UInt32
    SizeOfOptionalHeader::UInt16
    Characteristics::UInt16
end

function printfield(io::IO,string,fieldlength)
    print(io," "^max(fieldlength-length(string),0))
    print(io,string)
end
printentry(io::IO,header,values...) = (printfield(io,header,21);println(io," ",values...))

isrelocatable(h::COFFHeader) = (h.Characteristics & (IMAGE_FILE_DLL | IMAGE_FILE_EXECUTABLE_IMAGE)) == 0
isrelocatable(h::COFFHandle) = isrelocatable(readheader(h))

using Base.Dates

function show(io::IO,h::COFFHeader)
    printentry(io,"Machine",IMAGE_FILE_MACHINE[h.Machine])
    printentry(io,"NumberOfSections", h.NumberOfSections)
    printentry(io,"TimeDateStamp", Dates.DateTime(1970,1,1,0,0,0,0) + Dates.Second(h.TimeDateStamp))
    printentry(io,"PointerToSymbolTable", "0x",hex(h.PointerToSymbolTable))
    printentry(io,"NumberOfSymbols", h.NumberOfSymbols)
    printentry(io,"SizeOfOptionalHeader", h.SizeOfOptionalHeader)
    Characteristics = String[]
    for (k,v) in IMAGE_FILE_CHARACTERISTICS
        ((k&h.Characteristics) != 0) && push!(Characteristics, v)
    end
    printentry(io,"Characteristics",join(Characteristics,", "))
end

#
# Optional Header
#
@struct immutable OptionalHeaderStandard
    Magic::UInt16
    MajorLinkerVersion::UInt8
    MinorLinkerVersion::UInt8
    SizeOfCode::UInt32
    SizeOfInitializedData::UInt32
    SizeOfUninitializedData::UInt32
    AddressOfEntryPoint::UInt32
    BaseOfCode::UInt32
end

@struct immutable IMAGE_DATA_DIRECTORY
    VirtualAddress::UInt32
    Size::UInt32
end

@struct immutable DataDirectories
    ExportTable::IMAGE_DATA_DIRECTORY
    ImportTable::IMAGE_DATA_DIRECTORY
    ResourceTable::IMAGE_DATA_DIRECTORY
    ExceptionTable::IMAGE_DATA_DIRECTORY
    CertificateTable::IMAGE_DATA_DIRECTORY
    BaseRelocatioNTable::IMAGE_DATA_DIRECTORY
    Debug::IMAGE_DATA_DIRECTORY
    Architecture::IMAGE_DATA_DIRECTORY
    GlobalPtr::IMAGE_DATA_DIRECTORY
    TLSTable::IMAGE_DATA_DIRECTORY
    LoadConfigTable::IMAGE_DATA_DIRECTORY
    BoundImport::IMAGE_DATA_DIRECTORY
    IAT::IMAGE_DATA_DIRECTORY
    DelayImportDescriptor::IMAGE_DATA_DIRECTORY
    CLRRuntimeHeader::IMAGE_DATA_DIRECTORY
    Reserverd::IMAGE_DATA_DIRECTORY
end

module PE32
    using StructIO
    import ..OptionalHeaderStandard, ..DataDirectories

    @struct immutable OptionalHeaderWindows
        ImageBase::UInt32
        SectionAlignment::UInt32
        FileAlignment::UInt32
        MajorOperatingSystemVersion::UInt16
        MinorOperatingSystemVersion::UInt16
        MajorImageVersion::UInt16
        MinorImageVersion::UInt16
        MajorSubsystemVersion::UInt16
        MinorSubsystemVersion::UInt16
        Win32VersionValue::UInt32
        SizeOfImage::UInt32
        SizeOfHeaders::UInt32
        CheckSum::UInt32
        Subsystem::UInt16
        DllCharacteristics::UInt16
        SizeOfStackReserve::UInt32
        SizeOfStackCommit::UInt32
        SizeOfHeapReserve::UInt32
        SizeOfHeapCommit::UInt32
        LoaderFlag::UInt32
        NumberOfRvaAndSizes::UInt32
    end

    @struct immutable OptionalHeader
        standard::OptionalHeaderStandard
        BaseOfData::UInt32
        windows::OptionalHeaderWindows
        directories::DataDirectories
    end

end

module PE32Plus
    using StructIO
    import ..OptionalHeaderStandard, ..DataDirectories

    @struct immutable OptionalHeaderWindows
        ImageBase::UInt64
        SectionAlignment::UInt32
        FileAlignment::UInt32
        MajorOperatingSystemVersion::UInt16
        MinorOperatingSystemVersion::UInt16
        MajorImageVersion::UInt16
        MinorImageVersion::UInt16
        MajorSubsystemVersion::UInt16
        MinorSubsystemVersion::UInt16
        Win32VersionValue::UInt32
        SizeOfImage::UInt32
        SizeOfHeaders::UInt32
        CheckSum::UInt32
        Subsystem::UInt16
        DllCharacteristics::UInt16
        SizeOfStackReserve::UInt64
        SizeOfStackCommit::UInt64
        SizeOfHeapReserve::UInt64
        SizeOfHeapCommit::UInt64
        LoaderFlag::UInt32
        NumberOfRvaAndSizes::UInt32
    end

    @struct immutable OptionalHeader
        standard::OptionalHeaderStandard
        windows::OptionalHeaderWindows
        directories::DataDirectories
    end
end

# Section Table

@struct immutable tiny_fixed_string
    str::UInt64
end

import Base: bytestring, show, print

function bytestring(x::tiny_fixed_string)
    a8 = reinterpret(UInt8,[x.str])
    z = findfirst(a8,0)
    String(a8[1:(z == 0 ? length(a8) : z-1)])
end
show(io::IO,x::tiny_fixed_string) = show(io,bytestring(x))
print(io::IO,x::tiny_fixed_string) = print(io,bytestring(x))

==(x::tiny_fixed_string,y::String) = bytestring(x) == y
==(x::String,y::tiny_fixed_string) = y==x

*(a::String,b::tiny_fixed_string) = a*bytestring(b)

@struct immutable SectionHeader <: ObjFileBase.Section{COFFObjectHandle}
    Name::tiny_fixed_string
    VirtualSize::UInt32
    VirtualAddress::UInt32
    SizeOfRawData::UInt32
    PointerToRawData::UInt32
    PointerToRelocations::UInt32
    PointerToLinenumbers::UInt32
    NumberOfRelocations::UInt16
    NumberOfLinenumbers::UInt16
    Characteristics::UInt32
end

function sectionname(header::SectionHeader; strtab = nothing, errstrtab=true)
    name = bytestring(header.Name)
    if name[1] == '/'
        if strtab != nothing
            return strtab_lookup(strtab,parse(Int, name[2:end]))
        elseif errstrtab
            error("Section name refers to the strtab, but no strtab given")
        end
    end
    return name
end

function show(io::IO, header::SectionHeader; strtab = nothing)
    name = bytestring(header.Name)
    name2 = sectionname(header; strtab = strtab, errstrtab=false)
    printentry(io,"Name",name,name!=name2?" => "*name2:"")
    printentry(io,"VirtualSize", "0x", hex(header.VirtualSize))
    printentry(io,"VirtualAddress", "0x", hex(header.VirtualAddress))
    printentry(io,"SizeOfRawData", "0x", hex(header.SizeOfRawData))
    printentry(io,"PointerToRawData", "0x", hex(header.PointerToRawData))
    printentry(io,"PointerToRelocations", "0x", hex(header.PointerToRelocations))
    printentry(io,"PointerToLinenumbers", "0x", hex(header.PointerToLinenumbers))
    printentry(io,"NumberOfRelocations", "0x", hex(header.NumberOfRelocations))
    printentry(io,"NumberOfLinenumbers", "0x", hex(header.NumberOfLinenumbers))
    Characteristics = String[]
    for (k,v) in IMAGE_SCN_CHARACTERISTICS
        if k & IMAGE_SCN_ALIGN_MASK != 0
            continue
        end
        ((k&header.Characteristics) != 0) && push!(Characteristics, v)
    end
    if header.Characteristics & IMAGE_SCN_ALIGN_MASK != 0
        push!(Characteristics,
            IMAGE_SCN_CHARACTERISTICS[header.Characteristics & IMAGE_SCN_ALIGN_MASK])
    end
    printentry(io,"Characteristics",join(Characteristics,", "))
end

@struct immutable SymbolName
    name::UInt64
end

function show(io::IO, sname::SymbolName; strtab = nothing, showredirect=true)
    if sname.name & typemax(UInt32) == 0
        if strtab !== nothing
            if showredirect
                print(io, sname.name >> 32, " => ")
            end
            print(io,strtab_lookup(strtab,sname.name>>32))
        else
            print(io, "/", sname.name >> 32)
        end
    else
        print(io,bytestring(tiny_fixed_string(sname.name)))
    end
end

function symname(sname::SymbolName;  kwargs...)
    buf = IOBuffer()
    show(buf,sname; kwargs..., showredirect=false)
    takebuf_string(buf)
end

@struct immutable SymtabEntry <: ObjFileBase.SymtabEntry{COFFHandle}
    Name::SymbolName
    Value::UInt32
    SectionNumber::UInt16
    Type::UInt16
    StorageClass::UInt8
    NumberOfAuxSymbols::UInt8
end align_packed

const IMAGE_SYM_UNDEFINED = 0
function isundef(entry::SymtabEntry)
    entry.StorageClass in (
        IMAGE_SYM_CLASS_EXTERNAL_DEF,
        IMAGE_SYM_CLASS_UNDEFINED_LABEL,
        IMAGE_SYM_CLASS_UNDEFINED_STATIC) ||
    (entry.StorageClass == IMAGE_SYM_CLASS_EXTERNAL &&
        entry.SectionNumber == IMAGE_SYM_UNDEFINED)
end
isfunction(entry::SymtabEntry) = entry.Type == 0x20

symname(sname::SymtabEntry; kwargs...) = symname(sname.Name; kwargs...)

function show(io::IO, entry::SymtabEntry; strtab = nothing)
    print(io, "0x", hex(entry.Value, 8), " ")
    if entry.SectionNumber == 0
        printfield(io, "*UND*", 5)
    elseif entry.SectionNumber == (-1%UInt16)
        printfield(io, "*ABS*", 5)
    elseif entry.SectionNumber == (-2%UInt16)
        printfield(io, "*DBG*", 5)
    else
        printfield(io, dec(entry.SectionNumber), 5)
    end
    print(io, " ",hex(entry.Type, 4)," ")
    #print(io, IMAGE_SYM_CLASS[entry.StorageClass]," ")
    show(io, entry.Name; strtab = strtab)
end

@struct immutable RelocationEntry
    VirtualAddress::UInt32
    SymbolTableIndex::UInt32
    Type::UInt16
end align_packed

function show(io::IO, entry::RelocationEntry; machine=IMAGE_FILE_MACHINE_UNKNOWN, syms = northing, strtab=nothing)
    print(io, "0x", hex(entry.VirtualAddress,8), " ")
    if machine == IMAGE_FILE_MACHINE_UNKNOWN
        print(io,hex(entry.Type,4)," ")
    else
        printfield(io,MachineRelocationMap[machine][entry.Type],maximum(map(length,MachineRelocationMap[machine])))
    end
    printfield(io,"@"*string(dec(entry.SymbolTableIndex)),6)
    if syms !== nothing
        print(io," -> ",symname(syms[entry.SymbolTableIndex+1]; strtab = strtab))
    end
end

# # # Higer level interface

import Base: length, getindex, start, done, next

# # Sections
immutable Sections <: ObjFileBase.Sections{COFFHandle}
    h::COFFHandle
    num::UInt16
    offset::Int
    Sections(h::COFFHandle, num::UInt16, offset::Int) = new(h,num,offset)
    function Sections(handle::COFFHandle,header::COFFHeader=readheader(handle))
        Sections(handle, header.NumberOfSections, handle.header + sizeof(COFFHeader) + header.SizeOfOptionalHeader)
    end
end
ObjFileBase.handle(s::Sections) = s.h
ObjFileBase.Sections(h::COFFHandle) = Sections(h)
ObjFileBase.mangle_sname(h::COFFHandle, name) = string(".", name)

immutable SectionRef <: ObjFileBase.SectionRef{COFFHandle}
    handle::COFFHandle
    no::Int
    offset::Int
    header::SectionHeader
end
@Base.pure ObjFileBase.SectionRef{T<:COFFHandle}(::Type{T}) = SectionRef

ObjFileBase.handle(s::SectionRef) = s.handle
sectionname(ref::SectionRef) = sectionname(ref.header; strtab=strtab(ref.handle))
deref(ref::SectionRef) = ref.header

function show(io::IO,x::SectionRef)
println(io,"0x",hex(x.offset,8),": Section #",x.no)
show(io, x.header; strtab=strtab(x.handle))
end

length(s::Sections) = s.num
const SectionHeaderSize = sizeof(SectionHeader)
function getindex(s::Sections,n)
    if n < 1 || n > length(s)
        throw(BoundsError())
    end
    offset = s.offset + (n-1)*SectionHeaderSize
    seek(s.h,offset)
    SectionRef(s.h,n,offset,unpack(s.h, SectionHeader))
end

start(s::Sections) = 1
done(s::Sections,n) = n > length(s)
next(s::Sections,n) = (s[n],n+1)

# # Symbols
immutable Symbols  <: ObjFileBase.Symbols{COFFHandle}
    h::COFFHandle
    num::UInt32
    offset::Int
    Symbols(h::COFFHandle, num, offset) = new(h,num,offset)
    function Symbols(handle::COFFHandle,header::COFFHeader=readheader(handle))
        Symbols(handle, header.NumberOfSymbols, header.PointerToSymbolTable)
    end
end
ObjFileBase.Symbols(h::COFFHandle) = Symbols(h)
ObjFileBase.handle(s::Symbols) = s.h

# Special case until Base is fixed
function Base.findnext(testf::Function, A::Symbols, start::Integer)
    while !done(A, start)
        i = start
        val, start = next(A, i)
        if testf(val)
            return i
        end
    end
    return 0
end

immutable SymbolRef <: ObjFileBase.SymbolRef{COFFHandle}
    handle::COFFHandle
    num::UInt32
    offset::Int
    entry::SymtabEntry
end
symbolnum(ref::SymbolRef) = ref.num
deref(ref::SymbolRef) = ref.entry
symname(sym::SymbolRef; strtab=COFF.strtab(sym.handle), kwargs...) = symname(sym.entry; strtab=strtab, kwargs...)
isfunction(entry::SymbolRef) = isfunction(deref(entry))
Base.seekstart(entry::SymbolRef) = seek(entry.handle,
    sectionoffset(Sections(entry.handle)[entry.entry.SectionNumber])+entry.entry.Value)

function symbolvalue(entry::Union{SymtabEntry, SymbolRef}, sects)
    entry = deref(entry)
    Value = entry.Value
    if entry.SectionNumber != 0 && entry.SectionNumber != (-1%UInt16) &&
            entry.SectionNumber != (-2%UInt16)
        sec = sects[entry.SectionNumber]
        Value += deref(sec).VirtualAddress
    end
    Value
end

function show(io::IO,x::SymbolRef)
    print(io,'[')
    printfield(io,dec(x.num),5)
    print(io,"] ")
    show(io,x.entry; strtab=strtab(x.handle))
end

length(s::Symbols) = endof(s) # This is incorrect, but required due to quirks in Base
endof(s::Symbols) = s.num
const SymtabEntrySize = sizeof(SymtabEntry)
function getindex(s::Symbols,n)
    if n < 1 || n > endof(s)
        throw(BoundsError())
    end
    offset = s.offset + (n-1)*SymtabEntrySize
    seek(s.h,offset)
    SymbolRef(s.h,n,offset,unpack(s.h, SymtabEntry))
end

start(s::Symbols) = 1
done(s::Symbols,n) = n > endof(s)
next(s::Symbols,n) = (x=s[n];(x,n+x.entry.NumberOfAuxSymbols+1))
Base.iteratorsize(::Type{Symbols}) = Base.SizeUnknown()

# String table

immutable StrTab
    h::COFFHandle
    size::Int
    offset::Int
end

function StrTab(h::COFFHandle, header=readheader(h))
    offset = header.PointerToSymbolTable+header.NumberOfSymbols*SymtabEntrySize
    seek(h, offset)
    return StrTab(h,read(h,UInt32),offset)
end

function strtab(h::COFFHandle)
    if isdefined(h, :strtab)
        return h.strtab
    end
    h.strtab = StrTab(h)
end
ObjFileBase.StrTab(h::COFFHandle) = strtab(h)
ObjFileBase.StrTab(h::Symbols) = ObjFileBase.StrTab(handle(h))

function strtab_lookup(strtab::StrTab, offset)
    seek(strtab.h,offset+strtab.offset)
    # Strip trailing \0
    readuntil(strtab.h,'\0')[1:end-1]
end
#

const PEMAGIC = reinterpret(UInt32,UInt8['P','E','\0','\0'])[1]
const MZ = reinterpret(UInt16,UInt8['M','Z'])[1]
function readmeta(io::IO, ::Type{COFFHandle})
    start = position(io)
    if read(io,UInt16) == MZ
        # Get the PE Header offset
        seek(io, start+0x3c)
        off = read(io, UInt32)
        # PE File
        seek(io, start+off)
        read(io, UInt32) == PEMAGIC || throw(ObjFileBase.MagicMismatch("Invalid PE magic"))
    else
        seek(io,start)
        read(io, UInt32) == PEMAGIC || throw(ObjFileBase.MagicMismatch("Invalid PE magic"))
        seek(io, start)
    end
    COFFHandle(io,start,position(io))
end

readheader(h::COFFHandle) = (seek(h.io,h.header); unpack(h, COFFHeader))
function readoptheader(h::COFFHandle)
    seek(h.io,h.header + sizeof(COFFHeader))
    standard = unpack(h, OptionalHeaderStandard)
    if standard.Magic == 0x10b # PE32
        BaseOfData = read(h, UInt32)
        windows = unpack(h, PE32.OptionalHeaderWindows)
        dirs = unpack(h, DataDirectories)
        return PE32.OptionalHeader(standard, BaseOfData, windows, dirs)
    elseif standard.Magic == 0x20b # PE32Plus
        windows = unpack(h, PE32Plus.OptionalHeaderWindows)
        dirs = unpack(h, DataDirectories)
        return PE32Plus.OptionalHeader(standard, windows, dirs)
    else
        error("Unknown magic")
    end
end

function intptr(h::COFFHandle)
    header = readheader(h)
    header.Machine == IMAGE_FILE_MACHINE_AMD64 ? UInt64 :
        header.Machine == IMAGE_FILE_MACHINE_I386 ? UInt32 :
            error("Unknown Machine Type")
end

### Relocation support

immutable Relocations
    h::COFFHandle
    machine::Int
    sect::SectionHeader
end

immutable RelocationRef
    h::COFFHandle
    machine::Int
    reloc::RelocationEntry
end

show(io::IO, x::RelocationRef) = show(io,x.reloc; machine=x.machine, syms=Symbols(x.h), strtab=strtab(x.h))

Relocations(s::SectionRef) = Relocations(s.handle,readheader(s.handle).Machine,s.header)

length(s::Relocations) = s.sect.NumberOfRelocations
const RelocationEntrySize = sizeof(RelocationEntry)
function getindex(s::Relocations,n)
    if n < 1 || n > length(s)
        throw(BoundsError())
    end
    offset = s.sect.PointerToRelocations + (n-1)*RelocationEntrySize
    seek(s.h,offset)
    RelocationRef(s.h,s.machine,unpack(s.h, RelocationEntry))
end

start(s::Relocations) = 1
done(s::Relocations,n) = n > length(s)
next(s::Relocations,n) = (x=s[n];(x,n+1))

printtargetsymbol(io::IO,reloc::RelocationEntry, syms, strtab) = print(io,symname(syms[reloc.SymbolTableIndex+1]; strtab = strtab, showredirect = false))

function printRelocationInterpretation(io::IO, reloc::RelocationEntry, LocalValue::UInt64, machine, syms, sects, strtab)
    if machine == IMAGE_FILE_MACHINE_AMD64
        if reloc.Type == IMAGE_REL_AMD64_ABSOLUTE
            print(io,"0x",hex(LocalValue))
        elseif reloc.Type == IMAGE_REL_AMD64_ADDR64
            print(io,"(uint64_t) ")
            printtargetsymbol(io, reloc, syms, strtab)
            print(io," + 0x",hex(LocalValue,16))
        elseif reloc.Type == IMAGE_REL_AMD64_ADDR32
            print(io,"(uint32_t) ")
            printtargetsymbol(io, reloc, syms, strtab)
            print(io," + 0x",hex(LocalValue,8))
        elseif reloc.Type == IMAGE_REL_AMD64_ADDR32NB
            print(io,"(uint32_t) ")
            printtargetsymbol(io, reloc, syms, strtab)
            print(io," + 0x",hex(LocalValue,8))
            print(io," - ImageBase")
        elseif reloc.Type >= IMAGE_REL_AMD64_REL32 && reloc.Type <= IMAGE_REL_AMD64_REL32_5
            print(io,"(uint32_t) @pc-")
            printtargetsymbol(io, reloc, syms, strtab)
            print(io," + 0x",hex(LocalValue,8))
            add = reloc.Type-IMAGE_REL_AMD64_REL32
            print(io," + ",add)
        elseif reloc.Type == IMAGE_REL_AMD64_SECTION
            print(io,"(uint16_t) indexof(")
            printtargetsymbol(io, reloc, syms, strtab)
            print(io,")")
            LocalValue != 0 && print(io,"+",dec(LocalValue))
        elseif reloc.Type == IMAGE_REL_AMD64_SECREL
            print(io,"(uint32_t) ")
            printtargetsymbol(io, reloc, syms, strtab)
            print(io," - ")
            # Get the symbol's section
            sect = sects[syms[reloc.SymbolTableIndex].entry.SectionNumber]
            print(io,sectionname(sect))
        else
            error("Unsupported Relocations")
        end
    else
        error("Relocation Support not implemented for this Machine Type")
    end
end

function relocationLength(reloc::RelocationEntry)
    reloc.Type == IMAGE_REL_AMD64_ABSOLUTE ? 0 :
    reloc.Type == IMAGE_REL_AMD64_ADDR64 ? 8 :
    reloc.Type >= IMAGE_REL_AMD64_ADDR32 &&
    reloc.Type <= IMAGE_REL_AMD64_REL32_5 ? 4 :
    reloc.Type == IMAGE_REL_AMD64_SECTION ? 2 :
    reloc.Type == IMAGE_REL_AMD64_SECREL ? 4 :
    reloc.Type == IMAGE_REL_AMD64_SECREL7 ? 1 :
    error("Unknown relocation type")
end

function inspectRelocations(sect::SectionRef, relocs = Relocations(sect))
    data = readbytes(sect);
    handle = sect.handle
    header = readheader(handle)
    for x in relocs
      offset = x.reloc.VirtualAddress - sect.header.VirtualAddress
      size = COFF.relocationLength(x.reloc)
      # zext
      # + 1 for 1-indexed array
      Ld = data[offset+1:offset+size]
      Local = reinterpret(UInt64,vcat(Ld,zeros(UInt8,sizeof(UInt64)-size)))[1]
      print("*(",sectionname(sect),"+0x",hex(offset,8),") = ")
      COFF.printRelocationInterpretation(STDOUT, x.reloc, Local, header.Machine, Symbols(handle), Sections(handle), COFF.strtab(handle))
      println()
    end
end

import Base: readbytes

readbytes{T<:IO}(io::COFFHandle{T},sec::SectionHeader) = (seek(io,sec.PointerToRawData); readbytes(io, sec.SizeOfRawData))
readbytes(sec::SectionRef) = readbytes(sec.handle,sec.header)

function sectionsize(sec::SectionHeader)
    return sec.VirtualSize == 0 ?
        sec.SizeOfRawData : min(sec.SizeOfRawData, sec.VirtualSize)
end
sectionoffset(sec::SectionHeader) = sec.PointerToRawData

### .idata parsing
@struct immutable ImportDirectoryEntry
    ImportLookupTableRVA::UInt32
    Timestamp::UInt32
    ForwarderChain::UInt32
    NameRVA::UInt32
    ImportAddressTableRVA::UInt32
end

immutable ImportDirectoryEntries
    idata::SectionRef
end
# Does no bounds checking
function Base.getindex(entries::ImportDirectoryEntries, i)
    seek(entries.idata, (i-1)*sizeof(ImportDirectoryEntry))
    unpack(entries.idata.handle, ImportDirectoryEntry)
end
Base.start(it::ImportDirectoryEntries) = 1
Base.next(it::ImportDirectoryEntries,i) = (it[i], i+1)
function Base.done(it::ImportDirectoryEntries,i)
    x = it[i]
    x.ImportLookupTableRVA == 0 && x.Timestamp == 0 && x.ForwarderChain == 0 &&
        x.ImportAddressTableRVA == 0
end
Base.iteratorsize(it::ImportDirectoryEntries) = Base.SizeUnknown()

### DWARF support

using DWARF

function debugsections{T<:IO}(h::COFFHandle{T})
    sects = collect(Sections(h))
    snames = map(sectionname,sects)
    sections = Dict{String,SectionRef}()
    for i in 1:length(snames)
        # remove leading "."
        ind = findfirst(DWARF.DEBUG_SECTIONS,String(snames[i])[2:end])
        if ind != 0
            sections[DWARF.DEBUG_SECTIONS[ind]] = sects[i]
        end
    end
    ObjFileBase.DebugSections(h,sections)
end

include("mingw.jl")

end # module


text:
# Defaults for atomics

function to_modelmatrix(b, scale, offset, rotation)
    lift_node(scale, offset, rotation) do s, o, r
        q = Quaternion(r[4], r[1], r[2], r[3])
        transformationmatrix(o, s, q)
    end
end

@default function shared(scene, kw_args)
    visible = to_bool(visible)
    scale = to_scale(scale)
    offset = to_offset(offset)
    rotation = to_rotation(rotation)
    model = to_modelmatrix(scale, offset, rotation)
    camera = to_camera(camera)
    show = to_bool(show)
    light = to_static_vec(light)
end

# Note that this will create a function called surface_defaults
# This is not perfect, but for integrating this into the scene, it's the easiest to
# just have the default function name in the macro match the drawing function name.
@default function surface(scene, kw_args)
    x = to_surface(x)
    y = to_surface(y)
    z = to_surface(z)
    xor(
        begin # Colormap is first, so it will default to it
            colormap = to_colormap(colormap)
            # convert function should only have one argument right now, so we create this closure
            colornorm = ((s, colornorm) -> to_colornorm(s, colornorm, z))(colornorm)
        end,
        begin
            color = to_color(color)
        end,
        if (image,)
            image = to_image(image)
        end
    )
end

@default function lines(scene, kw_args)
    xor(
        begin
            positions = to_positions(positions)
        end,
        if (x, y, z)
            x = to_array(x)
            y = to_array(y)
            z = to_array(z)
            positions = to_positions((x, y, z))
        end,
        if (x, y)
            x = to_array(x)
            y = to_array(y)
            positions = to_positions((x, y))
        end
    )
    xor(
        begin
            color = to_color(color)
        end,
        begin
            colormap = to_colormap(colormap)
            intensity = to_intensity(intensity)
            colornorm = to_colornorm(colornorm, intensity)
        end
    )
    linewidth = to_float(linewidth)
    linestyle = to_linestyle(linestyle)
    drawover = to_bool(drawover)
end

@default function mesh(scene, kw_args)
    shading = to_bool(shading)
    attribute_id = to_attribut_id(attribute_id)
    color = to_color(color)

    xor(
        if (indices,)
            indices = to_faces(indices)
            xor(
                if (positions,)
                    positions = to_positions(positions)
                end,
                if (x, y, z)
                    x = to_array(x)
                    y = to_array(y)
                    z = to_array(z)
                    positions = to_positions((x, y, z))
                end
            )
            mesh = to_mesh(positions, indices, color, attribute_id)
        end,
        begin
            mesh = to_mesh(mesh)
        end
    )
end

@default function scatter(scene, kw_args)
    xor(
        begin
            positions = to_positions(positions)
        end,
        if (x, y, z)
            x = to_array(x)
            y = to_array(y)
            z = to_array(z)
            positions = to_positions((x, y, z))
        end,
        if (x, y)
            x = to_array(x)
            y = to_array(y)
            positions = to_positions((x, y))
        end
    )
    # Either you give a color, or a colormap.
    # For a colormap, you'll also need intensities
    xor(
        begin
            color = to_color(color)
        end,
        begin
            colormap = to_colormap(colormap)
            intensity = to_intensity(intensity)
            colornorm = to_colornorm(colornorm, intensity)
        end
    )
    marker = to_spritemarker(marker)

    strokecolor = to_color(strokecolor)
    strokewidth = to_float(strokewidth)

    glowcolor = to_color(glowcolor)
    glowwidth = to_float(glowwidth)

    markersize = to_markersize2d(markersize)

    rotations = to_rotations(rotations)
end

@default function meshscatter(scene, kw_args)
    xor(
        begin
            positions = to_positions(positions)
        end,
        if (x, y, z)
            x = to_array(x)
            y = to_array(y)
            z = to_array(z)
            positions = to_positions((x, y, z))
        end,
        if (x, y)
            x = to_array(x)
            y = to_array(y)
            positions = to_positions((x, y))
        end
    )
    # Either you give a color, or a colormap.
    # For a colormap, you'll also need intensities
    xor(
        begin
            color = to_color(color)
        end,
        begin
            colormap = to_colormap(colormap)
            intensity = to_intensity(intensity)
            colornorm = to_colornorm(colornorm, intensity)
        end
    )

    marker = to_mesh(marker)
    markersize = to_markersize3d(markersize)
    rotations = to_rotations(rotations)
end

@default function image(scene, kw_args)
    spatialorder = to_spatial_order(spatialorder)
    x = to_interval(x)
    y = to_interval(y)
    image = to_image(image)
end

@default function volume(scene, kw_args)
    volume = to_array(volume)
    xor(
        begin
            color = to_color(color)
        end,
        begin
            colormap = to_colormap(colormap)
            # convert function should only have one argument right now, so we create this closure
            colornorm = ((b, colornorm) -> to_colornorm(b, colornorm, volume))(colornorm)
        end
    )
    algorithm = to_volume_algorithm(algorithm)
    absorption = to_float(absorption)
    isovalue = to_float(isovalue)
    isorange = to_float(isorange)
end

@default function heatmap(scene, kw_args)
    linewidth = to_float(linewidth)
    levels = to_float(levels)
    heatmap = to_array(heatmap)

    colormap = to_colormap(colormap)
    # convert function should only have one argument right now, so we create this closure
    colornorm = ((b, colornorm) -> to_colornorm(b, colornorm, heatmap))(colornorm)
end

@default function axis(scene, kw_args)
    axisnames = to_text(axisnames)

    showticks = to_bool(showticks)
    tickfont2d = to_font(tickfont2d)
    tickfont3d = to_font(tickfont3d)
    showaxis = to_bool(showaxis)
    showgrid = to_bool(showgrid)

    scalefuncs = to_scalefunc(scalefuncs)
    gridcolors = to_color(gridcolors)
    gridthickness = to_3floats(gridthickness)
    axiscolors = to_color(axiscolors)
end

function expand_kwargs(scene, kw_args)
    # TODO get in all the shorthands from Plots.jl
    attributes = Dict{Symbol, Any}(kw_args)
    shared_defaults(scene, attributes)
end

const atomic_funcs = (
    :contour => """
        contour(x, y, z)
    Creates a contour plot of the plane spanning x::Vector, y::Vector, z::Matrix
    """,
    :image => """
        image(x, y, image) / image(image)
    Plots an image on range x, y (defaults to dimensions)
    """,
    # could be implemented via image, but might be optimized specifically by the backend
    :heatmap => """
        heatmap(x, y, values) / heatmap(values)
    Plots a image on heatmap x, y (defaults to dimensions)
    """,
    :volume => """
        volume(volume_data)
    Plots a volume
    """,
    # alternatively, mesh2d?
    # :poly => """
    # """,
    :surface => """
        surface(x, y, z)
    Plots a surface, where x y z are supposed to lie on a grid
    """,
    :lines => """
        lines(x, y, z) / lines(x, y) / lines(positions)
    Plots a connected line for each element in xyz/positions
    """,
    :linesegment => """
        linesegment(x, y, z) / linesegment(x, y) / linesegment(positions)
    Plots a line for each pair of points in xyz/positions

    ## Attributes:

    The same as for [`lines`](@ref)
    """,
    # alternatively, mesh3d? Or having only mesh instead of poly + mesh and figure out 2d/3d via dispatch
    :mesh => """
        mesh(x, y, z) / mesh(mesh_object) / mesh(x, y, z, faces) / mesh(xyz, faces)
    Plots a 3D mesh
    """,
    :scatter => """
        scatter(x, y, z) / scatter(x, y) / scatter(positions)
    Plots a marker for each element in xyz/positions
    """,
    :meshscatter => """
        meshscatter(x, y, z) / meshscatter(x, y) / meshscatter(positions)
    Plots a mesh for each element in xyz/positions
    """,
    # :text => """
    # """,
    # Doesn't really need to be an atomic, could be implemented via lines
    :wireframe => """
        wireframe(x, y, z) / wireframe(positions) / wireframe(mesh)
    Draws a wireframe either interpreted as a surface or mesh
    """,

    :legend => """
        legend(series, labels)
    creates a legend from an array of plots and labels
    """,

    :axis => """
        axis(xrange, yrange, [zrange])

    Creates a axis from a x,y,z ranges
    """
)

for (func, docstring) in atomic_funcs
    adoc = try
        f = eval(Symbol("$(func)_defaults"))
        sprint(x-> Markdown.plain(x, Docs.doc(f)))
    catch e
        ""
    end
    docstring = docstring * "\n\n ## Attributes:\n\n" * adoc
    @eval begin
        """
        $($(docstring))
        """
        function $func(a::T, args...; kw_args...) where T
            if T != Scene
                $func(get_global_scene(), a, args...; kw_args...)
            else
                # keyword argument signature should never contain the backend.
                # if so, it must come from the function defined below
                ts = join(typeof.(args), ", ")
                error("Signature $func($ts) is not implemented")
            end
        end
        function $func(scene::Scene, args...; kw_args...)
            $func(scene, args..., expand_kwargs(scene, kw_args))
        end
        export $func
    end
end

"""
Get's the function a scene got created with
"""
function plotfunction(scene::Scene)
    # TODO if we actually make new types per atomic primitive, we could do this performant and more elegantly
    scene.name == :scatter && return scatter
    scene.name == :lines && return lines
    scene.name == :linesegment && return linesegment
    scene.name == :mesh && return mesh
    scene.name == :poly && return poly
    scene.name == :surface && return surface
    scene.name == :wireframe && return wireframe
    scene.name == :image && return image
    scene.name == :heatmap && return heatmap
    scene.name == :volume && return volume
    scene.name == :meshscatter && return meshscatter
    error("$(scene.name) wasn't created by any plot command")
end

"""
Recreates a similar plot with new data
"""
function Base.similar(scene::Scene, newdata...; kw_args...)
    f = plotfunction(scene)
    p = parent(scene)
    attributes = expand_kwargs(scene, kw_args)
    attributes = merge(scene.data, attributes)
    for elem in (:positions, :x, :y, :z)
        delete!(attributes, elem)
    end
    f(p, newdata..., attributes)
end

for func in (:image, :heatmap, :lines, :surface)
    # Higher level atomic signatures
    @eval begin
        function $func(scene::Scene, data::AbstractMatrix, attributes::Dict)
            $func(scene, 1:size(data, 1), 1:size(data, 2), data, attributes)
        end
        function $func(scene::Scene, x::AbstractVector{T1}, y::AbstractVector{T2}, f::Function, attributes::Dict) where {T1, T2}
            if !applicable(f, x[1], y[1])
                error("You need to pass a function like f(x::$T1, y::$T2). Found: $f")
            end
            T = typeof(f(x[1], y[1]))
            z = similar(x, T, (length(x), length(y)))
            z .= f.(x, y')
            $func(scene, x, y, z, attributes)
        end
    end
end


text:
module QuantumLattices

include("Interfaces.jl")
include("Prerequisites/Prerequisites.jl")
include("Essentials/Essentials.jl")

using .Interfaces
using .Prerequisites
using .Essentials

# Interfaces
export id, value, rank, dimension
export ⊕, ⊗, ⋅, add!, sub!, mul!, div!
export expand, expand!, decompose, decompose!, permute

# Essentials
export dtype, kind, update!, reset!

# Essentials.QuantumOperators
export ID, QuantumOperator, OperatorProd, OperatorSum, Transformation, Numericalization, MatrixRepresentation, idtype, sequence, matrix, matrix!

# Essentials.QuantumNumbers
export AbelianNumber, AbelianNumbers, @abeliannumber, periods, Momentum, Momentum₁, Momentum₂, Momentum₃

# Essentials.Spatials
export distance, azimuthd, azimuth, polard, polar, volume, isparallel, isonline, isintratriangle, issubordinate, reciprocals, translate, rotate
export PID, CPID, Point, Bond, Lattice, SuperLattice, Cylinder, Bonds, Segment, BrillouinZone, ReciprocalZone, ReciprocalPath
export pidtype, rcoord, icoord, isintracell, bonds!, bonds, @line_str, @rectangle_str, @hexagon_str
export allbonds, zerothbonds, insidebonds, acrossbonds, intrabonds, interbonds

# Essentials.DegreesOfFreedom
export CompositeIID, CompositeInternal, Index, OID, Operator, Operators, IIDSpace, Hilbert
export Subscript, Subscripts, SubscriptsID, Coupling, Couplings, OIDToTuple, Table, Term, LaTeX, Boundary
export statistics, iidtype, ishermitian, ismodulatable, abbr, otype, latexformat, twist, plain, @subscript_str, @couplings

# Essentials.Frameworks
export Parameters, Generator, SimplifiedGenerator, Action, Engine, Assignment, Algorithm
export prepare!, register!, run!, rundependences!

# Essentials.QuantumSystems
## Canonical fermionic/bosonic systems
export majorana, annihilation, creation, flatex, blatex
export FID, Fock, FockCoupling, Onsite, Hopping, Pairing, Hubbard, InterOrbitalInterSpin, InterOrbitalIntraSpin, SpinFlip, PairHopping, Coulomb
export isnormalordered, @σ⁰_str, @σˣ_str, @σʸ_str, @σᶻ_str, @σ⁺_str, @σ⁻_str, @fc_str

## SU(2) spin systems
export slatex
export SID, Spin, SpinCoupling, SpinTerm, totalspin
export @heisenberg_str, @ising_str, @gamma_str, @dm_str, @sˣ_str, @sʸ_str, @sᶻ_str, @sc_str

## Phononic systems
export nlatex
export NID, Phonon, PhononCoupling, PhononKinetic, PhononPotential

## Magnon-Phonon coupled systems
export DMPhonon

end


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license

ndgrid(v::AbstractVector) = copy(v)

function ndgrid(v1::AbstractVector{T}, v2::AbstractVector{T}) where T
    m, n = length(v1), length(v2)
    v1 = reshape(v1, m, 1)
    v2 = reshape(v2, 1, n)
    (repeat(v1, 1, n), repeat(v2, m, 1))
end

function ndgrid_fill(a, v, s, snext)
    for j = 1:length(a)
        a[j] = v[div(rem(j-1, snext), s)+1]
    end
end

function ndgrid(vs::AbstractVector{T}...) where T
    n = length(vs)
    sz = map(length, vs)
    out = ntuple(i->Array{T}(uninitialized, sz), n)
    s = 1
    for i=1:n
        a = out[i]::Array
        v = vs[i]
        snext = s*size(a,i)
        ndgrid_fill(a, v, s, snext)
        s = snext
    end
    out
end