text:
module BoxChar

include("./../../constraints/geometric/constants.jl")
using JSON
using .Constants

using JSON

"""
json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox6/30box6.json")

add = sum(json["items_to_stack"][i]["weight"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*100*json["items_to_stack"][i]["depth"]*100*json["items_to_stack"][i]["height"]*100 for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox7/30box7.json")

add = sum(json["items_to_stack"][i]["weight"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*100*json["items_to_stack"][i]["depth"]*100*json["items_to_stack"][i]["height"]*100 for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox8/30box8.json")

add = sum(json["items_to_stack"][i]["weight"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*100*json["items_to_stack"][i]["depth"]*100*json["items_to_stack"][i]["height"]*100 for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox9/30box9.json")

add = sum(json["items_to_stack"][i]["weight"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*100*json["items_to_stack"][i]["depth"]*100*json["items_to_stack"][i]["height"]*100 for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox10/30box10.json")

add = sum(json["items_to_stack"][i]["weight"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*100*json["items_to_stack"][i]["depth"]*100*json["items_to_stack"][i]["height"]*100 for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

print('\n')
"""

"""
json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox1/30box1.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox2/30box2.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox3/30box3.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox4/30box4.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox5/30box5.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox6/30box6.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox7/30box7.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox8/30box8.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox9/30box9.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox10/30box10.json")

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])

print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox1/7box1.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox2/7box2.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox3/7box3.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox4/7box4.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox5/7box5.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox6/7box6.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox7/7box7.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox8/7box8.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox9/7box9.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments//Nbox10/7box10.json")

add = minimum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = sum(json["items_to_stack"][i]["weight"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*100*json["items_to_stack"][i]["depth"]*100*json["items_to_stack"][i]["height"]*100 for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = sum(json["items_to_stack"][i]["width"]*json["items_to_stack"][i]["depth"]*json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = maximum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["width"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["depth"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

add = minimum(json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"]))
print(add, '\n')

liste = [json["items_to_stack"][i]["height"] for i in 1:length(json["items_to_stack"])]
sort!(liste)
maks = liste[length(liste)]
andmaks = liste[length(liste)-1]
diff = maks - andmaks
diff = round(diff; digits=4)
print(diff, '\n')

wd = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["depth"]) for i in 1:length(json["items_to_stack"])]
wh = [abs(json["items_to_stack"][i]["width"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]
dh = [abs(json["items_to_stack"][i]["depth"] - json["items_to_stack"][i]["height"]) for i in 1:length(json["items_to_stack"])]

add = sum((wd[i] + wh[i] + dh[i])/3 for i in 1:length(json["items_to_stack"]))/length(json["items_to_stack"])
print(add, '\n')

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments/Nbox2/30box2.json")

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments/Nbox3/30box3.json")

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments/Nbox4/30box4.json")

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments/Nbox5/30box5.json")

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments/Nbox6/30box6.json")

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments/Nbox7/30box7.json")

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments/Nbox8/30box8.json")

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments/Nbox9/30box9.json")

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

json = JSON.parsefile("/Users/Toralf/Desktop/Master/Package/Master/src/data/experiments/Nbox10/30box10.json")

d = Dict()
for i in 1:length(json["items_to_stack"])
    tup = (json["items_to_stack"][i]["width"], json["items_to_stack"][i]["depth"], json["items_to_stack"][i]["height"])
    if get(d, tup, -1) == -1
        get!(d, tup, 1)
    else
        val = d[tup]
        #print(get(d, tup, -1), '\n')
        val = val + 1
        d[tup] = val
    end
end
print(sum(values(d))/length(d))
print('\n')

"""

end

ext: jl


text:
#
println("Dk) Test of the CoulombExcitation module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-CoulombExcitation.sum")
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(36.),
                        initialConfigs=[Configuration("1s^2 2s^2 2p^6")],
                        finalConfigs  =[Configuration("1s^2 2s^2 2p^5"), Configuration("1s^2 2s 2p^6") ],
                        process = Coulex(),
                        processSettings=CoulombExcitation.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], [3000., 4000.], false, true,
                                        LineSelection()) )

wb = perform(wa)
setDefaults("print summary: close", "")


text:
function get_output_str(var::Variable)
    output = "{DeepShiba.ShibaObject.Variable}:\n"
    if var.name != ""
        output *= "$(var.name) :\n"
    end

    output *= "data: $(var.data)\n"

    if (var.grad !== nothing)
        output *= "grad: $(var.grad.data)\n"
    end

    if (var.creator !== nothing)
        output *= "creator: $(typeof(var.creator))\n"
    else
        output *= "User-Defined\n"
    end
    return output
end

function Base.show(io::IO, var::Variable)
    print(io, var.data)
end

function Base.show(io::IO, ::MIME"text/plain", var::Variable)
    print(io, get_output_str(var))
end

function get_value_type(nested_var)
    if eltype(nested_var) <: Real
        return eltype(nested_var)
    else
        get_value_type(eltype(nested_var))
    end
end


text:
module Snopt

using SparseArrays

export snopta

const snoptlib = joinpath(dirname(@__FILE__), "../deps/src/libsnopt")

const codes = Dict(
1 => "Finished successfully: optimality conditions satisfied",
2 => "Finished successfully: feasible point found",
3 => "Finished successfully: requested accuracy could not be achieved",
11 => "The problem appears to be infeasible: infeasible linear constraints",
12 => "The problem appears to be infeasible: infeasible linear equalities",
13 => "The problem appears to be infeasible: nonlinear infeasibilities minimized",
14 => "The problem appears to be infeasible: infeasibilities minimized",
15 => "The problem appears to be infeasible: infeasible linear constraints in QP subproblem",
21 => "The problem appears to be unbounded: unbounded objective",
22 => "The problem appears to be unbounded: constraint violation limit reached",
31 => "Resource limit error: iteration limit reached",
32 => "Resource limit error: major iteration limit reached",
33 => "Resource limit error: the superbasics limit is too small",
41 => "Terminated after numerical difficulties: current point cannot be improved",
42 => "Terminated after numerical difficulties: singular basis",
43 => "Terminated after numerical difficulties: cannot satisfy the general constraints",
44 => "Terminated after numerical difficulties: ill-conditioned null-space basis",
51 => "Error in the user-supplied functions: incorrect objective derivatives",
52 => "Error in the user-supplied functions: incorrect constraint derivatives",
61 => "Undefined user-supplied functions: undefined function at the first feasible point",
62 => "Undefined user-supplied functions: undefined function at the initial point",
63 => "Undefined user-supplied functions: unable to proceed into undefined region",
71 => "User requested termination: terminated during function evaluation ",
74 => "User requested termination: terminated from monitor routine",
81 => "Insufficient storage allocated: work arrays must have at least 500 elements",
82 => "Insufficient storage allocated: not enough character storage",
83 => "Insufficient storage allocated: not enough integer storage",
84 => "Insufficient storage allocated: not enough real storage",
91 => "Input arguments out of range: invalid input argument",
92 => "Input arguments out of range: basis file dimensions do not match this problem",
141 => "System error: wrong number of basic variables",
142 => "System error: error in basis package"
)

const PRINTNUM = 18
const SUMNUM = 19

"""
    Names(prob, xnames, fnames)

Convenience to put custom names in output files.
Arguments to variables of same names in snOptA.
Use strings and vectors of strings.
"""
struct Names{TS}
    prob::TS
    xnames::Vector{TS}
    fnames::Vector{TS}
end

"""
    Names()

Default names (uses snopt defaults for xnames, fnames)
"""
Names() = Names("Opt Prob", [""], [""])

# Internal: truncater or pad string to 8 characters
function eightchar(name)
    n = length(name)
    if n >= 8
        return name[1:8]
    else
        return string(name, repeat(" ", 8-n))
    end
end

# Internal: struct of names in fortran format
struct NamesFStyle{TC}
    prob::Vector{TC}
    xnames::Vector{TC}
    fnames::Vector{TC}
end

# Internal: convert a list of names to fortran format
function process_list_of_names(names)
    allnames = ""
    for name in names
        allnames *= eightchar(name)
    end
    return Vector{Cuchar}(allnames)
end

# Internal: convert the names to fortran format
function processnames(names)

    prob = Vector{Cuchar}(eightchar(names.prob))

    xnames = process_list_of_names(names.xnames)
    fnames = process_list_of_names(names.fnames)

    return NamesFStyle(prob, xnames, fnames)
end

"""
    Start(start, ns, xstate, fstate, x, f, xmul, fmul)

A starting point for optimization.
Usually not used directly. Instead use ColdStart or WarmStart.
Arguments correspond to the same names in snOptA.
"""
struct Start{TI, TF}
    start::TI
    ns::TI
    xstate::Vector{TI}
    fstate::Vector{TI}
    x::Vector{TF}
    f::Vector{TF}
    xmul::Vector{TF}
    fmul::Vector{TF}
end

"""
    ColdStart(x0, nf)

A cold start.

# Arguments
- `x0::Vector{Float64}`: starting point
- `nf::Int64`: number of output functions (obj + constraints)
"""
function ColdStart(x0, nf)
    nx = length(x0)
    return Start(Cint(0), Cint(0),
        zeros(Cint, nx), zeros(Cint, nf),
        x0, zeros(nf), zeros(nx), zeros(nf)
        )
end

"""
    WarmStart(ns, xstate, fstate, x, f, xmul, fmul)

A warm start.  Arguments correspond to variables of same names in snOptA.
One of the outputs of snopta is a WarmStart object that can be reused as an input.
"""
WarmStart(ns, xstate, fstate, x, f, xmul, fmul) = Start(
    Cint(2), ns, xstate, fstate, x, f, xmul, fmul)

# internal use: workspace vectors
mutable struct Workspace{TI,TF,TS}
    lencw::TI
    cw::Vector{TS}
    leniw::TI
    iw::Vector{TI}
    lenrw::TI
    rw::Vector{TF}
end

# internal use: initialize workspace arrays from lengths
Workspace(lenc, leni, lenr) = Workspace(
    Cint(lenc), Array{Cuchar}(undef, lenc*8),
    Cint(leni), Array{Cint}(undef, leni),
    Cint(lenr), Array{Float64}(undef, lenr)
)

"""
    Outputs(gstar, iterations, major_iter, run_time, nInf, sInf, warm)

Outputs returned by the snopta function.

# Arguments
- `gstar::Vector{Float}`: constraints evaluated at the optimal point
- `iterations::Int`: total iteration count
- `major_iter::Int`: number of major iterations
- `run_time::Float`: solve time as reported by snopt
- `nInf::Int`: number of infeasibility constraints, see snopta docs
- `sInf::Float`: sum of infeasibility constraints, see snopta docs
- `warm::WarmStart`: a warm start object that can be used in a restart.
"""
struct Outputs{TF,TI,TW}
    gstar::Vector{TF}
    iterations::TI
    major_iter::TI
    run_time::TF
    nInf::TI
    sInf::TF
    warm::TW
end

# wrapper for snInit
function sninit(nx, nf)

    # temporary working arrays
    minlen = 500
    lencw = minlen
    leniw = minlen + 100*(nx + nf)
    lenrw = minlen + 200*(nx + nf)
    w = Workspace(lencw, leniw, lenrw)

    ccall( (:sninit_, snoptlib), Nothing,
        (Ref{Cint}, Ref{Cint}, Ptr{Cuchar}, Ref{Cint}, Ptr{Cint},
        Ref{Cint}, Ptr{Cdouble}, Ref{Cint}),
        PRINTNUM, SUMNUM, w.cw, w.lencw, w.iw,
        w.leniw, w.rw, w.lenrw)

    return w
end

# wrapper for openfiles. not defined with snopt, fortran file supplied in repo (from pyoptsparse)
function openfiles(printfile, sumfile)

    # open files for printing (not part of snopt distribution)
    printerr = Cint[0]
    sumerr = Cint[0]
    ccall( (:openfiles_, snoptlib), Nothing,
        (Ref{Cint}, Ref{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cuchar}, Ptr{Cuchar}),
        PRINTNUM, SUMNUM, printerr, sumerr, printfile, sumfile)

    if printerr[1] != 0
        @warn "failed to open print file"
    end
    if sumerr[1] != 0
        @warn "failed to open summary file"
    end

    return nothing
end

# wrapper for closefiles. not defined with snopt, fortran file supplied in repo (from pyoptsparse)
function closefiles()
    # close output files
    ccall( (:closefiles_, snoptlib), Nothing,
        (Ref{Cint}, Ref{Cint}),
        PRINTNUM, SUMNUM)

    return nothing
end

# wrapper for flushfiles. not defined with snopt, fortran file supplied in repo (from pyoptsparse)
function flushfiles()
    # flush output files to see progress
    ccall( (:flushfiles_, snoptlib), Nothing,
        (Ref{Cint}, Ref{Cint}),
        PRINTNUM, SUMNUM)

    return nothing
end

# wrapper for snSet, snSeti, snSetr
function setoptions(options, work)

    # --- set options ----
    errors = Cint[0]

    for key in keys(options)
        value = options[key]
        buffer = string(key, repeat(" ", 55-length(key)))  # buffer length is 55 so pad with space.

        if length(key) > 55
            @warn "invalid option, too long"
            continue
        end

        errors[1] = 0

        if typeof(value) == String

            value = string(value, repeat(" ", 72-length(value)))

            ccall( (:snset_, snoptlib), Nothing,
                (Ptr{Cuchar}, Ref{Cint}, Ref{Cint}, Ptr{Cint},
                Ptr{Cuchar}, Ref{Cint}, Ptr{Cint}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}),
                value, PRINTNUM, SUMNUM, errors,
                work.cw, work.lencw, work.iw, work.leniw, work.rw, work.lenrw)

        elseif isinteger(value)

            ccall( (:snseti_, snoptlib), Nothing,
                (Ptr{Cuchar}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ptr{Cint},
                Ptr{Cuchar}, Ref{Cint}, Ptr{Cint}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}),
                buffer, value, PRINTNUM, SUMNUM, errors,
                work.cw, work.lencw, work.iw, work.leniw, work.rw, work.lenrw)

        elseif isreal(value)

            ccall( (:snsetr_, snoptlib), Nothing,
                (Ptr{Cuchar}, Ref{Cdouble}, Ref{Cint}, Ref{Cint}, Ptr{Cint},
                Ptr{Cuchar}, Ref{Cint}, Ptr{Cint}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}),
                buffer, value, PRINTNUM, SUMNUM, errors,
                work.cw, work.lencw, work.iw, work.leniw, work.rw, work.lenrw)
        end

        if errors[1] > 0
            @warn errors[1], " errors encountered while setting options"
        end

    end

    return nothing
end

# wrapper for snMemA
function setmemory(INFO, nf, nx, nxname, nfname, neA, neG, work)

    mincw = Cint[0]
    miniw = Cint[0]
    minrw = Cint[0]

    # --- set memory requirements --- #
    ccall( (:snmema_, snoptlib), Nothing,
        (Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint},
        Ref{Cint}, Ref{Cint}, Ref{Cint},
        Ptr{Cuchar}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cdouble}, Ref{Cint}),
        INFO, nf, nx, nxname, nfname, neA, neG,
        mincw, miniw, minrw,
        work.cw, work.lencw, work.iw, work.leniw, work.rw, work.lenrw)

    if INFO[1] != 104
        @warn "error in snmema memory setting: ", INFO[1]
    end

    # --- resize arrays to match memory requirements
    if mincw[1] > work.lencw
        work.lencw = mincw[1]
        resize!(work.cw, work.lencw*8)
    end
    if miniw[1] > work.leniw
        work.leniw = miniw[1]
        resize!(work.iw, work.leniw)
    end
    if minrw[1] > work.lenrw
        work.lenrw = minrw[1]
        resize!(work.rw, work.lenrw)
    end

    memkey = ("Total character workspace", "Total integer   workspace",
        "Total real      workspace")
    memvalue = (work.lencw, work.leniw, work.lenrw)
    errors = Cint[0]
    for (key,value) in zip(memkey, memvalue)
        buffer = string(key, repeat(" ", 55-length(key)))  # buffer length is 55 so pad with space.
        errors[1] = 0
        ccall( (:snseti_, snoptlib), Nothing,
            (Ptr{Cuchar}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ptr{Cint},
            Ptr{Cuchar}, Ref{Cint}, Ptr{Cint}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}),
            buffer, value, PRINTNUM, SUMNUM, errors,
            work.cw, work.lencw, work.iw, work.leniw, work.rw, work.lenrw)
        if errors[1] > 0
            @warn errors[1], " error encountered while lengths in options from memory sizing"
        end
    end

    return nothing
end

"""
    f, fail = example!(g, df, dg, x, deriv)

The expected function signature for user functions.

# Arguments
- `g::Vector{Float64}`: (output) constraint vector, modified in place
- `df::Vector{Float64}`: (output) gradient vector, modified in place
- `dg::Vector{Float64}`: (output) constraint jacobian vector, modified in place.
    dgi/dxj in order corresponding to sparsity pattern provided to snopt
    (preferably column major order if dense)
- `x::Vector{Float64}`: (input) design variables, unmodified
- `deriv::Bool`: (input) if false snopt does not need derivatives that iteration so you can skip their computation.

# Returns
- `f::Float64`: objective value
- `fail::Bool`: true if function fails to compute at this x
"""
function example!(g, df, dg, x)
    return 0.0, false
end

# wrapper for usrfun (augmented with function pass-in)
function usrcallback(func!, status_::Ptr{Cint}, nx::Cint, x_::Ptr{Cdouble},
    needf::Cint, nf::Cint, f_::Ptr{Cdouble}, needG::Cint, ng::Cint,
    G_::Ptr{Cdouble}, cu::Ptr{Cuchar}, lencu::Cint, iu::Ptr{Cint},
    leniu::Cint, ru::Ptr{Cdouble}, lenru::Cint)

    # check if solution finished, no need to calculate more
    status = unsafe_load(status_)
    if status >= 2
        return
    end

    # unpack design variables
    x = unsafe_wrap(Array, x_, nx)

    # set functions
    f = unsafe_wrap(Array, f_, nf)
    G = unsafe_wrap(Array, G_, ng)
    f[1], fail = func!(@view(f[2:end]), @view(G[1:nx]), @view(G[nx+1:end]), x, needG > 0)

    # check if solutions fails
    if fail
        unsafe_store!(status_, -1, 1)
    end

end

# convenience function used internally to simplify passing in arguments associated with A (linear constraints).
# automatically computes iAfun, jAvar, A, lenA, neA.
function parseAmatrix(A)

    if typeof(A) <: SparseMatrixCSC
        # len = nnz(A)
        r, c, values = findnz(A)
        len = length(r)

        rows = zeros(Cint, len)
        rows .= r
        cols = zeros(Cint, len)
        cols .= c

    else  # dense
        if isempty(A)
            len = 0
            rows = Int32[1]
            cols = Int32[1]
            values = [0.0]
        else
            nf, nx = size(A)
            len = nf*nx
            rows = zeros(Cint, nx*nf)
            rows .= [i for i = 1:nf, j = 1:nx][:]
            cols = zeros(Cint, nx*nf)
            cols .= [j for i = 1:nf, j = 1:nx][:]
            values = A[:]
        end
    end

    return len, rows, cols, values
end

# commonly used convenience method to provide regular starting point (cold start)
function snopta(func!, x0::T, lx, ux, lg, ug, rows, cols,
    options=Dict(); A=[], names=Names(), objadd=0.0) where T<:Vector

    start = ColdStart(x0, length(lg)+1)

    return snopta(func!, start, lx, ux, lg, ug, rows, cols, options, A=A, names=names, objadd=objadd)
end

"""
    snopta(func!, x0, lx, ux, lg, ug, rows, cols,
        options=Dict(); A=[], names=Names(), objadd=0.0)

Main function call into snOptA.

# Arguments
- `func!::function`: follows function signature shown in example!
- `x0::Vector{Float64}` or `x0::WarmStart`: starting point
- `lx::Vector{Float64}`: lower bounds on x
- `ux::Vector{Float64}`: upper bounds on x
- `lg::Vector{Float64}`: lower bounds on g
- `ug::Vector{Float64}`: upper bounds on g
- `rows::Vector{Int64}`: sparsity pattern for constraint jacobian.  dg[k] corresponds to rows[k], cols[k]
- `cols::Vector{Int64}`: sparsity pattern for constraint jacobian.  dg[k] corresponds to rows[k], cols[k]
- `options::Dict`: dictionary of options (see Snopt docs)
- `A::Matrix` (if dense) or `SparseMatrixCSC`: linear constraints g += A*x
- `names::Names`: custom names for problem and variables for print file
- `objAdd::Float64`: adds a scalar to objective (see Snopt docs)

# Returns
- `xstar::Vector{Float64}`: optimal x
- `fstar::Vector{Float64}`: corresponding f
- `info::String`: termination message
- `out::Outputs`: various outputs
"""
function snopta(func!, start::Start, lx, ux, lg, ug, rows, cols,
    options=Dict(); A=[], names=Names(), objadd=0.0)

    # --- number of functions ----
    nx = length(start.x)
    ng = length(lg)
    nf = 1 + length(lg)
    lf = [0.0; lg]  # bounds on objective are irrelevant
    uf = [0.0; ug]

    # --- parse names -------
    nxname = length(names.xnames)
    if nxname != 1 && nxname != nx
        @warn "incorrect length for xnames"
        nxname = 1
    end
    nfname = length(names.fnames)
    if nfname != 1 && nfname != nf
        @warn "incorrect length for fnames"
        nfname = 1
    end

    if nxname == nx && nfname == nf
        names = processnames(names) # TODO: prob not type stable if names are added
    end

    # ---- parse linear constraints -------
    lenA, iAfun, jAvar, Aval = parseAmatrix(A)
    neA = lenA

    # ----- setup user function ---------------
    wrapper = function(status_::Ptr{Cint}, n::Cint, x_::Ptr{Cdouble},
        needf::Cint, nF::Cint, f_::Ptr{Cdouble}, needG::Cint, lenG::Cint,
        G_::Ptr{Cdouble}, cu::Ptr{Cuchar}, lencu::Cint, iu::Ptr{Cint},
        leniu::Cint, ru::Ptr{Cdouble}, lenru::Cint)

        usrcallback(func!, status_, n, x_, needf, nF, f_, needG, lenG,
            G_, cu, lencu, iu, leniu, ru, lenru)

        return nothing
    end

    # c wrapper to callback function
    usrfun = @cfunction($wrapper, Cvoid, (Ptr{Cint}, Ref{Cint}, Ptr{Cdouble},
        Ref{Cint}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}, Ref{Cint}, Ptr{Cdouble},
        Ptr{Cuchar}, Ref{Cint}, Ptr{Cint}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}))

    # ---- setup jacobian  sparsity pattern ------
    nsp = length(rows)
    lenG = nx + nsp
    neG = lenG
    iGfun = Array{Cint}(undef, lenG)
    jGvar = Array{Cint}(undef, lenG)

    objrow = 1

    # objective gradient (first row) assumed dense
    iGfun[1:nx] .= 1
    jGvar[1:nx] .= 1:nx

    # constraint jacobian
    for k = 1:nsp
        iGfun[nx + k] = rows[k] + 1  # adding one for objective row
        jGvar[nx + k] = cols[k]
    end

    # --- open files ------
    printfile = "snopt-print.out"
    sumfile = "snopt-summary.out"

    if haskey(options, "Print file")
        printfile = options["Print file"]
    end
    if haskey(options, "Summary file")
        sumfile = options["Summary file"]
    end
    openfiles(printfile, sumfile)

    # ----- initialize -------
    work = sninit(nx, nf)

    # --- set options ----
    setoptions(options, work)

    # ---- set memory requirements ------
    INFO = Cint[0]
    setmemory(INFO, nf, nx, nxname, nfname, neA, neG, work)

    # --- call snopta ----
    mincw = Cint[0]
    miniw = Cint[0]
    minrw = Cint[0]
    nInf = Cint[0]
    sInf = Cdouble[0]
    lencu = 1
    cu = Array{Cuchar}(undef, lencu*8)
    leniu = 1
    iu = Cint[0]
    lenru = 1
    ru = [0.0]
    ns = Cint[start.ns]

    ccall( (:snopta_, snoptlib), Nothing,
        (Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint},
        Ref{Cdouble}, Ref{Cint}, Ptr{Cuchar}, Ptr{Nothing},
        Ptr{Cint}, Ptr{Cint}, Ref{Cint}, Ref{Cint}, Ptr{Cdouble},
        Ptr{Cint}, Ptr{Cint}, Ref{Cint}, Ref{Cint},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cuchar}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cuchar},
        Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble},
        Ptr{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint},
        Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
        Ptr{Cuchar}, Ref{Cint}, Ptr{Cint}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint},
        Ptr{Cuchar}, Ref{Cint}, Ptr{Cint}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}),
        start.start, nf, nx, nxname, nfname,
        objadd, objrow, names.prob, usrfun,
        iAfun, jAvar, lenA, neA, Aval,
        iGfun, jGvar, lenG, neG,
        lx, ux, names.xnames, lf, uf, names.fnames,
        start.x, start.xstate, start.xmul, start.f, start.fstate, start.fmul,
        INFO, mincw, miniw, minrw,
        ns, nInf, sInf,
        cu, lencu, iu, leniu, ru, lenru,
        work.cw, work.lencw, work.iw, work.leniw, work.rw, work.lenrw)

    # close output files
    closefiles()

    # pack outputs
    warm = WarmStart(ns[1], start.xstate, start.fstate, start.x, start.f,
    start.xmul, start.fmul)

    out = Outputs(start.f[2:end], work.iw[421], work.iw[422], work.rw[462],
        nInf[1], sInf[1], warm)

    return start.x, start.f[1], codes[INFO[1]], out
end

end  # end module


text:
#
#  AdaptTimestep.jl
#  CollagenPolymer
#
#  Created by Christopher Revell on 06/07/2020.
#
#

module AdaptTimestep

# Import Julia packages
using LinearAlgebra
using StaticArrays

# Import local program modules

function adaptTimestep(F,W,magsF,magsW,σ)

    magsF .= norm.(F)
    maxF = maximum(magsF)

    magsW .= norm.(W)
    maxW = maximum(magsW)

    Δt = min(σ^2/(4.0*maxW^2),σ/(2.0*maxF))

    return Δt

end

export adaptTimestep

end

max_stars_repo_name: chris-revell/collagen-model-jl

max_issues_repo_name: chris-revell/collagen-model-jl

max_forks_repo_name: chris-revell/collagen-model-jl


text:
### Hastings & Powell 1991
### Model equations

### Define parameters and equations
f(u, a, b) = a*u/(1 + b*u)

function parameterized_model!(du, u, p, t)
    x, y, z = u
    a1, a2, b1, b2, d1, d2 = p
    du[1] = dx = x*(1-x) - f(x, a1, b1)*y
    du[2] = dy = f(x, a1, b1)*y - f(y, a2, b2)*z - d1*y
    du[3] = dz = f(y, a2, b2)*z - d2*z
end


text:
@testset "Foreign Exchange Currency" begin
    @testset "currency_exchange_rate" begin
        data = currency_exchange_rate("BTC", "USD")
        @test typeof(data) === Dict{String,Any}
        @test length(data) === 1
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "fx_intraday" begin
        data = fx_intraday("EUR", "USD")
        @test typeof(data) === Dict{String,Any}
        @test length(data) === 2
        sleep(TEST_SLEEP_TIME + 2*rand())
    end

    for f in (:fx_daily, :fx_weekly, :fx_monthly)[1:MAX_TESTS]
        @eval begin
            testname = string($f)
            @testset "$testname" begin
                data = $f("EUR", "USD")
                @test typeof(data) === Dict{String,Any}
                @test length(data) === 2
            end
        end
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

end


text:
#"""
#SiteType is a parameterized type which allows
#making Index tags into Julia types. Use cases
#include overloading functions such as `op`,
#`siteinds`, and `state` which generate custom
#operators, Index arrays, and IndexVals associated
#with Index objects having a certain tag.
#
#To make a SiteType type, you can use the string
#macro notation: `SiteType"MyTag"`
#
#To make an SiteType value or object, you can use
#the notation: `SiteType("MyTag")`
#"""
@eval struct SiteType{T}
  (f::Type{<:SiteType})() = $(Expr(:new, :f))
end

# Note that the complicated definition of
# SiteType above is a workaround for performance
# issues when creating parameterized types
# in Julia 1.4 and 1.5-beta. Ideally we
# can just use the following in the future:
# struct SiteType{T}
# end

SiteType(s::AbstractString) = SiteType{Tag(s)}()
SiteType(t::Tag) = SiteType{t}()
tag(::SiteType{T}) where {T} = T

macro SiteType_str(s)
  SiteType{Tag(s)}
end

# Keep TagType defined for backwards
# compatibility; will be deprecated later
const TagType = SiteType
macro TagType_str(s)
  TagType{Tag(s)}
end

#---------------------------------------
#
# op system
#
#---------------------------------------

#"""
#OpName is a parameterized type which allows
#making strings into Julia types for the purpose
#of representing operator names.
#The main use of OpName is overloading the
#`ITensors.op!` method which generates operators
#for indices with certain tags such as "S=1/2".
#
#To make a OpName type, you can use the string
#macro notation: `OpName"MyTag"`.
#
#To make an OpName value or object, you can use
#the notation: `OpName("myop")`
#"""
@eval struct OpName{Name}
  (f::Type{<:OpName})() = $(Expr(:new, :f))
end

# Note that the complicated definition of
# OpName above is a workaround for performance
# issues when creating parameterized types
# in Julia 1.4 and 1.5-beta. Ideally we
# can just use the following in the future:
# struct OpName{Name}
# end

OpName(s::AbstractString) = OpName{SmallString(s)}()
OpName(s::SmallString) = OpName{s}()
name(::OpName{N}) where {N} = N

macro OpName_str(s)
  OpName{SmallString(s)}
end

# Default implementations of op and op!
op(::OpName, ::SiteType, ::Index...; kwargs...) = nothing
op(::OpName, ::SiteType, ::SiteType,
   sitetypes_inds::Union{SiteType, Index}...; kwargs...) = nothing
op!(::ITensor, ::OpName, ::SiteType, ::Index...; kwargs...) = nothing
op!(::ITensor, ::OpName, ::SiteType, ::SiteType,
    sitetypes_inds::Union{SiteType, Index}...; kwargs...) = nothing

# Deprecated version, for backwards compatibility
op(::SiteType, ::Index, ::AbstractString; kwargs...) = nothing

function _sitetypes(ts::TagSet)
  Ntags = length(ts)
  return SiteType[SiteType(ts[n]) for n in 1:Ntags]
end

_sitetypes(i::Index) = _sitetypes(tags(i))

"""
    op(opname::String, s::Index; kwargs...)

Return an ITensor corresponding to the operator
named `opname` for the Index `s`. The operator
is constructed by calling an overload of either
the `op` or `op!` methods which take a `SiteType`
argument that corresponds to one of the tags of
the Index `s` and an `OpName"opname"` argument
that corresponds to the input operator name.

Operator names can be combined using the "*"
symbol, for example "S+*S-" or "Sz*Sz*Sz".
The result is an ITensor made by forming each operator
then contracting them together in a way corresponding
to the usual operator product or matrix multiplication.

The `op` system is used by the AutoMPO
system to convert operator names into ITensors,
and can be used directly such as for applying
operators to MPS.

# Example
```julia
s = Index(2,"Site,S=1/2")
Sz = op("Sz",s)
```
"""
function op(name::AbstractString,
            s::Index...;
            kwargs...)

  name = strip(name)

  # TODO: filter out only commons tags
  # if there are multiple indices
  commontags_s = commontags(s...)

  # Interpret operator names joined by *
  # as acting sequentially on the same site
  starpos = findfirst("*", name)
  if !isnothing(starpos)
    op1 = name[1:starpos.start-1]
    op2 = name[starpos.start+1:end]
    return product(op(op1, s...; kwargs...),
                   op(op2, s...; kwargs...))
  end

  common_stypes  = _sitetypes(commontags_s)
  push!(common_stypes,SiteType("Generic"))
  opn = OpName(name)

  #
  # Try calling a function of the form:
  #    op(::OpName, ::SiteType, ::Index; kwargs...)
  #
  for st in common_stypes
    res = op(opn, st, s...; kwargs...)
    if !isnothing(res)
      return res
    end
  end

  # otherwise try calling a function of the form:
  #    op!(::ITensor, ::OpName, ::SiteType, ::Index; kwargs...)
  #
  Op = emptyITensor(prime.(s)..., dag.(s)...)
  for st in common_stypes
    op!(Op, opn, st, s...; kwargs...)
    if !isempty(Op)
      return Op
    end
  end

  if length(s) > 1
    # No overloads for common tags found. It might be a
    # case of making an operator with mixed site types,
    # searching for overloads like:
    #   op(::OpName,
    #      ::SiteType...,
    #      ::Index...;
    #      kwargs...)
    #   op!(::ITensor, ::OpName,
    #       ::SiteType...,
    #       ::Index...;
    #       kwargs...)
    stypes = _sitetypes.(s)

    for st in Iterators.product(stypes...)
      res = op(opn, st..., s...; kwargs...)
      if !isnothing(res)
        return res
      end
    end

    Op = emptyITensor(prime.(s)..., dag.(s)...)
    for st in Iterators.product(stypes...)
      op!(Op, opn, st..., s...; kwargs...)
      if !isempty(Op)
        return Op
      end
    end
    error("Older op interface does not support multiple indices with mixed site types. You may want to overload `op(::OpName, ::SiteType..., ::Index...)` or `op!(::ITensor, ::OpName, ::SiteType..., ::Index...) for the operator \"$name\" and Index tags $(tags.(s)).")
  end

  #
  # otherwise try calling a function of the form:
  #   op(::SiteType, ::Index, ::AbstractString)
  #
  # (Note: this version is for backwards compatibility
  #  after version 0.1.10, and may be eventually
  #  deprecated)
  #
  for st in common_stypes
    res = op(st, s[1], name; kwargs...)
    if !isnothing(res)
      return res
    end
  end

  throw(ArgumentError("Overload of \"op\" or \"op!\" functions not found for operator name \"$name\" and Index tags: $(commontags_s))"))
end

# For backwards compatibility, version of `op`
# taking the arguments in the other order:
op(s::Index,
   opname::AbstractString;
   kwargs...) = op(opname, s; kwargs...)

"""
  op(opname::String,sites::Vector{<:Index},n::Int; kwargs...)

Return an ITensor corresponding to the operator
named `opname` for the n'th Index in the array
`sites`.
"""
op(opname::AbstractString,
   s::Vector{<:Index},
   ns::Vararg{Int, N};
   kwargs...) where {N} =
  op(opname, ntuple(n -> s[ns[n]], Val(N))...; kwargs...)

op(s::Vector{<:Index},
   opname::AbstractString,
   ns::Int...;
   kwargs...) =
  op(opname, s, ns...; kwargs...)

#---------------------------------------
#
# state system
#
#---------------------------------------

@eval struct StateName{Name}
  (f::Type{<:StateName})() = $(Expr(:new, :f))
end

StateName(s::AbstractString) = StateName{SmallString(s)}()
StateName(s::SmallString) = StateName{s}()
name(::StateName{N}) where {N} = N

macro StateName_str(s)
  StateName{SmallString(s)}
end

state(::SiteType, ::StateName) = nothing
state(::SiteType, ::AbstractString) = nothing

function state(s::Index,
               name::AbstractString)::IndexVal
  stypes  = _sitetypes(s)
  sname = StateName(name)

  # Try calling state(::SiteType"Tag",::StateName"Name")
  for st in stypes
    res = state(st,sname)
    !isnothing(res) && return s(res)
  end

  # Try calling state(::SiteType"Tag","Name")
  for st in stypes
    res = state(st,name)
    !isnothing(res) && return s(res)
  end

  throw(ArgumentError("Overload of \"state\" function not found for Index tags $(tags(s))"))
end

state(s::Index,n::Integer) = s[n]

state(sset::Vector{<:Index},j::Integer,st) = state(sset[j],st)

#---------------------------------------
#
# siteind system
#
#---------------------------------------

space(st::SiteType; kwargs...) = throw(MethodError("Overload of \"space\",\"siteind\", or \"siteinds\" functions not found for Index tag: $(tag(st))"))

function siteind(st::SiteType; addtags="", kwargs...)
  sp = space(st;kwargs...)
  return Index(sp,"Site,$(tag(st)),$addtags")
end

siteind(st::SiteType, n; kwargs...) = addtags(siteind(st; kwargs...),"n=$n")

siteind(tag::String; kwargs...) = siteind(SiteType(tag);kwargs...)

siteind(tag::String,n; kwargs...) = siteind(SiteType(tag),n;kwargs...)

# Special case of `siteind` where integer (dim) provided
# instead of a tag string
siteind(d::Integer,n::Integer; kwargs...) = Index(d,"Site,n=$n")

#---------------------------------------
#
# siteinds system
#
#---------------------------------------

siteinds(::SiteType, N; kwargs...) = nothing

function siteinds(tag::String,
                  N::Integer; kwargs...)
  st = SiteType(tag)

  si = siteinds(st,N;kwargs...)
  if !isnothing(si)
    return si
  end

  return [siteind(st,j; kwargs...) for j=1:N]
end

function siteinds(f::Function,
                  N::Integer; kwargs...)
  [siteind(f(n),n; kwargs...) for n=1:N]
end

# Special case of `siteinds` where integer (dim)
# provided instead of a tag string
function siteinds(d::Integer,
                  N::Integer; kwargs...)
  return [siteind(d,n; kwargs...) for n=1:N]
end

#---------------------------------------
#
# has_fermion_string system
#
#---------------------------------------

has_fermion_string(::OpName, ::SiteType) = nothing

function has_fermion_string(opname::AbstractString,
                            s::Index;
                            kwargs...)::Bool
  opname = strip(opname)
  Ntags = length(tags(s))
  stypes = _sitetypes(s)
  opn = OpName(opname)
  for st in stypes
    res = has_fermion_string(opn, st)
    !isnothing(res) && return res
  end
  return false
end


text:
using LinearAlgebra
using Test
using Random

using ProximalOperators
using ProximalAlgorithms

@testset "Lasso small (strongly convex, $T)" for T in [Float32, Float64]

    Random.seed!(777)
    dim = 5
    μf = T(1)
    Lf = T(10)

    x_star = convert(Vector{T}, 1.5 * rand(T, dim) .- 0.5)

    lam = (μf + Lf) / 2
    @test typeof(lam) == T

    D = Diagonal(sqrt(μf) .+ (sqrt(Lf) - sqrt(μf)) * rand(T, dim))
    D[1] = sqrt(μf)
    D[end] = sqrt(Lf)
    Q = qr(rand(T, (dim, dim))).Q
    A = Q * D * Q'
    b = A * x_star + lam * inv(A') * sign.(x_star)

    f = LeastSquares(A, b)
    h = NormL1(lam)

    TOL = T(1e-4)

    @testset "SFISTA" begin

        # SFISTA

        x0 =  A \ b
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.SFISTA(tol = TOL)
        y, it = solver(x0, f = f, h = h, Lf = Lf, μf = μf)
        @test eltype(y) == T
        @test norm(y - x_star) <= TOL
        @test it < 200
        @test x0 == x0_backup

    end

end


text:
module SympatricPopulation

export make_genome,
    make_diploid_genome,
    spawn_population,
    length,
    recombinate,
    reproduce,
    HybridPopulation,
    proportion_species1,
    crosspollination_sampling,
    plot_population_histories

using Distributions
using Plots
using Random
using Statistics
using StatsBase

import Base.print

# Constructors
"""
    HaploidGenome

Simulates a haploid genome as an array, each count representing a chromosome.
Each chromosome can be either `0` or `1` for two different alleles.
"""
struct HaploidGenome
    genome::Array{Int8}
end

"""
    DiploidGenome

Simulates a diploid genome as two `HaploidGenome` instances.
"""
struct DiploidGenome
    genome1::HaploidGenome
    genome2::HaploidGenome
end

"""
    HybridPopulation

Simulates a population `Nₑ` of individual `DiploidGenome` instances.
"""
mutable struct HybridPopulation
    population::Array
end

```
    PopulationHistory

Store history of a population in a 2d array.
```
mutable struct PopulationHistory
    diary::Array{Float64,2}
end

"""
    PopulationInTime

Simulate a population over time.
Storing properties of the population (computed by `measure_function`) in
an array of history (`PopulationHistory`)
"""
struct PopulationInTime
    population::HybridPopulation
    history::PopulationHistory
    measure_function::Function
end

"""
    VarianceAndFixationHistory

Store histories of important properties of a population.
Here we use the variance of genome assignment, loss of one or another species
and hybrid overshoot (a hybrid has at least a proportion of 0.25 of a stranger
species).
"""
mutable struct VarianceAndFixationHistory
    variances::Array{Float64, 2}
    fixations::Array{Float64, 3}
    hybrid_overshoot::Array{Float64, 2}
end

# Functions
"""
    make_genome(n_chromosomes::Int64, p_allele::Float64)

Build a `HaploidGenome` of a given number of chromosomes (`n_chromosomes`).
The assignment of each chromosome to one or the other allele is done
by a Bernoulli trial with probability `p_allele` for `1`.
"""
function make_genome(n_chromosomes::Int64, p_allele::Float64)
    HaploidGenome(
        convert(
            Array{Int64},
            rand(Bernoulli(p_allele), n_chromosomes)
            )
        )
end

"""
    make_diploid_genome(n_chromosomes::Int64, p_allele::Float64)

Build a `DiploidGenome` with `n_chromosomes` and `p_allele` as
2 `HaploidGenome` instances.
"""
function make_diploid_genome(n_chromosomes::Int64, p_allele::Float64)
    DiploidGenome(make_genome(n_chromosomes, p_allele),
        make_genome(n_chromosomes, p_allele))
end

"""
    spawn_population(Nₑ::Int64, n_chromosomes::Int64, p_allele::Float64)

Build a `DiploidGenome` with `n_chromosomes` and `p_allele` as
2 `HaploidGenome` instances.
"""
function spawn_population(Nₑ::Int64, n_chromosomes::Int64, p_allele::Float64)
    population = [
        make_diploid_genome(n_chromosomes, rand(Bernoulli(p_allele), 1)[1] ? 1. : 0.)
        for i in 1:Nₑ]
    HybridPopulation(population)
end

"""
    length(d_g::DiploidGenome)

Retrieve number of chromosomes of each `HaploidGenome` representing
the genome length.
Both genomes should be same length.
"""
function length(d_g::DiploidGenome)
    @assert size(d_g.genome1.genome)[1] == size(d_g.genome2.genome)[1]
    size(d_g.genome1.genome)[1]
end

"""
    length(h_p::HybridPopulation)

Retrieve number of indiviuals in a `HybridPopulation` representing
the genome length.
Both genomes should be same length.
"""
function length(h_p::HybridPopulation)
    size(h_p.population)[1]
end

Base.deepcopy(h_g::HaploidGenome) = begin
    HaploidGenome(deepcopy(h_g.genome))
end

Base.deepcopy(d_g::DiploidGenome) = begin
    DiploidGenome(deepcopy(d_g.genome1), deepcopy(d_g.genome2))
end

"""
    proportion_species1(d_g)

Retrieve proportion of genome assigned to the species encoded with 1.
"""
function proportion_species1(d_g::DiploidGenome)
    (sum(d_g.genome1.genome) + sum(d_g.genome2.genome)) / (2*length(d_g))
end

"""
    print(diploid_genome::DiploidGenome)

Print the chromosomes in given order for both genomes.
"""
function print(diploid_genome::DiploidGenome)
    println(" - Genome 1: ", diploid_genome.genome1.genome...)
    println(" - Genome 2: ", diploid_genome.genome2.genome...)
end

"""
    print(hybrid_population::HybridPopulation)

Print the genomes of specimen in a hybrid population.
"""
function print(hybrid_population::HybridPopulation)
    println("Hybrid population (",
        length(hybrid_population),
        " dipl. individuals, ",
        length(hybrid_population.population[1]),
        " chromosomes):")
    for (i, individual) in enumerate(hybrid_population.population)
        println("- Individual ", i, ":")
        print(individual)
    end
end

"""
    recombinate(d_g::DiploidGenome, p_recombination::Float64; shuffled_array::Array=nothing)

Make random recombination events X, with X ~ Ber(`p_recombination`).
"""
function recombinate(d_g::DiploidGenome, p_recombination::Float64; shuffled_array::Array=[])
    d_g_new = deepcopy(d_g)
    if shuffled_array == []
        genome_length = length(d_g_new)
	shuffled_array = shuffle(1:genome_length)
    end
    switch = rand(Bool)
    for i in shuffled_array
        if rand(Bernoulli(p_recombination), 1)[1]
            switch = !switch
        end
        if switch
            d_g_new.genome1.genome[i], d_g_new.genome2.genome[i] = d_g_new.genome2.genome[i], d_g_new.genome1.genome[i]
        end
    end
    d_g_new
end

"""
    reproduce(parent1::DiploidGenome, parent2::DiploidGenome, p_meioticdrive::Float64)

Simulate reproduction of two diplod genomes by taking on of each parents genomes.
`p_meioticdrive` defines the probability to choose the first genome of both parents
"""
function reproduce(parent1::DiploidGenome, parent2::DiploidGenome, p_meioticdrive::Float64)
    genome1 = rand(Bernoulli(p_meioticdrive), 1)[1] ?
        deepcopy(parent1.genome1) : deepcopy(parent1.genome2)
    genome2 = rand(Bernoulli(p_meioticdrive), 1)[1] ?
        deepcopy(parent2.genome1) : deepcopy(parent2.genome2)
    DiploidGenome(genome1, genome2)
end

"""
    store_history(
        population_history::PopulationHistory,
        population::HybridPopulation,
        measure_function::Function)

Store the history of a population evolving over time.
The abstraction level of history can be defined by using `measure_function`.
"""
function store_history(
    population_history::PopulationHistory, population::HybridPopulation,
    measure_function::Function, position::Int64)
    @assert size(population_history.diary)[1] == length(population)
    new_diary_entry = sort([measure_function(d_g) for d_g in population.population])
    population_history.diary[:, position] = new_diary_entry;
end

"""
    initiate_phistory(Nₑ::Int64, n_cycles::Int64)

Make a new `PopulationHistory` instance with a given population size `Nₑ`
and `n_cycles` generations.
"""
function initiate_phistory(Nₑ::Int64, n_cycles::Int64)
    PopulationHistory(zeros(Nₑ, n_cycles))
end

"""
    population_dt(
        Nₑ::Int64, n_cycles::Int64; n_chromosomes::Int64=10, n_turnover::Int64=1,
        p_allele::Float64=0.5, p_recombination::Float64=0.5, p_meioticdrive::Float64=0.5,
        p_crosspollination::Float64=0.5, homoz_lethality::Float64=0.1,
        measure_function::Function=proportion_species1, fitness_function::Function=genome_balance,
        sample_function::Function=random_parent_sampling)

Simulate a population over time.
History is abstracted by `measure_function` and stored.
"""
function population_dt(
    Nₑ::Int64, n_cycles::Int64; n_chromosomes::Int64=10, n_turnover::Int64=1,
    p_allele::Float64=0.5, p_recombination::Float64=0.5, p_meioticdrive::Float64=0.5,
    p_crosspollination::Float64=0.5, homoz_lethality::Float64=0.1,
    measure_function::Function=proportion_species1, fitness_function::Function=genome_balance,
    sample_function::Function=random_parent_sampling
    )
    @assert Nₑ >= 3
    population = spawn_population(Nₑ, n_chromosomes, p_allele)
    popu_history = initiate_phistory(Nₑ, n_cycles)
    store_history(popu_history, population, measure_function, 1)
    for i in 2:n_cycles
        population = die_and_spawn(Nₑ, population, n_turnover,
            p_recombination, p_meioticdrive, p_crosspollination,
            homoz_lethality, measure_function, fitness_function,
            sample_function)
        store_history(popu_history, population, measure_function, i)
    end
    PopulationInTime(population, popu_history, measure_function)
end

"""
    die_and_spawn(
        Nₑ::Int64, population::HybridPopulation, n_turnover::Int64,
        p_recombination::Float64, p_meioticdrive::Float64)

Let `n_turnover` individuals of the population die and spawn new ones.
Dieing individuals are chosen randomly. Newborns are product of reproduction
of two sampled specimen (no selfing), which had recombination event
`p_recombination` as well as possible meiotic drive preferring one or the
other genome `p_meioticdrive`.
"""
function die_and_spawn(
    Nₑ::Int64, population::HybridPopulation, n_turnover::Int64,
    p_recombination::Float64, p_meioticdrive::Float64, p_crosspollination::Float64,
    homoz_lethality::Float64, measure_function::Function,
    fitness_function::Function, sample_function::Function
    )
    new_population =  Array{DiploidGenome, 1}(undef, Nₑ)
    death_said = kill_the_hybrids(
        Nₑ, n_turnover, population, homoz_lethality, fitness_function)
    # chromosome_shuffle = shuffle(1:length(population.population[1]))
    for i in 1:Nₑ
        if i in death_said
            (parent1, parent2) = sample_function(population, measure_function,
                p_crosspollination)
	    parent1 = recombinate(parent1, p_recombination) # , shuffled_array=chromosome_shuffle)
	    parent2 = recombinate(parent2, p_recombination) # , shuffled_array=chromosome_shuffle)
            offspringenome = reproduce(parent1, parent2, p_meioticdrive)
            new_population[i] = offspringenome
        else
            new_population[i] = population.population[i]
        end
    end
    new_population = shuffle(new_population)
    HybridPopulation(new_population)
end

"""
    random_parent_sampling(
        population::HybridPopulation, m_f::Function, p_cp::Float64)

Sample random parent pair from population.
"""
function random_parent_sampling(
    population::HybridPopulation, m_f::Function, p_cp::Float64)
    StatsBase.sample(
                population.population, 2; replace=false)
    # if `replace == true`, one could allow selfing
end

"""
    crosspollination_sampling(
        population::HybridPopulation, measure_function::Function,
        p_crosspollination::Float64)

Sample parent pair from population like an average pollinator does.
`p_crosspollination` defines how much less the pollinators prefer the less
interesting species from the favorite species. Thus, it is always
0 ⪓ `p_crosspollination` ≤ 1.
"""
function crosspollination_sampling(
    population::HybridPopulation, measure_function::Function,
    p_crosspollination::Float64)
    if rand(Bernoulli(.5), 1)[1]
        weights = (1).- [measure_function(d_g) for d_g in population.population]
    else
        weights = [measure_function(d_g) for d_g in population.population]
    end
    weights = [p_crosspollination + (1 - p_crosspollination)*weight for weight in weights]
    StatsBase.sample(population.population, pweights(weights), 2; replace=false)
end

"""
    kill_the_hybrids(
        Nₑ::Int64, n_turnover::Int64, population::HybridPopulation,
        homoz_lethality::Float64, measure_function::Function)

Sample individuals to die with increased hybrid lethality.
`homoz_lethality` defines the probability of killing a homozygous
non-hybrid compared to a full hybrid.
Returns array of indices of individuals to kill.
"""
function kill_the_hybrids(
    Nₑ::Int64, n_turnover::Int64, population::HybridPopulation,
    homoz_lethality::Float64, measure_function::Function)
    ε = (1 - homoz_lethality) / .5
    weights = [measure_function(individual) * ε + homoz_lethality
        for individual in population.population]
    StatsBase.sample(1:Nₑ, pweights(weights), n_turnover; replace=false)
end

"""
    genome_balance(d_g::DiploidGenome)

Compute the proportion of heterozygous genome.
The return value is min({p₁, p₀}), where p₁, p₀ are the genomic
proportion of species 1 or 0 within one individual.
"""
function genome_balance(d_g::DiploidGenome)
    f(p) = p > .5 ? 1 - p : p
    p = f(proportion_species1(d_g))
    @assert p <= .5
    @assert p >= 0.
    p
end

"""
    uniform_kill(d_g::DiploidGenome)

Shape uniform distribution to kill without hybrid disadvantage.
"""
function uniform_kill(d_g::DiploidGenome)
    1
end

"""
    repeat_population_dt(
        Nₑ::Int64, n_cycles::Int64; n_chromosomes::Int64=10, n_turnover::Int64=1,
        n_iterations::Int64=10, p_allele::Float64=0.5, p_recombination::Float64=0.5,
        p_meioticdrive::Float64=0.5, p_crosspollination::Float64=0.5, homoz_lethality::Float64=0.1,
        measure_function::Function=proportion_species1, fitness_function::Function=genome_balance,
        sample_function::Function=random_parent_sampling
        )

Repeat the function `population_dt` for `n_iterations` times.
Return the mean of all population histories.
"""
function repeat_population_dt(
    Nₑ::Int64, n_cycles::Int64; n_chromosomes::Int64=10, n_turnover::Int64=1,
    n_iterations::Int64=10, p_allele::Float64=0.5, p_recombination::Float64=0.5,
    p_meioticdrive::Float64=0.5, p_crosspollination::Float64=0.5, homoz_lethality::Float64=0.1,
    measure_function::Function=proportion_species1, fitness_function::Function=genome_balance,
    sample_function::Function=random_parent_sampling
    )
    mean_population = nothing
    properties_history = VarianceAndFixationHistory(
        Array{Float64, 2}(undef, n_iterations, n_cycles),
        Array{Float64, 3}(undef, n_iterations, 2, n_cycles),
        Array{Float64, 2}(undef, n_iterations, n_cycles))
    for i in 1:n_iterations
        population = population_dt(
            Nₑ, n_cycles; n_chromosomes=n_chromosomes, n_turnover=n_turnover,
            p_allele=p_allele, p_recombination=p_recombination, p_meioticdrive=p_meioticdrive,
            p_crosspollination=p_crosspollination, homoz_lethality=homoz_lethality,
            measure_function=measure_function, fitness_function=fitness_function,
            sample_function=sample_function)
        if isnothing(mean_population)
            mean_population = deepcopy(population.history.diary)
        else
            mean_population = mean_population + population.history.diary
        end
        summarize_history(population.history, properties_history, i)
    end
    mean_population = (mean_population)./ n_iterations
    (mean_population, properties_history)
end

"""
    get_population_property(
        population_history::PopulationHistory, behavioural_function::Function)

Compute and return a given parameter for the current population
"""
function get_population_property(
    population_history::PopulationHistory, behavioural_function::Function)
    dropdims(mapslices(
        x -> behavioural_function(x), population_history.diary, dims=1),
        dims=1)
end

"""
    hybrid_overshoot(population_history_slice::Array)

Compute the frequency of hybrids withih a population.
"""
function hybrid_overshoot(population_history_slice::Array)
    HYBRID_CONTENT = .25  # must be < 0.5
    n_hybrids = sum((population_history_slice.>HYBRID_CONTENT).&(
            population_history_slice.<(1-HYBRID_CONTENT)))
    f_hybrids = n_hybrids / size(population_history_slice)[1]
    f_hybrids
end

"""
    lost_in_fixation(population_history_slice::Array, species1::Bool)

Compute the frequency of specimen of one species.
To be assigned to a species, an individual must share at least 0.9 of
the species genome.
"""
function lost_in_fixation(population_history_slice::Array, species1::Bool)
    HYBRID_CONTENT = .10  # must be < 0.5
    if species1
        n_homozygotes = sum((population_history_slice.>(1-HYBRID_CONTENT)))
    else
        n_homozygotes = sum((population_history_slice.<HYBRID_CONTENT))
    end
    f_homozygotes = n_homozygotes / size(population_history_slice)[1]
    f_homozygotes
end

function summarize_history(population_history::PopulationHistory,
    properties_history::VarianceAndFixationHistory, iteration_step::Int64)
    properties_history.variances[iteration_step, :] = get_population_property(
        population_history, var)
    properties_history.hybrid_overshoot[iteration_step, :] = get_population_property(
        population_history, hybrid_overshoot)
    loss_of_species = Array{Float64, 2}(undef, 2, size(population_history.diary)[2])
    for species in 0:1
        species1 = species==1
        loss_of_species[species+1, :] = get_population_property(
            population_history, x -> lost_in_fixation(x, species1))
    end
    properties_history.fixations[iteration_step, :, :] = loss_of_species
end

"""

Make different analytical plots of population histories.
"""
function plot_population_histories(
    Nₑ::Int64, n_cycles::Int64, pollination_range::Array;
    variance_plot::Bool=false, mean_heatmap::Bool=true, fixation_plot::Bool=false,
    hybridization_content_plot::Bool=false,
    n_chromosomes::Int64=10, n_turnover::Int64=1,
    n_iterations::Int64=10, p_allele::Float64=0.5, p_recombination::Float64=0.5,
    p_meioticdrive::Float64=0.5, homoz_lethality::Float64=0.1,
    measure_function::Function=proportion_species1, fitness_function::Function=genome_balance,
    sample_function::Function=random_parent_sampling, postfix::String="")
    HYBRID_OVERSHOOT = .1
    gr()
    vp = variance_plot ? plot() : nothing
    fp = fixation_plot ? plot() : nothing
    hcp = hybridization_content_plot ? plot() : nothing
    mhp_array = Array{Plots.Plot, 1}(undef, size(pollination_range)[2])
    for (j, i) in enumerate(pollination_range)
        population_history, properties_history = repeat_population_dt(
            Nₑ, n_cycles; n_chromosomes=n_chromosomes, n_turnover=n_turnover, p_allele=p_allele, p_recombination=p_recombination, p_meioticdrive=p_meioticdrive, n_iterations=n_iterations,
            p_crosspollination=i, homoz_lethality=homoz_lethality,
            measure_function=measure_function, fitness_function=fitness_function,
            sample_function=sample_function);
	variance_plot ? plot!(vp, (1:n_cycles), mapslices(sum, properties_history.hybrid_overshoot.>HYBRID_OVERSHOOT,
            dims=1)', linewidth=1, label=string(i)) : nothing
	hybridization_content_plot ? plot!(hcp, (1:n_cycles), mapslices(x -> mean(x), properties_history.variances, dims=1)',
            linewidth=.5, label=string(i)) : nothing
        fixation_plot ? plot!(fp, (1:n_cycles), dropdims(mapslices(mean,
            mapslices(x -> minimum(x), properties_history.fixations, dims=(2)),
                dims=1), dims=1)', linewidth=.5, label=string(i)) : nothing
        if mean_heatmap
            mhp_array[j] = heatmap(population_history, title=string(i), colorbar=:none)
        end
        println(string("step", i))
    end
    mhp = mean_heatmap ? plot(mhp_array..., size = (1200, 800), link = :y) : nothing
    mhp = mean_heatmap ? savefig(mhp, string("figures/heatmap_n_e", Nₑ , postfix, ".svg")) : nothing
    fixation_plot ? savefig(
        fp, string("figures/speciesextinction_n_e", Nₑ , postfix, ".svg")) : nothing
    variance_plot ? savefig(
        vp, string("figures/variance_n_e", Nₑ , postfix, ".svg")) : nothing
    hybridization_content_plot ? savefig(
        hcp, string("figures/hybrid_fraction_n_e", Nₑ , postfix, ".svg")) : nothing
end
# Main

end  # module

max_stars_repo_name: biothomme/Geosmin

max_issues_repo_name: biothomme/Geosmin

max_forks_repo_name: biothomme/Geosmin


text:
@testset "Minimum cost flow" begin

    # bipartite oriented + source & sink
    # source 5, sink 6, v1 {1, 2} v2 {3, 4}
    g = Graphs.DiGraph(6)
    Graphs.add_edge!(g, 5, 1)
    Graphs.add_edge!(g, 5, 2)
    Graphs.add_edge!(g, 3, 6)
    Graphs.add_edge!(g, 4, 6)
    Graphs.add_edge!(g, 1, 3)
    Graphs.add_edge!(g, 1, 4)
    Graphs.add_edge!(g, 2, 3)
    Graphs.add_edge!(g, 2, 4)
    cost = zeros(6, 6)
    cost[1, 3] = 10.0
    cost[1, 4] = 5.0
    cost[2, 3] = 2.0
    cost[2, 4] = 2.0
    # v2 -> sink have demand of one
    demand = spzeros(6, 6)
    demand[3, 6] = 1
    demand[4, 6] = 1
    capacity = ones(6, 6)

    o = GLPK.Optimizer
    node_demand = spzeros(Graphs.nv(g))
    flow = minimum_cost_flow(
        g,
        node_demand,
        capacity,
        cost;
        optimizer=o,
        edge_demand=demand,
        source_nodes=[5],
        sink_nodes=[6],
    )
    @test flow[5, 1] ≈ 1
    @test flow[5, 2] ≈ 1
    @test flow[3, 6] ≈ 1
    @test flow[4, 6] ≈ 1
    @test flow[1, 4] ≈ 1
    @test flow[2, 4] ≈ 0
    @test flow[2, 3] ≈ 1
    @test flow[1, 3] ≈ 0

    # flow conservation property
    for n1 in 1:4
        @test sum(flow[n1, nout] for nout in outneighbors(g, n1)) ≈
            sum(flow[nin, n1] for nin in inneighbors(g, n1))
    end

    #--- Same test without edge_demands & positive cost should result in null flow
    flow = minimum_cost_flow(
        g, node_demand, capacity, cost; optimizer=o, source_nodes=[5], sink_nodes=[6]
    )

    @test flow[5, 1] ≈ 0
    @test flow[5, 2] ≈ 0
    @test flow[3, 6] ≈ 0
    @test flow[4, 6] ≈ 0
    @test flow[1, 4] ≈ 0
    @test flow[2, 4] ≈ 0
    @test flow[2, 3] ≈ 0
    @test flow[1, 3] ≈ 0

    # flow conservation property
    for n1 in 1:4
        @test sum(flow[n1, :]) ≈ sum(flow[:, n1])
    end

    #run same test again with nodal demands
    node_demand = spzeros(6)
    node_demand[5] = -2
    node_demand[6] = 2
    #--- Same test without edge_demands
    flow = minimum_cost_flow(
        g, node_demand, capacity, cost; optimizer=o, source_nodes=[5], sink_nodes=[6]
    )
    @test flow[5, 1] ≈ 1
    @test flow[5, 2] ≈ 1
    @test flow[3, 6] ≈ 1
    @test flow[4, 6] ≈ 1
    @test flow[1, 4] ≈ 1
    @test flow[2, 4] ≈ 0
    @test flow[2, 3] ≈ 1
    @test flow[1, 3] ≈ 0
    @test sum(diag(flow)) ≈ 0

    # flow conservation property
    for n1 in 1:4
        @test sum(flow[n1, :]) ≈ sum(flow[:, n1])
    end

    # no demand => null flow
    d2 = spzeros(6, 6)
    flow = minimum_cost_flow(
        g,
        spzeros(Graphs.nv(g)),
        capacity,
        cost;
        optimizer=o,
        edge_demand=d2,
        source_nodes=[5],
        sink_nodes=[6],
    )
    for i in 1:6
        for j in 1:6
            @test flow[i, j] ≈ 0.0
        end
    end

    # graph without sink or source
    g = Graphs.DiGraph(6)
    # create circuit
    for s in 1:5
        Graphs.add_edge!(g, s, s + 1)
    end
    Graphs.add_edge!(g, 6, 1)
    Graphs.add_edge!(g, 2, 5) # shortcut

    capacity = ones(6, 6)
    demand = spzeros(6, 6)
    demand[1, 2] = 1
    costs = ones(6, 6)
    flow = minimum_cost_flow(
        g, spzeros(Graphs.nv(g)), capacity, costs; optimizer=o, edge_demand=demand
    )
    active_flows = [(1, 2), (2, 5), (5, 6), (6, 1)]
    for s in 1:6
        for t in 1:6
            if (s, t) in active_flows
                @test flow[s, t] ≈ 1.0
            else
                @test flow[s, t] ≈ 0.0
            end
        end
    end
    # flow conservation property
    for n1 in 1:6
        @test sum(flow[n1, :]) ≈ sum(flow[:, n1])
    end
    # higher short-circuit cost
    costs[2, 5] = 10.0
    flow = minimum_cost_flow(
        g, spzeros(Graphs.nv(g)), capacity, costs; optimizer=o, edge_demand=demand
    )
    active_flows = [(1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 1)]
    for s in 1:6
        for t in 1:6
            if (s, t) in active_flows
                @test flow[s, t] ≈ 1.0
            else
                @test flow[s, t] ≈ 0.0
            end
        end
    end
    # flow conservation property
    for n1 in 1:6
        @test sum(flow[n1, :]) ≈ sum(flow[:, n1])
    end
end


text:
# Canonical form of multivariate normal

### Generic types

immutable MvNormalCanon{P<:AbstractPDMat,V<:Union(Vector{Float64},ZeroVector{Float64})} <: AbstractMvNormal
    μ::V    # the mean vector
    h::V    # potential vector, i.e. inv(Σ) * μ
    J::P    # precision matrix, i.e. inv(Σ)
end

typealias FullNormalCanon MvNormalCanon{PDMat,Vector{Float64}}
typealias DiagNormalCanon MvNormalCanon{PDiagMat,Vector{Float64}}
typealias IsoNormalCanon  MvNormalCanon{ScalMat,Vector{Float64}}

typealias ZeroMeanFullNormalCanon MvNormalCanon{PDMat,ZeroVector{Float64}}
typealias ZeroMeanDiagNormalCanon MvNormalCanon{PDiagMat,ZeroVector{Float64}}
typealias ZeroMeanIsoNormalCanon  MvNormalCanon{ScalMat,ZeroVector{Float64}}

### Constructors

function MvNormalCanon{P<:AbstractPDMat}(μ::Vector{Float64}, h::Vector{Float64}, J::P)
    length(μ) == length(h) == dim(J) || throw(DimensionMismatch("Inconsistent argument dimensions"))
    MvNormalCanon{P,Vector{Float64}}(μ, h, J)
end

function MvNormalCanon{P<:AbstractPDMat}(J::P)
    z = ZeroVector(Float64, dim(J))
    MvNormalCanon{P,ZeroVector{Float64}}(z, z, J)
end

function MvNormalCanon{P<:AbstractPDMat}(h::Vector{Float64}, J::P)
    length(h) == dim(J) || throw(DimensionMismatch("Inconsistent argument dimensions"))
    MvNormalCanon{P,Vector{Float64}}(J \ h, h, J)
end

MvNormalCanon(h::Vector{Float64}, J::Matrix{Float64}) = MvNormalCanon(h, PDMat(J))
MvNormalCanon(h::Vector{Float64}, prec::Vector{Float64}) = MvNormalCanon(h, PDiagMat(prec))
MvNormalCanon(h::Vector{Float64}, prec::Float64) = MvNormalCanon(h, ScalMat(length(h), prec))

MvNormalCanon(J::Matrix{Float64}) = MvNormalCanon(PDMat(J))
MvNormalCanon(prec::Vector{Float64}) = MvNormalCanon(PDiagMat(prec))
MvNormalCanon(d::Int, prec::Float64) = MvNormalCanon(ScalMat(d, prec))

### Show

distrname(d::IsoNormalCanon) = "IsoNormalCanon"
distrname(d::DiagNormalCanon) = "DiagNormalCanon"
distrname(d::FullNormalCanon) = "FullNormalCanon"

distrname(d::ZeroMeanIsoNormalCanon) = "ZeroMeanIsoNormalCanon"
distrname(d::ZeroMeanDiagNormalCanon) = "ZeroMeanDiagormalCanon"
distrname(d::ZeroMeanFullNormalCanon) = "ZeroMeanFullNormalCanon"

### conversion between conventional form and canonical form

meanform{C,V}(d::MvNormalCanon{C,V}) = MvNormal{C,V}(d.μ, inv(d.J))

canonform{C}(d::MvNormal{C,Vector{Float64}}) = (J = inv(d.Σ); MvNormalCanon(d.μ, J * d.μ, J))
canonform{C}(d::MvNormal{C,ZeroVector{Float64}}) = MvNormalCanon(inv(d.Σ))

### Basic statistics

length(d::MvNormalCanon) = length(d.μ)
mean(d::MvNormalCanon) = convert(Vector{Float64}, d.μ)

var(d::MvNormalCanon) = diag(inv(d.J))
cov(d::MvNormalCanon) = full(inv(d.J))
invcov(d::MvNormalCanon) = full(d.J)
logdetcov(d::MvNormalCanon) = -logdet(d.J)

### Evaluation

sqmahal(d::MvNormalCanon, x::DenseVector{Float64}) = quad(d.J, x - d.μ)
sqmahal!(r::DenseVector{Float64}, d::MvNormalCanon, x::DenseMatrix{Float64}) = quad!(r, d.J, x .- d.μ)

# Sampling (for GenericMvNormal)

unwhiten_winv!(J::AbstractPDMat, x::DenseVecOrMat) = unwhiten!(inv(J), x)
unwhiten_winv!(J::PDiagMat, x::DenseVecOrMat) = whiten!(J, x)
unwhiten_winv!(J::ScalMat, x::DenseVecOrMat) = whiten!(J, x)

_rand!(d::MvNormalCanon, x::DenseVecOrMat) = add!(unwhiten_winv!(d.J, randn!(x)), d.μ)


text:
# CompilerInstance
CompilerInstance() = CompilerInstance(create_compiler_instance())

"""
    create_compiler_instance() -> CXCompilerInstance
Return a pointer to a `clang::CompilerInstance` object.
"""
function create_compiler_instance()
    status = Ref{CXInit_Error}(CXInit_NoError)
    instance = clang_CompilerInstance_create(status)
    @assert status[] == CXInit_NoError
    return instance
end

dispose(x::CompilerInstance) = clang_CompilerInstance_dispose(x)

# Diagnostics
function hasDiagnostics(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_hasDiagnostics(ci)
end

function getDiagnostics(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasDiagnostics(ci) "CompilerInstance has no diagnostics engine."
    return DiagnosticsEngine(clang_CompilerInstance_getDiagnostics(ci))
end

function setDiagnostics(ci::CompilerInstance, diag::DiagnosticsEngine)
    @check_ptrs ci diag
    return clang_CompilerInstance_setDiagnostics(ci, diag)
end

function getDiagnosticClient(ci::CompilerInstance)
    @check_ptrs ci
    return DiagnosticConsumer(clang_CompilerInstance_getDiagnosticClient(ci))
end

function createDiagnostics(ci::CompilerInstance,
                           client::DiagnosticConsumer=DiagnosticConsumer(C_NULL),
                           should_own_client=true)
    @check_ptrs ci
    return clang_CompilerInstance_createDiagnostics(ci, client, should_own_client)
end

function setShowPresumedLoc(ci::CompilerInstance, should_show::Bool=true)
    opt = getDiagnosticOpts(ci)
    return setShowPresumedLoc(opt, should_show)
end

function setShowColors(ci::CompilerInstance, should_show::Bool=true)
    opt = getDiagnosticOpts(ci)
    return setShowColors(opt, should_show)
end

# FileManager
function hasFileManager(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_hasFileManager(ci)
end

function getFileManager(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasFileManager(ci) "CompilerInstance has no file manager."
    return FileManager(clang_CompilerInstance_getFileManager(ci))
end

function setFileManager(ci::CompilerInstance, file_mgr::FileManager)
    @check_ptrs ci file_mgr
    return clang_CompilerInstance_setFileManager(ci, file_mgr)
end

function createFileManager(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_createFileManager(ci)
end

function getFileEntry(ci::CompilerInstance, filename::AbstractString, open_file::Bool=true)
    file_mgr = getFileManager(ci)
    return getFileEntry(file_mgr, filename; open_file)
end

# SourceManager
function hasSourceManager(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_hasSourceManager(ci)
end

function getSourceManager(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasSourceManager(ci) "CompilerInstance has no source manager."
    return SourceManager(clang_CompilerInstance_getSourceManager(ci))
end

function setSourceManager(ci::CompilerInstance, src_mgr::SourceManager)
    @check_ptrs ci src_mgr
    return clang_CompilerInstance_setSourceManager(ci, src_mgr)
end

function createSourceManager(ci::CompilerInstance, src_mgr::FileManager)
    @check_ptrs ci src_mgr
    @assert hasDiagnostics(ci) "CompilerInstance has no diagnostics."
    return clang_CompilerInstance_createSourceManager(ci, src_mgr)
end

function createSourceManager(ci::CompilerInstance)
    file_mgr = getFileManager(ci)
    return createSourceManager(ci, file_mgr)
end

function setMainFileID(ci::CompilerInstance, filename::AbstractString, open_file::Bool=true)
    file_entry = getFileEntry(ci, filename, open_file)
    src_mgr = getSourceManager(ci)
    return setMainFileID(src_mgr, file_entry)
end

"""
    get_main_file_id(ci::CompilerInstance) -> FileID
Return the main file ID.

This function allocates and one should call `dispose` to release the resources after using this object.
"""
getMainFileID(ci::CompilerInstance) = getMainFileID(getSourceManager(ci))

# Invocation
function hasInvocation(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_hasInvocation(ci)
end

function getInvocation(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasInvocation(ci) "CompilerInstance has no invocation"
    return CompilerInvocation(clang_CompilerInstance_getInvocation(ci))
end

function setInvocation(ci::CompilerInstance, cinv::CompilerInvocation)
    @check_ptrs ci cinv
    return clang_CompilerInstance_setInvocation(ci, cinv)
end

# Target
function hasTarget(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_hasTarget(ci)
end

function getTarget(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasTarget(ci) "CompilerInstance has no target"
    return TargetInfo(clang_CompilerInstance_getTarget(ci))
end

function setTargetAndLangOpts(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_setTargetAndLangOpts(ci)
end

function setTarget(ci::CompilerInstance, tgti::TargetInfo)
    @check_ptrs ci tgti
    return clang_CompilerInstance_setTarget(ci, tgti)
end

# Preprocessor
function hasPreprocessor(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_hasPreprocessor(ci)
end

function getPreprocessor(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasPreprocessor(ci) "CompilerInstance has no preprocessor"
    return Preprocessor(clang_CompilerInstance_getPreprocessor(ci))
end

function setPreprocessor(ci::CompilerInstance, pp::Preprocessor)
    @check_ptrs ci pp
    return clang_CompilerInstance_setPreprocessor(ci, pp)
end

function createPreprocessor(ci::CompilerInstance, kind=CXTranslationUnitKind_TU_Complete)
    @check_ptrs ci
    @assert hasTarget(ci) "CompilerInstance has no target."
    return clang_CompilerInstance_createPreprocessor(ci, kind)
end

# Sema
function hasSema(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_hasSema(ci)
end

function getSema(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasSema(ci) "CompilerInstance has no sema."
    return Sema(clang_CompilerInstance_getSema(ci))
end

function setSema(ci::CompilerInstance, sema::Sema)
    @check_ptrs ci sema
    return clang_CompilerInstance_setSema(ci, sema)
end

function createSema(ci::CompilerInstance, kind=CXTranslationUnitKind_TU_Complete)
    @check_ptrs ci
    @assert hasASTContext(ci) "CompilerInstance has no ASTContext."
    @assert hasASTConsumer(ci) "CompilerInstance has no ASTConsumer."
    return clang_CompilerInstance_createSema(ci, kind)
end

# ASTContext
function hasASTContext(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_hasASTContext(ci)
end

function getASTContext(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasASTContext(ci) "CompilerInstance has no AST context."
    return ASTContext(clang_CompilerInstance_getASTContext(ci))
end

function setASTContext(ci::CompilerInstance, ctx::ASTContext)
    @check_ptrs ci ctx
    return clang_CompilerInstance_setASTContext(ci, ctx)
end

function createASTContext(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasPreprocessor(ci) "CompilerInstance has no preprocessor."
    return clang_CompilerInstance_createASTContext(ci)
end

# ASTConsumer
function hasASTConsumer(ci::CompilerInstance)
    @check_ptrs ci
    return clang_CompilerInstance_hasASTConsumer(ci)
end

function getASTConsumer(ci::CompilerInstance)
    @check_ptrs ci
    @assert hasASTConsumer(ci) "CompilerInstance has no AST consumer."
    return ASTConsumer(clang_CompilerInstance_getASTConsumer(ci))
end

function setASTConsumer(ci::CompilerInstance, csr::AbstractASTConsumer)
    @check_ptrs ci csr
    return clang_CompilerInstance_setASTConsumer(ci, csr)
end

# CodeGenerator
function CreateLLVMCodeGen(ci::CompilerInstance, ctx::Context, mod_name::String="JLCC")
    @check_ptrs ci
    @assert ctx.ref != C_NULL "Context has a NULL pointer."
    return CodeGenerator(clang_CreateLLVMCodeGen(ci, ctx, mod_name))
end

# Actions
function ExecuteAction(ci::CompilerInstance, action::T) where {T<:AbstractFrontendAction}
    @check_ptrs ci action
    return clang_CompilerInstance_ExecuteAction(ci, action)
end

# Options
function getCodeGenOpts(ci::CompilerInstance)
    @check_ptrs ci
    return CodeGenOptions(clang_CompilerInstance_getCodeGenOpts(ci))
end

function getDiagnosticOpts(ci::CompilerInstance)
    @check_ptrs ci
    return DiagnosticOptions(clang_CompilerInstance_getDiagnosticOpts(ci))
end

function getFrontendOpts(ci::CompilerInstance)
    @check_ptrs ci
    return FrontendOptions(clang_CompilerInstance_getFrontendOpts(ci))
end

function getHeaderSearchOpts(ci::CompilerInstance)
    @check_ptrs ci
    return HeaderSearchOptions(clang_CompilerInstance_getHeaderSearchOpts(ci))
end

function getPreprocessorOpts(ci::CompilerInstance)
    @check_ptrs ci
    return PreprocessorOptions(clang_CompilerInstance_getPreprocessorOpts(ci))
end

function getTargetOpts(ci::CompilerInstance)
    @check_ptrs ci
    return TargetOptions(clang_CompilerInstance_getTargetOpts(ci))
end

function getLangOpts(ci::CompilerInstance)
    @check_ptrs ci
    return LangOptions(clang_CompilerInstance_getLangOpts(ci))
end

# PrintStats
function PrintStats(ci::CompilerInstance, ::Type{CodeGenOptions})
    opts = getCodeGenOpts(ci)
    return PrintStats(opts)
end

function PrintStats(ci::CompilerInstance, ::Type{DiagnosticOptions})
    opts = getDiagnosticOpts(ci)
    return PrintStats(opts)
end

function PrintStats(ci::CompilerInstance, ::Type{FrontendOptions})
    opts = getFrontendOpts(ci)
    return PrintStats(opts)
end

function PrintStats(ci::CompilerInstance, ::Type{HeaderSearchOptions})
    opts = getHeaderSearchOpts(ci)
    return PrintStats(opts)
end

function PrintStats(ci::CompilerInstance, ::Type{PreprocessorOptions})
    opts = getPreprocessorOpts(ci)
    return PrintStats(opts)
end

function PrintStats(ci::CompilerInstance, ::Type{TargetOptions})
    opts = getTargetOpts(ci)
    return PrintStats(opts)
end

function PrintStats(ci::CompilerInstance, ::Type{LangOptions})
    opts = getLangOpts(ci)
    return PrintStats(opts)
end

function PrintStats(ci::CompilerInstance, ::Type{FileManager})
    fm = getFileManager(ci)
    return PrintStats(fm)
end

function PrintStats(ci::CompilerInstance, ::Type{SourceManager})
    sm = getSourceManager(ci)
    return PrintStats(sm)
end

function PrintStats(ci::CompilerInstance, ::Type{HeaderSearch})
    pp = getPreprocessor(ci)
    opts = getHeaderSearchInfo(pp)
    return PrintStats(opts)
end

function PrintStats(ci::CompilerInstance, ::Type{Preprocessor})
    pp = getPreprocessor(ci)
    return PrintStats(pp)
end

function PrintStats(ci::CompilerInstance, ::Type{Sema})
    s = getSema(ci)
    return PrintStats(s)
end

function PrintStats(ci::CompilerInstance, ::Type{ASTContext})
    ctx = getASTContext(ci)
    return PrintStats(ctx)
end

function PrintStats(ci::CompilerInstance, ::Type{ASTConsumer})
    ctx = getASTConsumer(ci)
    return PrintStats(ctx)
end


text:
using TMI
using Documenter

DocMeta.setdocmeta!(TMI, :DocTestSetup, :(using TMI); recursive=true)

makedocs(;
    modules=[TMI],
    authors="G Jake Gebbie",
    repo="https://github.com/ggebbie/TMI.jl/blob/{commit}{path}#{line}",
    sitename="TMI.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ggebbie.github.io/TMI.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Top-level functions" => "top.md",
        "Configuration" => "config.md",
        "Grid" => "grid.md",
        "Boundary Conditions" => "boundaries.md",
        "Plots" => "plots.md",
        "Utilities" => "utils.md",
        "Older functions" => "legacy.md",
    ],
)

deploydocs(;
           repo="github.com/ggebbie/TMI.jl",
           devbranch="main",
)