text:
using ExperimentalSetup,Gtk.ShortNames, GtkReactive, Reactive, DataStructures
const ES = ExperimentalSetup
function create_log(f=string.(1:5))
    factors = [Factor(f[i], f[1:i]) for i in linearindices(f)]
    a = ES.Log(ES.Metadata(factors))
    push!(a, repeat(f[1,:], outer=[length(f)]), "comment")
    push!(a, f, "comment")
end

_unique_factor(factors::Vector{Factor}, fname::String) = all(fname ≠ y.name for y in factors)

function edit_metadata!(a::ExperimentalSetup.Log, builder)
    update = Signal(nothing)
    backup = deepcopy(a)
    foreach(update) do _
        empty!(builder["metadata.box"])
        for (i, f) in enumerate(a.md.factors)
            b = Box(:h)
            name = textbox(f.name)
            foreach(name, init=nothing) do fname
                if isempty(fname)
                    delete!(a, f)
                elseif _unique_factor(a.md.factors, fname)
                    a.md.factors[i] = Factor(fname, f.levels)
                end
                push!(update, nothing)
            end
            push!(b, name)
            push!(b, label(": "))
            for (j, l) in enumerate(f.levels)
                name = textbox(l)
                foreach(name, init=nothing) do lname
                    if isempty(lname)
                        delete!(a, f, l)
                    elseif lname ∉ f.levels
                        f.levels[j] = lname
                    end
                    push!(update, nothing)
                end
                push!(b, name)
            end
            add_level = textbox("")
            foreach(add_level, init=nothing) do lname
                lname ∉ f.levels && push!(f.levels, lname)
                push!(update, nothing)
            end
            push!(b, add_level)
            push!(builder["metadata.box"], b)
        end
        add_factor = textbox("")
        foreach(add_factor, init=nothing) do fname
            _unique_factor(a.md.factors, fname) && push!(a.md.factors, Factor(fname, ["level 1"]))
            push!(update, nothing)
        end
        push!(builder["metadata.box"], add_factor)
        showall(builder["metadata.box"])
    end

    ok = button(widget = builder["metadata.ok"]) do _
        populate_collect!(a.md.factors, builder)
        visible(builder["metadata.window"], false)
    end

    cancel = button(widget = builder["metadata.cancel"], init=nothing) do _
        a = deepcopy(backup)
        # push!(update, nothing)
        visible(builder["metadata.window"], false)
        nothing
    end

    return a
end

function populate_collect!(factors::Vector{Factor}, builder)
    empty!(builder["collect.run.label.box"])
    empty!(builder["collect.run.combobox"])
    for f in factors
        push!(builder["collect.run.label.box"], Label(f.name))
        cb = ComboBoxText()
        for level in f.levels
            push!(cb, level)
        end
        setproperty!(cb, :active, 0)
        push!(builder["collect.run.combobox"], cb)
    end
    return builder
end

function collect_run(builder, levels::Vector{Int}, comment::String)
    for (i, level) in enumerate(levels)
        setproperty!(builder["collect.run.combobox"][i], :active, level - 1)
    end
    setproperty!(builder["collect.run.comment"], :text, comment)
    canceled = false
    c = Condition()
    okh = signal_connect(builder["collect.run.ok"], :clicked) do _
        notify(c)
    end
    cancelh = signal_connect(builder["collect.run.cancel"], :clicked) do _
        canceled = true
        notify(c)
    end
    visible(builder["collect.run.window"], true)
    wait(c)
    visible(builder["collect.run.window"], false)
    levels_ = [Gtk.GLib.bytestring(Gtk.GAccessor.active_text(cb)) for cb in builder["collect.run.combobox"]]
    comment_ = getproperty(builder["collect.run.comment"], :text, String)
    return (levels_, comment_, canceled)
end

const run_row = """
<interface>
<requires lib="gtk+" version="3.20"/>
<object class="GtkImage" id="delete.image">
<property name="visible">True</property>
<property name="can_focus">False</property>
<property name="stock">gtk-delete</property>
</object>
<object class="GtkImage" id="edit.image">
<property name="visible">True</property>
<property name="can_focus">False</property>
<property name="stock">gtk-edit</property>
</object>
<object class="GtkBox" id="run.row.box">
<property name="visible">True</property>
<property name="can_focus">False</property>
<child>
<object class="GtkLabel" id="label.run">
<property name="visible">True</property>
<property name="can_focus">False</property>
<property name="label" translatable="yes">None</property>
<property name="ellipsize">middle</property>
<property name="width_chars">20</property>
<property name="max_width_chars">20</property>
<property name="xalign">0</property>
</object>
<packing>
<property name="expand">False</property>
<property name="fill">True</property>
<property name="position">0</property>
</packing>
</child>
<child>
<object class="GtkButton" id="delete.run">
<property name="visible">True</property>
<property name="can_focus">True</property>
<property name="receives_default">True</property>
<property name="tooltip_markup" translatable="yes"><b>Delete</b> this run!</property>
<property name="image">delete.image</property>
<property name="always_show_image">True</property>
</object>
<packing>
<property name="expand">False</property>
<property name="fill">True</property>
<property name="position">1</property>
</packing>
</child>
<child>
<object class="GtkButton" id="edit.run">
<property name="visible">True</property>
<property name="can_focus">True</property>
<property name="receives_default">True</property>
<property name="image">edit.image</property>
<property name="always_show_image">True</property>
</object>
<packing>
<property name="expand">False</property>
<property name="fill">True</property>
<property name="position">2</property>
</packing>
</child>
</object>
</interface>
"""

function edit_runs!(aᵗ::Signal{ExperimentalSetup.Log}, builder)
    foreach(aᵗ) do a
        empty!(builder["runs.box"])
        for (k, v) in a.reps
            levels, _ = ES.pop(a, k)
            l = string(join(levels, ","), ":", v.replicate)
            row = Builder(buffer = run_row)
            Gtk.GAccessor.text(row["label.run"], l)
            deleteh = button(widget = row["delete.run"])
            foreach(deleteh, init=nothing) do _
                delete!(a, k)
                push!(aᵗ, a)
                nothing
            end
            edith = button(widget = row["edit.run"])
            foreach(edith, init=nothing) do _
                levels, comment, canceled = collect_run(builder, a.md.setups[a.reps[k].setup].levels, a.reps[k].comment)
                canceled || replace!(a, k, levels, comment)
                push!(aᵗ, a)
                nothing
            end
            push!(builder["runs.box"], row["run.row.box"])
        end
        addrun = button("+")
        foreach(addrun, init=nothing) do _
            levels_int = a.last == 0 ? ones(Int, a.md.nfactors) : a.md.setups[a.reps[end].setup].levels
            levels, comment, canceled = collect_run(builder, levels_int, "")
            canceled || push!(a, levels, comment)
            push!(aᵗ, a)
            nothing
        end
        push!(builder["runs.box"], addrun)

        # push!(aᵗ, a)
        showall(builder["runs.box"])
    end

    c = Condition()
    editmetadatah = signal_connect(builder["edit.metadata"], :clicked) do _
        notify(c)
    end
    # visible(builder["main.window"], true)
    wait(c)
    # visible(builder["main.window"], false)

    return value(aᵗ)
end

function main(a::ExperimentalSetup.Log)
    builder = Builder(filename=joinpath("/home/yakir/.julia/v0.6/ExperimentalSetup/src", "main.glade"))
    aᵗ = Signal(a)
    edit_metadata!(aᵗ, builder)
    edit_runs!(aᵗ, builder)
    # bᵗ = map(x -> edit_metadata!(x, builder), aᵗ)
    # cᵗ = map(x -> edit_runs!(x, builder), bᵗ)
    # bind!(aᵗ, cᵗ)

    showall(builder["main.window"])
    showall(builder["metadata.window"])
    visible(builder["metadata.window"], false)
    showall(builder["collect.run.window"])
    visible(builder["collect.run.window"], false)
end

a = create_log()
main(a)
# factors = [Factor("a", ["1","1000000000"]), Factor("b", ["1","100000000000000"])]
# a = ES.Log(ES.Metadata(factors))
# push!(a, ["1","1"], "")
# push!(a, ["1000000000","100000000000000"], "")

# function that takes an a, sends it to the main window and waits, main window returns a and message, then this function decides what to do with the a: save? edit metadata? etc.. this function has no window.

# s1 = Signal(Void())
# h1 = signal_connect(builder["edit.metadata.menu"], :activate) do _
# push!(s1, nothing)
# end


text:
## load in the data
ENV["OMP_NUM_THREADS"] = 6
using PowerSpectra
using Healpix
using PyCall, PyPlot
using CSV, DataFrames, LinearAlgebra
# using BenchmarkTools
hp = pyimport("healpy")
nmt = pyimport("pymaster")

##
data_dir = "/home/zequnl/.julia/dev/PowerSpectra/test/data/"
# mask = readMapFromFITS(data_dir * "mask.fits", 1, Float64)
nside = 256
lmax = 3 * nside - 1

np = pyimport("numpy")
nawrapper_dir = "/home/zequnl/Projects/nawrapper/notebooks/"
cls = np.genfromtxt("$(nawrapper_dir)/data/example_cls.txt", unpack=true)
ell, cltt, clte, clee, clbb = cls[1,:], cls[2,:], cls[3,:], cls[4,:], cls[5,:]
nw = pyimport("nawrapper")
nltt = nw.get_Nl(theta_fwhm=20.0, sigma_T=150, l_max=lmax+1)
nlee = nltt ./ 100

n1 = hp.synfast(
        (nltt, nlee, nlee,
         np.zeros_like(cltt), np.zeros_like(cltt), np.zeros_like(cltt)),
        nside, verbose=false, pixwin=false, new=true)

##
flat_beam = SpectralVector(ones(3*nside))
flat_mask = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )
mask1 = readMapFromFITS("$(data_dir)/example_mask_1.fits", 1, Float64)
mask2 = readMapFromFITS("$(data_dir)/example_mask_2.fits", 1, Float64)

##
flat_beam = SpectralVector(ones(3*nside))
flat_mask = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )
m1 = CovField("143_hm1", mask1, mask1, flat_mask, flat_mask, flat_mask, flat_beam, flat_beam)
m2 = CovField("143_hm2", mask2, mask2, flat_mask, flat_mask, flat_mask, flat_beam, flat_beam)
workspace = SpectralWorkspace(m1, m2)

npix = nside2npix(nside)
map1_ = PolarizedHealpixMap{Float64, RingOrder, Vector{Float64}}(
    n1[1,:],
    n1[2,:],
    n1[3,:])
map2_ = deepcopy(map1_)

map1, map2 = deepcopy(map1_), deepcopy(map2_)
mask!(map1, mask1, mask1)
mask!(map2, mask2, mask2)
a1 = map2alm(map1)
a2 = map2alm(map2)

M_EE = mcm(workspace, "EE", "143_hm1", "143_hm2").parent
M_EB = mcm(workspace, "EB", "143_hm1", "143_hm2").parent
M = M_EE - M_EB * inv(M_EE) * M_EB
M[1,1], M[2,2] = 1.0, 1.0
# Cl_hat = alm2cl(a1[2], a2[2], lu(parent(M)), flat_beam, flat_beam)
Cl_hat = alm2cl(a1[2], a2[2], M)

##
(M₊ - M₋ * inv(M₊) * M₋) .- (M_EE - M_EB * inv(M_EE) * M_EB)

##
M₊ .- M_EE[3:end,3:end]
##

M₋ * inv(M₊) * M₋ .- M_EB[3:end,3:end] * inv(M_EE[3:end,3:end]) * M_EB[3:end,3:end]

##

# num_ell = size(M_EE,1)
# M22 = zeros(2num_ell, 2num_ell)

# M22[1:num_ell,1:num_ell] .= M_EE
# M22[num_ell+1:2num_ell,num_ell+1:2num_ell] .= M_EE
# M22[1:num_ell,num_ell+1:2num_ell] .= M_EB
# M22[num_ell+1:2num_ell,1:num_ell] .= M_EB

##
M22 = PowerSpectra.mcm22(workspace, m1, m2)
c_EE, c_BB = alm2cl((a1[2], a1[3]), (a2[2], a2[3]), M22)
# ĉ_EE = alm2cl(a1[2], a2[2])
# ĉ_BB = alm2cl(a1[3], a2[3])
# ctot = qr(M22, Val(true)) \ vcat(ĉ_EE, ĉ_BB)
# c_EE = ctot[1:num_ell]
# c_BB = ctot[num_ell+1:2num_ell];

##

##

map1, map2 = deepcopy(map1_), deepcopy(map2_)
b = nmt.NmtBin.from_nside_linear(nside, 1)
# f_1 = nmt.NmtField(mask1.pixels, [map1_.i.pixels])
# f_2 = nmt.NmtField(mask2.pixels, [map2_.i.pixels])
f_1 = nmt.NmtField(mask1.pixels, [map1.q.pixels, map1.u.pixels])
f_2 = nmt.NmtField(mask2.pixels, [map2.q.pixels, map2.u.pixels])

w = nmt.NmtWorkspace()
w.compute_coupling_matrix(f_1, f_2, b)
cl_coupled = nmt.compute_coupled_cell(f_1, f_2)
# cl_coupled[2:end, :] .= 0.0
cl_decoupled = w.decouple_cell(cl_coupled)

cl_coupled[2:end, :] .= 0.0
cl_decoupled_EEONLY = w.decouple_cell(cl_coupled)

##
clf()
plot(( c_EE[3:end] ./ cl_decoupled[1,:]), "-")
plot(( c_BB[3:end] ./ cl_decoupled[4,:]), "-")
xlim(0,150)
# plot(cl_decoupled[4,:])
# xlim(0,150)
# ylim(0,0.0001)
# plot(nlee)
# plot(cl_decoupled)
# yscale("log")
gcf()

##

##
clf()
# plot(cl_coupled[1,3:end] ./ Cl_hat[3:end] )

plot(aps_decoupled[1:end] ./ cl_decoupled_EEONLY[1,1:end], "-" )
# ylim(0,2)
xlim(0,100)
# ylim(-1e-2, 1e-2)
# plot(cl_nmt[1,:])
# plot(Cl_hat[3:end])
# xlim(0,100)

# yscale("log")
gcf()

##

M_EE = mcm(workspace, "EE", "143_hm1", "143_hm2")
M_EB = mcm(workspace, "EB", "143_hm1", "143_hm2")

##
# mcm_aps = parent(M)
M₊ = w.get_coupling_matrix()[1:4:4lmax+1, 1:4:4lmax+1][3:end,3:end]
M₋ = w.get_coupling_matrix()[1:4:4lmax+1, 4:4:4lmax+4][3:end,3:end]

MCM_nmt = M₊ - M₋ * inv(M₊) * M₋
##

# maximum(abs.(M₊ .- M_EE.parent[3:end,3:end]))
clf()
plt.imshow(M₊ ./ M_EE.parent[3:end,3:end])
gcf()

##
clf()
cl_bespoke = deepcopy(cl_coupled[1,3:end])
# linalg = pyimport("scipy.linalg")
# cl_bespoke = linalg.lu_solve(linalg.lu_factor(
#     mcm_nmt[3:end,3:end]
#     ), b=cl_bespoke)

@time ldiv!(lu(MCM_nmt), cl_bespoke)

# plot(cl_bespoke)
# plot( cl_decoupled_EEONLY[1,:])

# ylim(0,0.00004)
plot(cl_bespoke ./ cl_decoupled_EEONLY[1,:], "-")
# ylim(0.8,1.2)
xlim(0,30)
gcf()

##
clf()
plot(diag(M₋,4) ./ diag(M₊,4))
gcf()
##

ext: jl

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
using HttpServer
using WebSockets

const connections = Dict{Int,WebSocket}()
const usernames   = Dict{Int,String}()

function decodeMessage( msg )
    String(copy(msg))
end

wsh = WebSocketHandler() do req, client
    global connections
    @show connections[client.id] = client
    println("req is $req")
    notifyonline = "Connection from user number $(client.id) is now online."
    for (k,v) in connections
        if k != client.id
            try
                write(v, notifyonline)
            catch
                continue
            end
        end
    end
    while true
        try
            msg = read(client)
        catch
            telloffline = "User $(usernames[client.id]) disconnected."
            println(telloffline, "(The client id was $(client.id).)")
            delete!(connections, client.id)
            if haskey(usernames, client.id)
                delete!(usernames, client.id)
            end
            for (k,v) in connections
                try
                    write(v, telloffline)
                catch
                    continue
                end
            end
            return
        end
        msg = decodeMessage(msg)
        if startswith(msg, "setusername:")
            println("SETTING USERNAME: $msg")
            usernames[client.id] = msg[13:end]
            notifyusername = "User number $(client.id) chose $(usernames[client.id]) as name handle."
            for (k,v) in connections
                try
                    write(v, notifyusername)
                catch
                    println("Caught exception writing to user $k")
                    continue
                end
            end
        end
        if startswith(msg, "say:")
            println("EMITTING MESSAGE: $msg")
            for (k,v) in connections
                if k != client.id
                    try
                        write(v, (usernames[client.id] * ": " * msg[5:end]))
                    catch
                        println("Caught exception writing to user $k")
                        continue
                    end
                end
            end
        end
    end
end

onepage = readstring(Pkg.dir("WebSockets","examples","chat-client.html"))
httph = HttpHandler() do req::Request, res::Response
  Response(onepage)
end

server = Server(httph, wsh)
println("Chat server listening on 8000...")
run(server,8000)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Distributions, ImageCore, MappedArrays, Colors, FixedPointNumbers

function gauss(mean, var)
    function f0(I...)
        pdf(MvNormal(mean, sqrt(var)), [I...])
    end
    f0
end

function gauss(domain::Domain, var = 0.02)
    gauss([(size(domain)./2)...], var)
end

const tocolor = colorsigned(RGBA{N0f8}(174/255,55/255,14/255,1.),RGBA{N0f8}(.99,.99,.99,0.),RGBA{N0f8}(5/255,105/255,97/255,1.))

function toimage(state::State, scale = 2)
    mappedarray(i->tocolor(scalesigned(scale)(i)), state.current)
end

function toimage(img, scale = 2)
    mappedarray(i->tocolor(scalesigned(scale)(i)), img)
end

function dc_block(signal; R=0.995)
    cor_out = zeros(signal)
    cor_out[1] = signal[1]
    for i in 2:length(signal)
        cor_out[i] = signal[i] - signal[i-1] + R * cor_out[i-1]
    end
    cor_out
end

function make_sine(dt; A=1, f=220, duration=2., phase=0)
    Float64[A*sin(2π*f*t+phase) for t in 0:dt:duration]
end

function sin_response(; A=1, f=220, phase=0, length = 0.5, showvalues=true, duration=2., resource=CPUThreads((100,100,1)), gamma=0.05, fps=5, kw...)
    wave = UniformWave{3}(; kw...)
    signal = make_sine(wave.dt, A=A, f=f, phase=phase, duration = length)
    signal_response(signal, showvalues=showvalues, duration=duration, resource=resource, gamma=gamma, kw...)
end

function signal_response(signal; showvalues=true, fps=5, duration=2., resource=CPUThreads((100,100,1)), gamma=0.05, kw...)
    wave = UniformWave{3}(;kw...)
    speaker = PointSpeaker((3,1,2), signal)
    micro = PointMicrophone(3,7,2)
    sim = Simulator(wave, StoreSnapshots(1/fps), speaker, micro, ProgressDisplay(0.1, showvalues), resource=resource, duration=duration)
    domain = BoxDomain(6,8,4, gamma=gamma)
    backend = backend_init(sim.resource, domain, sim)
    state = state_init(backend, domain, sim);
    simulate!(state, backend, sim)
    backend_cleanup!(backend, state, sim)
    speaker.signal, micro.signal, sim
end

function impulse_response(; gamma=0., duration=2., resource=CPUThreads((100,100,1)), kw...)
    micro_in  = PointMicrophone(4,5,3)
    micro_out = PointMicrophone(1.7,2.8,0.9)
    sim = Simulator(UniformWave{3}(;kw...), micro_in, micro_out, ProgressDisplay(), resource=resource, duration=duration)
    domain = BoxDomain(6,8,4, gamma=gamma)
    backend = backend_init(sim.resource, domain, sim)
    state = state_init(backend, domain, sim);
    gridpos = floor.(Int, (4,5,3) ./ sim.wave.dx) .+ 1
    state.current[gridpos...] = sim.wave.c^2 * sim.wave.dt / sim.wave.dx^3
    simulate!(state, backend, sim)
    backend_cleanup!(backend, state, sim)
    micro_in.signal, micro_out.signal, sim
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ ccdc3522-c45c-11eb-152a-93999f4d6367
begin
    import Pkg
    # activate a clean environment
    Pkg.activate(mktempdir())

    Pkg.add([
        Pkg.PackageSpec(name="HTTP"),
        Pkg.PackageSpec(name="JSON3"),
        # ... keep adding your packages
    ])

    using HTTP
    using JSON3
    # ... place all usings and imports into this one cell
end

# ╔═╡ 20fc8448-c45f-11eb-0788-c9b8e5855477
begin
	j = JSON3.read(HTTP.get("https://pokeapi.co/api/v2/type/3").body)
#	j_dict = JSON.parse(j::AbstractString; dicttype=Dict, inttype=Int64)
end

# ╔═╡ c1d141bc-2b24-402e-adad-988b98d5b5b4
j

# ╔═╡ 22dd6dc2-c45f-11eb-098c-eb59381196f6
all_type3_pokemon_names = [x["pokemon"]["name"] for x in j[:pokemon]]

# ╔═╡ 460d14b8-c460-11eb-1879-9ba2f5fe8bfc
begin
	for x in j[:pokemon]
		poke_name = x["pokemon"]["name"]
		print(poke_name)
		poke_url = x["url"]
	end
end

# ╔═╡ Cell order:
# ╟─ccdc3522-c45c-11eb-152a-93999f4d6367
# ╠═20fc8448-c45f-11eb-0788-c9b8e5855477
# ╠═c1d141bc-2b24-402e-adad-988b98d5b5b4
# ╠═22dd6dc2-c45f-11eb-098c-eb59381196f6
# ╠═460d14b8-c460-11eb-1879-9ba2f5fe8bfc

max_stars_repo_name: alvarmaciel/paginas-pluto

max_issues_repo_name: alvarmaciel/paginas-pluto

max_forks_repo_name: alvarmaciel/paginas-pluto


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04pu(datfile, io=stdout)
    ZERO = 0.0e0
    ONE = 1.0e0
    TWO = 2.0e0
    NIN = 5
    NOUT = 6
    NMAX = 100
    LDA  = NMAX
    LDQG = NMAX
    LDRES  = NMAX
    LDU1 = NMAX
    LDU2 = NMAX
    LDWORK = 2*NMAX
    A = Array{Float64,2}(undef, LDA,NMAX)
    QG = Array{Float64,2}(undef, LDQG,NMAX+1)
    CS = Array{Float64,1}(undef, 2*NMAX)
    TAU = Array{Float64,1}(undef, NMAX)
    U1 = Array{Float64,2}(undef, LDU1,NMAX)
    U2 = Array{Float64,2}(undef, LDU2,NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    RES = Array{Float64,2}(undef, LDRES,3*NMAX+1)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    if ( N<=0 || N>NMAX )
        @error "invalid N=$N"
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    #FOREIGN.dlacpy!( 'All', N, N, A, LDA, RES(1,N+1), LDRES )
    RES[1:N,N+1:2N] .= A[1:N,1:N]

    vs = String[]
    _isz,_jsz = (N,N+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       QG[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 2

    #FOREIGN.dlacpy!( 'All', N, N+1, QG, LDQG, RES(1,2*N+1), LDRES )
    RES[1:N,2*N+1:3N+1] .= QG[1:N,1:N+1]

    close(f)
# interp call 3

    INFO = SLICOT.mb04pu!(N, 1, A, QG, CS, TAU, LDWORK)
    @test INFO == 0
    INFO == 0 || return

    # yes, this is from the original
    # INFO = 0
    # if ( INFO!=0 )
    # else
# interp call 4

    #FOREIGN.dlacpy!( 'Lower', N, N, A, LDA, U1, LDU1 )
    U1[1:N,1:N] .= tril(A[1:N,1:N])

# interp call 5

    #FOREIGN.dlacpy!( 'Lower', N, N, QG, LDQG, U2, LDU2 )
    U2[1:N,1:N] .= tril(QG[1:N,1:N])

# interp call 6

    INFO = SLICOT.mb04wp!(N, 1, U1, U2, CS, TAU)
    @test INFO == 0
    INFO == 0 || return

    if  N>2
# interp call 7

        # FOREIGN.dlaset!( 'Lower', N-2, N-2, ZERO, ZERO, A(3,1), LDA )
        triu!(view(A,2:N,1:N-2))

    end
    if  N>1
# interp call 8

        # FOREIGN.dlaset!( 'Lower', N-1, N-1, ZERO, ZERO, QG(2,1), LDQG )
        triu!(view(QG,1:N,1:N-1))

    end
# interp output 1
    println(io, "U1:")
    _nc = N
    _nr = N
    show(io, "text/plain", U1[1:_nr,1:_nc])
    println(io)
    println(io, "U2:")
    _nc = N
    _nr = N
    show(io, "text/plain", U2[1:_nr,1:_nc])
    println(io)

# unable to translate write statement:
#  write MA02JD( .FALSE., .FALSE., N, U1, LDU1, U2, LDU2, RES, LDRES )
    resnorm = SLICOT.ma02jd!(false,false,N,U1,U2)
    println(io, "U orth. residual norm = $resnorm")
    @test resnorm < 1e-12
# interp output 2
    println(io, "A:")
    _nc = N
    _nr = N
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

# interp output 3
    println(io, "QG:")
    _nc = N+1
    _nr = N
    show(io, "text/plain", QG[1:_nr,1:_nc])
    println(io)

# interp call 9

    #FOREIGN.dgemm!( 'No Transpose', 'No Transpose', N, N, N, ONE, U1, LDU1, A, LDA, ZERO, RES, LDRES )
    BLAS.gemm!('N','N',ONE,U1[1:N,1:N],A[1:N,1:N],ZERO,view(RES,1:N,1:N))

# interp call 10

    #FOREIGN.dgemm!( 'No Transpose', 'Transpose', N, N, N, -ONE, RES, LDRES, U1, LDU1, ONE, RES(1,N+1), LDRES )
    BLAS.gemm!('N','T',-ONE,RES[1:N,1:N],U1[1:N,1:N],ONE,view(RES,1:N,N+1:2N))

# interp call 11

    #FOREIGN.dgemm!( 'No Transpose', 'Transpose', N, N, N, ONE, U2, LDU2, A, LDA, ZERO, RES, LDRES )
    BLAS.gemm!('N','T',ONE,U2[1:N,1:N],A[1:N,1:N],ZERO,view(RES,1:N,1:N))

# interp call 12

    #FOREIGN.dgemm!( 'No Transpose', 'Transpose', N, N, N, ONE, RES, LDRES, U2, LDU2, ONE, RES(1,N+1), LDRES )
    BLAS.gemm!('N','T',ONE,RES[1:N,1:N],U2[1:N,1:N],ONE,view(RES,1:N,N+1:2N))

# interp call 13

    #FOREIGN.dsymm!( 'Right', 'Upper', N, N, ONE, QG(1,2), LDQG, U1, LDU1, ZERO, RES, LDRES )
    BLAS.symm!('R','U',ONE,QG[1:N,2:N+1],U1[1:N,1:N],ZERO,view(RES,1:N,1:N))

# interp call 14

    #FOREIGN.dgemm!( 'No Transpose', 'Transpose', N, N, N, -ONE, RES, LDRES, U2, LDU2, ONE, RES(1,N+1), LDRES )
    BLAS.gemm!('N','T',-ONE,RES[1:N,1:N],U2[1:N,1:N],ONE,view(RES,1:N,N+1:2N))

    RES[1:N,1:N] .= U2[1:N,1:N]
# interp call 15

    for i in 1:N
        # FOREIGN.dscal!( N, QG(I,I), RES(1,I), 1 )
        RES[1:N,i] .*= QG[i,i]
    end

# interp call 16

    #FOREIGN.dgemm!( 'No Transpose', 'Transpose', N, N, N, -ONE, RES, LDRES, U1, LDU1, ONE, RES(1,N+1), LDRES )
    BLAS.gemm!('N','T',-ONE,RES[1:N,1:N],U1[1:N,1:N],ONE,view(RES,1:N,N+1:2N))

# interp call 17

    #FOREIGN.dgemm!( 'No Transpose', 'No Transpose', N, N, N, ONE, U2, LDU2, A, LDA, ZERO, RES, LDRES )
    BLAS.gemm!('N','N',ONE,U2[1:N,1:N],A[1:N,1:N],ZERO,view(RES,1:N,1:N))

# interp call 18

    #FOREIGN.dsyr2k!( 'Lower', 'No Transpose', N, N, ONE, RES, LDRES, U1, LDU1, ONE, RES(1,2*N+1), LDRES )
    BLAS.syr2k!('L','N',ONE,RES[1:N,1:N],U1[1:N,1:N],ONE,view(RES,1:N,2*N+1:3N))

# interp call 19
    #FOREIGN.dscal!( N, ONE/TWO, QG(1,2), LDQG+1 )
    for i in 1:N
        QG[i,i+1] *= 1/2
    end

    RES[1:N,1:N] .= U2[1:N,1:N]
# interp call 20

    #FOREIGN.dtrmm!(  'Right', 'Upper' , 'No Transpose', 'Not unit', N, N, ONE, QG(1,2), LDQG, RES, LDRES )
    BLAS.trmm!('R','U','N','N',ONE,QG[1:N,2:N+1],view(RES,1:N,1:N))

# interp call 21

    #FOREIGN.dsyr2k!( 'Lower', 'No Transpose', N, N, ONE, RES, LDRES, U2, LDU2, ONE, RES(1,2*N+1), LDRES )
    BLAS.syr2k!('L','N',ONE,RES[1:N,1:N],U2[1:N,1:N],ONE,view(RES,1:N,2N+1:3N))

# interp call 22

    for i in 1:N
    #FOREIGN.dsyr!( 'Lower', N, -QG(I,I), U1(1,I), 1, RES(1,2*N+1), LDRES )
        BLAS.syr!('L',-QG[i,i],U1[1:N,i],view(RES,1:N,2*N+1:3N))
    end

# interp call 23

    #FOREIGN.dgemm!( 'No Transpose', 'No Transpose', N, N, N, ONE, U1, LDU1, A, LDA, ZERO, RES, LDRES )
    BLAS.gemm!('N','N',ONE,U1[1:N,1:N],A[1:N,1:N],ZERO,view(RES,1:N,1:N))

# interp call 24

    #FOREIGN.dsyr2k!( 'Upper', 'No Transpose', N, N, ONE, RES, LDRES, U2, LDU2, ONE, RES(1,2*N+2), LDRES )
    BLAS.syr2k!('U','N',ONE,RES[1:N,1:N],U2[1:N,1:N],ONE,view(RES,1:N,2*N+2:3*N+1))

    RES[1:N,1:N] .= U1[1:N,1:N]
# interp call 25

    #FOREIGN.dtrmm!(  'Right', 'Upper' , 'No Transpose', 'Not unit', N, N, ONE, QG(1,2), LDQG, RES, LDRES )
    BLAS.trmm!('R','U','N','N',ONE,QG[1:N,2:N+1],view(RES,1:N,1:N))

# interp call 26

    #FOREIGN.dsyr2k!( 'Upper', 'No Transpose', N, N, -ONE, RES, LDRES, U1, LDU1, ONE, RES(1,2*N+2), LDRES )
    BLAS.syr2k!('U','N',-ONE,RES[1:N,1:N],U1[1:N,1:N],ONE,view(RES,1:N,2*N+2:3*N+1))

# interp call 27
    for i in 1:N
        #FOREIGN.dsyr!( 'Upper', N, QG(I,I), U2(1,I), 1, RES(1,2*N+2), LDRES )
        BLAS.syr!('U',QG[i,i],U2[1:N,i],view(RES,1:N,2*N+2:3*N+1))
    end

# unable to translate write statement:
#  write MA02ID( 'Hamiltonian', 'Frobenius', N, RES(1,N+1), LDRES, RES(1,2*N+1), LDRES, DWORK )
    resnorm = SLICOT.ma02id!('H','F',N,view(RES,1:N,N+1:2*N),view(RES,1:N,2*N+1:3*N),DWORK)
    println(io, "residual norm = $resnorm")
    @test resnorm < 1e-12
end # run_mb04pu()


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    combinations, sortstablemax, sortstablemin

## combinations

struct Combinations{T}
    collection::Vector{T}
    indice::Vector{Int}
end

function Base.iterate(cmb::Combinations, state=nothing)
    lastidx = lastindex(cmb.collection)
    idcsize = length(cmb.indice)
    if cmb.indice[end] == lastidx
        if idcsize == 1
            return
        end
        n = 1
        while cmb.indice[end - n] == lastidx - n
            n += 1
            if n == idcsize
                return
            end
        end
        cmb.indice[idcsize - n] += 1
        for i in (idcsize - n + 1):idcsize
            cmb.indice[i] = cmb.indice[i - 1] + 1
        end
    else
        cmb.indice[end] += 1
    end
    return (cmb.collection[cmb.indice], state)
end

Base.eltype(::Type{Combinations{T}}) where {T} = Vector{T}
Base.IteratorSize(::Type{<:Combinations}) = Base.SizeUnknown()

"""
    combinations(iter, k::Int=2)

An iterator that yields `k`-combinations (card. k subsets of the collection)
"""
function combinations(iter, k::Int=2)
    col = collect(iter)
    T = eltype(col)
    if k < 0
        throw(DomainError(k, "k should not be negative"))
    elseif k > length(col)
        throw(ErrorException("$(k) is larger than the collection size"))
    elseif k == 0
        return (Int[],)
    else
        idc = collect(1:k)
        idc[end] -= 1
        return Combinations{T}(col, idc)
    end
end

function sortstablemax(iter; by=identity, kwargs...)
    isempty(iter) && return kwargs[:init]
    cmp(x, y) = ifelse(by(x) < by(y), y, x)
    return reduce(cmp, iter)
end

function sortstablemin(iter; by=identity, kwargs...)
    isempty(iter) && return kwargs[:init]
    cmp(x, y) = ifelse(by(x) > by(y), y, x)
    return reduce(cmp, iter)
end


text:
"""
    getkeypress()::Union{CharKey, ModifiedKey, SpecialKey, UnknownKey

Return the next key pressed.
"""
getkeypress()::Union{CharKey, ModifiedKey, SpecialKey, UnknownKey} = bytestokey(getkeybytes())

"""
    prompt(s::AbstractString)::String

Temporarily enter canonical mode and turn on echo to prompt for a string.
"""
function prompt(s)::String
    prev_mode = normalmode()
    print(s)
    s = readline(stdin)
    setmode(prev_mode)

    s
end

# Internal

function getkeybytes()
    c = read(stdin, UInt8)
    (c == 0x1b || c > 0xc0) && more(stdin) ? getkeybytes(c) : c
end

function getkeybytes(x::UInt8)
    c = read(stdin, UInt8)
    (x == 0x1b || x == 0xe2) && more(stdin) ? getkeybytes(x, c) : (x, c)
end

function getkeybytes(x::UInt8, y::UInt8)
    c = read(stdin, UInt8)
    x == 0x1b && y == 0x5b && more(stdin) ? getkeybytes(x, y, c) : (x, y, c)
end

function getkeybytes(x::UInt8, y::UInt8, z::UInt8)
    c = read(stdin, UInt8)
    more(stdin) ? getkeybytes(x, y, z, c) : (x, y, z, c)
end

function getkeybytes(x::UInt8, y::UInt8, z::UInt8, p::UInt8)
    c = read(stdin, UInt8)
    p != 0x7e && more(stdin) ? getkeybytes(x, y, z, p, c) : (x, y, z, p, c)
end

function getkeybytes(x::UInt8, y::UInt8, z::UInt8, p::UInt8, q::UInt8)
    c = read(stdin, UInt8)
    q != 0x7e && more(stdin) ? getkeybytes(x, y, z, p, q, c) : (x, y, z, p, q, c)
end

function getkeybytes(x::UInt8, y::UInt8, z::UInt8, p::UInt8, q::UInt8, r::UInt8)
    c = read(stdin, UInt8)
    (x, y, z, p, q, r, c)
end

function bytestokey(x::UInt8)::Union{CharKey, ModifiedKey, SpecialKey, UnknownKey}
    if ischar(x)
        CharKey(Char(x))
    else
        @match x begin
            0x00             => ctrl | CharKey(' ')
            0x09             => CharKey('\t')
            0x0a             => returnKey
            0x1b             => escape
            0x1d             => ctrl | CharKey('5')
            0x1e             => ctrl | CharKey('6')
            0x1f             => ctrl | CharKey('7')
            0x7f             => backspace
            _, if x < 27 end => ctrl | CharKey(Char(0x60 + x))
            _                => UnknownKey(x)
        end
    end
end

function bytestokey(xs::Tuple{UInt8, UInt8})::Union{CharKey, ModifiedKey, UnknownKey}
    @match xs begin
        (0x1b, y)    => alt | bytestokey(y)
        (0xc2, y)    => CharKey(Char(y))
        (0xc3, y)    => CharKey(Char(0x40 + y))
        (0xd0, 0xb6) => ctrl | CharKey('ö')
        (0xd1, 0x85) => ctrl | CharKey('å')
        (0xd1, 0x8b) => ctrl | CharKey('ß')
        (0xd1, 0x8d) => ctrl | CharKey('ä')
        _            => UnknownKey(xs)
    end
end

function bytestokey(xs::Tuple{UInt8, UInt8, UInt8})::Union{CharKey, ModifiedKey, SpecialKey, UnknownKey}
    @match xs begin
        (0x1b, 0x4f, 0x50) => f1
        (0x1b, 0x4f, 0x51) => f2
        (0x1b, 0x4f, 0x52) => f3
        (0x1b, 0x4f, 0x53) => f4
        (0x1b, 0x5b, 0x41) => arrowUp
        (0x1b, 0x5b, 0x42) => arrowDown
        (0x1b, 0x5b, 0x43) => arrowRight
        (0x1b, 0x5b, 0x44) => arrowLeft
        (0x1b, 0x5b, 0x46) => endKey
        (0x1b, 0x5b, 0x48) => home
        (0x1b, 0xc3, x)    => maybemod(alt, (0xc3, x), xs)
        (0x1b, 0xd0, x)    => maybemod(alt, (0xd0, x), xs)
        (0x1b, 0xd1, x)    => maybemod(alt, (0xd1, x), xs)
        (0xe2, 0x82, 0xac) => CharKey('€')
        _ => UnknownKey(xs)
    end
end

function bytestokey(xs::Tuple{UInt8, UInt8, UInt8, UInt8})::Union{SpecialKey, UnknownKey}
    @match xs begin
        (0x1b, 0x5b, 0x32, 0x7e) => ins
        (0x1b, 0x5b, 0x33, 0x7e) => del
        (0x1b, 0x5b, 0x35, 0x7e) => pgUp
        (0x1b, 0x5b, 0x36, 0x7e) => pgDn
        _                        => UnknownKey(xs)
    end
end

function bytestokey(xs::Tuple{UInt8, UInt8, UInt8, UInt8, UInt8})::Union{SpecialKey, UnknownKey}
    @match xs begin
        (0x1b, 0x5b, 0x31, 0x35, 0x7e) => f5
        (0x1b, 0x5b, 0x31, 0x37, 0x7e) => f6
        (0x1b, 0x5b, 0x31, 0x38, 0x7e) => f7
        (0x1b, 0x5b, 0x31, 0x39, 0x7e) => f8
        (0x1b, 0x5b, 0x32, 0x30, 0x7e) => f9
        (0x1b, 0x5b, 0x32, 0x31, 0x7e) => f10
        (0x1b, 0x5b, 0x32, 0x33, 0x7e) => f11
        (0x1b, 0x5b, 0x32, 0x34, 0x7e) => f12
        _                              => UnknownKey(xs)
    end
end

function bytestokey(xs::Tuple{UInt8, UInt8, UInt8, UInt8, UInt8, UInt8})::Union{ModifiedKey, UnknownKey}
    @match xs begin
        (0x1b, 0x5b, 0x31, 0x3b, 0x32, x)    => maybemod(shift, (0x1b, 0x5b, x), xs)
        (0x1b, 0x5b, 0x31, 0x3b, 0x33, x)    => maybemod(alt, (0x1b, 0x5b, x), xs)
        (0x1b, 0x5b, 0x31, 0x3b, 0x35, x)    => maybemod(ctrl, (0x1b, 0x5b, x), xs)
        (0x1b, 0x5b, 0x31, 0x3b, 0x37, x)    => maybemod(alt | ctrl, (0x1b, 0x5b, x), xs)
        (0x1b, 0x5b, x, 0x3b, 0x35, 0x7e)    => maybemod(ctrl, (0x1b, 0x5b, x, 0x7e), xs)
        (0x1b, 0x5b, 0x31, 0x3b, 0x36, 0x53) => ctrl | shift | f4
        (0x1b, 0x5b, 0x31, 0x3b, 0x36, 0x52) => ctrl | shift | f3
        (0x1b, 0x5b, 0x31, 0x3b, 0x36, 0x51) => ctrl | shift | f2
        (0x1b, 0x5b, 0x31, 0x3b, 0x36, 0x50) => ctrl | shift | f1
        (0x1b, 0x5b, 0x33, 0x3b, 0x32, 0x7e) => shift | del
        _                                    => UnknownKey(xs)
    end
end

function bytestokey(xs::Tuple{UInt8, UInt8, UInt8, UInt8, UInt8, UInt8, UInt8})::Union{ModifiedKey, UnknownKey}
    @match xs begin
        (0x1b, 0x5b, x, y, 0x3b, 0x35, 0x7e) => maybemod(ctrl, (0x1b, 0x5b, x, y, 0x7e), xs)
        (0x1b, 0x5b, x, y, 0x3b, 0x36, 0x7e) => maybemod(ctrl | shift, (0x1b, 0x5b, x, y, 0x7e), xs)
        _                                    => UnknownKey(xs)
    end
end

"""
    maybemod(mask, xs, ys)

Modify the key corresponding to `xs` if it is not unknown. Otherwise create
unknown key from `ys`.
"""
function maybemod(mask, xs, ys)
    k = bytestokey(xs)
    typeof(k) == UnknownKey ? UnknownKey(ys) : mask | k
end

ischar(c::UInt8) = 32 <= c <= 126 || 161 <= c <= 255

function more(stream)
    t = @task begin; peek(stream, UInt8); end
    schedule(t)
    sleep(0.001)
    istaskdone(t)
end


text:
# Grids

"""
	$(SIGNATURES)

Abstract grid type.
"""
abstract type AbstractGrid end

"""
	$(SIGNATURES)

Linear, equally spaced grid.
"""
struct LinearGrid{T <: AbstractFloat, I <: Integer} <: AbstractGrid
    xMin :: T
    xMax :: T
    nPoints :: I
end

"""
	$(SIGNATURES)

Grid with increasing or decreasing intervals. Interval(j) / interval(j-1) is a constant.
Constructed from:
- xMin, xMax: range
- nPoints: no of grid points
- ratio of interval(j) / interval(j-1)
"""
struct PowerSpacedGrid{T <: AbstractFloat, I <: Integer} <: AbstractGrid
    xMin :: T
    xMax :: T
    nPoints :: I
    intervalRatio :: T
end

"""
	$(SIGNATURES)

Grid with intervals that increase or decrease by constants.
"""
struct LinSpacedGrid{T <: AbstractFloat, I <: Integer} <: AbstractGrid
    xMin :: T
    xMax :: T
    nPoints :: I
    intervalDiff :: T
end

## -------------  Common functions

Base.eltype(g :: AbstractGrid) = Base.eltype(minimum(g));
Base.length(g :: AbstractGrid) = g.nPoints;
Base.minimum(g :: AbstractGrid) = g.xMin;
Base.maximum(g :: AbstractGrid) = g.xMax;
n_intervals(g :: AbstractGrid) = length(g) - 1;
intervals(g :: AbstractGrid) = diff(grid(g));
Base.firstindex(g :: AbstractGrid) = 1;
Base.lastindex(g :: AbstractGrid) = length(g);
Base.getindex(g :: AbstractGrid, idx) = grid(g)[idx];

function grid(g :: AbstractGrid)
    return minimum(g) .+ vcat(zero(eltype(g)), cumsum(intervals(g)))
end

function Base.iterate(g :: AbstractGrid, idx)
    if idx > length(g)
        return nothing
    else
        return g[idx], idx+1
    end
end

Base.iterate(g :: AbstractGrid) = Base.iterate(g, 1);

## --------------  Linear Grid

grid(g :: LinearGrid) = range(minimum(g), maximum(g), length = length(g));

function Base.getindex(g :: LinearGrid, idx)
    dx = (maximum(g) - minimum(g)) / n_intervals(g);
    return minimum(g) .+ dx .* (collect(idx) .- 1)
end

## --------------  Power Spaced Grid

function first_interval(g :: PowerSpacedGrid{T,I}) where {T, I}
    factor1 = (g.intervalRatio ^ (length(g) - 1) - one(T)) / (g.intervalRatio - 1);
    return (maximum(g) - minimum(g)) / factor1
end

function intervals(g :: PowerSpacedGrid{T,I}) where {T, I}
    intV = Vector{T}(undef, length(g) - 1);
    intV[1] = first_interval(g);
    for j = 2 : (length(g) - 1)
        intV[j] = intV[j-1] * g.intervalRatio;
    end
    return intV
end

# function grid(g :: PowerSpacedGrid{T,I}) where {T,I}
#     return minimum(g) .+ vcat(zero(eltype(g)), cumsum(intervals(g)))
# end

## ----------------  Linearly changing intervals

function first_interval(g :: LinSpacedGrid{T,I}) where {T, I}
    n = length(g);
    dX = maximum(g) - minimum(g);

    return dX / (n-1) - g.intervalDiff * (n-2) / 2.0;
end

interval(g :: LinSpacedGrid{T,I}, idx) where {T, I} =
    first_interval(g) .+ (idx .- 1) .* g.intervalDiff;

intervals(g :: LinSpacedGrid{T,I}) where {T, I} =
    interval(g, 1 : n_intervals(g))

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


text:
const SKIP_PM_VALIDATION = false

"""
    System

    A power system defined by fields for basepower, components, and forecasts.

    # Constructor
    ```julia
    System(basepower)
    System(components, forecasts, basepower)
    System(buses, generators, loads, branches, storage, basepower, forecasts, services, annex; kwargs...)
    System(buses, generators, loads, basepower; kwargs...)
    System(file; kwargs...)
    System(; buses, generators, loads, branches, storage, basepower, forecasts, services, annex, kwargs...)
    System(; kwargs...)
    ```

    # Arguments

    * `buses`::Vector{Bus} : an array of buses
    * `generators`::Vector{Generator} : an array of generators of (possibly) different types
    * `loads`::Vector{ElectricLoad} : an array of load specifications that includes timing of the loads
    * `branches`::Union{Nothing, Vector{Branch}} : an array of branches; may be `nothing`
    * `storage`::Union{Nothing, Vector{Storage}} : an array of storage devices; may be `nothing`
    * `basepower`::Float64 : the base power value for the system
    * `forecasts`::Union{Nothing, IS.Forecasts} : dictionary of forecasts
    * `services`::Union{Nothing, Vector{ <: Service}} : an array of services; may be `nothing`

    # Keyword arguments

    * `runchecks`::Bool : run available checks on input fields. If an error is found in a field, that component will not be added to the system and InvalidRange is thrown.
    * `configpath`::String : specify path to validation config file
    DOCTODO: any other keyword arguments? genmap_file, REGEX_FILE
"""
struct System <: PowerSystemType
    data::IS.SystemData
    basepower::Float64             # [MVA]
    runchecks::Bool
    internal::InfrastructureSystemsInternal

    function System(data, basepower, internal; kwargs...)
        runchecks = get(kwargs, :runchecks, true)
        sys = new(data, basepower, runchecks, internal)
    end
end

"""Construct an empty System. Useful for building a System while parsing raw data."""
function System(basepower; kwargs...)
    return System(_create_system_data_from_kwargs(; kwargs...), basepower)
end

function System(data, basepower; kwargs...)
    return System(data, basepower, InfrastructureSystemsInternal(); kwargs...)
end

"""System constructor when components are constructed externally."""
function System(buses::Vector{Bus},
                generators::Vector{<:Generator},
                loads::Vector{<:ElectricLoad},
                branches::Union{Nothing, Vector{<:Branch}},
                storage::Union{Nothing, Vector{<:Storage}},
                basepower::Float64,
                forecasts::Union{Nothing, IS.Forecasts},
                services::Union{Nothing, Vector{ <: Service}},
                annex::Union{Nothing,Dict}; kwargs...)

    data = _create_system_data_from_kwargs(; kwargs...)

    if isnothing(forecasts)
        forecasts = IS.Forecasts()
    end

    sys = System(data, basepower; kwargs...)

    arrays = [buses, generators, loads]
    if !isnothing(branches)
        push!(arrays, branches)
    end
    if !isnothing(storage)
        push!(arrays, storage)
    end
    if !isnothing(services)
        push!(arrays, services)
    end

    error_detected = false

    for component in Iterators.flatten(arrays)
        try
            add_component!(sys, component)
        catch e
            if isa(e, InvalidRange)
                error_detected = true
            else
                rethrow()
            end
        end
    end

    load_zones = isnothing(annex) ? nothing : get(annex, :LoadZones, nothing)
    if !isnothing(load_zones)
        for lz in load_zones
            try
                add_component!(sys, lz)
            catch e
                if isa(e, InvalidRange)
                    error_detected = true
                else
                    rethrow()
                end
            end
        end
    end

    runchecks = get(kwargs, :runchecks, true)

    if error_detected
        throw(InvalidRange("Invalid value(s) detected"))
    end

    if runchecks
        check!(sys)
    end

    return sys
end

"""System constructor without nothing-able arguments."""
function System(buses::Vector{Bus},
                generators::Vector{<:Generator},
                loads::Vector{<:ElectricLoad},
                basepower::Float64; kwargs...)
    return System(buses, generators, loads, nothing, nothing, basepower, nothing, nothing, nothing; kwargs...)
end

"""System constructor with keyword arguments."""
function System(; basepower=100.0,
                buses,
                generators,
                loads,
                branches,
                storage,
                forecasts,
                services,
                annex,
                kwargs...)
    return System(buses, generators, loads, branches, storage, basepower, forecasts, services, annex; kwargs...)
end

"""Constructs a non-functional System for demo purposes."""
function System(::Nothing; buses=[Bus(nothing)],
                generators=[ThermalStandard(nothing), RenewableFix(nothing)],
                loads=[PowerLoad(nothing)],
                branches=nothing,
                storage=nothing,
                basepower=100.0,
                forecasts = nothing,
                services = nothing,
                annex = nothing,
                kwargs...)
    return System(buses, generators, loads, branches, storage, basepower, forecasts, services, annex; kwargs...)
end

"""
    to_json(sys::System, filename::AbstractString)

Serializes a system to a JSON string.
"""
function to_json(sys::System, filename::AbstractString)
    return IS.to_json(sys, filename)
end

"""
    to_json(io::IO, sys::System)

Serializes a system an IO stream in JSON.
"""
function to_json(io::IO, sys::System)
    return IS.to_json(io, sys)
end

"""Constructs a System from a JSON file."""
function System(filename::String)
    sys = IS.from_json(System, filename)
    check!(sys)
    return sys
end

"""
    add_component!(sys::System, component::T; kwargs...) where T <: Component

Add a component to the system.

Throws ArgumentError if the component's name is already stored for its concrete type.

Throws InvalidRange if any of the component's field values are outside of defined valid
range.

# Examples
```julia
sys = System(100.0)

# Add a single component.
add_component!(sys, bus)

# Add many at once.
buses = [bus1, bus2, bus3]
generators = [gen1, gen2, gen3]
foreach(x -> add_component!(sys, x), Iterators.flatten((buses, generators)))
```
"""
function add_component!(sys::System, component::T; kwargs...) where T <: Component
    if T <: Branch
        arc = get_arc(component)
        check_bus(sys, get_from(arc), arc)
        check_bus(sys, get_to(arc), arc)
    elseif Bus in fieldtypes(T)
        check_bus(sys, get_bus(component), component)
    end

    if sys.runchecks && !validate_struct(sys, component)
        throw(InvalidValue("Invalid value for $(component)"))
    end

    IS.add_component!(sys.data, component; kwargs...)
end

"""
    add_forecasts!(
                   sys::System,
                   metadata::Union{AbstractString, Vector{IS.TimeseriesFileMetadata}};
                   resolution=nothing,
                  )

Adds forecasts from a metadata file or metadata descriptors.

# Arguments
- `sys::System`: system
- `metadata::Union{AbstractString, Vector{IS.TimeseriesFileMetadata}}`: metdata filename
  that includes an array of IS.TimeseriesFileMetadata instances or a vector.
- `resolution::DateTime.Period=nothing`: skip forecast that don't match this resolution.
"""
function add_forecasts!(
                        sys::System,
                        metadata::Union{AbstractString, Vector{IS.TimeseriesFileMetadata}};
                        resolution=nothing)
    forecasts = Vector{Forecast}()
    for forecast in make_forecasts(sys, metadata; resolution=resolution)
        component = IS.get_component(forecast)
        if component isa LoadZones
            uuids = Set([IS.get_uuid(x) for x in get_buses(component)])
            for component_ in (load for load in get_components(ElectricLoad, sys)
                              if get_bus(load) |> IS.get_uuid in uuids)
                if forecast isa Deterministic
                    forecast_ = Deterministic(component_,
                                                 IS.get_label(forecast),
                                                 IS.get_timeseries(forecast))
                # TODO: others
                else
                    @assert false
                end
                push!(forecasts, forecast_)
            end
        else
            push!(forecasts, forecast)
        end
    end

    add_forecasts!(sys, forecasts)
end

"""
    iterate_components(sys::System)

Iterates over all components.

# Examples
```julia
for component in iterate_components(sys)
    @show component
end
```

See also: [`get_components`](@ref)
"""
function iterate_components(sys::System)
    return IS.iterate_components(sys.data)
end

"""
    remove_components!(::Type{T}, sys::System) where T <: Component

Remove all components of type T from the system.

Throws ArgumentError if the type is not stored.
"""
function remove_components!(::Type{T}, sys::System) where T <: Component
    return IS.remove_components!(T, sys.data)
end

"""
    remove_component!(sys::System, component::T) where T <: Component

Remove a component from the system by its value.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(sys::System, component)
    return IS.remove_component!(sys.data, component)
end

"""
    remove_component!(
                      ::Type{T},
                      sys::System,
                      name::AbstractString,
                      ) where T <: Component

Remove a component from the system by its name.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(
                           ::Type{T},
                           sys::System,
                           name::AbstractString,
                          ) where T <: Component
    return IS.remove_component!(T, sys.data, name)
end

"""
    get_component(
                  ::Type{T},
                  sys::System,
                  name::AbstractString
                 )::Union{T, Nothing} where {T <: Component}

Get the component of concrete type T with name. Returns nothing if no component matches.

See [`get_components_by_name`](@ref) if the concrete type is unknown.

Throws ArgumentError if T is not a concrete type.
"""
function get_component(::Type{T}, sys::System, name::AbstractString) where T <: Component
    return IS.get_component(T, sys.data, name)
end

"""
    get_components(
                   ::Type{T},
                   sys::System,
                  )::FlattenIteratorWrapper{T} where {T <: Component}

Returns an iterator of components. T can be concrete or abstract.
Call collect on the result if an array is desired.

# Examples
```julia
iter = PowerSystems.get_components(ThermalStandard, sys)
iter = PowerSystems.get_components(Generator, sys)
generators = PowerSystems.get_components(Generator, sys) |> collect
generators = collect(PowerSystems.get_components(Generator, sys))
```

See also: [`iterate_components`](@ref)
"""
function get_components(
                        ::Type{T},
                        sys::System,
                       )::FlattenIteratorWrapper{T} where {T <: Component}
    return IS.get_components(T, sys.data)
end

"""
    get_components_by_name(
                           ::Type{T},
                           sys::System,
                           name::AbstractString
                          )::Vector{T} where {T <: Component}

Get the components of abstract type T with name. Note that PowerSystems enforces unique
names on each concrete type but not across concrete types.

See [`get_component`](@ref) if the concrete type is known.

Throws ArgumentError if T is not an abstract type.
"""
function get_components_by_name(
                                ::Type{T},
                                sys::System,
                                name::AbstractString
                               )::Vector{T} where {T <: Component}
    return IS.get_components_by_name(T, sys.data, name)
end

"""
    get_component_forecasts(
                            ::Type{T},
                            sys::System,
                            initial_time::Dates.DateTime,
                            ) where T <: Component

Get the forecasts of a component of type T with initial_time.
    The resulting container can contain Forecasts of dissimilar types.

Throws ArgumentError if T is not a concrete type.

See also: [`get_component`](@ref)
"""
function get_component_forecasts(
                                 ::Type{T},
                                 sys::System,
                                 initial_time::Dates.DateTime,
                                ) where T <: Component
    return IS.get_component_forecasts(T, sys.data, initial_time)
end

"""
    add_forecast!(sys::System, forecast)

Adds forecast to the system.

# Arguments
- `sys::System`: system
- `forecast`: Any object of subtype forecast

Throws ArgumentError if the forecast's component is not stored in the system.

"""
function add_forecast!(sys::System, forecast)
    return IS.add_forecast!(sys.data, forecast)
end

"""
    add_forecast!(
                  sys::System,
                  filename::AbstractString,
                  component::Component,
                  label::AbstractString,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

Add a forecast to a system from a CSV file.

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of
scaling_factor.
"""
function add_forecast!(
                       sys::System,
                       filename::AbstractString,
                       component::Component,
                       label::AbstractString,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    return IS.add_forecast!(sys.data, filename, component, label, scaling_factor)
end

"""
    add_forecast!(
                  sys::System,
                  ta::TimeSeries.TimeArray,
                  component,
                  label,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

Add a forecast to a system from a TimeSeries.TimeArray.

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of
scaling_factor.
"""
function add_forecast!(
                       sys::System,
                       ta::TimeSeries.TimeArray,
                       component,
                       label,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    return IS.add_forecast!(sys.data, ta, component, label, scaling_factor)
end

"""
    add_forecast!(
                  sys::System,
                  df::DataFrames.DataFrame,
                  component,
                  label,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

Add a forecast to a system from a DataFrames.DataFrame.

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of
scaling_factor.
"""
function add_forecast!(
                       sys::System,
                       df::DataFrames.DataFrame,
                       component,
                       label,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    return IS.add_forecast!(sys.data, df, component, label, scaling_factor)
end

"""
    add_forecasts!(sys::System, forecasts)

Add forecasts to the system.

# Arguments
- `sys::System`: system
- `forecasts`: iterable (array, iterator, etc.) of Forecast values

Throws DataFormatError if
- A component-label pair is not unique within a forecast array.
- A forecast has a different resolution than others.
- A forecast has a different horizon than others.

Throws ArgumentError if the forecast's component is not stored in the system.

"""
function add_forecasts!(sys::System, forecasts)
    return IS.add_forecasts!(sys.data, forecasts)
end

"""
    make_forecasts(sys::System, metadata_file::AbstractString; resolution=nothing)

Return a vector of forecasts from a metadata file.

# Arguments
- `data::SystemData`: system
- `metadata_file::AbstractString`: path to metadata file
- `resolution::{Nothing, Dates.Period}`: skip any forecasts that don't match this resolution

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of what the file
should contain.
"""
function make_forecasts(sys::System, metadata_file::AbstractString; resolution=nothing)
    return IS.make_forecasts(sys.data, metadata_file, PowerSystems; resolution=resolution)
end

"""
    make_forecasts(data::SystemData, timeseries_metadata::Vector{TimeseriesFileMetadata};
                   resolution=nothing)

Return a vector of forecasts from a vector of TimeseriesFileMetadata values.

# Arguments
- `data::SystemData`: system
- `timeseries_metadata::Vector{TimeseriesFileMetadata}`: metadata values
- `resolution::{Nothing, Dates.Period}`: skip any forecasts that don't match this resolution
"""
function make_forecasts(sys::System, metadata::Vector{IS.TimeseriesFileMetadata};
                        resolution=nothing)
    return IS.make_forecasts(sys.data, metadata, PowerSystems; resolution=resolution)
end

"""
    get_forecasts(::Type{T}, sys::System, initial_time::Dates.DateTime)

Return an iterator of forecasts. T can be concrete or abstract.

Call collect on the result if an array is desired.

This method is fast and efficient because it returns an iterator to existing vectors.

# Examples
```julia
iter = PowerSystems.get_forecasts(Deterministic{RenewableFix}, sys, initial_time)
iter = PowerSystems.get_forecasts(Forecast, sys, initial_time)
forecasts = PowerSystems.get_forecasts(Forecast, sys, initial_time) |> collect
forecasts = collect(PowerSystems.get_forecasts(Forecast, sys))
```

See also: [`iterate_forecasts`](@ref)
"""
function get_forecasts(
                       ::Type{T},
                       sys::System,
                       initial_time::Dates.DateTime,
                      )::FlattenIteratorWrapper{T} where T <: Forecast
    return IS.get_forecasts(T, sys.data, initial_time)
end

"""
    get_forecasts(
                  ::Type{T},
                  sys::System,
                  initial_time::Dates.DateTime,
                  components_iterator,
                  label::Union{String, Nothing}=nothing,
                 )::Vector{Forecast}

# Arguments
- `forecasts::Forecasts`: system
- `initial_time::Dates.DateTime`: time designator for the forecast
- `components_iter`: iterable (array, iterator, etc.) of Component values
- `label::Union{String, Nothing}`: forecast label or nothing

Return forecasts that match the components and label.

This method is slower than the first version because it has to compare components and label
as well as build a new vector.

Throws ArgumentError if eltype(components_iterator) is a concrete type and no forecast is
found for a component.
"""
function get_forecasts(
                       ::Type{T},
                       sys::System,
                       initial_time::Dates.DateTime,
                       components_iterator,
                       label::Union{String, Nothing}=nothing,
                      )::Vector{T} where T <: Forecast
    return IS.get_forecasts(T, sys.data, initial_time, components_iterator, label)
end

"""
    iterate_forecasts(sys::System)

Iterates over all forecasts in order of initial time.

# Examples
```julia
for forecast in iterate_forecasts(sys)
    @show forecast
end
```

See also: [`get_forecasts`](@ref)
"""
function iterate_forecasts(sys::System)
    return IS.iterate_forecasts(sys.data)
end

"""
    remove_forecast(sys::System, forecast::Forecast)

Remove the forecast from the system.

Throws ArgumentError if the forecast is not stored.
"""
function remove_forecast!(sys::System, forecast::Forecast)
    return IS.remove_forecast!(sys.data, forecast)
end

"""
    clear_forecasts!(sys::System)

Remove all forecasts from the system.
"""
function clear_forecasts!(sys::System)
    return IS.clear_forecasts!(sys.data)
end

"""
    split_forecasts!(
                     sys::System,
                     forecasts::FlattenIteratorWrapper{T}, # must be an iterable
                     interval::Dates.Period,
                     horizon::Int,
                    ) where T <: Forecast

Replaces system forecasts with a set of forecasts by incrementing through an iterable
set of forecasts by interval and horizon.

"""
function split_forecasts!(
                          sys::System,
                          forecasts::FlattenIteratorWrapper{T}, # must be an iterable
                          interval::Dates.Period,
                          horizon::Int,
                         ) where T <: Forecast
    return IS.split_forecasts!(sys.data, forecasts, interval, horizon)
end

"""
    get_forecast_initial_times(sys::System)::Vector{Dates.DateTime}

Return sorted forecast initial times.

"""
function get_forecast_initial_times(sys::System)::Vector{Dates.DateTime}
    return IS.get_forecast_initial_times(sys.data)
end

"""
    get_forecasts_horizon(sys::System)

Return the horizon for all forecasts.
"""
function get_forecasts_horizon(sys::System)
    return IS.get_forecasts_horizon(sys.data)
end

"""
    get_forecasts_initial_time(sys::System)

Return the earliest initial_time for a forecast.
"""
function get_forecasts_initial_time(sys::System)
    return IS.get_forecasts_initial_time(sys.data)
end

"""
    get_forecasts_interval(sys::System)

Return the interval for all forecasts.
"""
function get_forecasts_interval(sys::System)
    return IS.get_forecasts_interval(sys.data)
end

"""
    get_forecasts_resolution(sys::System)

Return the resolution for all forecasts.
"""
function get_forecasts_resolution(sys::System)
    return IS.get_forecasts_resolution(sys.data)
end

"""
    validate_struct(sys::System, value::PowerSystemType)

Validates an instance of a PowerSystemType against System data.
Returns true if the instance is valid.

Users implementing this function for custom types should consider implementing
InfrastructureSystems.validate_struct instead if the validation logic only requires data
contained within the instance.
"""
function validate_struct(sys::System, value::PowerSystemType)::Bool
    return true
end

function check!(sys::System)
    buses = get_components(Bus, sys)
    slack_bus_check(buses)
    buscheck(buses)
end

function JSON2.read(io::IO, ::Type{System})
    raw = JSON2.read(io, NamedTuple)
    sys = System(float(raw.basepower); runchecks=raw.runchecks)
    component_cache = Dict{Base.UUID, Component}()

    # Buses and Arcs are encoded as UUIDs.
    composite_components = [Bus]
    for composite_component in composite_components
        for component in IS.get_components_raw(IS.SystemData, composite_component, raw.data)
            comp = IS.convert_type(composite_component, component)
            add_component!(sys, comp)
            component_cache[IS.get_uuid(comp)] = comp
        end
    end

    # Skip Services this round because they have Devices.
    for c_type_sym in IS.get_component_types_raw(IS.SystemData, raw.data)
        c_type = getfield(PowerSystems, Symbol(IS.strip_module_names(string(c_type_sym))))
        (c_type in composite_components || c_type <: Service) && continue
        for component in IS.get_components_raw(IS.SystemData, c_type, raw.data)
            comp = IS.convert_type(c_type, component, component_cache)
            add_component!(sys, comp)
            component_cache[IS.get_uuid(comp)] = comp
        end
    end

    # Now get the Services.
    for c_type_sym in IS.get_component_types_raw(IS.SystemData, raw.data)
        c_type = getfield(PowerSystems, Symbol(IS.strip_module_names(string(c_type_sym))))
        if c_type <: Service
            for component in IS.get_components_raw(IS.SystemData, c_type, raw.data)
                comp = IS.convert_type(c_type, component, component_cache)
                add_component!(sys, comp)
                component_cache[IS.get_uuid(comp)] = comp
            end
        end
    end

    # Now get the Forecasts, which have reference to components.
    IS.convert_forecasts!(sys.data, raw.data, component_cache)
    return sys
end

function JSON2.write(io::IO, component::T) where T <: Component
    return JSON2.write(io, encode_for_json(component))
end

function JSON2.write(component::T) where T <: Component
    return JSON2.write(encode_for_json(component))
end

"""
Encode composite buses as UUIDs.
"""
function encode_for_json(component::T) where T <: Component
    fields = fieldnames(T)
    vals = []

    for name in fields
        val = getfield(component, name)
        if val isa Bus
            push!(vals, IS.get_uuid(val))
        else
            push!(vals, val)
        end
    end

    return NamedTuple{fields}(vals)
end

function IS.convert_type(
                         ::Type{T},
                         data::NamedTuple,
                         component_cache::Dict,
                        ) where T <: Component
    @debug T data
    values = []
    for (fieldname, fieldtype)  in zip(fieldnames(T), fieldtypes(T))
        val = getfield(data, fieldname)
        if fieldtype <: Bus
            uuid = Base.UUID(val.value)
            bus = component_cache[uuid]
            push!(values, bus)
        elseif fieldtype <: Component
            # Recurse.
            push!(values, IS.convert_type(fieldtype, val, component_cache))
        else
            obj = IS.convert_type(fieldtype, val)
            push!(values, obj)
        end
    end

    return T(values...)
end

function get_bus(sys::System, bus_number::Int)
    for bus in get_components(Bus, sys)
        if bus.number == bus_number
            return bus
        end
    end

    return nothing
end

"""
    get_buses(sys::System, bus_numbers::Set{Int})

Return all buses values with bus_numbers.
"""
function get_buses(sys::System, bus_numbers::Set{Int})
    buses = Vector{Bus}()
    for bus in get_components(Bus, sys)
        if bus.number in bus_numbers
            push!(buses, bus)
        end
    end

    return buses
end

"""
Throws ArgumentError if the bus is not stored in the system.
"""
function check_bus(sys::System, bus::Bus, component::Component)
    name = get_name(bus)
    if isnothing(get_component(Bus, sys, name))
        throw(ArgumentError("$component has bus $name that is not stored in the system"))
    end
end

function IS.compare_values(x::System, y::System)::Bool
    match = true

    if !IS.compare_values(x.data, y.data)
        @debug "SystemData values do not match"
        match = false
    end

    if x.basepower != y.basepower
        @debug "basepower does not match" x.basepower y.basepower
        match = false
    end

    return match
end

function _create_system_data_from_kwargs(; kwargs...)
    validation_descriptor_file = nothing
    runchecks = get(kwargs, :runchecks, true)
    if runchecks
        validation_descriptor_file = get(kwargs, :configpath,
                                         POWER_SYSTEM_STRUCT_DESCRIPTOR_FILE)
    end

    return IS.SystemData(; validation_descriptor_file=validation_descriptor_file)
end


text:
@with_kw mutable struct AdagradParameters{T <: AbstractFloat} <: AbstractPolicyParameters
    ϵ::T = 1e-6
end

struct Adagrad <: AbstractSmoothing
    params::AdagradParameters{Float64}
    grms::Vector{Float64}

    function (::Type{Adagrad})(; kwargs...)
        return new(AdagradParameters(; kwargs...),Vector{Float64}())
    end
end

function initialize!(adagrad::Adagrad,x₀::Vector{Float64})
    resize!(adagrad.grms,length(x₀))
    adagrad.grms .= zeros(length(x₀))
end

function smooth!(adagrad::Adagrad,klocal::Integer,kglobal::Integer,x::AbstractVector,gprev::AbstractVector,gcurr::AbstractVector)
    @unpack ϵ = adagrad.params
    adagrad.grms .+= gprev .* gprev
    gcurr .= gprev ./ (.√adagrad.grms .+ ϵ)
end


text:
@testset "LinkedList" begin

    @testset "basic tests" begin
        l = list(2, 3)
        @test length(l) == 2
        @test eltype(l) == Int
        @test eltype(typeof(l)) == Int
    end

    @testset "l0" begin
        l0 = nil(Char)
        @test length(l0) == 0
        @test l0 == nil(Char)
        @test l0 == nil()
        @test sprint(show,l0) == "nil(Char)"
    end

    @testset "l1" begin
        l1 = nil()
        @test length(l1) == 0
        @test l1 == nil()
        @test l1 == nil(Int)
        @test sprint(show,l1) == "nil()"
        @test typeof(list()) === typeof(l1)
        @test copy(l1) == l1
        @test map((x) -> x*2,l1) == l1
    end

    @testset "l2" begin
        l1 = nil()
        l2 = cons(1, l1)
        @test length(l2) == 1
        @test head(l2) == 1
        @test l2 == cons(1, l1)
        @test l2 == list(1)
        @test sprint(show,l2) == "list(1)"
        @test cat(l2) == l2
    end

    @testset "l3" begin
        l3 = list(2, 3)
        @test isa(l3, Cons{Int})
        @test length(l3) == 2
        @test head(l3) == 2
        @test head(tail(l3)) == 3
        @test l3 == list(2,3)
        @test collect(l3) == [2; 3]
        @test collect(copy(l3)) == [2; 3]
        @test sprint(show,l3) == "list(2, 3)"
    end

    @testset "l4" begin
        l1 = nil()
        l2 = cons(1, l1)
        l3 = list(2, 3)
        l4 = cat(l1, l2, l3)
        @test length(l4) == 3
        @test l4 == list(1, 2, 3)
        @test collect(l4) == [1; 2; 3]
        @test collect(copy(l4)) == [1; 2; 3]
        @test sprint(show,l4) == "list(1, 2, 3)"
    end

    @testset "l5" begin
        l4 = list(1, 2, 3)
        l5 = map((x) -> x*2, l4)
        @test isa(l5, Cons{Int})
        @test collect(l5) == [2; 4; 6]
    end

    @testset "l5b" begin
        l5 = list(2, 4, 6)
        l5b = map((x) -> "$x", l5)
        @test isa(l5b, Cons{String})
        @test collect(l5b) == ["2"; "4"; "6"]
    end

    @testset "l5" begin
        l5 = list(2, 4, 6)
        l6 = filter((x) -> x < 6, l5)
        @test length(l6) == 2
        @test l6.head == 2
        @test l6.tail.head == 4
        @test collect(l6) == [2, 4]
    end

    @testset "l7" begin
        l6 = list(2, 4)
        l7 = reverse(l6)
        @test length(l7) == 2
        @test l7.head == 4
        @test l7.tail.head == 2
    end

    @testset "l8" begin
        l5b = list("2", "4", "6")
        l6 = list(2, 4)
        l8 = cat(l5b, l6)
        @test collect(l8) == ["2"; "4"; "6"; 2; 4]
    end

    @testset "l9" begin
        l9 = cat(list(1, 2), list(3.0, 4.0))
        @test isa(l9, Cons{Real})
        @test collect(l9) == [1; 2; 3.0; 4.0]
    end

    @testset "l10" begin
        l10 = list(2, 4, 5.6, 10.5)
        @test collect(l10) == [2; 4; 5.6; 10.5]
    end

    @testset "test identity map" begin
        ex = :(a+b+2 * 2 + foo(2))
        l11 = list(ex.args...)
        @test collect(map(x->x, l11)) == collect(l11)
    end

end # @testset LinkedList


text:
# Copyright 2019 Martin Holters
# See accompanying license file.

function wrap_nleq_expr(nn, nq, expr)
    Base.depwarn("nonlinear_eq should be given as a function, not an expression", :Element)

    res_symbs = [gensym("res_$n") for n in 1:nn]
    J_symbs = [gensym("J_$(m)_$n") for m in 1:nn, n in 1:nq]

    function rewrite_refs(expr::Expr)
        if expr.head == :ref
            if expr.args[1] === :res
                return res_symbs[expr.args[2]]
            elseif expr.args[1] === :J
                return J_symbs[expr.args[2], expr.args[3]]
            end
            return expr
        else
            return Expr(expr.head, rewrite_refs.(expr.args)...)
        end
    end

    rewrite_refs(x::Any) = x

    expr = rewrite_refs(expr)

    return eval(quote
        @inline function (q)
            $(nn > 0 || nq > 0 ? Expr(:local, res_symbs..., J_symbs...) : nothing)
            $(expr)
            res = SVector{$nn}($(res_symbs...))
            J = SMatrix{$nn,$nq}($(J_symbs...))
            return (res, J)
        end
    end)
end

function ports_from_old_pins(pins)
    ports = [pins[2i-1] => pins[2i] for i in 1:length(pins)÷2]
    Base.depwarn("`Element(..., pins=$(repr(pins)))` is depreated, use `Element(..., ports=$(repr(ports, context=:typeinfo=>typeof(ports))))` instead", :Element)
    return ports
end

Base.@deprecate(composite_element(circ::Circuit, pins::Vector{<:Pair}),
    composite_element(circ,  pinmap=Dict(pins...)))


text:
using SCS
const SCS_PRIMAL_FACTORY =
    MOI.OptimizerWithAttributes(SCS.Optimizer, MOI.Silent() => true)
const SCS_DUAL_FACTORY = dual_optimizer(SCS_PRIMAL_FACTORY)
const SCS_PRIMAL_OPT = MOI.instantiate(SCS_PRIMAL_FACTORY)
const SCS_DUAL_OPT = MOI.instantiate(SCS_DUAL_FACTORY)

push!(primal_power_cone_factory, SCS_PRIMAL_FACTORY)
push!(dual_power_cone_factory, SCS_DUAL_FACTORY)
push!(dual_power_cone_optimizer, SCS_DUAL_OPT)
push!(primal_power_cone_optimizer, SCS_PRIMAL_OPT)

@testset "SCS Exponential Cone problems" begin
    list_of_exp_problems = [exp1_test, exp2_test]
    test_strong_duality(list_of_exp_problems, SCS_PRIMAL_FACTORY)
end

@testset "SCS SOC problems" begin
    list_of_soc_problems = [soc3_test, soc4_test, soc5_test, soc6_test]
    test_strong_duality(list_of_soc_problems, SCS_PRIMAL_FACTORY)
    list_of_soc_problems = [soc1_test, soc2_test]
    test_strong_duality(
        list_of_soc_problems,
        SCS_PRIMAL_FACTORY,
        atol = 1e-3,
        rtol = 1e-3,
    )
end

@testset "SCS RotatedSOC problems" begin
    list_of_rsoc_problems = [rsoc1_test, rsoc2_test, rsoc3_test]
    test_strong_duality(list_of_rsoc_problems, SCS_PRIMAL_FACTORY)
    list_of_rsoc_problems = [rsoc4_test]
    test_strong_duality(
        list_of_rsoc_problems,
        SCS_PRIMAL_FACTORY,
        atol = 1e-3,
        rtol = 1e-3,
    )
end

@testset "SCS SDP triangle problems" begin
    list_of_sdp_triang_problems =
        [sdpt1_test, sdpt2_test, sdpt3_test, sdpt4_test]
    test_strong_duality(list_of_sdp_triang_problems, SCS_PRIMAL_FACTORY)
end

@testset "SCS Power Cone problems" begin
    list_of_pow_problems = [pow1_test, pow2_test]
    test_strong_duality(list_of_pow_problems, SCS_PRIMAL_FACTORY)
end


text:
# Tests for Cost DEA Models
@testset "CostDEAModel" begin

    ## Test Cost DEA Model with Cooper et al. (2007)
    # Test agains book results
    X = [3 2; 1 3; 4 6]
    Y = [3; 5; 6]
    W = [4 2; 4 2; 4 2]

    deacostcooper = deacost(X, Y, W, rts = :CRS)

    @test typeof(deacostcooper) == CostDEAModel

    @test nobs(deacostcooper) == 3
    @test ninputs(deacostcooper) == 2
    @test noutputs(deacostcooper) == 1
    @test ismonetary(deacostcooper) == false

    @test efficiency(deacostcooper, :Economic)   ≈ [0.375; 1; 0.429] atol = 1e-3
    @test efficiency(deacostcooper, :Technical)  ≈ [0.9  ; 1; 0.6  ] atol = 1e-3
    @test efficiency(deacostcooper, :Allocative) ≈ [0.417; 1; 0.714] atol = 1e-3

    ## Test Cost DEA Model with Zofío and Prieto (2006) data.
    # Test agains results in R
    X = [5 3; 2 4; 4 2; 4 8; 7 9]
    Y = [7 4; 10 8; 8 10; 5 4; 3 6]
    W = [2 1; 2 1; 2 1; 2 1; 2 1.0]

    # Cost CRS
    deacostcrs = deacost(X, Y, W, rts = :CRS)
    @test efficiency(deacostcrs) ≈ [0.4307692308;
                                    1.000;
                                    1.000;
                                    0.250;
                                    0.2608695652]
    @test efficiency(deacostcrs, :Technical) ≈ [0.6364;
                                1.000;
                                1.000;
                                0.250;
                                0.2609]  atol = 1e-3
    @test efficiency(deacostcrs, :Allocative) ≈ [0.6769230769;
                                1.000;
                                1.000;
                                1.000;
                                1.000]

    @test efficiency(deacost(targets(deacostcrs, :X), targets(deacostcrs, :Y), W, rts = :CRS)) ≈ ones(5)

    # Cost VRS
    deacostvrs = deacost(X, Y, W, rts = :VRS)
    @test efficiency(deacostvrs) ≈ [0.6153846154;
                                    1.000;
                                    1.000;
                                    0.500;
                                    0.3478260870]
    @test efficiency(deacostvrs, :Technical) ≈ [0.750;
                                1.000;
                                1.000;
                                0.500 ;
                                0.375]  atol = 1e-3
    @test efficiency(deacostvrs, :Allocative) ≈ [0.8205128205;
                                1.000;
                                1.000;
                                1.000;
                                0.9275362319]

    @test efficiency(deacost(targets(deacostvrs, :X), targets(deacostvrs, :Y), W, rts = :VRS)) ≈ ones(5)

    # Check defaults
    @test efficiency(deacost(X, Y, W)) == efficiency(deacostvrs)
    @test efficiency(deacostvrs, :Economic) == efficiency(deacostvrs)

    # Print
    show(IOBuffer(), deacostcooper)

    # Test errors
    @test_throws DimensionMismatch deacost([1; 2 ; 3], [4 ; 5], [1; 1; 1]) #  Different number of observations
    @test_throws DimensionMismatch deacost([1; 2; 3], [4; 5; 6], [1; 2; 3; 4]) # Different number of observation in prices
    @test_throws DimensionMismatch deacost([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3]) # Different number of input prices and inputs
    @test_throws ArgumentError deacost([1; 2; 3], [4; 5; 6], [1; 2; 3], rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError deacost([1; 2; 3], [4; 5; 6], [1; 2; 3], dispos = :Error) # Invalid disposability
    @test_throws ArgumentError normfactor(deacost(X, Y, W)) # CostDEAModel has no normalization factor

    # ------------------
    # Weak Disposability Tests
    # ------------------

    X = [1; 2; 3; 2; 4]
    Y = [2; 3; 4; 1; 3]
    W = [1; 1; 1; 1; 1]

    deacostStrong = deacost(X, Y, W, dispos = :Strong)
    @test efficiency(deacostStrong, :Economic) ≈ [1.0; 1.0; 1.0; 0.5; 0.5]
    @test efficiency(deacostStrong, :Technical) ≈ [1.0; 1.0; 1.0; 0.5; 0.5]
    @test efficiency(deacostStrong, :Allocative) ≈ [1.0; 1.0; 1.0; 1.0; 1.0]

    deacostWeak = deacost(X, Y, W, dispos = :Weak)
    @test efficiency(deacostWeak, :Economic) ≈ [1.0; 1.0; 1.0; 1.0; 0.5]
    @test efficiency(deacostWeak, :Technical) ≈ [1.0; 1.0; 1.0; 1.0; 0.5]
    @test efficiency(deacostWeak, :Allocative) ≈ [1.0; 1.0; 1.0; 1.0; 1.0]

    # ------------------
    # Test Vector and Matrix inputs and outputs
    # ------------------
    # Tests against results in R

    # Inputs is Matrix, Outputs is Vector
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1]

    @test efficiency(deacost(X, Y, W)) ≈ [1; 0.8; 0.8; 0.5714285714; 0.4; 0.5714285714; 0.5714285714; 0.4166666667]

    # Inputs is Vector, Output is Matrix
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]
    W = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deacost(X, Y, W)) ≈ [1; 1; 1; 1; 1; 1; 1; 1]

    # Inputs is Vector, Output is Vector
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]
    W = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deacost(X, Y, W)) ≈ [1; 1; 1; 1; 0.5; 0.4761904762; 0.8571428571; 0.2843710157]

end


text:
using Test
using AutoDiff.Backward

@testset "f(x,y,z) = x + 2*y/z" begin
    f(x, y, z) = x + 2.0*y/z
    g = gradient(f)

    @test all(isapprox.(g(1.0, 2.0, 3.0), [1.0, 2.0/3.0, -4.0/9.0]))
end

@testset "one-arg functions" begin
    f(x) = 2.0*x + 3.0
    g = gradient(f)

    @test isapprox(g(1.0), 2.0)
end

@testset "vector-mode functions" begin
    f(xs) = sqrt(xs[1]*xs[1] + xs[2]*xs[2] + xs[3]*xs[3])
    g = gradient(f)
    xs = randn(3)
    @test all(isapprox.(g(xs), [xs[1]/f(xs), xs[2]/f(xs), xs[3]/f(xs)]))
end

@testset "array-mode functions" begin
    f(xs) = sqrt(sum(xs.*xs))
    g = gradient(f)
    xs = randn(3, 5, 7)
    r = f(xs)
    @test all(isapprox.(g(xs), xs./r))
    @test size(g(xs))==(3,5,7)
end

@testset "reuse variable" begin
    f(x, y, z) = x - 3.0*x/(y*z) + z*y
    g = gradient(f)

    x = randn()
    y = randn()
    z = randn()

    @test all(isapprox.(g(x, y, z), [1.0 - 3.0/(y*z), 3.0*x/(y*y*z) + z, 3.0*x/(y*z*z) + y]))
end

@testset "unary minus" begin
    f(x,y,z) = -x + (-y)*z
    g = gradient(f)
    x,y,z = randn(3)
    @test all(isapprox.(g(x,y,z), [-1.0, -z, -y]))
end

@testset "exp" begin
    f(x,y,z) = x*exp(y*z)
    g = gradient(f)
    x,y,z = randn(3)
    ff = f(x,y,z)
    @test all(isapprox.(g(x,y,z), [ff/x, z*ff, y*ff]))
end

@testset "sqrt" begin
    f(x,y,z) = sqrt(x*x + y*y + z*z)
    f_and_g = value_and_gradient(f)
    x,y,z = randn(3)

    ff, gg = f_and_g(x,y,z)

    @test all(isapprox.(gg, [x/ff, y/ff, z/ff]))
end

@testset "trig" begin
    @testset "sin" begin
        f(x,y,z) = x*sin(y*z)
        g = gradient(f)
        x,y,z = randn(3)
        @test all(isapprox.(g(x,y,z), [sin(y*z), x*z*cos(y*z), x*y*cos(y*z)]))
    end

    @testset "cos" begin
        f(x,y,z) = x*cos(y*z)
        g = gradient(f)
        x,y,z = randn(3)
        @test all(isapprox.(g(x,y,z), [cos(y*z), -z*x*sin(y*z), -x*y*sin(y*z)]))
    end

    @testset "tan" begin
        f(x,y,z) = x*tan(y*z)
        ff(x,y,z) = x*sin(y*z)/cos(y*z)
        g = gradient(f)
        gg = gradient(ff)
        x,y,z = randn(3)
        @test all(isapprox.(g(x,y,z), gg(x,y,z)))
    end

    @testset "sec" begin
        f(x,y,z) = x*sec(y*z)
        ff(x,y,z) = x/cos(y*z)
        g = gradient(f)
        gg = gradient(ff)
        x,y,z = randn(3)
        @test all(isapprox.(g(x,y,z), gg(x,y,z)))
    end

    @testset "csc" begin
        f(x,y,z) = x*csc(y*z)
        ff(x,y,z) = x/sin(y*z)
        g = gradient(f)
        gg = gradient(ff)
        x,y,z = randn(3)
        @test all(isapprox.(g(x,y,z), gg(x,y,z)))
    end

    @testset "cot" begin
        f(x,y,z) = x*cot(y*z)
        ff(x,y,z) = x/tan(y*z)
        g = gradient(f)
        gg = gradient(ff)
        x,y,z = randn(3)
        @test all(isapprox.(g(x,y,z), gg(x,y,z)))
    end

    @testset "atan" begin
        f(x,y,z) = atan(fff(x,y,z))
        ff(x,y,z) = tan(f(x,y,z))
        fff(x,y,z) = sqrt(x*x+y*y+z*z)
        g = gradient(f)
        gg = gradient(ff)
        ggg = gradient(fff)
        x,y,z = randn(3)
        @test all(isapprox.(ggg(x,y,z), gg(x,y,z)))
    end

    @testset "atan2" begin
        f(theta) = atan(sin(theta),cos(theta))
        g = gradient(f)
        theta = 2*pi*randn()
        @test isapprox(g(theta), 1.0)
    end
end

@testset "sum like functions" begin
    @testset "sum" begin
        f(x, y, z) = sum([x,y,z])
        g = gradient(f)
        x,y,z = randn(3)
        @test all(isapprox.(g(x,y,z), ones(3)))
    end
end

@testset "log" begin
    f(x,y,z) = log(x*y*z)
    g = gradient(f)
    x,y,z = abs.(randn(3))
    @test all(isapprox.(g(x,y,z), [1.0/x, 1.0/y, 1.0/z]))
end

@testset "log1p" begin
    f(x,y,z) = log1p(x*y*z)
    g = gradient(f)
    x,y,z = 1e-3.*randn(3)
    @test all(isapprox.(g(x,y,z), [y*z/(1+x*y*z), x*z/(1+x*y*z), x*y/(1+x*y*z)]))
end

@testset "logsumexp" begin
    f(x,y,z) = logsumexp(x, logsumexp(y, z))
    f_exact(x,y,z) = log(exp(x) + exp(y) + exp(z))
    g = gradient(f)
    g_exact = gradient(f_exact)
    gg = gradient(logsumexp) # Array version

    x,y,z = randn(3)

    @test all(isapprox.(g(x,y,z), g_exact(x,y,z)))
    @test all(isapprox.(gg([x,y,z]), g_exact(x,y,z)))
end

@testset "intermediate reuse" begin
    function f(x,y,z)
        q = x*y + cos(z*x)
        x*sin(q) + z*q*q
    end
    g = gradient(f)
    x,y,z = randn(3)

    @test all(isapprox.(g(x,y,z),
                        [(2.0*z*(x*y + cos(x*z)) + x*cos(x*y + cos(x*z)))*(y-z*sin(z*x)) + sin(x*y + cos(x*z)),
                         x*(2.0*z*(x*y + cos(x*z)) + x*cos(x*y + cos(x*z))),
                         (x*y + cos(x*z))^2 - x*(2.0*z*(x*y + cos(x*z)) + x*(cos(x*y + cos(x*z))))*sin(x*z)]))
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ a341a594-e99b-11eb-376f-af0d7f593c65
begin
    using HTTP, Markdown
end

# ╔═╡ 55bf5b39-c01d-4696-97e2-5273f14a5ef1
begin
	function get_hh( hno :: Int )
		r  = HTTP.request("GET", "http://192.168.1.2:8002/hhld/$hno")
		s = String(r.body)
		@eval @md_str $s
	end

end

# ╔═╡ 75121b89-4e25-4f85-aaef-6de8de0aa18d
get_hh( 4000 )

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"

[compat]
HTTP = "~0.9.12"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "c6a1fff2fd4b1da29d3dccaffb1e1001244d844e"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.12"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ 9ad89f76-bf0f-4626-96b2-c26a1fcc47f0

# ╔═╡ Cell order:
# ╠═a341a594-e99b-11eb-376f-af0d7f593c65
# ╠═55bf5b39-c01d-4696-97e2-5273f14a5ef1
# ╟─75121b89-4e25-4f85-aaef-6de8de0aa18d
# ╠═0946f3fb-afcd-41fe-8496-eb4083894720
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002
# ╠═9ad89f76-bf0f-4626-96b2-c26a1fcc47f0

max_stars_repo_path: scripts/pluto_get_hh backup 1.jl

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_path: scripts/pluto_get_hh backup 1.jl

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_path: scripts/pluto_get_hh backup 1.jl

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function myreverse!(A)
    a = firstindex(A)
    b = lastindex(A)
    count = 0
    while (a < b)
        A[a + count], A[b + count] = A[b + count], A[a + count]
        count += 1
    end
end


text:
# Do `include("abstract.jl")` on two Julia versions and then run this script.
# You will likely have to change the path for `sigstable` below.
# Typically this should be run using the in-development version of Julia (or whatever "latest" is in your comparison)

using PyPlot

function parseline(line)
    m = match(r"(.*) => (.*)", line)
    sigstr, count = m.captures[1], parse(Int, m.captures[2])
    sig = try
        ex = Meta.parse(sigstr)
        eval(ex)
    catch
        return nothing
    end
    return sig, count
end

function parsedata(filename)
    lines = readlines(filename)
    sigcount = IdDict{Any,Int}()
    for line in lines
        ret = parseline(line)
        ret === nothing && continue
        sig, count = ret
        sigcount[sig] = count
    end
    return sigcount
end

function split_comparable(sigc1, sigc2)
    c1, c2, sigs = Int[], Int[], Any[]
    for (sig, c) in sigc1
        push!(sigs, sig)
        push!(c1, sigc1[sig])
        push!(c2, get(sigc2, sig, 0))
    end
    for (sig, c) in sigc2
        if !haskey(sigc1, sig)
            push!(sigs, sig)
            push!(c1, 0)
            push!(c2, c)
        end
    end
    return sigs, c1, c2
end

function tally0(c1, c2)
    nz1, nz2 = 0, 0
    for (a1, a2) in zip(c1, c2)
        a1 == a2 == 0 && continue
        a1 == 0 && (nz1 += 1)
        a2 == 0 && (nz2 += 1)
    end
    return nz1, nz2
end

sigmaster = parsedata("/tmp/methdata_$VERSION.log")
sigstable, stablever = parsedata("/tmp/methdata_1.5.1-pre.29.log"), "1.5"
# sigstable, stablever = parsedata("/tmp/methdata_1.4.2.log"), "1.4"

sigs, c1, c2 = split_comparable(sigstable, sigmaster)
nz1, nz2 = tally0(c1, c2)
println("$stablever has $nz1 with no backedges, master has $nz2")
mx1, mx2 = maximum(c1), maximum(c2)
isexported(sig) = (ft = Base.unwrap_unionall(sig).parameters[1]; isdefined(Main, ft.name.mt.name))
colors = [isexported(sig) ? "magenta" : "green" for sig in sigs]

function on_click(event)
    x, y = event.xdata, event.ydata
    normsqrdist(pr) = ((pr[1]-x)/mx1)^2 + ((pr[2]-y)/mx2)^2
    idx = argmin(normsqrdist.(zip(c1, c2)))
    println(sigs[idx])
end
begin
    hfig, axs = plt.subplots(2, 1)
    plt.subplots_adjust(hspace=0.3)
    logedges = LinRange(0, log10(max(mx1, mx2)+2), 30)
    ax = axs[1]
    ax.hist(log10.(c1 .+ 1), bins=logedges)
    ax.set_xlabel("log₁₀(# backedges + 1), $stablever")
    ax.set_ylabel("# 'at risk' signatures")
    ax = axs[2]
    ax.hist(log10.(c2 .+ 1), bins=logedges)
    ax.set_xlabel("log₁₀(# backedges + 1), 1.6")
    ax.set_ylabel("# 'at risk' signatures")

    display(hfig)
    fig, ax = plt.subplots()
    ax.scatter(c1 .+ 1, c2 .+ 1, c=colors)  # + 1 for the log-scaling
    ax.set_xlabel("# backedges + 1, $stablever")
    ax.set_ylabel("# backedges + 1, 1.6")
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_aspect("equal")
    fig.canvas.callbacks.connect("button_press_event", on_click)
    fig
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    inchi, inchikey

using libinchi_jll

mutable struct inchi_Output
    szInChI::Cstring
    szAuxInfo::Cstring
    szMessage::Cstring
    szLog::Cstring
    inchi_Output() = new(C_NULL, C_NULL, C_NULL, C_NULL)
end

"""
    inchi(molblock::String) -> String
    inchi(mol::GraphMol) -> String

Generate InChI string from molblock string or molecule
"""
function inchi(molblock::String)
    output = inchi_Output()
    ccall(
        (:MakeINCHIFromMolfileText, libinchi),
        Int32,
        (Cstring, Cstring, Ref{inchi_Output}),
        molblock, "-W60", output)
    if output.szInChI == C_NULL
        println("InChI error with $(molblock)")
        if output.szMessage != C_NULL
            println("message: ", unsafe_string(output.szMessage))
        end
        if output.szLog != C_NULL
            println("log: ", unsafe_string(output.szLog))
        end
        res = nothing
    else
        res = unsafe_string(output.szInChI)
    end
    # Free string buffers allocated by MakeINCHIFromMolfileText
    ccall(
        (:FreeINCHI, libinchi),
        Cvoid,
        (Ref{inchi_Output},),
        output
    )
    return res
end

inchi(mol::GraphMol) = inchi(printv2mol(mol))

"""
    inchikey(inchi::String) -> String
    inchikey(mol::GraphMol) -> String

Generate InChI key from InChI string or molecule
"""
function inchikey(inchi::Union{String,Nothing})
    inchi === nothing && return nothing
    ikeybuf = pointer(Vector{UInt8}(undef, 256))
    # TODO: need extra buffer?
    xtra1buf = pointer(Vector{UInt8}(undef, 256))
    xtra2buf = pointer(Vector{UInt8}(undef, 256))
    ccall(
        (:GetINCHIKeyFromINCHI, libinchi),
        Int32,
        (Cstring, Int32, Int32, Cstring, Cstring, Cstring),
        inchi, 1, 1, ikeybuf, xtra1buf, xtra2buf)
    return unsafe_string(ikeybuf)
end

inchikey(mol::GraphMol) = inchikey(inchi(mol))


text:
include("infdecomp_base.jl")
include("expcone.jl")
include("graddesc.jl")
include("infdecomp.jl")
include("read_pdf.jl")
include("dot_try.jl")
using InfDecomp
# solvers_list = [:Mosek, :Ipopt, :Knitro_Ip, :Knitro_IpCG, :Knitro_AS, :Knitro_SQP, :Cvxopt]

#solvers_list = [:Mosek]

tmpFloatDatatype = Float64

# if length(ARGS) != 2
#     print(
#     """Use w/ two arguments: solver & instance.
#     Solver:  one of ECOS_L, Mosek, GD.
#     instance: the name of a file in ../PDFs/Data/
#     """)
#     exit(1)
# end

if length(ARGS) != 3
    print(
    """Use w/ three arguments: solver, instance, stats_file.
    Solver:  one of ECOS_L, Mosek, GD.
    instance: the name of a file in ../PDFs/Data/
    stats_file: the name of the file the statistics should be written in.
    """)
    exit(1)
end

@show ARGS

solver = nothing

if     uppercase(ARGS[1]) == "ECOS_L"
    solver = :ECOS_L
elseif uppercase(ARGS[1]) == "MOSEK"
    solver = :Mosek
elseif uppercase(ARGS[1]) == "SCS_D"
    solver = :SCS_D
elseif uppercase(ARGS[1]) == "GD"
    solver = :My_GradDesc
else
    println("I don't recognize the solver (only ECOS_L, SCS_D, Mosek, GD are allowed)")
    exit(1)
end

instance = ARGS[2]

# test(solver, "AND-0.1-10000.dens", "AND-0.1-10000.dens",  tmpFloatDatatype)

# stats_file = "dot_stats.csv"
stats_file = ARGS[3]

# open(stats_file, "a") do ffile
#     for i in fieldnames(InfDecomp.Solution_and_Stats)
#         print(ffile,"\t",i,", ")
#     end
#     println(ffile,"")
# end

println("########################################################################################################################################################################################################")
println("vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv")
println("Starting computations of the distribution: ", instance," with ", solver)
println("")
test(solver, instance, instance,  tmpFloatDatatype, stats_file=stats_file)
println("")
println("Finishing computations of the distribution: ", instance," with ", solver)
println("^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^")

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
using TrajectoryOptimization
using Test
using LinearAlgebra
using Random
using StaticArrays
using SparseArrays
using ForwardDiff
using RobotDynamics
using BenchmarkTools
const TO = TrajectoryOptimization

include("test_models.jl")

@testset "Costs" begin
    include("cost_tests.jl")
    include("objective_tests.jl")
end

@testset "Constraints" begin
    include("constraint_tests.jl")
    include("dynamics_constraints.jl")
    include("constraint_list.jl")
    include("constraint_sets.jl")
end

@testset "Problems" begin
    include("problems_tests.jl")
end

# @testset "Utils" begin
# end

@testset "NLP" begin
    include("nlp_tests.jl")
    include("moi_test.jl")
end

using NBInclude
@testset "Examples" begin
    # @test_nowarn include(joinpath(@__DIR__, "..", "examples", "quickstart.jl"))
    # @nbinclude(joinpath(@__DIR__, "..", "examples", "Cartpole.ipynb"); softscope=true)
    # @nbinclude(joinpath(@__DIR__, "..", "examples", "Quadrotor.ipynb"); softscope=true)
end


text:
mutable struct Node
    depth
    is_leaf
    prediction
    feature_index
    thresh
    left
    right
end

function Node(depth, is_leaf, prediction)
    Node(depth, is_leaf, prediction, nothing, nothing, nothing, nothing)
end
struct WeightedDecisionStump
    root_node

    function WeightedDecisionStump(X, y, weights)
        node = Node(1, false, nothing)
        feature_index, thresh, left_prediction = split(X, y, weights)
        node.feature_index = feature_index
        node.thresh = thresh
        node.left = Node(2, true, left_prediction)
        node.right = Node(2, true, -left_prediction)

        new(node)
    end
end

function split(X, y, weights)
    best_feature = 1
    best_thresh = minimum(X[:, 1])
    best_left_prediction = 1
    error_val = Inf

    for (feature_index, feature_array) in enumerate(eachcol(X))
        sorted_indices = sortperm(feature_array)
        X_sorted = X[sorted_indices, :]
        y_sorted = y[sorted_indices]
        weights_sorted = weights[sorted_indices]
        predictions = zeros(length(y_sorted))

        for left_prediction = (-1, 1)
            for i = 1:(length(y_sorted) - 1)
                predictions[1:i] .= left_prediction
                predictions[i+1:end] .= -left_prediction
                temp_error_val = misclassification_rate(predictions, y_sorted, weights_sorted)
                if temp_error_val < error_val
                    error_val = temp_error_val
                    best_feature = feature_index
                    best_thresh = (X_sorted[i, feature_index] + X_sorted[i+1, feature_index]) / 2
                    best_left_prediction = left_prediction
                end
            end
        end
    end

    return best_feature, best_thresh, best_left_prediction
end

function misclassification_rate(predictions, y, weights)
    sum((predictions .!= y) .* weights)
end

struct AdaBoostModel
    num_models
    models
    alphas

    function AdaBoostModel(X, y, num_models)
        # AdaBoost model requires labels being -1 or 1
        @assert Set(y) == Set([-1, 1])
        weights = fill(1.0 / length(y), length(y))
        models = WeightedDecisionStump[]
        alphas = Float64[]

        for _ = 1:num_models
            model = WeightedDecisionStump(X, y, weights)
            predictions = predict(model, X)
            mistakes = predictions .!= y
            error = sum(weights .* mistakes) / sum(weights)
            alpha = error == 0.0 ? 1.0 : log((1-error)/error)
            weights = weights .* exp.(alpha .* mistakes)
            weights = weights ./ sum(weights)

            push!(models, model)
            push!(alphas, alpha)
        end

        new(num_models, models, alphas)
    end
end

function predict(model::WeightedDecisionStump, X)
    root_node = model.root_node
    feature_index = root_node.feature_index
    thresh = root_node.thresh
    predictions = zeros(size(X)[1])

    left_indices = X[:, feature_index] .< thresh
    predictions[left_indices] .= root_node.left.prediction
    predictions[.!left_indices] .= root_node.right.prediction
    return predictions
end

function predict(model::AdaBoostModel, X)
    models = model.models
    alphas = model.alphas
    predictions = zeros(size(X)[1])
    votes = zeros(size(X)[1])

    for (alpha, model) in zip(alphas, models)
        predictions = predict(model, X)
        votes .+= alpha.*predictions
    end

    return sign.(sum(votes, dims=2))
end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia