text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")

using JuMP
using Cbc
using Plots
using CSV
using DataFrames
using GLPK

##

# Define the time interval and timesteps
t_s = 1
t_f = 168
offset = 0 #equals basic example offset = 6531 / pv and wind data
timesteps = t_s:t_f

#pv, wind and demand data
#demand = CSV.read("../data/demand_Industriepark.csv", DataFrame)[timesteps, 1]
#wind = CSV.read("../data/wind_Karholz.csv", DataFrame)[timesteps, 1]
#pv = CSV.read("../data/pv_Halle18.csv", DataFrame)[timesteps, 1]
pv = CSV.read("../basic_example_normalized.csv", DataFrame)[timesteps.+offset, 3]
wind = CSV.read("../basic_example_normalized.csv", DataFrame)[timesteps.+offset, 4]
demand = CSV.read("../basic_example_normalized.csv", DataFrame)[timesteps.+offset, 2]

#heatsystem demand assumption
heatdemand = zeros(Int, length(timesteps))
fill!(heatdemand, 0)

#Global system parameter
c_i = .30
c_o = .05
c_w = 2000.
c_pv = 700.
c_storage = 600.
#investment_budget = 1000000.

#adding heat storage parameters
c_heatstorage = 200.
c_heatpump = 10.
COP = 2. #coefficient of perfomance

time_scale = 10. *365*24/length(timesteps)
storage_scale = 1

m = Model(GLPK.Optimizer)

#varables declaration
@variable(m, 0 <= u_pv)
@variable(m, 0 <= u_w)
@variable(m, 0 <= u_storage)
@variable(m, gci[t in timesteps]>=0)
@variable(m, gco[t in timesteps]>=0)
@variable(m, storage[t in timesteps])

#variable declaration heatsystem
@variable(m, 0 <= u_heatpump)
@variable(m, 0 <= u_heatstorage)
@variable(m, heatstorage[t in timesteps])
@variable(m, 0 <= heatpumpflow[t in timesteps])

#energy balance
@constraint(m,[t in timesteps], gci[t]-gco[t]+u_pv*pv[t]+u_w*wind[t]-demand[t]+storage[t] - heatpumpflow[t]/COP ==0) #add - heatpumpflow[t]/COP

#objective function
@objective(m, Min, c_pv/time_scale*u_pv + c_w/time_scale*u_w + u_storage*c_storage/time_scale + c_i*sum(gci) - c_o*sum(gco) + u_heatstorage*c_heatstorage/time_scale + u_heatpump* c_heatpump/time_scale )  #add + u_heatstorage*c_heatstorage/time_scale + u_heatpump* c_heatpump/time_scale

#budget constraint
#@constraint(m, c_pv*u_pv + c_w*u_w + c_storage*u_storage + c_heatstorage*u_heatstorage + c_heatpump*u_heatpump  <= investment_budget )

#battery in the example
@constraint(m, [t in timesteps], -0.5*u_storage*storage_scale <= sum(storage[1:t]))
@constraint(m, [t in timesteps], sum(storage[1:t]) <= 0.5*u_storage*storage_scale)
@constraint(m, sum(storage) == 0)

#heat balance
@constraint(m, [t in timesteps], heatpumpflow[t] + heatstorage[t] - heatdemand[t] == 0) #heatpumpflow can be either positive and negative
@constraint(m, [t in timesteps], heatpumpflow[t] <= u_heatpump) #necessary size of heatpump

#heatstorage constraints like in battary storage
@constraint(m, [t in timesteps], 0 <= sum(heatstorage[1:t])+ 0.5 * u_heatstorage)
@constraint(m, [t in timesteps], sum(heatstorage[1:t]) <= 0.5 * u_heatstorage)
@constraint(m, [t in timesteps], sum(heatstorage) == 0)

print(m)

status = optimize!(m)

println("Objective value: ", getobjectivevalue(m))
println("u_pv = ", getvalue(u_pv))
println("u_w = ", getvalue(u_w))
println("u_storage = ", getvalue(u_storage))
println("u_heatstorage = ", getvalue(u_heatstorage))
println("u_heatpump = ", getvalue(u_heatpump))
println("gci = ", value.(gci))
println("gco = ", value.(gco))

## Additional test code

##midday without ancillary service

#data supply
#demand = [1164.728, 1182.104, 1194.048]
#pv = [0.0391025641025641, 0.019230769230769232, 0.023717948717948717]
#wind = [0.9772296015180265, 0.9276462338090916, 0.9802821549377114]

##Utility functions
#=function clamped_sin(x)
    return sin(x-8)>0 ? 3. .*sin(x-8) : 0.
end=#
#pv = [clamped_sin(x/24. * 2. * pi) for x in timesteps]

#timesteps = 1:24
#pv = [clamped_sin(x/24. * 2. * pi) for x in timesteps]

##
#defining battary constraint
#@constraint(m, [t in timesteps], -0.5*u_b*b_scale <= sum(b[1:t]))
#@constraint(m, [t in timesteps], sum(b[1:t]) <= 0.5*u_b*b_scale)
#@constraint(m, b[1]==b[length(timesteps)])

#defining heatpump constraint, with binary vaiable z
#@constraint(m, [t in timesteps], heatdemand[t] <= 1000 * (1 - z[t]) )
#@constraint(m, [t in timesteps], -heatpumpflow[t] * 1.1 + heatdemand[t]  == 1000*z[t])

#Energybalance without heat storage:  -heatdemand[t]* Coefficient of Performance (COP) from heatpump
#Energybalance with heat storage: -(heatdemand[t]+heatstorage[t]) -> heatpumpflow[t]


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

using Markdown
using InteractiveUtils

# ╔═╡ b88059f0-f793-11ea-0031-23a60382d51a
using Pkg

# ╔═╡ b8809a8c-f793-11ea-382c-5763b1d11ca5
begin
	using Optim
	using StatisticalRethinking
end

# ╔═╡ fdf51ec0-f794-11ea-327a-33958a9cd52a
md"## Intro-logpdf-2s.jl"

# ╔═╡ 3330c924-f793-11ea-28b7-6d61030db6f9
md"##### This script shows clip-04-26-29s.jl using Optim and a loglik function."

# ╔═╡ b8812010-f793-11ea-19d0-61be51252ae5
md"## snippet 4.26"

# ╔═╡ 5d5852d4-34cf-11eb-2618-4b83ea072b1e
begin
	df = DataFrame(CSV.read(sr_datadir("Howell1.csv"), DataFrame; delim=';'))
	df2 = filter(row -> row[:age] >= 18, df)
end;

# ╔═╡ b88ee9c0-f793-11ea-08c9-e9b05a6e3990
md"## snippet 4.27"

# ╔═╡ 4433d922-34cf-11eb-379a-cfe1dbd939cb
# Our first model:

stan4_1 = "
  μ ~ Normal(178,20) # prior
  σ ~ Uniform(0, 50) # prior
  height ~ Normal(μ, σ) # likelihood
";

# ╔═╡ b898fb18-f793-11ea-202f-672f3b026e96
md"## snippet 4.28"

# ╔═╡ b8998394-f793-11ea-021c-bf822c791396
function loglik(x)
  ll = 0.0
  ll += log(pdf(Normal(178, 20), x[1]))
  ll += log(pdf(Uniform(0, 50), x[2]))
  ll += sum(logpdf.(Normal(x[1], x[2]), df2.height))
  -ll
end

# ╔═╡ b8a3c502-f793-11ea-34f8-eddc53070e3b
md"## snippet 4.29"

# ╔═╡ 9e16251c-34cf-11eb-2af5-21ac06afd361
# Start values

begin
	lower = [0.0, 0.0]
	upper = [250.0, 50.0]
	x0 = [170.0, 10.0]
end

# ╔═╡ b8bcaa18-f793-11ea-3e04-0fa9333a1435
res = optimize(loglik, lower, upper, x0)

# ╔═╡ b8bd8988-f793-11ea-0f93-4bf863950fe3
Optim.minimizer(res)

# ╔═╡ b8cab216-f793-11ea-1880-9dc5d6a84e79
# Our second model:

stan4_2 = "
  μ ~ Normal(178, 0.1) # prior
  σ ~ Uniform(0, 50) # prior
  height ~ Normal(μ, σ) # likelihood
";

# ╔═╡ 9c55acca-f794-11ea-26ae-c505f54777f2
function loglik2(x)
  ll = 0.0
  ll += log(pdf(Normal(178, 0.1), x[1]))
  ll += log(pdf(Uniform(0, 50), x[2]))
  ll += sum(logpdf.(Normal(x[1], x[2]), df2.height))
  -ll
end

# ╔═╡ b8dd1d66-f793-11ea-1210-99abfb4919c9
x1 = [178.0, 40.0] # Initial values can't be outside the the box.

# ╔═╡ b8ecbb40-f793-11ea-219a-4d4d89afbede
res2 = optimize(loglik2, lower, upper, x1)

# ╔═╡ d9e99bd8-34ce-11eb-1779-231c74d97cd1
Optim.minimizer(res2)

# ╔═╡ b8f4b8ec-f793-11ea-1f5e-6d1225f2da1a
md"##### Notice the increase of σ."

# ╔═╡ b8fc6e14-f793-11ea-275c-35a406bdaefc
md"## End of intro-logpdfs-2.jl"

# ╔═╡ Cell order:
# ╟─fdf51ec0-f794-11ea-327a-33958a9cd52a
# ╟─3330c924-f793-11ea-28b7-6d61030db6f9
# ╠═b88059f0-f793-11ea-0031-23a60382d51a
# ╠═b8809a8c-f793-11ea-382c-5763b1d11ca5
# ╟─b8812010-f793-11ea-19d0-61be51252ae5
# ╠═5d5852d4-34cf-11eb-2618-4b83ea072b1e
# ╠═b88ee9c0-f793-11ea-08c9-e9b05a6e3990
# ╠═4433d922-34cf-11eb-379a-cfe1dbd939cb
# ╟─b898fb18-f793-11ea-202f-672f3b026e96
# ╠═b8998394-f793-11ea-021c-bf822c791396
# ╟─b8a3c502-f793-11ea-34f8-eddc53070e3b
# ╠═9e16251c-34cf-11eb-2af5-21ac06afd361
# ╠═b8bcaa18-f793-11ea-3e04-0fa9333a1435
# ╠═b8bd8988-f793-11ea-0f93-4bf863950fe3
# ╠═b8cab216-f793-11ea-1880-9dc5d6a84e79
# ╠═9c55acca-f794-11ea-26ae-c505f54777f2
# ╠═b8dd1d66-f793-11ea-1210-99abfb4919c9
# ╠═b8ecbb40-f793-11ea-219a-4d4d89afbede
# ╠═d9e99bd8-34ce-11eb-1779-231c74d97cd1
# ╟─b8f4b8ec-f793-11ea-1f5e-6d1225f2da1a
# ╟─b8fc6e14-f793-11ea-275c-35a406bdaefc

ext: jl


text:
function Inverse()
    return Inverse(())
end

function value(obj::Inverse, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::Inverse, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
type Car
    l::Float64
    acc::Float64
    v::Float64
    vmax::Float64
    pos::Float64
    over::Bool
end

type NS
    ρ::Float64
    p::Float64
    q::Float64
    L::Float64
    dt::Float64
    lsafe::Float64
end

function setup(;l=1.0, lsafe=7.0,acc=1.0, vmax=5.0, ρ=0.1, p=0.15, q=0.0, L=1000.0, dt=1.0)
    N = round(Int, ρ*L/(l+lsafe))
    NS(ρ,p,q,L,dt,lsafe), [Car(l,acc,0.0,vmax,pos,false) for pos in sort(sample(l:(l+lsafe):L,N,replace=false))]
end

function forward!(ns::NS, car::Car, observ::Float64)
    ispass = false
    if car.pos <= observ && car.pos + ns.dt*car.v > observ
        ispass = true
    end
    car.pos += ns.dt*car.v
    if car.pos > ns.L
        car.pos -= ns.L
    end
    ispass
end
function isin!(ns::NS, car::Car, observ1::Float64,observ2::Float64)
    isin = false
    if observ1 < car.pos <= observ2
        isin = true
    end
    isin
end

function onestep!(ns::NS, cars::Vector{Car}, cars1::Vector{Car}, ns_cars_idx::Vector{Int},
    nsos_cars_idx::Vector{Int}, idx::Vector{Int})
    overcnt = 0
    passcnt = 0
    incnt=0
    N = length(cars)
    resize!(ns_cars_idx,0)
    resize!(nsos_cars_idx,0)
    for i=2:N-1
        if rand() < ns.q
            push!(nsos_cars_idx, i)
        else
            push!(ns_cars_idx, i)
        end
    end
    append!(ns_cars_idx, [1;N])

    sort!(nsos_cars_idx, rev=true)

    for i in ns_cars_idx
        ip1 = i<N ? i+1 : 1
        ip2 = i<N-1 ? i+2 : (i<N ? 1 : 2)
        if cars[i].v + ns.dt*cars[i].acc > cars[i].vmax
            cars[i].v = cars[i].vmax
        else
            cars[i].v += ns.dt*cars[i].acc
        end
        d = cars[ip1].pos - cars[i].pos - cars[ip1].l-ns.lsafe
        if d<0
            d += ns.L
        end
        if d <= cars[i].v
            cars[i].v = d
        end
        if rand()<ns.p && cars[i].v >= ns.dt*cars[i].acc
            cars[i].v -= ns.dt*cars[i].acc
        end
    end

    for i in nsos_cars_idx
        ip1 = i<N ? i+1 : 1
        ip2 = i<N-1 ? i+2 : (i<N ? 1 : 2)
        ip3 = i<N-2 ? i+3 : (i<N-1 ? 1 : 2)
        if cars[i].v + ns.dt*cars[i].acc > cars[i].vmax
            cars[i].v = cars[i].vmax
        else
            cars[i].v += ns.dt*cars[i].acc
        end
        d = cars[ip1].pos - cars[i].pos - cars[ip1].l-ns.lsafe
        if cars[ip2].over
          d1 = cars[ip3].pos + ns.dt*cars[ip3].v - cars[ip1].pos - ns.dt*cars[ip1].v-cars[ip3].l-ns.lsafe
        else
          d1 = cars[ip2].pos + ns.dt*cars[ip2].v - cars[ip1].pos - ns.dt*cars[ip1].v-cars[ip2].l-ns.lsafe
        end
        if d<0
            d += ns.L
        end
        if d1<0
            d1 += ns.L
        end
        if cars[ip1].over
            cars[i].v = min(d + ns.dt*cars[ip1].v-cars[ip2].l, cars[i].v)
            if rand()<ns.p && cars[i].v >= ns.dt*cars[i].acc
                cars[i].v -= ns.dt*cars[i].acc
            end
        elseif d + cars[i].l + cars[ip1].l + ns.dt*cars[ip1].v <= ns.dt*cars[i].v && d1 >= cars[i].l
            cars[i].v = d + cars[i].l + cars[ip1].l + ns.dt*cars[ip1].v
            cars[i].over = true
            overcnt += 1
        else
            cars[i].v=min(d+ns.dt*cars[ip1].v,cars[i].v)
            if rand()<ns.p && cars[i].v >= ns.dt*cars[i].acc
                cars[i].v -= ns.dt*cars[i].acc
            end
        end
    end

    vsum = 0.0
    for i in eachindex(cars)
        if isin!(ns, cars[i], ns.L/2, ns.L/2+0.1*ns.L)
            incnt += cars[i].l
        end
        if forward!(ns, cars[i], ns.L/2)
            passcnt += 1
            vsum += cars[i].v
        end
    end
    for i in eachindex(cars)
        cars1[i] = cars[i]
    end
    carpos = [cars[i].pos for i in eachindex(cars)]
    idx = sortperm(carpos)
    for (i, pos) in enumerate(idx)
        cars[i] = cars1[pos]
        cars[i].over = false
    end
    overcnt, passcnt,incnt, vsum
end

function spacetime(file, ns, cars; pre_T=1000,T=100)
    ns_cars_idx = Int[]
    nsos_cars_idx = Int[]
    sizehint!(ns_cars_idx, length(cars))
    sizehint!(nsos_cars_idx, length(cars))
    idx = Array(Int, length(cars))
    cars1 = copy(cars)
    Road = fill(-1, T, round(Int,ns.L))
    for i=1:pre_T
        onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
    end
    for t=1:T
        for car in cars
            Road[t, round(Int,car.pos)] = car.v
        end
        onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
    end
    f = open(file, "w")
    for i=1:size(Road,1)
        for j=1:size(Road,2)
            v = Road[i,j] < 0 ? "*" : string(Road[i,j])
            write(f, v)
        end
        write(f, '\n')
    end
    close(f)
end
function trajectory(ns, cars; pre_T=1000,T=100)
    ns_cars_idx = Int[]
    nsos_cars_idx = Int[]
    sizehint!(ns_cars_idx, length(cars))
    sizehint!(nsos_cars_idx, length(cars))
    idx = Array(Int, length(cars))
    cars1 = copy(cars)
    Road = fill(-1, T, round(Int,ns.L))
    for i=1:pre_T
        onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
    end
    for t=1:T
        for car in cars
          carend=car.pos-car.l+1
          if carend > 0
            Road[t, (round(Int,carend)):(round(Int,car.pos))] = car.v
          else
            Road[t,1:(round(Int,car.pos))] = car.v
            Road[t,(round(Int,carend+ns.L)):(round(Int,ns.L))] = car.v
          end
        end
        onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
    end
    Road
end

function draw(ns::NS,car::Car, t=1, T=100)
    (context(),rectangle(car.pos/ns.L, t/T, car.l/ns.L, 1/ns.L),fill(RGB(1-car.v/car.vmax,0,0)))
end

function trajectory!(ns, cars; T=100)
    allcars = Any[]
    ns_cars_idx = Int[]
    nsos_cars_idx = Int[]
    sizehint!(ns_cars_idx, length(cars))
    sizehint!(nsos_cars_idx, length(cars))
    idx = Array(Int, length(cars))
    cars1 = copy(cars)
      for t=1:T
          for car in cars
              push!(allcars, draw(ns,car, t, T))
          end
          onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
      end
      compose(context(), allcars...)
  end

function flowrate(ns, cars; pre_T=1000, T=2000)
    passcntsum = 0
    overcntsum = 0
    ns_cars_idx = Int[]
    nsos_cars_idx = Int[]
    sizehint!(ns_cars_idx, length(cars))
    sizehint!(nsos_cars_idx, length(cars))
    idx = Array(Int, length(cars))
    cars1 = copy(cars)
    for i=1:pre_T
        onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
    end
    for i=1:T
        #a, b = onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
        overcnt, passcnt,incnt, vsum=onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
        overcntsum += overcnt
        passcntsum += passcnt
    end
    overcntsum/T/length(cars), passcntsum/T
end

function flux!(ns, cars; pre_T=1000, T=100*60)
    density = Float64[]
    flux = Float64[]
    velocity = Float64[]
    passcnt = 0
    overcnt = 0
    ns_cars_idx = Int[]
    nsos_cars_idx = Int[]
    sizehint!(ns_cars_idx, length(cars))
    sizehint!(nsos_cars_idx, length(cars))
    idx = Array(Int, length(cars))
    cars1 = copy(cars)
    for i=1:pre_T
        onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
    end
    fluxsum = 0.0
    vsumsum = 0.0
    for i=1:T
        #a, b, c = onestep1!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
        overcnt, passcnt,incnt, vsum=onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
        fluxsum += passcnt
        vsumsum += vsum
        if i%60 == 0
            push!(density,fluxsum*fluxsum/vsumsum)
            push!(flux,fluxsum)
            push!(velocity,vsumsum/fluxsum)
            fluxsum = 0.0
            vsumsum = 0.0
        end
    end
    density, flux, velocity
end

function timeheadway(ns, cars;pre_T=1000, T=2000)
    free=Float64[]
    sync=Float64[]
    jams=Float64[]
    passcnt = 0
    overcnt = 0
    N = length(cars)
    ns_cars_idx = Int[]
    nsos_cars_idx = Int[]
    sizehint!(ns_cars_idx, length(cars))
    sizehint!(nsos_cars_idx, length(cars))
    idx = Array(Int, length(cars))
    cars1 = copy(cars)
    for i=1:pre_T
        onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
    end
    for t=1:T
        overcnt, passcnt,incnt, vsum=onestep!(ns, cars, cars1, ns_cars_idx, nsos_cars_idx, idx)
        density=incnt/(0.1*ns.L)
        for i in eachindex(cars)
            if isin!(ns, cars[i], ns.L/2,ns.L/2+50)
                ip1 = i<N ? i+1 : 1
                d = cars[ip1].pos - cars[i].pos
                if d<0
                    d += ns.L
                end
                if cars[i].v > 0
                    τ=d/cars[i].v
                    if 0<τ<10
                      density<=0.12 ? push!(free,τ) : (density >0.35 ? push!(jams,τ) : push!(sync,τ))
                      #ns.ρ<=0.1 ? push!(free,τ) : (ns.ρ >0.5 ? push!(jams,τ) : push!(sync,τ))
                    end

                end
            end
        end
    end
    free,sync,jams
end

#ns, cars = setup(p=0.25,q=0.5,ρ=0.6,vmax=5,l=1, L=1000)
#timeheadway("free1.csv","sync.csv1","jams.csv", ns, cars;pre_T=1000,T=10000)
#trajectory(ns, cars; T=100)
#ns, cars = setup(ρ=0.2, q=0.3, p=0.25, vmax=25, l=5)
#density,flux,velocity=flux!(ns, cars; pre_T=1000, T=100*60)
#using PyPlot
#plot(density,flux,"o")


text:
# https://adventofcode.com/2020/day/18

function parse_expression(chars::Vector{Char})
    expression = []
    while !isempty(chars)
        ch = popfirst!(chars)

        if ch == ' '
            continue
        elseif isdigit(ch)
            digit = parse(Int, ch)
            if !isempty(expression) && isa(expression[end], Number)
                expression[end] = 10 * expression[end] + digit
            else
                push!(expression, digit)
            end
        elseif ch == '*'
            push!(expression, :*)
        elseif ch == '+'
            push!(expression, :+)
        elseif ch == '('
            partial, chars = parse_expression(chars)
            push!(expression, partial)
        elseif ch == ')'
            return expression, chars
        else
            error("invalid symbol: '$(ch)'")
        end
    end
    return expression, chars
end

function parse_expression(string::AbstractString)
    expression, _ = parse_expression(collect(string))
    return expression
end

@assert parse_expression("1 + 2") == [1, :+, 2]
@assert parse_expression("1 + (2 * 3)") == [1, :+, [2, :*, 3]]

function evaluate(expression)
    result = 0
    op = (+)
    for x in expression
        if isa(x, Symbol)
            op = (x == :+) ? (+) : (*)
            continue
        elseif isa(x, Vector)
            x = evaluate(x)
        end
        result = op(result, x)
    end
    return result
end

parse_eval(expr) = evaluate(parse_expression(expr))

@assert parse_eval("1 + (2 * 3) + (4 * (5 + 6))") == 51
@assert parse_eval("2 * 3 + (4 * 5)") == 26
@assert parse_eval("5 + (8 * 3 + 9 + 3 * 4 * 3)") == 437
@assert parse_eval("5 * 9 * (7 * 3 * 3 + 9 * 3 + (8 + 6 * 4))") == 12240
@assert parse_eval("((2 + 4 * 9) * (6 + 9 * 8 + 6) + 6) + 2 + 4 * 2") == 13632

function part1(input)
    return reduce(+, map(parse_eval, split(input, '\n', keepempty=false)))
end

function bracket(expr)
    new_expr = []
    i = 1
    while i <= length(expr)
        x = expr[i]
        if isa(x, Vector)
            x = bracket(x)
        end

        if x == :+
            y = expr[i + 1]
            if isa(y, Vector)
                y = bracket(y)
            end

            new_expr[end] = [new_expr[end], :+, y]
            i += 2
        else
            push!(new_expr, x)
            i += 1
        end
    end
    return length(new_expr) == 1 ? first(new_expr) : new_expr
end

@assert bracket([2, :*, 3, :+, 4]) == [2, :*, [3, :+, 4]]
@assert bracket([2, :*, [3, :*, 4], :+, 5]) == [2, :*, [[3, :*, 4], :+, 5]]
@assert bracket([2, :+, 3, :+, 4, :*, 5]) == [[[2, :+, 3], :+, 4], :*, 5]
@assert bracket([5, :+, [8, :*, 3, :+, 9, :+, 3, :*, 4, :*, 3]]) == [5, :+, [8, :*, [[3, :+, 9], :+, 3], :*, 4, :*, 3]]

parse_bracket_eval(expr) = evaluate(bracket(parse_expression(expr)))

@assert parse_bracket_eval("1 + (2 * 3) + (4 * (5 + 6))") == 51
@assert parse_bracket_eval("2 * 3 + (4 * 5)") == 46
@assert parse_bracket_eval("5 + (8 * 3 + 9 + 3 * 4 * 3)") == 1445
@assert parse_bracket_eval("5 * 9 * (7 * 3 * 3 + 9 * 3 + (8 + 6 * 4))") == 669060
@assert parse_bracket_eval("((2 + 4 * 9) * (6 + 9 * 8 + 6) + 6) + 2 + 4 * 2") == 23340

function part2(input)
    return reduce(+, map(parse_bracket_eval, split(input, '\n', keepempty=false)))
end

test = read("data/day-18.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 1408133923393
@assert result2 == 314455761823725

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
module Multiplexers

function stack(io::IO,msg::Vector{UInt8})
    frontbytes = reinterpret(UInt8,Int16[length(msg)])
    item = UInt8[frontbytes...,msg...]
    write(io,item)
end

function unstack(io::IO)
    sizebytes = [read(io,UInt8),read(io,UInt8)]
    size = reinterpret(Int16,sizebytes)[1]

    msg = UInt8[]
    for i in 1:size
        push!(msg,read(io,UInt8))
    end
    return msg
end

function unstack(io::IOBuffer)
    bytes = take!(io)
    size = reinterpret(Int16,bytes[1:2])[1]
    msg = bytes[3:size+2]
    if length(bytes)>size+2
        write(io,bytes[size+3:end])
    end
    return msg
end

struct Line <: IO
    socket::IO
    ch::Channel{UInt8}
    n::UInt8
end

Line(socket,n) = Line(socket,Channel{UInt8}(Inf),UInt8(n))

import Base.write
write(line::Line,msg::UInt8) = stack(line.socket,UInt8[line.n,msg])
write(line::Line,msg::Vector{UInt8}) = stack(line.socket,UInt8[line.n,msg...])
write(line::Line,msg::String) = stack(line.socket,UInt8[line.n,msg...])

import Base.read
read(line::Line,x::Type{UInt8}) = take!(line.ch)

"""
Takes multiple input lines and routes them to a single line.
"""
function route(lines::Vector{Line},socket::IO)
    while true
        data = unstack(socket)
        if data==UInt8[0]
            stack(socket,UInt8[0])
            return
        else
            n,msg = Int(data[1]),data[2:end]

            if typeof(msg)==UInt8
                put!(lines[n].ch,msg)
            else
                bmsg = take!(IOBuffer(msg))
                for byte in bmsg #[end:-1:1]
                    put!(lines[n].ch,byte)
                end
            end
        end
    end
end

### Some high end interface
struct Multiplexer{T<:IO}
    socket::IO
    lines::Vector{T}
    daemon::Task
end

function Multiplexer(socket::IO,N::Integer)
    lines = Line[Line(socket,i) for i in 1:N]
    daemon = @async route(lines,socket)
    Multiplexer(socket,lines,daemon)
end

import Base.wait
wait(mux::Multiplexer) = wait(mux.daemon)

import Base.close
function close(mux::Multiplexer)
    stack(mux.socket,UInt8[0])
    wait(mux)
end

"""
A function which one uses to forward forward traffic from multiple sockets into one socket by multiplexing.
"""
function forward(ios::Vector{IO},socket::IO)
    mux = Multiplexer(socket,length(ios))

    tasks = []

    for (line,io) in zip(mux.lines,ios)
        task1 = @async while true
            byte = read(line,UInt8)
            write(io,byte)
        end
        push!(tasks,task1)

        task2 = @async while true
            byte = read(io,UInt8)
            write(line,byte) # One needs to think about some kind of buffering (traffic now is multiplied by 2)
        end
        push!(tasks,task2)
    end

    wait(mux)

    byte = unstack(socket)
    @assert byte==UInt8[0] "$byte is not a termination byte"

    for t in tasks
        @async Base.throwto(t,InterruptException())
    end
end

function Multiplexer(socket::IO,ios::Vector{IO})
    daemon = @async forward(ios,socket)
    Multiplexer(socket,ios,daemon)
end

export Multiplexer, Line, read, write

end # module


text:
function get_chain(G::GAP.GapObj)
  @assert GAP.Globals.IsSolvable(G)
  if GAP.Globals.IsAbelian(G)
    return [GAP.Globals.TrivialSubgroup(G), G]
  end
  lN = GAP.Globals.CharacteristicSubgroups(G)
  acceptable = trues(length(lN))
  orders = [GAP.Globals.Size(lN[i]) for i = 1:length(lN)]
  for i = 1:length(lN)
    if !acceptable[i]
      continue
    end
    for j = i+1:length(lN)
      if orders[j] == orders[i]
        id_i = GAP.Globals.IdGroup(lN[i])
        id_j = GAP.Globals.IdGroup(lN[j])
        if id_i == id_j
          id_fi = GAP.Globals.IdGroup(GAP.Globals.FactorGroup(G, lN[i]))
          id_fj = GAP.Globals.IdGroup(GAP.Globals.FactorGroup(G, lN[j]))
          if id_fi == id_fj
            acceptable[i] = false
            acceptable[j] == false
            break
          end
        end
      end
    end
  end
  subs_for_series = [lN[i] for i = 1:length(lN) if acceptable[i]]
  res = [GAP.Globals.TrivialSubgroup(G)]
  while res[end] != G
    contain = [x for x in subs_for_series if GAP.Globals.IsSubset(x, res[end])]
    abelian_subgroups = [x for x in contain if GAP.Globals.IsAbelian(GAP.Globals.FactorGroup(x, res[end]))]
    j = 0
    ord = 0
    for i = 1:length(abelian_subgroups)
      x = abelian_subgroups[i]
      sx = GAP.Globals.Size(x)
      if sx > ord
        j = i
        ord = sx
      end
    end
    push!(res, abelian_subgroups[j])
  end
  return res
end

function get_chain(G::GAP.GapObj, H::GAP.GapObj)
  lN = GAP.Globals.NormalSubgroups(G)
  possible_subs = []
  idH = GAP.Globals.IdGroup(H)
  for i = 1:length(lN)
    if GAP.Globals.IdGroup(GAP.Globals.FactorGroup(G, lN[i])) == idH
      push!(possible_subs, lN[i])
    end
  end
  if isempty(possible_subs)
    error("The field can't be extended!")
  end
  res = possible_subs[1]
  if length(possible_subs) > 1
    csubs = [GAP.Globals.IdGroup(x) for x in possible_subs]
    best_occ = length(csubs)
    for i = 1:length(csubs)
      occ = 0
      for j = 1:length(csubs)
        if csubs[i] == csubs[j]
          occ += 1
        end
      end
      if occ == 1
        res = possible_subs[i]
        break
      elseif occ < best_occ
        res = possible_subs[i]
      end
    end
  end
  mp = GAP.Globals.NaturalHomomorphismByNormalSubgroup(G, res)
  K = GAP.Globals.Kernel(mp)
  return get_chain(K)
end


text:
#######################################################################
# Diagonalize methods
#   (All but LinearAlgebraPackage `@require` some package to be loaded)
#######################################################################
abstract type AbstractDiagonalizeMethod end

struct HamiltonianBlochFunctor{H<:Union{Hamiltonian,ParametricHamiltonian},A<:AbstractMatrix,M}
    h::H
    matrix::A
    mapping::M
end

(f::HamiltonianBlochFunctor)(vertex) = bloch!(f.matrix, f.h, map_phiparams(f.mapping, vertex))

struct Diagonalizer{M<:AbstractDiagonalizeMethod,F,O<:Union{OrbitalStructure,Missing}}
    method::M
    matrixf::F         # functor or function matrixf(φs) that produces matrices to be diagonalized
    orbstruct::O       # store structure of original Hamiltonian if available (to allow unflattening eigenstates)
    perm::Vector{Int}  # reusable permutation vector
    perm´::Vector{Int} # reusable permutation vector
end

struct NoUnflatten end

"""
    diagonalizer(h::Union{Hamiltonian,ParametricHamiltonian}; method = LinearAlgebraPackage(), mapping = missing)

Build a `d::Diagonalizer` object that, when called as `d(φs)` , uses the specified
diagonalization `method` to produce the sorted eigenpairs `(εs, ψs)` of `h` at Bloch
phases/parameters given by `mapping`. See `bandstructure` for further details.

A 0D Hamiltonian `h` also supports `d = diagonalizer(h)`. In this case `d` can be called
with no arguments and gives the same information as `spectrum`, `d() ≈ Tuple(spectrum(h))`.

# Examples
```jldoctest
julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(1));

julia> d = diagonalizer(h)
Diagonalizer with method : LinearAlgebraPackage{NamedTuple{(), Tuple{}}}

julia> d((0, 0)) |> first
2-element Vector{Float64}:
 -3.0
  3.0

julia> h = wrap(h); d = diagonalizer(h);

julia> d() .≈ Tuple(spectrum(h))
(true, true)
```

# See also
    `bandstructure`, `spectrum`
"""
function diagonalizer(h::Union{Hamiltonian,ParametricHamiltonian}; method = LinearAlgebraPackage(), mapping = missing)
    matrix = similarmatrix(h, method_matrixtype(method, h))
    matrixf = HamiltonianBlochFunctor(h, matrix, mapping)
    perm = Vector{Int}(undef, size(matrix, 2))
    orbstruct = parent(h).orbstruct
    return Diagonalizer(method, matrixf, orbstruct, perm, copy(perm))
end

function diagonalizer(matrixf::Function, dimh; method = LinearAlgebraPackage())
    perm = Vector{Int}(undef, dimh)
    return Diagonalizer(method, matrixf, missing, perm, copy(perm))
end

@inline function (d::Diagonalizer)(φs, ::NoUnflatten)
    ϵ, ψ = diagonalize(d.matrixf(φs), d.method)
    issorted(ϵ, by = real) || sorteigs!(ϵ, ψ, d)
    fixphase!(ψ)
    return ϵ, ψ
end

@inline (d::Diagonalizer)(n::NoUnflatten) = d((), n)

function (d::Diagonalizer)(φs)
    ϵ, ψ = d(φs, NoUnflatten())
    ψ´ = unflatten_orbitals_or_reinterpret(ψ, d.orbstruct)
    return ϵ, ψ´
end

@inline (d::Diagonalizer)() = d(())

function sorteigs!(ϵ::AbstractVector, ψ::AbstractMatrix, d)
    p, p´ = d.perm, d.perm´
    resize!(p, length(ϵ))
    resize!(p´, length(ϵ))
    sortperm!(p, ϵ, by = real, alg = Base.DEFAULT_UNSTABLE)
    Base.permute!!(ϵ, copy!(p´, p))
    Base.permutecols!!(ψ, copy!(p´, p))
    return ϵ, ψ
end

fixphase!(ψ::AbstractArray{<:Real}) = ψ
function fixphase!(ψ::AbstractArray{<:Complex})
    for col in eachcol(ψ)
        col ./= cis(angle(mean(col)))
    end
    return ψ
end

function Base.show(io::IO, d::Diagonalizer)
    print(io, "Diagonalizer with method : $(summary(d.method))")
end

## Diagonalize methods ##

checkloaded(package::Symbol) = isdefined(Main, package) ||
    throw(ArgumentError("Package $package not loaded, need to be `using $package`."))

## LinearAlgebra ##
struct LinearAlgebraPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
    kw::K
end

LinearAlgebraPackage(; kw...) = LinearAlgebraPackage(values(kw))

function diagonalize(matrix, method::LinearAlgebraPackage)
    ϵ, ψ = eigen!(matrix; (method.kw)...)
    return ϵ, ψ
end

## Arpack ##
struct ArpackPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
    kw::K
end

ArpackPackage(; kw...) = (checkloaded(:Arpack); ArpackPackage(values(kw)))

function diagonalize(matrix, method::ArpackPackage)
    ϵ, ψ = Main.Arpack.eigs(matrix; (method.kw)...)
    return ϵ, ψ
end

## ArnoldiMethod ##
struct ArnoldiMethodPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
    kw::K
end

ArnoldiMethodPackage(; kw...) = (checkloaded(:ArnoldiMethod); ArnoldiMethodPackage(values(kw)))

function diagonalize(matrix, method::ArnoldiMethodPackage)
    p = first(Main.ArnoldiMethod.partialschur(matrix; (method.kw)...))
    ϵ, ψ = p.eigenvalues, p.Q
    return ϵ, ψ
end

## IterativeSolvers ##

struct KrylovKitPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
    kw::K
end

KrylovKitPackage(; kw...) = (checkloaded(:KrylovKit); KrylovKitPackage(values(kw)))

function diagonalize(matrix::AbstractMatrix{M}, method::KrylovKitPackage) where {M}
    ishermitian(matrix) || throw(ArgumentError("Only Hermitian matrices supported with KrylovKitPackage for the moment"))
    origin = get(method.kw, :origin, 0.0)
    howmany = get(method.kw, :howmany, 1)
    kw´ = Base.structdiff(method.kw, NamedTuple{(:origin,:howmany)}) # Remove origin option
    lmap = shiftandinvert(matrix, origin)
    T = eltypevec(matrix)
    n = size(matrix, 2)
    x0 = rand(T, n)
    ϵ, ψ, _ = Main.KrylovKit.eigsolve(x -> lmap * x, x0, howmany, :LM; kw´...)

    ϵ´ = invertandshift(ϵ, origin)
    resize!(ϵ´, howmany)

    dimh = size(matrix, 2)
    ψ´ = Matrix{T}(undef, dimh, howmany)
    for i in 1:howmany
        copyslice!(ψ´, CartesianIndices((1:dimh, i:i)), ψ[i], CartesianIndices((1:dimh,)))
    end

    return ϵ´, ψ´
end

### matrix/vector types

similarmatrix(h, method::AbstractDiagonalizeMethod) = similarmatrix(h, method_matrixtype(method, h))

method_matrixtype(::LinearAlgebraPackage, h) = Matrix{blockeltype(h)}
method_matrixtype(::AbstractDiagonalizeMethod, h) = flatten

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
#Function for calculating Rydberg blockade shifts
#24/07/2017

function BlockadeShift(atom,nn,ll,jj,mj)

  pceV = 2*π*rpceV #Planck's constant in eV
  ZZ, spin, eneigval, params = GetAtomParams(atom,nn,ll,jj)

  #Create ranges for θ and R (R in atomic units)
  θ = collect(0:0.01:(π/2))
  RRSI = (collect(4:0.01:10))*1e-6
  RR = RRSI/bohrrad

  #Set defect threshold
  defectmax = 100e9 #Maximum energy defect in Hz
  entol = defectmax*pceV/atomenergy #Converted to atomic units
  ntol = 4 #Maximum change in n

  #Determine single particle state energy
  singen = eneigval

  #Create set of allowable interacting state quantum numbers
  nvec = rectpulse((nn-ntol):(nn+ntol),ll+2)
  lcands = collect((ll-1):2:(ll+1))
  lcands = filter(lcands->lcands≥0,lcands)
  if lcands[1]==0
    smalllvec = [0; rectpulse(lcands[2:end],2)]
    jmaker = [0.5; repmat([-0.5, 0.5],1,length(lcands[2:end]))]
  else
    smalllvec = rectpulse(lcands,2)
    jmaker = repmat([-0.5, 0.5],1,length(lcands))
  end
  jmaker = reshape(jmaker,length(smalllvec),1)
  smalljvec = smalllvec+jmaker
  lvec = repmat(smalllvec,1,Int(length(nvec)/length(smalllvec)))
  jvec = repmat(smalljvec,1,Int(length(nvec)/length(smalljvec)))
  lvec = reshape(lvec,length(nvec),1)
  jvec = reshape(jvec,length(nvec),1)
  truncspace1 = hcat(nvec, lvec, jvec)'
  pairs1 = combvechack(truncspace1,truncspace1)
  Sp1 = size(pairs1)

  #Check which of these states satisfy the infinite separation energy defect condition and selection rules for j

  pindex = zeros(Sp1[2])
  infdefects = zeros(Sp1[2])

  for kk = 1:Sp1[2]
    energy1 = envalfunc(atom,pairs1[1,kk],pairs1[2,kk],pairs1[3,kk])
    energy2 = envalfunc(atom,pairs1[4,kk],pairs1[5,kk],pairs1[6,kk])
    infdefects[kk] = energy1 + energy2 - 2*singen #Energy defect
    jvec = collect((jj-1):(jj+1))
    j1check = length(find(jvec -> jvec==pairs1[3,kk],jvec))
    j2check = length(find(jvec -> jvec==pairs1[6,kk],jvec))
    if abs(infdefects[kk])≤entol&&j1check==1&&j2check==1
      pindex[kk] = 1
    end
  end
  qvals = find(pindex -> pindex==1,pindex)
  pairs2 = vcat(pairs1[:,qvals],infdefects[qvals]')
  Lp2 = Int(length(infdefects[qvals]))

  matel2part = zeros(length(θ),Lp2) #Array to store matrix elements in

  #Call function to calculate matrix elements
  for kk = 1:Lp2
    matel2part[:,kk] = matrixel2p(atom,nn,ll,jj,mj,pairs2[1,kk],pairs2[2,kk],pairs2[3,kk],pairs2[4,kk],pairs2[5,kk],pairs2[6,kk],pairs2[7,kk],θ)
  end

  #Compute blockade shift in atomic units
  summation = -sum(matel2part,2)
  maxs = maximum(abs(summation))
  kindex = find(summation -> abs(summation)==maxs,summation) #Finds index of maximum blockade shift

  blockadeshiftau = summation[kindex]./(RR.^6)

  RRmesh = repmat(RR,1,length(summation))'
  summationmesh = repmat(summation,1,length(RR))
  blockadeshiftmesh = summationmesh./(RRmesh.^6)

  #Convert to GHz
  encon = (atomenergy/pceV)*1e-9 #Factor from atomic energy units to GHz
  blockadeshiftGHz = blockadeshiftau*encon
  blockadeshiftmeshGHz = blockadeshiftmesh*encon

  #Curve fitting
  C_6val= blockadeshiftGHz[end]*(RRSI[end]^6) #GHz/m⁶

  return RRSI, θ, blockadeshiftmeshGHz, C_6val
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
"""
$(TYPEDEF)

Introduce direct observations on all dimensions of a Pose2 variable:

Example:
--------
```julia
PriorPose2( MvNormal([10; 10; pi/6.0], Matrix(Diagonal([0.1;0.1;0.05].^2))) )
```
"""
mutable struct PriorPose2{T} <: IncrementalInference.FunctorSingleton  where {T <: IncrementalInference.SamplableBelief}
    Z::T
    PriorPose2{T}() where T = new{T}()
    PriorPose2{T}(x::T) where {T <: IncrementalInference.SamplableBelief}  = new{T}(x)
end
PriorPose2(x::T) where {T <: IncrementalInference.SamplableBelief} = PriorPose2{T}(x)

function getSample(p2::PriorPose2, N::Int=1)
  return (rand(p2.Z,N), )
end

## Serialization support

"""
$(TYPEDEF)
"""
mutable struct PackedPriorPose2  <: IncrementalInference.PackedInferenceType
    str::String
    PackedPriorPose2() = new()
    PackedPriorPose2(x::String) = new(x)
end
function convert(::Type{PackedPriorPose2}, d::PriorPose2)
  return PackedPriorPose2(string(d.Z))
end
function convert(::Type{PriorPose2}, d::PackedPriorPose2)
  distr = extractdistribution(d.str)
  return PriorPose2{typeof(distr)}(distr)
end

## NOTE likely deprecated comparitors, see DFG compareFields, compareAll instead
function compare(a::PriorPose2,b::PriorPose2; tol::Float64=1e-10)
  TP = true
  TP = TP && norm(a.Z.μ-b.Z.μ) < tol
  TP = TP && norm(a.Z.Σ.mat-b.Z.Σ.mat) < tol
  return TP
end


text:
using DataFrames
#using UrlDownload
using CSV, HTTP
import Downloads
using Test
using TerminalPager
using PrettyTables
using Formatting

# See learn_Downloads.jl for how to download

"""
This is not supported in 1.0

See learn_empty_constructor()
"""
function unsupported_DataFrame()
    @test_throws ArgumentError DataFrame(
            [String, Float64, String],
            ["name", "score", "note"]
    )
end

function learn_constructor()

    # Like this one the most
    df = DataFrame( # Pair constructor
        "a" => 1:2,
        "b" => 0  , # auto fill
    )
    # 2×2 DataFrame
    #  Row │ a      b
    #      │ Int64  Int64
    # ─────┼──────────────
    #    1 │     1      0
    #    2 │     2      0
    # DataFrame([:a => 1:2, :b => 13]) # vector of Pairs constructor
    @test [1,2] == df[!,:a]
    @test [0,0] == df[!,:b]

    df = DataFrame((a=[1, 2], b=[3, 4])) # Tables.jl table constructor
    # 2×2 DataFrame
    #  Row │ a      b
    #      │ Int64  Int64
    # ─────┼──────────────
    #    1 │     1      3
    #    2 │     2      4
    @test [1,2] == df[!,:a]
    @test [3,4] == df[!,:b]

    df = DataFrame([(a=1, b=0), (a=2, b=0)]) # Tables.jl table constructor
    # 2×2 DataFrame
    #  Row │ a      b
    #      │ Int64  Int64
    # ─────┼──────────────
    #    1 │     1      0
    #    2 │     2      0
    @test [1,2] == df[!,:a]
    @test [0,0] == df[!,:b]

    df = DataFrame([:a => 1:2, :b => 0]) # vector of Pairs constructor
    # 2×2 DataFrame
    #  Row │ a      b
    #      │ Int64  Int64
    # ─────┼──────────────
    #    1 │     1      0
    #    2 │     2      0
    @test [1,2] == df[!,:a]
    @test [0,0] == df[!,:b]

    df = DataFrame(Dict(:a => 1:2, :b => 0)) # dictionary constructor
    #  2×2 DataFrame
    #   Row │ a      b
    #       │ Int64  Int64
    #  ─────┼──────────────
    #     1 │     1      0
    #     2 │     2      0
    @test [1,2] == df[!,:a]
    @test [0,0] == df[!,:b]

    df = DataFrame(a=1:2, b=0) # keyword argument constructor
    #  2×2 DataFrame
    #   Row │ a      b
    #       │ Int64  Int64
    #  ─────┼──────────────
    #     1 │     1      0
    #     2 │     2      0
    @test [1,2] == df[!,:a]
    @test [0,0] == df[!,:b]

    df = DataFrame([[1, 2], [0, 0]], [:a, :b]) # vector of vectors constructor
    #  2×2 DataFrame
    #   Row │ a      b
    #       │ Int64  Int64
    #  ─────┼──────────────
    #     1 │     1      0
    #     2 │     2      0
    @test [1,2] == df[!,:a]
    @test [0,0] == df[!,:b]

    df = DataFrame([[1, 2], [0, 0]], ["a", "b"]) # vector of vectors constructor
    #  2×2 DataFrame
    #   Row │ a      b
    #       │ Int64  Int64
    #  ─────┼──────────────
    #     1 │     1      0
    #     2 │     2      0
    @test [1,2] == df[!,:a]
    @test [0,0] == df[!,:b]
    @test [1,2] == df[!,"a"]
    @test [0,0] == df[!,"b"]

    df = DataFrame([1 0; 2 0], :auto) # matrix constructor
    #  2×2 DataFrame
    #   Row │ x1     x2
    #       │ Int64  Int64
    #  ─────┼──────────────
    #     1 │     1      0
    #     2 │     2      0
    @test [1,2] == df[!,:x1]
    @test [0,0] == df[!,:x2]
end

"""
Create a empty DataFrame then fill row-by-row.
"""
function learn_empty_constructor()
    # Create empty data frame
    df = DataFrame(
        [Float64[], String[], Int64[]],  # column data type
        ["cond", "name", "count"]        # column names
    )

    # Add rows
    push!(df, [13.0, "H-1", 5])
    push!(df, [17.0, "H-2", 7])

    # Test
    @test [ 13.0,  17.0] == df[!, "cond"]
    @test ["H-1", "H-2"] == df[!, "name"]
    @test [    5,     7] == df[!, :count]
end

function learn_csv()
    url = "https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data-original"

    # Download to a file
    tmp = Downloads.download(url)
    csv = CSV.File(
        tmp,
        delim = ' ',
        ignorerepeated = true,
        header = [
            :mpg, :cylinders, :displacement, :horsepower,
            :weight, :acceleration, :year, :origin, :name,
        ],
        missingstring = "NA",
    )
    rm(tmp)
    df = DataFrame(csv)

    return df
end

"""

https://stackoverflow.com/questions/51240161/how-to-insert-a-column-in-a-julia-dataframe-at-specific-position-without-referr
"""
function learn_add_column()
    df = DataFrame(
        "name"  => ["Liam", "Sophie", "Jacob"],
        "score" => [8.0, 8.0, 12.0],
    )

    insertcols!(
        df,
        3,
        "note" => ["4.17 below the mean", "4.17 below the mean", "0.17 below the mean"],
        makeunique = true
    )

    @test ["4.17 below the mean", "4.17 below the mean", "0.17 below the mean"] == df[!,"note"]
end

function learn_get_column()
    df = DataFrame(
        "name"  => ["Liam", "Sophie", "Jacob"],
        "score" => [8.0, 8.0, 12.0],
        "note"  => ["4.17 below the mean", "4.17 below the mean", "0.17 below the mean"]
    )

    @test ["Liam", "Sophie", "Jacob"] == df[!, "name"] # by ref
    @test ["Liam", "Sophie", "Jacob"] == df[:, "name"] # by value
end

function learn_groupby_1()
    mdf = DataFrame(
        [
            [1,        2,    3],        # A
            [true, false, true],        # B
            [1,    nothing,  3],        # C
            [nothing, nothing, nothing] # D
        ],
        ["A", "B", "C", "D"]
    )

    # group by :B
    gdf = groupby(mdf, "B")
    df = gdf[(true,)] # group of B == true
    #pager(df)
    df[:,"D"] = df[!,"A"] ./ df[!,"C"]
    #pager(mdf)
end

function learn_hcat() # https://dataframes.juliadata.org/stable/lib/functions/#Base.hcat
    df1 = DataFrame(A=1:3, B=1:3)
    df2 = DataFrame(A=4:6, B=4:6)

    df3 = hcat(df1, df2, makeunique=true) # copycols = true
    @test df3.A !== df1.A
    df3 = hcat(df1, df2, makeunique=true, copycols=false) # by ref so ...
    @test df3.A === df1.A

    @test_throws ArgumentError hcat(df1, df2, makeunique=false) # Duplicate variable names

    df1 = DataFrame(A=1:3, B=1:3)
    df2 = DataFrame(A=4:6,        C=4:6, D=5:7)
    df3 = hcat(df1, df2[!, [2, 3]])
    @test df3.A == df1.A
    @test df3.C == df2.C
end

function learn_row_iteration()
    df = simple_table()
    for row in eachrow(df)
        # @show typeof(row)
        # @show row["A"]
    end
end

function learn_set_value()
    df = simple_table()
    @test 1.0 == df[1,"A"]
    df[1,"A"] = 13.0 # How efficient is this?
    @test 13.0 == df[1,"A"]

    # Row is referenced to internal data
    df = simple_table()
    row = df[1,:]
    @test 1 == row["A"]
    row["A"] = 13.0
    @test 13 == row["A"]
    @test 13 == df[1,"A"] # referenced
end

function learn_sort()
    mdf = DataFrame(
        [
            Int64[1, 2, 3], # A
            Int64[3, 2, 1], # B
            Int64[3, 1, 2], # C
        ],
        ["A", "B", "C"]
    )

    df = sort(mdf, ["B"])
    @test [1,2,3] == df[!,"B"]
    @test [2,1,3] == df[!,"C"]
end

function learn_transform_1()
    #pager(simple_table())
    # 10×3 DataFrame
    # Row │ A        B        C
    #     │ Float64  Float64  Float64
    #─────┼───────────────────────────
    #   1 │     1.0     11.0    101.0
    #   2 │     2.0     12.0    102.0
    #   3 │     3.0     13.0    103.0
    #   4 │     4.0     14.0    104.0
    #   5 │     5.0     15.0    105.0
    #   6 │     6.0     16.0    106.0
    #   7 │     7.0     17.0    107.0
    #   8 │     8.0     18.0    108.0
    #   9 │     9.0     19.0    109.0
    #  10 │    10.0     20.0    110.0

    # Assign same value back to the same column
    df = simple_table()
    transform!(df, "A" => x -> x, renamecols=false )
    @test simple_table() == df

    # Assign same value to a new column
    # Notice the () around the anonymous function
    df = simple_table()
    transform!(df, "A" => (x -> x) => "AA")
    @test df[!,"A"] == df[!,:AA]

    # Create A2 = 2*A
    df = simple_table()
    transform!(df, "A" => (x -> x*2) => "A2")
    @test df[!,"A"] .* 2 == df[!, "A2"]

    # Create ABC = A + B + C
    # Notice the () around ((a,b,c) -> a+b+c), this is a must
    df = simple_table()
    transform!(df, ["A", "B", "C"] => ((a,b,c) -> a+b+c) => "ABC")
    @test df[!, "A"] + df[!, "B"] + df[!, "C"] == df[!,"ABC"]

    # Same but use a function
    df = simple_table()
    transform!(df, ["A", "B", "C"] =>
        function(a,b,c)
            return a+b+c
        end
        => "ABC"
    )
    @test df[!, "A"] + df[!, "B"] + df[!, "C"] == df[!,"ABC"]

    # Create D = min(A,B)
    df = simple_table()
    df[:,"A2"] = [1.2, 1.1, 3.0, 4.3, 5.0, 5.9, 7.2, 7.8, 9.3, 9.4]
    transform!(df, ["A", "A2"] => ((a,b) -> min.(a,b)) => "D")
    @show df
    transform!(df, ["A", "A2"] => (ByRow((a,b)->min(a,b))) => "D")
    @show df
    transform!(df, ["A"] => ((a) -> min.(a,5)) => "D")
    @show df
end

function learn_filter_1()
    #pager(simple_table())
    # 10×3 DataFrame
    # Row │ A        B        C
    #     │ Float64  Float64  Float64
    #─────┼───────────────────────────
    #   1 │     1.0     11.0    101.0
    #   2 │     2.0     12.0    102.0
    #   3 │     3.0     13.0    103.0
    #   4 │     4.0     14.0    104.0
    #   5 │     5.0     15.0    105.0
    #   6 │     6.0     16.0    106.0
    #   7 │     7.0     17.0    107.0
    #   8 │     8.0     18.0    108.0
    #   9 │     9.0     19.0    109.0
    #  10 │    10.0     20.0    110.0

    df = simple_table()
    filter!(row -> row["A"] <=5, df)
    @test (5,3) == size(df)

    # Same but use a function, use this for longer function
    df = simple_table()
    filter!(function(row)
        return row["A"] <= 5
    end, df)
    @test (5,3) == size(df)

    # Keep 3 <= A <= 7
    df = simple_table()
    filter!(row -> 3 <= row["A"] <= 7, df)
    @test (5,3) == size(df)
    @test [3,4,5,6,7] == df[!,"A"]

    # Remove 3 <= A <= 7
    df = simple_table()
    filter!(row -> !(3 <= row["A"] <= 7), df)
    @test (5,3) == size(df)
    @test [1,2,8,9,10] == df[!,"A"]
    @test [11,12,18,19,20] == df[!,"B"]
    @test [101,102,108,109,110] == df[!,"C"]

end
"""

Printing and formatting with PrettyTables
"""
function learn_formatting()
    formatter = function(v, i, j)
        if j in [2]
            return sprintf1("%'7.0f", v)
        else
            return v
        end
    end

    df = DataFrame(Dict(
        "a" => Float64[1000000, 2000000],
        "b" => Float64[10000000, 20000000],
    ))

    pretty_table(df, formatters=formatter)

end

function simple_table()
    df = DataFrame(
        [1.0:1:10, 11.0:1:20, 101.0:1:110],
        ["A", "B", "C"]
    )

    return df
end

@testset "Base" begin
    unsupported_DataFrame()

    learn_add_column()
    learn_get_column()

    learn_constructor()
    learn_empty_constructor()
    #learn_csv()

    learn_groupby_1()
    learn_hcat()

    learn_row_iteration()
    learn_set_value()
    learn_sort()
    learn_transform_1()
    learn_filter_1()
end

#df = learn_csv()

nothing


text:
#jl #! format: off
# # Simulations with TAMU data and [PowerSimulations.jl](https://github.com/NREL/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# This is a basic simulation example using the [TAMU Cases](https://electricgrids.engr.tamu.edu/).

# ## Dependencies
using SIIPExamples
pkgpath = pkgdir(SIIPExamples)
using Dates

# ### Modeling Packages
using PowerSystems
using PowerSimulations
using PowerSystemCaseBuilder

# ### Optimization packages
# For this simple example, we can use the HiGHS solver with a relatively relaxed tolerance.
using HiGHS # mip solver
solver = optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.5)

# ### Create a `System` from TAMU data
# We have included some of the TAMU cases (with truncated time series info)
# in the PowerSystemCaseBuilder testing, so we can just use that.
sys = build_system(PSYTestSystems, "tamu_ACTIVSg2000_sys")
transform_single_time_series!(sys, 24, Dates.Hour(24));

# ## Run a PCM
# note that the TAMU data doesn't contain startup and shutdown costs, or minimum up/down
# time limits, so a UC problem merely respects minmum generation levels.

sim_folder = mkpath(joinpath(pkgpath, "TAMU-sim"))
models = SimulationModels(
    decision_models = [
        DecisionModel(template_unit_commitment(), sys, name = "UC", optimizer = solver),
    ],
)
# ### Define and build a simulation
sim = Simulation(
    name = "TAMU-test",
    steps = 3,
    models = models,
    sequence = SimulationSequence(
        models = models,
        ini_cond_chronology = InterProblemChronology(),
    ),
    simulation_folder = sim_folder,
)

build!(sim)

# ### Execute the simulation
execute!(sim)

# ### Load and analyze results
sim_results = SimulationResults(sim);
uc_results = get_decision_problem_results(sim_results, "UC")
read_realized_variable(uc_results, "OnVariable__ThermalStandard")

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
mutable struct CmdDependency
    exec::Base.Cmd
    test_args::Base.Cmd
	validate_success::Bool
    validate_stdout::Function
    validate_stderr::Function
    exit_when_fail::Bool
end

"""
# Struct

	mutable struct CmdDependency
		exec::Base.Cmd
		test_args::Base.Cmd
		validate_success::Bool
		validate_stdout::Function
		validate_stderr::Function
		exit_when_fail::Bool
	end

# Methods

	CmdDependency(;
		exec::Base.Cmd=``,
		test_args::Base.Cmd=``,
		validate_success::Bool=false,
		validate_stdout::Function=do_nothing,
		validate_stderr::Function=do_nothing,
		exit_when_fail::Bool=true
	)

Create Command Dependency (`CmdDependency`).

# Arguments

- `exec::AbstractCmd`: the command to call the dependency.

- `test_args::AbstractCmd`: for testing purposes, the command to be appended to `exec`.

- `validate_success::Bool`: when checking the dependency, whether to validate the exit code == 0.

- `validate_stdout::Function`: a function takes standard out as `String` and return the validation result as `::Bool`.

- `validate_stderr::Function`: a function takes standard error as `String` and return the validation result as `::Bool`.

- `exit_when_fail::Bool`: if validation fails, whether to throw error and exit.

# Example

	julia = CmdDependency(
		exec = Base.julia_cmd(),
		test_args = `--version`,
		validate_success = true,
		validate_stdout = x -> occursin(r"^julia version", x),
		validate_stderr = do_nothing,
		exit_when_fail = true
	)

	check_dependency(julia)

"""
function CmdDependency(;
	exec::Cmd=``,
	test_args::Cmd=``,
	validate_success::Bool=false,
	validate_stdout::Function=do_nothing,
	validate_stderr::Function=do_nothing,
	exit_when_fail::Bool=true
	)
	CmdDependency(exec, test_args, validate_success, validate_stdout, validate_stderr, exit_when_fail)
end

"""
	readall(cmd::Cmd) -> (standard_out::String, standard_err::String, success::Bool)
"""
function readall(cmd::Base.AbstractCmd)
    out = Pipe()
    err = Pipe()

	process = try
	    run(pipeline(cmd, stdout=out, stderr=err))
	catch
		nothing
	end
    close(out.in)
    close(err.in)

    return (
        read(out, String),
        read(err, String),
        !isnothing(process)
    )
end

"""
	check_dependency(p::CmdDependency) -> Bool

Check `CmdDependency` by evaluating:

	`\$(p.exec) \$(p.test_args)`

If success, return `true`.

If fail, return `false`, or throw DependencyError when `p.exit_when_fail` set to `true`.
"""
function check_dependency(p::CmdDependency)
	if isempty(p.test_args)
		@warn "Cannot check dependency for test_args not set: $p"
		return true
	end
    out, err, success = readall(`$p $(p.test_args)`)
	has_dependency = true

	if p.validate_success && !success
		@goto dependency_error
	end

	try
		res = p.validate_stdout(out)
		res === false && error("DependencyError: invalid: $p: the `validate_stdout` function returns false.")
	catch e
		rethrow(e)
		@goto dependency_error
	end

	try
		res = p.validate_stderr(err)
		res === false && error("DependencyError: invalid: $p: the `validate_stdout` function returns false.")
	catch e
		rethrow(e)
		@goto dependency_error
	end

	return true

	@label dependency_error
    @error timestamp() * "DependencyError: invalid: $p" CHECK_ARGS=p.test_args
	println(stderr, "Dependency check stdout:")
	println(stderr, out)
	println(stderr, "Dependency check stderr:")
	println(stderr, err)
	if p.exit_when_fail
		error("DependencyError: invalid: $p")
	end
	return false
end

# Interpolation in Cmd
# It allows CmdDependency to be interpolated in `$p`.
Base.arg_gen(p::CmdDependency) = Base.arg_gen(p.exec)

"""
	check_dependency_file(path::Union{AbstractString,Cmd}; exit_when_false=true) -> Bool

Checke whether a file exists. Return `::Bool`.
"""
function check_dependency_file(path::AbstractString; exit_when_false=true)
	has_dependency = isfile(path)
	if exit_when_false && !has_dependency
		error("DependencyError: not a valid file: $path")
	end
	has_dependency
end
check_dependency_file(path::Cmd; exit_when_false=true) = check_dependency_file(str(path); exit_when_false=exit_when_false)

"""
	check_dependency_dir(path::Union{AbstractString,Cmd}; exit_when_false=true) -> Bool

Checke whether a directory exists. Return `::Bool`.
"""
function check_dependency_dir(path::AbstractString; exit_when_false=true)
	has_dependency = isdir(path)
	if exit_when_false && !has_dependency
		error("DependencyError: not a valid directory: $path")
	end
	has_dependency
end
check_dependency_dir(path::Cmd; exit_when_false=true) = check_dependency_file(str(path); exit_when_false=exit_when_false)


text:
module SeaPearlExtras

using SeaPearl
using Plots
using GraphPlot
using CSV
using BSON
using DataFrames

include("metrics/metrics.jl")
include("staterepresentation/representation.jl")

end # module


text:
# SAM Record
# ==========

type SAMRecord
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # indexes
    qname::UnitRange{Int}
    flag::UnitRange{Int}
    rname::UnitRange{Int}
    pos::UnitRange{Int}
    mapq::UnitRange{Int}
    cigar::UnitRange{Int}
    rnext::UnitRange{Int}
    pnext::UnitRange{Int}
    tlen::UnitRange{Int}
    seq::UnitRange{Int}
    qual::UnitRange{Int}
    fields::Vector{UnitRange{Int}}
end

"""
    SAMRecord()

Create an unfilled `SAMRecord` object.
"""
function SAMRecord()
    return SAMRecord(
        UInt8[], 1:0,
        # qname-mapq
        1:0, 1:0, 1:0, 1:0, 1:0,
        # cigar-seq
        1:0, 1:0, 1:0, 1:0, 1:0,
        # qual and fields
        1:0, UnitRange{Int}[])
end

"""
    SAMRecord(data::Vector{UInt8})

Create a `SAMRecord` object from `data` containing a SAM record.
This function verifies the format and indexes fields for accessors.
Note that the ownership of `data` is transferred to a new `SAMRecord` object.
"""
function SAMRecord(data::Vector{UInt8})
    return convert(SAMRecord, data)
end

function Base.convert(::Type{SAMRecord}, data::Vector{UInt8})
    record = SAMRecord(
        data, 1:0,
        # qname-mapq
        1:0, 1:0, 1:0, 1:0, 1:0,
        # cigar-seq
        1:0, 1:0, 1:0, 1:0, 1:0,
        # qual and fields
        1:0, UnitRange{Int}[])
    index!(record)
    return record
end

"""
    SAMRecord(data::Vector{UInt8})

Create a `SAMRecord` object from `data` containing a SAM record.
This function verifies the format and indexes fields for accessors.
"""
function SAMRecord(str::AbstractString)
    return convert(SAMRecord, str)
end

function Base.convert(::Type{SAMRecord}, str::AbstractString)
    return SAMRecord(convert(Vector{UInt8}, str))
end

function initialize!(record::SAMRecord)
    record.filled = 1:0
    record.qname = 1:0
    record.flag = 1:0
    record.rname = 1:0
    record.pos = 1:0
    record.mapq = 1:0
    record.cigar = 1:0
    record.rnext = 1:0
    record.pnext = 1:0
    record.tlen = 1:0
    record.seq = 1:0
    record.qual = 1:0
    empty!(record.fields)
    return record
end

function isfilled(record::SAMRecord)
    return !isempty(record.filled)
end

function datarange(record::SAMRecord)
    return record.filled
end

function checkfilled(record::SAMRecord)
    if !isfilled(record)
        throw(ArgumentError("unfilled SAM record"))
    end
end

function Base.show(io::IO, record::SAMRecord)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "    sequence name: ", seqname(record))
        println(io, "             flag: ", flag(record))
        println(io, "        reference: ", refname(record))
        println(io, "         position: ", leftposition(record))
        println(io, "  mapping quality: ", mappingquality(record))
        println(io, "            CIGAR: ", cigar(record))
        println(io, "   next reference: ", nextrefname(record))
        println(io, "    next position: ", nextleftposition(record))
        println(io, "  template length: ", templatelength(record))
        println(io, "         sequence: ", sequence(String, record))
        println(io, "   base qualities: ", qualities(String, record))
          print(io, "  optional fields:")
        for field in record.fields
            print(io, ' ', String(record.data[field]))
        end
    else
        print(io, " <not filled>")
    end
end

function Base.print(io::IO, record::SAMRecord)
    write(io, record)
    return nothing
end

function Base.write(io::IO, record::SAMRecord)
    checkfilled(record)
    r = datarange(record)
    return unsafe_write(io, pointer(record.data, first(r)), length(r))
end

function Base.:(==)(record1::SAMRecord, record2::SAMRecord)
    if isfilled(record1) == isfilled(record2) == true
        r1 = datarange(record1)
        r2 = datarange(record2)
        return length(r1) == length(r2) && memcmp(pointer(record1.data, first(r1)), pointer(record2.data, first(r2)), length(r1)) == 0
    else
        return isfilled(record1) == isfilled(record2) == false
    end
end

# TODO
function Base.isless(rec1::SAMRecord, rec2::SAMRecord)
    # compared by the left-most position of an alignment
    if rec1.name == rec2.name
        return isless(rec1.pos, rec2.pos)
    else
        return isless(rec1.name, rec2.name)
    end
end

function Base.copy(record::SAMRecord)
    return SAMRecord(
        copy(record.data),
        record.filled,
        record.qname,
        record.flag,
        record.rname,
        record.pos,
        record.mapq,
        record.cigar,
        record.rnext,
        record.pnext,
        record.tlen,
        record.seq,
        record.qual,
        copy(record.fields))
end

function ismapped(record::SAMRecord)
    return isfilled(record) && (flag(record) & SAM_FLAG_UNMAP == 0)
end

function Bio.seqname(record::SAMRecord)
    checkfilled(record)
    if ismissing(record, record.qname)
        missingerror(:seqname)
    end
    return String(record.data[record.qname])
end

function Bio.hasseqname(record::SAMRecord)
    return isfilled(record) && !ismissing(record, record.qname)
end

function flag(record::SAMRecord)
    checkfilled(record)
    return unsafe_parse_decimal(UInt16, record.data, record.flag)
end

function hasflag(record::SAMRecord)
    return isfilled(record)
end

function refname(record::SAMRecord)
    checkfilled(record)
    if ismissing(record, record.rname)
        missingerror(:refname)
    end
    return String(record.data[record.rname])
end

function hasrefname(record::SAMRecord)
    return isfilled(record) && !ismissing(record, record.rname)
end

function Bio.leftposition(record::SAMRecord)
    checkfilled(record)
    pos = unsafe_parse_decimal(Int, record.data, record.pos)
    if pos == 0
        missingerror(:leftposition)
    end
    return pos
end

function Bio.hasleftposition(record::SAMRecord)
    return isfilled(record) && (length(record.pos) != 1 || record.data[first(record.pos)] != UInt8('0'))
end

function rightposition(record::SAMRecord)
    return leftposition(record) + alignment_length(record) - 1
end

function mappingquality(record::SAMRecord)
    checkfilled(record)
    qual = unsafe_parse_decimal(UInt8, record.data, record.mapq)
    if qual == 0xff
        missingerror(:mappingquality)
    end
    return qual
end

function hasmappingquality(record::SAMRecord)
    return isfilled(record) && unsafe_parse_decimal(UInt8, record.data, record.mapq) != 0xff
end

function cigar(record::SAMRecord)
    checkfilled(record)
    if ismissing(record, record.cigar)
        missingerror(:cigar)
    end
    return String(record.data[record.cigar])
end

function hascigar(record::SAMRecord)
    return isfilled(record) && !ismissing(record, record.cigar)
end

function nextrefname(record::SAMRecord)
    checkfilled(record)
    if ismissing(record, record.rnext)
        missingerror(:nextrefname)
    end
    return String(record.data[record.rnext])
end

function hasnextrefname(record::SAMRecord)
    return isfilled(record) && !ismissing(record, record.rnext)
end

function nextleftposition(record::SAMRecord)
    checkfilled(record)
    pos = unsafe_parse_decimal(Int, record.data, record.pnext)
    if pos == 0
        missingerror(:nextleftposition)
    end
    return pos
end

function hasnextleftposition(record::SAMRecord)
    return isfilled(record) && (length(record.pnext) != 1 || record.data[first(record.pnext)] != UInt8('0'))
end

function templatelength(record::SAMRecord)
    checkfilled(record)
    len = unsafe_parse_decimal(Int, record.data, record.tlen)
    if len == 0
        missingerror(:templatelength)
    end
    return len
end

function hastemplatelength(record::SAMRecord)
    return isfilled(record) && (length(record.tlen) != 1 || record.data[first(record.tlen)] != UInt8('0'))
end

function sequence(record::SAMRecord)
    checkfilled(record)
    if ismissing(record, record.seq)
        missingerror(:sequence)
    end
    seqlen = length(record.seq)
    ret = Bio.Seq.DNASequence(seqlen)
    Bio.Seq.encode_copy!(ret, 1, record.data, first(record.seq), seqlen)
    return ret
end

function hassequence(record::SAMRecord)
    return isfilled(record) && !ismissing(record, record.seq)
end

function sequence(::Type{String}, record::SAMRecord)
    checkfilled(record)
    return String(record.data[record.seq])
end

function qualities(record::SAMRecord)
    checkfilled(record)
    if ismissing(record, record.qual)
        missingerror(:qualities)
    end
    qual = record.data[record.qual]
    for i in 1:endof(qual)
        @inbounds qual[i] -= 33
    end
    return qual
end

function hasqualities(record::SAMRecord)
    return isfilled(record) && !ismissing(record, record.qual)
end

function qualities(::Type{String}, record::SAMRecord)
    checkfilled(record)
    return String(record.data[record.qual])
end

function seqlength(record::SAMRecord)
    checkfilled(record)
    if ismissing(record, record.seq)
        missingerror(:seqlength)
    end
    return length(record.seq)
end

function alignment(record::SAMRecord)
    if ismapped(record)
        return Alignment(cigar(record), 1, leftposition(record))
    else
        return Alignment(AlignmentAnchor[])
    end
end

function Base.haskey(record::SAMRecord, tag::AbstractString)
    return findtag(record, tag) > 0
end

function Base.keys(record::SAMRecord)
    checkfilled(record)
    return [String(record.data[first(f):first(f)+1]) for f in record.fields]
end

function Base.values(record::SAMRecord)
    return [record[k] for k in keys(record)]
end

function Base.getindex(record::SAMRecord, tag::AbstractString)
    i = findtag(record, tag)
    if i == 0
        throw(KeyError(tag))
    end
    field = record.fields[i]
    # data type
    typ = record.data[first(field)+3]
    lo = first(field) + 5
    if i == endof(record.fields)
        hi = last(field)
    else
        hi = first(record.fields[i+1]) - 2
    end
    if typ == UInt8('A')
        @assert lo == hi
        return Char(record.data[lo])
    elseif typ == UInt8('i')
        return unsafe_parse_decimal(Int, record.data, lo:hi)
    elseif typ == UInt8('f')
        # TODO: Call a C function directly for speed?
        return parse(Float32, SubString(record.data[lo:hi]))
    elseif typ == UInt8('Z')
        return String(record.data[lo:hi])
    elseif typ == UInt8('H')
        return parse_hexarray(record.data, lo:hi)
    elseif typ == UInt8('B')
        return parse_typedarray(record.data, lo:hi)
    else
        throw(ArgumentError("type code '$(Char(typ))' is not defined"))
    end
end

function findtag(record::SAMRecord, tag::AbstractString)
    checkfilled(record)
    if sizeof(tag) != 2
        return 0
    end
    t1, t2 = UInt8(tag[1]), UInt8(tag[2])
    for (i, field) in enumerate(record.fields)
        p = first(field)
        if record.data[p] == t1 && record.data[p+1] == t2
            return i
        end
    end
    return 0
end

function optional_fields(record::SAMRecord)
    return Dict(k => record[k] for k in keys(record))
end

# Return the length of alignment.
function alignment_length(record::SAMRecord)
    if length(record.cigar) == 1 && record.data[first(record.cigar)] == UInt8('*')
        return 0
    end
    ret = 0
    len = 0  # operation length
    for i in record.cigar
        c = record.data[i]
        if c ∈ UInt8('0'):UInt8('9')
            len = len * 10 + (c - UInt8('0'))
        else
            op = convert(Operation, Char(c))
            if ismatchop(op) || isdeleteop(op)
                ret += len
                len = 0
            end
        end
    end
    return ret
end

# r"[0-9]+" must match `data[range]`.
function unsafe_parse_decimal{T<:Unsigned}(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int})
    x = zero(T)
    @inbounds for i in range
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return x
end

# r"[-+]?[0-9]+" must match `data[range]`.
function unsafe_parse_decimal{T<:Signed}(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int})
    lo = first(range)
    if data[lo] == UInt8('-')
        sign = T(-1)
        lo += 1
    elseif data[lo] == UInt8('+')
        sign = T(+1)
        lo += 1
    else
        sign = T(+1)
    end
    x = zero(T)
    @inbounds for i in lo:last(range)
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return sign * x
end

function parse_hexarray(data::Vector{UInt8}, range::UnitRange{Int})
    @assert iseven(length(range))
    ret = Vector{UInt8}(length(range) >> 1)
    byte2hex(b) = b ∈ 0x30:0x39 ? (b - 0x30) : b ∈ 0x41:0x46 ? (b - 0x41 + 0x0A) : error("not in [0-9A-F]")
    j = 1
    for i in first(range):2:last(range)-1
        ret[j] = (byte2hex(data[range[i]]) << 4) | byte2hex(data[range[i+1]])
        j += 1
    end
    return ret
end

function parse_typedarray(data::Vector{UInt8}, range::UnitRange{Int})
    # format: [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+
    t = data[first(range)]
    xs = split(String(data[first(range)+2:last(range)]))
    if t == UInt8('c')
        return [parse(Int8, x) for x in xs]
    elseif t == UInt8('C')
        return [parse(UInt8, x) for x in xs]
    elseif t == UInt8('s')
        return [parse(Int16, x) for x in xs]
    elseif t == UInt8('S')
        return [parse(UInt16, x) for x in xs]
    elseif t == UInt8('i')
        return [parse(Int32, x) for x in xs]
    elseif t == UInt8('I')
        return [parse(UInt32, x) for x in xs]
    elseif t == UInt8('f')
        return [parse(Float32, x) for x in xs]
    else
        throw(ArgumentError("type code '$(Char(t))' is not defined"))
    end
end

function ismissing(record::SAMRecord, range::UnitRange{Int})
    return length(range) == 1 && record.data[first(range)] == UInt8('*')
end

function memcmp(p1::Ptr, p2::Ptr, n::Integer)
    return ccall(:memcmp, Cint, (Ptr{Void}, Ptr{Void}, Csize_t), p1, p2, n)
end


text:
import Base.-

-(A::Assoc,B::Assoc) = minus(A::Assoc,B::Assoc)
function minus(A::Assoc,B::Assoc)
    #Check A,B, if string => Logical
    At = A
    Bt = B
    if(A.val != [1.0])
        At = logical(A)
    end
    if(B.val !=[1.0])
        Bt = logical(B)
    end
    ## A*B operation
    ABrow = sortedunion(At.row,Bt.row)
    ABcol = sortedunion(At.col,Bt.col)
    if (size(ABrow,1) == 0 || size(ABcol,1) == 0)
        return Assoc([1],[1],0,(+))
    end

    ABA =  spzeros(length(ABrow),length(ABcol))
    Arow = searchsortedmapping(At.row,ABrow)
    Acol = searchsortedmapping(At.col,ABcol)
    Brow = searchsortedmapping(Bt.row,ABrow)
    Bcol = searchsortedmapping(Bt.col,ABcol)
    ABA[Arow,Acol] += At.A
    ABA[Brow,Bcol] -= Bt.A
    AB = Assoc(ABrow,ABcol, Array{Union{AbstractString,Number}}([1.0]), ABA)

    return condense(AB) #Incase of negation
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
module SplittingPool

# scan into a sorted list of free buffers, splitting buffers along the way

import ..@pool_timeit, ..actual_alloc, ..actual_free

using DataStructures

using CUDAdrv

## tunables

# pool boundaries
const SMALL_CUTOFF   = 2^20    # 1 MiB
const LARGE_CUTOFF   = 2^25    # 32 MiB

# memory size rounding
const SMALL_ROUNDOFF = 2^9     # 512 bytes
const LARGE_ROUNDOFF = 2^17    # 128 KiB
const HUGE_ROUNDOFF  = 2^20    # 1 MiB

# maximum overhead (unused space) when selecting a buffer
# small and large buffers with large overhead will be split, but huge buffers never are
const SMALL_OVERHEAD = typemax(Int)
const LARGE_OVERHEAD = typemax(Int)
const HUGE_OVERHEAD  = 0

## block of memory

using Printf

@enum BlockState begin
    AVAILABLE
    ALLOCATED
    FREED
end

const block_id = Ref(UInt(0))

# TODO: it would be nice if this could be immutable, since that's what OrderedSet requires
mutable struct Block
    buf::Mem.Buffer     # base allocation
    sz::Integer         # size into it
    off::Integer        # offset into it

    state::BlockState
    prev::Union{Nothing,Block}
    next::Union{Nothing,Block}

    id::UInt

    Block(buf, sz=sizeof(buf), off=0, state=AVAILABLE, prev=nothing, next=nothing) =
        new(buf, sz, off, state, prev, next, block_id[]+=1)
end

Base.sizeof(block::Block) = block.sz
Base.pointer(block::Block) = pointer(block.buf) + block.off

Base.convert(::Type{Mem.Buffer}, block::Block) = similar(block.buf, pointer(block), sizeof(block))

iswhole(block::Block) = block.prev === nothing && block.next === nothing

## block utilities

function Base.show(io::IO, block::Block)
    fields = [@sprintf("#%d", block.id)]
    push!(fields, @sprintf("%s at %p", Base.format_bytes(sizeof(block)), pointer(block)))
    push!(fields, "$(block.state)")
    block.prev !== nothing && push!(fields, @sprintf("prev=Block(#%d)", block.prev.id))
    block.next !== nothing && push!(fields, @sprintf("next=Block(#%d)", block.next.id))

    print(io, "Block(", join(fields, ", "), ")")
end

# split a block at size `sz`, returning the newly created block
function split!(block, sz)
    @assert sz < block.sz
    split = Block(block.buf, sizeof(block) - sz, block.off + sz)
    block.sz = sz

    # update links
    split.prev = block
    split.next = block.next
    if block.next !== nothing
        block.next.prev = split
    end
    block.next = split

    return split
end

# merge a sequence of blocks `blocks`
function merge!(head, blocks...)
    for block in blocks
        @assert head.next === block
        head.sz += block.sz

        # update links
        tail = block.next
        head.next = tail
        if tail !== nothing
            tail.prev = head
        end
    end

    return head
end

function actual_free(block::Block)
    @assert iswhole(block) "Cannot free a split block"
    if block.state == ALLOCATED
        error("Cannot free an allocated block")
    elseif block.state == FREED
        error("Double-free")
    else
        actual_free(block.buf)
        block.state = FREED
    end
    return
end

## pooling

using Base.Threads
const pool_lock = SpinLock()   # protect against deletion from freelists

function scan!(blocks, sz, max_overhead=typemax(Int))
    max_sz = max(sz + max_overhead, max_overhead)   # protect against overflow
    lock(pool_lock) do
        for block in blocks
            if sz <= sizeof(block) <= max_sz
                delete!(blocks, block)
                return block
            end
        end
        return nothing
    end
end

function incremental_compact!(blocks)
    # we mutate the list of blocks, so take a copy
    blocks = Set(blocks)

    compacted = 0
    lock(pool_lock) do
        while !isempty(blocks)
            block = pop!(blocks)
            szclass = size_class(sizeof(block))
            available = (available_small, available_large, available_huge)[szclass]

            # get the first unallocated node in a chain
            head = block
            while head.prev !== nothing && head.prev.state == AVAILABLE
                head = head.prev
            end

            # construct a chain of unallocated blocks
            chain = [head]
            let block = head
                while block.next !== nothing && block.next.state == AVAILABLE
                    block = block.next
                    @assert block in available
                    push!(chain, block)
                end
            end

            # compact the chain into a single block
            if length(chain) > 1
                for block in chain
                    delete!(available, block)
                    delete!(blocks, block)
                end
                block = merge!(chain...)
                @assert !in(block, available) "Collision in the available memory pool"
                push!(available, block)
                compacted += length(chain) - 1
            end
        end
    end
    return compacted
end

# TODO: partial reclaim on ordered list?
function reclaim!(blocks, sz)
    freed = 0
    candidates = []
    lock(pool_lock) do
        # mark non-split blocks
        for block in blocks
            if iswhole(block)
                push!(candidates, block)
            end
        end

        # free them
        for block in candidates
            delete!(blocks, block)
            freed += sizeof(block)
            actual_free(block)
        end
    end

    return freed
end

function repopulate!(blocks)
    lock(pool_lock) do
        for block in blocks
            szclass = size_class(sizeof(block))
            available = (available_small, available_large, available_huge)[szclass]
            @assert !in(block, available) "Collision in the available memory pool"
            push!(available, block)
        end
    end
end

const SMALL = 1
const LARGE = 2
const HUGE  = 3

# sorted containers need unique keys, which the size of a block isn't.
# mix in the block address to keep the key sortable, but unique.
# the size is shifted 24 bits, and as many identifier bits are
# mixed in, supporting 16777216 unique allocations of up to 1 TiB.
unique_sizeof(block::Block) = (UInt64(sizeof(block))<<24) | (UInt64(block.id) & (2<<24-1))
const UniqueIncreasingSize = Base.By(unique_sizeof)

const available_small = SortedSet{Block}(UniqueIncreasingSize)
const available_large = SortedSet{Block}(UniqueIncreasingSize)
const available_huge  = SortedSet{Block}(UniqueIncreasingSize)
const allocated = Dict{Mem.Buffer,Block}()
const freed = Vector{Block}()

function size_class(sz)
    if sz <= SMALL_CUTOFF
        SMALL
    elseif SMALL_CUTOFF < sz <= LARGE_CUTOFF
        LARGE
    else
        HUGE
    end
end

function pool_alloc(sz)
    szclass = size_class(sz)

    # round of the block size
    req_sz = sz
    roundoff = (SMALL_ROUNDOFF, LARGE_ROUNDOFF, HUGE_ROUNDOFF)[szclass]
    sz = cld(sz, roundoff) * roundoff
    szclass = size_class(sz)

    # select a pool
    available = (available_small, available_large, available_huge)[szclass]

    # determine the maximum scan overhead
    max_overhead = (SMALL_OVERHEAD, LARGE_OVERHEAD, HUGE_OVERHEAD)[szclass]

    block = nothing
    for phase in 1:3
        if phase == 2
           @pool_timeit "$phase.0 gc(false)" GC.gc(false)
        elseif phase == 3
           @pool_timeit "$phase.0 gc(true)" GC.gc(true)
        end

        if !isempty(freed)
            @pool_timeit "$phase.1 repopulate + compact" begin
                # `freed` may be modified concurrently, so take a copy
                pending = copy(freed)
                empty!(freed)

                repopulate!(pending)
                incremental_compact!(pending)
            end
        end

        @pool_timeit "$phase.2 scan" begin
            block = scan!(available, sz, max_overhead)
        end
        block === nothing || break

        @pool_timeit "$phase.3 alloc" begin
            buf = actual_alloc(sz)
            block = buf === nothing ? nothing : Block(buf)
        end
        block === nothing || break

        @pool_timeit "$phase.4 reclaim + alloc" begin
            reclaim!(available_small, sz)
            reclaim!(available_large, sz)
            reclaim!(available_huge, sz)
            buf = actual_alloc(sz)
            block = buf === nothing ? nothing : Block(buf)
        end
        block === nothing || break
    end

    if block !== nothing
        @assert size_class(sizeof(block)) == szclass "Cannot satisfy an allocation in pool $szclass with a buffer in pool $(size_class(sizeof(block)))"

        # split the block if that would yield one from the same pool
        # (i.e. don't split off chunks that would be way smaller)
        # TODO: avoid creating unaligned blocks here (doesn't currently seem to happen
        #       because of the roundoff, but we should take care anyway)
        remainder = sizeof(block) - sz
        if szclass != HUGE && remainder > 0 && size_class(remainder) == szclass
            split = split!(block, sz)
            push!(available, split)
        end

        block.state = ALLOCATED
    end

    return block
end

function pool_free(block)
    block.state = AVAILABLE
    push!(freed, block)
end

## interface

init() = return

function deinit()
    @assert isempty(allocated) "Cannot deinitialize memory pool with outstanding allocations"

    repopulate!(pending)
    incremental_compact!(freed)
    empty!(freed)

    for available in (available_small, available_large, available_huge)
        while !isempty(available)
            block = pop!(available)
            actual_free(block)
        end
    end

    return
end

function alloc(sz)
    block = pool_alloc(sz)
    if block !== nothing
        buf = convert(Mem.Buffer, block)
        @assert !haskey(allocated, buf)
        allocated[buf] = block
        return buf
    else
        return nothing
    end
end

function free(buf)
    block = allocated[buf]
    delete!(allocated, buf)
    pool_free(block)
    return
end

used_memory() = mapreduce(sizeof, +, values(allocated); init=0)

cached_memory() = mapreduce(sizeof, +, union(available_small, available_large, available_huge); init=0)

function dump()
    println("Allocated buffers: $((Base.format_bytes(used_memory())))")
    for block in sort(collect(values(allocated)); by=sizeof)
        println(" - ", block)
    end

    println("Available, but fragmented buffers: $((Base.format_bytes(cached_memory())))")
    for block in available_small
        println(" - small ", block)
    end
    for block in available_large
        println(" - large ", block)
    end
    for block in available_huge
        println(" - huge ", block)
    end
end

end


text:
function modeljs(model, x, result::WebIO.Observable)
    code, weights = FluxJS.compile(model, x)
    code = join(split(code, "\n")[1:end-2],"\n") # remove fetch statement
    code = WebIO.JSString(code)

    JSExpr.@js (() -> begin
        $(code)
        JSExpr.@var weights = $weights.map(e -> begin
            tf.tensor(e)
        end)
        JSExpr.@var x = tf.tensor($x)
        model.setWeights(weights)
        window.model = model
        $result[] = model(x).dataSync()
    end)
end

function Base.take!(o::Channel,timeout::Int)
    taken = false
    @async begin
        sleep(timeout)
        !taken && put!(o, nothing)
    end
    s = take!(o)
    taken = true
    return s
end

function compare(res, out)
    out == nothing && return false
    resjs = Array{eltype(res),1}(undef, length(out))
    for i in keys(out)
        resjs[Main.parse(Int, i) + 1] = out[i]
    end
    return all(x -> x, abs.(res .- resjs) .< 10.0^(-5))
end

function testjs(w, model, x)
    s = Scope()
    r = Observable(s, "result", Dict())
    output = Channel{Any}(1)
    mjs = modeljs(model, x, r)
    onimport(s, mjs)
    on(r) do o
        put!(output, o)
    end
    Blink.body!(w, s)
    res = [data(model(x))...]
    @test compare(res, take!(output, 5))
end

loadseq(w) = nothing

function loadseq_(w, config, cb)
    s = Scope(imports=config["files"])
    obs = Observable(s, "loaded", false)
    onimport(s, JSExpr.@js ()->begin
        JSExpr.@var names = $(config["names"])
        JSExpr.@var modules = arguments
        names.forEach((e, i) -> begin
            window[e] = modules[i]
        end)
        $obs[] = true
    end)
    Blink.body!(w, s)

    b = Channel{Any}(1)
    on(obs) do x
        put!(b, 1)
    end
    s = take!(b, 3) # timeout of 3 secs
    s == nothing && throw("files not loading")
    cb()
end

function loadseq(w, config, args...)
    loadseq_(w, config, () -> loadseq(w, args...))
end

function setupWindow()
    w = Blink.Window(Dict(:show => false))
    libs = [
        "https://unpkg.com/bson@2.0.8/browser_build/bson.js",
        "//cdnjs.cloudflare.com/ajax/libs/tensorflow/0.11.7/tf.js"]
    fluxjs = [normpath("$(@__DIR__)/../lib/flux.js")]
    loadseq(w, Dict("files"=>libs, "names"=>["BSON", "tf"]), Dict("files"=>fluxjs, "names"=>["fluxjs"]))
    return w
end


text:
module Gmsh

const deps_path = abspath(joinpath(dirname(@__FILE__), "..", "deps"))

if !isfile(joinpath(deps_path, "deps.jl"))
    error("Gmsh not installed properly, run Pkg.build(\"Gmsh\"), restart Julia and try again")
end
include(joinpath(deps_path, "deps.jl"))
include(joinpath(deps_path, "usr", "lib", "gmsh.jl"))
global const api = gmsh

function __init__()
    # Always check your dependencies that live in `deps.jl`
    check_deps()
end

end # module


text:
function readinput(filename)
    map(group -> filter(!isempty, split(group, "")), readlines(filename))
end

function countadjacentoccupied_p1(values, dd, row, col, m, n)
    occupied = 0
    for (x, y) in dd
        r, c = row + x, col + y
        if !(r in 1:m && c in 1:n)
            continue
        end
        occupied += values[r][c] == "#"
    end
    return occupied
end

function countadjacentoccupied_p2(values, dd, row, col, m, n)
    occupied = 0
    for (x, y) in dd
        r, c = row + x, col + y
        while (valid = (r in 1:m && c in 1:n)) && values[r][c] == "."
            r, c = r + x, c + y
        end
        occupied += (valid && values[r][c] == "#")
    end
    return occupied
end

function countoccupied(values, m, n)
    return sum(values[r][c] == "#" for r in 1:m for c in 1:n)
end

function shuffleseats(original, countadjacentoccupied, threshold)
    changed = deepcopy(original)
    for i in 1:m
        for j in 1:n
            if original[i][j] == "L" && countadjacentoccupied(original, dd, i, j, m, n) == 0
                changed[i][j] = "#"
            elseif original[i][j] == "#" && countadjacentoccupied(original, dd, i, j, m, n) >= threshold
                changed[i][j] = "L"
            end
        end
    end
    return changed
end

function convergeseats(seats, countadjacentoccupied, threshold)
    while true
        next = shuffleseats(seats, countadjacentoccupied, threshold)
        if next == seats
            return next
        end
        seats = next
    end
end

dd = [(1, 0),(-1, 0),(0, 1),(0, -1),(1, 1),(1, -1),(-1, 1),(-1, -1)]

seats = readinput("d11.in")

m = length(seats)
n = length(seats[1])

println("P1: ", countoccupied(convergeseats(seats, countadjacentoccupied_p1, 4), m, n))
println("P2: ", countoccupied(convergeseats(seats, countadjacentoccupied_p2, 5), m, n))

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
# # include this first to help with crashing??
# try
#   @eval using Gtk
# catch err
#   warn("Gtk not loaded. err: $err")
# end

using VisualRegressionTests
using ExamplePlots

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

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

using Plots, FactCheck
# import Images, ImageMagick

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

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

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

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

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

#   win = Gtk.GtkWindowLeaf("Should we make this the new reference image?")
#   push!(win, Gtk.GtkFrameLeaf(imgbox))

#   showall(win)

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

#   destroy(win)
# end

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

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

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

  # ensure consistent results
  srand(1234)

  # test function
  func = (fn, idx) -> begin
    map(eval, example.exprs)
    png(fn)
  end

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

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

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

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

  # the test
  vtest = VisualTest(func, reffn, idx)
  test_images(vtest, popup=popup, sigma=sigma, eps=eps)

  # try

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

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

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

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

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

  #   if isinteractive()

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

  #   else

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

  # end
end

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


text:
using Test

using Qecsim.GenericModels
using Qecsim.Model
using Qecsim.BasicModels:FiveQubitCode
using Qecsim.PauliTools:to_pauli
using Random:MersenneTwister

# utility test function
function test_simpleerrormodel(
    error_model::SimpleErrorModel,
    p::Real,
    expected_label::AbstractString,
    expected_probability_distribution::NTuple{4,Real},
    excluded_paulis::AbstractString=""
)
    # label
    @test label(error_model) == expected_label
    # probability_distribution
    @test probability_distribution(error_model, p) == expected_probability_distribution
    # generate
    code = FiveQubitCode()
    error = generate(error_model, code, p)
    error_pauli = to_pauli(error)
    @test length(error_pauli) == nkd(code)[1]
    for pauli in excluded_paulis
        @test count(==(pauli), error_pauli) == 0
    end
    # generate (seeded)
    error1 = generate(error_model, code, p, MersenneTwister(13))
    error2 = generate(error_model, code, p, MersenneTwister(13))
    @test error1 == error2
end

@testset "BitFlipErrorModel" begin
    error_model = BitFlipErrorModel()
    for p in [0.1, 1//10]
        test_simpleerrormodel(error_model, p, "Bit-flip", (1 - p, p, 0, 0), "YZ")
    end
end

@testset "BitPhaseFlipErrorModel" begin
    error_model = BitPhaseFlipErrorModel()
    for p in [0.1, 1//10]
        test_simpleerrormodel(error_model, p, "Bit-phase-flip", (1 - p, 0, p, 0), "XZ")
    end
end

@testset "DepolarizingErrorModel" begin
    error_model = DepolarizingErrorModel()
    for p in [0.1, 1//10]
        test_simpleerrormodel(error_model, p, "Depolarizing", (1 - p, p / 3, p / 3, p / 3))
    end
end

@testset "PhaseFlipErrorModel" begin
    error_model = PhaseFlipErrorModel()
    for p in [0.1, 1//10]
        test_simpleerrormodel(error_model, p, "Phase-flip", (1 - p, 0, 0, p), "XY")
    end
end


text:
abstract type AbstractObject{T <: SSDFloat} end

@inline function in(pt::AbstractCoordinatePoint{T}, c::AbstractObject{T})::Bool where {T}
    return in(pt, c.geometry)
end

# function in(p::AbstractCoordinatePoint{T}, c::AbstractObject{T}, rs::Vector{T}) where T
#     rv = false
#     for g in c.geometry_positive
#         if typeof(g) in []#[ConeMantle{T}]
#             in(p, g, rs) ? rv = true : nothing
#         else
#             (p in g) ? rv = true : nothing
#         end
#     end
#     return rv
# end

function in(p::AbstractCoordinatePoint{T}, v::AbstractVector{<:AbstractObject{T}}) where T
    rv = false
    for object in v
        if p in object
            rv = true
        end
    end
    return rv
end

function find_closest_gridpoint(point::CylindricalPoint{T}, grid::CylindricalGrid{T}) where T
    return Int[searchsortednearest( grid[:r].ticks, point.r), searchsortednearest(grid[:φ].ticks, point.φ), searchsortednearest(grid[:z].ticks, point.z)]
end
function find_closest_gridpoint(point::CartesianPoint{T}, grid::CylindricalGrid{T}) where T
    find_closest_gridpoint(CylindricalPoint(point),grid)
end

function find_closest_gridpoint(point::CartesianPoint{T}, grid::CartesianGrid{T}) where T
    return Int[searchsortednearest( grid[:x].ticks, point.x), searchsortednearest(grid[:y].ticks, point.y), searchsortednearest(grid[:z].ticks, point.z)]
end
function find_closest_gridpoint(point::CylindricalPoint{T}, grid::CartesianGrid{T}) where T
    find_closest_gridpoint(CartesianPoint(point),grid)
end

function paint_object(object::AbstractObject, grid::CylindricalGrid{T}) where T
    stepsize::Vector{T}= [minimum(diff(grid[:r].ticks)), IntervalSets.width(grid[:φ].interval) == 0.0 ? 0.05236 : minimum(diff(grid[:φ].ticks)), minimum(diff(grid[:z].ticks))]
    samples = filter(x-> x in object.geometry, vcat([sample(g, stepsize) for g in object.geometry_positive]...))
    object_gridpoints = unique!([find_closest_gridpoint(sample_point,grid) for sample_point in samples])
    return object_gridpoints
end
function paint_object(object::AbstractObject{T}, grid::CylindricalGrid{T}, φ::T) where T
    closest_φ_idx=searchsortednearest(grid[:φ].ticks, φ)
    stepsize::Vector{T}= [minimum(diff(grid[:r].ticks)), IntervalSets.width(grid[:φ].interval) == 0.0 ? 0.05236 : minimum(diff(grid[:φ].ticks)), minimum(diff(grid[:z].ticks))]
    samples = filter(x-> x in object.geometry, vcat([sample(g, stepsize) for g in object.geometry_positive]...))
    object_gridpoints = unique!([find_closest_gridpoint(sample_point,grid) for sample_point in samples])
    return filter(x -> x[2]==closest_φ_idx, object_gridpoints)
end

function paint_object(object::AbstractObject{T}, grid::CartesianGrid{T}) where T
    stepsize::Vector{T}= [minimum(diff(grid[:x].ticks)), minimum(diff(grid[:y].ticks)), minimum(diff(grid[:z].ticks))]
    samples = filter(x-> x in object.geometry, vcat([sample(g, stepsize) for g in object.geometry_positive]...))
    object_gridpoints = unique!([find_closest_gridpoint(sample_point,grid) for sample_point in samples])
    return object_gridpoints
end


text:
@testset "'SampleSpace' ........................ " begin

  ss = SampleSpace(collect(5:10), 25)

  ss = SampleSpace(5:10, 25)

  @test interimsamplesizerange(ss) == collect(5:10)

  @test maxsamplesize(ss) == 25

  @test all(possible.(5:10, ss))

  @test all(.!possible.([-3, 0, 4, 11, 20], ss))

  @test isgroupsequential(ss) == false

  ss = SampleSpace(5:10, 25, n2min = 10)
  ss = SampleSpace(5:10, 25, maxnfact = 5)
  ss = SampleSpace(5:10, 25, nmincont = 10)
  ss = SampleSpace(5:10, 25, maxvariables = 10000)

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
using StaticArrays, SparseArrays
const UP, DOWN = 1, 2
const Float = Float64
# const Operator = SparseMatrixCSC{Float,Int64}
const Operator = Matrix{Float}


text:
export LogNormal, momentMatching, logMomentMatching

"""
Description:

    A log-normal node with location-scale parameterization:

    f(out,m,s) = logN(out|m, s) = 1/out (2π s)^{-1/2} exp(-1/(2s) (log(out) - m)^2))

Interfaces:

    1. out
    2. m (location)
    3. s (squared scale)

Construction:

    LogNormal(out, m, s, id=:some_id)
"""
mutable struct LogNormal <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function LogNormal(out, m, s; id=generateId(LogNormal))
        @ensureVariables(out, m, s)
        self = new(id, Array{Interface}(undef, 3), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:m] = self.interfaces[2] = associate!(Interface(self), m)
        self.i[:s] = self.interfaces[3] = associate!(Interface(self), s)

        return self
    end
end

slug(::Type{LogNormal}) = "log𝒩"

format(dist::ProbabilityDistribution{Univariate, LogNormal}) = "$(slug(LogNormal))(m=$(format(dist.params[:m])), s=$(format(dist.params[:s])))"

ProbabilityDistribution(::Type{Univariate}, ::Type{LogNormal}; m::Float64=1.0, s::Float64=1.0) = ProbabilityDistribution{Univariate, LogNormal}(Dict(:m=>m, :s=>s))
ProbabilityDistribution(::Type{LogNormal}; m::Float64=1.0, s::Float64=1.0) = ProbabilityDistribution{Univariate, LogNormal}(Dict(:m=>m, :s=>s))

dims(dist::ProbabilityDistribution{Univariate, LogNormal}) = 1

vague(::Type{LogNormal}) = ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=huge)

unsafeMean(dist::ProbabilityDistribution{Univariate, LogNormal}) = exp(dist.params[:m] + 0.5*dist.params[:s])
unsafeLogMean(dist::ProbabilityDistribution{Univariate, LogNormal}) = dist.params[:m]

unsafeVar(dist::ProbabilityDistribution{Univariate, LogNormal}) = (exp(dist.params[:s]) - 1.0)*exp(2.0*dist.params[:m] + dist.params[:s])
unsafeLogVar(dist::ProbabilityDistribution{Univariate, LogNormal}) = dist.params[:s]

unsafeCov(dist::ProbabilityDistribution{Univariate, LogNormal}) = unsafeVar(dist)
unsafeLogCov(dist::ProbabilityDistribution{Univariate, LogNormal}) = dist.params[:s]

logPdf(dist::ProbabilityDistribution{Univariate, LogNormal},x) = -0.5*(log(2pi)+log(dist.params[:s])) -log(x) -0.5*(log(x)-dist.params[:m])^2/dist.params[:s]
isProper(dist::ProbabilityDistribution{Univariate, LogNormal}) = (dist.params[:s] > 0.0)

function sample(dist::ProbabilityDistribution{Univariate, LogNormal})
    return exp(dist.params[:m]+sqrt(dist.params[:s])*randn())
end

"""
Gamma approximation to the log-normal distribution using Laplace's method
"""
laplace(::Type{Gamma}, dist::ProbabilityDistribution{Univariate, LogNormal}) = ProbabilityDistribution(Univariate, Gamma, a=1/dist.params[:s], b=1/dist.params[:s]*exp(-dist.params[:m]))

@symmetrical function prod!(x::ProbabilityDistribution{Univariate, LogNormal},
                            y::ProbabilityDistribution{Univariate, Gamma},
                            z::ProbabilityDistribution{Univariate, Gamma}=ProbabilityDistribution(Univariate, Gamma, a=1.0, b=1.0))

    x_approx = laplace(Gamma, x)
    z.params[:a] = x_approx.params[:a] + y.params[:a] - 1.0
    z.params[:b] = x_approx.params[:b] + y.params[:b]

    return z
end

@symmetrical function prod!(x::ProbabilityDistribution{Univariate, LogNormal},
                            y::ProbabilityDistribution{Univariate, PointMass},
                            z::ProbabilityDistribution{Univariate, PointMass}=ProbabilityDistribution(Univariate, PointMass, m=0.0))

    (y.params[:m] > 0.0) || error("PointMass location $(y.params[:m]) should be positive")
    z.params[:m] = y.params[:m]

    return z
end

# Entropy functional
function differentialEntropy(dist::ProbabilityDistribution{Univariate, LogNormal})
    0.5*log(dist.params[:s]) +
    dist.params[:m] + 0.5 +
    0.5*log(2*pi)
end

# Average energy functional
function averageEnergy(::Type{LogNormal}, marg_out::ProbabilityDistribution{Univariate}, marg_m::ProbabilityDistribution{Univariate}, marg_s::ProbabilityDistribution{Univariate})
    unsafeLogMean(marg_out) +
    0.5*log(2*pi) +
    0.5*unsafeLogMean(marg_s) +
    0.5*unsafeInverseMean(marg_s)*( unsafeCov(marg_m) + unsafeLogCov(marg_out) + (unsafeLogMean(marg_out) - unsafeMean(marg_m))^2 )
end