text:
using OrdinaryDiffEq, ODE, ODEInterfaceDiffEq, LSODA, Sundials, DiffEqDevTools

f = (du,u,p,t) -> begin
  @inbounds begin
  x = view(u,1:7)   # x
  y = view(u,8:14)  # y
  v = view(u,15:21) # x′
  w = view(u,22:28) # y′
  du[1:7] .= v
  du[8:14].= w
  for i in 15:28
    du[i] = zero(u[1])
  end
  for i=1:7,j=1:7
    if i != j
      r = ((x[i]-x[j])^2 + (y[i] - y[j])^2)^(3/2)
      du[14+i] += j*(x[j] - x[i])/r
      du[21+i] += j*(y[j] - y[i])/r
    end
  end
  end
end

prob = ODEProblem(f,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],(0.0,3.0))

abstols = 1.0 ./ 10.0 .^ (6:9)
reltols = 1.0 ./ 10.0 .^ (3:6);

using Plots; gr()

sol = solve(prob,Vern8(),abstol=1/10^12,reltol=1/10^10,maxiters=1000000)
test_sol = TestSolution(sol);
plot(sol)

#setups = [Dict(:alg=>ode45())]
#wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=10000)
#plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>dopri5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=10000)
plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100,maxiters=10000,error_estimate=:L2,dense_errors=true)
plot(wp)

#setups = [Dict(:alg=>ode78())]
#wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
#plot(wp)

setups = [Dict(:alg=>DP8())
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>dop853())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)

setups = [Dict(:alg=>odex())
          Dict(:alg=>Vern7())
          Dict(:alg=>CVODE_Adams())
          Dict(:alg=>lsoda())
          Dict(:alg=>Vern6())
          Dict(:alg=>Tsit5())
          Dict(:alg=>ddeabm())
          Dict(:alg=>ARKODE(Sundials.Explicit(),order=6))
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=20)
plot(wp)

setups = [Dict(:alg=>DP8())
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100,maxiters=1000,error_estimate=:L2,dense_errors=true)
plot(wp)

setups = [Dict(:alg=>Tsit5())
          Dict(:alg=>Vern9())
          Dict(:alg=>VCABM())
          Dict(:alg=>AitkenNeville(min_order=1, max_order=9, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointDeuflhard(min_order=1, max_order=9, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))]
solnames = ["Tsit5","Vern9","VCABM","AitkenNeville","Midpoint Deuflhard","Midpoint Hairer Wanner"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>ExtrapolationMidpointDeuflhard(min_order=1, max_order=9, init_order=9, threading=false))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=false))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, sequence = :romberg, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, sequence = :bulirsch, threading=true))]
solnames = ["Deuflhard","No threads","standard","Romberg","Bulirsch"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=5, max_order=11, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=15, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=5, max_order=7, init_order=6, threading=true))]
solnames = ["1","2","3","4","5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using Documenter
using Bootstrap

makedocs(
    modules = [Bootstrap],
    format = :html,
    sitename = "Bootstrap.jl",
    doctest = true,
    checkdocs = :exports,
    linkcheck = true,
    html_disable_git = true, # disable source and edit links to github
    html_canonical = "https://juliangehring.github.io/Bootstrap.jl/stable/",
    pages = [
        "Home" => "index.md",
        "Library" => "library.md",
        "Changelog" => "NEWS.md",
        "License" => "LICENSE.md"
    ]
)

deploydocs(
    repo = "github.com/juliangehring/Bootstrap.jl.git",
    julia = "1.0",
    target = "build",
    deps   = nothing,
    make   = nothing
)

ext: jl


text:
cdf_inters = []

struct GeneralizedKolmogorov #<: Distribution{Univariate, Continuous}
    β::Float64 #value in the range (0,1]
end

length(d::GeneralizedKolmogorov) = 1

struct GeneralizedBrownianBridge #<: Distribution{Univariate, Continuous}
    a::Float64
    T::Float64
end

length(d::GeneralizedBrownianBridge) = 1

"""
#cdf of the supremum over [0,T] of absolute value of generalized brownian bridge with B(T) = a
"""
function cdf(d::GeneralizedBrownianBridge, x::Float64)
    x ≤ abs(d.a) && return 0.0
    total, j = 1.0, 1
    ek(k) = exp(-2k*x*(k*x-d.a)/d.T)
    while true
        summand = ek(2j) + ek(-2j) - ek(2j-1) - ek(1-2j)
        total += summand
        if abs(summand) < 10^-8
            return total
        end
        if total < 10^-6
            return 0.0
        end
        j +=1
    end
end

"""
cdf of the bubble statistic
"""
function cdf(d::GeneralizedKolmogorov, x::Float64; precomputed_interpolation = true)
    β = d.β

    if precomputed_interpolation
        global cdf_inters
        if isempty(cdf_inters)
            @info "Precomputing integrals"
            cdf_inters = make_cdfs() #carry out a pre-compute if first time
        end

        index = round((β-β_min) / β_step, digits = 6) + 1
        if isinteger(index)
            return cdf_inters[Int(index)](x)
        else
            index_A, index_B = floor(Int,index), ceil(Int,index)
            α = index - index_A
            val_A, val_B = cdf_inters[index_A](x), cdf_inters[index_B](x)
            return (1-α)*val_A + α*val_B
        end

    else #Do not use a precomputed interpolation
        if β == 1.0
            return cdf(GeneralizedBrownianBridge(0,1),x)
        else
            integral, _  = quadgk(0.0,sqrt(β/(1-β))*x) do z
                                                            cdf(GeneralizedBrownianBridge(sqrt((1-β)/β)*z, (1-β)/β),x) *
                                                            cdf(GeneralizedBrownianBridge(sqrt((1-β)*β)*z, β),x) *
                                                            pdf(Normal(),z)
                                                        end
            return 2integral
        end
    end
end

pdf(d::GeneralizedKolmogorov, x::Float64 ; precomputed_interpolation = true) = derivative((x)->cdf(d,x; precomputed_interpolation), x)

quantile(d::GeneralizedKolmogorov, u::Float64 ; precomputed_interpolation = true) = find_zero((x)->(cdf(d, x; precomputed_interpolation)-u), 1/d.β); #initial condition at 1/β seems to work

#QQQQ - problem with broadcasting
cdf(d::GeneralizedKolmogorov, xv::AbstractVector{Float64}) = [cdf(d,x) for x in xv]
pdf(d::GeneralizedKolmogorov, xv::AbstractVector{Float64}) = [pdf(d,x) for x in xv];
quantile(d::GeneralizedKolmogorov, uv::AbstractVector{Float64}) = [quantile(d,u) for u in uv];


text:
#!usr/bin/env julia

include("call_cc.jl")

@show @call_cc cc begin
    5*6+cc(7)
end # ⇒ 9

block(content::Function, except::Function) =
    (@call_cc success begin
        except(@call_cc fail begin
            success(content(fail))
        end)
    end)

block(
    raise -> begin
        function safe_div(a, b)
            if b == 0
                raise("divide by zero")
            else
                div(a, b)
            end
        end
        @show safe_div(3/4)
        @show safe_div(9/6)
        @show safe_div(2/0)
    end,
    errmsg -> begin
        print_with_color(:red, errmsg)
    end
)

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
using GPUBenchmarks, Plots, Colors, BenchmarkTools, UnicodeFun

#########################
# Parameters

# lol we can't just use 0.0, because Plots.jl errors then
missing_val = 0.0 + eps(Float64)
text2pix = 1.0
ywidth = 2.5
window_size = (800, 500)

nice_colors = map([
    (0, 116, 217), # Blue
    (127, 219, 255), # Aqua
    (133,  20,  75), # Maroon
    (240,  18, 190), # Fuchsia
    (46, 204,  64), # Green
    (177,  13, 201), # Purple
    (57, 204, 204), # Teal
    (61, 153, 112), # Olive
    (255, 220,   0), # Yellow
    (255,  65,  54), # Red
    (1, 255, 112), # Lime
    (255, 133,  27), # Orange
    (0,  31,  63), # Navy
    (221, 221, 221), # Silver
    (17,  17,  17), # Black
    (170, 170, 170), # Gray
]) do x
    RGB(map(v-> v / 255, x)...)
end
nice_colors_hash = ["#a6cee3","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6","#6a3d9a","#1f78b4"]
nice_colors = parse.(Colorant, nice_colors_hash)

#######################
# helpers

function get_time(x::BenchmarkTools.Trial)
    time = minimum(x).time
    isinf(time) ? missing_val : time / 10^9
end

function judged_push!(benchset, benchmark, name)
    get_time(benchmark) == missing_val && return
    for benches in benchset
        judgment = judge(minimum(first(benches)[2]), minimum(benchmark))
        if judgment.time == :invariant
            push!(benches, (name, benchmark))
            return
        end
    end
    # we have a new unique benchresult and can open a new group
    push!(benchset, [(name, benchmark)])
    return
end
function get_log_n(N)
    # we should only use log10  or log2 for now!
    isinteger(log10(N)) && return string("10", to_superscript(Int(log10(N))))
    ispow2(N) && return string("2", to_superscript(Int(log2(N))))
    string(N)
end

function prettytime(t)
    if t < 1e3
        1, "ns"
    elseif t < 1e6
        1e3, "μs"
    elseif t < 1e9
        1e6, "ms"
    else
        1e9, "s"
    end
end
function smart_round(t)
    if t < 0.0001
        round(t, 9)
    elseif t < 0.1
        round(t, 4)
    else
        round(t, 1)
    end
end

get_trial(x::BenchmarkTools.Trial) = x
function get_trial(x)
    BenchmarkTools.Trial(
        BenchmarkTools.Parameters(),
        fill(Inf, 1000),
        fill(Inf, 1000),
        typemax(Int),
        typemax(Int)
    )
end

grouptime(x) = minimum(last(first(x))).time
function speedup!(benchset)
    sort!(benchset, by = grouptime)
    slowest = grouptime(benchset[end])
    map(x-> slowest / grouptime(x), benchset)
end

rect(w, h, x, y) = Shape(x + [0,w,w,0], y + [0,0,h,h])
function plot_speedup!(p, position, label, color)
    annotation = text(label, 9, RGB(0.2, 0.2, 0.2), :right, "helvetica")
    ps = annotation.font.pointsize
    w = (10 * ps) / text2pix
    shape = rect(-w, ps * ywidth, (position)...)
    plot!(p, shape, linewidth = 0, linecolor = RGBA(0,0,0,0), color = color, m = (color, stroke(0)))
    annotate!(p, [((position .+ (-4, 10text2pix))..., annotation)])
    position .- (w + (5 * text2pix), 0)
end

function plot_label!(p, position, label, color)
    shape = rect(10, 9 * ywidth, (position .- (12, 0))...)
    plot!(p, shape, linewidth = 0, color = color, linecolor = RGBA(1,1,1,0.2))
    position .- ((15text2pix), 0)
end

function plot_benchset(p, position, wstart, benchset, label_colors, speed_cmap)
    speedups = speedup!(benchset)
    abs_times = map(x-> minimum(x[1][2]).time, benchset)
    iterator = zip(reverse(speedups), speed_cmap, reverse(benchset), reverse(abs_times))
    for (speedup, scolor, benches, abs_time) in iterator
        position = plot_speedup!(p, position, prettytime(abs_time), scolor)
        position = plot_speedup!(p, position, @sprintf("%8.1fx", speedup), scolor)
        for (name, bench) in benches
            position = plot_label!(p, position, name, label_colors[name])
        end
        position = (wstart, position[2] + 11 * ywidth + 5)
    end
    position
end

function plot_samples(suite, baseline, devices, legend_colors)
    Nmax = maximum(unique(map(x-> x.N, suite)))
    sort_dev(dev) = minimum(first(filter(x-> x.device == dev && x.N == Nmax, suite)).benchmark).time
    devices = sort(devices, by = sort_dev)
    max_time = maximum(map(x-> minimum(x.benchmark).time, suite))
    divisor, unit = prettytime(max_time)
    str = IOBuffer()
    Ns = map(x-> x.N, baseline)
    Ns = [Ns[1], Ns[end]]
    print(str, "| device |")
    for n in Ns
        print(str, " N = ", get_log_n(n), " |")
    end
    print(str, "\n| --- | ")
    for i = 1:(length(Ns))
        print(str, " :---: |")
    end
    for device in devices
        c = legend_colors[device]
        colorhex = bytes2hex(reinterpret(UInt8, [c.r, c.g, c.b]))
        label = device_label(device)
        print(str, "\n| ", "![$label](https://placehold.it/15/$(colorhex)/000000?text=+) ", label, " | ")
        for n in Ns
            bench = filter(x-> x.N == n && x.device == device, suite)[1].benchmark
            basetime = minimum(filter(x-> x.N == n, baseline)[1].benchmark).time
            t = minimum(bench).time
            speedup = basetime / t
            print(str, " `", smart_round(t / divisor), " ", unit, "` `", smart_round(speedup), "x` |")
        end
    end
    String(take!(str))
end
function github_url(isimage, name...)
    str = joinpath(
        "https://github.com/JuliaGPU/GPUBenchmarks.jl/blob/master/",
        name...,
        isimage ? "?raw=true" : ""
    )
    # there is a better way to do this in HTTParser or so... keep forgetting where and how
    replace(str, " ", "%20")
end

##########################################
# plotting code

gr(size = window_size)

md_io = open(GPUBenchmarks.dir("results", string("results.md")), "w")

println(md_io, """
# GPU Benchmarks

Packages benchmarked:

[CuArrays](https://github.com/FluxML/CuArrays.jl) appears as: **cuarrays**

[GPUArrays](https://github.com/JuliaGPU/CLArrays.jl) appears as: **clarrays**

[ArrayFire](https://github.com/gaika/ArrayFire.jl) appears as: **arrayfire cl**, **arrayfire cu**

Julia Base Arrays appear as: **julia base**

Hardware used for GPU: **GTX 950**

Hardware used for Julia single and multithreaded backends: **Intel® Core™ i7-6700 CPU @ 3.40GHz × 4**

Julia's Array implementation is used as a baseline for performance and precision.
So the baseline is in most cases the maximum single threaded performance with SIMD acceleration.
The mean difference in the precision compared to the Julia baseline is plotted as points, with the size of difference corelating with point size.

---

""")

function device_label(device)
    # TODO rename devices in GPUArrays and GPUBenchmarks
    str = replace(string(device), "_", " ")
    if device == "opencl"
        "gpuarrays cl"
    elseif device == "cudanative"
        "gpuarrays cudanative"
    elseif device == "julia"
        "gpuarrays threaded"
    else
        str
    end
end
db = GPUBenchmarks.get_database()

most_current = filter(x-> x.timestamp == GPUBenchmarks.last_time_stamp(), GPUBenchmarks.get_database())

names = unique(GPUBenchmarks.device.(most_current))
prevts = sort(unique(GPUBenchmarks.timestamp.(db)))[1]

prevbs = filter(x-> x.timestamp == prevts, db)
namesprev = unique(GPUBenchmarks.device.(prevbs))
to_merge_names = namesprev[[1, 4, 6, 7]]
benchprev = filter(x-> x.device in to_merge_names, prevbs)
cu_pde = filter(x-> x.device == "cudanative" && x.name == "PDE", prevbs)

most_current = vcat(cu_pde, benchprev, most_current)
using GPUBenchmarks: codepath, name

most_current = map(most_current) do bench
    if bench.device == "cudanative"
        return BenchResult(bench, device = "cuarrays")
    else
        bench
    end
end

codepaths = unique(codepath.(most_current))

for elem in most_current
    t = get_time(elem.benchmark)
    if length(t) == 14
        println(elem)
    end
end

xx = filter(x-> x.device == "clarrays_gpu" && x.name == "Juliaset", most_current)
new_most_current = filter(x-> !(x.device == "clarrays_gpu" && x.name == "Juliaset"), most_current)

newxx = BenchResult[]
for N in ns
    push!(newxx, first(filter(x-> x.N == N, xx)))
end
most_current  = vcat(newxx, new_most_current)

for code_path in codepaths
    println(code_path)
    suites = unique(name.(filter(x-> codepath(x) == code_path, most_current)))
    mod = include(code_path)
    jl_name = basename(code_path)
    file_name, ext = splitext(jl_name)
    println(md_io, "### ", titlecase(file_name))
    println(md_io, mod.description)
    for suitename in suites
        suite = filter(x-> name(x) == suitename, most_current)
        println(md_io, "#### ", titlecase(suitename))
        i = 1
        legend_colors = Dict()
        main_plot = plot(
            xaxis = ("Problem size N", :log10), yaxis = "Speedup",
            legend = :topleft,
            background_color_legend = RGBA(1, 1, 1, 0.6),
            top_margin = 0,
            foreground_color_grid = RGB(0.6, 0.6, 0.6),
            axiscolor = RGB(0.2, 0.2, 0.2),
            markerstrokewidth = 0,
        );
        devices = unique(map(x-> x.device, suite))
        baseline = sort(filter(x-> x.device == "julia_base", suite), by = (x)-> x.N)
        base_times, Ns = map(x-> x.benchmark, baseline), map(x-> x.N, baseline)
        base_times = get_time.(base_times)
        for device in devices
            println(device, " ", suitename)
            device_benches = sort(filter(x-> x.device == device, suite), by = (x)-> x.N)
            if length(device_benches) == 14
                device_benches = device_benches[1:2:length(device_benches)]
            end
            times, Ns = map(x-> x.benchmark, device_benches), map(x-> x.N, device_benches)
            meandiff = map(x-> x.meandiffrence, device_benches) .* 3000.0

            println(">", length(times))
            times = base_times ./ get_time.(times)
            color = nice_colors[i]
            legend_colors[device] = color
            error_cmap = linspace(colorant"#E53A15", colorant"#AAE500", length(Ns))
            plot!(main_plot, Ns, times, line = (1, 0.4, color), m = (color, 5, stroke(0)), label = device_label(device))
            i += 1
        end
        legend_str = plot_samples(suite, baseline, devices, legend_colors)

        layout = @layout [
            a{0.5h}
            a{0.5w} a{0.5w}
        ]
        plot(main_plot, size = window_size ./ (1, 2))
        plotbase = GPUBenchmarks.dir("results", "plots")
        isdir(plotbase) || mkdir(plotbase)
        pngpath = joinpath(plotbase, suitename * ".png")
        println(pngpath)
        savefig(pngpath)
        println(pngpath)
        img_url = github_url(true, split(pngpath, Base.Filesystem.path_separator)[end-2:end]...)

        code_url = github_url(false, "benchmark", jl_name)
        println(md_io, "[![$suitename]($img_url)]($code_url)")
        println(md_io)
        println(md_io, legend_str)
        println(md_io)
        println(md_io, "[code]($code_url)")
        println(md_io)
        println(md_io, "___")
        println(md_io)
    end
end
close(md_io)


text:
using OrdinaryDiffEq, Zygote
using DiffEqSensitivity, Test, ForwardDiff

abstol = 1e-12
reltol = 1e-12
savingtimes = 0.5

function test_continuous_callback(cb, g, dg!; only_backsolve=false)
  function fiip(du, u, p, t)
    du[1] = u[2]
    du[2] = -p[1]
  end
  function foop(u, p, t)
    dx = u[2]
    dy = -p[1]
    [dx, dy]
  end

  u0 = [5.0, 0.0]
  tspan = (0.0, 2.5)
  p = [9.8, 0.8]

  prob = ODEProblem(fiip, u0, tspan, p)
  proboop = ODEProblem(fiip, u0, tspan, p)

  sol1 = solve(prob, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes)
  sol2 = solve(prob, Tsit5(), u0=u0, p=p, abstol=abstol, reltol=reltol, saveat=savingtimes)

  if cb.save_positions == [1, 1]
    @test length(sol1.t) != length(sol2.t)
  else
    @test length(sol1.t) == length(sol2.t)
  end

  du01, dp1 = @time Zygote.gradient(
    (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=BacksolveAdjoint())),
    u0, p)

  du01b, dp1b = Zygote.gradient(
    (u0, p) -> g(solve(proboop, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=BacksolveAdjoint())),
    u0, p)

  du01c, dp1c = Zygote.gradient(
    (u0, p) -> g(solve(proboop, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=BacksolveAdjoint(checkpointing=false))),
    u0, p)

  if !only_backsolve
    @test_broken du02, dp2 = @time Zygote.gradient(
      (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=ReverseDiffAdjoint())), u0, p)

    du03, dp3 = @time Zygote.gradient(
      (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=InterpolatingAdjoint(checkpointing=true))),
      u0, p)

    du03c, dp3c = Zygote.gradient(
      (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=InterpolatingAdjoint(checkpointing=false))),
      u0, p)

    du04, dp4 = @time Zygote.gradient(
      (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=QuadratureAdjoint())),
      u0, p)
  end
  dstuff = @time ForwardDiff.gradient(
    (θ) -> g(solve(prob, Tsit5(), u0=θ[1:2], p=θ[3:4], callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes)),
    [u0; p])

  @info dstuff

  @test du01 ≈ dstuff[1:2]
  @test dp1 ≈ dstuff[3:4]
  @test du01b ≈ dstuff[1:2]
  @test dp1b ≈ dstuff[3:4]
  @test du01c ≈ dstuff[1:2]
  @test dp1c ≈ dstuff[3:4]
  if !only_backsolve
    @test_broken du01 ≈ du02
    @test du01 ≈ du03 rtol = 1e-7
    @test du01 ≈ du03c rtol = 1e-7
    @test du03 ≈ du03c
    @test du01 ≈ du04
    @test_broken dp1 ≈ dp2
    @test dp1 ≈ dp3
    @test dp1 ≈ dp3c
    @test dp3 ≈ dp3c
    @test dp1 ≈ dp4 rtol = 1e-7

    @test_broken du02 ≈ dstuff[1:2]
    @test_broken dp2 ≈ dstuff[3:4]
  end

  cb2 = DiffEqSensitivity.track_callbacks(CallbackSet(cb), prob.tspan[1], prob.u0, prob.p, BacksolveAdjoint())
  sol_track = solve(prob, Tsit5(), u0=u0, p=p, callback=cb2, abstol=abstol, reltol=reltol, saveat=savingtimes)

  adj_prob = ODEAdjointProblem(sol_track, BacksolveAdjoint(), dg!, sol_track.t, nothing,
    callback=cb2,
    abstol=abstol, reltol=reltol)
  adj_sol = solve(adj_prob, Tsit5(), abstol=abstol, reltol=reltol)
  @test du01 ≈ -adj_sol[1:2, end]
  @test dp1 ≈ adj_sol[3:4, end]

  # adj_prob = ODEAdjointProblem(sol_track,InterpolatingAdjoint(),dg!,sol_track.t,nothing,
  # 						 callback = cb2,
  # 						 abstol=abstol,reltol=reltol)
  # adj_sol = solve(adj_prob, Tsit5(), abstol=abstol,reltol=reltol)
  #
  # @test du01 ≈ -adj_sol[1:2,end]
  # @test dp1 ≈ adj_sol[3:6,end]
end

println("Continuous Callbacks")
@testset "Continuous callbacks" begin
  @testset "simple loss function bouncing ball" begin
    g(sol) = sum(sol)
    function dg!(out, u, p, t, i)
      (out .= -1)
    end

    @testset "callbacks with no effect" begin
      condition(u, t, integrator) = u[1] # Event when event_f(u,t) == 0
      affect!(integrator) = (integrator.u[2] += 0)
      cb = ContinuousCallback(condition, affect!, save_positions=(false, false))
      test_continuous_callback(cb, g, dg!)
    end
    @testset "callbacks with no effect except saving the state" begin
      condition(u, t, integrator) = u[1]
      affect!(integrator) = (integrator.u[2] += 0)
      cb = ContinuousCallback(condition, affect!, save_positions=(true, true))
      test_continuous_callback(cb, g, dg!)
    end
    @testset "+= callback" begin
      condition(u, t, integrator) = u[1]
      affect!(integrator) = (integrator.u[2] += 50.0)
      cb = ContinuousCallback(condition, affect!, save_positions=(true, true))
      test_continuous_callback(cb, g, dg!)
    end
    @testset "= callback with parameter dependence and save" begin
      condition(u, t, integrator) = u[1]
      affect!(integrator) = (integrator.u[2] = -integrator.p[2] * integrator.u[2])
      cb = ContinuousCallback(condition, affect!, save_positions=(true, true))
      test_continuous_callback(cb, g, dg!)
    end
    @testset "= callback with parameter dependence but without save" begin
      condition(u, t, integrator) = u[1]
      affect!(integrator) = (integrator.u[2] = -integrator.p[2] * integrator.u[2])
      cb = ContinuousCallback(condition, affect!, save_positions=(false, false))
      test_continuous_callback(cb, g, dg!; only_backsolve=true)
    end
    @testset "= callback with non-linear affect" begin
      condition(u, t, integrator) = u[1]
      affect!(integrator) = (integrator.u[2] = integrator.u[2]^2)
      cb = ContinuousCallback(condition, affect!, save_positions=(true, true))
      test_continuous_callback(cb, g, dg!)
    end
    @testset "= callback with terminate" begin
      condition(u, t, integrator) = u[1]
      affect!(integrator) = (integrator.u[2] = -integrator.p[2] * integrator.u[2]; terminate!(integrator))
      cb = ContinuousCallback(condition, affect!, save_positions=(true, true))
      test_continuous_callback(cb, g, dg!; only_backsolve=true)
    end
  end
  @testset "MSE loss function bouncing-ball like" begin
    g(u) = sum((1.0 .- u) .^ 2) ./ 2
    dg!(out, u, p, t, i) = (out .= 1.0 .- u)
    condition(u, t, integrator) = u[1]
    @testset "callback with non-linear affect" begin
      function affect!(integrator)
        integrator.u[1] += 3.0
        integrator.u[2] = integrator.u[2]^2
      end
      cb = ContinuousCallback(condition, affect!, save_positions=(true, true))
      test_continuous_callback(cb, g, dg!)
    end
    @testset "callback with non-linear affect and terminate" begin
      function affect!(integrator)
        integrator.u[1] += 3.0
        integrator.u[2] = integrator.u[2]^2
        terminate!(integrator)
      end
      cb = ContinuousCallback(condition, affect!, save_positions=(true, true))
      test_continuous_callback(cb, g, dg!; only_backsolve=true)
    end
  end
  @testset "MSE loss function free particle" begin
    g(u) = sum((1.0 .- u) .^ 2) ./ 2
    dg!(out, u, p, t, i) = (out .= 1.0 .- u)
    function fiip(du, u, p, t)
      du[1] = u[2]
      du[2] = 0
    end
    function foop(u, p, t)
      dx = u[2]
      dy = 0
      [dx, dy]
    end

    u0 = [5.0, -1.0]
    p = [0.0, 0.0]
    tspan = (0.0, 2.0)

    prob = ODEProblem(fiip, u0, tspan, p)
    proboop = ODEProblem(fiip, u0, tspan, p)

    condition(u, t, integrator) = u[1] # Event when event_f(u,t) == 0
    affect!(integrator) = (integrator.u[2] = -integrator.u[2])
    cb = ContinuousCallback(condition, affect!)

    du01, dp1 = Zygote.gradient(
      (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=BacksolveAdjoint())),
      u0, p)

    dstuff = @time ForwardDiff.gradient(
      (θ) -> g(solve(prob, Tsit5(), u0=θ[1:2], p=θ[3:4], callback=cb, abstol=abstol, reltol=reltol, saveat=savingtimes)),
      [u0; p])

    @info dstuff

    @test du01 ≈ dstuff[1:2]
    @test dp1 ≈ dstuff[3:4]
  end
end


text:
using Documenter

using Pkg
Pkg.activate("../")
using LinRegOutliers

push!(LOAD_PATH,"../src/")
makedocs(
	sitename="LinRegOutliers",
	pages=[
		"datasets.md",
		"types.md",
		"diagnostics.md",
		"algorithms.md"
	]
)

max_stars_repo_name: angeris/LinRegOutliers

max_issues_repo_name: angeris/LinRegOutliers

max_forks_repo_name: angeris/LinRegOutliers


text:
#params for RvLineList.jl

#module params
#export max_spectra_to_use
#export inst, data_path, target_subdir, fits_target_str
#export espresso_mask_filename, orders_to_use, norm_type, ccf_mid_velocity

max_spectra_to_use = 1000000
if max_spectra_to_use < 200
   @warn "param.in setting max_spectra_to_use to " * max_spectra_to_use
end

#expres_data_path = "/home/awise/data/expres/"
#neid_data_path="/home/awise/data/neid/"
#harpsn_data_path="/home/awise/data/harps-n/"
#harps_data_path="/home/awise/data/harps/"

pipeline_output_path_ebf11 = "/gpfs/group/ebf11/default/ebf11/neid_solar/data/pipeline/v1.1/outputs/ebf11/test1/"
#pipeline_output_path_ebf11 = "/home/awise/data/neid/solar/"

pipeline_output_path_afw5465 = "/gpfs/group/ebf11/default/ebf11/neid_solar/data/pipeline/v1.1/outputs/afw5465/test1/"
#pipeline_output_path_afw5465 = "/home/awise/data/neid/solar/"

neid_data_path = "/gpfs/group/ebf11/default/ebf11/neid_solar/data/pipeline/v1.1/L2/"

output_dir = joinpath(pwd(),"outputs")
VALD_output=true #whether or not to carry VALD line data through to the final mask

inst = NEID
#target_subdir="101501"

local fits_target_str
local espresso_mask_filename
local ccf_mid_velocity
local data_path

if inst==EXPRES
   @assert @isdefined expres_data_path
   data_path = joinpath(expres_data_path,target_subdir)
   fits_target_str = target_subdir
   if fits_target_str == "101501"
      espresso_mask_filename = "G9.espresso.mas"
      ccf_mid_velocity = -5.0e3
   end
   if fits_target_str == "10700"
      espresso_mask_filename = "G8.espresso.mas"
      ccf_mid_velocity = -16640.0
   end
   if fits_target_str == "34411"  #G1?
      espresso_mask_filename = "G2.espresso.mas"
      ccf_mid_velocity = 66500.0
   end
   if fits_target_str == "26965" # K0?
      espresso_mask_filename = "G9.espresso.mas"
      ccf_mid_velocity = -40320.0
   end
end

if inst==NEID #we are using solar data
   fits_target_str="neid"
   @assert @isdefined neid_data_path
   data_path = neid_data_path
   espresso_mask_filename = "G2.espresso.mas"
   ccf_mid_velocity = 0.0
end

if inst==HARPS #the target is Alpha Centauri B
   @assert @isdefined harps_data_path
   data_path = joinpath(harps_data_path,target_subdir)
   espresso_mask_filename = "K2.espresso.mas"
   ccf_mid_velocity = -22000.0
end

if inst==HARPSN #we are using solar data
   @assert @isdefined harpsn_data_path
   data_path = joinpath(harpsn_data_path,target_subdir)
   espresso_mask_filename = "G2.espresso.mas"
   ccf_mid_velocity = 0.0
end

#data set params
#EXPRES_excalibur_only = true #whether or not to only include excalibur wavelengths in expres data - this parameter has been replaced with orders_to_use
#orders_to_use = 43:72 #EXPRES excalibur range
#orders_to_use = sort(sample(43:72,10,replace=false)) #Ten random orders from the EXPRES excalibur range
#order_to_use = 12:83 #EXPRES all useable orders
#orders_to_use = orders_to_use_default(NEID2D())
#orders_to_use = 4:118 # all NEID orders with wavelength values
orders_to_use = 17:115 # NEID orders with minimal NaNs in extracted orders from daily flux averages

norm_type = :continuum #normalization type. Options: :raw, :continuum, :blaze

# VALD mask params

# lineprops = lineprops_101501 #which VALD line list to use in mask creation - not works in param file yet because mask names do not use this param as they should
#maskWavelengths = string("Reiners",target) #keyword for mask wavelengths - should be e.g. Reiners or Reiners101501
allowBlends = 0 #number of blends allowed with each line, e.g. 0 for single lines only, 1 for doublets only, [0,1] for singlets and doublets
overlap_cutoff = 1e-5 #distance between lines required for them to be categorized as a blend, expressed as a fraction of the speed of light
depth_cutoff=0.05
iron1Only = "all" # which species to use - options are "Fe1", "nonFe1", "all"
badLineFilter = "none" #which mask to use to filter out bad line - only lines within ~ 3 km/s of one of these mask lines will be kept
rejectTelluricSlope = 2000.0 #derivative of spectrum required for telluric line rejection - a value of 0 turns off telluric rejection
nbin = 1 #number of mask subdivisions to make
#bin_n = 1 #which subdivision to use - one is first
binParam = :depth #which mask parameter to bin multiple masks by - only important for subdividing masks - can be :depth or :lambda
depthPercentile = true # if binning by depth, use sum of depths for bin weights

# Empirical mask params
quant = "90" #quantile for stability of fit params for empirical masks
#min_frac_converged = "90" #minimum fraction of the line fits that converged in the dataset for the empirical line to be used
line_width_50 = 7392.0 #NEID solar line width for ESPRESSO G2 mask, mask_scale_factor = 2.0, other ccf params default values

# RV params
#lsf_width = 2.0e3 #estimated line width for telluric avoidance
#mask_type = :gaussian #mask shape. Options: :tophat, :gaussian, :supergaussian, :halfcos
#mask_scale_factor = 4.0 #approximate number of pixels in the width scale of the mask
#fwtf = 2.0 #fraction of the FWHM of the CCF to fit
#RV_fit_type = :gaussian #type of fit to use on CCF to measure RV. Options: :quadratic, :gaussian

#global tophap_ccf_mask_scale_factor=1.6

#end

max_stars_repo_name: alexander-wise/RvLineList

max_issues_repo_name: alexander-wise/RvLineList

max_forks_repo_name: alexander-wise/RvLineList


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "htslib"
version = v"1.14"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/samtools/htslib/releases/download/$(version.major).$(version.minor)/htslib-$(version.major).$(version.minor).tar.bz2",
                  "ed221b8f52f4812f810eebe0cc56cd8355a5c9d21c62d142ac05ad0da147935f")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/htslib-*
export CPPFLAGS="-I${includedir}"
export LDFLAGS=-L${libdir}
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install

# Delete static library
rm "${prefix}/lib/libhts.a"

# On Windows, product files are renamed so that BB can find them.
if [[ ${target} == *-mingw32 ]]; then
    # Add .exe extension to executables.
    for exefile in bgzip tabix htsfile; do
        mv "${bindir}/${exefile}" "${bindir}/${exefile}${exeext}"
    done
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# NOTE: Configuring i686-w64-mingw32 fails due to 'unable to find the recv()
# function' error. So we skip it for now.
platforms = supported_platforms(; experimental=true, exclude=Sys.iswindows)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libhts", "hts"], :libhts),
    ExecutableProduct("bgzip", :bgzip),
    ExecutableProduct("tabix", :tabix),
    ExecutableProduct("htsfile", :htsfile),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("Bzip2_jll"; compat="1.0.8"),
    Dependency("XZ_jll"),
    Dependency("LibCURL_jll"),
    Dependency("OpenSSL_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
export legendrePoly

function legendrePoly(::Type{T},order::Int) where {T}
    order<0 && error("order must be >= 0")
    P = [Polynomial{T}(1), Polynomial{T}([0,1])]
    x = Polynomial{T}([0,1])
    order < 2 && return P[order+1]
    for i = 1:order-1
        Ptmp = ( (2*i + 1 ) * x * P[2] - i * P[1]) / (i+1)
        P[1] = P[2]
        P[2] = Ptmp
    end
    return P[2]
end

legendrePoly(order::T) where {T} = legendrePoly(T,convert(Int,order))


text:
function azeotrope_pressure(model::EoSModel, T; v0 =nothing)
    TYPE = promote_type(eltype(T),eltype([0.5,0.5]))
#     lb_v = lb_volume(model,x)
    ts = T_scales(model,[0.5,0.5])
    pmix = p_scale(model,[0.5,0.5])
    if v0 === nothing
        v0 = x0_bubble_pressure(model,T,[0.5,0.5])
    end
    len = length(v0[1:end-1])
    #xcache = zeros(eltype(x0),len)
    Fcache = zeros(eltype(v0[1:end-1]),len)
    f! = (F,z) -> Obj_az_pressure(model, F, T, exp10(z[1]), exp10(z[2]), z[3:end],z[3:end],ts,pmix)
    r  =Solvers.nlsolve(f!,v0[1:end-1],LineSearch(Newton()))
    sol = Solvers.x_sol(r)
    v_l = exp10(sol[1])
    v_v = exp10(sol[2])
    y = FractionVector(sol[3:end])
    P_sat = pressure(model,v_l,T,y)
    return (P_sat, v_l, v_v, y)
end

function Obj_az_pressure(model::EoSModel, F, T, v_l, v_v, x, y,ts,ps)
    x   = FractionVector(x)
    y   = FractionVector(y) #julia magic, check misc.jl
    μ_l = VT_chemical_potential(model,v_l,T,x)
    μ_v = VT_chemical_potential(model,v_v,T,y)
    p_l = pressure(model,v_l,T,x)
    p_v = pressure(model,v_v,T,y)
    for i in 1:length(x)
        F[i] = (μ_l[i]-μ_v[i])/(R̄*ts[i])
    end
    F[end] = (p_l-p_v)/ps
    return F
end


text:
""
function run_hpf_iv(file, model_type::Type, optimizer; kwargs...)
    return _PMs.run_model(file, model_type, optimizer, build_hpf_iv; ref_extensions=[ref_add_xfmr!],  solution_processors=[ _HPM.sol_data_model!], multinetwork=true, kwargs...)
end

""
function build_hpf_iv(pm::AbstractPowerModel)
    bounded = false

    for (n, network) in _PMs.nws(pm)
        _PMs.variable_bus_voltage(pm, nw=n, bounded=bounded)
        _PMs.variable_bus_voltage_magnitude_sqr(pm, nw=n, bounded=bounded)
        variable_transformer_voltage(pm, nw=n, bounded=bounded)

        _PMs.variable_branch_current(pm, nw=n, bounded=bounded)
        _PMs.variable_dcline_current(pm, nw=n)
        variable_transformer_current(pm, nw=n, bounded=bounded)

        variable_load_current(pm, nw=n, bounded=bounded)
        variable_gen_current(pm, nw=n, bounded=bounded)

    end

    for (n, network) in _PMs.nws(pm)
        for i in _PMs.ids(pm, :ref_buses, nw=n)
            constraint_ref_bus(pm, i, nw=n)
        end

        for i in _PMs.ids(pm, :bus, nw=n)
            constraint_current_balance(pm, i, nw=n)
            constraint_vm_auxiliary_variable(pm, i, nw=n)
            # constraint_voltage_harmonics_relative_magnitude(pm, i, nw=n)
        end

        for g in _PMs.ids(pm, :gen, nw=n)
            # _PMs.constraint_gen_active_bounds(pm, g, nw=n)
            # _PMs.constraint_gen_reactive_bounds(pm, g, nw=n)
        end

        for i in _PMs.ids(pm, :load, nw=n)
            constraint_load_power(pm, i, nw=n)
        end

        for b in _PMs.ids(pm, :branch, nw=n)
            _PMs.constraint_current_from(pm, b, nw=n)
            _PMs.constraint_current_to(pm, b, nw=n)

            _PMs.constraint_voltage_drop(pm, b, nw=n)

            # _PMs.constraint_current_limit(pm, b, nw=n)
        end

        for t in _PMs.ids(pm, :xfmr, nw=n)
            constraint_transformer_core_excitation(pm, t, nw=n)
            constraint_transformer_core_voltage_drop(pm, t, nw=n)
            constraint_transformer_core_voltage_balance(pm, t, nw=n)
            constraint_transformer_core_current_balance(pm, t, nw=n)

            constraint_transformer_winding_config(pm, t, nw=n)
            constraint_transformer_winding_current_balance(pm, t, nw=n)
        end

        # for d in _PMs.ids(pm, :dcline, nw=n)
        #     _PMs.constraint_dcline_power_losses(pm, d, nw=n)
        # end
    end

    #constraints across harmonics
    fundamental = 1
    # for i in _PMs.ids(pm, :bus, nw=fundamental)
    #     constraint_voltage_magnitude_rms(pm, i)
    #     constraint_voltage_thd(pm, i, fundamental=fundamental)
    # end

    # for b in _PMs.ids(pm, :branch, nw=fundamental)
    #     constraint_current_limit_rms(pm, b)
    # end

    for g in _PMs.ids(pm, :gen, nw=fundamental)
        constraint_active_filter(pm, g, fundamental=fundamental)
    end

    _PMs.objective_min_fuel_and_flow_cost(pm)
    # objective_distortion_minimization(pm)
end


text:
using FastRefine
using LazySets
import FastRefine: forward_network, forward_affine_map, ishull

nnet = read_nnet("nnet/test3.nnet")

delta = 0.05

solver = DimTree(delta)

in_hyper = Hyperrectangle(fill(1.0, 3), fill(1.0, 3))
out_hyper = Hyperrectangle(fill(0.0, 2), fill(10.0, 2))
problem = Problem(nnet, in_hyper, out_hyper)

file = open("results/group5.txt", "a")
print(file, "Test Result of Group 5:\n\n")

#solver4

time4 = 0

for i = 1:1
    timed_result =@timed solve(solver, problem)
    print(file, "DimTree - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time4 += timed_result.time
end

print(file, "Average time: " * string(time4/1) * " s\n\n")

close(file)


text:
#!/usr/bin/env julia

using Documenter, Logging, NarrativeTest, Test
using FunSQL

ENV["LINES"] = "24"
ENV["COLUMNS"] = "80"

if isempty(ARGS)

    @testset "FunSQL" begin

    @info "Running doctests..."
    DocMeta.setdocmeta!(
        FunSQL,
        :DocTestSetup,
        quote
            using FunSQL:
                SQLTable,
                Agg, Append, As, Asc, Bind, Define, Desc, From, Fun, Get,
                Group, Highlight, Iterate, Join, LeftJoin, Limit, Order,
                Partition, Select, Sort, Var, Where, With,
                AGG, AS, ASC, CASE, CTE, DESC, FROM, FUN, GROUP, HAVING, ID, JOIN,
                KW, LIMIT, LIT, OP, ORDER, PARTITION, SELECT, SORT, UNION, VAR,
                WHERE, WINDOW, WITH,
                render
            using Dates
        end)
    with_logger(Logging.ConsoleLogger(stderr, Logging.Warn)) do
        doctest(FunSQL)
    end

    @info "Running narrative tests..."
    NarrativeTest.testset(joinpath(@__DIR__, "../docs/src"))

    end

else
    NarrativeTest.testset(ARGS)
end


text:
# ---
# title: 337. House Robber III
# id: problem337
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Tree, Depth-first Search
# link: <https://leetcode.com/problems/house-robber-iii/description/>
# hidden: true
# ---
#
# The thief has found himself a new place for his thievery again. There is only
# one entrance to this area, called the "root." Besides the root, each house has
# one and only one parent house. After a tour, the smart thief realized that
# "all houses in this place forms a binary tree". It will automatically contact
# the police if two directly-linked houses were broken into on the same night.
#
# Determine the maximum amount of money the thief can rob tonight without
# alerting the police.
#
# **Example 1:**
#
#
#
#     Input: [3,2,3,null,3,null,1]
#
#          3
#         / \
#        2   3
#         \   \
#          3   1
#
#     Output: 7
#     Explanation:  Maximum amount of money the thief can rob = 3 + 3 + 1 = **7**.
#
# **Example 2:**
#
#
#
#     Input: [3,4,5,1,3,null,1]
#
#          3
#         / \
#        4   5
#       / \   \
#      1   3   1
#
#     Output: 9
#     Explanation:  Maximum amount of money the thief can rob = 4 + 5 = **9**.
#
#
#
## @lc code=start
using LeetCode

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


text:
mutable struct DeployKey <: GitHubType
    id::Union{Int, Nothing}
    key::Union{String, Nothing}
    url::Union{HTTP.URI, Nothing}
    title::Union{String, Nothing}
    verified::Union{Bool, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    read_only::Union{Bool, Nothing}
end

DeployKey(data::Dict) = json2github(DeployKey, data)

namefield(key::DeployKey) = key.id

@api_default function deploykey(api::GitHubAPI, repo, deploykey_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/keys/$(name(deploykey_obj))"; options...)
    return DeployKey(result)
end

@api_default function deploykeys(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/keys"; options...)
    return map(DeployKey, results), page_data
end

@api_default function create_deploykey(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/keys"; options...)
    return DeployKey(result)
end

@api_default function delete_deploykey(api::GitHubAPI, repo, item; options...)
    return gh_delete(api, "/repos/$(name(repo))/keys/$(name(item))"; options...)
end


text:
function write_job_header(f, job::Job, environment)
    environment !== nothing && write(f, environment)
    job_header = job.header
    for line in job_header
        if occursin("\n", line)
            write(f, line)
        else
            write(f, line * "\n")
        end
    end
    modules = AbstractString[]
    for e in map(x->x.exec, job.calculations)
        for m in e.modules
            if !(m ∈ modules)
                push!(modules, m)
            end
        end
    end
    if !isempty(modules)
        write(f, """module load $(join(modules, " "))\n""")
    end
end

# function writetojob(f, job, calculations::Vector{Calculation{Abinit}}, environment; kwargs...)
#     abinit_jobfiles = write_abi_datasets(calculations, job.dir; kwargs...)
#     abifiles = String[]
#     num_abi = 0
#     for (filename, pseudos, runcommand) in abinit_jobfiles
#         push!(abifiles, filename)
#         file, ext = splitext(filename)
#         write(f,
#               "$runcommand << !EOF\n$filename\n$(file * ".out")\n$(job.name * "_Xi$num_abi")\n$(job.name * "_Xo$num_abi")\n$(job.name * "_Xx$num_abi")\n")
#         for pp in pseudos
#             write(f, "$pp\n")
#         end
#         write(f, "!EOF\n")
#         num_abi += 1
#     end
#     return abifiles
# end

function writetojob(f, job, calculations::Vector{Calculation{Elk}}, environment; kwargs...)
    save(calculations, job.structure, job.dir; kwargs...)
    should_run = any(map(x -> x.run, calculations))
    if !should_run
        write(f, "#")
    end
    writeexec(f, calculations[1].exec, environment)
    write(f, "< $(calculations[1].infile) > $(calculations[1].outfile)\n")
    return calculations
end

function writeexec(f, exec::Exec, environment)
    if exec.parallel
        @assert environment !== nothing "Exec $(exec.exec) flagged as parallel but no valid environment was supplied."
        write(f, environment.MPI_command * " ")
    end
    write(f, string(exec) * " ")
end

function writetojob(f, job, calculation::Calculation, environment; kwargs...)
    filename   = calculation.infile
    should_run = calculation.run
    save(calculation, job.structure, joinpath(job, calculation.infile); kwargs...)
    if !should_run
        write(f, "#")
    end
    writeexec(f, calculation.exec, environment)
    write(f, "< $filename > $(calculation.outfile)\n")
    return (calculation,)
end

function writetojob(f, job, _calculation::Calculation{Wannier90}, environment; kwargs...)
    filename   = _calculation.infile
    should_run = _calculation.run
    id         = findfirst(isequal(_calculation), job.calculations)
    seedname   = _calculation.name

    nscf_calc = getfirst(x -> Calculations.isnscf(x), job.calculations)
    if nscf_calc !== nothing
        runexec = nscf_calc.exec
        # For elk the setup necessary for the wan_calc needs to be done before writing the wan calculation
        # because it's inside elk.in
        if eltype(nscf_calc) == QE
            pw2wancalculation = qe_generate_pw2wancalculation(_calculation, nscf_calc,
                                                              runexec)
            preprocess = pop!(_calculation, :preprocess, false)
            wannier_plot = pop!(_calculation, :wannier_plot, nothing)

            if !preprocess || !should_run
                write(f, "#")
            end
            writeexec(f, _calculation.exec, environment)
            write(f,
                  "-pp $filename > $(_calculation.outfile)\n")

            save(_calculation, job.structure, joinpath(job, _calculation.infile); kwargs...)
            writetojob(f, job, pw2wancalculation, environment; kwargs...)
            _calculation[:preprocess] = preprocess
            wannier_plot !== nothing && (_calculation[:wannier_plot] = wannier_plot)
        elseif eltype(nscf_calc) == Elk
            pw2wancalculation = job["elk2wannier"]
        end
    end

    if !should_run
        write(f, "#")
    end
    writeexec(f, _calculation.exec, environment)
    write(f, "$filename > $(_calculation.outfile)\n")
    return (_calculation,)
end

#TODO: This should take scratch dir of server

function write_job_preamble(f, job::Job) end
function write_job_postamble(f, job::Job) end
# function write_job_preamble(f, job::Job)
#     if !isempty(job.server_dir)
#         dir = splitpath(job.dir)[end]
#         write(f, "cp -r $(job.dir) $(job.server_dir) \n")
#         write(f, "cd $(joinpath(job.server_dir, dir)) \n")
#     end
# end

# function write_job_postamble(f, job::Job)
#     if !isempty(job.server_dir)
#         dir = splitpath(job.dir)[end]
#         write(f, "cp -r $(joinpath(job.server_dir, dir)) $(splitdir(job.dir)[1])\n")
#         write(f, "rm -r $(joinpath(job.server_dir, dir))\n")
#     end
# end

"""
    write_job_files(job::Job; kwargs...)

Writes all the calculation files and job file that are linked to a Job.
Kwargs will be passed down to various writetojob functions.
"""
function write_job_files(job::Job, environment; kwargs...)

    if any(x ->eltype(x) == QE, job.calculations)
        for a in unique(job.structure.atoms)
            write(joinpath(job, "$(a.element.symbol).UPF"), a.pseudo)
        end
    end

    open(joinpath(job, "job.tt"), "w") do f
        write(f, "#!/bin/bash\n")
        write(f, "#SBATCH -J $(job.name) \n")
        write_job_header(f, job, environment)
        write_job_preamble(f, job)
        written_calculations = Calculation[]
        abicalculations = filter(x -> eltype(x) == Abinit, job.calculations)
        !isempty(abicalculations) && writetojob(f, job, abicalculations, environment; kwargs...)
        elkcalculations = filter(x -> eltype(x) == Elk, job.calculations)
        !isempty(elkcalculations) &&
            append!(written_calculations, writetojob(f, job, elkcalculations, environment; kwargs...))

        for i in job.calculations
            if i ∉ written_calculations
                append!(written_calculations, writetojob(f, job, i, environment; kwargs...))
            end
        end
        return write_job_postamble(f, job)
    end
end

function read_job_line(line)
    line = strip(line)
    line = replace(line, ['>', '<'] => " ")
    if line[1] == '#'
        run = false
        line = strip(line[2:end], '#')
    else
        run = true
    end
    spl = strip_split(line)

    calculation = spl[end-1]
    output = spl[end]
    spl = spl[1:end-2]
    exec_and_flags = Pair{String,Vector{SubString}}[]
    #TODO This is not really nice, we don't handle execs that are unparseable...
    #     Not sure how we can actually do this
    for s in spl
        d, e = splitdir(s)
        if e ∈ Calculations.allexecs()
            push!(exec_and_flags, s => SubString[])
        elseif !isempty(exec_and_flags)
            # else
            push!(last(exec_and_flags[end]), s)
        end
    end

    execs = Exec[]
    for (e, flags) in exec_and_flags
        dir, efile = splitdir(e)
        dir = replace(dir, "~" => homedir())
        if occursin("pw2wannier90", efile)
            continue
        end
        if occursin("mpi", e)
            push!(execs, Exec(exec=efile, dir=dir, flags=Calculations.parse_mpi_flags(flags)))
        elseif efile == "wannier90.x"
            push!(execs, Exec(exec=efile, dir=dir, flags=Calculations.parse_wan_execflags(flags)))
        elseif any(occursin.(Calculations.QE_EXECS, (efile,)))
            push!(execs, Exec(exec=efile, dir=dir, flags=Calculations.parse_qe_execflags(flags)))
        elseif any(occursin.(Calculations.ELK_EXECS, (efile,)))
            calculation = "elk.in"
            output = "elk.out"
            push!(execs, Exec(; exec = efile, dir = dir))
        else
            push!(execs, Exec(exec=efile, dir=dir, flags=Calculations.parse_generic_flags(flags)))
        end
    end
    return execs, calculation, output, run
end

# TODO: make this work again
# function read_job_filenames(job_file::String)
#     calculation_files = String[]
#     output_files = String[]
#     open(job_file, "r") do f
#         readline(f)
#         while !eof(f)
#             line = readline(f)
#             if isempty(line)
#                 continue
#             end
#             if occursin(".x", line)
#                 runcommand, exec, calculation, output, run = read_job_line(line)
#                 !in(calculation,  calculation_files)  && push!(calculation_files,  calculation)
#                 !in(output, output_files) && push!(output_files, output)
#             end
#         end
#     end
#     return calculation_files, output_files
# end

function calculationparser(exec::Exec)
    if Calculations.is_qe_exec(exec)
        qe_read_calculation
    elseif Calculations.is_wannier_exec(exec)
        wan_read_calculation
    elseif Calculations.is_elk_exec(exec)
        elk_read_calculation
    end
end

function read_job_script(job_file::String)
    name = ""
    dir = splitdir(job_file)[1]
    header = Vector{String}()
    scratch_dir = ""
    calcs = NamedTuple{(:exec, :infile, :outfile, :run), Tuple{Exec, String, String, Bool}}[]
    open(job_file, "r") do f
        readline(f)
        while !eof(f)
            line = readline(f)

            isempty(line) && continue

            if occursin("#SBATCH", line)
                if occursin("-J", line) || occursin("--job-name", line)
                    name = split(split(line, "=")[end])[end]
                else
                    push!(header, line)
                end
            elseif occursin("cp -r", line) && isempty(scratch_dir)
                scratch_dir = split(line)[end]
            elseif Calculations.has_parseable_exec(line)
                execs, infile, outfile, run = read_job_line(line)
                inpath = joinpath(dir, infile)
                if ispath(inpath)
                    exec = getfirst(Calculations.isparseable, execs)
                    if exec !== nothing
                        if length(execs) == 1
                            exec.parallel = false
                        end
                        push!(calcs, (exec = exec, infile = joinpath(dir, infile), outfile = joinpath(dir, outfile), run = run))
                    end
                end
            else
                push!(header, line)
            end
        end
    end

    module_line_ids = findall(x -> occursin("module", x), header)

    if module_line_ids !== nothing
        module_lines = header[module_line_ids]
        deleteat!(header, module_line_ids)
        modules = map(x -> split(x)[end], module_lines)
    end
    #TODO cleanup
    known_es = Calculations.load_execs()
    execs = filter(x -> !any(y -> y.dir == x.dir && y.exec == x.exec, values(known_es)), unique(map(x->x.exec, calcs)))
    for e in execs
        i = 1
        runnable = Calculations.isrunnable(e)
        while length(e.modules) < length(modules) && !runnable
            push!(e.modules, modules[i])
            runnable = Calculations.isrunnable(e)
            i+=1
        end
        if runnable
            e.name = """$(e.exec)@$(join(splitpath(e.dir)[2:end], "_"))"""
            Calculations.maybe_register(e)
        end
    end
    for c in calcs
        e1 = c.exec
        for e in [execs; known_es]
            if e.dir == e1.dir && e.exec == e1.exec
                e1.modules = e.modules
            end
        end
    end

    runtime = Jobs.environment_from_jobscript(job_file)
    ename = Jobs.environment_name(runtime)
    if ename !== nothing
        environment = ename
        #TODO only works with slurm!!!
        deleteat!(header, findall(x -> x[1:7] == "#SBATCH", header))
        deleteat!(header, findall(x -> x[1:6] == "export", header))
    else
        environment = ""
    end
    return (name, calcs, header, environment)
end

function parse_calculations(calcs)
    structures = Structure[]
    outcalcs = Calculation[]
    for calc in calcs
        c = calculationparser(calc.exec)(calc.infile; exec = calc.exec, infile = splitpath(calc.infile)[end], outfile = splitpath(calc.outfile)[end], run = calc.run)
        if c[2] !== nothing
            push!(structures, c[2])
        end
        if c[1] !== nothing
            push!(outcalcs, c[1])
        end
    end
    if !isempty(structures)
        structure = Structures.mergestructures(structures)
    else
        structure = Structure()
        @warn "No valid structures could be read from calculation files."
    end
    return (calculations = outcalcs, structure=structure)
end


text:
using Test
using CUnion

@testset "CUnion tests" begin
    @union struct U
        x::UInt16
        y::UInt32
        struct z
            a::UInt8
            b::UInt8
            c::UInt16
        end
    end
    u = U(0xabcd0123)
    @test u.x === 0x0123
    @test u.y === 0xabcd0123
    @test u.z.a === 0x23
    @test u.z.b === 0x01
    @test u.z.c === 0xabcd
    @test isprimitivetype(U)
    @test reinterpret(UInt32, u) === 0xabcd0123
    @test U(Dict(unionfields(U))[:z](0x23, 0x01, 0xabcd)) === u
end #testset


text:
"""
This function does nothing if the method is not specified for the project's buildphase.
"""
function retire_project!(project::P) where P <: Project{<: BuildPhase}
    return project
end

"""
This function converts Existing, Planned and Queued Projects to Retired buildphase.
"""
function retire_project!(project::P) where P <: Union{Project{Existing}, Project{Planned}, Project{Queue}}
    return convert(Project{Retired}, project)
end

"""
This function does nothing if PSY System is not defined.
"""
function remove_system_component!(sys::Nothing,
                                  project::P,
                                  ) where P <: Project{<: BuildPhase}
    return
end

"""
This function removes Existing projects from the list of Devices in PSY System.
"""
function remove_system_component!(sys::PSY.System,
                                   project::P,
                                  ) where P <: Project{Existing}
    component = PSY.get_components_by_name(PSY.Device, sys, get_name(project))[1]
    PSY.remove_component!(sys, component)
end

"""
This function does nothing if the method is not specified for the project's type and buildphase.
"""
function remove_renewable_gen_data!(project::P,
                                     simulation_dir::String,
                                     ) where P <: Project{<:BuildPhase}
    return
 end

"""
This function removes the PSY System timeseries data for Existing RenewableGen projects.
"""
  function remove_renewable_gen_data!(project::RenewableGenEMIS{Existing},
                                     simulation_dir::String,
                                     )
        load_n_vg_df =  read_data(joinpath(simulation_dir, "timeseries_data_files", "Net Load Data", "load_n_vg_data.csv"))
        DataFrames.select!(load_n_vg_df, DataFrames.Not(get_name(project)))
        write_data(joinpath(simulation_dir, "timeseries_data_files", "Net Load Data"), "load_n_vg_data.csv", load_n_vg_df)
    return
 end

 """
This function removes the future realized profits of projects which are being retired.
"""
function remove_future_profits!(project::Project, iteration_year::Int64)
    products = get_products(project)
    finance_data = get_finance_data(project)
    profit_array_length =  size(get_realized_profit(finance_data), 2)
    for year in iteration_year:profit_array_length
        for product in get_name.(products)
            set_realized_profit!(finance_data, product, year, 0.0)
        end
    end
    return
end

"""
This function does nothing if the method is not specified for the project's buildphase.
"""
function take_retirement_decision(project::P,
                                  index::Int64,
                                  projects::Vector{<: Project{<: BuildPhase}},
                                  scenario_data::Vector{Scenario},
                                  market_prices::MarketPrices,
                                  risk_preference::R,
                                  sys::Union{Nothing, PSY.System},
                                  simulation_dir::String,
                                  iteration_year::Int64,
                                  yearly_horizon::Int64,
                                  simulation_years::Int64,
                                  rep_hour_weight::Vector{Float64},
                                  capacity_forward_years::Int64,
                                  solver::JuMP.MOI.OptimizerWithAttributes) where {P <: Project{<: BuildPhase}, R <: RiskPreference}

    return
end

"""
This function decides whether or not to retire Existing and Planned projects.
"""
function take_retirement_decision(project::P,
                                  index::Int64,
                                  projects::Vector{<: Project{<: BuildPhase}},
                                  scenario_data::Vector{Scenario},
                                  market_prices::MarketPrices,
                                  risk_preference::R,
                                  sys::Union{Nothing, PSY.System},
                                  simulation_dir::String,
                                  iteration_year::Int64,
                                  yearly_horizon::Int64,
                                  simulation_years::Int64,
                                  rep_hour_weight::Vector{Float64},
                                  capacity_forward_years::Int64,
                                  solver::JuMP.MOI.OptimizerWithAttributes) where {P <: Union{Project{Existing}, Project{Planned}}, R <: RiskPreference}

    finance_data = get_finance_data(project)
    queue_cost = get_queue_cost(finance_data)

    if  get_construction_year(project) <= simulation_years

        update_project_utility!(project,
                              scenario_data,
                              market_prices,
                              risk_preference,
                              rep_hour_weight,
                              iteration_year,
                              yearly_horizon,
                              queue_cost,
                              capacity_forward_years,
                              solver)

        project_utility = get_expected_utility(finance_data)[iteration_year]

        annual_revenue = 0.0
        for scenario in scenario_data
            annual_revenue += get_probability(scenario) *
                             sum(get_scenario_profit(finance_data)[get_name(scenario)][iteration_year][:, iteration_year])
        end

        annual_profit = annual_revenue - get_fixed_OM_cost(finance_data)

        if annual_profit < -1 && project_utility < 0
            set_retirement_year!(projects[index], iteration_year)
            projects[index] = retire_project!(project)
            remove_system_component!(sys, project)
            remove_renewable_gen_data!(project, simulation_dir)
            remove_future_profits!(project, iteration_year)
        end
    end

    return

end

"""
This function decides whether or not to retire Queued projects.
"""
function take_retirement_decision(project::P,
                                  index::Int64,
                                  projects::Vector{<: Project{<: BuildPhase}},
                                  scenario_data::Vector{Scenario},
                                  market_prices::MarketPrices,
                                  risk_preference::R,
                                  sys::Union{Nothing, PSY.System},
                                  simulation_dir::String,
                                  iteration_year::Int64,
                                  yearly_horizon::Int64,
                                  simulation_years::Int64,
                                  rep_hour_weight::Vector{Float64},
                                  capacity_forward_years::Int64,
                                  solver::JuMP.MOI.OptimizerWithAttributes) where {P <: Project{Queue}, R <: RiskPreference}

    finance_data = get_finance_data(project)
    queue_cost = get_queue_cost(finance_data)

    if  get_construction_year(project) <= simulation_years

        update_project_utility!(project,
                              scenario_data,
                              market_prices,
                              risk_preference,
                              rep_hour_weight,
                              iteration_year,
                              yearly_horizon,
                              queue_cost,
                              capacity_forward_years,
                              solver)

        project_utility = get_expected_utility(finance_data)[iteration_year]

        if project_utility < 0
            set_retirement_year!(projects[index], iteration_year)
            projects[index] = retire_project!(project)
            remove_future_profits!(project, iteration_year)
        end
    end

    return

end

"""
This function makes the retirement decisions each year.
Returns nothing.
"""
function retire_unprofitable!(investor::Investor,
                              sys::Union{Nothing, PSY.System},
                              simulation_dir::String,
                              iteration_year::Int64,
                              yearly_horizon::Int64,
                              simulation_years::Int64,
                              capacity_forward_years::Int64,
                              solver::JuMP.MOI.OptimizerWithAttributes)

    scenario_data = get_scenario_data(get_forecast(investor))
    market_prices = get_market_prices(investor)
    risk_preference = get_risk_preference(investor)
    projects = get_projects(investor)
        for i in 1:length(projects)
        take_retirement_decision(projects[i],
                                 i,
                                 projects,
                                 scenario_data,
                                 market_prices,
                                 risk_preference,
                                 sys,
                                 simulation_dir,
                                 iteration_year,
                                 yearly_horizon,
                                 simulation_years,
                                 get_rep_hour_weight(investor),
                                 capacity_forward_years,
                                 solver)

        end
    return
end


text:
using Documenter, PhysicalConstants, Unitful

## Generate list of constants
open(joinpath(@__DIR__, "src", "constants.md"), "w") do io
    print(io, """
              ## List of Set of Constants

              !!! note

                  Each dataset listed below exports by default only the full
                  long names of the constants.  Short aliases are provided for
                  convenience, but they are not exported, to avoid polluting
                  the main namespace with dozens of short-named variables.
                  Users can to import the short names of the variables they
                  use most frequently, as shown in the examples above.
              """
          )
    for set in (PhysicalConstants.CODATA2014, PhysicalConstants.CODATA2018)
        println(io)
        println(io, "### ", nameof(set))
        println(io)
        constants = names(set)
        others = setdiff(names(set, all = true), constants)
        println(io, "| Long name | Short | Value | Unit |")
        println(io, "| --------- | ----- | ----- | ---- |")
        for constant in constants
            c = getfield(set, constant)
            if c isa PhysicalConstants.PhysicalConstant
                sym = others[findall(x -> c === getfield(set, x), others)][1]
                println(io, "| `", constant, "` | `", sym, "` | ", ustrip(float(c)), " | ",
                        unit(c) == Unitful.NoUnits ? "" : "`$(unit(c))`", " |")
            end
        end
    end
end

## Build docs
makedocs(
    modules = [PhysicalConstants],
    sitename = "PhysicalConstants",
    pages    = [
        "Introduction" => "index.md",
        "Usage" => "usage.md",
        "List of Constants" => "constants.md",
        "API Reference" => "reference.md",
    ]
)

### Deploy docs
deploydocs(
    repo = "github.com/JuliaPhysics/PhysicalConstants.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
    push_preview = true,
)

max_issues_repo_name: JuliaPhysics/Constants

max_forks_repo_name: JuliaPhysics/Constants


text:
using StaticArrays

@enum GlassType MODEL MIL AGF OTHER AIR

"""
Object identifying a glass, containing a type (e.g. `MODEL`, `MIL`, `OTHER` or `AGF`) depending on how the glass is defined, and an integer ID.
Air is `AIR:0`, others are on the form `AGF:N`, for example.
"""
struct GlassID
    type::GlassType
    num::Int
end

Base.show(io::IO, a::GlassID) = print(io, "$(string(a.type)):$(a.num)")
type(a::GlassID) = a.type
num(a::GlassID) = a.num

ismodel(a::GlassID) = a.type === MODEL
isMIL(a::GlassID) = a.type === MIL
isAGF(a::GlassID) = a.type === AGF
isother(a::GlassID) = a.type === OTHER
isair(a::GlassID) = a.type === AIR

"""
Abstract type encapsulating all glasses.
"""
abstract type AbstractGlass end

"""
Stores all attributes relating to a glass type specified in an .AGF glass catalog.

Never used directly, instead created using catalog glasses, e.g. `GlassCat.SCHOTT.N_BK7`.

In order to prevent type ambiguities in OpticSim.jl we can't have this type paramaterized.
"""
struct Glass <: AbstractGlass
    ID::GlassID
    dispform::Int
    C1::Float64
    C2::Float64
    C3::Float64
    C4::Float64
    C5::Float64
    C6::Float64
    C7::Float64
    C8::Float64
    C9::Float64
    C10::Float64
    λmin::Float64
    λmax::Float64
    D₀::Float64
    D₁::Float64
    D₂::Float64
    E₀::Float64
    E₁::Float64
    λₜₖ::Float64
    temp::Float64
    ΔPgF::Float64
    PR::Float64
    relcost::Float64
    TCE::Float64
    CR::Float64
    status::Int
    SR::Float64
    transmission::Union{Nothing,SVector{100,SVector{3,Float64}}} # could use Vector on 1.5.X, but on 1.6 this causes allocs so we need to use a SVector...
    transmissionN::Int
    Nd::Float64
    AR::Float64
    FR::Float64
    exclude_sub::Int
    Vd::Float64
    ignore_thermal_exp::Int
    p::Float64
    meltfreq::Int

    function Glass(ID::GlassID, dispform, C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, λmin, λmax, D₀, D₁, D₂, E₀, E₁, λₜₖ, temp, ΔPgF, PR, relcost, TCE, CR, status, SR, transmission, Nd, AR, FR, exclude_sub, Vd, ignore_thermal_exp, p, meltfreq)
        # need a constructor to massage the transmission data
        if transmission === nothing
            transmission_s = nothing
            transmissionN = -1
        else
            fill = [SVector(0.0, 0.0, 0.0) for _ in 1:(100 - length(transmission))]
            transmission_s = SVector{100,SVector{3,Float64}}(transmission..., fill...)
            transmissionN = length(transmission)
        end
        return new(ID, dispform, C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, λmin, λmax, D₀, D₁, D₂, E₀, E₁, λₜₖ, temp, ΔPgF, PR, relcost, TCE, CR, status, SR, transmission_s, transmissionN, Nd, AR, FR, exclude_sub, Vd, ignore_thermal_exp, p, meltfreq)
    end
end

"""
    glassid(g::AbstractGlass) -> GlassID

Get the ID of the glass, see [`GlassID`](@ref).
"""
glassid(g::Glass) = g.ID
"""
    glassname(g::Union{AbstractGlass,GlassID})

Get the name (including catalog) of the glass, or glass with this ID.
"""
glassname(g::Glass) = glassname(g.ID)
function glassname(ID::GlassID)
    if ismodel(ID)
        return "GlassCat.ModelGlass.$(num(ID))"
    elseif isMIL(ID)
        return "GlassCat.GlassFromMIL.$(num(ID))"
    elseif isair(ID)
        return "GlassCat.Air"
    elseif isother(ID)
        return OTHER_GLASS_NAMES[num(ID)]
    else
        return AGF_GLASS_NAMES[num(ID)]
    end
end

function Base.show(io::IO, g::Glass)
    print(io, glassname(g))
end

struct AirType <: AbstractGlass end
Base.show(io::IO, ::AirType) = print(io, "Air")
glassid(::AirType) = GlassID(AIR, 0)
glassname(::AirType) = "GlassCat.Air"

"""
    isair(a) -> Bool

Tests if `a` is Air.
"""
isair(::AirType) = true
isair(::AbstractGlass) = false

"""
    glassforid(ID::GlassID)

Get the glass for a given ID.
"""
function glassforid(ID::GlassID)
    if ismodel(ID)
        return MODEL_GLASSES[num(ID)]::Glass
    elseif isMIL(ID)
        return MIL_GLASSES[num(ID)]::Glass
    elseif isair(ID)
        return GlassCat.Air
    elseif isother(ID)
        return OTHER_GLASSES[num(ID)]::Glass
    else
        return AGF_GLASSES[num(ID)]::Glass
    end
end

"""
    info([io::IO], glass::AbstractGlass)

Print out all data associated with `glass` in an easily readable format.

# Examples
```julia-repl
julia> info(GlassCat.RPO.IG4)
ID:                                                AGF:52
Dispersion formula:                                Schott (1)
Dispersion formula coefficients:
     a₀:                                           6.91189161
     a₁:                                           -0.000787956404
     a₂:                                           -4.22296071
     a₃:                                           142.900646
     a₄:                                           -1812.32748
     a₅:                                           7766.33028
Valid wavelengths:                                 3.0μm to 12.0μm
Reference temperature:                              20.0°C
Thermal ΔRI coefficients:
     D₀:                                           3.24e-5
     D₁:                                           0.0
     D₂:                                           0.0
     E₀:                                           0.0
     E₁:                                           0.0
     λₜₖ:                                          0.0
TCE (÷1e-6):                                       20.4
Ignore thermal expansion:                          false
Density (p):                                       4.47g/m³
ΔPgF:                                              0.0
RI at sodium D-Line (587nm):                       1.0
Abbe Number:                                       0.0
Cost relative to N_BK7:                              ?
Status:                                            Standard (0)
Melt frequency:                                    0
Exclude substitution:                              false
```
"""
function info(io::IO, glass::Glass)
    D = glass.dispform

    println(io, "$(rpad("ID:", 50)) $(glass.ID)")

    if (D == -2)
        println(io, "$(rpad("Dispersion formula:", 50)) Cauchy (-2)")
    elseif (D == -1)
        println(io, "$(rpad("Dispersion formula:", 50)) Fitted for model/MIL glass")
    else
        println(io, "$(rpad("Dispersion formula:", 50)) $(DISPFORM_NAMES[D]) ($D)")
    end
    println(io, "Dispersion formula coefficients:")
    if (D == -2) # Cauchy
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C:", 45)) $(glass.C3)")
        println(io, "     $(rpad("D:", 45)) $(glass.C4)")
        println(io, "     $(rpad("E:", 45)) $(glass.C5)")
        println(io, "     $(rpad("F:", 45)) $(glass.C6)")
    elseif (D == -1)
        println(io, "     $(rpad("C₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("C₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("C₃:", 45)) $(glass.C4)")
    elseif (D == 1) # Schott
        println(io, "     $(rpad("a₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("a₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("a₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("a₃:", 45)) $(glass.C4)")
        println(io, "     $(rpad("a₄:", 45)) $(glass.C5)")
        println(io, "     $(rpad("a₅:", 45)) $(glass.C6)")
    elseif (D == 2)  # Sellmeier1
        println(io, "     $(rpad("K₁:", 45)) $(glass.C1)")
        println(io, "     $(rpad("L₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("K₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("L₂:", 45)) $(glass.C4)")
        println(io, "     $(rpad("K₃:", 45)) $(glass.C5)")
        println(io, "     $(rpad("L₃:", 45)) $(glass.C6)")
    elseif (D == 3)  # Herzberger
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C:", 45)) $(glass.C3)")
        println(io, "     $(rpad("D:", 45)) $(glass.C4)")
        println(io, "     $(rpad("E:", 45)) $(glass.C5)")
        println(io, "     $(rpad("F:", 45)) $(glass.C6)")
    elseif (D == 4)  # Sellmeier2
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("λ₁:", 45)) $(glass.C3)")
        println(io, "     $(rpad("B₂:", 45)) $(glass.C4)")
        println(io, "     $(rpad("λ₂:", 45)) $(glass.C5)")
    elseif (D == 5)  # Conrady
        println(io, "     $(rpad("n₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("A:", 45)) $(glass.C2)")
        println(io, "     $(rpad("B:", 45)) $(glass.C3)")
    elseif (D == 6)  # Sellmeier3
        println(io, "     $(rpad("K₁:", 45)) $(glass.C1)")
        println(io, "     $(rpad("L₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("K₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("L₂:", 45)) $(glass.C4)")
        println(io, "     $(rpad("K₃:", 45)) $(glass.C5)")
        println(io, "     $(rpad("L₃:", 45)) $(glass.C6)")
        println(io, "     $(rpad("K₄:", 45)) $(glass.C7)")
        println(io, "     $(rpad("L₄:", 45)) $(glass.C8)")
    elseif (D == 7) || (D == 8)  # HandbookOfOptics1/2
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C:", 45)) $(glass.C3)")
        println(io, "     $(rpad("D:", 45)) $(glass.C4)")
    elseif (D == 9)  # Sellmeier4
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C:", 45)) $(glass.C3)")
        println(io, "     $(rpad("D:", 45)) $(glass.C4)")
        println(io, "     $(rpad("E:", 45)) $(glass.C5)")
    elseif (D == 10) || (D == 12) # Extended1/2
        println(io, "     $(rpad("a₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("a₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("a₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("a₃:", 45)) $(glass.C4)")
        println(io, "     $(rpad("a₄:", 45)) $(glass.C5)")
        println(io, "     $(rpad("a₅:", 45)) $(glass.C6)")
        println(io, "     $(rpad("a₆:", 45)) $(glass.C7)")
        println(io, "     $(rpad("a₇:", 45)) $(glass.C8)")
    elseif (D == 11)  # Sellmeier5
        println(io, "     $(rpad("K₁:", 45)) $(glass.C1)")
        println(io, "     $(rpad("L₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("K₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("L₂:", 45)) $(glass.C4)")
        println(io, "     $(rpad("K₃:", 45)) $(glass.C5)")
        println(io, "     $(rpad("L₃:", 45)) $(glass.C6)")
        println(io, "     $(rpad("K₄:", 45)) $(glass.C7)")
        println(io, "     $(rpad("L₄:", 45)) $(glass.C8)")
        println(io, "     $(rpad("K₅:", 45)) $(glass.C9)")
        println(io, "     $(rpad("L₅:", 45)) $(glass.C10)")
    elseif (D == 13)  # Extended3
        println(io, "     $(rpad("a₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("a₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("a₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("a₃:", 45)) $(glass.C4)")
        println(io, "     $(rpad("a₄:", 45)) $(glass.C5)")
        println(io, "     $(rpad("a₅:", 45)) $(glass.C6)")
        println(io, "     $(rpad("a₆:", 45)) $(glass.C7)")
        println(io, "     $(rpad("a₇:", 45)) $(glass.C8)")
        println(io, "     $(rpad("a₈:", 45)) $(glass.C9)")
    else
        println(io, "     INVALID DISPERSION FORMULA!!")
    end
    println(io, "$(rpad("Valid wavelengths:", 50)) $(glass.λmin)μm to $(glass.λmax)μm")
    println(io, "$(rpad("Reference temperature:", 50)) $(glass.temp)°C")

    if !isnan(glass.D₀) && (glass.D₀ != 0 || glass.D₁ != 0 || glass.D₂ != 0 || glass.E₀ != 0 || glass.E₁ != 0)
        println(io, "Thermal ΔRI coefficients:")
        println(io, "     $(rpad("D₀:", 45)) $(glass.D₀)")
        println(io, "     $(rpad("D₁:", 45)) $(glass.D₁)")
        println(io, "     $(rpad("D₂:", 45)) $(glass.D₂)")
        println(io, "     $(rpad("E₀:", 45)) $(glass.E₀)")
        println(io, "     $(rpad("E₁:", 45)) $(glass.E₁)")
        println(io, "     $(rpad("λₜₖ:", 45)) $(glass.λₜₖ)")
    end

    println(io, "$(rpad("TCE (÷1e-6):", 50)) $(glass.TCE)")
    println(io, "$(rpad("Ignore thermal expansion:", 50)) $(glass.ignore_thermal_exp == 1)")

    println(io, "$(rpad("Density (p):", 50)) $(glass.p)g/m³")
    println(io, "$(rpad("ΔPgF:", 50)) $(glass.ΔPgF)")

    println(io, "$(rpad("RI at sodium D-Line (587nm):", 50)) $(glass.Nd)")
    println(io, "$(rpad("Abbe Number:", 50)) $(glass.Vd)")

    println(io, "$(rpad("Cost relative to N_BK7:", 50)) $(glass.relcost == -1 ? "?" : glass.relcost)")

    if glass.CR != -1 || glass.FR != -1 || glass.SR != -1 || glass.AR != -1 || glass.PR != -1
        println(io, "Environmental resistance:")
        println(io, "     $(rpad("Climate (CR):", 45)) $(glass.CR == -1 ? "?" : glass.CR)")
        println(io, "     $(rpad("Stain (FR):", 45)) $(glass.FR == -1 ? "?" : glass.FR)")
        println(io, "     $(rpad("Acid (SR):", 45)) $(glass.SR == -1 ? "?" : glass.SR)")
        println(io, "     $(rpad("Alkaline (AR):", 45)) $(glass.AR == -1 ? "?" : glass.AR)")
        println(io, "     $(rpad("Phosphate (PR):", 45)) $(glass.PR == -1 ? "?" : glass.PR)")
    end

    println(io, "$(rpad("Status:", 50)) $(STATUS[glass.status + 1]) ($(glass.status))")
    println(io, "$(rpad("Melt frequency:", 50)) $(glass.meltfreq == -1 ? "?" : glass.meltfreq)")
    println(io, "$(rpad("Exclude substitution:", 50)) $(glass.exclude_sub == 1)")

    if glass.transmission !== nothing
        println(io, "Transmission data:")
        println(io, "$(lpad("Wavelength", 15))$(lpad("Transmission", 15))$(lpad("Thickness", 15))")
        for i in 1:(glass.transmissionN)
            λ, t, τ = glass.transmission[i]
            println(io, "$(lpad("$(λ)μm", 15))$(lpad(t, 15))$(lpad("$(τ)mm", 15))")
        end
    end
end

function info(io::IO, ::AirType)
    println(io, "GlassCat.Air")
    println(io, "Material representing air, RI is always 1.0 at system temperature and pressure, absorption is always 0.0.")
end

info(g::AbstractGlass) = info(stdout, g)

function docstring(io::IO, glass::Glass)
    D = glass.dispform

    println(io, "**ID**$(glass.ID)\n")

    if (D == -2)
        println(io, "**Dispersion formula**Cauchy (-2)\n")
    elseif (D == -1)
        println(io, "**Dispersion formula**Fitted for model/MIL glass\n")
    else
        println(io, "**Dispersion formula**$(DISPFORM_NAMES[D]) ($D)\n")
    end
    println(io, "**Dispersion formula coefficients**\n")
    println(io, "|Coefficient|Value|")
    println(io, "|:----------|-----|")
    if (D == -2) # Cauchy
        println(io, "|A |$(glass.C1)|")
        println(io, "|B |$(glass.C2)|")
        println(io, "|C |$(glass.C3)|")
        println(io, "|D |$(glass.C4)|")
        println(io, "|E |$(glass.C5)|")
        println(io, "|F |$(glass.C6)|")
    elseif (D == -1)
        println(io, "|C₀|$(glass.C1)|")
        println(io, "|C₁|$(glass.C2)|")
        println(io, "|C₂|$(glass.C3)|")
        println(io, "|C₃|$(glass.C4)|")
    elseif (D == 1) # Schott
        println(io, "|a₀|$(glass.C1)|")
        println(io, "|a₁|$(glass.C2)|")
        println(io, "|a₂|$(glass.C3)|")
        println(io, "|a₃|$(glass.C4)|")
        println(io, "|a₄|$(glass.C5)|")
        println(io, "|a₅|$(glass.C6)|")
    elseif (D == 2)  # Sellmeier1
        println(io, "|K₁|$(glass.C1)|")
        println(io, "|L₁|$(glass.C2)|")
        println(io, "|K₂|$(glass.C3)|")
        println(io, "|L₂|$(glass.C4)|")
        println(io, "|K₃|$(glass.C5)|")
        println(io, "|L₃|$(glass.C6)|")
    elseif (D == 3)  # Herzberger
        println(io, "|A |$(glass.C1)|")
        println(io, "|B |$(glass.C2)|")
        println(io, "|C |$(glass.C3)|")
        println(io, "|D |$(glass.C4)|")
        println(io, "|E |$(glass.C5)|")
        println(io, "|F |$(glass.C6)|")
    elseif (D == 4)  # Sellmeier2
        println(io, "|A |$(glass.C1)|")
        println(io, "|B₁|$(glass.C2)|")
        println(io, "|λ₁|$(glass.C3)|")
        println(io, "|B₂|$(glass.C4)|")
        println(io, "|λ₂|$(glass.C5)|")
    elseif (D == 5)  # Conrady
        println(io, "|n₀|$(glass.C1)|")
        println(io, "|A |$(glass.C2)|")
        println(io, "|B |$(glass.C3)|")
    elseif (D == 6)  # Sellmeier3
        println(io, "|K₁|$(glass.C1)|")
        println(io, "|L₁|$(glass.C2)|")
        println(io, "|K₂|$(glass.C3)|")
        println(io, "|L₂|$(glass.C4)|")
        println(io, "|K₃|$(glass.C5)|")
        println(io, "|L₃|$(glass.C6)|")
        println(io, "|K₄|$(glass.C7)|")
        println(io, "|L₄|$(glass.C8)|")
    elseif (D == 7) || (D == 8)  # HandbookOfOptics1/2
        println(io, "|A |$(glass.C1)|")
        println(io, "|B |$(glass.C2)|")
        println(io, "|C |$(glass.C3)|")
        println(io, "|D |$(glass.C4)|")
    elseif (D == 9)  # Sellmeier4
        println(io, "|A |$(glass.C1)|")
        println(io, "|B |$(glass.C2)|")
        println(io, "|C |$(glass.C3)|")
        println(io, "|D |$(glass.C4)|")
        println(io, "|E |$(glass.C5)|")
    elseif (D == 10) || (D == 12) # Extended1/2
        println(io, "|a₀|$(glass.C1)|")
        println(io, "|a₁|$(glass.C2)|")
        println(io, "|a₂|$(glass.C3)|")
        println(io, "|a₃|$(glass.C4)|")
        println(io, "|a₄|$(glass.C5)|")
        println(io, "|a₅|$(glass.C6)|")
        println(io, "|a₆|$(glass.C7)|")
        println(io, "|a₇|$(glass.C8)|")
    elseif (D == 11)  # Sellmeier5
        println(io, "|K₁|$(glass.C1)|")
        println(io, "|L₁|$(glass.C2)|")
        println(io, "|K₂|$(glass.C3)|")
        println(io, "|L₂|$(glass.C4)|")
        println(io, "|K₃|$(glass.C5)|")
        println(io, "|L₃|$(glass.C6)|")
        println(io, "|K₄|$(glass.C7)|")
        println(io, "|L₄|$(glass.C8)|")
        println(io, "|K₅|$(glass.C9)|")
        println(io, "|L₅|$(glass.C10)|")
    elseif (D == 13)  # Extended3
        println(io, "|a₀|$(glass.C1)|")
        println(io, "|a₁|$(glass.C2)|")
        println(io, "|a₂|$(glass.C3)|")
        println(io, "|a₃|$(glass.C4)|")
        println(io, "|a₄|$(glass.C5)|")
        println(io, "|a₅|$(glass.C6)|")
        println(io, "|a₆|$(glass.C7)|")
        println(io, "|a₇|$(glass.C8)|")
        println(io, "|a₈|$(glass.C9)|")
    else
        println(io, "     INVALID DISPERSION FORMULA!!\n")
    end
    println(io)

    println(io, "**Valid wavelengths**$(glass.λmin)μm to $(glass.λmax)μm\n")
    println(io, "**Reference temperature**$(glass.temp)°C\n")

    if !isnan(glass.D₀) && (glass.D₀ != 0 || glass.D₁ != 0 || glass.D₂ != 0 || glass.E₀ != 0 || glass.E₁ != 0)
        println(io, "**Thermal ΔRI coefficients**\n")
        println(io, "|Coefficient                 |Value      |")
        println(io, "|:---------------------------|-----------|")
        println(io, "|D₀                          |$(glass.D₀)|")
        println(io, "|D₁                          |$(glass.D₁)|")
        println(io, "|D₂                          |$(glass.D₂)|")
        println(io, "|E₀                          |$(glass.E₀)|")
        println(io, "|E₁                          |$(glass.E₁)|")
        println(io, "|λₜₖ                          |$(glass.λₜₖ)|")
        println(io)
    end

    println(io, "**TCE (÷1e-6)**$(glass.TCE)\n")
    println(io, "**Ignore thermal expansion**$(glass.ignore_thermal_exp == 1)\n")

    println(io, "**Density (p)**$(glass.p)g/m³\n")
    println(io, "**ΔPgF**$(glass.ΔPgF)\n")

    println(io, "**RI at sodium D-Line (587nm)**$(glass.Nd)\n")
    println(io, "**Abbe Number**$(glass.Vd)\n")

    println(io, "**Cost relative to N_BK7**$(glass.relcost == -1 ? "?" : glass.relcost)\n")

    if glass.CR != -1 || glass.FR != -1 || glass.SR != -1 || glass.AR != -1 || glass.PR != -1
        println(io, "**Environmental resistance**\n")
        println(io, "|Environment   |Resistance                        |")
        println(io, "|:-------------|----------------------------------|")
        println(io, "|Climate (CR)  |$(glass.CR == -1 ? "?" : glass.CR)|")
        println(io, "|Stain (FR)    |$(glass.FR == -1 ? "?" : glass.FR)|")
        println(io, "|Acid (SR)     |$(glass.SR == -1 ? "?" : glass.SR)|")
        println(io, "|Alkaline (AR) |$(glass.AR == -1 ? "?" : glass.AR)|")
        println(io, "|Phosphate (PR)|$(glass.PR == -1 ? "?" : glass.PR)|")
        println(io)
    end

    println(io, "**Status**$(STATUS[glass.status + 1]) ($(glass.status))\n")
    println(io, "**Melt frequency**$(glass.meltfreq == -1 ? "?" : glass.meltfreq)\n")
    println(io, "**Exclude substitution**$(glass.exclude_sub == 1)\n")

    if glass.transmission !== nothing
        println(io, "**Transmission data**\n")
        println(io, "|Wavelength|Transmission|Thickness|")
        println(io, "|----------|------------|---------|")
        for i in 1:(glass.transmissionN)
            λ, t, τ = glass.transmission[i]
            println(io, "|$(λ)μm|$t|$(τ)mm|")
        end
    end
end