text:
export SimulatedRobot, SimulatedRobotParams

Base.@kwdef struct SimulatedRobotParams <: DeviceParams
  defaultVelocity::Vector{typeof(1.0u"mm/s")} = [10,10,10]u"mm/s"
  axisRange::Vector{Vector{typeof(1.0u"mm")}} = [[0,500],[0,400],[0,250]]u"mm"
  namedPositions::Dict{String, Vector{typeof(1.0u"mm")}} = Dict("origin" => [0,0,0]u"mm")
  movementOrder::String = "default"
  coordinateSystem::ScannerCoordinateSystem = ScannerCoordinateSystem(3)
end

SimulatedRobotParams(dict::Dict) = params_from_dict(SimulatedRobotParams, prepareRobotDict(dict))
Base.@kwdef mutable struct SimulatedRobot <: Robot
"Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::SimulatedRobotParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}

  state::RobotState = INIT
  referenced::Bool = false
  position::Vector{typeof(1.0u"mm")} = [0,0,0]u"mm"
  connected::Bool = false
end

Base.close(rob::SimulatedRobot) = nothing

_isReferenced(rob::SimulatedRobot) = rob.referenced

dof(rob::SimulatedRobot) = 3
_getPosition(rob::SimulatedRobot) = rob.position
axisRange(rob::SimulatedRobot) = rob.params.axisRange
movementOrder(rob::SimulatedRobot) = rob.params.movementOrder

function _moveAbs(rob::SimulatedRobot, pos::Vector{<:Unitful.Length}, speed::Union{Vector{<:Unitful.Velocity},Nothing})
  @debug "SimulatedRobot: Moving to position $("["*join([string(x) for x in pos], ", ")*"]")."
  vel = speed!==nothing ? speed : rob.params.defaultVelocity
  projected_time = maximum(abs.(rob.position-pos)./vel)
  sleep(ustrip(u"s",projected_time))
  rob.position = pos
  @debug "SimulatedRobot: Movement completed"
end

function _moveRel(rob::SimulatedRobot, dist::Vector{<:Unitful.Length}, speed::Union{Vector{<:Unitful.Velocity},Nothing})
  @debug "SimulatedRobot: Moving distance $("["*join([string(x) for x in dist], ", ")*"]")"
  vel = speed!==nothing ? speed : rob.params.defaultVelocity
  projected_time = maximum(abs.(dist)./vel)
  sleep(ustrip(u"s",projected_time))
  rob.position += dist
  @debug "SimulatedRobot: Movement completed"
end

function _enable(rob::SimulatedRobot)
  @debug "SimulatedRobot enabled"
end

function _disable(rob::SimulatedRobot)
  @debug "SimulatedRobot disabled"
end

function _reset(rob::SimulatedRobot)
  @debug "SimulatedRobot reset"
  rob.connected = false
end

function _setup(rob::SimulatedRobot)
  @debug "SimulatedRobot setup"
  rob.connected = true
  sleep(0.5)
  @debug "SimulatedRobot: Finished setup"
end

function _doReferenceDrive(rob::SimulatedRobot)
  @debug "SimulatedRobot: Doing reference drive"
  sleep(2)
  @debug "SimulatedRobot: Reference drive complete"
  rob.referenced = true
end


text:
using Distributions
import Distributions: ContinuousUnivariateDistribution, pdf, logpdf

struct Improper{T} <: ContinuousUnivariateDistribution
    l::T
    u::T
    function Improper{T}(l::Number, u::Number) where T<:Number
        @assert any((!isfinite(l), !isfinite(u))) "At least one of the bounds must be Inf"
        @assert l < u "Lower bound must be smaller than upper bound"
        new{Float64}(float(l), float(u))
    end
end

Improper() = Improper{Float64}(-Inf, +Inf)
Improper(a, b) = Improper{Float64}(a, b)

function insupport(d::Improper, x::Real)
    (x >= d.l && x <= d.u) ? true : false
end

Distributions.pdf(d::Improper, x::Real) = insupport(d, x) ? 1.0 : 0.0
Distributions.logpdf(d::Improper, x::Real) = insupport(d, x) ? 0.0 : Inf
Base.minimum(d::Improper) = d.l
Base.maximum(d::Improper) = d.u
export Improper

ext: jl

max_stars_repo_name: gragusa/ProductDistributions

max_issues_repo_name: gragusa/ProductDistributions

max_forks_repo_name: gragusa/ProductDistributions


text:
abstract type AbstractPardisoLU{Tv,Ti} <: AbstractLUFactorization{Tv,Ti} end

mutable struct PardisoLU{Tv, Ti} <: AbstractPardisoLU{Tv,Ti}
    A::Union{ExtendableSparseMatrix{Tv,Ti},Nothing}
    ps::Pardiso.PardisoSolver
    phash::UInt64
end

function PardisoLU{Tv,Ti}(;iparm=nothing,dparm=nothing,mtype=nothing) where {Tv,Ti}
    fact=PardisoLU{Tv,Ti}(nothing,Pardiso.PardisoSolver(),0)
    default_initialize!(fact,iparm,dparm,mtype)
end

"""
```
PardisoLU(;valuetype=Float64,
           indextype=Int64,
           iparm::Vector,
           dparm::Vector,
           mtype::Int)

PardisoLU(matrix; iparm,dparm,mtype)
```

LU factorization based on pardiso. For using this, you need to issue `using Pardiso`
and have the pardiso library from  [pardiso-project.org](https://pardiso-project.org)
[installed](https://github.com/JuliaSparse/Pardiso.jl#pardiso-60).

The optional keyword arguments `mtype`, `iparm`  and `dparm` are
(Pardiso internal parameters)[https://github.com/JuliaSparse/Pardiso.jl#readme].

Forsetting them, one can also access the `PardisoSolver` e.g. like
```
using Pardiso
plu=PardisoLU()
Pardiso.set_iparm!(plu.ps,5,13.0)
```
"""
PardisoLU(;valuetype::Type=Float64, indextype::Type=Int64, kwargs...)=PardisoLU{valuetype,indextype}(;kwargs...)

#############################################################################################
mutable struct MKLPardisoLU{Tv, Ti} <: AbstractPardisoLU{Tv,Ti}
    A::Union{ExtendableSparseMatrix{Tv,Ti},Nothing}
    ps::Pardiso.MKLPardisoSolver
    phash::UInt64
end

function MKLPardisoLU{Tv,Ti}(;iparm=nothing,mtype=nothing) where {Tv,Ti}
    fact=MKLPardisoLU{Tv,Ti}(nothing,Pardiso.MKLPardisoSolver(),0)
    default_initialize!(fact, iparm,nothing,mtype)
end

"""
```
MKLPardisoLU(;valuetype=Float64,
           indextype=Int64,
           iparm::Vector,
           mtype::Int)

MKLPardisoLU(matrix; iparm, mtype)
```

LU factorization based on pardiso. For using this, you need to issue `using Pardiso`.
This version  uses the early 2000's fork in Intel's MKL library.

The optional keyword arguments `mtype` and `iparm`  are
(Pardiso internal parameters)[https://github.com/JuliaSparse/Pardiso.jl#readme].

For setting them you can also access the `PardisoSolver` e.g. like
```
using Pardiso
plu=MKLPardisoLU()
Pardiso.set_iparm!(plu.ps,5,13.0)
```
"""
MKLPardisoLU(;valuetype::Type=Float64, indextype::Type=Int64,kwargs...)=MKLPardisoLU{valuetype,indextype}(;kwargs...)

##########################################################################################
function default_initialize!(fact::AbstractPardisoLU{Tv,Ti}, iparm,dparm,mtype) where {Tv, Ti}
    if !isnothing(mtype)
        my_mtype=mtype
    elseif Tv<:Complex
        my_mtype=Pardiso.COMPLEX_NONSYM
    else
        my_mtype=Pardiso.REAL_NONSYM
    end

    Pardiso.set_matrixtype!(fact.ps,Pardiso.REAL_NONSYM)

    if !isnothing(iparm)
        for i=1:min(length(iparm),length(fact.ps.iparm))
            Pardiso.set_iparm!(fact.ps,i,iparm[i])
        end
    end

    if !isnothing(dparm)
        for i=1:min(length(dparm),length(fact.ps.dparm))
            Pardiso.set_dparm!(fact.ps,i,dparm[i])
        end
    end
    fact
end

function update!(lufact::AbstractPardisoLU{Tv,Ti}) where {Tv, Ti}
    ps=lufact.ps
    flush!(lufact.A)
    Acsc=lufact.A.cscmatrix
    if lufact.phash!=lufact.A.phash
        Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
        Pardiso.pardiso(ps, Tv[], Acsc, Tv[])
        Pardiso.set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT)
        lufact.phash=lufact.A.phash
    else
        Pardiso.set_phase!(ps, Pardiso.NUM_FACT)
    end
    Pardiso.fix_iparm!(ps, :N)
    Pardiso.pardiso(ps, Tv[], Acsc, Tv[])
    lufact
end

function LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, lufact::AbstractPardisoLU, v::AbstractArray{T,1} where T)
    ps=lufact.ps
    Acsc=lufact.A.cscmatrix
    Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
    Pardiso.pardiso(ps, u, Acsc, v)
    u
end

LinearAlgebra.ldiv!(fact::AbstractPardisoLU, v::AbstractArray{T,1} where T)=ldiv!(v,fact,copy(v))

@eval begin
    @makefrommatrix PardisoLU
    @makefrommatrix MKLPardisoLU
end


text:
# main.jl
using ProgressMeter
# Include helper functions
include("deps.jl")

# in_path = "images/sukii.jpg"
in_path = "images/lawrence.jpg"
# Constants:
shape = "Ellipse"
num_shape = 2500
pick_from_n = 12
max_age_of_shape = 15

in_img = loadFloatImage(in_path)

m, n = size(in_img)
bounds = [m, n]

# Initialize the background color to the average of the image.
mean_color = mean(in_img)
canvas = fill(mean_color, size(in_img))

@showprogress 1 "Fitting Shapes..." for i = 1:num_shape
avail_shapes = Array{primitive}(pick_from_n)
avail_canvases = fill(canvas, pick_from_n)
avail_score = Array{Float64}(pick_from_n)
# Generate N random shapes, pick the best one, and start hill climbing
    for j = 1:pick_from_n
        avail_shapes[j] = genRandShape(shape, bounds)
        avail_canvases[j] = applyMask(avail_canvases[j], avail_shapes[j], in_img)
        avail_score[j] = MSE(in_img, avail_canvases[j])
    end
    best_error, s_index = findmin(avail_score)
    # extract the "best" shape
    cur_shape = avail_shapes[s_index]
    # cur_canvas = avail_canvases[s_index]
    # Mutate + Hill Climb
    canvas = hillClimb(canvas, cur_shape, in_img, max_age_of_shape)

end
# imshow(canvas)
save("images/2500Ellipse.png", canvas)

max_stars_repo_name: mcbaron/primitive

max_issues_repo_name: mcbaron/primitive

max_forks_repo_name: mcbaron/primitive


text:
function transform(subject_number, magic, loop_size)
    v = 1
    for i in 1:loop_size
        v *= subject_number
        v %= magic
    end
    return v
end

function find_loop_size(subject_number, magic, target)
    v = 1
    i = 1
    while true
        v *= subject_number
        v %= magic
        v == target && return i
        i += 1
    end
end

function part1()
    magic = 20201227
    card_loop_size = find_loop_size(7, magic, 12092626)
    door_loop_size = find_loop_size(7, magic, 4707356)
    transform(12092626, magic, door_loop_size)
end

# No part2 :-) I had to finish Day 20 to get the final golden star!

# ---- Iteration ----

struct CryptoTransform{T <: Integer}
    subject::T
    magic::T
end

function Base.iterate(ct::CryptoTransform, state = nothing)
    if state === nothing
        return 1, 1
    else
        state = state * ct.subject % magic
        return state, state
    end
end

Base.IteratorSize(::Type{CryptoTransform}) = Base.IsInfinite()

transform_iter(ct::CryptoTransform, loop_size::Int) =
    first(x for (i,x) in enumerate(ct) if i == loop_size)

find_loop_size_iter(ct::CryptoTransform, public_key::Int) =
    first(i for (i,x) in enumerate(ct) if x == public_key)

function part1_iter()
    magic = 20201227
    ct = CryptoTransform(7, magic)
    card_loop_size = find_loop_size_iter(ct, 12092626)
    door_loop_size = find_loop_size_iter(ct, 4707356)

    search = CryptoTransform(12092626, magic)
    transform_iter(search, door_loop_size)
end

# ---- Functional style ----

using Lazy

# reducer
calc(subject::Int, magic::Int) = (x, _) -> subject * x % magic

# apply reducer many times
trans(subject::Int, magic::Int, n::Int) = foldl(calc(subject, magic), 1:n; init = 1)

# recursive until target is hit
@rec function loopsize(subject::Int, magic::Int, target::Int, x::Int, n::Int)
    v = calc(subject, magic)(x, n)
    v == target ? n : loopsize(subject, magic, target, v, n + 1)
end

function part1_lazy()
    magic = 20201227
    card_loop_size = loopsize(7, magic, 12092626, 1, 1)
    door_loop_size = loopsize(7, magic, 4707356, 1, 1)
    trans(12092626, magic, door_loop_size)
end

# ---- Using channel (credit: Doug) ----

function transform_subject_number(subjnum = 7, buffer_size = 10000)
    x = 1
    Channel{Int}(buffer_size) do c
        while true
            x = (x * subjnum) % 20201227
            put!(c, x)
        end
    end
end

function loopsize_doug(key)
    return first(i for (i, v) = enumerate(transform_subject_number())
                        if v == key)
end

function transform_doug(key, loop_size)
    return first(v for (i, v) = enumerate(transform_subject_number(key))
                        if i == loop_size)
end

function part1_doug()
    card_loop_size = loopsize_doug(12092626)
    door_loop_size = loopsize_doug(4707356)
    transform_doug(12092626, door_loop_size)
end


text:
"The set of all natural numbers."
struct NaturalNumbers <: Domain{Int}
end

in(x, d::NaturalNumbers) = isinteger(x) && (x >= 0)
in(x::AbstractArray, d::NaturalNumbers) = false
==(::NaturalNumbers, ::NaturalNumbers) = true

approx_in(x::Real, d::NaturalNumbers, tol) =
    (abs(x-round(x)) < tol) && (round(Int, x) ∈ d)

"The set of all integers."
struct Integers <: Domain{Int}
end

in(x, d::Integers) = isinteger(x)
in(x::AbstractArray, d::Integers) = false
==(::Integers, ::Integers) = true

approx_in(x::Real, d::Integers, tol) =
    (abs(x-round(x)) < tol) && (round(Int, x) ∈ d)

"The set of all real numbers."
struct RealNumbers <: Domain{Float64}
end

in(x::Number, d::RealNumbers) = isreal(x)
# isreal also allows real arrays, we want to disallow that here:
in(x::AbstractArray, d::RealNumbers) = false

approx_in(x::Complex, d::RealNumbers, tol) = (imag(x) < tol) && (real(x) ∈ d)
==(::RealNumbers, ::RealNumbers) = true

"The set of all rationals."
struct Rationals <: Domain{Rational{Int}}
end

in(x, d::Rationals) = x ∈ Integers()
in(x::Rational, d::Rationals) = true
==(::Rationals, ::Rationals) = true

"The set of all complex numbers whose real and imaginary parts are real numbers."
struct ComplexNumbers <: Domain{Complex{Float64}}
end

in(x::Complex{T}, d::ComplexNumbers) where {T} =
    isreal(real(x)) && isreal(imag(x))
in(x::Complex{T}, d::ComplexNumbers) where {T<:Real} = true
in(x, d::ComplexNumbers) = x ∈ RealNumbers()
==(::ComplexNumbers, ::ComplexNumbers) = true

"The set of natural numbers."
const ℕ = NaturalNumbers()
"The set of integers."
const ℤ = Integers()
"The set of rational numbers."
const ℚ = Rationals()
"The set of real numbers."
const ℝ = RealNumbers()
"The set of complex numbers."
const ℂ = ComplexNumbers()

"The space ℝ^1."
const ℝ1 = VcatDomain(ℝ)
"The space ℝ^2."
const ℝ2 = VcatDomain(ℝ, ℝ)
"The space ℝ^3."
const ℝ3 = VcatDomain(ℝ, ℝ, ℝ)
"The space ℝ^4."
const ℝ4 = VcatDomain(ℝ, ℝ, ℝ, ℝ)


text:
isnotfinite(x::Float64) = reinterpret(UInt64,x) & 0x7ff0000000000000 === 0x7ff0000000000000
isqnan(x::Float64) = reinterpret(UInt64,x) & 0x7ff8000000000000 === 0x7ff8000000000000
isainf(x::Float64) = reinterpret(UInt64,x) & 0x7ff8000000000000 === 0x7ff0000000000000

isnotfinite(x::Float32) = reinterpret(UInt32,x) & 0x7f800000 === 0x7f800000
isqnan(x::Float32) = reinterpret(UInt32,x)      & 0x7fc00000 === 0x7fc00000
isainf(x::Float32) = reinterpret(UInt32,x)      & 0x7fc00000 === 0x7f800000

isnotfinite(x::DoubleFloat{T}) where {T<:IEEEFloat} = isnan(LO(x))
isqnan(x::DoubleFloat{T}) where {T<:IEEEFloat} = isqnan(HI(x))
isainf(x::DoubleFloat{T}) where {T<:IEEEFloat} = isainf(HI(x))

isnotfinite(x::AbstractFloat) = !isfinite(x)
isqnan(x::AbstractFloat) = isnan(x)
isainf(x::AbstractFloat) = isinf(x)

include("ops/op_fp_dd.jl")
include("ops/op_fpfp_dd.jl")
include("ops/op_fpdd_dd.jl")
include("ops/op_dd_dd.jl")
include("ops/op_ddfp_dd.jl")
include("ops/op_ddsi_dd.jl")
include("ops/op_dddd_dd.jl")
include("ops/op_fp_db.jl")
include("ops/op_db_db.jl")
include("ops/op_fpfp_db.jl")
include("ops/op_fpdb_db.jl")
include("ops/op_dbfp_db.jl")
include("ops/op_dbsi_db.jl")
include("ops/op_dbdb_db.jl")
include("ops/arith.jl")


text:
#
# hierarchicalsolve
#

# Classical matrix case

\{S<:AbstractMatrix,U<:LowRankMatrix}(H::HierarchicalMatrix{S,U},f::AbstractVecOrMat) = hierarchicalsolve(H,f)

hierarchicalsolve(H::AbstractMatrix,f::AbstractVecOrMat) = H\f

function hierarchicalsolve{S<:AbstractMatrix,U<:LowRankMatrix}(H::HierarchicalMatrix{S,U},f::AbstractVecOrMat)
    T,nf = promote_type(eltype(H),eltype(f)),size(f,2)

    # Pre-compute Factorization

    !isfactored(H) && factorize!(H)

    # Partition HierarchicalMatrix

    (H11,H22),(H21,H12) = partition(H)

    # Off-diagonal low-rank matrix assembly

    U12,V12 = H12.U,H12.V
    U21,V21 = H21.U,H21.V

    # Partition the right-hand side

    (f1,f2) = (f[1:size(H11,1),:],f[1+size(H11,1):end,:])

    # Solve recursively

    H22f2,H11f1 = hierarchicalsolve(H22,f2),hierarchicalsolve(H11,f1)

    # Compute pivots

    v12,v21 = computepivots(V12,V21,H11f1,H22f2,H.factorization,nf)

    # Solve again with updated right-hand sides

    RHS1 = f1-U12*v12
    RHS2 = f2-U21*v21
    reshape([hierarchicalsolve(H11,RHS1);hierarchicalsolve(H22,RHS2)],size(f))
end

function factorize!{S<:AbstractMatrix,U<:LowRankMatrix}(H::HierarchicalMatrix{S,U})
    # Partition HierarchicalMatrix

    (H11,H22),(H21,H12) = partition(H)

    # Off-diagonal low-rank matrix assembly

    U12,V12 = H12.U,H12.V
    U21,V21 = H21.U,H21.V

    # Solve recursively

    H22U21,H11U12 = hierarchicalsolve(H22,U21),hierarchicalsolve(H11,U12)

    # Compute A

    r1,r2 = size(V12,2),size(V21,2)
    for i=1:r1,j=1:r2
        H.A[i,j+r1] += dot(V12[:,i],H22U21[:,j])
        H.A[j+r2,i] += dot(V21[:,j],H11U12[:,i])
    end

    # Compute factorization

    H.factorization = pivotldufact(H.A,r1,r2)#lufact(H.A)
    H.factored = true
end

function computepivots{V1,V2,A1,A2,T}(V12::Matrix{V1},V21::Matrix{V2},H11f1::Matrix{A1},H22f2::Matrix{A2},A::Factorization{T},nf::Int)
    P = promote_type(V1,V2,A1,A2)
    r1,r2 = size(V12,2),size(V21,2)
    b = zeros(P,r1+r2,nf)
    for i=1:nf
        for j=1:r1
            b[j,i] += dot(V12[:,j],H22f2[:,i])
        end
        for j=1:r2
            b[j+r1,i] += dot(V21[:,j],H11f1[:,i])
        end
    end
    vc = A\b
    vc[1:r1,1:nf],vc[r1+1:r1+r2,1:nf]
end

function computepivots{T}(V12::Matrix{T},V21::Matrix{T},H11f1::Vector{T},H22f2::Vector{T},A::PivotLDU{T},nf::Int)
    r1,r2 = size(V12,2),size(V21,2)
    b1,b2 = zeros(T,r1),zeros(T,r2)
    for j=1:r1
        b1[j] += dot(V12[:,j],H22f2)
    end
    for j=1:r2
        b2[j] += dot(V21[:,j],H11f1)
    end
    A_ldiv_B1B2!(A,b1,b2)
end

function computepivots{T}(V12::Matrix{T},V21::Matrix{T},H11f1::Matrix{T},H22f2::Matrix{T},A::PivotLDU{T},nf::Int)
    r1,r2 = size(V12,2),size(V21,2)
    b1,b2 = zeros(T,r1,nf),zeros(T,r2,nf)
    for i=1:nf
        for j=1:r1
            b1[j,i] += dot(V12[:,j],H22f2[:,i])
        end
        for j=1:r2
            b2[j,i] += dot(V21[:,j],H11f1[:,i])
        end
    end
    A_ldiv_B1B2!(A,b1,b2)
end


text:
#    _|_|    _|                    _|                                      _|      _|_|_|    _|_|_|    _|_|_|_|
#  _|    _|  _|_|_|      _|_|_|  _|_|_|_|  _|  _|_|    _|_|_|    _|_|_|  _|_|_|_|  _|    _|  _|    _|  _|
#  _|_|_|_|  _|    _|  _|_|        _|      _|_|      _|    _|  _|          _|      _|_|_|    _|    _|  _|_|_|
#  _|    _|  _|    _|      _|_|    _|      _|        _|    _|  _|          _|      _|        _|    _|  _|
#  _|    _|  _|_|_|    _|_|_|        _|_|  _|          _|_|_|    _|_|_|      _|_|  _|        _|_|_|    _|

abstract type AbstractPDF{N} <: AbstractFunctionWithParameters end  # N is the dimension

# calls
function (d::AbstractPDF)(x; p=freepars(d))
    normalization = normalizationintegral(d; p)
    if normalization ≈ 0.0
        println("Error: normalization ≈ 0!")
        normalization = 1.0
    end
    allp = p+fixedpars(d)
    return func(d,x; p=allp) / normalization
end

# 1 dims
function normalizationintegral(d::AbstractPDF{1}; p=freepars(d))
    allp = p+fixedpars(d)
    quadgk(x->func(d, x; p=allp), lims(d)...)[1]
end
#
# 2 dims
# function normalizationintegral(d::AbstractPDF{2}; p=freepars(d))
#     xmap = x->mapx_to_unit(x,lims(d))
#     curhe((x,f)->f[1]=func(d, xmap(x); p=p), lims(d)...)[1]
# end

"""
noparsnormf(d::AbstractPDF; p=pars(d))

Returns a single-argument lambda-function with parameters fixed to `p` and normalization computed.
"""
noparsnormf(d::AbstractPDF; p=pars(d)) = (ns=normalizationintegral(d;p=p); (x;kw...)->func(d,x;p=p)/ns)

# by default, I assume that the fields "lims" is present
lims(d::AbstractPDF) = getfield(d, :lims)
updatevalueorflag(d::AbstractPDF, s::Symbol, isfree::Bool, v=getproperty(pars(d),s)) =
    typeof(d)(updatevalueorflag(d.p, s, isfree, v), d.lims)

# other methods
function integral(d::AbstractPDF{1}, lims; p=freepars(d))
    allpars = p+fixedpars(d)
    quadgk(x->func(d,x; p=allpars), lims...)[1] / normalizationintegral(d; p=allpars)
end

####################################################################

struct Normalized{T<:AbstractFunctionWithParameters,L} <: AbstractPDF{1}
    lineshape::T
    lims::L
end
lineshape(d::Normalized) = getfield(d, :lineshape)

# two methods to be defined
import Base: getproperty
getproperty(d::Normalized, sym::Symbol) = sym==:p ? pars(lineshape(d)) : getfield(d, sym)
func(d::Normalized, x::NumberOrTuple; p=pars(d)) = func(lineshape(d), x; p)
pars(d::Normalized, isfree::Bool) = pars(lineshape(d), isfree)
updatevalueorflag(d::Normalized, s::Symbol, isfree::Bool, v=getproperty(pars(d),s)) =
    Normalized(updatevalueorflag(lineshape(d), s, isfree, v), d.lims)

# short cuts
# 1 argument
Normalized(f;p,lims) = Normalized(FunctionWithParameters(f; p), lims)

######################################################################

fixedshapepdf(f, lims) = Normalized((x;p)->f(x); lims=lims, p=∅)

######################################################################

#    _|                                                    _|      _|_|
#  _|_|_|_|  _|    _|  _|_|_|      _|_|    _|_|_|      _|_|_|    _|
#    _|      _|    _|  _|    _|  _|_|_|_|  _|    _|  _|    _|  _|_|_|_|
#    _|      _|    _|  _|    _|  _|        _|    _|  _|    _|    _|
#      _|_|    _|_|_|  _|_|_|      _|_|_|  _|_|_|      _|_|_|    _|
#                  _|  _|                  _|
#              _|_|    _|                  _|

"""
    @makepdftype MyPDF(x;p) = unnormdensity(x, p.a, p.b)

    Expected form of the expression is `f(x;p)` on the left
"""
macro makepdftype(ex)
    #
    fpx = ex.args[1].args
    name = fpx[1]
    p = fpx[2].args[1]
    (p != :p) && error("expected format f(x;p) = ... " )
    x = fpx[3]
    body = ex.args[2]
    #
    quote
        struct $name{T,N} <: AbstractPDF{1}
            p::T
            lims::N
        end
        $(esc(name))(;p,lims) = $(esc(name))(p, lims)

        import AlgebraPDF: func, pars, updatevalueorflag
        #
        $(esc(:func))(d::$(esc(name)), $(esc(x))::NumberOrTuple; p=$(esc(:pars))(d)) = $(esc(body))
    end
end


text:
using Plots
plotly()

using simulator.fresnelltools
using simulator.materials.spesifics

θs = 0:1e-3:π/2
θ = 0.82
λs = 300e-9:1e-9:1200e-9
λ = 797e-9
ds = 1e-9:1e-9:100e-9
d = 53e-9

n_gold = Au
n_fiber = 1.4
n_air = 1

R_te(λ, θ, d) = abs((ThreeLayerSystem(n_fiber, n_gold, n_air, λ, θ, d)[1] * [1, 0])[2])^2
R_tm(λ, θ, d) = abs((ThreeLayerSystem(n_fiber, n_gold, n_air, λ, θ, d)[2] * [1, 0])[2])^2
R_ratio(λ, θ, d) = R_te(λ, θ, d) / R_tm(λ, θ, d)

function Inside(λ, θ, d)
    border_1_te, border_1_tm, θg = FresnellBoundrary(n_fiber, n_gold(λ), θ)
    border_2_te, border_2_tm, θa = FresnellBoundrary(n_gold(λ), n_air, θg)

    bulk_fiber(d) = FresnellSlab(n_fiber, 2*π/λ, d, θ)
    bulk_gold(d) = FresnellSlab(n_gold(λ), 2*π/λ, d, θg)
    bulk_air(d) = FresnellSlab(n_air, 2*π/λ, d, θa)

    border_1_te_f(d) = border_1_te
    border_1_tm_f(d) = border_1_tm
    border_2_te_f(d) = border_2_te
    border_2_tm_f(d) = border_2_tm

    search = (Up, Un) -> abs(Up + Un)^2

    u_te, d_te = Interrogator(
        [bulk_fiber, border_1_te_f, bulk_gold, border_2_te_f, bulk_air],
        [20e-9, 0, d, 0, 20e-9],
        1e-10,
        search)

    u_tm, d_tm = Interrogator(
        [bulk_fiber, border_1_tm_f, bulk_gold, border_2_tm_f, bulk_air],
        [20e-9, 0, d, 0, 20e-9],
        1e-10,
        search)

    d_te, u_te, u_tm
end

# plot(λs, R_ratio.(λs, θ, d), label="θ=0.856", title="TE Reflection / TM Reflection at d=53nm")
# plot(θs, (θ) -> (R_ratio.(λ, θ + 1e-9im, d)))
# plot(ds, R_ratio.(λ, θ, ds))

# plot(λs, R_te.(λs, θ, d), label="te", title="Reflection at θ=0.82, d=53nm", legend=:bottomright)
# plot!(λs, R_tm.(λs, θ, d), label="tm", xlabel="wavelength", ylabel="reflected power")

# ds, u_te, u_tm = Inside(650e-9, θ, d)
# ds = ds .* 1e9
# _, u_tep, u_tmp = Inside(786e-9, θ, d)
#
# plot(ds, u_te, label="te normal", title="Gold thin film, power distribution",
#     xlabel="distance in nm", ylabel="all propagating power",
#     legend=:topleft)
# plot!(ds, u_tm, label="tm normal")
# plot!(ds, u_tep, label="te plasmone")
# plot!(ds, u_tmp, label="tm plasmone")


text:
using Test
import UUIDs
import Pluto: PlutoRunner, Notebook, WorkspaceManager, Cell, ServerSession, ClientSession, update_run!

@testset "Macro analysis" begin
    🍭 = ServerSession()
    🍭.options.evaluation.workspace_use_distributed = false

    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient

    @testset "Base macro call" begin
        notebook = Notebook([
            Cell("@enum Fruit 🍎 🍐"),
            Cell("my_fruit = 🍎"),
            Cell("jam(fruit::Fruit) = cook(fruit)"),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test cell(1) |> noerror
        @test [:🍎, :🍐] ⊆ notebook.topology.nodes[cell(1)].definitions
        @test :Fruit ∈ notebook.topology.nodes[cell(1)].funcdefs_without_signatures
        @test Symbol("@enum") ∈ notebook.topology.nodes[cell(1)].references

        @test cell(2) |> noerror
        @test :🍎 ∈ notebook.topology.nodes[cell(2)].references

        @test cell(3) |> noerror
        @test :Fruit ∈ notebook.topology.nodes[cell(3)].references
    end

    @testset "User defined macro 1" begin
        notebook = Notebook([
            Cell("""macro my_assign(sym, val)
              :(\$(esc(sym)) = \$(val))
            end"""),
            Cell("@my_assign x 1+1"),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test :x ∈ notebook.topology.nodes[cell(2)].definitions
        @test Symbol("@my_assign") ∈ notebook.topology.nodes[cell(2)].references

        update_run!(🍭, notebook, notebook.cells)

        # Works on second time because of old workspace
        @test :x ∈ notebook.topology.nodes[cell(2)].definitions
        @test Symbol("@my_assign") ∈ notebook.topology.nodes[cell(2)].references
    end

    @testset "User defined macro 2" begin
        notebook = Notebook([
            Cell("@my_identity(f(123))"),
            Cell(""),
            Cell(""),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        setcode(cell(2), """macro my_identity(expr)
            esc(expr)
        end""")
        update_run!(🍭, notebook, cell(2))

        setcode(cell(3), "f(x) = x")
        update_run!(🍭, notebook, cell(3))

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror
        @test cell(1).output.body == "123"

        update_run!(🍭, notebook, cell(1))
        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror
    end

    @testset "User defined macro 3" begin
        notebook = Notebook([
            Cell("""
            macro mymap()
                quote
                    [1, 2, 3] .|> sqrt
                end
            end
            """),
            Cell("@mymap")
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test cell(1) |> noerror
        @test cell(2) |> noerror

        update_run!(🍭, notebook, cell(1))

        @test cell(2) |> noerror
    end

    @testset "User defined macro 4" begin
        notebook = Notebook([
            Cell("""macro my_assign(ex)
                esc(ex)
            end"""),
            Cell("@macroexpand @my_assign 1+1"),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test Symbol("@my_assign") ∈ notebook.topology.nodes[cell(2)].references
    end

    @testset "User defined macro 5" begin
        notebook = Notebook([
            Cell("""macro dynamic_values(ex)
                [:a, :b, :c]
            end"""),
            Cell("myarray = @dynamic_values()"),
        ])
        references(idx) = notebook.topology.nodes[notebook.cells[idx]].references

        update_run!(🍭, notebook, notebook.cells)

        @test :a ∉ references(2)
        @test :b ∉ references(2)
        @test :c ∉ references(2)
    end

    @testset "User defined macro 6" begin
        notebook = Notebook([
            Cell("""macro my_macro()
                esc(:(y + x))
            end"""),
            Cell("""function my_function()
                @my_macro()
            end"""),
            Cell("my_function()"),
            Cell("x = 1"),
            Cell("y = 2"),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test [Symbol("@my_macro"), :x, :y] ⊆ notebook.topology.nodes[cell(2)].references
        @test cell(3).output.body == "3"
    end

    @testset "Function docs" begin
        notebook = Notebook([
            Cell("""
                "my function doc"
                f(x) = 2x
            """),
            Cell("f"),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test :f ∈ notebook.topology.nodes[cell(1)].funcdefs_without_signatures
        @test :f ∈ notebook.topology.nodes[cell(2)].references
    end

    @testset "Expr sanitization" begin
        struct A; end
        f(x) = x
        unserializable_expr = :($(f)(A(), A[A(), A(), A()], PlutoRunner, PlutoRunner.sanitize_expr))

        get_expr_types(other) = typeof(other)
        get_expr_types(ex::Expr) = get_expr_types.(ex.args)

        flatten(x, acc=[]) = push!(acc, x)
        function flatten(arr::AbstractVector, acc=[]) foreach(x -> flatten(x, acc), arr); acc end

        sanitized_expr = PlutoRunner.sanitize_expr(unserializable_expr)
        types = sanitized_expr |> get_expr_types |> flatten |> Set

        # Checks that no fancy type is part of the serialized expression
        @test Set([Nothing, Symbol, QuoteNode]) == types
    end

    @testset "Macrodef cells not root of run" begin
        notebook = Notebook([
            Cell(""),
            Cell(""),
            Cell(""),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test all(noerror, notebook.cells)

        setcode(cell(1), raw"""
            macro test(sym)
                esc(:($sym = true))
            end
        """)
        update_run!(🍭, notebook, cell(1))

        setcode(cell(2), "x")
        setcode(cell(3), "@test x")
        update_run!(🍭, notebook, notebook.cells[2:3])

        @test cell(2).output.body == "true"
        @test all(noerror, notebook.cells)
    end

    @testset "Reverse order" begin
        notebook = Notebook([Cell() for _ in 1:3])
        cell(idx) = notebook.cells[idx]
        update_run!(🍭, notebook, notebook.cells)

        setcode(cell(1), "x")
        update_run!(🍭, notebook, cell(1))

        @test cell(1).errored == true

        setcode(cell(2), "@bar x")
        update_run!(🍭, notebook, cell(2))

        @test cell(1).errored == true
        @test cell(2).errored == true

        setcode(cell(3), raw"""macro bar(sym)
            esc(:($sym = "yay"))
        end""")
        update_run!(🍭, notebook, cell(3))

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror
        @test cell(1).output.body == "\"yay\""
    end

    @testset "@a defines @b" begin
        notebook = Notebook([Cell() for _ in 1:4])
        cell(idx) = notebook.cells[idx]
        update_run!(🍭, notebook, notebook.cells)

        setcode(cell(1), "x")
        update_run!(🍭, notebook, cell(1))

        @test cell(1).errored == true

        setcode(cell(3), "@a()")
        setcode(cell(2), raw"""macro a()
            quote
                macro b(sym)
                    esc(:($sym = 42))
                end
            end |> esc
        end""")
        update_run!(🍭, notebook, notebook.cells[2:3])

        @test cell(1).errored == true
        @test cell(2) |> noerror
        @test cell(3) |> noerror

        setcode(cell(4), "@b x")
        update_run!(🍭, notebook, cell(4))

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror
        @test cell(4) |> noerror
        @test cell(1).output.body == "42"
    end

    @testset "Redefines macro with new SymbolsState" begin
        notebook = Notebook(Cell.([
            "@b x",
            raw"""macro b(sym)
                esc(:($sym = 42))
            end""",
            "x",
            "y",
        ]))
        cell(idx) = notebook.cells[idx]
        update_run!(🍭, notebook, notebook.cells)

        @test cell(3).output.body == "42"
        @test cell(4).errored == true

        setcode(cell(2), """macro b(_)
            esc(:(y = 42))
        end""")
        update_run!(🍭, notebook, cell(2))

        @test cell(4).output.body == "42"
        @test cell(3).errored == true

        notebook = Notebook(Cell.([
            "@b x",
            raw"""macro b(sym)
                esc(:($sym = 42))
            end""",
            "x",
            "y",
        ]))
        update_run!(🍭, notebook, notebook.cells)

        @test cell(3).output.body == "42"
        @test cell(4).errored == true

        setcode(cell(2), """macro b(_)
            esc(:(y = 42))
        end""")
        update_run!(🍭, notebook, [cell(1), cell(2)])

        # Cell 4 is executed even because cell(1) is in the root
        # of the reactive run because the expansion is done with the new version
        # of the macro in the new workspace because of the current_roots parameter of `resolve_topology`.
        # See Run.jl#resolve_topology.
        @test cell(4).output.body == "42"
        @test cell(3).errored == true
    end

    @testset "Reactive macro update does not invalidate the macro calls" begin
        notebook = Notebook(Cell.([
            raw"""macro b(sym)
                if z > 40
                    esc(:($sym = $z))
                else
                    esc(:(y = $z))
                end
            end""",
            "z = 42",
            "@b(x)",
            "x",
            "y",
        ]))
        cell(idx) = notebook.cells[idx]
        update_run!(🍭, notebook, notebook.cells)

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror
        @test cell(4) |> noerror
        @test cell(5).errored == true

        setcode(cell(2), "z = 39")

        # running only 2, running all cells here works however
        update_run!(🍭, notebook, cell(2))

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror
        @test cell(4).output.body != "42"
        @test cell(4).errored == true
        @test cell(5) |> noerror
    end

    @testset "Explicitely running macrocalls updates the reactive node" begin
        notebook = Notebook(Cell.([
            "@b()",
            "ref = Ref{Int}(0)",
            raw"""macro b()
                ex = if iseven(ref[])
                    :(x = 10)
                else
                    :(y = 10)
                end |> esc
                ref[] += 1
                ex
            end""",
            "x",
            "y",
        ]))
        cell(i) = notebook.cells[i]
        update_run!(🍭, notebook, notebook.cells)

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror
        @test cell(4) |> noerror
        @test cell(5).errored == true

        update_run!(🍭, notebook, cell(1))

        @test cell(4).errored == true
        @test cell(5) |> noerror
    end

    @testset "Implicitely running macrocalls updates the reactive node" begin
        notebook = Notebook(Cell.([
            "updater; @b()",
            "ref = Ref{Int}(0)",
            raw"""macro b()
                ex = if iseven(ref[])
                    :(x = 10)
                else
                    :(y = 10)
                end |> esc
                ref[] += 1
                ex
            end""",
            "x",
            "y",
            "updater = 1",
        ]))
        cell(i) = notebook.cells[i]
        update_run!(🍭, notebook, notebook.cells)

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror
        @test cell(4) |> noerror
        output_1 = cell(4).output.body
        @test cell(5).errored == true
        @test cell(6) |> noerror

        setcode(cell(6), "updater = 2")
        update_run!(🍭, notebook, cell(6))

        # the output of cell 4 has not changed since the underlying computer
        # has not been regenerated. To update the reactive node and macrocall
        # an explicit run of @b() must be done.
        @test cell(4).output.body == output_1
        @test cell(5).errored == true
    end

    @testset "Weird behavior" begin
        # https://github.com/fonsp/Pluto.jl/issues/1591

        notebook = Notebook(Cell.([
            "macro huh(_) throw(\"Fail!\") end",
            "huh(e) = e",
            "@huh(z)",
            "z = 101010",
        ]))
        cell(idx) = notebook.cells[idx]
        update_run!(🍭, notebook, notebook.cells)

        @test cell(3).errored == true

        setcode(cell(3), "huh(z)")
        update_run!(🍭, notebook, cell(3))

        @test cell(3) |> noerror
        @test cell(3).output.body == "101010"

        setcode(cell(4), "z = 1234")
        update_run!(🍭, notebook, cell(4))

        @test cell(3) |> noerror
        @test cell(3).output.body == "1234"
    end

    @testset "Cell failing first not re-run?" begin
        notebook = Notebook(Cell.([
            "x",
            "@b x",
            raw"macro b(sym) esc(:($sym = 42)) end",
        ]))
        update_run!(🍭, notebook, notebook.cells)

        # CELL 1 "x" was run first and failed because the definition
        # of x was not yet found. However, it was not run re-run when the definition of
        # x ("@b(x)") was run. Should it? Maybe set a higher precedence to cells that define
        # macros inside the notebook.
        @test_broken noerror(notebook.cells[1]; verbose=false)
    end

    @testset "@a defines @b initial loading" begin
        notebook = Notebook(Cell.([
            "x",
            "@b x",
            "@a",
            raw"""macro a()
                quote
                    macro b(sym)
                        esc(:($sym = 42))
                    end
                end |> esc
            end"""
        ]))
        cell(idx) = notebook.cells[idx]
        update_run!(🍭, notebook, notebook.cells)

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror
        @test cell(4) |> noerror
        @test cell(1).output.body == "42"
    end

    @testset "Macro with long compile time gets function wrapped" begin
        ms = 1e-3
        ns = 1e-9
        sleep_time = 40ms

        notebook = Notebook(Cell.([
            "updater; @b()",
            """macro b()
                x = rand()
                sleep($sleep_time)
                :(1+\$x)
            end""",
            "updater = :slow",
        ]))
        cell(idx) = notebook.cells[idx]
        update_run!(🍭, notebook, notebook.cells)

        @test noerror(cell(1))
        runtime = cell(1).runtime*ns
        output_1 = cell(1).output.body
        @test sleep_time <= runtime

        setcode(cell(3), "updater = :fast")
        update_run!(🍭, notebook, cell(3))

        @test noerror(cell(1))
        runtime = cell(1).runtime*ns
        @test runtime < sleep_time # no recompilation!

        # output is the same because no new compilation happens
        @test output_1 == cell(1).output.body

        # force recompilation by explicitely running the cell
        update_run!(🍭, notebook, cell(1))

        @test cell(1) |> noerror
        @test output_1 != cell(1).output.body
        output_3 = cell(1).output.body

        setcode(cell(1), "@b()") # changing code generates a new 💻
        update_run!(🍭, notebook, cell(1))

        @test cell(1) |> noerror
        @test output_3 != cell(1).output.body
    end

    @testset "Macro Prefix" begin
        🍭.options.evaluation.workspace_use_distributed = true

        notebook = Notebook(Cell.([
            "@sprintf \"answer = %d\" x",
            "x = y+1",
            raw"""
            macro def(sym)
                esc(:($sym=41))
            end
            """,
            "@def y",
            "import Printf: @sprintf",
        ]))
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, cell(1))

        @test cell(1).errored == true

        update_run!(🍭, notebook, cell(5))

        @test occursin("UndefVarError: x", cell(1).output.body[:msg])

        update_run!(🍭, notebook, cell(3))
        update_run!(🍭, notebook, cell(2))
        update_run!(🍭, notebook, cell(4))

        @test cell(1) |> noerror

        WorkspaceManager.unmake_workspace((🍭, notebook))
        🍭.options.evaluation.workspace_use_distributed = false
    end

    @testset "Package macro 1" begin
        notebook = Notebook([
            Cell("using Dates"),
            Cell("df = dateformat\"Y-m-d\""),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, cell(2))

        @test cell(2).errored == true
        @test occursinerror("UndefVarError: @dateformat_str", cell(2)) == true

        update_run!(🍭, notebook, notebook.cells)

        @test cell(1) |> noerror
        @test cell(2) |> noerror

        notebook = Notebook([
            Cell("using Dates"),
            Cell("df = dateformat\"Y-m-d\""),
        ])
        update_run!(🍭, notebook, notebook.cells)

        @test cell(1) |> noerror
        @test cell(2) |> noerror
    end

    @testset "Package macro 2" begin
        🍭.options.evaluation.workspace_use_distributed = true

        notebook = Notebook([
            Cell("z = x^2 + y"),
            Cell("@variables x y"),
            Cell("""
            begin
                import Pkg
                Pkg.activate(mktempdir())
                Pkg.add(Pkg.PackageSpec(name="Symbolics", version="1"))
                import Symbolics: @variables
            end
            """),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells[1:2])

        @test cell(1).errored == true
        @test cell(2).errored == true

        update_run!(🍭, notebook, cell(3))

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(2) |> noerror

        update_run!(🍭, notebook, notebook.cells)

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror

        update_run!(🍭, notebook, cell(2))

        @test cell(1) |> noerror
        @test cell(2) |> noerror

        setcode(cell(2), "@variables 🐰 y")
        update_run!(🍭, notebook, cell(2))

        @test cell(1).errored
        @test cell(2) |> noerror

        setcode(cell(1), "z = 🐰^2 + y")
        update_run!(🍭, notebook, cell(1))

        @test cell(1) |> noerror
        @test cell(2) |> noerror

        WorkspaceManager.unmake_workspace((🍭, notebook))

        🍭.options.evaluation.workspace_use_distributed = false
    end

    @testset "Previous workspace for unknowns" begin
        notebook = Notebook([
            Cell("""macro my_identity(expr)
              expr
            end"""),
            Cell("(@__MODULE__, (@my_identity 1 + 1))"),
            Cell("@__MODULE__"),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, cell(1))
        update_run!(🍭, notebook, notebook.cells[2:end])

        @test cell(1) |> noerror
        @test cell(2) |> noerror
        @test cell(3) |> noerror

        module_from_cell2 = cell(2).output.body[:elements][1][2][1]
        module_from_cell3 = cell(3).output.body

        @test module_from_cell2 == module_from_cell3
    end

    @testset "Definitions" begin
        notebook = Notebook([
            Cell("""macro my_assign(sym, val)
                :(\$(esc(sym)) = \$(val))
            end"""),
            Cell("c = :hello"),
            Cell("@my_assign b c"),
            Cell("b"),
        ])
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)
        update_run!(🍭, notebook, notebook.cells)

        @test ":hello" == cell(3).output.body
        @test ":hello" == cell(4).output.body
        @test :b ∈ notebook.topology.nodes[cell(3)].definitions
        @test [:c, Symbol("@my_assign")] ⊆ notebook.topology.nodes[cell(3)].references

        setcode(notebook.cells[2], "c = :world")
        update_run!(🍭, notebook, cell(2))

        @test ":world" == cell(3).output.body
        @test ":world" == cell(4).output.body
    end

    @testset "Macros using import" begin
        notebook = Notebook(Cell.([
            """
            @option "option_a" struct OptionA
                option_value::option_type
            end
            """,
            "option_type = String",
            "import Configurations: @option",
        ]))
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test :option_type ∈ notebook.topology.nodes[cell(1)].references
        @test cell(1) |> noerror
    end

    @testset "GlobalRefs in macros should be respected" begin
        notebook = Notebook(Cell.([
            """
            macro identity(expr)
                expr
            end
            """,
            """
            x = 20
            """,
            """
            let x = 10
                @identity(x)
            end
            """,
        ]))
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test all(cell.([1,2,3]) .|> noerror)
        @test cell(3).output.body == "20"
    end

    @testset "GlobalRefs shouldn't break unreached undefined references" begin
        notebook = Notebook(Cell.([
            """
            macro get_x_but_actually_not()
                quote
                    if false
                        x
                    else
                        :this_should_be_returned
                    end
                end
            end
            """,
            """
            @get_x_but_actually_not()
            """,
        ]))
        cell(idx) = notebook.cells[idx]

        update_run!(🍭, notebook, notebook.cells)

        @test all(cell.([1,2]) .|> noerror)
        @test cell(2).output.body == ":this_should_be_returned"
    end
end


text:
#**************************************************************************************
# test_GeoEfficiency.jl package
# ============================== part of the GeoEfficiency.jl package.
#
#
#
#**************************************************************************************

using GeoEfficiency
using Base.Test


text:
"""
    Featurization

A Featurization collectively stores a set of FeatureDescriptors (and any other
supporting attributes that it may further require), and defines how to featurize
an Atoms object by using the encoding schemes defined by the FeatureDescriptors
stored.
"""
module Featurization

using ..ChemistryFeaturization.AbstractType: AbstractFeaturization

import ..ChemistryFeaturization.encodable_elements
encodable_elements(fzn::AbstractFeaturization) = throw(MethodError(encodable_elements, fzn))
export encodable_elements

import ..ChemistryFeaturization.encode
encode(fzn::AbstractFeaturization, object_to_be_encoded) = throw(MethodError(encode, fzn))
export encode

import ..ChemistryFeaturization.decode
decode(fzn::AbstractFeaturization, encoded_feature) = throw(MethodError(decode, fzn))
export decode

include("graphnodefeaturization.jl")
export GraphNodeFeaturization, encode, decode

include("weavefeaturization.jl")
export WeaveFeaturization

end


text:
#=
LeapFrogger.jl
Jacob Child
Flow Lab Onboarding Project
April 23rd, 2022

Equations
V(vector) = (Gamma(vector) cross r(vector))/(2*pi*r^2)

Pseudocode- Using the reference image and the right hand rule, compute the velocities for each
vortex "point" individually. The bottom left point is 1, then go clockwise.
For Pt1- compute the effects of pts 2,3,4 on it, add them together for Pt1TotVeloc, then using a
small time step, compute its new location and repeat the process :||
All points have the same process
Things to watch- don't update the position of any points before all the calculations for
that time step are done.

Coordinate system:
(0; 0; 0) = [x; y; z] = inbetween and in line with pts 1 and 2
Pt 1 = btm left, the rest follow clockwise.
Righthand rule = out of page = positive, into page = negative.
=#

#needed packages
using LinearAlgebra
using Plots
#using DataFrame

#Inputs/Variables
d = 1 #distance between vortex ring top and bottom points, and between rings
G = 1   #strength of vorticity = Gamma, same for all rings/points
dt = .01 #time step size
x = collect(0:dt:480) #gives a time vector from 0-xseconds with .01 second step size
plot()  #to clear the plot

#p1 = zeros(lenght(x), 3)
#Point 1
Pt1Gamma = [0; 0; -G]
Pt1Location = [0; -d/2; 0]
x1 = []
y1 = []

#Point 2
Pt2Gamma = [0; 0; G]
Pt2Location = [0; d/2; 0]
x2 = []
y2 = []

#Point 3
Pt3Gamma = [0; 0; G]
Pt3Location = [d; d/2; 0]
x3 = []
y3 = []

#Point 4
Pt4Gamma = [0; 0; -G]
Pt4Location = [d; -d/2; 0]
x4 = []
y4 = []

#before putting everything in functions, practice for pt 1
#2 on 1
# R21 = Pt1Location - Pt2Location #make sure positive vs negative is correct!
# println(R21)
# MagnitudeR21 = sqrt(R21[1]^2 + R21[2]^2 + R21[3]^2) #the magnitude of the distance btwn the 2 pts
# V21Inf = cross(Pt2Gamma, R21) / (2*pi*MagnitudeR21) #The influence pt2 has on 1 for velocity
# println(V21Inf)
# dt = .01    #time step
# Pt1dp = V21Inf * dt
# Pt1Location = Pt1Location + Pt1dp
# println(Pt1Location)

function v_influence(Location1, Location2, Gamma)
    R = Location1 - Location2 #how to make sure this is right positive vs negative?
    MagnitudeR = sqrt(R[1]^2 + R[2]^2 + R[3]^2) #the magnitude of the distance btwn the 2 pts
    VInf = cross(Gamma, R) / (2*pi*MagnitudeR^2) #The influence pt2 has on 1 for velocity
    return VInf
end

function total_influence(VAInf,VBInf, VCInf)
    VTotInf = VAInf + VBInf +VCInf
    return VTotInf

end

function update_location(Location, VInfTot, dt)
    dp = VInfTot*dt
    NewLocation = Location + dp
    return NewLocation
end

###Looping Function

for i in x #i=1:length(x)

    #Point 1
    V21Inf = v_influence(Pt1Location, Pt2Location, Pt2Gamma)
    V31Inf = v_influence(Pt1Location, Pt3Location, Pt3Gamma)
    V41Inf = v_influence(Pt1Location, Pt4Location, Pt4Gamma)
    V1InfTot = total_influence(V21Inf, V31Inf, V41Inf)
    #println(V21Inf)

    #Point 2
    V12Inf = v_influence(Pt2Location, Pt1Location, Pt1Gamma)
    V32Inf = v_influence(Pt2Location, Pt3Location, Pt3Gamma)
    V42Inf = v_influence(Pt2Location, Pt4Location, Pt4Gamma)
    V2InfTot = total_influence(V12Inf, V32Inf, V42Inf)

    #Point 3
    V13Inf = v_influence(Pt3Location, Pt1Location, Pt1Gamma)
    V23Inf = v_influence(Pt3Location, Pt2Location, Pt2Gamma)
    V43Inf = v_influence(Pt3Location, Pt4Location, Pt4Gamma)
    V3InfTot = total_influence(V13Inf, V23Inf, V43Inf)

    #Point 4
    V14Inf = v_influence(Pt4Location, Pt1Location, Pt1Gamma)
    V24Inf = v_influence(Pt4Location, Pt2Location, Pt2Gamma)
    V34Inf = v_influence(Pt4Location, Pt3Location, Pt3Gamma)
    V4InfTot = total_influence(V14Inf, V24Inf, V34Inf)

    #Update Location of All points
    global Pt1Location = update_location(Pt1Location, V1InfTot, dt) #Make sure to update the positions last!
    # p1[i,:] = update_location(Pt1Location, V1InfTot, dt)
    global Pt2Location = update_location(Pt2Location, V2InfTot, dt)
    global Pt3Location = update_location(Pt3Location, V3InfTot, dt)
    global Pt4Location = update_location(Pt4Location, V4InfTot, dt)
    global x1 = push!(x1,Pt1Location[1]) #creates a series of x and y values for each point
    global y1 = push!(y1,Pt1Location[2])
    global x2 = push!(x2,Pt2Location[1])
    global y2 = push!(y2,Pt2Location[2])
    global x3 = push!(x3,Pt3Location[1])
    global y3 = push!(y3,Pt3Location[2])
    global x4 = push!(x4,Pt4Location[1])
    global y4 = push!(y4,Pt4Location[2])
    #x1,y1,_ = update_location(Pt1Location, V1InfTot, dt)
    # #debugging
    # A = V12Inf[1]   #x direction velocity influence of pt 3 on pt 2
    # #B = V12Inf[2]   #y direction velocity influence of pt 3 on pt2
    # display(plot!([A], legend = false, m = 3))

end
#plots
#Vortex Ring 1 (back)
display(plot!(x1,y1,  label = false, color = :blue, xlabel="Position (m)", ylabel = "Position (m)")) #plot only does vectors, not floats etc, number values are plot commands
display(plot!(x2,y2,  label = false, color = :blue, title = "Leap Frogging Vortex Ring Positions Over Time"))

#Vortex Ring 2 (front)
display(plot!(x3,y3,  label = false, color = :green))
display(plot!(x4,y4,  label = false, color = :green))

# plt = plot(xaxis="X distance (m)", yaxis="Y distance (m)")
# plot!(x1, y1)
# display(plt)

#println(Pt1Location, Pt2Location, Pt3Location, Pt4Location)


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)
module FreeType

using Compat

isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")) || error("FreeType was not build properly. Please run Pkg.build(\"FreeType\")")
include(joinpath("..", "deps", "deps.jl"))

function __init__()
    dl_ptr = Libdl.dlopen_e(freetype)
	dl_ptr == C_NULL && error("FreeType library could not be loaded correctly. Please check library path ($freetype), or run Pkg.build(\"FreeType\")")
	Libdl.dlclose(dl_ptr)
end

typealias FT_Int16 Int16
typealias FT_UInt16 UInt16
typealias FT_Int32 Cint
typealias FT_UInt32 UInt32
typealias FT_Fast Cint
typealias FT_UFast UInt32
abstract FT_MemoryRec_
typealias FT_Memory Ptr{FT_MemoryRec_}
typealias FT_Alloc_Func Ptr{Void}
typealias FT_Free_Func Ptr{Void}
typealias FT_Realloc_Func Ptr{Void}
type FT_StreamDesc
    _FT_StreamDesc::Ptr{Void}
end
typealias FT_Stream_IoFunc Ptr{Void}
typealias FT_Stream_CloseFunc Ptr{Void}
type FT_StreamRec
    base::Ptr{Cuchar}
    size::Culong
    pos::Culong
    descriptor::FT_StreamDesc
    pathname::FT_StreamDesc
    read::FT_Stream_IoFunc
    close::FT_Stream_CloseFunc
    memory::FT_Memory
    cursor::Ptr{Cuchar}
    limit::Ptr{Cuchar}
end
typealias FT_Stream Ptr{FT_StreamRec}
typealias FT_Pos Clong
immutable FT_Vector
    x::FT_Pos
    y::FT_Pos
end
immutable FT_BBox
    xMin::FT_Pos
    yMin::FT_Pos
    xMax::FT_Pos
    yMax::FT_Pos
end

const FT_LOAD_DEFAULT = @compat Int32(0)
const FT_LOAD_TARGET_NORMAL = @compat Int32(0)
const FT_LOAD_NO_SCALE = @compat Int32(1 << 0)
const FT_LOAD_NO_HINTING = @compat Int32(1 << 1)
const FT_LOAD_RENDER = @compat Int32(1 << 2)
const FT_LOAD_NO_BITMAP = @compat Int32(1 << 3)
const FT_LOAD_VERTICAL_LAYOUT = @compat Int32(1 << 4)
const FT_LOAD_FORCE_AUTOHINT = @compat Int32(1 << 5)
const FT_LOAD_CROP_BITMAP = @compat Int32(1 << 6)
const FT_LOAD_PEDANTIC = @compat Int32(1 << 7)
const FT_LOAD_ADVANCE_ONLY = @compat Int32(1 << 8)
const FT_LOAD_IGNORE_GLOBAL_ADVANCE_WIDTH = @compat Int32(1 << 9)
const FT_LOAD_NO_RECURSE = @compat Int32(1 << 10)
const FT_LOAD_IGNORE_TRANSFORM = @compat Int32(1 << 11)
const FT_LOAD_MONOCHROME = @compat Int32(1 << 12)
const FT_LOAD_LINEAR_DESIGN = @compat Int32(1 << 13)
const FT_LOAD_SBITS_ONLY = @compat Int32(1 << 14)
const FT_LOAD_NO_AUTOHINT = @compat Int32(1 << 15)
const FT_LOAD_TARGET_LIGHT = @compat Int32(1 << 16)
const FT_LOAD_TARGET_MONO = @compat Int32(1 << 17)
const FT_LOAD_TARGET_LCD = @compat Int32(196608)
const FT_LOAD_TARGET_LCD_V = @compat Int32(1 << 18)
const FT_LOAD_COLOR = @compat Int32(1 << 20)

const FT_FACE_FLAG_SCALABLE         = @compat Int32(1 <<  0)
const FT_FACE_FLAG_FIXED_SIZES      = @compat Int32(1 <<  1)
const FT_FACE_FLAG_FIXED_WIDTH      = @compat Int32(1 <<  2)
const FT_FACE_FLAG_SFNT             = @compat Int32(1 <<  3)
const FT_FACE_FLAG_HORIZONTAL       = @compat Int32(1 <<  4)
const FT_FACE_FLAG_VERTICAL         = @compat Int32(1 <<  5)
const FT_FACE_FLAG_KERNING          = @compat Int32(1 <<  6)
const FT_FACE_FLAG_FAST_GLYPHS      = @compat Int32(1 <<  7)
const FT_FACE_FLAG_MULTIPLE_MASTERS = @compat Int32(1 <<  8)
const FT_FACE_FLAG_GLYPH_NAMES      = @compat Int32(1 <<  9)
const FT_FACE_FLAG_EXTERNAL_STREAM  = @compat Int32(1 << 10)
const FT_FACE_FLAG_HINTER           = @compat Int32(1 << 11)
const FT_FACE_FLAG_CID_KEYED        = @compat Int32(1 << 12)
const FT_FACE_FLAG_TRICKY           = @compat Int32(1 << 13)
const FT_FACE_FLAG_COLOR            = @compat Int32(1 << 14)

# begin enum FT_Pixel_Mode
typealias FT_Pixel_Mode UInt32
const FT_PIXEL_MODE_NONE = @compat UInt32(0)
const FT_PIXEL_MODE_MONO = @compat UInt32(1)
const FT_PIXEL_MODE_GRAY = @compat UInt32(2)
const FT_PIXEL_MODE_GRAY2 = @compat UInt32(3)
const FT_PIXEL_MODE_GRAY4 = @compat UInt32(4)
const FT_PIXEL_MODE_LCD = @compat UInt32(5)
const FT_PIXEL_MODE_LCD_V = @compat UInt32(6)
const FT_PIXEL_MODE_BGRA = @compat UInt32(7)
const FT_PIXEL_MODE_MAX = @compat UInt32(8)
# end enum FT_Pixel_Mode
immutable FT_Bitmap
    rows::Cint
    width::Cint
    pitch::Cint
    buffer::Ptr{Cuchar}
    num_grays::Int16
    pixel_mode::UInt8
    palette_mode::UInt8
    palette::Ptr{Void}
end
immutable FT_Outline
    n_contours::Int16
    n_points::Int16
    points::Ptr{FT_Vector}
    tags::Ptr{UInt8}
    contours::Ptr{Int16}
    flags::Cint
end
typealias FT_Outline_MoveToFunc Ptr{Void}
typealias FT_Outline_LineToFunc Ptr{Void}
typealias FT_Outline_ConicToFunc Ptr{Void}
typealias FT_Outline_CubicToFunc Ptr{Void}
type FT_Outline_Funcs
    move_to::FT_Outline_MoveToFunc
    line_to::FT_Outline_LineToFunc
    conic_to::FT_Outline_ConicToFunc
    cubic_to::FT_Outline_CubicToFunc
    shift::Cint
    delta::FT_Pos
end
# begin enum FT_Glyph_Format
typealias FT_Glyph_Format UInt32
const FT_GLYPH_FORMAT_NONE = @compat UInt32(0)
const FT_GLYPH_FORMAT_COMPOSITE = @compat UInt32(1668246896)
const FT_GLYPH_FORMAT_BITMAP = @compat UInt32(1651078259)
const FT_GLYPH_FORMAT_OUTLINE = @compat UInt32(1869968492)
const FT_GLYPH_FORMAT_PLOTTER = @compat UInt32(1886154612)
# end enum FT_Glyph_Format
abstract FT_RasterRec_
typealias FT_Raster Ptr{FT_RasterRec_}
type FT_Span
    x::Int16
    len::UInt16
    coverage::Cuchar
end
typealias FT_SpanFunc Ptr{Void}
typealias FT_Raster_BitTest_Func Ptr{Void}
typealias FT_Raster_BitSet_Func Ptr{Void}
type FT_Raster_Params
    target::Ptr{FT_Bitmap}
    source::Ptr{Void}
    flags::Cint
    gray_spans::FT_SpanFunc
    black_spans::FT_SpanFunc
    bit_test::FT_Raster_BitTest_Func
    bit_set::FT_Raster_BitSet_Func
    user::Ptr{Void}
    clip_box::FT_BBox
end
typealias FT_Raster_NewFunc Ptr{Void}
typealias FT_Raster_DoneFunc Ptr{Void}
typealias FT_Raster_ResetFunc Ptr{Void}
typealias FT_Raster_SetModeFunc Ptr{Void}
typealias FT_Raster_RenderFunc Ptr{Void}
type FT_Raster_Funcs
    glyph_format::FT_Glyph_Format
    raster_new::FT_Raster_NewFunc
    raster_reset::FT_Raster_ResetFunc
    raster_set_mode::FT_Raster_SetModeFunc
    raster_render::FT_Raster_RenderFunc
    raster_done::FT_Raster_DoneFunc
end
typealias FT_Bool Cuchar
typealias FT_FWord Int16
typealias FT_UFWord UInt16
typealias FT_Char UInt8
typealias FT_Byte Cuchar
typealias FT_Bytes Ptr{FT_Byte}
typealias FT_Tag FT_UInt32
typealias FT_String UInt8
typealias FT_Short Int16
typealias FT_UShort UInt16
typealias FT_Int Cint
typealias FT_UInt UInt32
typealias FT_Long Clong
typealias FT_ULong Culong
typealias FT_F2Dot14 Int16
typealias FT_F26Dot6 Clong
typealias FT_Fixed Clong
typealias FT_Error Cint
typealias FT_Pointer Ptr{Void}
typealias FT_Offset Csize_t
typealias FT_PtrDist Cptrdiff_t
type FT_UnitVector
    x::FT_F2Dot14
    y::FT_F2Dot14
end
type FT_Matrix
    xx::FT_Fixed
    xy::FT_Fixed
    yx::FT_Fixed
    yy::FT_Fixed
end
type FT_Data
    pointer::Ptr{FT_Byte}
    length::FT_Int
end
typealias FT_Generic_Finalizer Ptr{Void}
immutable FT_Generic
    data::Ptr{Void}
    finalizer::FT_Generic_Finalizer
end
abstract FT_ListNodeRec_
typealias FT_ListNode Ptr{FT_ListNodeRec_}
type FT_ListNodeRec <: FT_ListNodeRec_
    prev::FT_ListNode
    next::FT_ListNode
    data::Ptr{Void}
end
immutable FT_ListRec
    head::FT_ListNode
    tail::FT_ListNode
end
typealias FT_List Ptr{FT_ListRec}
# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 UInt32
const FT_Mod_Err_Base = @compat UInt32(0)
const FT_Mod_Err_Autofit = @compat UInt32(0)
const FT_Mod_Err_BDF = @compat UInt32(0)
const FT_Mod_Err_Bzip2 = @compat UInt32(0)
const FT_Mod_Err_Cache = @compat UInt32(0)
const FT_Mod_Err_CFF = @compat UInt32(0)
const FT_Mod_Err_CID = @compat UInt32(0)
const FT_Mod_Err_Gzip = @compat UInt32(0)
const FT_Mod_Err_LZW = @compat UInt32(0)
const FT_Mod_Err_OTvalid = @compat UInt32(0)
const FT_Mod_Err_PCF = @compat UInt32(0)
const FT_Mod_Err_PFR = @compat UInt32(0)
const FT_Mod_Err_PSaux = @compat UInt32(0)
const FT_Mod_Err_PShinter = @compat UInt32(0)
const FT_Mod_Err_PSnames = @compat UInt32(0)
const FT_Mod_Err_Raster = @compat UInt32(0)
const FT_Mod_Err_SFNT = @compat UInt32(0)
const FT_Mod_Err_Smooth = @compat UInt32(0)
const FT_Mod_Err_TrueType = @compat UInt32(0)
const FT_Mod_Err_Type1 = @compat UInt32(0)
const FT_Mod_Err_Type42 = @compat UInt32(0)
const FT_Mod_Err_Winfonts = @compat UInt32(0)
const FT_Mod_Err_GXvalid = @compat UInt32(0)
const FT_Mod_Err_Max = @compat UInt32(1)
# end enum ANONYMOUS_1
immutable FT_Glyph_Metrics
    width::FT_Pos
    height::FT_Pos
    horiBearingX::FT_Pos
    horiBearingY::FT_Pos
    horiAdvance::FT_Pos
    vertBearingX::FT_Pos
    vertBearingY::FT_Pos
    vertAdvance::FT_Pos
end
type FT_Bitmap_Size
    height::FT_Short
    width::FT_Short
    size::FT_Pos
    x_ppem::FT_Pos
    y_ppem::FT_Pos
end
abstract FT_LibraryRec_
typealias FT_Library Ptr{FT_LibraryRec_}
abstract FT_ModuleRec_
typealias FT_Module Ptr{FT_ModuleRec_}
abstract FT_DriverRec_
typealias FT_Driver Ptr{FT_DriverRec_}
abstract FT_RendererRec_
typealias FT_Renderer Ptr{FT_RendererRec_}
abstract FT_SizeRec_
typealias FT_Size Ptr{FT_SizeRec_}
abstract FT_CharMapRec_
typealias FT_CharMap Ptr{FT_CharMapRec_}
# begin enum FT_Encoding
typealias FT_Encoding UInt32
const FT_ENCODING_NONE = @compat UInt32(0)
const FT_ENCODING_MS_SYMBOL = @compat UInt32(1937337698)
const FT_ENCODING_UNICODE = @compat UInt32(1970170211)
const FT_ENCODING_SJIS = @compat UInt32(1936353651)
const FT_ENCODING_GB2312 = @compat UInt32(1734484000)
const FT_ENCODING_BIG5 = @compat UInt32(1651074869)
const FT_ENCODING_WANSUNG = @compat UInt32(2002873971)
const FT_ENCODING_JOHAB = @compat UInt32(1785686113)
const FT_ENCODING_MS_SJIS = @compat UInt32(1936353651)
const FT_ENCODING_MS_GB2312 = @compat UInt32(1734484000)
const FT_ENCODING_MS_BIG5 = @compat UInt32(1651074869)
const FT_ENCODING_MS_WANSUNG = @compat UInt32(2002873971)
const FT_ENCODING_MS_JOHAB = @compat UInt32(1785686113)
const FT_ENCODING_ADOBE_STANDARD = @compat UInt32(1094995778)
const FT_ENCODING_ADOBE_EXPERT = @compat UInt32(1094992453)
const FT_ENCODING_ADOBE_CUSTOM = @compat UInt32(1094992451)
const FT_ENCODING_ADOBE_LATIN_1 = @compat UInt32(1818326065)
const FT_ENCODING_OLD_LATIN_2 = @compat UInt32(1818326066)
const FT_ENCODING_APPLE_ROMAN = @compat UInt32(1634889070)
# end enum FT_Encoding
abstract FT_Face_InternalRec_
typealias FT_Face_Internal Ptr{FT_Face_InternalRec_}
abstract FT_Size_InternalRec_
typealias FT_Size_Internal Ptr{FT_Size_InternalRec_}
immutable FT_Size_Metrics
    x_ppem::FT_UShort
    y_ppem::FT_UShort
    x_scale::FT_Fixed
    y_scale::FT_Fixed
    ascender::FT_Pos
    descender::FT_Pos
    height::FT_Pos
    max_advance::FT_Pos
end
abstract FT_SubGlyphRec_
typealias FT_SubGlyph Ptr{FT_SubGlyphRec_}
abstract FT_Slot_InternalRec_
typealias FT_Slot_Internal Ptr{FT_Slot_InternalRec_}
type FT_GlyphSlotRec
    library::FT_Library
    face::Ptr{Void}
    next::Ptr{FT_GlyphSlotRec}
    reserved::FT_UInt
    generic::FT_Generic
    metrics::FT_Glyph_Metrics
    linearHoriAdvance::FT_Fixed
    linearVertAdvance::FT_Fixed
    advance::FT_Vector
    format::FT_Glyph_Format
    bitmap::FT_Bitmap
    bitmap_left::FT_Int
    bitmap_top::FT_Int
    outline::FT_Outline
    num_subglyphs::FT_UInt
    subglyphs::FT_SubGlyph
    control_data::Ptr{Void}
    control_len::Clong
    lsb_delta::FT_Pos
    rsb_delta::FT_Pos
    other::Ptr{Void}
    internal::FT_Slot_Internal
end
typealias FT_GlyphSlot Ptr{FT_GlyphSlotRec}
type FT_FaceRec
    num_faces::FT_Long
    face_index::FT_Long
    face_flags::FT_Long
    style_flags::FT_Long
    num_glyphs::FT_Long
    family_name::Ptr{FT_String}
    style_name::Ptr{FT_String}
    num_fixed_sizes::FT_Int
    available_sizes::Ptr{FT_Bitmap_Size}
    num_charmaps::FT_Int
    charmaps::Ptr{FT_CharMap}
    generic::FT_Generic
    bbox::FT_BBox
    units_per_EM::FT_UShort
    ascender::FT_Short
    descender::FT_Short
    height::FT_Short
    max_advance_width::FT_Short
    max_advance_height::FT_Short
    underline_position::FT_Short
    underline_thickness::FT_Short
    glyph::FT_GlyphSlot
    size::FT_Size
    charmap::FT_CharMap
    driver::FT_Driver
    memory::FT_Memory
    stream::FT_Stream
    sizes_list::FT_ListRec
    autohint::FT_Generic
    extensions::Ptr{Void}
    internal::FT_Face_Internal
end
typealias FT_Face Ptr{FT_FaceRec}
type FT_CharMapRec <: FT_CharMapRec_
    face::FT_Face
    encoding::FT_Encoding
    platform_id::FT_UShort
    encoding_id::FT_UShort
end
type FT_SizeRec <: FT_SizeRec_
    face::FT_Face
    generic::FT_Generic
    metrics::FT_Size_Metrics
    internal::FT_Size_Internal
end
type FT_Parameter
    tag::FT_ULong
    data::FT_Pointer
end
type FT_Open_Args
    flags::FT_UInt
    memory_base::Ptr{FT_Byte}
    memory_size::FT_Long
    pathname::Ptr{FT_String}
    stream::FT_Stream
    driver::FT_Module
    num_params::FT_Int
    params::Ptr{FT_Parameter}
end
# begin enum FT_Size_Request_Type
typealias FT_Size_Request_Type UInt32
const FT_SIZE_REQUEST_TYPE_NOMINAL = @compat UInt32(0)
const FT_SIZE_REQUEST_TYPE_REAL_DIM = @compat UInt32(1)
const FT_SIZE_REQUEST_TYPE_BBOX = @compat UInt32(2)
const FT_SIZE_REQUEST_TYPE_CELL = @compat UInt32(3)
const FT_SIZE_REQUEST_TYPE_SCALES = @compat UInt32(4)
const FT_SIZE_REQUEST_TYPE_MAX = @compat UInt32(5)
# end enum FT_Size_Request_Type
type FT_Size_RequestRec
    _type::FT_Size_Request_Type
    width::FT_Long
    height::FT_Long
    horiResolution::FT_UInt
    vertResolution::FT_UInt
end
typealias FT_Size_Request Ptr{FT_Size_RequestRec}
# begin enum FT_Render_Mode
typealias FT_Render_Mode UInt32
const FT_RENDER_MODE_NORMAL = @compat UInt32(0)
const FT_RENDER_MODE_LIGHT = @compat UInt32(1)
const FT_RENDER_MODE_MONO = @compat UInt32(2)
const FT_RENDER_MODE_LCD = @compat UInt32(3)
const FT_RENDER_MODE_LCD_V = @compat UInt32(4)
const FT_RENDER_MODE_MAX = @compat UInt32(5)
# end enum FT_Render_Mode
# begin enum FT_Kerning_Mode
typealias FT_Kerning_Mode UInt32
const FT_KERNING_DEFAULT = @compat UInt32(0)
const FT_KERNING_UNFITTED = @compat UInt32(1)
const FT_KERNING_UNSCALED = @compat UInt32(2)
# end enum FT_Kerning_Mode
function FT_Init_FreeType(alibrary)
    ccall((:FT_Init_FreeType,freetype),FT_Error,(Ptr{FT_Library},),alibrary)
end
function FT_Done_FreeType(library::FT_Library)
    ccall((:FT_Done_FreeType,freetype),FT_Error,(FT_Library,),library)
end
function FT_New_Face(library::FT_Library,filepathname::AbstractString,face_index::Integer, aface)
    ccall((:FT_New_Face,freetype),FT_Error,(FT_Library,Ptr{UInt8},FT_Long,Ptr{FT_Face}),library,filepathname,face_index,aface)
end
function FT_New_Memory_Face(library::FT_Library,file_base::Ptr{FT_Byte},file_size::FT_Long,face_index::FT_Long,aface::Ptr{FT_Face})
    ccall((:FT_New_Memory_Face,freetype),FT_Error,(FT_Library,Ptr{FT_Byte},FT_Long,FT_Long,Ptr{FT_Face}),library,file_base,file_size,face_index,aface)
end
function FT_Open_Face(library::FT_Library,args::Ptr{FT_Open_Args},face_index::FT_Long,aface::Ptr{FT_Face})
    ccall((:FT_Open_Face,freetype),FT_Error,(FT_Library,Ptr{FT_Open_Args},FT_Long,Ptr{FT_Face}),library,args,face_index,aface)
end
function FT_Attach_File(face::FT_Face,filepathname::AbstractString)
    ccall((:FT_Attach_File,freetype),FT_Error,(FT_Face,Ptr{UInt8}),face,filepathname)
end
function FT_Attach_Stream(face::FT_Face,parameters::Ptr{FT_Open_Args})
    ccall((:FT_Attach_Stream,freetype),FT_Error,(FT_Face,Ptr{FT_Open_Args}),face,parameters)
end
function FT_Reference_Face(face::FT_Face)
    ccall((:FT_Reference_Face,freetype),FT_Error,(FT_Face,),face)
end
function FT_Done_Face(face::FT_Face)
    ccall((:FT_Done_Face,freetype),FT_Error,(FT_Face,),face)
end
function FT_Select_Size(face::FT_Face,strike_index::FT_Int)
    ccall((:FT_Select_Size,freetype),FT_Error,(FT_Face,FT_Int),face,strike_index)
end
function FT_Request_Size(face::FT_Face,req::FT_Size_Request)
    ccall((:FT_Request_Size,freetype),FT_Error,(FT_Face,FT_Size_Request),face,req)
end
function FT_Set_Char_Size(face::FT_Face, char_width, char_height, horz_resolution, vert_resolution)
    ccall((:FT_Set_Char_Size,freetype),FT_Error,(FT_Face,FT_F26Dot6,FT_F26Dot6,FT_UInt,FT_UInt),face,char_width,char_height,horz_resolution,vert_resolution)
end
function FT_Set_Pixel_Sizes(face::FT_Face,pixel_width::FT_UInt,pixel_height::FT_UInt)
    ccall((:FT_Set_Pixel_Sizes,freetype),FT_Error,(FT_Face,FT_UInt,FT_UInt),face,pixel_width,pixel_height)
end
function FT_Load_Glyph(face::FT_Face,glyph_index::FT_UInt,load_flags::FT_Int32)
    ccall((:FT_Load_Glyph,freetype),FT_Error,(FT_Face,FT_UInt,FT_Int32),face,glyph_index,load_flags)
end
function FT_Load_Char(face::FT_Face,char_code::Char,load_flags::FT_Int32)
    ccall((:FT_Load_Char,freetype),FT_Error,(FT_Face,FT_ULong,FT_Int32),face,char_code,load_flags)
end
function FT_Set_Transform(face::FT_Face,matrix::Ptr{FT_Matrix},delta::Ptr{FT_Vector})
    ccall((:FT_Set_Transform,freetype),Void,(FT_Face,Ptr{FT_Matrix},Ptr{FT_Vector}),face,matrix,delta)
end
function FT_Render_Glyph(slot::FT_GlyphSlot,render_mode::FT_Render_Mode)
    ccall((:FT_Render_Glyph,freetype),FT_Error,(FT_GlyphSlot,FT_Render_Mode),slot,render_mode)
end
function FT_Get_Kerning(face::FT_Face,left_glyph::FT_UInt,right_glyph::FT_UInt,kern_mode::FT_UInt, akerning)
    ccall((:FT_Get_Kerning,freetype),FT_Error,(FT_Face,FT_UInt,FT_UInt,FT_UInt,Ptr{FT_Vector}),face,left_glyph,right_glyph,kern_mode,akerning)
end
function FT_Get_Track_Kerning(face::FT_Face,point_size::FT_Fixed,degree::FT_Int,akerning::Ptr{FT_Fixed})
    ccall((:FT_Get_Track_Kerning,freetype),FT_Error,(FT_Face,FT_Fixed,FT_Int,Ptr{FT_Fixed}),face,point_size,degree,akerning)
end
function FT_Get_Glyph_Name(face::FT_Face,glyph_index::FT_UInt,buffer::FT_Pointer,buffer_max::FT_UInt)
    ccall((:FT_Get_Glyph_Name,freetype),FT_Error,(FT_Face,FT_UInt,FT_Pointer,FT_UInt),face,glyph_index,buffer,buffer_max)
end
function FT_Get_Postscript_Name(face::FT_Face)
    ccall((:FT_Get_Postscript_Name,freetype),Ptr{UInt8},(FT_Face,),face)
end
function FT_Select_Charmap(face::FT_Face,encoding::FT_Encoding)
    ccall((:FT_Select_Charmap,freetype),FT_Error,(FT_Face,FT_Encoding),face,encoding)
end
function FT_Set_Charmap(face::FT_Face,charmap::FT_CharMap)
    ccall((:FT_Set_Charmap,freetype),FT_Error,(FT_Face,FT_CharMap),face,charmap)
end
function FT_Get_Charmap_Index(charmap::FT_CharMap)
    ccall((:FT_Get_Charmap_Index,freetype),FT_Int,(FT_CharMap,),charmap)
end
function FT_Get_Char_Index(face::FT_Face,charcode::Char)
    ccall((:FT_Get_Char_Index,freetype),FT_UInt,(FT_Face,FT_ULong),face,charcode)
end
function FT_Get_First_Char(face::FT_Face,agindex::Ptr{FT_UInt})
    ccall((:FT_Get_First_Char,freetype),FT_ULong,(FT_Face,Ptr{FT_UInt}),face,agindex)
end
function FT_Get_Next_Char(face::FT_Face,char_code::FT_ULong,agindex::Ptr{FT_UInt})
    ccall((:FT_Get_Next_Char,freetype),FT_ULong,(FT_Face,FT_ULong,Ptr{FT_UInt}),face,char_code,agindex)
end
function FT_Get_Name_Index(face::FT_Face,glyph_name::Ptr{FT_String})
    ccall((:FT_Get_Name_Index,freetype),FT_UInt,(FT_Face,Ptr{FT_String}),face,glyph_name)
end
function FT_Get_SubGlyph_Info(glyph::FT_GlyphSlot,sub_index::FT_UInt,p_index::Ptr{FT_Int},p_flags::Ptr{FT_UInt},p_arg1::Ptr{FT_Int},p_arg2::Ptr{FT_Int},p_transform::Ptr{FT_Matrix})
    ccall((:FT_Get_SubGlyph_Info,freetype),FT_Error,(FT_GlyphSlot,FT_UInt,Ptr{FT_Int},Ptr{FT_UInt},Ptr{FT_Int},Ptr{FT_Int},Ptr{FT_Matrix}),glyph,sub_index,p_index,p_flags,p_arg1,p_arg2,p_transform)
end
function FT_Get_FSType_Flags(face::FT_Face)
    ccall((:FT_Get_FSType_Flags,freetype),FT_UShort,(FT_Face,),face)
end
function FT_Face_GetCharVariantIndex(face::FT_Face,charcode::FT_ULong,variantSelector::FT_ULong)
    ccall((:FT_Face_GetCharVariantIndex,freetype),FT_UInt,(FT_Face,FT_ULong,FT_ULong),face,charcode,variantSelector)
end
function FT_Face_GetCharVariantIsDefault(face::FT_Face,charcode::FT_ULong,variantSelector::FT_ULong)
    ccall((:FT_Face_GetCharVariantIsDefault,freetype),FT_Int,(FT_Face,FT_ULong,FT_ULong),face,charcode,variantSelector)
end
function FT_Face_GetVariantSelectors(face::FT_Face)
    ccall((:FT_Face_GetVariantSelectors,freetype),Ptr{FT_UInt32},(FT_Face,),face)
end
function FT_Face_GetVariantsOfChar(face::FT_Face,charcode::FT_ULong)
    ccall((:FT_Face_GetVariantsOfChar,freetype),Ptr{FT_UInt32},(FT_Face,FT_ULong),face,charcode)
end
function FT_Face_GetCharsOfVariant(face::FT_Face,variantSelector::FT_ULong)
    ccall((:FT_Face_GetCharsOfVariant,freetype),Ptr{FT_UInt32},(FT_Face,FT_ULong),face,variantSelector)
end
function FT_MulDiv(a::FT_Long,b::FT_Long,c::FT_Long)
    ccall((:FT_MulDiv,freetype),FT_Long,(FT_Long,FT_Long,FT_Long),a,b,c)
end
function FT_DivFix(a::FT_Long,b::FT_Long)
    ccall((:FT_DivFix,freetype),FT_Long,(FT_Long,FT_Long),a,b)
end
function FT_RoundFix(a::FT_Fixed)
    ccall((:FT_RoundFix,freetype),FT_Fixed,(FT_Fixed,),a)
end
function FT_CeilFix(a::FT_Fixed)
    ccall((:FT_CeilFix,freetype),FT_Fixed,(FT_Fixed,),a)
end
function FT_FloorFix(a::FT_Fixed)
    ccall((:FT_FloorFix,freetype),FT_Fixed,(FT_Fixed,),a)
end
function FT_Vector_Transform(vec::Ptr{FT_Vector},matrix::Ptr{FT_Matrix})
    ccall((:FT_Vector_Transform,freetype),Void,(Ptr{FT_Vector},Ptr{FT_Matrix}),vec,matrix)
end
function FT_Library_Version(library::FT_Library,amajor::Ptr{FT_Int},aminor::Ptr{FT_Int},apatch::Ptr{FT_Int})
    ccall((:FT_Library_Version,freetype),Void,(FT_Library,Ptr{FT_Int},Ptr{FT_Int},Ptr{FT_Int}),library,amajor,aminor,apatch)
end
function FT_Face_CheckTrueTypePatents(face::FT_Face)
    ccall((:FT_Face_CheckTrueTypePatents,freetype),FT_Bool,(FT_Face,),face)
end
function FT_Face_SetUnpatentedHinting(face::FT_Face,value::FT_Bool)
    ccall((:FT_Face_SetUnpatentedHinting,freetype),FT_Bool,(FT_Face,FT_Bool),face,value)
end

export FT_Init_FreeType
export FT_Done_FreeType
export FT_New_Face
export FT_New_Memory_Face
export FT_Open_Face
export FT_Attach_File
export FT_Attach_Stream
export FT_Reference_Face
export FT_Done_Face
export FT_Select_Size
export FT_Request_Size
export FT_Set_Char_Size
export FT_Set_Pixel_Sizes
export FT_Load_Glyph
export FT_Load_Char
export FT_Set_Transform
export FT_Render_Glyph
export FT_Get_Kerning
export FT_Get_Track_Kerning
export FT_Get_Glyph_Name
export FT_Get_Postscript_Name
export FT_Select_Charmap
export FT_Set_Charmap
export FT_Get_Charmap_Index
export FT_Get_Char_Index
export FT_Get_First_Char
export FT_Get_Next_Char
export FT_Get_Name_Index
export FT_Get_SubGlyph_Info
export FT_Get_FSType_Flags
export FT_Face_GetCharVariantIndex
export FT_Face_GetCharVariantIsDefault
export FT_Face_GetVariantSelectors
export FT_Face_GetVariantsOfChar
export FT_Face_GetCharsOfVariant
export FT_MulDiv
export FT_DivFix
export FT_RoundFix
export FT_CeilFix
export FT_FloorFix
export FT_Vector_Transform
export FT_Library_Version
export FT_Face_CheckTrueTypePatents
export FT_Face_SetUnpatentedHinting

export FT_Int16
export FT_UInt16
export FT_Int32
export FT_UInt32
export FT_Fast
export FT_UFast
export FT_Memory
export FT_Alloc_Func
export FT_Free_Func
export FT_Realloc_Func
export FT_Stream
export FT_Stream_IoFunc
export FT_Stream_CloseFunc
export FT_Pos
export FT_Pixel_Mode
export FT_Outline_MoveToFunc
export FT_Outline_LineToFunc
export FT_Outline_ConicToFunc
export FT_Outline_CubicToFunc
export FT_Glyph_Format
export FT_Raster
export FT_SpanFunc
export FT_Raster_BitTest_Func
export FT_Raster_BitSet_Func
export FT_Raster_NewFunc
export FT_Raster_DoneFunc
export FT_Raster_ResetFunc
export FT_Raster_SetModeFunc
export FT_Raster_RenderFunc
export FT_Bool
export FT_FWord
export FT_UFWord
export FT_Char
export FT_Byte
export FT_Bytes
export FT_Tag
export FT_String
export FT_Short
export FT_UShort
export FT_Int
export FT_UInt
export FT_Long
export FT_ULong
export FT_F2Dot14
export FT_F26Dot6
export FT_Fixed
export FT_Error
export FT_Pointer
export FT_Offset
export FT_PtrDist
export FT_Generic_Finalizer
export FT_ListNode
export FT_List
export FT_Library
export FT_Module
export FT_Driver
export FT_Renderer
export FT_Face
export FT_Size
export FT_GlyphSlot
export FT_CharMap
export FT_Encoding
export FT_Face_Internal
export FT_Size_Internal
export FT_SubGlyph
export FT_Slot_Internal
export FT_Size_Request_Type
export FT_Size_Request
export FT_Render_Mode
export FT_Kerning_Mode

export FT_LOAD_DEFAULT
export FT_LOAD_TARGET_NORMAL
export FT_LOAD_NO_SCALE
export FT_LOAD_NO_HINTING
export FT_LOAD_RENDER
export FT_LOAD_NO_BITMAP
export FT_LOAD_VERTICAL_LAYOUT
export FT_LOAD_FORCE_AUTOHINT
export FT_LOAD_CROP_BITMAP
export FT_LOAD_PEDANTIC
export FT_LOAD_ADVANCE_ONLY
export FT_LOAD_IGNORE_GLOBAL_ADVANCE_WIDTH
export FT_LOAD_NO_RECURSE
export FT_LOAD_IGNORE_TRANSFORM
export FT_LOAD_MONOCHROME
export FT_LOAD_LINEAR_DESIGN
export FT_LOAD_SBITS_ONLY
export FT_LOAD_NO_AUTOHINT
export FT_LOAD_TARGET_LIGHT
export FT_LOAD_TARGET_MONO
export FT_LOAD_TARGET_LCD
export FT_LOAD_TARGET_LCD_V
export FT_LOAD_COLOR

export FT_FACE_FLAG_SCALABLE
export FT_FACE_FLAG_FIXED_SIZES
export FT_FACE_FLAG_FIXED_WIDTH
export FT_FACE_FLAG_SFNT
export FT_FACE_FLAG_HORIZONTAL
export FT_FACE_FLAG_VERTICAL
export FT_FACE_FLAG_KERNING
export FT_FACE_FLAG_FAST_GLYPHS
export FT_FACE_FLAG_MULTIPLE_MASTERS
export FT_FACE_FLAG_GLYPH_NAMES
export FT_FACE_FLAG_EXTERNAL_STREAM
export FT_FACE_FLAG_HINTER
export FT_FACE_FLAG_CID_KEYED
export FT_FACE_FLAG_TRICKY
export FT_FACE_FLAG_COLOR

export FT_RENDER_MODE_NORMAL
export FT_RENDER_MODE_LIGHT
export FT_RENDER_MODE_MONO
export FT_RENDER_MODE_LCD
export FT_RENDER_MODE_LCD_V
export FT_RENDER_MODE_MAX

export FT_PIXEL_MODE_NONE
export FT_PIXEL_MODE_MONO
export FT_PIXEL_MODE_GRAY
export FT_PIXEL_MODE_GRAY2
export FT_PIXEL_MODE_GRAY4
export FT_PIXEL_MODE_LCD
export FT_PIXEL_MODE_LCD_V
export FT_PIXEL_MODE_BGRA
export FT_PIXEL_MODE_MAX

end


text:
"""Load from VTU
Obviously we want to load the line mesh. Meshes of triangles/tetrahedrons
are loaded in as a graph of edges.
"""
function read_curve_vtu(path::AbstractString)
    last(splitext(path)) == ".vtu" || throw("$(path) is not a VTU file")
    # The hierchy is doc
    #                  piece
    #                    points
    #                    cells
    xroot = root(parse_file(path))
    grid = first(child_elements(xroot))
    piece = first(child_elements(grid))

    npoints = parse(Int, attribute(piece, "NumberOfPoints"))

    points, cells = child_elements(piece)
    # The data is X0 Y0 Z0  X1 Y1 Z1
    points_data = split(strip(content(points)), "  ")
    # Identify dimensionality from offets
    segments, offsets, types = child_elements(cells)

    # At this point we have one of line/tri/tet mesh. Line mesh is easiest
    # as we already have all the edges for tri and tet they need to be
    # made manually
    mesh_types = Set(split(strip(content(types)), " "))
    @assert length(mesh_types) == 1
    mesh_type = pop!(mesh_types)

    # A line mesh
    if mesh_type == "3"
        dim = parse(Int, string(first(content(offsets))))
        # Index data is i0 j0  i1 j1
        segments_data = split(strip(content(segments)), "  ")

        nsegments = parse(Int, attribute(piece, "NumberOfCells"))
p        # Finally read in
        segments = Vector{Tuple{Int, Int}}(nsegments)
        for (i, seg) in enumerate(segments_data)
            # One based indexing
            segments[i] = Tuple((1+parse(Int, v)) for v in split(seg))
        end

        points = Vector{SVector{dim, Float64}}(npoints)
        for (i, point) in enumerate(points_data)
            points[i] = SVector((parse(Float64, v) for v in split(point)[1:dim])...)
        end

        return (points, segments)
    end

    # A simplex mesh of triangles or tets
    @assert mesh_type == "5" || mesh_type == "10"

    nvertices = (mesh_type == "5") ? 3 : 4
    # Cell local
    edge_indices = collect(Combinatorics.combinations(collect(1:nvertices), 2))
    # Cells are made of indicices - need to split to edges
    cell_data = split(strip(content(segments)), "  ")

    # Ath this point we don't really know about the segments counts
    segments = Set{Tuple{Int, Int}}()
    for cell in imap(s -> map(x -> parse(Int, x), split(s, " ")), cell_data)
        # Break cell to edges
        for (i0, i1) in edge_indices
            # One based indexing
            v0, v1 = cell[i0]+1, cell[i1]+1
            push!(segments, (v0 < v1) ? (v0, v1) : (v1, v0))
        end
    end
    # Unique can be converted
    segments = collect(segments)

    # We load points assuming 3d
    points = Vector{SVector{3, Float64}}(npoints)
    for (i, point) in enumerate(points_data)
        points[i] = SVector((parse(Float64, v) for v in split(point))...)
    end

    # For tet the 3d assumption is correct and we're done
    nvertices == 4 && return (points, segments)
    # For triangle if it is in 3d (not z == 0 plane) 3d was right
    norm(map(last, points), 1) > 0 && return (points, segments)
    # Otherwise we project to 2d
    points2d = Vector{SVector{2, Float64}}(npoints)
    for (row, point) in enumerate(points)
        points2d[row] = SVector(point[1], point[2])
    end

    (points2d, segments)
end

"""Write curve function"""
function write(path::AbstractString, f::CurveFunction)
    name, ext = splitext(path)
    ext == ".vtu" || throw("$(path) is not a VTU file")

    curve = f.curve
    cells = Vector{MeshCell}(length(curve.segments))
    for (i, seg) in enumerate(curve.segments)
        cells[i] = MeshCell(VTKCellTypes.VTK_LINE, [seg[1], seg[2]])
    end

    points = zeros((geometrical_dim(f), length(curve.points)))
    for (col, point) in enumerate(curve.points)
        for (row, value) in enumerate(point)
            points[row, col] = value
        end
    end

    vtk_file = vtk_grid(name, points, cells)
    topological_dim(f) == 0 && vtk_point_data(vtk_file, f.data, "point_values")
    topological_dim(f) == 1 && vtk_cell_data(vtk_file, f.data, "edge_values")

    vtk_save(vtk_file)
end

"""Write curve"""
function write{D, T}(path::AbstractString, curve::Curve{D, T})
    name, ext = splitext(path)
    ext == ".vtu" || throw("$(path) is not a VTU file")

    cells = Vector{MeshCell}(length(curve.segments))
    for (i, seg) in enumerate(curve.segments)
        cells[i] = MeshCell(VTKCellTypes.VTK_LINE, [seg[1], seg[2]])
    end

    points = zeros((D, length(curve.points)))
    for (col, point) in enumerate(curve.points)
        for (row, value) in enumerate(point)
            points[row, col] = value
        end
    end

    vtk_file = vtk_grid(name, points, cells)
    vtk_save(vtk_file)
end


text:
"""
create (n x n) identity matrix

# Arguments
`T`: type of the matrix elements

`n`: dimensions
"""
function eye(T, n::Int)
  return Matrix{T}(LinearAlgebra.I, n, n)
end

"""
left order a binary matrix (Z)

convenience function for:
```julia
julia> sortslices(Z, dims=2, rev=true)
```
"""
function leftOrder(Z::Matrix{T}) where T
  return sortslices(Z, dims=2, rev=true)
end

"""
Generate a simple Z matrix. For simulation
studies. J, K are integers and dimensions
of the (binary) Z matrix. J must be
a multiple of K.
"""
function genSimpleZ(J::Int, K::Int)
  g = div(J, K)
  @assert g == J / K

  return kron(eye(Int, K), fill(1,g))
end

"""
Assert that a Z matrix is valid
"""
function isValidZ(Z)
  hasRepeatedColumns = size(unique(Z, dims=2), 2) < size(Z, 2)
  hasColOfZero =  any(sum(Z, dims=1) .== 0)
  return  !(hasRepeatedColumns || hasRepeatedColumns)
end

"""
Generate a (J x K) Z matrix for simulation studies.

`prob1` is the desired proportion of 1's in the matrix.
"""
function genZ(J::Int, K::Int, prob1::Float64)::Matrix{Int}
  @assert 0 < prob1 < 1
  Z = Int.(rand(J, K) .> prob1)
  Z = sortslices(Z, dims=1, rev=true)
  Z = leftOrder(Z)

  # if size(unique(Z, dims=2), 2) < K || any(sum(Z, dims=1) .== 0) || all(sum(Z, dims=2) .== 0)
  if isValidZ(Z)
    return Z
  else
    return genZ(J, K, prob1)
  end
end

#= Test
Z = Int.(randn(3,5) .> 0)
=#

function genData(J::Int, N::Vector{Int}, K::Int, L::Dict{Int, Int};
                 useSimpleZ::Bool=true, prob1::Float64=.6,
                 sortLambda::Bool=false, propMissingScale::Float64=0.7)

  I = length(N)
  Z = useSimpleZ ? genSimpleZ(J, K) : genZ(J, K, prob1)

  genData(J=J, N=N, K=K, L=L, Z=Z,
          beta=[-9.2, -2.3], # linear missing mechanism
          sig2=fill(0.1, I),
          mus=Dict(0=>collect(range(-5, length=L[0], stop=-1)),
                   1=>collect(range(1, length=L[1], stop=5))),
          a_W=[float(i) for i in 1:K],
          a_eta=Dict([ z => [float(l) for l in 1:L[z]] for z in 0:1 ]),
          sortLambda=sortLambda, propMissingScale=propMissingScale)
end

function genData(; J::Int, N::Vector{Int}, K::Int, L::Dict{Int, Int},
                 Z::Matrix{Int}, beta::Vector{Float64},
                 sig2::Vector{Float64}, mus::Dict{Int, Vector{Float64}},
                 a_W::Vector{Float64}, a_eta::Dict{Int, Vector{Float64}},
                 sortLambda::Bool=false, propMissingScale::Float64=0.7,
                 eps=zeros(length(N)))

  # Check eps
  @assert length(eps) == length(N)
  @assert all(0 .<= eps .<= 1)

  # Check Z dimensions
  @assert ncol(Z) == K && nrow(Z) == J

  # Check N dimensions
  @assert all(N .> 0)

  I = length(N)

  # Check sig2 dimensions
  @assert all(sig2 .> 0)
  @assert length(sig2) == I

  # Check mus dimensions
  @assert all(mus[0] .< 0) && all(mus[1] .> 0)
  @assert length(mus[0]) == L[0] && length(mus[1]) == L[1]

  # Check a_W dimensions
  @assert length(a_W) == K
  @assert all(a_W .> 0)

  # Check a_eta dimensions
  @assert length(a_eta[0]) == L[0]
  @assert length(a_eta[1]) == L[1]
  @assert all(a_eta[0] .> 0)
  @assert all(a_eta[1] .> 0)

  # Simulate W
  W = zeros(I, K)
  for i in 1:I
    W[i,:] = rand(Dirichlet(Random.shuffle(a_W)))
  end

  # Simulate lambda
  lam = [begin
           p = [eps[i]; (1 - eps[i]) * W[i, :]]
           rand(Categorical(p), N[i]) .- 1
         end for i in 1:I]
  if sortLambda
    lam = [sort(lami) for lami in lam]
  end

  # Simulate eta
  eta = Dict([ z => zeros(Float64, I, J, L[z]) for z in 0:1 ])
  for i in 1:I
    for j in 1:J
      eta[0][i, j, :] = rand(Dirichlet(Random.shuffle(a_eta[0])))
      eta[1][i, j, :] = rand(Dirichlet(Random.shuffle(a_eta[1])))
    end
  end

  # Simulate gam
  gam = [zeros(Int, Ni, J) for Ni in N]
  for i in 1:I
    for j in 1:J
      for n in 1:N[i]
        k = lam[i][n]
        if k > 0
          z = Z[j, k]
          gam[i][n,j] = rand(Categorical(eta[z][i, j, :]))
        else
          gam[i][n,j] = 0
        end
      end
    end
  end

  # Generate y
  #y = [zeros(Union{Float64, Missing}, Ni, J) for Ni in N]
  y = [zeros(Float64, Ni, J) for Ni in N]
  y_complete = [zeros(Float64, Ni, J) for Ni in N]

  function z_get(i, n, j)
    return Z[j, lam[i][n]]
  end

  function mu_get(i, n, j)
    l = gam[i][n, j]
    if l > 0
      z = z_get(i, n, j)
      out = mus[z][l]
    else
      out = 0.0
    end
    return out
  end

  for i in 1:I
    for j in 1:J
      for n in 1:N[i]
        sig_i = lam[i][n] > 0 ? sqrt(sig2[i]) : 3.0
        y_complete[i][n, j] = rand(Normal(mu_get(i, n, j), sig_i))
      end

      # Set some to be missing
      p_miss = [prob_miss(y_complete[i][n, j], beta) for n in 1:N[i]]
      prop_missing = rand() * propMissingScale * sum(W[i,:] .* (1 .- Z[j,:]))
      num_missing = Int(round(N[i] * prop_missing))
      idx_missing = Distributions.wsample(1:N[i], p_miss, num_missing, replace=false)
      y[i][:, j] .= y_complete[i][:, j] .+ 0
      y[i][idx_missing, j] .= NaN
    end
  end

  return Dict(:y=>y, :y_complete=>y_complete, :Z=>Z, :W=>W,
              :eta=>eta, :mus=>mus, :sig2=>sig2, :lam=>lam, :gam=>gam,
              :beta=>beta, :eps=>eps)
end # genData

#precompile(genData, (String, Int, Vector{Int}, Int, Int, Bool, Float64, Bool, Float64))


text:
"""
    fullonelevel

Compute the mesh topology structures that correspond to the full one-level
downward and upward adjacency graph. Memory is used to store the downward incidence
relations (3, 2), (2, 1), (1, 0), and the upward incidence relations (0, 1), (1,
2), and (2, 3) and the locations of the vertices.
"""
module fullonelevel
using StaticArrays
using MeshCore: nvertices, nshapes, manifdim, attribute
using MeshCore: ir_skeleton, ir_bbyfacets, ir_bbyridges, ir_transpose
using MeshSteward: vtkwrite, T4block
using BenchmarkTools
using Test

include("usedbytes.jl")

function test()
    ts = time()
    n = 3
    membytes = 0; summembytes = 0
    @info "Initial (3, 0)"
    @time connectivity = T4block(1.0, 2.0, 3.0, n*7, n*9, n*10, :a; intbytes = 4)
    ir30 = connectivity
    @show "($(manifdim(ir30.left)), $(manifdim(ir30.right)))"
    @show (nshapes(ir30.left), nshapes(ir30.right))
    @show membytes = usedbytes(ir30._v)
    geom = attribute(ir30.right, "geom")
    @show membytes = usedbytes(geom.v)
    summembytes += membytes
    @show time() - ts

    ts = time()
    @info "Skeleton: facets. (2, 0)"
    @time ir20 = ir_skeleton(ir30)
    @show "($(manifdim(ir20.left)), $(manifdim(ir20.right)))"
    @show (nshapes(ir20.left), nshapes(ir20.right))
    @show membytes = usedbytes(ir20._v)
    summembytes += membytes

    @info "Bounded-by facets. (3, 2)"
    @time ir32 = ir_bbyfacets(ir30, ir20)
    @show "($(manifdim(ir32.left)), $(manifdim(ir32.right)))"
    @show (nshapes(ir32.left), nshapes(ir32.right))
    @show membytes = usedbytes(ir32._v)
    summembytes += membytes

    @info "Skeleton: ridges. (1, 0)"
    @time ir10 = ir_skeleton(ir20)
    @show "($(manifdim(ir10.left)), $(manifdim(ir10.right)))"
    @show (nshapes(ir10.left), nshapes(ir10.right))
    @show membytes = usedbytes(ir10._v)
    summembytes += membytes

    @info "Bounded-by facets. (2, 1)"
    @time ir21 = ir_bbyfacets(ir20, ir10)
    @show "($(manifdim(ir21.left)), $(manifdim(ir21.right)))"
    @show (nshapes(ir21.left), nshapes(ir21.right))
    @show membytes = usedbytes(ir21._v)
    summembytes += membytes

    @info "Transpose. (0, 1)"
    @time tr01 = ir_transpose(ir10)
    @show "($(manifdim(tr01.left)), $(manifdim(tr01.right)))"
    @show (nshapes(tr01.left), nshapes(tr01.right))
    @show membytes = usedbytes(tr01._v)
    summembytes += membytes

    @info "Transpose. (1, 2)"
    @time tr12 = ir_transpose(ir21)
    @show "($(manifdim(tr12.left)), $(manifdim(tr12.right)))"
    @show (nshapes(tr12.left), nshapes(tr12.right))
    @show membytes = usedbytes(tr12._v)
    summembytes += membytes

    @info "Transpose. (2, 3)"
    @time tr23 = ir_transpose(ir32)
    @show "($(manifdim(tr23.left)), $(manifdim(tr23.right)))"
    @show (nshapes(tr23.left), nshapes(tr23.right))
    @show membytes = usedbytes(tr23._v)
    summembytes += membytes

    # Print the total number of megabytes used to store the database
    @show summembytes/2^20
    @show time() - ts

    # vtkwrite("speedtest1", connectivity)
    true
end
end
using .fullonelevel
# using BenchmarkTools
# @btime fullonelevel.test()
fullonelevel.test()


text:
module TDD
include("graph.jl")
include("interval.jl")
end


text:
module DocMaker
import Reproducible, Random

Random.seed!(1234)

function Reproducible.render(o::Reproducible.CodeBlock, r::Val{:thing}; kw...)
    """
    ```
    The first line of your code was $(strip((o.rows[1].input)))
    ```
    ```
    The output of that line was $(o.rows[1].repldisplay)
    ```
    """
end

Reproducible.build(joinpath(@__DIR__(), "index.md"); toc=true)
end


text:
"""
    RatingCurve(G::Vector{Gauging}, a::Real, b::Real, c::Real)

Construct an object of type RatingCurve
"""
struct RatingCurve
    gauging::Vector{Gauging}
    a::Real
    b::Real
    c::Real
    RatingCurve(gauging::Vector{Gauging}, a::T,b::T,c::T) where {T<:Real} = new(gauging,a,b,c)
end

RatingCurve(gauging::Vector{Gauging}, a::Real, b::Real, c::Real) = RatingCurve(gauging, promote(a, b, c)...)
RatingCurve(gauging::Vector{Gauging}, a::Integer, b::Integer, c::Integer) = RatingCurve(gauging, float(a), float(b), float(c))
RatingCurve() = RatingCurve(Gauging[], 0.0, 0.0, 0.0)

Base.Broadcast.broadcastable(obj::RatingCurve) = Ref(obj)

"""
    function bic(rc::RatingCurve)

BIC of the rating curve model.
"""
function bic(rc::RatingCurve)

    n = length(rc.gauging)
    p = 3

    @assert n>3

    σ̂² = RatingCurves.var(rc)

    return n*log(σ̂²[]) + p*log(n)

end

"""
    cint(rc::RatingCurve; nboot::Int=1000, α::Real=.05)

Rating curve parameter confidence intervals of level `1-α` obtained by a bootstrap sample of size `nboot`.
"""
function cint(rc::RatingCurve; nboot::Int=100, α::Real=.05)

    a = Vector{Float64}(undef, nboot)
    b = Vector{Float64}(undef, nboot)
    c = Vector{Float64}(undef, nboot)

    for i in 1:nboot

        Gᵢ = RatingCurves.bootstrap(rc.gauging)
        rcᵢ = rcfit(Gᵢ)

        a[i] = rcᵢ.a
        b[i] = rcᵢ.b
        c[i] = rcᵢ.c

    end

    M = hcat(
            quantile(a, [α/2, 1-α/2]),
            quantile(b, [α/2, 1-α/2]),
            quantile(c, [α/2, 1-α/2])
            )

    return M

end

"""
    discharge(rc::RatingCurve, h::Real)

Compute the estimated discharge at level `h` with the rating curve `rc`.
"""
function discharge(rc::RatingCurve, h::Real)

    y = logdischarge(rc, h)

    return exp(y)

end

"""
    level(rc::RatingCurve, q::Real)

Estimate the level corresponding to the discharge `q` with the rating curve `rc`.
"""
function level(rc::RatingCurve, q::Real)

    @assert q>0

    h = rc.b + (q/rc.a)^(1/rc.c)

    return h

end

"""
    logdischarge(rc::RatingCurve, h::Real)

Compute the log of the estimated discharge at level `h` with the rating curve `rc`.
"""
function logdischarge(rc::RatingCurve, h::Real)

    if (h>rc.b) & (rc.a>0)
        y = log(rc.a) + rc.c*log(h-rc.b)
    else
        y = -Inf
    end

    return y

end

@doc raw"""
    pint(rc::RatingCurve, level::Real, α::Real=0.05, rtol::Real=.05)

`1-α` confidence interval of the estimated discharge at level `h` with the rating curve `rc`.

See also [`pintlog`](@ref)

## Details

Estimate the `1-α` confidence interval of the estimated discharge corresponding to the rating curve `rc`at level `h`
with the relative error of the discharge measurement `rtol`.

### Measurement error

It is assumed that ``(q_i \pm \operatorname{rtol} \times q_i)`` contains the true discharge 95% of the time. Assuming a Gaussian distribution for the
discharge measurement, the discharge ``Q_i`` distribution is given as follows:

`` Q_i \sim \mathcal{N} \left\{ q_i,  \left(\frac{\operatorname{rtol} q_i}{1.96}\right)^2 \right\}.``

The discharge standard variation is therefore ``\tau_i = \frac{\operatorname{rtol} q_i}{1.96}.``

### Residual error

The estimation error ``σ²`` corresponding to the point on the rating curve is estimated with the sum of squared residuals in the log space.
Let ``qᵢ`` be the observed discharge and ``q̂ᵢ`` be the corresponding discharge estimation. For the ``n>3`` gaugings, the
sum of squared residuals is defined as:

``\operatorname{SSE} = \sum_{i=1}^n ( \log qᵢ - \log q̂ᵢ)^2.``

The variance estimation of the log residuals are therefore:

``σ̂ₑ² = \frac{1}{n-3} \operatorname{SSE}.``

The linear regression prediction variance of the estimated discharge at level `h₀` is given by:

``σ̂² = x₀^⊤ (X^\top X)^{-1} x₀``

where ``x₀^⊤ = [1 \log (h₀ - b)]`` and $X$ is the structure matrix of the regression model.

!!! note
    The uncertainty on $b$ is neglected in this expression.

### Log discharge prediction error

The log prediction variance is assumed to be the sum of the estimation variance and the discharge measurement variance in the log space.
It is then assumed that the prediction error in the log space is Gaussian in order to compute the confidence interval bounds.

### Discharge prediction error

The confidence interval bounds in the original space are given by the exponential of the bounds in the log space.

"""
function pint(rc::RatingCurve, h₀::Real, α::Real=0.05, rtol::Real=.05)

    res = pintlog(rc, h₀, α, rtol)

    return exp.(res)

end

"""
    pintlog(rc::RatingCurve, h₀::Real, α::Real=0.05, rtol::Real=.05)

`1-α` confidence interval of the estimated log discharge at level `h₀` with the rating curve `rc`.

### Details

`rtol` represents the relative uncertainty of the dishcarge so that the true discharge is included in the interval `q ± 1.96*rtol` 95% of the time
"""
function pintlog(rc::RatingCurve, h₀::Real, α::Real=0.05, rtol::Real=.05)

    @assert 0<α<1
    @assert 0<rtol<1
    @assert h₀>rc.b

    # Residuals variance in the log space
    σ̂ₑ² = RatingCurves.var(rc)[]

    # Prediction variance given by the linear regression (uncertainty on b is not taken into account)
    G = rc.gauging
    x = log.(level.(rc.gauging) .- rc.b)
    X = hcat(ones(length(x)), x)
    x₀ = [1, log(h₀ - rc.b)]
    σ̂² = σ̂ₑ² * (1+x₀'/(X'X)*x₀)

    # Relative discharge variance in the log space
    τ² = (rtol/1.96)^2

    pd = Normal(logdischarge(rc, h₀), sqrt(σ̂² + τ²))

    lower = quantile(pd, α/2)
    upper = quantile(pd, 1-α/2)

    return [lower, upper]

end

"""
    Base.show(io::IO, obj::EVA)
Override of the show function for the objects of type RatingCurve.
"""
function Base.show(io::IO, obj::RatingCurve)

    showRatingCurve(io, obj)

end

function showRatingCurve(io::IO, obj::RatingCurve, prefix::String = "")

    println(io, prefix, "RatingCurve")
    println(io, prefix ,"   ", "data: ", typeof(obj.gauging), "[", length(obj.gauging), "]")
    println(io, prefix, "   ", "Parameters:")
    println(io, prefix, "      ", "a = ", obj.a)
    println(io, prefix, "      ", "b = ", obj.b)
    println(io, prefix, "      ", "c = ", obj.c)

end

"""
    sse(rc::RatingCurve)

Compute the sum of the squares of residuals between the curve and the gauging in the log space.
"""
function sse(rc::RatingCurve)

    h = level.(rc.gauging)
    q = discharge.(rc.gauging)

    y = log.(q)

    ŷ = logdischarge.(rc, h)

    SSE = sum((y - ŷ).^2)

    return [SSE]

end

"""
    var(rc::RatingCurve)

Estimate the variance of the errors in the log space.
"""
function var(rc::RatingCurve)

    n = length(rc.gauging)
    p = 3

    SSE = RatingCurves.sse(rc)

    σ̂² = SSE./(n-p)

    return σ̂²

end


text:
using Plots
plotly()

function tsplot(model::RMB.InventoryModel,ka)
    T=size(ka,1)
    dd,dc,ds,dy,dr=model.dims
    ddd=[ind2sub((dd,dc,ds),ka[t,1])[1]-1 for t=1:T]
    cc=[model.wholesaleprices[ind2sub((dd,dc,ds),ka[t,1])[2]] for t=1:T]
    ss=[ind2sub((dd,dc,ds),ka[t,1])[3]-1 for t=1:T]
    p=plot(layout = (3,1),size=(1000,400))
    plot!(p[1],1:T,ss,label="stocks")
    plot!(p[2],1:T,cc,label="wholesale prices",l=:green)
    plot!(p[3],1:T,ddd,label="demand",l=:orange)
    display(p)
end