text:
########################################################################################
# Solution
########################################################################################
abstract type AbstractProbODESolution{T,N,S} <: DiffEqBase.AbstractODESolution{T,N,S} end
mutable struct ProbODESolution{
    T, N, uType, puType, uType2, DType, tType, rateType, xType, diffType, llType, P, A, IType, DE
} <: AbstractProbODESolution{T,N,uType}
    u::uType
    pu::puType
    u_analytic::uType2
    errors::DType
    t::tType
    k::rateType
    x::xType
    diffusions::diffType
    log_likelihood::llType
    prob::P
    alg::A
    interp::IType
    dense::Bool
    tslocation::Int
    destats::DE
    retcode::Symbol
end
ProbODESolution{T, N}(
    u, pu, u_analytic, errors, t, k, x, diffusions, log_likelihood, prob, alg, interp,
    dense, tslocation, destats, retcode) where {T, N} =
        ProbODESolution{
            T, N, typeof(u), typeof(pu), typeof(u_analytic), typeof(errors), typeof(t),
            typeof(k), typeof(x), typeof(diffusions), typeof(log_likelihood), typeof(prob),
            typeof(alg), typeof(interp), typeof(destats)}(
                u, pu, u_analytic, errors, t, k, x, diffusions, log_likelihood, prob, alg,
                interp, dense, tslocation, destats, retcode)

function DiffEqBase.solution_new_retcode(sol::ProbODESolution{T,N}, retcode) where {T,N}
    ProbODESolution{T, N}(sol.u, sol.pu, sol.u_analytic, sol.errors, sol.t, sol.k, sol.x,
                          sol.diffusions, sol.log_likelihood, sol.prob, sol.alg, sol.interp,
                          sol.dense, sol.tslocation, sol.destats, retcode)
end

# Used to build the initial empty solution in OrdinaryDiffEq.__init
function DiffEqBase.build_solution(
    prob::DiffEqBase.AbstractODEProblem,
    alg::GaussianODEFilter,
    t, u;
    k = nothing,
    retcode = :Default,
    destats = nothing,
    dense = true,
    kwargs...)

    T = eltype(eltype(u))
    N = length((size(prob.u0)..., length(u)))

    d = length(prob.u0)
    uEltype = eltype(prob.u0)
    cov = PSDMatrix(LowerTriangular(zeros(uEltype, d, d)))
    # cov = zeros(uEltype, d, d)
    pu = StructArray{Gaussian{Vector{eltype(prob.u0)}, typeof(cov)}}(undef, 0)
    x = StructArray{Gaussian{Vector{eltype(prob.u0)}, typeof(cov)}}(undef, 0)

    interp = GaussianODEFilterPosterior(alg, prob.u0)

    if DiffEqBase.has_analytic(prob.f)
        u_analytic = Vector{typeof(prob.u0)}()
        errors = Dict{Symbol,real(eltype(prob.u0))}()
    else
        u_analytic = nothing
        errors = nothing
    end

    ll = zero(uEltype)
    return ProbODESolution{T, N}(
        u, pu, u_analytic, errors, t, [], x, [], ll, prob, alg, interp, dense, 0, destats, retcode,
    )
end

function DiffEqBase.build_solution(sol::ProbODESolution{T,N}, u_analytic, errors) where {T,N}
    return ProbODESolution{T, N}(
        sol.u, sol.pu, u_analytic, errors, sol.t, sol.k, sol.x, sol.diffusions,
        sol.log_likelihood, sol.prob,
        sol.alg, sol.interp, sol.dense, sol.tslocation, sol.destats, sol.retcode,
    )
end

########################################################################################
# Compat with classic ODE solutions, to enable analysis with DiffEqDevTools.jl
########################################################################################
mutable struct MeanProbODESolution{
    T, N, uType, uType2, DType, tType, rateType, P, A, IType, DE, PSolType
} <: DiffEqBase.AbstractODESolution{T,N,uType}
    u::uType
    u_analytic::uType2
    errors::DType
    t::tType
    k::rateType
    prob::P
    alg::A
    interp::IType
    dense::Bool
    tslocation::Int
    destats::DE
    retcode::Symbol
    probsol::PSolType
end
function mean(sol::ProbODESolution{T,N}) where {T,N}
    return MeanProbODESolution{
        T, N, typeof(sol.u), typeof(sol.u_analytic), typeof(sol.errors),
        typeof(sol.t), typeof(sol.k), typeof(sol.prob), typeof(sol.alg),
        typeof(sol.interp), typeof(sol.dense), typeof(sol.tslocation),
        typeof(sol.destats), typeof(sol.retcode), typeof(sol)
    }(
        sol.u, sol.u_analytic, sol.errors, sol.t, sol.k, sol.prob,
        sol.alg, sol.interp, sol.dense, sol.tslocation, sol.destats,
        sol.retcode, sol)
end
(sol::MeanProbODESolution)(t::Real, deriv::Val{N}=Val(0)) where {N} =
    mean(sol.probsol(t, deriv))
(sol::MeanProbODESolution)(t::AbstractVector, deriv=Val(0)) =
    DiffEqArray(mean(sol.probsol(t, deriv).u), t)

########################################################################################
# Dense Output
########################################################################################
abstract type AbstractODEFilterPosterior <: DiffEqBase.AbstractDiffEqInterpolation end
struct GaussianODEFilterPosterior{SPType, AType, QType, PType} <: AbstractODEFilterPosterior
    d::Int
    q::Int
    SolProj::SPType
    A::AType
    Q::QType
    Precond::PType
    smooth::Bool
end
function GaussianODEFilterPosterior(alg, u0)
    uElType = eltype(u0)
    d = length(u0)
    q = alg.order

    Proj(deriv) = kron([i==(deriv+1) ? 1 : 0 for i in 1:q+1]', diagm(0 => ones(d)))
    SolProj = Proj(0)

    A, Q = ibm(d, q, uElType)
    Precond = preconditioner(d, q)
    GaussianODEFilterPosterior(
        d, q, SolProj, A, Q, Precond, alg.smooth)
end
DiffEqBase.interp_summary(interp::GaussianODEFilterPosterior) = "Gaussian ODE Filter Posterior"

function (posterior::GaussianODEFilterPosterior)(tval::Real, t, x, diffusions)
    @unpack A, Q, d, q, Precond = posterior

    if tval < t[1]
        error("Invalid t<t0")
    end
    if tval in t
        idx = sum(t .<= tval)
        @assert t[idx] == tval
        return x[idx]
    end

    idx = sum(t .<= tval)
    prev_t = t[idx]
    prev_rv = x[idx]
    diffusion = diffusions[minimum((idx, end))]

    # Extrapolate
    h1 = tval - prev_t
    P = Precond(h1)
    PI = inv(P)
    Qh = apply_diffusion(Q, diffusion)
    goal_pred = predict(P * prev_rv, A, Qh)
    goal_pred = PI * goal_pred

    if !posterior.smooth || tval >= t[end]
        return goal_pred
    end

    next_t = t[idx+1]
    next_smoothed = x[idx+1]

    # Smooth
    h2 = next_t - tval
    P = Precond(h2)
    PI = inv(P)
    goal_pred = P * goal_pred
    next_smoothed = P * next_smoothed
    Qh = apply_diffusion(Q, diffusion)

    goal_smoothed, _ = smooth(goal_pred, next_smoothed, A, Qh)

    return PI * goal_smoothed
end
function (sol::ProbODESolution)(t::Real, deriv::Val{N}=Val(0)) where {N}
    @unpack q, d = sol.interp
    return sol.interp.SolProj * sol.interp(t, sol.t, sol.x, sol.diffusions)
end
(sol::ProbODESolution)(t::AbstractVector, deriv=Val(0)) = DiffEqArray(StructArray(sol.(t, deriv)), t)

########################################################################################
# Plotting
########################################################################################
@recipe function f(sol::AbstractProbODESolution; ribbon_width=1.96)
    times = range(sol.t[1], sol.t[end], length=1000)
    dense_post = sol(times).u
    values = stack(mean(dense_post))
    stds = stack(std(dense_post))
    ribbon --> ribbon_width * stds
    xguide --> "t"
    yguide --> "u(t)"
    label --> hcat(["u$(i)(t)" for i in 1:length(sol.u[1])]...)
    return times, values
end

########################################################################################
# Sampling from a solution
########################################################################################
"""Helper function to sample from our covariances, which often have a "cross" of zeros
For the 0-cov entries the outcome of the sampling is deterministic!"""
function _rand(x::Gaussian, n::Int=1)
    m, C = x.μ, x.Σ
    @assert C isa PSDMatrix

    sample = m .+ C.L*randn(length(m), n)
    return sample
end

function sample_back(x_curr::Gaussian, x_next_sample::AbstractVector, Ah::AbstractMatrix, Qh::AbstractMatrix, PI=I)
    m_p, P_p = Ah*x_curr.μ, Ah*x_curr.Σ*Ah' + Qh
    P_p_inv = inv(Symmetric(P_p))
    Gain = x_curr.Σ * Ah' * P_p_inv

    m = x_curr.μ + Gain * (x_next_sample - m_p)

    P = X_A_Xt(x_curr.Σ, (I - Gain*Ah)) + X_A_Xt(Qh, Gain)

    assert_nonnegative_diagonal(P)
    return Gaussian(m, P)
end

function sample(sol::ProbODESolution, n::Int=1)
    sample(sol.t, sol.x, sol.diffusions, sol.t, sol.interp, n)
end
function sample(ts, xs, diffusions, difftimes, posterior, n::Int=1)

    @unpack A, Q, d, q, Precond = posterior
    E0 = kron([i==1 ? 1 : 0 for i in 1:q+1]', diagm(0 => ones(d)))
    dim = d*(q+1)

    x = xs[end]
    sample = _rand(x, n)
    @assert size(sample) == (dim, n)

    sample_path = zeros(length(ts), dim, n)
    sample_path[end, :, :] .= sample
    # @info "final value and samples" x.μ sample sample_path[end, :]

    for i in length(xs)-1:-1:1
        dt = ts[i+1] - ts[i]

        i_diffusion = sum(difftimes .<= ts[i])
        diffusion = diffusions[i_diffusion]

        Qh = apply_diffusion(Q, diffusion)
        P = Precond(dt)
        PI = inv(P)

        for j in 1:n
            sample_p = P*sample_path[i+1, :, j]
            x_prev_p = P*xs[i]

            prev_sample_p = sample_back(x_prev_p, sample_p, A, Qh, PI)

            # sample_path[i, :, j] .= PI*prev_sample_p.μ
            sample_path[i, :, j] .= PI*_rand(prev_sample_p)[:]
        end
    end

    return sample_path[:, 1:d, :]
end
function dense_sample(sol::ProbODESolution, n::Int=1)
    times = range(sol.t[1], sol.t[end], length=1000)
    states = StructArray([sol.interp(t, sol.t, sol.x, sol.diffusions) for t in times])

    sample(times, states, sol.diffusions, sol.t, sol.interp, n), times
end


text:
@testset "check_result.jl" begin
    @testset "_check_add!!_behavour" begin
        check = ChainRulesTestUtils._check_add!!_behavour

        check(10.0, 2.0)
        check(11.0, Zero())
        check([10.0, 20.0],  @thunk([2.0, 0.0]))

        check(12.0, InplaceableThunk(@thunk(2.0), X̄ -> error("Should not have in-placed")))

        check([10.0, 20.0], InplaceableThunk(
            @thunk([2.0, 0.0]),
            X̄ -> (X̄[1] += 2.0; X̄)
        ))
        @test fails(()->check([10.0, 20.0], InplaceableThunk(
            @thunk([2.0, 0.0]),
            X̄ -> (X̄[1] += 3.0; X̄),
        )))
    end
end

ext: jl


text:
function check_bipartiteness(A)
  return true
end

function check_unipartiteness(A)
  size(A,1) == size(A,2) || throw(ArgumentError("The matrix for a unipartite network must be square (yours has size $(size(A)))"))
end

function check_bipartiteness(A, T, B)
  check_bipartiteness(A)
  check_name_vector(T)
  check_name_vector(B)
  eltype(T) == eltype(B) || throw(ArgumentError("The species in a bipartite network levels must have the same type (you gave $(eltype(T)) and $(eltype(B)))"))
  length(unique(vcat(T,B))) == sum(size(A)) || throw(ArgumentError("Species names cannot be shared between the levels of a bipartite network"))
  length(T) == size(A,1) || throw(ArgumentError("The length of top-level species names must match matrix size ($(length(T)) v. $(size(A,1)))"))
  length(B) == size(A,2) || throw(ArgumentError("The length of bottom-level species names must match matrix size ($(length(B)) v. $(size(A,2)))"))
end

function check_unipartiteness(A, S)
  check_unipartiteness(A)
  check_name_vector(S)
  length(S) == size(A,1) || throw(ArgumentError("The length of species names must match matrix size ($(length(S)) v. $(size(A,1)))"))
end

function check_name_vector(N)
  length(N) == length(unique(N)) || throw(ArgumentError("The names of species must be unique"))
end

function check_probability_values(A)
  minimum(A) >= zero(eltype(A)) || throw(ArgumentError("The probabilities must be at least 0.0"))
  maximum(A) <= one(eltype(A)) || throw(ArgumentError("The probabilities must be at most 1.0"))
end


text:
# new Sandshark example
# add more julia processes
# nprocs() < 7 ? addprocs(7-nprocs()) : nothing

using Caesar, RoME, KernelDensityEstimate, IncrementalInference
using Interpolations
using Distributions

using RoMEPlotting
using Gadfly, DataFrames
using ProgressMeter

const TU = TransformUtils

Gadfly.set_default_plot_size(35cm,25cm)

include(joinpath(Pkg.dir("Caesar"), "examples", "marine", "auv", "Sandshark","Plotting.jl"))

# datadir = joinpath(ENV["HOME"],"data","sandshark","sample_wombat_2018_09_07","processed","extracted")
datadir = joinpath(ENV["HOME"],"data","sandshark","full_wombat_2018_07_09","extracted")
matcheddir = joinpath(datadir, "matchedfilter", "particles")
beamdir = joinpath(datadir, "beamformer", "particles")

function loaddircsvs(datadir)
  # https://docs.julialang.org/en/v0.6.1/stdlib/file/#Base.Filesystem.walkdir
  datadict = Dict{Int, Array{Float64}}()
  for (root, dirs, files) in walkdir(datadir)
    # println("Files in $root")
    for file in files
      # println(joinpath(root, file)) # path to files
      data = readdlm(joinpath(root, file),',')
      datadict[parse(Int,split(file,'.')[1])/1000000] = data
    end
  end
  return datadict
end

rangedata = loaddircsvs(matcheddir)
azidata = loaddircsvs(beamdir)
timestamps = intersect(sort(collect(keys(rangedata))), sort(collect(keys(azidata))))

# NAV data
navdata = Dict{Int, Vector{Float64}}()
navfile = readdlm(joinpath(datadir, "nav_data.csv"))
for row in navfile
    s = split(row, ",")
    id = round(Int, 1000*parse(s[1]))
    # round(Int, 1000 * parse(s[1])) = 1531153292381
    navdata[id] = parse.(s)
end
navkeys = sort(collect(keys(navdata)))
# NAV colums are X,Y = 7,8
# lat,long = 9,10
# time,pitch,roll,heading,speed,[Something], internal_x,internal_y,internal_lat,internal_long, yaw_rad

# LBL data - note the timestamps need to be exported as float in future.
lbldata = Dict{Int,  Vector{Float64}}()
lblfile = readdlm(joinpath(datadir, "lbl.csv"))
for row in lblfile
    s = split(row, ",")
    id = round(Int, 1000*parse(s[1]))
    @show s
    if s[2] != "NaN"
        lbldata[id] = parse.(s)
    end
end
lblkeys = sort(collect(keys(lbldata)))

# GET Y = north,  X = East,  Heading along +Y clockwise [0,360)]
# east = Float64[]
# north = Float64[]
# heading = Float64[]

# WANT X = North,  Y = West,  Yaw is right and rule from +X (0) towards +Y pi/2, [-pi,pi)
# so the drawPoses picture will look flipped from Nicks picture
# remember theta = atan2(y,x)    # this is right hand rule
X = Float64[]
Y = Float64[]
yaw = Float64[]
for id in navkeys
  # push!(east, getindex(navdata[id],7)) # x-column csv
  # push!(north, getindex(navdata[id],8)) # y-column csv
  # push!(heading, getindex(navdata[id],4))
  # push!(yaw, TU.wrapRad(-deg2rad(getindex(navdata[id],4))))

  push!(X, getindex(navdata[id],8) )
  push!(Y, -getindex(navdata[id],7) )
  push!(yaw, TU.wrapRad(-deg2rad(getindex(navdata[id],4))) )  # rotation about +Z
end

lblX = Float64[]
lblY = Float64[]
for id in lblkeys
    push!(lblX, getindex(lbldata[id],3) )
    push!(lblY, -getindex(lbldata[id],2) )
end

# Build interpolators for x, y, yaw
interp_x = LinearInterpolation(navkeys, X)
interp_y = LinearInterpolation(navkeys, Y)
interp_yaw = LinearInterpolation(navkeys, yaw)

## Caching factors
ppbrDict = Dict{Int, Pose2Point2Range}()
odoDict = Dict{Int, Pose2Pose2}()
NAV = Dict{Int, Vector{Float64}}()
# Step: Selecting a subset for processing and build up a cache of the factors.
epochs = timestamps[50:2:100]
lastepoch = 0
for ep in epochs
  if lastepoch != 0
    # @show interp_yaw(ep)
    deltaAng = interp_yaw(ep) - interp_yaw(lastepoch)

    wXi = TU.SE2([interp_x(lastepoch);interp_y(lastepoch);interp_yaw(lastepoch)])
    wXj = TU.SE2([interp_x(ep);interp_y(ep);interp_yaw(ep)])
    iDXj = se2vee(wXi\wXj)
    NAV[ep] = iDXj
    # println("$(iDXj[1]), $(iDXj[2]), $(iDXj[3])")

    odoDict[ep] = Pose2Pose2(MvNormal(NAV[ep], diagm([0.1;0.1;0.005].^2)))
  end
  rangepts = rangedata[ep][:]
  rangeprob = kde!(rangepts)
  # azipts = azidata[ep][:,1]
  # aziprob = kde!(azipts)

  # prep the factor functions
  ppbrDict[ep] = Pose2Point2Range(rangeprob)
  lastepoch = ep
end

## Step: Building the factor graph
fg = initfg()
l1Uncertainty = 0.1
nonParamStep = 1 # Number of poses between nonparametric acoustic factors
# Add a central beacon with a prior
addVariable!(fg, :l1, Point2)
addFactor!(fg, [:l1], IIF.Prior( MvNormal([0.6; -16], l1Uncertainty^2*eye(2)) ))
index = 0
for ep in epochs
    curvar = Symbol("x$index")
    addVariable!(fg, curvar, Pose2)

    # xi -> l1 - nonparametric factor
    if (index+1) % nonParamStep == 0
        info(" - Adding $curvar->l1 factor...")
        # addFactor!(fg, [curvar; (index < length(epochs)/2 ? :l1 : :l2)], ppbrDict[ep])
        addFactor!(fg, [curvar; :l1], ppbrDict[ep])
    end

    if ep != epochs[1]
      # Odo factor x(i-1) -> xi
      info(" - Adding x$(index-1)->$curvar odo factor")
      addFactor!(fg, [Symbol("x$(index-1)"); curvar], odoDict[ep])
    else
      # Prior to the first pose location (a "GPS" prior)
      initLoc = [interp_x(ep);interp_y(ep);interp_yaw(ep)]
      info(" - Adding a initial location prior at $curvar, $initLoc")
      addFactor!(fg, [curvar], IIF.Prior( MvNormal(initLoc, diagm([0.1;0.1;0.05].^2)) ))
    end
    # Heading partial prior
    info(" - Adding heading prior on $curvar")
    addFactor!(fg, [curvar], RoME.PartialPriorYawPose2(Normal(interp_yaw(ep), deg2rad(3))))
    index+=1
end

# writeGraphPdf(fg, engine="dot")

ensureAllInitialized!(fg)
solveTree!(fg)

# drawPosesLandms(fg)
drawPosesLandmarksAndOdo(fg, ppbrDict, navkeys, X, Y, lblkeys, lblX, lblY, "Multi-modal ISAM (bearing+range measurement every $nonParamStep poses, l1 uncertainty $l1Uncertainty)")


text:
using YAML, DataFrames, CSV, Plots
using Statistics
using Streamfall

HERE = @__DIR__
DATA_PATH = joinpath(HERE, "../test/data/hymod/")

# Load and generate stream network
network = YAML.load_file(joinpath(DATA_PATH, "hymod_network.yml"))
sn = create_network("HyMod Network", network)

# Load climate data
date_format = "YYYY-mm-dd"
obs_data = CSV.File(joinpath(DATA_PATH, "leaf_river_data.csv"),
                    comment="#",
                    dateformat=date_format) |> DataFrame

hist_streamflow = obs_data[:, ["Date", "leaf_river_outflow"]]
climate_data = obs_data[:, ["Date", "leaf_river_P", "leaf_river_ET"]]
climate = Climate(climate_data, "_P", "_ET")

run_basin!(sn, climate)

nid, node = sn["leaf_river"]
@info mean(node.outflow)


text:
module TestMPBWrapper

using Ipopt
using LinearAlgebra
using MathProgBase
using SparseArrays
using Test

function test_mpb()
    solver = IpoptSolver()
    sol = linprog([-1, 0], [2 1], '<', 1.5, solver)
    @test sol.status == :Optimal
    @test abs(sol.objval - -0.75) <= 1e-6
    @test norm(sol.sol - [0.75, 0.0]) <= 1e-6

    sol = linprog([-1, 0], sparse([2 1]), '<', 1.5, solver)
    @test sol.status == :Optimal
    @test abs(sol.objval - -0.75) <= 1e-6
    @test norm(sol.sol - [0.75, 0.0]) <= 1e-6

    # test infeasible problem:
    # min x
    # s.t. 2x+y <= -1
    # x,y >= 0
    sol = linprog([1, 0], [2 1], '<', -1, solver)
    #    @test sol.status == :Infeasible

    # test unbounded problem:
    # min -x-y
    # s.t. -x+2y <= 0
    # x,y >= 0
    return sol = linprog([-1, -1], [-1 2], '<', [0], solver)
    #    @test sol.status == :Unbounded
end

function mathprogbase_file(file::String)
    return joinpath(dirname(dirname(pathof(MathProgBase))), "test", file)
end

include(mathprogbase_file("nlp.jl"))
include(mathprogbase_file("quadprog.jl"))

function test_mathprogbase_nlp()
    nlptest(IpoptSolver())
    nlptest_nohessian(IpoptSolver())
    convexnlptest(IpoptSolver())
    return rosenbrocktest(IpoptSolver())
end

function test_mathprogbase_quadprog()
    quadprogtest(IpoptSolver())
    return qpdualtest(IpoptSolver())
end

function test_restoration_options()
    # Test retoration only options
    #
    # Warm start with infeasible solution, force restoration on initial iteration.
    # But limit to 0 iterations. Forces :UserLimit exit.
    m = MathProgBase.NonlinearModel(
        IpoptSolver(start_with_resto = "yes", resto_max_iter = 0),
    )
    l = [1, 1, 1, 1]
    u = [5, 5, 5, 5]
    lb = [25, 40]
    ub = [Inf, 40]
    MathProgBase.loadproblem!(m, 4, 2, l, u, lb, ub, :Min, HS071())
    MathProgBase.setwarmstart!(m, [0, 15, 15, 11])
    MathProgBase.optimize!(m)
    @test MathProgBase.status(m) == :UserLimit
end

end  # TestMPBWrapper

runtests(TestMPBWrapper)


text:
include("synonym.jl")
include("location.jl")
include("exon.jl")

@doc """ build index for genes
""" ->
function build_index()
    info("Building index now")
    @sync begin
        @spawn build_index_synonym()
        build_index_geneloc()
        build_index_exon()
    end
    info("Building index done")
    true
end

@doc """ load index for synonym genes
""" ->
function load_index()
    if !isfile(id_genes_fl) || !isfile(gene_id_fl)  ||
        !isfile(chr_sted_gene_fl) || !isfile(gene_chrsted_fl)
        build_index()
    end

    global id_genes = load(id_genes_fl, "id_genes_dict")
    global gene_id  = load(gene_id_fl,  "gene_id_dict")
    global chr_sted_gene = load(chr_sted_gene_fl, "chr_sted_gene_dict")
    global gene_chrsted  = load(gene_chrsted_fl, "gene_chrsted_dict")
    #gene_id,id_genes,chr_sted_gene,gene_chrsted

    nothing
end


text:
module webkit

  export tag, page, tag_table;

  mutable struct tag;
    name::Symbol;
    pairs::Bool;
    attrib::NamedTuple;
    attach::Array{tag};
    text::String;
    end;

  mutable struct page;
    doctype::tag;
    head::tag;
    body::tag;
    end;

## def

  tag(text::String) = 			  tag(:empty, false, NamedTuple(), [], text);
  tag(name::Symbol) = 			  tag(name, false, NamedTuple(), [], "");
  tag(name::Symbol, pairs::Bool) = 	  tag(name, pairs, NamedTuple(), [], "");
  tag(name::Symbol, text::String) = 	  tag(name, true, NamedTuple(), [], text);

  tag(name::Symbol, attach::Array{tag}) = tag(name, true, NamedTuple(), attach, "");
  tag(name::Symbol, attach::tag) = 	  tag(name, true, NamedTuple(), [attach], "");
  tag(name::Symbol, attrib::NamedTuple) = tag(name, true, attrib, [], "");

  tag(parent::tag, child::tag) = 	  push!(parent.attach, child);
  tag(parent::tag, childr::Array{tag}) =  append!(parent.attach, childr);

  schema_patch(z::Array{NamedTuple{Z, T}, N}) where{Z, T, N} = Z;
  schema_patch(z::Array{NamedTuple{Z}}) where{Z} = Z;

  function page(;doctype::tag = tag("<!DOCTYPE html><html>"), head::tag = tag(:head, true), body::tag = tag(:body, "Hello Webkit!"))
    page(doctype, head, body);
    end;

  function (x::page)(;generate = false)
    if(generate)
      io = IOBuffer();
      writetags(io, x.doctype);
      writetags(io, x.head);
      writetags(io, x.body);
      write(io, "</html>");
      return( String(take!(io)));
      end;
    nothing;
    end;

##

# low level api

function writetags(io::IO, xtag::tag)

  if(xtag.name !== :empty)
    write(io, "<$(xtag.name)");
    for (k, v) in pairs(xtag.attrib)
      write(io, " $k = \"$v\" ");
      end
    !xtag.pairs && write(io, " /");
    write(io, ">");
    end;

  write(io, xtag.text);
  for x in xtag.attach
    writetags(io, x);
    end;

  xtag.pairs && write(io, "</$(xtag.name)>");
  xtag.pairs && length(xtag.attach)>1 && write(io, "\n");
  end

function tag_table(table, caption="")

  xtable = tag(:table, true);
  (caption !== "") && tag(xtable, tag(:caption, caption));
  xlabels = hasmethod(schema_patch, (typeof(table), )) ? schema_patch(table) : [];
  thead = tag(:thead, true);
  [tag(thead, tag(:th, "$x")) for x in xlabels];
  tbody = tag(:tbody, true);

  for (id, xrow) in enumerate(table)
    local trow = tag(:tr, true);
    for elem in values(xrow)
      tag(trow, tag(:td, "$elem"));
      end;
    tag(tbody, trow);
    end;

  tag(xtable, thead);
  tag(xtable, tbody);
  return xtable;
  end;

  end # module


text:
# Find the Reverse Complement of a String

const NUCLEOTIDE_COMPLEMENTS = Dict('A' => 'T', 'C' => 'G', 'T' => 'A', 'G' => 'C')

calc_reverse_complement(genome) = join(map(n -> NUCLEOTIDE_COMPLEMENTS[n], Iterators.reverse(genome)), "")

function main()
    dna = readline()
    result = calc_reverse_complement(dna)
    println(result)
end

main()


text:
using Test
using Example2018

# Polynomials form a ring
# 1. addition is associative
# 2. additive identity
# 3. additive inverse
# 4. addition commutes

@testset "Polynomials form a ring" begin

	z = Polynomial([0])

	for n in 1:100

		deg = rand(0:100)
		p = Polynomial(rand(-1000000:1000000, deg+1))

		@test p + z == p
		@test z + p == p
		@test p + (-p) == z

	end

end

# Scaler Multiplication
# 1. A scaler multiplied by a polynomial gives a polynomial
# 2. A polynomial multiplied by zero gives zero

@testset "Scaler Multiplication" begin

	z = Polynomial([0])

	for n in 1:100

		deg = rand(0:100)
		sclr = rand()
		p = Polynomial(rand(-1000000:1000000, deg+1))

		@test poly(sclr*p) == true
		@test 0*p == z

	end

end

# Derivative Test
# 1. The derivative of a polynomial must be a polynomial
# 2. If the polynomial has degree n, the derivative must have degree (n-1)

@testset "Derivative Test" begin

	for n in 1:100

		deg = rand(0:100)
		p = Polynomial(rand(-1000000:1000000, deg+1))
		dp = deriv(p)

		@test poly(dp) == true
		@test length(dp.coeffs) == deg

	end

end


text:
const log2π = log(2π)
const _lgamma_cs = (0.99999999999980993,
                    676.5203681218851,
                    -1259.1392167224028,
                    771.32342877765313,
                    -176.61502916214059,
                    12.507343278686905,
                    -0.13857109526572012,
                    9.9843695780195716e-6,
                    1.5056327351493116e-7)

macro cufunc(ex)
    def = MacroTools.splitdef(ex)
    f = def[:name]
    def[:name] = Symbol(:cu, f)
    def[:body] = CuArrays.replace_device(def[:body])
    push!(CuArrays._cufuncs, f)
    quote
        $(esc(ex))
        $(esc(MacroTools.combinedef(def)))
        CuArrays.cufunc(::typeof($(esc(f)))) = $(esc(def[:name]))
    end
end

"""
logΓ function, computed via Lanczos approximation.
"""
function lgamma end

@cufunc function lgamma(x::T) where T<:Real
    s = zero(T)
    for k in 9:-1:2
        s += T(_lgamma_cs[k]) / (x + k - 2)
    end
    s += T(_lgamma_cs[1])
    t = x + T(13)/2
    T(log2π)/2 + (x-one(T)/2)log(t) - t + log(s)
end


text:
module PhaseRec
using FFTW
using LinearAlgebra
using Statistics
using Logging
using Base.Iterators

include("reconstruction.jl")
export phaserec, two_point

end # module


text:
@testset "chainrules" begin
    x = randn()
    z = x + exp(randn())
    y = z + exp(randn())
    test_frule(DistributionsAD.uniformlogpdf, x, y, z)
    test_rrule(DistributionsAD.uniformlogpdf, x, y, z)

    # StatsFuns: https://github.com/JuliaStats/StatsFuns.jl/pull/106
    x = exp(randn())
    y = exp(randn())
    z = logistic(randn())
    test_frule(StatsFuns.betalogpdf, x, y, z)
    test_rrule(StatsFuns.betalogpdf, x, y, z)

    x = exp(randn())
    y = exp(randn())
    z = exp(randn())
    test_frule(StatsFuns.gammalogpdf, x, y, z)
    test_rrule(StatsFuns.gammalogpdf, x, y, z)

    x = exp(randn())
    y = exp(randn())
    test_frule(StatsFuns.chisqlogpdf, x, y)
    test_rrule(StatsFuns.chisqlogpdf, x, y)

    x = exp(randn())
    y = exp(randn())
    z = exp(randn())
    test_frule(StatsFuns.fdistlogpdf, x, y, z)
    test_rrule(StatsFuns.fdistlogpdf, x, y, z)

    x = exp(randn())
    y = randn()
    test_frule(StatsFuns.tdistlogpdf, x, y)
    test_rrule(StatsFuns.tdistlogpdf, x, y)

    # TODO: Re-enable if https://github.com/JuliaMath/SpecialFunctions.jl/pull/325 is fixed
    # use `BigFloat` to avoid Rmath implementation in finite differencing check
    # (returns `NaN` for non-integer values)
    #n = rand(1:100)
    #x = BigFloat(n)
    #y = big(logistic(randn()))
    #z = BigFloat(rand(1:n))
    #test_frule(StatsFuns.binomlogpdf, x, y, z)
    #test_rrule(StatsFuns.binomlogpdf, x, y, z)

    x = big(exp(randn()))
    y = BigFloat(rand(1:100))
    test_frule(StatsFuns.poislogpdf, x, y)
    test_rrule(StatsFuns.poislogpdf, x, y)

    _, pb = rrule(StatsFuns.poislogpdf, 0.0, 0.0)
    _, x̄1, _ = pb(1)
    @test x̄1 == -1
    _, pb = rrule(StatsFuns.poislogpdf, 0.0, 1.0)
    _, x̄1, _ = pb(1)
    @test x̄1 == Inf

    # PoissonBinomial
    test_rrule(Distributions.poissonbinomial_pdf_fft, rand(50))
    if isdefined(Distributions, :poissonbinomial_pdf)
        test_rrule(Distributions.poissonbinomial_pdf, rand(50))
    end
end


text:
function convolute(a,b)
    len_a = length(a)
    len_b = length(b)
    len_c = len_a + len_b -1
    T_B = zeros(Float64, len_c, len_a)
    c = zeros(Float64, len_c)

    #Initialize T_B
    for j=1:len_a
        for i=j:len_c
            k = i-j+1
            if(k<=len_b)
                T_B[i,j] = b[k]
            else
                T_B[i,j] = 0
            end
        end
    end

    #Multiply T(B) with a to get c
    for i=1:len_c
        for j=1:1
            sum = 0
            for k=1:len_a
                sum+=T_B[i,k]*a[k]
            end
            c[i,j] = sum
        end
    end

    return c
end

u = Array{Int64}(undef, 50)
for i=1:50
    x = rand(Int,1)
    if(x[1]%2==0)
        u[i]=1
    else
        u[i]=-1
    end
end

c = [1,0.7,-0.3]
y = convolute(u, c)

using Plots
pyplot()
plot([y])

h = [0.9, -0.5, 0.5, -0.4, 0.3, -0.3, 0.2, -0.1]
z = convolute(y, h)
plot([z])

using DSP
a=[1]
b=[1,2]
c = conv(a,b)
print(c)


text:
function rbegin_func(arbiter_ptr::Ptr{Void}, space_ptr::Ptr{Void}, data::Ptr{Void})
    data = unsafe_pointer_to_objref(data)
    begin_func = data[1][1]
    userdata = data[2]
    arbiter = Arbiter(arbiter_ptr)
    space = Space(space_ptr)

    begin_func(arbiter, space, userdata)
end
function rpresolve_func(arbiter_ptr::Ptr{Void}, space_ptr::Ptr{Void}, data::Ptr{Void})
    data = unsafe_pointer_to_objref(data)
    presolve_func = data[1][2]
    userdata = data[2]
    arbiter = Arbiter(arbiter_ptr)
    space = Space(space_ptr)

    presolve_func(arbiter, space, userdata)
end
function rpostsolve_func(arbiter_ptr::Ptr{Void}, space_ptr::Ptr{Void}, data::Ptr{Void})
    data = unsafe_pointer_to_objref(data)
    postsolve_func = data[1][3]
    userdata = data[2]
    arbiter = Arbiter(arbiter_ptr)
    space = Space(space_ptr)

    postsolve_func(arbiter, space, userdata)
    return nothing
end
function rseparate_func(arbiter_ptr::Ptr{Void}, space_ptr::Ptr{Void}, data::Ptr{Void})
    data = unsafe_pointer_to_objref(data)
    separate_func = data[1][4]
    userdata = data[2]
    arbiter = Arbiter(arbiter_ptr)
    space = Space(space_ptr)

    separate_func(arbiter, space, userdata)
    return nothing
end

function CollisionHandler(typeA::Integer, typeB::Integer, begin_func::Function, presolve_func::Function,
                          postsolve_func::Function, separate_func::Function, userdata::Array{Any})

    callbacks = [begin_func, presolve_func, postsolve_func, separate_func]
    real_data = pointer_from_objref([callbacks, userdata])
    c_begin = cfunction(rbegin_func, Bool, (Ptr{Void}, Ptr{Void}, Ptr{Void}))
    c_presolve = cfunction(rpresolve_func, Bool, (Ptr{Void}, Ptr{Void}, Ptr{Void}))
    c_postsolve = cfunction(rpostsolve_func, Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}))
    c_separate = cfunction(rseparate_func, Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}))
    ccall(dlsym(libchipmunkjl, :cpCollisionHandler_create), Ptr{Void},
         (Uint32, Uint32, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}),
          typeA, typeB, c_begin, c_presolve, c_postsolve, c_separate, real_data)
end


text:
module RayTraceEllipsoids

using CoordinateTransformations
using StaticArrays: SVector, SDiagonal
using LinearAlgebra: normalize, ⋅, norm
# import Optim
# using OffsetArrays
# using WeightedOnlineStats
# using AngleBetweenVectors

export V3, Ray, EmptyRay, TraceRay, DeadRay, AbstractRay, Ellipsoid, Interface, OpticUnit, raytrace!, trace!, Medium, NonRetina, Retina, Interface, Glued, Active, Light, Cylinder, cspheroid, Signal

const ENERGYTHRESHOLD = 0.01

"""
    V3 = SVector{3,Float64}

Points and directions are just a point in 3D space best described as a static `SVector`. `V3` is an alias for that.
"""
const V3 = SVector{3, Float64}

abstract type AbstractRay end

struct DeadRay <: AbstractRay end

"""
    Ray(orig::V3, dir::V3, int::Float64)

The main Ray type with a ray origin, `orig`, a direction, `dir`, and an intensity `int`. The ray's direction gets normalized.
"""
struct Ray <: AbstractRay
    orig::V3
    dir::V3
    int::Float64
end
Ray(o, d) = Ray(o, d, 1.0)
# Ray() = Ray(zero(V3), zero(V3))

"""
    EmptyRay(orig::V3, dir::V3, int::Float64)

A Ray type for plotting. Same as `Ray` but doesn't get registered in retinas. Mainly used for plotting.
"""
struct EmptyRay <: AbstractRay
    orig::V3
    dir::V3
    int::Float64
end
EmptyRay(o, d) = EmptyRay(o, d, 1.0)

"""
    TraceRay(orig::V3, dir::V3, int::Float64)

A Ray type for plotting. Same as `Ray` but doesn't get registered in retinas. Mainly used for plotting.
"""
struct TraceRay <: AbstractRay
    orig::V3
    dir::V3
end

# TraceRay() = TraceRay(zero(V3), zero(V3))

# struct PSFRay <: AbstractRay
    # orig::V3
    # dir::V3
    # int::Float64
    # PSFRay(o::V3, d::V3, int::Float64) = new(o, d, int)
# end
# PSFRay() = PSFRay(zero(V3), zero(V3))

"""
    Ellipsoid(c::V3, r::V3, dir::V3, open::Float64)

    An ellipsoid with a center, `c`, and radii, `r`, as well as a direction (gets automatically normalized), `dir`, and an opening angle, `open`, creating a dome (or window), that the ellipsoid is defined in. Note that `open` is the angle between the dome's edge and the direction of the dome (so actually half the opening angle) and is defined in **some angular units** (using UnitfulAngles, for example: u"°").

`Ellipsoid` has 6 additional fields all relating to various spatial transformations that convert the ellipsoid to a unit-sphere and back again. These are all `CoordinateTransformations`.
```
"""

struct Ellipsoid
    positive::Bool # is the (normalized) dome in the positive end of the z-axis?
    center::Translation{V3} # translate the ellipsoid to zero
    scale::LinearMap{SDiagonal{3,Float64}} # scale to a unit-sphere
    center_scale::AffineMap{SDiagonal{3,Float64},V3} # translate and scale to a unit-sphere
    uncenter_unscale::AffineMap{SDiagonal{3,Float64},V3}

    function Ellipsoid(c::V3, r::V3, positive::Bool)
        uncenter = Translation(c)
        unscale = LinearMap(SDiagonal(r))
        center = inv(uncenter)
        scale = inv(unscale)
        center_scale = scale∘center
        uncenter_unscale = inv(center_scale)
        new(positive, center, scale, center_scale, uncenter_unscale)
    end

end

cspheroid(cz::Real, rxy::Real, rz::Real, h::Bool) = Ellipsoid(V3(0.0, 0.0, cz), V3(rxy, rxy, rz), h)

isgoodz(oz, l, dz, positive::Bool) = oz > -l*dz ? positive : !positive

"""
    distance(orig::V3, dir::V3)

Return the two distances between a point with origin, `orig`, and direction, `dir`, and the two (potentially identical) intersection points with a unit-sphere.
"""
function distance(orig::V3, dir::V3, positive::Bool)
    b = -orig⋅dir
    disc = b^2 - orig⋅orig + 1
    if disc ≥ 0
        d = sqrt(disc)
        t2 = b + d
        if t2 ≥ 0
            t1 = b - d
            if t1 > 0
                isgoodz(orig[3], t1, dir[3], positive) && return t1
                isgoodz(orig[3], t2, dir[3], positive) && return t2
            else
                isgoodz(orig[3], t2, dir[3], positive) && return t2
            end
        end
    end
    nothing
end

function distance(r::AbstractRay, s::Ellipsoid)
    orig = s.center_scale(r.orig) # transform the ray's origin according to the ellipsoid
    dir = normalize(s.scale(r.dir)) # scale the ray's direction as well
    l = distance(orig, dir, s.positive)
    untransform(l, orig, dir, r.orig, s.uncenter_unscale)
end
untransform(l::Nothing, orig, dir, rorig, uncenter_unscale) = l
function untransform(l::Float64, orig, dir, rorig, uncenter_unscale)
    o = orig + l*dir
    p = uncenter_unscale(o)
    norm(p - rorig)
end

struct Cylinder
    scale::LinearMap{SDiagonal{3,Float64}}
    center_scale::AffineMap{SDiagonal{3,Float64},V3}
    uncenter_unscale::AffineMap{SDiagonal{3,Float64},V3}

    function Cylinder(rxy::Real, z1::Real, z2::Real)
        cz = (z1 + z2)/2
        uncenter = Translation(V3(0.0, 0.0, cz))
        @assert z1 < z2 "cylinder z1 seems bigger than z2"
        rz = (z2 - z1)/2
        unscale = LinearMap(SDiagonal(rxy, rxy, rz))
        center = inv(uncenter)
        scale = inv(unscale)
        center_scale = scale∘center
        uncenter_unscale = inv(center_scale)
        new(scale, center_scale, uncenter_unscale)
    end

end

function _getl(orig::V3, l::Float64, dir::V3, ucs, rorig::V3)
    pᵗ = orig + l*dir
    p2 = ucs(pᵗ)
    norm(p2 - rorig)
end

function incylinder(r::Ray, s::Cylinder, L::Float64)
    orig = s.center_scale(r.orig)
    dir = normalize(s.scale(r.dir))
    a = dir[1]^2 + dir[2]^2
    if a == 0
        if orig[1]^2 + orig[2]^2 < 1
            return L
        else
            return 0.0
        end
    end
    b = 2orig[1]*dir[1] + 2orig[2]*dir[2]
    c = orig[1]^2 + orig[2]^2 - 1
    Δ = b^2 - 4a*c
    if Δ ≤ 0
        return 0.0
    end
    sqrtΔ = sqrt(Δ)
    l1 = (-b - sqrtΔ)/2a
    l2 = (-b + sqrtΔ)/2a
    l1, l2 = l1 < l2 ? (l1, l2) : (l2, l1)
    if l2 < 0
        return 0.0
    end
    if l1 < 0
        l = _getl(orig, l2, dir, s.uncenter_unscale, r.orig)
        if 0 < l ≤ L
            return l
        else
            return L
        end
    else
        l1 = _getl(orig, l1, dir, s.uncenter_unscale, r.orig)
        if l1 > L
            return 0.0
        end
        l2 = _getl(orig, l2, dir, s.uncenter_unscale, r.orig)
        if l2 > L
            return L - l1
        else
            return l2 - l1
        end
    end
end

abstract type Medium end

struct NonRetina <: Medium
    absorption_coefficient::Float64
end

mutable struct Signal
    photoreceptor::Float64
    retina::Float64
    Signal() = new(0.0, 0.0)
end

struct Retina <: Medium
    cylinder::Cylinder
    absorption_coefficient::Float64
    signal::Signal
    # photons::Base.RefValue{Int}
    # pixels::OffsetArray{Float64}
    # n::Base.RefValue{Int}
    # var::WeightedMean{Float64}
    # σ²::WeightedVariance{Float64}
    # data::Vector{Float64}
    # weight::Vector{Float64}
    # nodalz::Float64
    function Retina(cylinder::Cylinder, absorption_coefficient::Float64)
        # new(absorption_coefficient, WeightedMean{Float64}(), Float64[], Float64[], nodalz)
        # pixels = OffsetArray{Float64}(undef, -400:400, -400:400)
        # fill!(pixels, 0.0)
        new(cylinder, absorption_coefficient, Signal())#, WeightedVariance{Float64}())
    end
end

const LOG10EXP = 1/(log10(exp(1)))

function absorption(l::Real, absorption_coefficient::Real)
    kl = absorption_coefficient*l
    kl/(LOG10EXP + kl)
end
# absorption(l::Float64, absorption_coefficient::Float64) = exp(-l*absorption_coefficient)

function register!(m::Retina, r::Ray, l, w)
    linside = incylinder(r, m.cylinder, l)
    m.signal.photoreceptor += r.int*absorption(linside, m.absorption_coefficient)
    m.signal.retina += w
    # if linside > 0
    #     m.photons[] += 1
    # end
    # fit!(m.σ², r.orig[1], w/2)
    # fit!(m.σ², r.orig[2], w/2)
    nothing
end

register!(m, r, l, w) = nothing

abstract type Interface end

struct Glued <: Interface
    passthrough::Float64
end

"""
    Interface(normal::AffineMap, n::Float64)

Build an optical interface from a AffineMap that transforms a point on the ellipsoid of the interface to the normal at that point, `normal`, and the refractive index ratio between the inside and the outside of the ellipsoid, `n`.
"""
struct Active <: Interface
    passthrough::Float64
    normal::AffineMap{SDiagonal{3, Float64}, V3}
    n::Float64
    n²::Float64
    Active(pass::Float64, normal::AffineMap{SDiagonal{3, Float64}, V3}, n::Float64) = new(pass, normal, n, n^2)
end

refract(n, dir, a, b, N) = n*dir + (n*a - sqrt(1 - b))*N

reflect(dir, a, N) = dir + 2a*N

"""
    newdirection(r::Ray, i::Interface)

Refract or reflect a ray with an interface. Update the direction of the ray. Returns if event failed, which is always `true` (see `raytrace!` for details why that is so).
"""
function newdirection(orig, dir, i::Active)
    N = normalize(i.normal(orig))
    a = -dir⋅N
    b = i.n²*(1 - a^2)
    dir = if b ≤ 1
        refract(i.n, dir, a, b, N)
    else
        reflect(dir, a, N)
    end
    normalize(dir)
end

newdirection(orig, dir, i::Glued) = dir

"""
    OpticUnit(surface::Ellipsoid, pointin::Bool, n::Float64, register::Bool, name::String)

Build an optical unit from a surface. `pointin` dictates if the normal should be pointing in or out. `n` is the refractive index. `register` indicates whether this unit should register intersection points (thus functions as a retina).
"""
struct OpticUnit{T <: Medium}
    surface::Ellipsoid
    medium::T
    interface::Interface

    function OpticUnit(surface::Ellipsoid, pointin::Bool, n::Float64, medium::T, passthrough::Float64) where {T<:Medium}
        interface = if n ≠ 1
            i = Float64((-1)^pointin)
            dir = LinearMap(SDiagonal(i, i, i))
            normal = dir∘surface.scale∘surface.scale∘surface.center
            Active(passthrough, normal, n)
        else
            Glued(passthrough)
        end
        new{T}(surface, medium, interface)
    end
end

trace!(c::OpticUnit, r::DeadRay) = r

"""
    trace!(c::OpticUnit, r::Ray)

Advance a ray to the intersection point with an ellipsoid. If the intersection was successful, bend the ray. Updates the ray accordingly. Returns intersection failure.
"""
function trace!(c::OpticUnit, r::AbstractRay)
    l = distance(r, c.surface)
    mutate(l, r, c)
end

newlocation(orig, l, dir) = orig + l*dir

mutate(l::Nothing, r, c) = DeadRay()
function mutate(l::Float64, r::TraceRay, c)
    orig = newlocation(r.orig, l, r.dir)
    dir = newdirection(orig, r.dir, c.interface)
    TraceRay(orig, dir)
end
function mutate(l::Float64, r::AbstractRay, c)
    w = absorbed(r, l, c.medium)
    int = newintensity(r, w)
    _mutate(int, w, r, l, c)
end
_mutate(int::Nothing, w, r, l, c) = DeadRay()
function _mutate(int, w, r, l, c)
    register!(c.medium, r, l, w)
    orig = newlocation(r.orig, l, r.dir)
    int = passed(c.interface.passthrough, int)
    __mutate(int, orig, r, c.interface)
end

function __mutate(int::Float64, orig, r, interface)
    dir = newdirection(orig, r.dir, interface)
    typeof(r)(orig, dir, int)
end

__mutate(int::Nothing, orig, r, interface) = DeadRay()

function passed(passthrough, int)
    int *= passthrough
    int < ENERGYTHRESHOLD ?  nothing : int
end

absorbed(r, l, m) = r.int*absorption(l, m.absorption_coefficient)

function newintensity(r, w)
    int = r.int - w
    int < ENERGYTHRESHOLD ?  nothing : int
end

raytrace!(ous, r::DeadRay) = r
function raytrace!(ous, r::AbstractRay)
    @inbounds for ou in ous
        r = trace!(ou, r)
    end
    r
end

# light

struct Light
    cosa::Float64
    orig::V3
    rotm::LinearMap{RotY{Float64}}
    nom::Float64
    radius²::Float64

end
function Light(distance::Real, aperture::Real, rc, θ::Real, nodalz::Real)
    nodal = V3(0, 0, nodalz)
    radius = aperture/2
    h = rc.rz*sqrt(1 - (radius/rc.rxy)^2) + rc.cz
    vs = (V3(0, radius, h), V3(radius, 0, h), V3(-radius, 0, h))
    rotm = LinearMap(RotY{Float64}(θ))
    tran = recenter(rotm, nodal)
    orig = tran(V3(0, 0, nodalz + distance))
    v1 = normalize(nodal - orig)
    cosa = Inf
    for vi in vs
        _cosa = getcosa(vi, orig, v1)
        if _cosa < cosa
            cosa = _cosa
        end
    end
    nom = h - orig[3]
    Light(cosa, orig, rotm, nom, radius^2)
end
#=function Light(distance::Real, aperture::Real, rc, θ::Real, γ::Real, nodalz::Real)

    =##=nodal = V3(0, 0, nodalz)
    rotm = LinearMap(RotZY{Float64}(γ, θ))
    tran = recenter(rotm, nodal)
    orig = tran(V3(0, 0, nodalz + distance))

    v1 = normalize(nodal - orig)
    radius = aperture/2
    h = rc.rz*sqrt(1 - (radius/rc.rxy)^2) + rc.cz
    res = Optim.optimize(ξ -> ξ2v(ξ, radius, h, orig, v1), 0, π)
    cosa = Optim.minimum(res)
    acosd(cosa)=##=

    nodal = V3(0, 0, nodalz)
    radius = aperture/2
    h = rc.rz*sqrt(1 - (radius/rc.rxy)^2) + rc.cz
    vs = (V3(0, radius, h), V3(radius, 0, h), V3(-radius, 0, h))
    roty = LinearMap(RotY{Float64}(θ))
    trany = recenter(roty, nodal)
    origy = trany(V3(0, 0, nodalz + distance))
    v1 = normalize(nodal - origy)
    # cosa = minimum(getcosa(vi, origy, v1) for vi in vs)
    cosa = Inf
    for vi in vs
        _cosa = getcosa(vi, origy, v1)
        if _cosa < cosa
            cosa = _cosa
        end
    end
    rotm = LinearMap(RotZY{Float64}(γ, θ))
    tran = recenter(rotm, nodal)
    orig = tran(V3(0, 0, nodalz + distance))
    nom = h - orig[3]
    new(cosa, orig, rotm, nom, radius^2)
end=#

function getcosa(vi, origy, v1)
    vin = normalize(vi - origy)
    vin⋅v1
end

#=function ξ2v(ξ, radius, h, orig, v1)
sinξ, cosξ = sincos(ξ)
v = V3(radius*cosξ, radius*sinξ, h) - orig
normalize(v)⋅v1
end=#

#=function Light(distance::Real, aperture::Real, rxy::Real, rz::Real, cz::Real, θ::Real)
radius = aperture/2
h = rz*sqrt(1 - (radius/rxy)^2)
nodal = V3(0, 0, cz + h)
rotm = LinearMap(RotY{Float64}(θ))
tran = recenter(rotm, nodal)
orig = tran(nodal + V3(0, 0,distance))
v1 = nodal - orig
v2 = V3(radius, 0, h + cz) - orig
cosalpha1 = normalize(v1)⋅normalize(v2)
v2 = V3(-radius, 0, h + cz) - orig
cosalpha2 = normalize(v1)⋅normalize(v2)
cosalpha3 = cos(atan(radius/distance))
cosa = min(cosalpha1, cosalpha2, cosalpha3)
z = cz + h
nom = z - orig[3]
Light(cosa, orig, rotm, nom, nodal, radius + 1e-7) # adding some fuzz factor to the radius
end=#

function (l::Light)(::Type{T}, b1::Real, b2::Real) where {T <: AbstractRay}
    z = b1*(l.cosa - 1) - l.cosa
    k = sqrt(1 - z^2)
    theta = 2b2
    dir = l.rotm(V3(k*cospi(theta), k*sinpi(theta), z))
    t = l.nom/dir[3]
    p0 = l.orig + t*dir
    if p0[1]^2 + p0[2]^2 ≤ l.radius²
        T(l.orig, dir)
    else
        DeadRay()
    end
end

function (l::Light)(::Type{T}, b::Real) where {T <: AbstractRay}
    α = acos(l.cosa)
    a = 2α*(b - 0.5)
    dir = l.rotm(normalize(V3(sin(a), 0.0, -cos(a))))
    t = l.nom/dir[3]
    p0 = l.orig + t*dir
    if p0[1]^2 + p0[2]^2 ≤ l.radius²
        T(l.orig, dir)
    else
        DeadRay()
    end
end

(l::Light)(::Type{T}) where {T <: AbstractRay} = l(T, rand(), rand())

function raytrace!(ous, l::Light; n=10^5) # after some testing I found that 10^5 rays is enough for stable stats
    ni = Ref(0)
    @inbounds for _ in 1:n
        r = l(Ray)
        addone!(ni, r)
        raytrace!(ous, r)
    end
    ni[]
end

addone!(ni, r::DeadRay) = nothing
function addone!(ni, r::AbstractRay)
    ni[] += 1
    nothing
end

end #module


text:
using Test

@testset "Barnes-Hut algorithm" begin
    include("../src/barnes_hut.jl")

    @testset "space fitting" begin
        quad = Quad(-2, -2, 2, 2)
        @test fit(quad, [1, 1]) == :ne
        @test fit(quad, [-1, 1]) == :nw
        @test fit(quad, [1, -1]) == :se
        @test fit(quad, [-1, -1]) == :sw

        @test fit(quad, [0, 0]) == :ne
        @test fit(quad, [1, 0]) == :ne
        @test fit(quad, [0, 1]) == :ne
        @test_throws DomainError fit(quad, [3, 3])

        @test quadrant_from_symbol(quad, :ne) == Quad(0,0,2,2)
        @test quadrant_from_symbol(quad, :nw) == Quad(-2,0,0,2)
        @test quadrant_from_symbol(quad, :se) == Quad(0,-2,2,0)
        @test quadrant_from_symbol(quad, :sw) == Quad(-2,-2,0,0)
    end
    @testset "tree creation" begin

    end
end


text:
export
    Unconstrained,
    Polygonal,
    nextboundary

abstract type Domain end

"""
    Unconstrained

Signature of domains without boundaries.
"""
struct Unconstrained <: Domain end

"""
    Polygonal

Encapsulate a domain with polygonal boundaries. The boundaries are defined with
the normals to the facets and the intercepts of the corresponding hyperplanes.
"""
struct Polygonal <: Domain
    normals::Matrix{Real} # dim CxD where C=#constraints, D=#features
    intercepts::Vector{Real} # dim C
end

# -----------------------------------------------------------------------------

"""
    nextboundary(ud::Unconstrained, x, v)

Return (NaN, NaN) corresponding to the unconstrained case.
"""
function nextboundary(ud::Unconstrained,
                      x::Vector{<:Real}, v::Vector{<:Real})
    return (NaN, NaN)
end

"""
    nextboundary(pd::Polygonal, x, v)

Return the time of next boundary hit along the current trajectory when the
domain is polygonal and return the normal of the corresponding boundary.
The current point is `x` and the velocity `v`.
"""
function nextboundary(pd::Polygonal,
                      x::Vector{<:Real}, v::Vector{<:Real})
    # hitting time along trajectory (x+tv) for a boundary (normal,intercept) is
    # t = intercept/(normal dot v) - (x dot normal)/(normal dot v)
    nsv  = pd.normals * v
    hits = (pd.intercepts - pd.normals * x) ./ nsv
    # hard threshold times with nsv ~ 0 (near-parallel case),
    # remove negative times and times ~ 0 (for numerical stability)
    mask = map(|, abs.(nsv) .< 1e-10, hits .< 1e-10)
    for (i, e) ∈ enumerate(mask)
        e && (hits[i] = Inf)
    end
    # get time of hit + index of corresponding boundary
    (t_hit, j) = findmin(hits)
    # return time of hit + normal vector to boundary
    (t_hit, pd.normals[j, :])
end


text:
# # include this first to help with crashing??
# using Gtk

# don't let pyplot use a gui... it'll crash
# note: Agg will set gui -> :none in PyPlot
ENV["MPLBACKEND"] = "Agg"
try
  @eval import PyPlot
end

# macro test_approx_eq_sigma_eps(A, B, sigma, eps)

include("../docs/example_generation.jl")

# # make and display one plot
# function test_examples(pkg::Symbol, idx::Int; debug = true)
#   Plots._debugMode.on = debug
#   println("Testing plot: $pkg:$idx:$(examples[idx].header)")
#   backend(pkg)
#   backend()
#   map(eval, examples[idx].exprs)
#   plt = current()
#   gui(plt)
#   plt
# end

using Plots, FactCheck
import Images, ImageMagick

# if !isdefined(ImageMagick, :init_deps)
#   function ImageMagick.init_deps()
#     ccall((:MagickWandGenesis,libwand), Void, ())
#   end
# end

function makeImageWidget(fn)
  img = Gtk.GtkImageLeaf(fn)
  vbox = Gtk.GtkBoxLeaf(:v)
  push!(vbox, Gtk.GtkLabelLeaf(fn))
  push!(vbox, img)
  show(img)
  vbox
end

function replaceReferenceImage(tmpfn, reffn)
  cmd = `cp $tmpfn $reffn`
  run(cmd)
  info("Replaced reference image with: $cmd")
end

"Show a Gtk popup with both images and a confirmation whether we should replace the new image with the old one"
function compareToReferenceImage(tmpfn, reffn)
  @eval import Gtk

  # add the images
  imgbox = Gtk.GtkBoxLeaf(:h)
  push!(imgbox, makeImageWidget(tmpfn))
  push!(imgbox, makeImageWidget(reffn))

  # add the buttons
  # doNothingButton = Gtk.GtkButtonLeaf("Skip")
  # replaceReferenceButton = Gtk.GtkButtonLeaf("Replace reference image")
  # btnbox = Gtk.GtkButtonBoxLeaf(:h)
  # push!(btnbox, doNothingButton)
  # push!(btnbox, replaceReferenceButton)

  # # create the window
  # box = Gtk.GtkBoxLeaf(:v)
  # push!(box, imgbox)
  # push!(box, btnbox)
  win = Gtk.GtkWindowLeaf("Should we make this the new reference image?")
  push!(win, Gtk.GtkFrameLeaf(imgbox))

  showall(win)

  # now ask the question
  if Gtk.ask_dialog("Should we make this the new reference image?", "No", "Yes")
    replaceReferenceImage(tmpfn, reffn)
  end

  destroy(win)

  # # we'll wait on this condition
  # c = Condition()
  # Gtk.on_signal_destroy((x...) -> notify(c), win)

  # Gtk.signal_connect(replaceReferenceButton, "clicked") do widget
  #   replaceReferenceImage(tmpfn, reffn)
  #   notify(c)
  # end

  # Gtk.signal_connect(doNothingButton, "clicked") do widget
  #   notify(c)
  # end

  # # wait until a button is clicked, then close the window
  # Gtk.showall(win)
  # wait(c)
  # Gtk.destroy(win)
end

# TODO: use julia's Condition type and the wait() and notify() functions to initialize a Window, then wait() on a condition that
#       is referenced in a button press callback (the button clicked callback will call notify() on that condition)

function image_comparison_tests(pkg::Symbol, idx::Int; debug = false, sigma = [1,1], eps = 1e-2)

  # first
  Plots._debugMode.on = debug
  info("Testing plot: $pkg:$idx:$(PlotExamples.examples[idx].header)")
  backend(pkg)
  backend()

  # ensure consistent results
  srand(1234)

  # run the example
  map(eval, PlotExamples.examples[idx].exprs)

  # save the png
  tmpfn = tempname() * ".png"
  png(tmpfn)

  # load the saved png
  tmpimg = Images.load(tmpfn)

  # reference image location
  # refdir = joinpath(Pkg.dir("Plots"), "test", "refimg", "v$(VERSION.major).$(VERSION.minor)", string(pkg))
  refdir = joinpath(Pkg.dir("Plots"), "test", "refimg", string(pkg))
  try
    run(`mkdir -p $refdir`)
  catch err
    display(err)
  end
  reffn = joinpath(refdir, "ref$idx.png")

  try

    # info("Comparing $tmpfn to reference $reffn")

    # load the reference image
    refimg = Images.load(reffn)

    # run the comparison test... a difference will throw an error
    # NOTE: sigma is a 2-length vector with x/y values for the number of pixels
    #       to blur together when comparing images
    diffpct = Images.test_approx_eq_sigma_eps(tmpimg, refimg, sigma, eps)

    # we passed!
    info("Reference image $reffn matches.  Difference: $diffpct")
    return true

  catch err
    warn("Image did not match reference image $reffn. err: $err")
    showerror(Base.STDERR, err)

    if isinteractive()

      # if we're in interactive mode, open a popup and give us a chance to examine the images
      warn("Should we make this the new reference image?")
      compareToReferenceImage(tmpfn, reffn)
      return

    else

      # if we rejected the image, or if we're in automated tests, throw the error
      rethrow(err)
    end

  end
end

function image_comparison_tests(pkg::Symbol; skip = [], debug = false, sigma = [1,1], eps = 1e-2)
  for i in 1:length(PlotExamples.examples)
    i in skip && continue
    @fact image_comparison_tests(pkg, i, debug=debug, sigma=sigma, eps=eps) --> true
  end
end


text:
using EEGTest
using Base.Test

# write your own tests here
include("eegtests.jl")


text:
## Basic weight functions

export LebesgueMeasure,
    Lebesgue,
    LebesgueDomain,
    lebesguemeasure,
    dx,
    LegendreWeight,
    JacobiWeight,
    ChebyshevWeight,
    ChebyshevTWeight,
    ChebyshevUWeight,
    LaguerreWeight,
    HermiteWeight,
    GaussianWeight,
    DiracWeight,
    point,
    jacobi_α,
    jacobi_β,
    laguerre_α

####################
# Lebesgue measures
####################

"Supertype of Lebesgue measures"
abstract type LebesgueMeasure{T} <: Weight{T} end

unsafe_weightfun(μ::LebesgueMeasure, x) = one(codomaintype(μ))

islebesguemeasure(m::Measure) = false
islebesguemeasure(m::LebesgueMeasure) = true

==(μ1::LebesgueMeasure, μ2::LebesgueMeasure) = support(μ1) == support(μ2)

"The Lebesgue measure on the space `FullSpace{T}`."
struct Lebesgue{T} <: LebesgueMeasure{T}
end

Lebesgue() = Lebesgue{Float64}()
similar(μ::Lebesgue, ::Type{T}) where {T} = Lebesgue{T}()

Base.show(io::IO, μ::Lebesgue) = print(io, "dx")

"The Lebesgue measure on the unit interval `[0,1]`."
struct LebesgueUnit{T} <: LebesgueMeasure{T}
end

LebesgueUnit() = LebesgueUnit{Float64}()
similar(μ::LebesgueUnit, ::Type{T}) where {T <: Real} = Lebesgue{T}()
support(μ::LebesgueUnit{T}) where {T} = UnitInterval{T}()

isnormalized(μ::LebesgueUnit) = true

Base.show(io::IO, μ::LebesgueUnit) = print(io, "dx(0..1)")

"Lebesgue measure supported on a general domain."
struct LebesgueDomain{T} <: LebesgueMeasure{T}
    domain  ::  Domain{T}
end

similar(μ::LebesgueDomain, ::Type{T}) where {T} = LebesgueDomain{T}(μ.domain)
support(m::LebesgueDomain) = m.domain

Base.show(io::IO, mime::MIME"text/plain", μ::LebesgueDomain) = composite_show(io, mime, μ)
Display.displaystencil(μ::LebesgueDomain) = ["dx(", support(μ), ")"]

###############
# Dirac weight
###############

"A continuous Dirac measure at a point"
struct DiracWeight{T} <: Weight{T}
    point   ::  T
end

similar(μ::DiracWeight, ::Type{T}) where {T} = DiracWeight{T}(μ.point)

point(μ::DiracWeight) = μ.point
support(μ::DiracWeight) = Point(μ.point)
isnormalized(μ::DiracWeight) = true
unsafe_weightfun(μ::DiracWeight, x) = convert(codomaintype(μ), Inf)

Base.show(io::IO, μ::DiracWeight) = iszero(point(μ)) ? print(io, "δ(x)dx") : print(io, "δ(x-$(repr(point(μ))))dx")

###################################
# Classical orthogonal polynomials
###################################

"The Legendre weight is a Lebesgue measure on the interval `[-1,1]`."
struct LegendreWeight{T} <: LebesgueMeasure{T}
end
LegendreWeight() = LegendreWeight{Float64}()

similar(μ::LegendreWeight, ::Type{T}) where {T <: Real} = LegendreWeight{T}()
support(μ::LegendreWeight{T}) where {T} = ChebyshevInterval{T}()

Base.show(io::IO, μ::LegendreWeight) = print(io, "dx(-1..1)  (Legendre)")
Display.object_parentheses(μ::LegendreWeight) = true

"Supertype of Jacobi weights."
abstract type AbstractJacobiWeight{T} <: Weight{T} end

==(μ1::AbstractJacobiWeight, μ2::AbstractJacobiWeight) =
    jacobi_α(μ1) == jacobi_α(μ2) && jacobi_β(μ1) == jacobi_β(μ2)

"The Jacobi weight on the interval `[-1,1]`."
struct JacobiWeight{T} <: AbstractJacobiWeight{T}
    α   ::  T
    β   ::  T

    JacobiWeight{T}(α = zero(T), β = zero(T)) where {T} = new(α, β)
end
JacobiWeight() = JacobiWeight{Float64}()
JacobiWeight(α, β) = JacobiWeight(promote(α, β)...)
JacobiWeight(α::T, β::T) where {T<:AbstractFloat} = JacobiWeight{T}(α, β)
JacobiWeight(α::N, β::N) where {N<:Number} = JacobiWeight(float(α), float(β))

jacobi_weightfun(x, α, β) = (1+x)^α * (1-x)^β

similar(μ::JacobiWeight, ::Type{T}) where {T <: Real} = JacobiWeight{T}(μ.α, μ.β)
support(μ::JacobiWeight{T}) where {T} = ChebyshevInterval{T}()
unsafe_weightfun(μ::JacobiWeight, x) = jacobi_weightfun(x, μ.α, μ.β)

jacobi_α(μ::JacobiWeight) = μ.α
jacobi_β(μ::JacobiWeight) = μ.β

Base.show(io::IO, μ::JacobiWeight) = print(io, "(1+x)^$(jacobi_α(μ))(1-x)^$(jacobi_β(μ)) dx  (Jacobi)")
Display.object_parentheses(μ::JacobiWeight) = true

"""
The `Chebyshev` or `ChebyshevT` weight is the measure on `[-1,1]` with the
Chebyshev weight function `w(x) = 1/√(1-x^2)`.
"""
struct ChebyshevTWeight{T} <: AbstractJacobiWeight{T}
end
ChebyshevTWeight() = ChebyshevTWeight{Float64}()

const ChebyshevWeight = ChebyshevTWeight

chebyshev_weight_firstkind(x) = 1/sqrt(1-x^2)

jacobi_α(μ::ChebyshevTWeight{T}) where {T} = -one(T)/2
jacobi_β(μ::ChebyshevTWeight{T}) where {T} = -one(T)/2

similar(μ::ChebyshevTWeight, ::Type{T}) where {T <: Real} = ChebyshevTWeight{T}()
support(μ::ChebyshevTWeight{T}) where {T} = ChebyshevInterval{T}()
unsafe_weightfun(μ::ChebyshevTWeight, x) = chebyshev_weight_firstkind(x)

weightfunction(μ::ChebyshevTWeight) = chebyshev_weight_firstkind

Base.show(io::IO, μ::ChebyshevTWeight) = print(io, "1/√(1+x)^2 dx  (ChebyshevT)")
Display.object_parentheses(μ::ChebyshevTWeight) = true

"""
The ChebyshevU weight is the measure on `[-1,1]` with the Chebyshev weight
function of the second kind `w(x) = √(1-x^2).`
"""
struct ChebyshevUWeight{T} <: AbstractJacobiWeight{T}
end
ChebyshevUWeight() = ChebyshevUWeight{Float64}()

chebyshev_weight_secondkind(x) = sqrt(1-x^2)

jacobi_α(μ::ChebyshevUWeight{T}) where {T} = one(T)/2
jacobi_β(μ::ChebyshevUWeight{T}) where {T} = one(T)/2

similar(μ::ChebyshevUWeight, ::Type{T}) where {T <: Real} = ChebyshevUWeight{T}()
support(μ::ChebyshevUWeight{T}) where {T} = ChebyshevInterval{T}()
unsafe_weightfun(μ::ChebyshevUWeight, x) = chebyshev_weight_secondkind(x)

weightfunction(μ::ChebyshevUWeight) = chebyshev_weight_secondkind

Base.show(io::IO, μ::ChebyshevUWeight) = print(io, "√(1+x)^2 dx  (ChebyshevU)")
Display.object_parentheses(μ::ChebyshevUWeight) = true

convert(::Type{JacobiWeight}, μ::LegendreWeight{T}) where {T} =
    JacobiWeight{T}(0, 0)
convert(::Type{JacobiWeight}, μ::ChebyshevTWeight{T}) where {T} =
    JacobiWeight{T}(-one(T)/2, -one(T)/2)
convert(::Type{JacobiWeight}, μ::ChebyshevUWeight{T}) where {T} =
    JacobiWeight{T}(one(T)/2, one(T)/2)

function convert(::Type{ChebyshevTWeight}, μ::JacobiWeight{T}) where {T}
    (jacobi_α(μ) ≈ -one(T)/2 && jacobi_β(μ) ≈ -one(T)/2) || throw(InexactError(:convert, ChebyshevTWeight, μ))
    ChebyshevTWeight{T}()
end

function convert(::Type{ChebyshevUWeight}, μ::JacobiWeight{T}) where {T}
    (jacobi_α(μ) ≈ one(T)/2 && jacobi_β(μ) ≈ one(T)/2) || throw(InexactError(:convert, ChebyshevUWeight, μ))
    ChebyshevUWeight{T}()
end

function convert(::Type{LegendreWeight}, μ::JacobiWeight{T}) where {T}
    (μ.α ≈ 0 && μ.β ≈ 0) || throw(InexactError(:convert, LegendreWeight, μ))
    LegendreWeight{T}()
end

jacobi_α(μ::LegendreWeight{T}) where {T} = zero(T)
jacobi_β(μ::LegendreWeight{T}) where {T} = zero(T)

==(μ1::AbstractJacobiWeight, μ2::LegendreWeight) =
    jacobi_α(μ1) == jacobi_α(μ2) && jacobi_β(μ1) == jacobi_β(μ2)
==(μ1::LegendreWeight, μ2::AbstractJacobiWeight) =
    jacobi_α(μ1) == jacobi_α(μ2) && jacobi_β(μ1) == jacobi_β(μ2)

"The generalised Laguerre measure on the halfline `[0,∞)`."
struct LaguerreWeight{T} <: Weight{T}
    α   ::  T

    LaguerreWeight{T}(α = zero(T)) where {T} = new(α)
end
LaguerreWeight() = LaguerreWeight{Float64}()
LaguerreWeight(α::T) where {T<:AbstractFloat} = LaguerreWeight{T}(α)
LaguerreWeight(α) = LaguerreWeight(float(α))

laguerre_weightfun(x, α) = exp(-x) * x^α

similar(μ::LaguerreWeight, ::Type{T}) where {T <: Real} = LaguerreWeight{T}(μ.α)
support(μ::LaguerreWeight{T}) where {T} = HalfLine{T}()
isnormalized(m::LaguerreWeight) = m.α == 0

unsafe_weightfun(μ::LaguerreWeight, x) = laguerre_weightfun(x, μ.α)

laguerre_α(μ::LaguerreWeight) = μ.α

Base.show(io::IO, μ::LaguerreWeight) =
    laguerre_α(μ) == 0 ? print(io, "exp(-x)dx  (Laguerre)") : print(io, "x^$(laguerre_α(μ))exp(-x)dx  (Laguerre)")
Display.object_parentheses(μ::LaguerreWeight) = true

"The Hermite measure with weight exp(-x^2) on the real line."
struct HermiteWeight{T} <: Weight{T}
end
HermiteWeight() = HermiteWeight{Float64}()

hermite_weightfun(x) = exp(-x^2)

similar(μ::HermiteWeight, ::Type{T}) where {T <: Real} = HermiteWeight{T}()
unsafe_weightfun(μ::HermiteWeight, x) = hermite_weightfun(x)

weightfunction(μ::HermiteWeight) = hermite_weightfun

Base.show(io::IO, μ::HermiteWeight) = print(io, "exp(-x^2)dx  (Hermite)")
Display.object_parentheses(μ::HermiteWeight) = true

"The Gaussian measure with weight exp(-|x|^2/2)."
struct GaussianWeight{T} <: Weight{T}
end
GaussianWeight() = GaussianWeight{Float64}()

gaussian_weightfun(x, ::Type{T} = prectype(x)) where {T} =
    1/(2*convert(T, pi))^(length(x)/2) * exp(-norm(x)^2)

similar(μ::GaussianWeight, ::Type{T}) where {T} = GaussianWeight{T}()
isnormalized(μ::GaussianWeight) = true
unsafe_weightfun(μ::GaussianWeight, x) = gaussian_weightfun(x, prectype(μ))

Base.show(io::IO, μ::GaussianWeight{T}) where {T<:StaticTypes} =
    print(io, "1/(2π)^($(euclideandimension(T)/2)) exp(-|x|^2)dx  (Gauss)")
Base.show(io::IO, μ::GaussianWeight) = print(io, "C(d) exp(-|x|^2)dx  (Gauss)")
Display.object_parentheses(μ::GaussianWeight) = true

"The Lebesgue measure associated with the given domain"
lebesguemeasure(domain::UnitInterval{T}) where {T} = LebesgueUnit{T}()
lebesguemeasure(domain::ChebyshevInterval{T}) where {T} = LegendreWeight{T}()
lebesguemeasure(domain::FullSpace{T}) where {T} = Lebesgue{T}()
lebesguemeasure(domain::Domain{T}) where {T} = LebesgueDomain{T}(domain)

dx(domain::Domain) = lebesguemeasure(domain)


text:
"""
PoroTrack(args...)  -> PyPlot Plot

Return a plot of Density-Neutron logs in a single track.

<br> The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Depth|Mandatory| -- | Array{Number,1}| 1D Array of Depths|
|Density|Mandatory| -- | Array{Number,1}|1D Array of Density Values|
|Neutron|Mandatory| -- | Array{Number,1}|1D Array of Neutron Values|
|DepthFrom|Optional| minimum(Depth) | Number |Set the Depth upper limit to plot|
|DepthTo|Optional| maximum(Depth) | Number |  Set the Depth lower limit to plot|
|denax_color|Optional| "darkred" |denax_color="color" |  Set the Density curve color
|ntr_color|Optional| "darkgreen" |ntr_color="color" |  Set the Neutron curve Color|
|Dtick|Optional|true|Dtick=Bool|Set the Depth ticks On/Off|
|LimeSync|Optionalfalse|LimeSync=true|Sync the plot for Limestone|
|Fill|Optional|true|Fill=Bool| Fill when crossed Neutron-Density|
"""
function PoroTrack(Depth,Density,Neutron;DepthFrom=minimum(Depth),DepthTo=maximum(Depth),
                    denax_color="darkred",ntr_color="darkblue",Dtick=false,LimeSync=false, Fill=true)

 if LimeSync==true
        d=2.71
    else
        d=2.65
    end

    m=(d-1.9)/(0-0.45)
    b=-m*0.45+1.9
    denlim=-0.15*m.+b

denax=gca()
denax.plot(Density,Depth,color=denax_color)
denax.set_ylim([DepthTo,DepthFrom])
denax.set_xlabel("Density")
xticks(range(1.9,stop=denlim,length=4))
denax.set_xlim([1.9,denlim])
denax.tick_params("both",labelsize=6)
denax.grid(true)
denax.set_yticks(range(DepthFrom, stop=DepthTo, length=51),minor=true)
denax.set_yticks(range(DepthFrom, stop=DepthTo, length=11))

ntrax=denax.twiny()
ntrax.plot(Neutron,Depth,color=ntr_color)
xlabel("Neutron")
xticks(range(0.45,stop=-0.15,length=4))
xlim(0.45,-0.15)
ntrax.tick_params("both",labelsize=6)
ntrax.grid(true)
setp(denax.get_yticklabels(),visible=Dtick)
 #Convert the Neutron values to Density Values in order to fill the cross Density-Neutron
#When the track is callibrated for sandstone use m=-1.666667 and b=2.65
#When the track is callibrated for limestone use m=-1.8 and b=2.71

NtrToDen=Neutron.*m.+b
if Fill==true
denax.fill_betweenx(Depth,Density,NtrToDen,where=(Density .< NtrToDen),color="red")
end
    end


text:
using DiffGPFlux
using Measurements
using Plots
using Random: randperm
using Printf

# build training and testing data,
# for X array along the 1st dimension is the feature,
# and along the 2nd dimension is the data samples,
# y array is required to be a Vector.
X = reshape(Array(range(0,2π,length=100)), 1, 100)
y = dropdims(sin.(X) + 0.01 .* randn(1, 100), dims=1)
∂y = dropdims(cos.(X) + 0.01 .* randn(1, 100), dims=1)

index = randperm(100)
train_index = index[1:20]
test_index = index[21:end]
X_train, y_train, ∂y_train = X[:, train_index], y[train_index], ∂y[train_index]
X_test, y_test, ∂y_test = X[:, test_index], y[test_index], ∂y[test_index]

γ0 = seed_duals([0.5], Float64, 1)
β = 0.1
GP = GaussProcess(γ0, β, rbf, ∇ₓrbf)
θ = [GP.γ]
indexes = [[1]]

function loss(gp, X̄, ȳ, ∂ȳ, x)
    ∂ₓf = x -> ∇ₓpredict(gp, X̄, ȳ, x)
    ∂ₓy = dropdims(∂ₓf(X̄), dims=1)
    return mean((∂ₓy-∂ȳ).^2)
end

opt = ADAM(0.01)

plt = plot(1, marker=2)
@gif for i in 1:300
    L = loss(GP, X_train, y_train, ∂y_train, X_train)
    if i%25 == 1
        @printf("Step %d, train error %.4f, test error %.4f\n",
                i, L.value,
                (loss(GP, X_test, y_test, ∂y_test, X_test)).value)
    end
    push!(plt, L.value)
    θ̄_array = [g for g in L.partials]
    θ̄ = grads(θ̄_array, indexes)
    update!(opt, θ, θ̄, 1)
    end every 10

@show(GP.γ)

γ_values = [GP.γ[1].value]
non_dual_GP = GaussProcess(γ_values, β, rbf, ∇ₓrbf)
pred_y, pred_σ = predict(non_dual_GP, X_train, y_train, X)
pred_∂y = ∇ₓpredict(non_dual_GP, X_train, y_train, X)
plot(X', pred_y.±pred_σ, color=:blue)
plot!(X', y, seriestype=:scatter, color=:red)

plot(X', pred_∂y', color=:green)
plot!(X', ∂y, seriestype=:scatter, color=:red)

max_stars_repo_name: HamletWantToCode/DiffGPFlux

max_issues_repo_name: HamletWantToCode/DiffGPFlux


text:
module MathLink

export mlinitialize, mldeinitialize, mlopen, mlclose, mlactivate,
mlerror, mlerrormessage, mlclearerror, mlflush, mlready,
mlnewpacket, mlendpacket, mlnextpacket,
mlgetnext, mlgettype, mlputnext, mlgetnextraw, mlgetrawtype,
mlput, mlget, MLFunction, mltesthead, @mlput,
mlloopbackopen, mltransferexpression, mltransfertoendofloopbacklink,
mlcreatemark, mlseektomark, mldestroymark,
mlgetlinkedenvidstring, mlsetenvidstring, mllinkname, mltolinkid, mlfromlinkid

import Base.convert

include("setup.jl")
include("typesconsts.jl")
include("mlfn.jl")
include("mlextra.jl")
include("macro.jl")

# module
end


text:
using FileCmp
using FileCmp: Info, files_equal, got_eof, bytes_read
using Test

@testset "FileCmp.jl" begin
    write_file = function(fname::String, data)
        open(fname, "w") do io
            write(io, data)
        end
    end
    dir = mktempdir()
    nbytes = 1000
    nbytes_trunc = 500
    v1 = rand(UInt8, nbytes)
    path = fname -> joinpath(dir, fname)
    write_file(path("test1"), v1)
    write_file(path("test2"), @view v1[1:nbytes_trunc])
    v2 = rand(UInt8, nbytes)
    write_file(path("test3"), v2)
    cp(path("test1"), path("test4"))
    v3 = rand(UInt8, nbytes_trunc)
    write_file(path("test5"), v3)

    @test filecmp(path("test1"), path("test1"))
    _info = filecmp(path("test1"), path("test1"); info=true)
    @test files_equal(_info)
    @test got_eof(_info) == 0
    @test bytes_read(_info) == 0  # files are the same
    @test ! filecmp(path("test1"), path("test2"))

    _info = filecmp(path("test1"), path("test2"); info=true)
    @test isa(_info, Info)
    @test _info._byte_index == 0
    @test _info._size_cmp == 1
    @test got_eof(_info) == 1
    @test bytes_read(_info) == nbytes_trunc
    @test filecmp(path("test1"), path("test2"); limit=filesize(path("test2")))
    @test ! filecmp(path("test1"), path("test2"); limit=(filesize(path("test2")) + 1))

    @test ! filecmp(path("test1"), path("test3"))

    _info = filecmp(path("test1"), path("test3"); info=true)
    @test isa(_info, Info)
    @test _info._byte_index != 0
    @test _info._size_cmp == 0
    @test ! files_equal(_info)
    @test got_eof(_info) == 0
    @test filecmp(path("test1"), path("test4"))
    @test files_equal(filecmp(path("test1"), path("test4"); info=true))
    @test ! filecmp(path("test1"), path("test5"))

    _info = filecmp(path("test1"), path("test5"); info=true)
    @test isa(_info, Info)
    @test _info._byte_index != 0
    @test _info._size_cmp == 1
    @test ! files_equal(_info)
    @test got_eof(_info) == 0

    @test_throws SystemError filecmp("sdfsdf", "sdfdd")
    @test_throws SystemError filecmp("sdfsdf", path("test1"))
    @test_throws SystemError filecmp(path("test1"), "sdfd")
end


text:
using Compose
using Test

include("misc.jl")
@testset "SVG Correctness Tests" begin
    include("svg.jl")
end
include("immerse.jl")

# Run the examples
cd(joinpath(@__DIR__, "output"))

exampledir = joinpath(@__DIR__, "examples")
for ex in readdir(exampledir)
    endswith(ex, ".jl") || continue
    file = joinpath(exampledir, ex)
    expr = quote
        module $(Symbol(replace(ex, "."=>"_")))
        using Compose
        using Colors
        using Random
        Random.seed!(1) #Needed so that SVG uuid is reproducible
        include($file)
        end
    end
    expr.head = :toplevel
    eval(expr)
end


text:
# -----------------------------------------------
# Generalized Linear Models routines
# -----------------------------------------------

function Base.show(io::IO, fitresult::AbstractFittingResult)
    # println(io, "£$(money.pounds).$(money.shillings)s.$(money.pence)d")
    f_ = "-------------"
    #s_ = "  "
    println("\n# $f_ RegResult [ $(typeof(fitresult.callertype)) ] $f_ #\n")
    #println(" Beta: $(fitresult.beta)\n")
    #present(fitresult.beta)
    println(" N of obs.: \t$(@sprintf("%8d", fitresult.N))")
    println(" MSE:       \t$(@sprintf("%8f", fitresult.mse))\n")
    println("$f_$f_$f_$f_$f_$f_$f_")
    println("     Y     |      β         σ [S.E.]       t          P > |t|      [95% Conf. Interval] ")
    println("$f_$f_$f_$f_$f_$f_$f_")
    for i in 1:length(fitresult.beta)
        print("    X[$i]   |")
        print(" $(@sprintf("%10f", fitresult.beta[i]))  ")
        print(" $(@sprintf("%10f", sqrt(fitresult.variance[i,i])))  ")
        print(" $(@sprintf("%8.3f", fitresult.tstat[i]))    ")
        print(" $(@sprintf("%8.3f", fitresult.tsig[i]))      ")
        print(" $(@sprintf("%8f", fitresult.confint[i][1]))   $(@sprintf("%8f", fitresult.confint[i][2]))")
        println("")
    end
    println("$f_$f_$f_$f_$f_$f_$f_")
    #print(" Variance: ")
    #present(fitresult.variance)
    try
        robust = fitresult.robust
        println("\n Robustness: $(robust)")
    catch end
    try
        converged = fitresult.converged
        println("\n Converged: $(converged)")
    catch end
    try
        iterations = fitresult.iterations
        if iterations > 0 println("\n [ Iterations: $(iterations) ]") end
    catch end
    println("\n# $f_ EndResult [ $(typeof(fitresult.callertype)) ] $f_ #")
end

function Base.show(io::IO, fitresult::glmfitresult)
    # println(io, "£$(money.pounds).$(money.shillings)s.$(money.pence)d")
    f_ = "-------------"
    #s_ = "  "
    println("\n# $f_ RegResult [ $(typeof(fitresult.callertype)) ] $f_ #\n")
    #println(" Beta: $(fitresult.beta)\n")
    #present(fitresult.beta)
    println(" N of obs.: \t$(@sprintf("%8d", fitresult.N))")
    if typeof(fitresult.callertype) == Logit
        println(" D:        \t$(@sprintf("%8f", fitresult.fit))\n")
    else
        println(" MSE:       \t$(@sprintf("%8f", fitresult.fit))\n")
    end
    println("$f_$f_$f_$f_$f_$f_$f_")
    println("     Y     |      β         σ [S.E.]       t          P > |t|      [95% Conf. Interval] ")
    println("$f_$f_$f_$f_$f_$f_$f_")
    for i in 1:length(fitresult.beta)
        print("    X[$i]   |")
        print(" $(@sprintf("%10f", fitresult.beta[i]))  ")
        print(" $(@sprintf("%10f", sqrt(fitresult.variance[i,i])))  ")
        print(" $(@sprintf("%8.3f", fitresult.tstat[i]))    ")
        print(" $(@sprintf("%8.3f", fitresult.tsig[i]))      ")
        print(" $(@sprintf("%8f", fitresult.confint[i][1]))   $(@sprintf("%8f", fitresult.confint[i][2]))")
        println("")
    end
    println("$f_$f_$f_$f_$f_$f_$f_")
    #print(" Variance: ")
    #present(fitresult.variance)
    try
        robust = fitresult.robust
        println("\n Robustness: $(robust)")
    catch end
    try
        converged = fitresult.converged
        println("\n Converged: $(converged)")
    catch end
    try
        iterations = fitresult.iterations
        if iterations > 0 println("\n [ Iterations: $(iterations) ]") end
    catch end
    println("\n# $f_ EndResult [ $(typeof(fitresult.callertype)) ] $f_ #")
end


text:
##############################################################################
#
# Test problem
#
#############################################################################

function rosenbrock(x)

    n = 6;

    # Initializations
    f = 0

    evenIdx = 2:2:n
    oddIdx  = 1:2:(n-1)

    f1  = x[evenIdx] .- x[oddIdx].^2
    f2  = 1 .- x[oddIdx]

    # Function
    f   = sum( f1.^2 .+ f2.^2 )

    return f
end


text:
module TheNumberLine

using Luxor

export NumberLinePlot, NumberLineExpression
export aSlider, markers, arrows

include("main_elements.jl")
include("standard_number_line.jl")

end


text:
using Statistics

# redondeamos
function Muestra(v)
    vr = map(x->round(x, digits=3),v)

    d= 13
    cnt = 1
    l = length(v)
    while(cnt+d < l)
        println("& $(vr[cnt:cnt+d]) \\\\")
        cnt += d
    end
    println("& $(vr[cnt:end])")
    println("Media de valores = $(mean(vr))")
    println("desviación típica = $(std(vr))")
end

L = [
    [ 84.28571428571429	0.0	42.142857142857146	63.054625 ],
    [ 80.0	0.0	40.0	55.254541999999994 ],
    [  80.0	0.0	40.0	44.4975 ],
    [  90.0	0.0	45.0	62.065083 ],
    [  85.5072463768116	0.0	42.7536231884058	55.16675 ],
    [  83.95859213250517	0.0	41.97929606625259	56.0077],
    [  4.193591874396594	0.0	2.096795937198297	7.418055484420731],
]

A =[
    [0.5887686685732444, 0.9531744518589133, 0.5866821215009963, 0.782117254528122, 0.5623068027317195, 0.6564704480457578, 0.40674526334604383, 0.6092219347224108, 0.4552812202097233, 0.6414910458940486, 0.41585795996186825, 0.5132057742067271, 0.5858436606291706, 0.6432692307692305, 0.6143725947110118, 0.5923336234647213, 0.45750931623190894, 0.48422240758394236, 0.6892960616086276, 0.5370658067206863, 0.40613346043851245, 0.526081775604095, 0.5556662151087918, 0.6505324028274908, 0.821819643441515, 0.8320472522093112, 0.7271819722487023, 0.6688067310482029, 0.7362589190282229, 0.8981491427113681, 0.6142272547606314, 0.3433301557038448, 0.5020524892511816, 0.5551902109151572, 0.5174156010650538, 0.5355033067206862, 0.697372370849644, 0.5904585163074899, 0.5369783391589873, 0.7698220527486493, 0.5924077352580687, 0.7089259612329201, 0.7378260986318287, 0.9025683535734513],
    [0.03908944540742903, 0.10470510872715894, 0.009605606600198751, 0.08750347902131649, 0.000432002224942483, 0.03444562533261263, 0.016393063448458048, 0.0645436367671232, 0.07204338176315458, 0.09040109104441989, 0.09431232580640499, 0.030365669146988485, 0.08458274193229523, 0.19941856241644776, 0.02117555900999259, 0.008230673946275514, 0.08049213538408105, 0.06439523762061847, 0.09463394198220276, 0.0590562452596907, 0.05284602408406678, 0.011556501072216433, 0.1574868850767859, 0.10903078639232158, 0.1794055372461562, 0.12248525685757713, 0.21406928452373764, 0.1288926474824063, 0.013741269778967032, 0.04075031910619335, 0.1331650773668219, 0.01864100520804161, 0.0866785752255467, 0.018528845429171395, 0.06707799659994318, 0.05818278120597983, 0.09091750441647033, 0.001632348144195952, 0.04774751790610398, 0.02834875186917607, 0.060173587813287854, 0.03672956591258585, 0.11973504684026441, 0.023423090873852793]
    ]
for l in A
    Muestra(l)
    println("----")
end


text:
# abstract data type for CCDData
"""
    AbstractCCDData{T}

Supertype for `CCDData` based on `AbstractMatrix` interface.

Every subtype of this type should have methods `data` and `hdr` defined.
"""
abstract type AbstractCCDData{T} <: AbstractMatrix{T} end

# extending the AbstractMatrix interface
Base.size(ccd::AbstractCCDData) = size(data(ccd))
Base.getindex(ccd::AbstractCCDData, inds...) = getindex(data(ccd), inds...) # default fallback for operations on Array
Base.setindex!(ccd::AbstractCCDData, v, inds...) = setindex!(data(ccd), v, inds...) # default fallback for operations on Array
Base.getindex(ccd::AbstractCCDData, inds::AbstractString...) = getindex(hdr(ccd), inds...) # accesing header using strings
Base.setindex!(ccd::AbstractCCDData, v, inds::AbstractString...) = setindex!(hdr(ccd), v, inds...) # modifying header using strings
Base.getindex(ccd::AbstractCCDData, inds::Symbol...) = getindex(ccd, string.(inds)...) # accessing header using symbol
Base.setindex!(ccd::AbstractCCDData, v, inds::Symbol...) = setindex!(ccd, v, string.(inds)...) # modifying header using Symbol
Base.promote_rule(::Type{AbstractCCDData{T}}, ::Type{AbstractCCDData{V}}) where {T,V} = AbstractCCDData{promote_type{T,V}}

# custom data type to hold ImageHDU
"""
    CCDData(data::AbstractMatrix, [hdr::FITSHeader])

Struct to store `ImageHDU`, derived from [`AbstractCCDData`](@ref).

`CCDData` acts like a matrix with a header associated.

```julia
ccd = CCDData(zeros(4, 4))

ccd[1]
```
This accesses the 1st element in matrix associated with `ccd`.

```
ccd["SIMPLE"]
```
One can also access the header directly from `ccd`, the key can be `Symbol` as well.

```
ccd[:SIMPLE] = false
```
Header values can be directly modified from `ccd`.

One can perform arithmetic operations on it as well:

```julia
ccd1 = CCDData(zeros(4, 4))

ccd2 = CCDData(ones(4, 4))

sum_ccd1 = ccd1 + ccd2

sum_ccd2 = ccd2 + ccd1
```

`sum_ccd1` has the header of `ccd1` whereas `sum_ccd2` has the header `ccd2`.

If header is not provided in the `CCDData` constructor, `default_header` is used to generate the header.
"""
struct CCDData{T,M<:AbstractMatrix{T}} <: AbstractCCDData{T}
    data::M
    hdr::FITSHeader
end

"""
    CCDData(hdu::ImageHDU)

Loads `ImageHDU` as `CCDData`.
"""
CCDData(hdu::ImageHDU) = CCDData(getdata(hdu), read_header(hdu))

CCDData(data) = CCDData(data, default_header(data))

"""
    CCDData(path::AbstractString; hdu = 1)

Loads HDU from `hdu` index in `FITS` file at `path` as `CCDData`.
"""
CCDData(path::AbstractString; hdu = 1) = FITS(f -> CCDData(f[hdu]), path)

CCDData{T}(data, hdr) where {T<:Number} = CCDData(T.(data), hdr)
data(ccd::CCDData) = ccd.data
hdr(ccd::CCDData) = ccd.hdr
Base.copy(ccd::CCDData) = CCDData(copy(data(ccd)), deepcopy(hdr(ccd)))
Base.convert(::Type{CCDData{T}}, ccd::CCDData{V}) where {T,V} = CCDData{T}(data(ccd), hdr(ccd))
Base.view(ccd::CCDData, inds...) = CCDData(view(data(ccd), inds...), hdr(ccd))
Base.selectdim(ccd::CCDData, d::Integer, idxs) = CCDData(selectdim(data(ccd), d, idxs), hdr(ccd))
# broadcast mechanics
Base.BroadcastStyle(::Type{<:CCDData}) = Broadcast.ArrayStyle{CCDData}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{CCDData}}, ::Type{T}) where T
    ccd = find_ccd(bc)
    CCDData(similar(data(ccd), T, axes(bc)), hdr(ccd))
end
"`A = find_ccd(As)` returns the first CCDData among the arguments."
find_ccd(bc::Base.Broadcast.Broadcasted) = find_ccd(bc.args)
find_ccd(args::Tuple) = find_ccd(find_ccd(args[1]), Base.tail(args))
find_ccd(x) = x
find_ccd(::Tuple{}) = nothing
find_ccd(a::CCDData, rest) = a
find_ccd(::Any, rest) = find_ccd(rest)


text:
# ---
# title: 637. Average of Levels in Binary Tree
# id: problem637
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/average-of-levels-in-binary-tree/description/>
# hidden: true
# ---
#
# Given a non-empty binary tree, return the average value of the nodes on each
# level in the form of an array.
#
# **Example 1:**
#
#
#
#     Input:
#         3
#        / \
#       9  20
#         /  \
#        15   7
#     Output: [3, 14.5, 11]
#     Explanation:
#     The average value of nodes on level 0 is 3,  on level 1 is 14.5, and on level 2 is 11. Hence return [3, 14.5, 11].
#
#
# **Note:**
#
#   1. The range of node's value is in the range of 32-bit signed integer.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export SetProperty, CLOSED, OPEN
export AbstractVertexSet, Neighbors, Mirrors, Vertex, UnionOf
export NeighborCover, MirrorCover, NearestNeighbors
export vertices, neighbors, neighborcover

# Vertex Set Notations
abstract type AbstractVertexSet end

abstract type SetProperty end
struct CLOSED <: SetProperty end
struct OPEN <: SetProperty end

"""
    Neighbors{SP<:SetProperty, NTH} <: AbstractVertexSet

Neighbors, if `SP` is `CLOSED`, then this set contains vertex itself.
"""
struct Neighbors{SP<:SetProperty, NTH} <: AbstractVertexSet
    i::Int
end

const NearestNeighbors{SP} = Neighbors{SP, 1}
const NeighborCover = Neighbors{CLOSED, 1}

"""
    Mirrors{SP<:SetProperty} <: AbstractVertexSet

Mirrors and itself. A vertex `w` is a mirror of `v` if `w ∈ N²(v)` and `N[v]\\N(w)` is a clique.
"""
struct Mirrors{SP<:SetProperty} <: AbstractVertexSet
    i::Int
end
const MirrorCover{SP} = Mirrors{CLOSED}

struct Vertex<:AbstractVertexSet
    i::Int
end

"""
    UnionOf{TA, TB}<:AbstractVertexSet

A union of ...
"""
struct UnionOf{TA, TB}<:AbstractVertexSet
    A::TA
    B::TB
end
Base.:∪(A::AbstractVertexSet, B::AbstractVertexSet) = UnionOf(A, B)

"""
    neighborcover(eg::EliminateGraph, i::Int)

Get neighbors of vertex `i`, including itself.
"""
function neighborcover end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()

@testset "NormInfinity" begin
    bridged_mock = MOIB.Constraint.NormInfinity{Float64}(mock)

    MOIT.basic_constraint_tests(bridged_mock, config,
        include = [(F, MOI.NormInfinityCone) for F in [
            MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}, MOI.VectorQuadraticFunction{Float64}
        ]])

    @testset "norminf1test" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
            [1.0, 0.5, 1.0],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0.0, 1.0, 0.0, 0.0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [-1]])

        MOIT.norminf1vtest(bridged_mock, config)
        MOIT.norminf1ftest(bridged_mock, config)

        var_names = ["x", "y", "z"]

        @testset "Test mock model" begin
            MOI.set(mock, MOI.VariableName(), MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
            nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
            @test length(nonneg) == 1
            MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")
            zeros = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Zeros}())
            @test length(zeros) == 2
            MOI.set(mock, MOI.ConstraintName(), zeros[1], "x_eq")
            MOI.set(mock, MOI.ConstraintName(), zeros[2], "y_eq")

            s = """
            variables: x, y, z
            nonneg: [x + -1.0y, x + -1.0z, x + y, x + z] in MathOptInterface.Nonnegatives(4)
            x_eq: [-1.0 + x] in MathOptInterface.Zeros(1)
            y_eq: [-0.5 + y] in MathOptInterface.Zeros(1)
            maxobjective: y + z
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(mock, model, var_names, ["nonneg", "x_eq", "y_eq"])
        end

        @testset "Test bridged model" begin
            MOI.set(bridged_mock, MOI.VariableName(), MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
            norminf = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormInfinityCone}())
            @test length(norminf) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), norminf[1], "norminf")
            zeros = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Zeros}())
            @test length(zeros) == 2
            MOI.set(bridged_mock, MOI.ConstraintName(), zeros[1], "x_eq")
            MOI.set(bridged_mock, MOI.ConstraintName(), zeros[2], "y_eq")

            s = """
            variables: x, y, z
            norminf: [1.0x, y, z] in MathOptInterface.NormInfinityCone(3)
            x_eq: [-1.0 + x] in MathOptInterface.Zeros(1)
            y_eq: [-0.5 + y] in MathOptInterface.Zeros(1)
            maxobjective: y + z
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(bridged_mock, model, var_names, ["norminf", "x_eq", "y_eq"])
        end

        ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormInfinityCone}()))
        nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())

        @testset "$attr" for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
            @test MOI.supports(bridged_mock, attr, typeof(ci))
            value = [4.0, 1.0, -2.0]
            MOI.set(bridged_mock, attr, ci, value)
            @test MOI.get(bridged_mock, attr, ci) ≈ value
            if attr isa MOI.ConstraintPrimalStart
                nonneg_value = Float64[3, 6, 5, 2]
            else
                nonneg_value = [0.25, 2.25, 1.25, 0.25]
            end
            @test MOI.get(mock, attr, nonneg[1]) ≈ nonneg_value
        end

        test_delete_bridge(bridged_mock, ci, 3, (
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),))
    end

    @testset "norminf3test" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
            [2, -1, -1, -1],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [vcat(fill(inv(3), 3), zeros(3)), fill(inv(3), 3)])
        MOIT.norminf3test(bridged_mock, config)

        var_names = ["x", "y1", "y2", "y3"]

        @testset "Test mock model" begin
            MOI.set(mock, MOI.VariableName(), MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
            nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
            @test length(nonneg) == 2
            MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg1")
            MOI.set(mock, MOI.ConstraintName(), nonneg[2], "nonneg2")

            s = """
            variables: x, y1, y2, y3
            nonneg1: [x + -1.0y1 + -3.0, x + -1.0y2 + -3.0, x + -1.0y3 + -3.0, x + y1 + 1.0, x + y2 + 1.0, x + y3 + 1.0] in MathOptInterface.Nonnegatives(6)
            nonneg2: [y1 + 1.0, y2 + 1.0, y3 + 1.0] in MathOptInterface.Nonnegatives(3)
            minobjective: x
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(mock, model, var_names, ["nonneg1", "nonneg2"])
        end

        @testset "Test bridged model" begin
            MOI.set(bridged_mock, MOI.VariableName(), MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
            norminf = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormInfinityCone}())
            @test length(norminf) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), norminf[1], "norminf")
            nonneg = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
            @test length(nonneg) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), nonneg[1], "nonneg")

            s = """
            variables: x, y1, y2, y3
            norminf: [x + -1.0, y1 + 2.0, y2 + 2.0, y3 + 2.0] in MathOptInterface.NormInfinityCone(4)
            nonneg: [y1 + 1.0, y2 + 1.0, y3 + 1.0] in MathOptInterface.Nonnegatives(3)
            minobjective: x
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(bridged_mock, model, var_names, ["norminf", "nonneg"])
        end

        ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormInfinityCone}()))
        nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())

        @testset "$attr" for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
            @test MOI.supports(bridged_mock, attr, typeof(ci))
            value = Float64[5, 1, -2, 3]
            MOI.set(bridged_mock, attr, ci, value)
            @test MOI.get(bridged_mock, attr, ci) ≈ value
            if attr isa MOI.ConstraintPrimalStart
                nonneg_value = Float64[4, 7, 2, 6, 3, 8]
            else
                nonneg_value = [-1, 11, -1, 5, -1, 17] / 6
            end
            @test MOI.get(mock, attr, nonneg[1]) ≈ nonneg_value
        end

        test_delete_bridge(bridged_mock, ci, 4, (
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 1),))
    end
end

@testset "NormOne" begin
    bridged_mock = MOIB.Constraint.NormOne{Float64}(mock)

    MOIT.basic_constraint_tests(bridged_mock, config,
        include = [(F, MOI.NormOneCone) for F in [
            MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}, MOI.VectorQuadraticFunction{Float64}
        ]])

    @testset "norminf1test" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
            [1.0, 0.5, 0.5, 0.5, 0.5],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1.0, 1.0, 1.0, 0.0, 0.0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [0]])

        MOIT.normone1vtest(bridged_mock, config)
        MOIT.normone1ftest(bridged_mock, config)

        var_names = ["x", "y", "z"]
        MOI.set(bridged_mock, MOI.VariableName(), MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)

        nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
        u, v = MOI.get(mock, MOI.ListOfVariableIndices())[4:5]
        @testset "Test mock model" begin
            MOI.set(mock, MOI.VariableName(), u, "u")
            MOI.set(mock, MOI.VariableName(), v, "v")
            @test length(nonneg) == 1
            MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")
            zeros = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Zeros}())
            @test length(zeros) == 2
            MOI.set(mock, MOI.ConstraintName(), zeros[1], "x_eq")
            MOI.set(mock, MOI.ConstraintName(), zeros[2], "y_eq")

            s = """
            variables: x, y, z, u, v
            nonneg: [x + -1.0u + -1.0v, u + -1.0y, v + -1.0z, u + y, v + z] in MathOptInterface.Nonnegatives(5)
            x_eq: [-1.0 + x] in MathOptInterface.Zeros(1)
            y_eq: [-0.5 + y] in MathOptInterface.Zeros(1)
            maxobjective: y + z
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(mock, model, [var_names; "u"; "v"], ["nonneg", "x_eq", "y_eq"])
        end

        @testset "Test bridged model" begin
            normone = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormOneCone}())
            @test length(normone) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), normone[1], "normone")
            zeros = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Zeros}())
            @test length(zeros) == 2
            MOI.set(bridged_mock, MOI.ConstraintName(), zeros[1], "x_eq")
            MOI.set(bridged_mock, MOI.ConstraintName(), zeros[2], "y_eq")

            s = """
            variables: x, y, z
            normone: [1.0x, y, z] in MathOptInterface.NormOneCone(3)
            x_eq: [-1.0 + x] in MathOptInterface.Zeros(1)
            y_eq: [-0.5 + y] in MathOptInterface.Zeros(1)
            maxobjective: y + z
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(bridged_mock, model, var_names, ["normone", "x_eq", "y_eq"])
        end

        ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormOneCone}()))

        @testset "$attr" for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
            @test MOI.supports(bridged_mock, attr, typeof(ci))
            value = [4.0, 1.0, -2.0]
            MOI.set(bridged_mock, attr, ci, value)
            @test MOI.get(bridged_mock, attr, ci) ≈ value
            if attr isa MOI.ConstraintPrimalStart
                @test MOI.get(mock, MOI.VariablePrimalStart(), u) == 1
                @test MOI.get(mock, MOI.VariablePrimalStart(), v) == 2
                @test MOI.get(mock, attr, nonneg[1]) == Float64[1, 0, 4, 2, 0]
            else
                @test MOI.get(mock, attr, nonneg[1]) == Float64[4, 0, 2, 1, 0]
            end
        end

        test_delete_bridge(bridged_mock, ci, 3, (
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),))
    end

    @testset "normone3test" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
            [4, -1, -1, -1, 1, 1, 1],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [vcat(ones(4), zeros(3)), ones(3)])
        MOIT.normone3test(bridged_mock, config)

        var_names = ["x", "y1", "y2", "y3"]

        @testset "Test mock model" begin
            var_names_all = vcat(var_names, "z1", "z2", "z3")
            MOI.set(mock, MOI.VariableName(), MOI.get(mock, MOI.ListOfVariableIndices()), var_names_all)
            nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
            @test length(nonneg) == 2
            MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg1")
            MOI.set(mock, MOI.ConstraintName(), nonneg[2], "nonneg2")

            s = """
            variables: x, y1, y2, y3, z1, z2, z3
            nonneg1: [x + -1.0 + -1.0z1 + -1.0z2 + -1.0z3, z1 + -1.0y1 + -2.0, z2 + -1.0y2 + -2.0, z3 + -1.0y3 + -2.0, z1 + y1 + 2.0, z2 + y2 + 2.0, z3 + y3 + 2.0] in MathOptInterface.Nonnegatives(7)
            nonneg2: [y1 + 1.0, y2 + 1.0, y3 + 1.0] in MathOptInterface.Nonnegatives(3)
            minobjective: x
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(mock, model, var_names_all, ["nonneg1", "nonneg2"])
        end

        @testset "Test bridged model" begin
            MOI.set(bridged_mock, MOI.VariableName(), MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
            normone = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormOneCone}())
            @test length(normone) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), normone[1], "normone")
            nonneg = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
            @test length(nonneg) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), nonneg[1], "nonneg")

            s = """
            variables: x, y1, y2, y3
            normone: [x + -1.0, y1 + 2.0, y2 + 2.0, y3 + 2.0] in MathOptInterface.NormOneCone(4)
            nonneg: [y1 + 1.0, y2 + 1.0, y3 + 1.0] in MathOptInterface.Nonnegatives(3)
            minobjective: x
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(bridged_mock, model, var_names, ["normone", "nonneg"])
        end

        ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormOneCone}()))
        nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())

        @testset "$attr" for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
            @test MOI.supports(bridged_mock, attr, typeof(ci))
            value = (attr isa MOI.ConstraintPrimalStart ? vcat(3, ones(3)) : vcat(1, fill(-1, 3)))
            MOI.set(bridged_mock, attr, ci, value)
            @test MOI.get(bridged_mock, attr, ci) ≈ value
            nonneg_value = (attr isa MOI.ConstraintPrimalStart ? vcat(zeros(4), fill(2.0, 3)) : vcat(ones(4), zeros(3)))
            @test MOI.get(mock, attr, nonneg[1]) ≈ nonneg_value
        end

        test_delete_bridge(bridged_mock, ci, 4, (
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 1),))
    end
end