text:
using Gtk, Graphics

global flag=0
const dir = ( (-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1) )

mutable struct Board
    table::Matrix{Int}
    turn::Int # 1-> white, 0->black
    n::Int # linear size of board

    function Board()
        n = 8
        t = ones(Int, n, n) .* Int(-1)
        a = div(n, 2);
        t[a,a] = t[a+1, a+1] = 1
        t[a,a+1] = t[a+1, a] = 0
        new(t, 1, n)
    end
end

function flip!(b::Matrix{Int}, x::Int, y::Int, h::Int, v::Int, bound::Int, turn::Int)
    ( !(0 < x <= bound) || !(0 < y <= bound) ) && return
    b[x,y] == turn && return (global flag=1; nothing)
    flip!(b, x+h,y+v, h, v, bound, turn)
    flag == 1 && (b[x,y] = turn; nothing)
end
flip!(b::Board, x::Int, y::Int, h::Int, v::Int) = flip!(b.table, x, y, h, v, b.n, b.turn)

function initialize!(c::GtkCanvas)
    @guarded draw(c) do widget
        n = 8
        ctx = getgc(c)
        h = height(c)
        w = width(c)
        space = h/n
        r = space/2

        rectangle(ctx, 0, 0, w, h)
        set_source_rgb(ctx, 0, 0.6, 0)
        fill(ctx)

        for x in 1:n
            move_to(ctx, x*w/n, 0)
            line_to(ctx, x*w/n, h)
        end
        for y in 1:7
            move_to(ctx, 0, y*w/n)
            line_to(ctx, w, y*w/n)
        end
        set_source_rgb(ctx, 0, 0, 0)
        stroke(ctx)
        reveal(widget)
    end
    show(c)
end

function draw_board!(c::GtkCanvas, b::Board)
    table = b.table
    n = 8
    for row in 1:n
        for col in 1:n
            color = table[row,col]
            if color != -1
                @guarded draw(c) do widget
                    ctx = getgc(c)
                    h = height(c)
                    w = width(c)
                    space = h/n
                    r = space/2

                    set_source_rgb(ctx, color, color, color)
                    arc(ctx, (row-1)*space+r, (col-1)*space+r, 0.75*r, 0, 2pi)
                    fill(ctx)
                    reveal(widget)
                end
            end
        end
    end
    show(c)
end

b = Board()
c = @GtkCanvas()
n = 400
bnum = sum(b.table .== 0)
wnum = sum(b.table .== 1)
win = GtkWindow(c,
        "Reversi: White $(lpad(wnum, 2)), Black: $(lpad(bnum, 2)), Turn: " * (b.turn == 1 ? "White" : "Black"),
        n, n)
initialize!(c)
draw_board!(c, b)

c.mouse.button1press = @guarded (widget, event) -> begin
    n = 8
    turn = b.turn
    ctx = getgc(widget)
    set_source_rgb(ctx, turn, turn, turn)
    h = height(c)
    space = h/n
    r = space/2
    x, y = Int(div(event.x, space)), Int(div(event.y, space))
    x += 1; y += 1
    b.table[x, y] = b.turn
    for (h,v) in dir
        if 0 < x+h <= n && 0 < y+v <= n
            b.table[x+h, y+v] == -1 && continue
            flip!(b, x+h, y+v, h, v)
            global flag=0
        end
    end
    initialize!(c)
    draw_board!(c, b)
    draw_board!(c, b)
    reveal(widget)
    b.turn = xor(1, b.turn)
    bnum = sum(b.table .== 0)
    wnum = sum(b.table .== 1)
    set_gtk_property!(win, :title, "Reversi: White $(lpad(wnum, 2)), Black: $(lpad(bnum, 2)), Turn: " * (b.turn == 1 ? "White" : "Black"))
end
initialize!(c)
draw_board!(c,b)

if !isinteractive()
    cond = Condition()
    signal_connect(win, :destroy) do widget
        notify(cond)
    end
    wait(cond)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
# AbstractBandedMatrix must implement

# A BlockBandedMatrix is a BlockMatrix, but is not a BandedMatrix
abstract type AbstractBlockBandedMatrix{T} <: AbstractBlockMatrix{T} end

"""
    blockbandwidths(A)

Returns a tuple containing the upper and lower blockbandwidth of `A`.
"""
blockbandwidths(A::AbstractMatrix) = (nblocks(A,1)-1 , nblocks(A,2)-1)

"""
    blockbandwidth(A,i)

Returns the lower blockbandwidth (`i==1`) or the upper blockbandwidth (`i==2`).
"""
blockbandwidth(A::AbstractMatrix, k::Integer) = blockbandwidths(A)[k]

"""
    bandrange(A)

Returns the range `-blockbandwidth(A,1):blockbandwidth(A,2)`.
"""
blockbandrange(A::AbstractMatrix) = -blockbandwidth(A,1):blockbandwidth(A,2)

# start/stop indices of the i-th column/row, bounded by actual matrix size
@inline blockcolstart(A::AbstractVecOrMat, i::Integer) = Block(max(i-colblockbandwidth(A,2)[i], 1))
@inline  blockcolstop(A::AbstractVecOrMat, i::Integer) = Block(max(min(i+colblockbandwidth(A,1)[i], nblocks(A, 1)), 0))
@inline blockrowstart(A::AbstractVecOrMat, i::Integer) = Block(max(i-rowblockbandwidth(A,1)[i], 1))
@inline  blockrowstop(A::AbstractVecOrMat, i::Integer) = Block(max(min(i+rowblockbandwidth(A,2)[i], nblocks(A, 2)), 0))

for Func in (:blockcolstart, :blockcolstop, :blockrowstart, :blockrowstop)
    @eval $Func(A, i::Block{1}) = $Func(A, Int(i))
end

@inline blockcolrange(A::AbstractVecOrMat, i) = blockcolstart(A,i):blockcolstop(A,i)
@inline blockrowrange(A::AbstractVecOrMat, i) = blockrowstart(A,i):blockrowstop(A,i)

# length of i-the column/row
@inline blockcollength(A::AbstractVecOrMat, i) = max(Int(blockcolstop(A, i)) - Int(blockcolstart(A, i)) + 1, 0)
@inline blockrowlength(A::AbstractVecOrMat, i) = max(Int(blockrowstop(A, i)) - Int(blockrowstart(A, i)) + 1, 0)

# this gives the block bandwidth in each block column/row
@inline colblockbandwidths(A::AbstractMatrix) = Fill.(blockbandwidths(A), nblocks(A,2))
@inline rowblockbandwidths(A::AbstractMatrix) = Fill.(blockbandwidths(A), nblocks(A,1))

@inline colblockbandwidth(bs, i::Int) = colblockbandwidths(bs)[i]
@inline rowblockbandwidth(bs, i::Int) = rowblockbandwidths(bs)[i]

"""
    isblockbanded(A)

returns true if a matrix implements the block banded interface.
"""
isblockbanded(::AbstractBlockBandedMatrix) = true
isblockbanded(_) = false

# override bandwidth(A,k) for each AbstractBandedMatrix
# override inbands_getindex(A,k,j)

# return id of first empty diagonal intersected along row k
function _firstblockdiagrow(A::AbstractMatrix, k::Int)
    a, b = blockrowstart(A, k), blockrowstop(A, k)
    c = a == 1 ? b+1 : a-1
    c-k
end

# return id of first empty diagonal intersected along column j
function _firstblockdiagcol(A::AbstractMatrix, j::Int)
    a, b = blockcolstart(A, j), blockcolstop(A, j)
    r = a == 1 ? b+1 : a-1
    j-r
end

## BlockSlice1 is a conveneience for views

const BlockSlice1 = BlockSlice{Block{1,Int}}

######################################
#  RaggedMatrix interface
######################################

@inline function colstart(A::AbstractBlockBandedMatrix, i::Integer)
    bs = A.block_sizes.block_sizes
    bs.cumul_sizes[1][Int(blockcolstart(A, _find_block(bs, 2, i)[1]))]
end
@inline function colstop(A::AbstractBlockBandedMatrix, i::Integer)
    bs = A.block_sizes.block_sizes
    bs.cumul_sizes[1][Int(blockcolstop(A, _find_block(bs, 2, i)[1]))+1]-1
end
@inline function rowstart(A::AbstractBlockBandedMatrix, i::Integer)
    bs = A.block_sizes.block_sizes
    bs.cumul_sizes[2][Int(blockrowstart(A, _find_block(bs, 1, i)[1]))]
end
@inline function rowstop(A::AbstractBlockBandedMatrix, i::Integer)
    bs = A.block_sizes.block_sizes
    bs.cumul_sizes[2][Int(blockrowstop(A, _find_block(bs, 1, i)[1]))+1]-1
end

# default implementation loops over all indices, including zeros
function fill!(A::AbstractBlockBandedMatrix, val::Any)
  iszero(val) || throw(BandError(A))
  fill!(A.data, val)
  A
end

ext: jl


text:
title = "Algorithms"
weight = 20


text:
module MonteCarloMarkovKernels
using StaticArrays
using LinearAlgebra
using Random
import Statistics.mean

include("simulateChain.jl")
include("randomWalkMetropolis.jl")
include("adaptiveMetropolis.jl")
include("batchMeans.jl")
include("spectralVariance.jl")
include("visualize.jl")

export simulateChain!, simulateChain,
  simulateChainProgress!, simulateChainProgress,
  makeRWMKernel, makeAMKernel,
  estimateBM, estimateSV,
  kde

end


text:
#=
function Base.show(io::IO, op_problem::OperationsProblem)
    println(io, "Operation Model")
end
=#

function _organize_model(
    val::Dict{Symbol, T},
    field::Symbol,
    io::IO,
) where {T <: Union{DeviceModel, ServiceModel}}
    println(io, "  $(field): ")
    for (i, ix) in val
        println(io, "      $(i):")
        for inner_field in fieldnames(T)
            inner_field == :services && continue
            value = getfield(val[i], Symbol(inner_field))

            if !isnothing(value)
                println(io, "        $(inner_field) = $value")
            end
        end
    end

end

"""
    Base.show(io::IO, ::MIME"text/plain", op_problem::OperationsProblem)

This function goes through the fields in OperationsProblem and then in OperationsProblemTemplate,
if the field contains a Device model dictionary, it calls organize_device_model() &
prints the data by field, key, value. If the field is not a Device model dictionary,
and a value exists for that field it prints the value.

"""
function Base.show(io::IO, m::MIME"text/plain", op_problem::OperationsProblem)
    show(io, m, op_problem.template)
end

function Base.show(io::IO, ::MIME"text/plain", template::OperationsProblemTemplate)
    println(io, "\nOperations Problem Specification")
    println(io, "============================================\n")

    for field in fieldnames(OperationsProblemTemplate)
        val = getfield(template, Symbol(field))
        if typeof(val) <: Dict{Symbol, <:Union{DeviceModel, ServiceModel}}
            println(io, "============================================")
            _organize_model(val, field, io)
        else
            if !isnothing(val)
                println(io, "  $(field):  $(val)")
            else
                println(io, "no data")
            end
        end
    end
    println(io, "============================================")
end

function Base.show(io::IO, psi_container::PSIContainer)
    println(io, "PSIContainer()")
end
#=
function Base.show(io::IO, op_problem::SimulationSequence)
    # Here ASCII Art
    println(io, "SimulationSequence()")
end
=#
function Base.show(io::IO, sim::Simulation)
    println(io, "Simulation()")
end
#=
function Base.show(io::IO, results::OperationsProblemResults)
    println(io, "Results Model")
 end

=#

function Base.show(io::IO, ::MIME"text/plain", results::IS.Results)
    println(io, "\nResults")
    println(io, "===============\n")

    for (k, v) in results.variable_values
        time = DataFrames.DataFrame(Time = results.time_stamp[!, :Range])
        if size(time, 1) == size(v, 1)
            var = hcat(time, v)
        else
            var = v
        end
        println(io, "$(k)")
        println(io, "==================")
        println(io, "$(var)\n")
    end
    println(io, "Optimizer Log")
    println(io, "-------------")
    for (k, v) in results.optimizer_log
        if !isnothing(v)
            println(io, "        $(k) = $(v)")
        end
    end
    println(io, "\n")
    for (k, v) in results.total_cost
        println(io, "Total Cost: $(k) = $(v)")
    end
end

function Base.show(io::IO, ::MIME"text/html", results::IS.Results)
    println(io, "
Results
")
    for (k, v) in results.variable_values
        time = DataFrames.DataFrame(Time = results.time_stamp[!, :Range])
        if size(time, 1) == size(v, 1)
            var = hcat(time, v)
        else
            var = v
        end
        println(io, "$(k)")
        show(io, MIME"text/html"(), var)
    end
    println(io, "
Optimizer Log
")
    for (k, v) in results.optimizer_log
        if !isnothing(v)
            println(io, "
        $(k) = $(v)
")
        end
    end
    println(io, "\n")
    for (k, v) in results.total_cost
        println(io, "
Total Cost: $(v)
")
    end
end

function Base.show(io::IO, stage::Stage)
    println(io, "Stage()")
end

function Base.show(io::IO, ::MIME"text/html", services::Dict{Symbol, PSI.ServiceModel})
    println(io, "
Services
")
    for (k, v) in services
        println(io, "
$(k)
")
        println(io, "
$(v)
")
    end
end

function Base.show(io::IO, ::MIME"text/html", sim_results::SimulationResultsReference)
    println(io, "
Simulation Results Reference
")
    println(io, "
Results Folder: $(sim_results.results_folder)
")
    println(io, "
Reference Tables
")
    for (k, v) in sim_results.ref
        println(io, "
$(k)
")
        for (i, x) in v
            println(io, "
$(i): dataframe size $(size(x))
")
        end
    end
    for (k, v) in sim_results.chronologies
        println(io, "
$(k)
")
        println(io, "
time length: $(v)
")
    end
end

function Base.show(io::IO, ::MIME"text/plain", sim_results::SimulationResultsReference)
    println(io, "Simulation Results Reference\n")
    println(io, "Results Folder: $(sim_results.results_folder)\n")
    println(io, "Reference Tables\n")
    for (k, v) in sim_results.ref
        println(io, "$(k)\n")
        try
            for (i, x) in v
                println(io, "$(i): dataframe size $(size(x))\n")
            end
        catch
            print(io, "Custom Stage")
        end
    end
    for (k, v) in sim_results.chronologies
        println(io, "$(k)\n")
        println(io, "time length: $(v)\n")
    end
end

function _count_stages(sequence::Array)
    stages = Dict{Int64, Int64}()
    stage = 1
    count = 0
    for i in 1:length(sequence)
        if sequence[i] == stage
            count += 1
            if i == length(sequence)
                stages[stage] = count
            end
        else
            stages[stage] = count
            stage += 1
            count = 1
        end
    end
    return stages
end

function _print_feedforward(io::IO, feed_forward::Dict, to::Array, from::Any)
    for (keys, sync) in feed_forward
        period = sync.periods
        stage1 = keys[1]
        stage2 = keys[2]
        spaces = " "^(length(stage2) + 2)
        dashes = "-"^(length(stage2) + 2)
        if period <= 12
            times = period
            line5 = string("└─$stage2 "^times, "($period) to : $to")
        else
            times = 12
            line5 = string("└─$stage2 "^times, "... (x$period) to : $to")
        end
        if times == 1
            line1 = "$stage1--┐ from : $from"
            println(io, "$line1\n$spaces|\n$spaces$line5\n")
        else
            if times == 2
                line3 = string("┌", string(dashes, "┤"))
                spacing = 0
            elseif times == 3
                line3 = string("┌", string(dashes, "┼"), string(dashes, "┐"))
                spacing = 0
            elseif iseven(times)
                spacing = (Int(times / 2) - 2)
                line3 = string(
                    "┌",
                    string(dashes, "┬")^spacing,
                    "----",
                    "┼",
                    string(dashes, "┬")^(spacing + 1),
                    "----┐",
                )
            else
                spacing = Int((times / 2) - 1.5)
                line3 = string(
                    "┌",
                    string(dashes, "┬")^spacing,
                    "----",
                    "┼",
                    string(dashes, "┬")^(spacing),
                    "----┐",
                )
            end
            line1 = string("     "^(spacing), " $stage1--┐ from : $from")
            line2 = string("     "^(spacing), " "^length(stage1), "   |")
            line4 = string("|", string(spaces, "|")^(times - 2), "    |")
            println(io, "$line1\n$line2\n$line3\n$line4\n$line4\n$line5\n")
        end
    end
end
function _print_inter_stages(io::IO, stages::Dict{Int64, Int64})
    list = sort!(collect(keys(stages)))
    for i in list
        num = stages[i]
        total = length(list)
        if total > 5
            println(io, "Too many stages to print.")
        else
            if length("$num") == 1
                num = " (x0$num)"
            elseif length("$num") == 3
                num = " ($num)"
            elseif length("$num") == 4
                num = "($num)"
            else
                num = " (x$num)"
            end
            if i == 1
                if total == 2
                    println(io, "$i")
                else
                    println(io, "$i\n|")
                end
            else
                N = 2^(i - 2)
                space_count = 10 * (2^(total - i))
                if i == total
                    print1 = "|             ┌----/"^(N - 1)
                    print2 = "|             |     "^(N - 1)
                    print3 = "$i --> $i ...$num   "^N
                    println(io, "$print1|\n$print2|\n$print3")
                else
                    spaces = " "^(space_count - 11)
                    up = Int(space_count / 2)
                    print = string("$i ", " "^(space_count - 4), "  $i ...$num", spaces)^N
                    println(io, "$print")
                    if i !== total - 1
                        indent1 = string(
                            string("|", " "^(up + 7), "┌", "-"^(up - 10), "/")^(2 * N - 1),
                            "|",
                        )
                        indent2 = string("|", " "^(up + 7), "|", " "^(up - 9))^(2 * N - 1)
                        println(io, "$indent1\n$indent2|")
                    end
                end
            end
        end
    end
end

function _print_intra_stages(io::IO, stages::Dict{Int64, Int64})
    list = sort!(collect(keys(stages)))
    for i in list
        num = stages[i]
        total = length(list)
        if total > 5
            println(io, "Too many stages to print.")
        else
            if length("$num") == 1
                num = " (x0$num)"
            elseif length("$num") == 3
                num = " ($num)"
            elseif length("$num") == 4
                num = "($num)"
            else
                num = " (x$num)"
            end
            if i == 1
                println(io, "$i\n")
            else
                N = 2^(i - 2)
                space_count = 10 * (2^(total - i))
                if i == total
                    print = "$i --> $i ...$num   "^N
                    println(io, "$print\n\n")
                else
                    spaces = " "^(space_count - 11)
                    print = string("$i ", "-"^(space_count - 4), "> $i ...$num", spaces)^N
                    indent = string(" "^(space_count))^(N * 2)
                    println(io, "$print\n$indent")
                end
            end
        end
    end
end

function Base.show(io::IO, sequence::SimulationSequence)
    stages = _count_stages(sequence.execution_order)
    println(io, "Feed Forward Chronology")
    println(io, "-----------------------\n")
    to = []
    from = String("")
    for (k, v) in sequence.feedforward
        println(io, "$(k[1]): $(typeof(v)) -> $(k[3])\n")
        to = String.(v.affected_variables)
        if isa(v, SemiContinuousFF)
            from = String.(v.binary_from_stage)
        elseif isa(v, RangeFF)
            from = String.([v.variable_from_stage_ub, v.variable_from_stage_lb])
        else
            from = String.(v.variable_from_stage)
        end
        _print_feedforward(io, sequence.feedforward_chronologies, to, from)
    end
    println(io, "Initial Condition Chronology")
    println(io, "----------------------------\n")
    if sequence.ini_cond_chronology == IntraStageChronology()
        _print_intra_stages(io, stages)
    elseif sequence.ini_cond_chronology == InterStageChronology()
        _print_inter_stages(io, stages)
    end
end


text:
function Thermostat!(sys::System,comm::MPI.Comm,root::Int64)
    for thermostat in sys.thermostat
        if thermostat!=NoThermostat()
            apply_thermostat!(sys,thermostat,comm,root)
        end
    end
end

function apply_thermostat!(sys::System,thermost::BerendsenNVT,comm::MPI.Comm,root::Int64)
    if sys.current_temp[1]==0.0
        sys.current_temp[1]=1e-8
    end
    scalT  = sqrt(1.0+sys.dt*(sys.temp/sys.current_temp[1]-1.0)/thermost.τ)
    apply_thermostat_particle!.(sys.atoms.vels,sys.atoms.forces,sys.atoms.masses,scalT,sys.dt,Ref(thermost))
    #@show sys.forces
end

function apply_thermostat_particle!(vels::Array{Float64,1},
                                    forces::Array{Float64,1},
                                    masses::Float64,
                                    scalT::Float64,
                                    dt::Float64,
                                    thermost::BerendsenNVT)
    vels[1] = vels[1] * scalT
    vels[2] = vels[2] * scalT
    vels[3] = vels[3] * scalT
end

function apply_thermostat!(sys::System,thermost::LangevinNVT,comm::MPI.Comm,root::Int64)
    coeff=sqrt(6.0 * thermost.gamma * sys.temp*thermost.invdt)
    apply_thermostat_particle!.(sys.atoms.vels,sys.atoms.forces,sys.atoms.masses,coeff,sys.dt,Ref(thermost))
end

function apply_thermostat_particle!(vels::Array{Float64,1},
                                    forces::Array{Float64,1},
                                    masses::Float64,
                                    coeff::Float64,
                                    dt::Float64,
                                    thermost::LangevinNVT)
    r_force_x = RandomNumber()*coeff - thermost.gamma*vels[1]
    r_force_y = RandomNumber()*coeff - thermost.gamma*vels[2]
    r_force_z = RandomNumber()*coeff - thermost.gamma*vels[3]
    vels[1] += 0.5 * dt * r_force_x / masses
    vels[2] += 0.5 * dt * r_force_y / masses
    vels[3] += 0.5 * dt * r_force_z / masses
    forces[1] += r_force_x
    forces[2] += r_force_y
    forces[3] += r_force_z
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# for 1d and 2d histograms
function split_smallest(histogram, intervals)
    intervals = sort(intervals)

    hist = deepcopy(histogram)
    hists_weights = Array{Array{Real}}(undef, length(intervals))
    dims = size(hist.weights)

    weights = vec(hist.weights)
    totalweights = sum(weights)
    weight_ids = sortperm(weights, rev=true) # starting with highest weight

    for i in 1:length(intervals)
        hists_weights[i] = zeros(length(weights))
    end

    rel_weights = weights/totalweights

    for (i, intv) in enumerate(intervals)
        sum_weights = 0.
        for w in weight_ids
            if sum_weights < intv
                hists_weights[i][w] = hist.weights[w]
                sum_weights += rel_weights[w]
            else
                break
            end
        end
    end

    hists = Array{StatsBase.Histogram}(undef, length(intervals))

    for i in 1:length(intervals)
        hists[i] = deepcopy(hist)
        hists[i].weights = reshape(hists_weights[i], dims)
    end

    realintervals = calculate_hist_percentage(hists, hist)

    return reverse(hists), histogram, reverse(realintervals)
end

# for 1d histograms
function split_central(histogram, intervals)
    intervals = sort(intervals)
    intervals = (1 .-intervals)/2

    hist = deepcopy(histogram)
    hists = Array{StatsBase.Histogram}(undef, length(intervals))

    for i in 1:length(intervals)
        hists[i] = deepcopy(hist)
    end

    totalweights = sum(hist.weights)

    rel_weights = hist.weights/totalweights

    for (i, intv) in enumerate(intervals)
        sum_left = 0.
        sum_right = 0.

        for l in 1:length(hist.weights)
           if sum_left + rel_weights[l] < intv
                sum_left = sum_left + rel_weights[l]
                hists[i].weights[l] = 0
           else
                break
           end
        end

        for r in length(hist.weights):-1:1
           if sum_right + rel_weights[r] < intv
                sum_right = sum_right + rel_weights[r]
                hists[i].weights[r] = 0
           else
                break
           end
        end
    end

    realintervals = calculate_hist_percentage(hists, hist)

    return reverse(hists), histogram, reverse(realintervals)
end

function calculate_hist_percentage(hists, hist)
# calculate percentage really enclosed inside the intervals (possible differences for to large bins)
    totalweights = sum(hist.weights)
    realintervals = zeros(length(hists))

    for (i, h) in enumerate(hists)
        realintervals[i] = sum(h.weights)/totalweights
    end

    return realintervals
end

function bin_centers(h)
    m = length(h.edges[1])-1
    n = length(h.edges[2])-1

    x = Vector{Real}(undef, m)
    y = Vector{Real}(undef, n)

    for i in 1:m
        x[i] = h.edges[1][i]+0.5*(h.edges[1][i+1]-h.edges[1][i])
    end

    for j in 1:n
        y[j] =  h.edges[2][j]+0.5*(h.edges[2][j+1]-h.edges[2][j])
    end

    return x, y, h.weights
end

function calculate_levels(weights, intervals)
    intervals = sort(intervals)
    levels = Vector{Real}(undef, length(intervals)+1)

    w = vec(weights)
    w = sort(w, rev=true)

    weight_ids = sortperm(w, rev=true);
    summedw = sum(w)

    sumw = 0.0

    for w_id in weight_ids
        if(sumw <= intervals[end])
            i = findfirst(x -> sumw <= x, intervals)
            levels[i] = w[w_id]
            sumw += w[w_id]/summedw
        end
    end

    levels[end] = 1.1*summedw

    return sort(levels)
end


text:
# Julia wrapper for header: /usr/local/include/lz4.h
# docstrings copied from /usr/local/include/lz4.h

function check_compression_error(ret::Integer, function_name::AbstractString)
    ret == 0 && throw(LZ4Exception(function_name, "Compression failed."))
    return ret
end

function check_decompression_error(ret::Integer, function_name::AbstractString)
    ret < 0 && throw(LZ4Exception(function_name, "Decompression failed."))
    return ret
end

function check_initialized(stream::Ptr{LZ4_stream_t})
    if stream == Ptr{LZ4_stream_t}(C_NULL)
        throw(LZ4Exception("LZ4_stream_t", "Uninitialized compression stream"))
    end
end

function check_initialized(stream::Ptr{LZ4_streamDecode_t})
    if stream == Ptr{LZ4_streamDecode_t}(C_NULL)
        throw(LZ4Exception("LZ4_streamDecode_t", "Uninitialized decompression stream"))
    end
end

"""
    LZ4_compress_fast(src, dst, srcsize, dstcapacity=1)

Compresses `srcsize` bytes from buffer `src`
into already allocated `dst` buffer of size `dstcapacity`.
Compression is guaranteed to succeed if `dstcapacity` >= LZ4_compressBound(srcsize).
It also runs faster, so it's a recommended setting.
If the function cannot compress `src` into a limited `dst` budget,
compression stops *immediately*, and the function result is zero.
As a consequence, `dst` content is not valid.

Also allows to select an `acceleration` factor.
The larger the acceleration value, the faster the algorithm, but also the lesser the compression.
It's a trade-off. It can be fine tuned, with each successive value providing roughly +~3% to speed.
Values <= 0 will be replaced by the default which is 1.

This function never writes outside `dst` buffer, nor read outside `source` buffer.
    srcsize : supported max value is LZ4_MAX_INPUT_VALUE
    dstcapacity : full or partial size of buffer `dst` (which must be already allocated)

Returns the number of bytes written into buffer `dst` (necessarily <= dstcapacity)
"""
function LZ4_compress_fast(src, dst, srcsize, dstcapacity, acceleration=1)
    ret = ccall((:LZ4_compress_fast, liblz4), Cint, (Ptr{UInt8}, Ptr{UInt8}, Cint, Cint, Cint), src, dst, srcsize, dstcapacity, acceleration)
    check_compression_error(ret, "LZ4_compress_fast")
end

"""
    LZ4_compress_destSize(src, dst, srcsize, dstcapacity)

Reverse the logic : compresses as much data as possible from `src` buffer
into already allocated buffer `dst` of size `dstcapacity`.
This function either compresses the entire `src` content into `dst` if it's large enough,
or fill `dst` buffer completely with as much data as possible from `src`.
    *srcsize : will be modified to indicate how many bytes where read from `src` to fill `dst`.
                  New value is necessarily <= old value.

Returns number of bytes written into `dst` (necessarily <= dstcapacity)
"""
function LZ4_compress_destSize(src, dst, srcsize, dstcapacity)
    ret = ccall((:LZ4_compress_destSize, liblz4), Cint, (Ptr{UInt8}, Ptr{UInt8}, Ptr{Cint}, Cint), src, dst, srcsize, dstcapacity)
    check_compression_error(ret, "LZ4_compress_destSize")
end

"""
    LZ4_compressBound(inputsize)

Provides the maximum size that LZ4 compression may output in a "worst case" scenario (input data not compressible)
This function is primarily useful for memory allocation purposes (destination buffer size).
Macro LZ4_COMPRESSBOUND() is also provided for compilation-time evaluation (stack memory allocation for example).
Note that LZ4_compress_default() compress faster when dest buffer size is >= LZ4_compressBound(srcsize)
    inputsize  : max supported value is LZ4_MAX_INPUT_SIZE

Returns maximum output size in a "worst case" scenario or 0, if input size is too large ( > LZ4_MAX_INPUT_SIZE)
"""
function LZ4_compressBound(inputsize)
    ccall((:LZ4_compressBound, liblz4), Cint, (Cint,), inputsize)
end

"""
    LZ4_createStream()

Will allocate and initialize an `LZ4_stream_t` structure.
"""
function LZ4_createStream()
    str = ccall((:LZ4_createStream, liblz4), Ptr{LZ4_stream_t}, ())
    if str == C_NULL # Could not allocate memory
        throw(OutOfMemoryError())
    end
    return str
end

"""
    LZ4_freeStream(streamptr::Ptr{LZ4_stream_t})

Releases memory allocated by LZ4_createStream.
"""
function LZ4_freeStream(streamptr::Ptr{LZ4_stream_t})
    ccall((:LZ4_freeStream, liblz4), Cint, (Ptr{LZ4_stream_t},), streamptr)
end

"""
    LZ4_resetStream(streamptr)

An LZ4_stream_t structure can be allocated once and re-used multiple times.
Use this function to start compressing a new stream.
"""
function LZ4_resetStream(streamptr::Ptr{LZ4_stream_t})
    check_initialized(streamptr)
    ccall((:LZ4_resetStream, liblz4), Cvoid, (Ptr{LZ4_stream_t},), streamptr)
end

"""
    LZ4_compress_fast_continue(streamptr, src, dst, srcsize, dstcapacity, acceleration)

Compress content into `src` using data from previously compressed blocks, improving compression ratio.
`dst` buffer must be already allocated.
If dstcapacity >= LZ4_compressBound(srcsize), compression is guaranteed to succeed, and runs faster.

Important : Up to 64KB of previously compressed data is assumed to remain present and unmodified in memory !
Special 1 : If input buffer is a double-buffer, it can have any size, including < 64 KB.
Special 2 : If input buffer is a ring-buffer, it can have any size, including < 64 KB.

Returns size of compressed block.
After an error, the stream status is invalid, it can only be reset or freed.
"""
function LZ4_compress_fast_continue(streamptr::Ptr{LZ4_stream_t}, src, dst, srcsize, dstcapacity, acceleration=1)
    check_initialized(streamptr)
    ret = ccall((:LZ4_compress_fast_continue, liblz4), Cint, (Ptr{LZ4_stream_t}, Ptr{UInt8}, Ptr{UInt8}, Cint, Cint, Cint), streamptr, src, dst, srcsize, dstcapacity, acceleration)
    check_compression_error(ret, "LZ4_compress_fast_continue")
end

"""
    LZ4_decompress_safe(src, dst, cmpsize, dstcapacity)

Decompresses `cmpsize` bytes from buffer `src` into `dst`.
If destination buffer is not large enough, decoding will stop and output an error.
If the source stream is detected malformed, the function will stop decoding and error.
This function is protected against buffer overflow exploits, including malicious data packets.
It never writes outside output buffer, nor reads outside input buffer.

cmpsize : is the exact complete size of the compressed block.
dstcapacity : is the size of destination buffer, which must be already allocated.

Returns the number of bytes decompressed into destination buffer (necessarily <= dstcapacity)
"""
function LZ4_decompress_safe(src, dst, cmpsize, dstcapacity)
    ret = ccall((:LZ4_decompress_safe, liblz4), Cint, (Ptr{UInt8}, Ptr{UInt8}, Cint, Cint), src, dst, cmpsize, dstcapacity)
    check_decompression_error(ret, "LZ4_decompress_safe")
end

"""
    LZ4_createStreamDecode()

These decoding functions work the same as
creation / destruction of streaming decompression tracking structure.
A tracking structure can be re-used multiple times sequentially.
"""
function LZ4_createStreamDecode()
    str = ccall((:LZ4_createStreamDecode, liblz4), Ptr{LZ4_streamDecode_t}, ())
    if str == C_NULL # Could not allocate memory
        throw(OutOfMemoryError())
    end
    return str
end

"""
    LZ4_freeStreamDecode(streamptr)

These decoding functions work the same as
creation / destruction of streaming decompression tracking structure.
A tracking structure can be re-used multiple times sequentially.
"""
function LZ4_freeStreamDecode(streamptr::Ptr{LZ4_streamDecode_t})
    ccall((:LZ4_freeStreamDecode, liblz4), Cint, (Ptr{LZ4_streamDecode_t},), streamptr)
end

"""
    LZ4_setStreamDecode(streamptr, dictionary, dictsize)

An LZ4_streamDecode_t structure can be allocated once and re-used multiple times.
Use this function to start decompression of a new stream of blocks.
A dictionary can optionally be set. Use NULL or size 0 for a simple reset order.
Returns 1 if OK
"""
function LZ4_setStreamDecode(streamptr::Ptr{LZ4_streamDecode_t}, dictionary=C_NULL, dictsize=0)
    check_initialized(streamptr)
    ret = ccall((:LZ4_setStreamDecode, liblz4), Cint, (Ptr{LZ4_streamDecode_t}, Ptr{UInt8}, Cint), streamptr, dictionary, dictsize)
    ret == 0 && throw(LZ4Exception("LZ4_setStreamDecode", "Stream reinitialization failed."))
    return ret
end

"""
    LZ4_decompress_safe_continue(streamptr, src, dst, srcsize, dstcapacity)

These decoding functions allow decompression of consecutive blocks in "streaming" mode.
A block is an unsplittable entity, it must be presented entirely to a decompression function.
Decompression functions only accept one block at a time.
Previously decoded blocks *must* remain available at the memory position where they were decoded (up to 64 KB).

Special : if application sets a ring buffer for decompression, it must respect one of the following conditions :
- Exactly same size as encoding buffer, with same update rule (block boundaries at same positions)
  In which case, the decoding & encoding ring buffer can have any size, including very small ones ( < 64 KB).
- Larger than encoding buffer, by a minimum of maxBlockSize more bytes.
  maxBlockSize is implementation dependent. It's the maximum size of any single block.
  In which case, encoding and decoding buffers do not need to be synchronized,
  and encoding ring buffer can have any size, including small ones ( < 64 KB).
- _At least_ 64 KB + 8 bytes + maxBlockSize.
  In which case, encoding and decoding buffers do not need to be synchronized,
  and encoding ring buffer can have any size, including larger than decoding buffer.
Whenever these conditions are not possible, save the last 64KB of decoded data into a safe buffer,
and indicate where it is saved using `LZ4_setStreamDecode()` before decompressing next block.
"""
function LZ4_decompress_safe_continue(streamptr::Ptr{LZ4_streamDecode_t}, src, dst, srcsize, dstcapacity)
    check_initialized(streamptr)
    ret = ccall((:LZ4_decompress_safe_continue, liblz4), Cint, (Ptr{LZ4_streamDecode_t}, Ptr{UInt8}, Ptr{UInt8}, Cint, Cint), streamptr, src, dst, srcsize, dstcapacity)
    check_decompression_error(ret, "LZ4_decompress_safe_continue")
end


text:
immutable xc_func_info_type
    number::Cint
    kind::Cint
    name::Cstring
    familly::Cint
    refs::Cstring
    flags::Cint
    min_dens::Cdouble
    min_grad::Cdouble
    min_tau::Cdouble
    min_zeta::Cdouble
    init::Ptr{Void}
    end_::Ptr{Void}
    lda::Ptr{Void}
    gga::Ptr{Void}
    mgga::Ptr{Void}
end

immutable xc_func_type
    info::Ptr{xc_func_info_type}
    nspin::Cint
    n_func_aux::Cint
    func_aux::Ptr{xc_func_type}
    mix_coef::Ptr{Cdouble}
    cam_omega::Cdouble
    cam_alpha::Cdouble
    cam_beta::Cdouble
    func::Cint
    n_rho::Cint
    n_sigma::Cint
    n_tau::Cint
    n_lapl::Cint
    n_zk::Cint
    n_vrho::Cint
    n_vsigma::Cint
    n_vtau::Cint
    n_vlapl::Cint
    n_v2rho2::Cint
    n_v2sigma2::Cint
    n_v2tau2::Cint
    n_v2lapl2::Cint
    n_v2rhosigma::Cint
    n_v2rhotau::Cint
    n_v2rholapl::Cint
    n_v2sigmatau::Cint
    n_v2sigmalapl::Cint
    n_v2lapltau::Cint
    n_v3rho3::Cint
    n_v3rho2sigma::Cint
    n_v3rhosigma2::Cint
    n_v3sigma3::Cint
    params::Ptr{Void}
    xc_func_type() = new()
end

const XC_UNPOLARIZED          = 1
const XC_POLARIZED            = 2
const XC_LDA_X                = 1   #  Exchange
const XC_LDA_C_PW_MOD         = 13  #  Perdew & Wang (Modified)

type XC
    __xc__::Ref{xc_func_type}

    function XC(id::Int, spin::Int)
        __xc__ = Ref{xc_func_type}()
        ccall((:xc_func_init,libxc), Void, (Ref{xc_func_type}, Cint, Cint), __xc__, id, spin)
        xc = new(__xc__)
        finalizer(xc, xc -> ccall((:xc_func_end,libxc), Void, (Ref{xc_func_type},), xc.__xc__))
        xc
    end
end

function lda_pot!(xc::XC, rho::AbstractArray{Float64}, v::AbstractArray{Float64})
    @assert length(rho) == length(v)
    ccall((:xc_lda_vxc,libxc), Void,
        (Ref{xc_func_type}, Cint, Ptr{Cdouble}, Ptr{Cdouble}),
        xc.__xc__, length(rho), rho, v)
end

function lda_pot(xc::XC, rho::AbstractArray{Float64})
    v = similar(rho)
    lda_pot!(xc, rho, v)
    v
end

function lda_en_density!(xc::XC, rho::AbstractArray{Float64}, ε::AbstractArray{Float64})
    @assert length(rho) == length(ε)
    ccall((:xc_lda_exc,libxc), Void,
    	(Ref{xc_func_type}, Cint, Ptr{Cdouble}, Ptr{Cdouble}),
    	xc.__xc__, length(rho), rho, ε)
end

function lda_en_density(xc::XC, rho::AbstractArray{Float64})
    ε = similar(rho)
    lda_en_density!(xc, rho, ε)
    ε
end

function lda_en_density_pot!(xc::XC, rho::AbstractArray{Float64},
		ε::AbstractArray{Float64}, v::AbstractArray{Float64})
    @assert length(rho) == length(ε) == length(v)
    ccall((:xc_lda_exc_vxc,libxc), Void, (Ref{xc_func_type}, Cint, Ptr{Cdouble},
    	Ptr{Cdouble}, Ptr{Cdouble}), xc.__xc__, length(rho), rho, ε, v)
end

immutable ldav
    __xc__::XC
end

immutable ldae
    __xc__::XC
end

Base.call(xc::ldav,rho,v) = lda_pot!(xc.__xc__,rho,v)
Base.call(xc::ldav,rho) = lda_pot(xc.__xc__,rho)
Base.call(xc::ldae,rho,e) = lda_en_density!(xc.__xc__,rho,e)
Base.call(xc::ldae,rho) = lda_en_density(xc.__xc__,rho)

immutable lda
	__xc__::XC
    v::ldav
    e::ldae
end

function lda(id,spin)
	__xc__ = XC(id,spin)
	lda(__xc__,ldav(__xc__),ldae(__xc__))
end

Base.call(xc::lda,rho,e,v) = lda_en_density_pot!(xc.__xc__,rho,e,v)

doc"""
`ldaxc` is the pw-lda exchange-correlation functional [1,2].

calls the functionals `XC_LDA_X`, `XC_LDA_C_PW_MOD` in `libxc`.

spin-polarized and spin-unpolarized versions
can be created via `ldaxcpol()` and `ldaxcunpol()`

e.g.
```julia
xc = ldaxcunpol()

vx = xc.x.v(rho)  # exchange potential
ex = xc.x.e(rho)  # exchange energy density
vc = xc.c.v(rho)  # correlation potential
ec = xc.c.e(rho)  # correlation energy density

# in-place versions
xc.x.v(rho,vx)
xc.x.e(rho,ex)
xc.c.v(rho,vc)
xc.c.e(rho,ec)

# more efficient to compute both together
xc.x(rho,ex,vx)
xc.c(rho,ec,vc)

# compute all components
xc(rho,ex,vx,ec,vc)
```

[1] J. P. Perdew and Y. Wang. Accurate and simple analytic representation
of the electron-gas correlation energy. Phys. Rev. B, 45:13244–13249, 1992.

[2] D. M. Ceperley and B. J. Alder. Ground state of the electron gas by
a stochastic method. Phys. Rev. Lett., 45:566–569, 1980.
"""
immutable ldaxc
    x::lda
    c::lda
end

function Base.call(xc::ldaxc,rho,ex,vx,ec,vc)
	xc.x(rho,ex,vx)
	xc.c(rho,ec,vc)
end

ldaxunpol() = lda(XC_LDA_X,XC_UNPOLARIZED)
ldaxpol() = lda(XC_LDA_X,XC_POLARIZED)
ldacunpol() = lda(XC_LDA_C_PW_MOD,XC_UNPOLARIZED)
ldacpol() = lda(XC_LDA_C_PW_MOD,XC_POLARIZED)

ldaxcunpol() = ldaxc(ldaxunpol(),ldacunpol())
ldaxcpol() = ldaxc(ldaxpol(),ldacpol())

Base.show(io::IO, xc::ldav) = print(io, summary(xc))
Base.show(io::IO, xc::ldae) = print(io, summary(xc))
Base.show(io::IO, xc::lda) = print(io, summary(xc))
Base.show(io::IO, xc::ldaxc) = print(io, summary(xc))
Base.show(io::IO, xc::XC) = print(io, summary(xc))

# const XC_UNPOLARIZED          = 1
# const XC_POLARIZED            = 2

# const XC_NON_RELATIVISTIC     = 0
# const XC_RELATIVISTIC         = 1

# const XC_EXCHANGE             = 0
# const XC_CORRELATION          = 1
# const XC_EXCHANGE_CORRELATION = 2
# const XC_KINETIC              = 3

# const XC_FAMILY_UNKNOWN       = -1
# const XC_FAMILY_LDA           = 1
# const XC_FAMILY_GGA           = 2
# const XC_FAMILY_MGGA          = 4
# const XC_FAMILY_LCA           = 8
# const XC_FAMILY_OEP           = 16
# const XC_FAMILY_HYB_GGA       = 32
# const XC_FAMILY_HYB_MGGA      = 64

# # flags that can be used in info.flags
# const XC_FLAGS_HAVE_EXC       = ( 1 <<  0)  #    1
# const XC_FLAGS_HAVE_VXC       = ( 1 <<  1)  #    2
# const XC_FLAGS_HAVE_FXC       = ( 1 <<  2)  #    4
# const XC_FLAGS_HAVE_KXC       = ( 1 <<  3)  #    8
# const XC_FLAGS_HAVE_LXC       = ( 1 <<  4)  #   16
# const XC_FLAGS_1D             = ( 1 <<  5)  #   32
# const XC_FLAGS_2D             = ( 1 <<  6)  #   64
# const XC_FLAGS_3D             = ( 1 <<  7)  #  128
# const XC_FLAGS_HYB_CAM        = ( 1 <<  8)  #  256
# const XC_FLAGS_STABLE         = ( 1 <<  9)  #  512
# const XC_FLAGS_DEVELOPMENT    = ( 1 << 10)  # 1024

# const XC_TAU_EXPLICIT         = 0
# const XC_TAU_EXPANSION        = 1

# # These are old names keep for compatibility, and that should disappear soon
# const XC_GGA_XC_LB            = 160
# const XC_GGA_K_ABSR1          = 506
# const XC_GGA_K_ABSR2          = 507

# # from xc_funcs.h

# const XC_LDA_X                       = 1  #  Exchange
# const XC_LDA_C_WIGNER                = 2  #  Wigner parametrization
# const XC_LDA_C_RPA                   = 3  #  Random Phase Approximation
# const XC_LDA_C_HL                    = 4  #  Hedin & Lundqvist
# const XC_LDA_C_GL                    = 5  #  Gunnarson & Lundqvist
# const XC_LDA_C_XALPHA                = 6  #  Slater Xalpha
# const XC_LDA_C_VWN                   = 7  #  Vosko, Wilk, & Nussair (5
# const XC_LDA_C_VWN_RPA               = 8  #  Vosko, Wilk, & Nussair (RPA
# const XC_LDA_C_PZ                    = 9  #  Perdew & Zunger
# const XC_LDA_C_PZ_MOD                = 10  #  Perdew & Zunger (Modified
# const XC_LDA_C_OB_PZ                 = 11  #  Ortiz & Ballone (PZ
# const XC_LDA_C_PW                    = 12  #  Perdew & Wang
# const XC_LDA_C_PW_MOD                = 13  #  Perdew & Wang (Modified
# const XC_LDA_C_OB_PW                 = 14  #  Ortiz & Ballone (PW
# const XC_LDA_C_2D_AMGB               = 15  #  Attaccalite et al
# const XC_LDA_C_2D_PRM                = 16  #  Pittalis, Rasanen & Marques correlation in 2D
# const XC_LDA_C_vBH                   = 17  #  von Barth & Hedin
# const XC_LDA_C_1D_CSC                = 18  #  Casula, Sorella, and Senatore 1D correlation
# const XC_LDA_X_2D                    = 19  #  Exchange in 2D
# const XC_LDA_XC_TETER93              = 20  #  Teter 93 parametrization
# const XC_LDA_X_1D                    = 21  #  Exchange in 1D
# const XC_LDA_C_ML1                   = 22  #  Modified LSD (version 1) of Proynov and Salahub
# const XC_LDA_C_ML2                   = 23  #  Modified LSD (version 2) of Proynov and Salahub
# const XC_LDA_C_GOMBAS                = 24  #  Gombas parametrization
# const XC_LDA_C_PW_RPA                = 25  #  Perdew & Wang fit of the RPA
# const XC_LDA_C_1D_LOOS               = 26  #  P-F Loos correlation LDA
# const XC_LDA_C_RC04                  = 27  #  Ragot-Cortona
# const XC_LDA_C_VWN_1                 = 28  #  Vosko, Wilk, & Nussair (1
# const XC_LDA_C_VWN_2                 = 29  #  Vosko, Wilk, & Nussair (2
# const XC_LDA_C_VWN_3                 = 30  #  Vosko, Wilk, & Nussair (3
# const XC_LDA_C_VWN_4                 = 31  #  Vosko, Wilk, & Nussair (4
# const XC_LDA_K_TF                    = 50  #  Thomas-Fermi kinetic energy functional
# const XC_LDA_K_LP                    = 51  #  Lee and Parr Gaussian ansatz
# const XC_GGA_C_Q2D                   = 47  #  Chiodo et al
# const XC_GGA_X_Q2D                   = 48  #  Chiodo et al
# const XC_GGA_X_PBE_MOL               = 49  #  Del Campo, Gazquez, Trickey and Vela (PBE-like
# const XC_GGA_K_TFVW                  = 52  #  Thomas-Fermi plus von Weiszaecker correction
# const XC_GGA_K_REVAPBEINT            = 53  #  interpolated version of REVAPBE
# const XC_GGA_K_APBEINT               = 54  #  interpolated version of APBE
# const XC_GGA_K_REVAPBE               = 55  #  revised APBE
# const XC_GGA_X_AK13                  = 56  #  Armiento & Kuemmel 2013
# const XC_GGA_K_MEYER                 = 57  #  Meyer,  Wang, and Young
# const XC_GGA_X_LV_RPW86              = 58  #  Berland and Hyldgaard
# const XC_GGA_X_PBE_TCA               = 59  #  PBE revised by Tognetti et al
# const XC_GGA_X_PBEINT                = 60  #  PBE for hybrid interfaces
# const XC_GGA_C_ZPBEINT               = 61  #  spin-dependent gradient correction to PBEint
# const XC_GGA_C_PBEINT                = 62  #  PBE for hybrid interfaces
# const XC_GGA_C_ZPBESOL               = 63  #  spin-dependent gradient correction to PBEsol
# const XC_GGA_XC_OPBE_D               = 65  #  oPBE_D functional of Goerigk and Grimme
# const XC_GGA_XC_OPWLYP_D             = 66  #  oPWLYP-D functional of Goerigk and Grimme
# const XC_GGA_XC_OBLYP_D              = 67  #  oBLYP-D functional of Goerigk and Grimme
# const XC_GGA_X_VMT84_GE              = 68  #  VMT{8,4} with constraint satisfaction with mu           = mu_GE
# const XC_GGA_X_VMT84_PBE             = 69  #  VMT{8,4} with constraint satisfaction with mu           = mu_PBE
# const XC_GGA_X_VMT_GE                = 70  #  Vela, Medel, and Trickey with mu           = mu_GE
# const XC_GGA_X_VMT_PBE               = 71  #  Vela, Medel, and Trickey with mu           = mu_PBE
# const XC_GGA_C_N12_SX                = 79  #  N12-SX functional from Minnesota
# const XC_GGA_C_N12                   = 80  #  N12 functional from Minnesota
# const XC_GGA_X_N12                   = 82  #  N12 functional from Minnesota
# const XC_GGA_C_VPBE                  = 83  #  variant PBE
# const XC_GGA_C_OP_XALPHA             = 84  #  one-parameter progressive functional (XALPHA version
# const XC_GGA_C_OP_G96                = 85  #  one-parameter progressive functional (G96 version
# const XC_GGA_C_OP_PBE                = 86  #  one-parameter progressive functional (PBE version
# const XC_GGA_C_OP_B88                = 87  #  one-parameter progressive functional (B88 version
# const XC_GGA_C_FT97                  = 88  #  Filatov & Thiel correlation
# const XC_GGA_C_SPBE                  = 89  #  PBE correlation to be used with the SSB exchange
# const XC_GGA_X_SSB_SW                = 90  #  Swarta, Sola and Bickelhaupt correction to PBE
# const XC_GGA_X_SSB                   = 91  #  Swarta, Sola and Bickelhaupt
# const XC_GGA_X_SSB_D                 = 92  #  Swarta, Sola and Bickelhaupt dispersion
# const XC_GGA_XC_HCTH_407P            = 93  #  HCTH/407
# const XC_GGA_XC_HCTH_P76             = 94  #  HCTH p          =7/6
# const XC_GGA_XC_HCTH_P14             = 95  #  HCTH p          =1/4
# const XC_GGA_XC_B97_GGA1             = 96  #  Becke 97 GGA-1
# const XC_GGA_XC_HCTH_A               = 97  #  HCTH-A
# const XC_GGA_X_BPCCAC                = 98  #  BPCCAC (GRAC for the energy
# const XC_GGA_C_REVTCA                = 99  #  Tognetti, Cortona, Adamo (revised
# const XC_GGA_C_TCA                   = 100  #  Tognetti, Cortona, Adamo
# const XC_GGA_X_PBE                   = 101  #  Perdew, Burke & Ernzerhof exchange
# const XC_GGA_X_PBE_R                 = 102  #  Perdew, Burke & Ernzerhof exchange (revised
# const XC_GGA_X_B86                   = 103  #  Becke 86 Xalfa,beta,gamma
# const XC_GGA_X_HERMAN                = 104  #  Herman et al original GGA
# const XC_GGA_X_B86_MGC               = 105  #  Becke 86 Xalfa,beta,gamma (with mod. grad. correction
# const XC_GGA_X_B88                   = 106  #  Becke 88
# const XC_GGA_X_G96                   = 107  #  Gill 96
# const XC_GGA_X_PW86                  = 108  #  Perdew & Wang 86
# const XC_GGA_X_PW91                  = 109  #  Perdew & Wang 91
# const XC_GGA_X_OPTX                  = 110  #  Handy & Cohen OPTX 01
# const XC_GGA_X_DK87_R1               = 111  #  dePristo & Kress 87 (version R1
# const XC_GGA_X_DK87_R2               = 112  #  dePristo & Kress 87 (version R2
# const XC_GGA_X_LG93                  = 113  #  Lacks & Gordon 93
# const XC_GGA_X_FT97_A                = 114  #  Filatov & Thiel 97 (version A
# const XC_GGA_X_FT97_B                = 115  #  Filatov & Thiel 97 (version B
# const XC_GGA_X_PBE_SOL               = 116  #  Perdew, Burke & Ernzerhof exchange (solids
# const XC_GGA_X_RPBE                  = 117  #  Hammer, Hansen & Norskov (PBE-like
# const XC_GGA_X_WC                    = 118  #  Wu & Cohen
# const XC_GGA_X_MPW91                 = 119  #  Modified form of PW91 by Adamo & Barone
# const XC_GGA_X_AM05                  = 120  #  Armiento & Mattsson 05 exchange
# const XC_GGA_X_PBEA                  = 121  #  Madsen (PBE-like
# const XC_GGA_X_MPBE                  = 122  #  Adamo & Barone modification to PBE
# const XC_GGA_X_XPBE                  = 123  #  xPBE reparametrization by Xu & Goddard
# const XC_GGA_X_2D_B86_MGC            = 124  #  Becke 86 MGC for 2D systems
# const XC_GGA_X_BAYESIAN              = 125  #  Bayesian best fit for the enhancement factor
# const XC_GGA_X_PBE_JSJR              = 126  #  JSJR reparametrization by Pedroza, Silva & Capelle
# const XC_GGA_X_2D_B88                = 127  #  Becke 88 in 2D
# const XC_GGA_X_2D_B86                = 128  #  Becke 86 Xalfa,beta,gamma
# const XC_GGA_X_2D_PBE                = 129  #  Perdew, Burke & Ernzerhof exchange in 2D
# const XC_GGA_C_PBE                   = 130  #  Perdew, Burke & Ernzerhof correlation
# const XC_GGA_C_LYP                   = 131  #  Lee, Yang & Parr
# const XC_GGA_C_P86                   = 132  #  Perdew 86
# const XC_GGA_C_PBE_SOL               = 133  #  Perdew, Burke & Ernzerhof correlation SOL
# const XC_GGA_C_PW91                  = 134  #  Perdew & Wang 91
# const XC_GGA_C_AM05                  = 135  #  Armiento & Mattsson 05 correlation
# const XC_GGA_C_XPBE                  = 136  #  xPBE reparametrization by Xu & Goddard
# const XC_GGA_C_LM                    = 137  #  Langreth and Mehl correlation
# const XC_GGA_C_PBE_JRGX              = 138  #  JRGX reparametrization by Pedroza, Silva & Capelle
# const XC_GGA_X_OPTB88_VDW            = 139  #  Becke 88 reoptimized to be used with vdW functional of Dion et al
# const XC_GGA_X_PBEK1_VDW             = 140  #  PBE reparametrization for vdW
# const XC_GGA_X_OPTPBE_VDW            = 141  #  PBE reparametrization for vdW
# const XC_GGA_X_RGE2                  = 142  #  Regularized PBE
# const XC_GGA_C_RGE2                  = 143  #  Regularized PBE
# const XC_GGA_X_RPW86                 = 144  #  refitted Perdew & Wang 86
# const XC_GGA_X_KT1                   = 145  #  Keal and Tozer version 1
# const XC_GGA_XC_KT2                  = 146  #  Keal and Tozer version 2
# const XC_GGA_C_WL                    = 147  #  Wilson & Levy
# const XC_GGA_C_WI                    = 148  #  Wilson & Ivanov
# const XC_GGA_X_MB88                  = 149  #  Modified Becke 88 for proton transfer
# const XC_GGA_X_SOGGA                 = 150  #  Second-order generalized gradient approximation
# const XC_GGA_X_SOGGA11               = 151  #  Second-order generalized gradient approximation 2011
# const XC_GGA_C_SOGGA11               = 152  #  Second-order generalized gradient approximation 2011
# const XC_GGA_C_WI0                   = 153  #  Wilson & Ivanov initial version
# const XC_GGA_XC_TH1                  = 154  #  Tozer and Handy v. 1
# const XC_GGA_XC_TH2                  = 155  #  Tozer and Handy v. 2
# const XC_GGA_XC_TH3                  = 156  #  Tozer and Handy v. 3
# const XC_GGA_XC_TH4                  = 157  #  Tozer and Handy v. 4
# const XC_GGA_X_C09X                  = 158  #  C09x to be used with the VdW of Rutgers-Chalmers
# const XC_GGA_C_SOGGA11_X             = 159  #  To be used with hyb_gga_x_SOGGA11-X
# const XC_GGA_X_LB                    = 160  #  van Leeuwen & Baerends
# const XC_GGA_XC_HCTH_93              = 161  #  HCTH functional fitted to  93 molecules
# const XC_GGA_XC_HCTH_120             = 162  #  HCTH functional fitted to 120 molecules
# const XC_GGA_XC_HCTH_147             = 163  #  HCTH functional fitted to 147 molecules
# const XC_GGA_XC_HCTH_407             = 164  #  HCTH functional fitted to 407 molecules
# const XC_GGA_XC_EDF1                 = 165  #  Empirical functionals from Adamson, Gill, and Pople
# const XC_GGA_XC_XLYP                 = 166  #  XLYP functional
# const XC_GGA_XC_B97                  = 167  #  Becke 97
# const XC_GGA_XC_B97_1                = 168  #  Becke 97-1
# const XC_GGA_XC_B97_2                = 169  #  Becke 97-2
# const XC_GGA_XC_B97_D                = 170  #  Grimme functional to be used with C6 vdW term
# const XC_GGA_XC_B97_K                = 171  #  Boese-Martin for Kinetics
# const XC_GGA_XC_B97_3                = 172  #  Becke 97-3
# const XC_GGA_XC_PBE1W                = 173  #  Functionals fitted for water
# const XC_GGA_XC_MPWLYP1W             = 174  #  Functionals fitted for water
# const XC_GGA_XC_PBELYP1W             = 175  #  Functionals fitted for water
# const XC_GGA_XC_SB98_1a              = 176  #  Schmider-Becke 98 parameterization 1a
# const XC_GGA_XC_SB98_1b              = 177  #  Schmider-Becke 98 parameterization 1b
# const XC_GGA_XC_SB98_1c              = 178  #  Schmider-Becke 98 parameterization 1c
# const XC_GGA_XC_SB98_2a              = 179  #  Schmider-Becke 98 parameterization 2a
# const XC_GGA_XC_SB98_2b              = 180  #  Schmider-Becke 98 parameterization 2b
# const XC_GGA_XC_SB98_2c              = 181  #  Schmider-Becke 98 parameterization 2c
# const XC_GGA_X_LBM                   = 182  #  van Leeuwen & Baerends modified
# const XC_GGA_X_OL2                   = 183  #  Exchange form based on Ou-Yang and Levy v.2
# const XC_GGA_X_APBE                  = 184  #  mu fixed from the semiclassical neutral atom
# const XC_GGA_K_APBE                  = 185  #  mu fixed from the semiclassical neutral atom
# const XC_GGA_C_APBE                  = 186  #  mu fixed from the semiclassical neutral atom
# const XC_GGA_K_TW1                   = 187  #  Tran and Wesolowski set 1 (Table II
# const XC_GGA_K_TW2                   = 188  #  Tran and Wesolowski set 2 (Table II
# const XC_GGA_K_TW3                   = 189  #  Tran and Wesolowski set 3 (Table II
# const XC_GGA_K_TW4                   = 190  #  Tran and Wesolowski set 4 (Table II
# const XC_GGA_X_HTBS                  = 191  #  Haas, Tran, Blaha, and Schwarz
# const XC_GGA_X_AIRY                  = 192  #  Constantin et al based on the Airy gas
# const XC_GGA_X_LAG                   = 193  #  Local Airy Gas
# const XC_GGA_XC_MOHLYP               = 194  #  Functional for organometallic chemistry
# const XC_GGA_XC_MOHLYP2              = 195  #  Functional for barrier heights
# const XC_GGA_XC_TH_FL                = 196  #  Tozer and Handy v. FL
# const XC_GGA_XC_TH_FC                = 197  #  Tozer and Handy v. FC
# const XC_GGA_XC_TH_FCFO              = 198  #  Tozer and Handy v. FCFO
# const XC_GGA_XC_TH_FCO               = 199  #  Tozer and Handy v. FCO
# const XC_GGA_C_OPTC                  = 200  #  Optimized correlation functional of Cohen and Handy
# const XC_GGA_K_VW                    = 500  #  von Weiszaecker functional
# const XC_GGA_K_GE2                   = 501  #  Second-order gradient expansion (l           = 1/9
# const XC_GGA_K_GOLDEN                = 502  #  TF-lambda-vW form by Golden (l           = 13/45
# const XC_GGA_K_YT65                  = 503  #  TF-lambda-vW form by Yonei and Tomishima (l           = 1/5
# const XC_GGA_K_BALTIN                = 504  #  TF-lambda-vW form by Baltin (l           = 5/9
# const XC_GGA_K_LIEB                  = 505  #  TF-lambda-vW form by Lieb (l           = 0.185909191
# const XC_GGA_K_ABSP1                 = 506  #  gamma-TFvW form by Acharya et al [g           = 1 - 1.412/N^(1/3
# const XC_GGA_K_ABSP2                 = 507  #  gamma-TFvW form by Acharya et al [g           = 1 - 1.332/N^(1/3
# const XC_GGA_K_GR                    = 508  #  gamma-TFvW form by Gázquez and Robles
# const XC_GGA_K_LUDENA                = 509  #  gamma-TFvW form by Ludeña
# const XC_GGA_K_GP85                  = 510  #  gamma-TFvW form by Ghosh and Parr
# const XC_GGA_K_PEARSON               = 511  #  Pearson
# const XC_GGA_K_OL1                   = 512  #  Ou-Yang and Levy v.1
# const XC_GGA_K_OL2                   = 513  #  Ou-Yang and Levy v.2
# const XC_GGA_K_FR_B88                = 514  #  Fuentealba & Reyes (B88 version
# const XC_GGA_K_FR_PW86               = 515  #  Fuentealba & Reyes (PW86 version
# const XC_GGA_K_DK                    = 516  #  DePristo and Kress
# const XC_GGA_K_PERDEW                = 517  #  Perdew
# const XC_GGA_K_VSK                   = 518  #  Vitos, Skriver, and Kollar
# const XC_GGA_K_VJKS                  = 519  #  Vitos, Johansson, Kollar, and Skriver
# const XC_GGA_K_ERNZERHOF             = 520  #  Ernzerhof
# const XC_GGA_K_LC94                  = 521  #  Lembarki & Chermette
# const XC_GGA_K_LLP                   = 522  #  Lee, Lee & Parr
# const XC_GGA_K_THAKKAR               = 523  #  Thakkar 1992
# const XC_GGA_X_WPBEH                 = 524  #  short-range version of the PBE
# const XC_GGA_X_HJS_PBE               = 525  #  HJS screened exchange PBE version
# const XC_GGA_X_HJS_PBE_SOL           = 526  #  HJS screened exchange PBE_SOL version
# const XC_GGA_X_HJS_B88               = 527  #  HJS screened exchange B88 version
# const XC_GGA_X_HJS_B97X              = 528  #  HJS screened exchange B97x version
# const XC_GGA_X_ITYH                  = 529  #  short-range recipe for exchange GGA functionals
# const XC_GGA_X_SFAT                  = 530  #  short-range recipe for exchange GGA functionals
# const XC_HYB_GGA_X_N12_SX            = 81  #  N12-SX functional from Minnesota
# const XC_HYB_GGA_XC_B3PW91           = 401  #  The original (ACM) hybrid of Becke
# const XC_HYB_GGA_XC_B3LYP            = 402  #  The (in)famous B3LYP
# const XC_HYB_GGA_XC_B3P86            = 403  #  Perdew 86 hybrid similar to B3PW91
# const XC_HYB_GGA_XC_O3LYP            = 404  #  hybrid using the optx functional
# const XC_HYB_GGA_XC_mPW1K            = 405  #  mixture of mPW91 and PW91 optimized for kinetics
# const XC_HYB_GGA_XC_PBEH             = 406  #  aka PBE0 or PBE1PBE
# const XC_HYB_GGA_XC_B97              = 407  #  Becke 97
# const XC_HYB_GGA_XC_B97_1            = 408  #  Becke 97-1
# const XC_HYB_GGA_XC_B97_2            = 410  #  Becke 97-2
# const XC_HYB_GGA_XC_X3LYP            = 411  #  maybe the best hybrid
# const XC_HYB_GGA_XC_B1WC             = 412  #  Becke 1-parameter mixture of WC and PBE
# const XC_HYB_GGA_XC_B97_K            = 413  #  Boese-Martin for Kinetics
# const XC_HYB_GGA_XC_B97_3            = 414  #  Becke 97-3
# const XC_HYB_GGA_XC_MPW3PW           = 415  #  mixture with the mPW functional
# const XC_HYB_GGA_XC_B1LYP            = 416  #  Becke 1-parameter mixture of B88 and LYP
# const XC_HYB_GGA_XC_B1PW91           = 417  #  Becke 1-parameter mixture of B88 and PW91
# const XC_HYB_GGA_XC_mPW1PW           = 418  #  Becke 1-parameter mixture of mPW91 and PW91
# const XC_HYB_GGA_XC_MPW3LYP          = 419  #  mixture of mPW and LYP
# const XC_HYB_GGA_XC_SB98_1a          = 420  #  Schmider-Becke 98 parameterization 1a
# const XC_HYB_GGA_XC_SB98_1b          = 421  #  Schmider-Becke 98 parameterization 1b
# const XC_HYB_GGA_XC_SB98_1c          = 422  #  Schmider-Becke 98 parameterization 1c
# const XC_HYB_GGA_XC_SB98_2a          = 423  #  Schmider-Becke 98 parameterization 2a
# const XC_HYB_GGA_XC_SB98_2b          = 424  #  Schmider-Becke 98 parameterization 2b
# const XC_HYB_GGA_XC_SB98_2c          = 425  #  Schmider-Becke 98 parameterization 2c
# const XC_HYB_GGA_X_SOGGA11_X         = 426  #  Hybrid based on SOGGA11 form
# const XC_HYB_GGA_XC_HSE03            = 427  #  the 2003 version of the screened hybrid HSE
# const XC_HYB_GGA_XC_HSE06            = 428  #  the 2006 version of the screened hybrid HSE
# const XC_HYB_GGA_XC_HJS_PBE          = 429  #  HJS hybrid screened exchange PBE version
# const XC_HYB_GGA_XC_HJS_PBE_SOL      = 430  #  HJS hybrid screened exchange PBE_SOL version
# const XC_HYB_GGA_XC_HJS_B88          = 431  #  HJS hybrid screened exchange B88 version
# const XC_HYB_GGA_XC_HJS_B97X         = 432  #  HJS hybrid screened exchange B97x version
# const XC_HYB_GGA_XC_CAM_B3LYP        = 433  #  CAM version of B3LYP
# const XC_HYB_GGA_XC_TUNED_CAM_B3LYP  = 434  #  CAM version of B3LYP tunes for excitations
# const XC_HYB_GGA_XC_BHANDH           = 435  #  Becke half-and-half
# const XC_HYB_GGA_XC_BHANDHLYP        = 436  #  Becke half-and-half with B88 exchange
# const XC_HYB_GGA_XC_MB3LYP_RC04      = 437  #  B3LYP with RC04 LDA
# const XC_HYB_GGA_XC_MPWLYP1M         = 453  #  MPW with 1 par. for metals/LYP
# const XC_HYB_GGA_XC_REVB3LYP         = 454  #  Revised B3LYP
# const XC_HYB_GGA_XC_CAMY_BLYP        = 455  #  BLYP with yukawa screening
# const XC_HYB_GGA_XC_PBE0_13          = 456  #  PBE0-1/3
# const XC_MGGA_XC_OTPSS_D             = 64  #  oTPSS_D functional of Goerigk and Grimme
# const XC_MGGA_C_CS                   = 72  #  Colle and Salvetti
# const XC_MGGA_C_MN12_SX              = 73  #  MN12-SX functional of Minnesota
# const XC_MGGA_C_MN12_L               = 74  #  MN12-L functional of Minnesota
# const XC_MGGA_C_M11_L                = 75  #  M11-L functional of Minnesota
# const XC_MGGA_C_M11                  = 76  #  M11 functional of Minnesota
# const XC_MGGA_C_M08_SO               = 77  #  M08-SO functional of Minnesota
# const XC_MGGA_C_M08_HX               = 78  #  M08-HX functional of Minnesota
# const XC_MGGA_X_LTA                  = 201  #  Local tau approximation of Ernzerhof & Scuseria
# const XC_MGGA_X_TPSS                 = 202  #  Perdew, Tao, Staroverov & Scuseria exchange
# const XC_MGGA_X_M06_L                = 203  #  M06-Local functional of Minnesota
# const XC_MGGA_X_GVT4                 = 204  #  GVT4 from Van Voorhis and Scuseria
# const XC_MGGA_X_TAU_HCTH             = 205  #  tau-HCTH from Boese and Handy
# const XC_MGGA_X_BR89                 = 206  #  Becke-Roussel 89
# const XC_MGGA_X_BJ06                 = 207  #  Becke & Johnson correction to Becke-Roussel 89
# const XC_MGGA_X_TB09                 = 208  #  Tran & Blaha correction to Becke & Johnson
# const XC_MGGA_X_RPP09                = 209  #  Rasanen, Pittalis, and Proetto correction to Becke & Johnson
# const XC_MGGA_X_2D_PRHG07            = 210  #  Pittalis, Rasanen, Helbig, Gross Exchange Functional
# const XC_MGGA_X_2D_PRHG07_PRP10      = 211  #  PRGH07 with PRP10 correction
# const XC_MGGA_X_REVTPSS              = 212  #  revised Perdew, Tao, Staroverov & Scuseria exchange
# const XC_MGGA_X_PKZB                 = 213  #  Perdew, Kurth, Zupan, and Blaha
# const XC_MGGA_X_M05                  = 214  #  M05 functional of Minnesota
# const XC_MGGA_X_M05_2X               = 215  #  M05-2X functional of Minnesota
# const XC_MGGA_X_M06_HF               = 216  #  M06-HF functional of Minnesota
# const XC_MGGA_X_M06                  = 217  #  M06 functional of Minnesota
# const XC_MGGA_X_M06_2X               = 218  #  M06-2X functional of Minnesota
# const XC_MGGA_X_M08_HX               = 219  #  M08-HX functional of Minnesota
# const XC_MGGA_X_M08_SO               = 220  #  M08-SO functional of Minnesota
# const XC_MGGA_X_MS0                  = 221  #  MS exchange of Sun, Xiao, and Ruzsinszky
# const XC_MGGA_X_MS1                  = 222  #  MS1 exchange of Sun, et al
# const XC_MGGA_X_MS2                  = 223  #  MS2 exchange of Sun, et al
# const XC_MGGA_X_MS2H                 = 224  #  MS2 hybrid exchange of Sun, et al
# const XC_MGGA_X_M11_L                = 226  #  M11-L functional of Minnesota
# const XC_MGGA_X_MN12_L               = 227  #  MN12-L functional from Minnesota
# const XC_MGGA_X_MN12_SX              = 228  #  MN12-SX functional from Minnesota
# const XC_MGGA_C_CC06                 = 229  #  Cancio and Chou 2006
# const XC_MGGA_X_MK00                 = 230  #  Exchange for accurate virtual orbital energies
# const XC_MGGA_C_TPSS                 = 231  #  Perdew, Tao, Staroverov & Scuseria correlation
# const XC_MGGA_C_VSXC                 = 232  #  VSxc from Van Voorhis and Scuseria (correlation part
# const XC_MGGA_C_M06_L                = 233  #  M06-Local functional of Minnesota
# const XC_MGGA_C_M06_HF               = 234  #  M06-HF functional of Minnesota
# const XC_MGGA_C_M06                  = 235  #  M06 functional of Minnesota
# const XC_MGGA_C_M06_2X               = 236  #  M06-2X functional of Minnesota
# const XC_MGGA_C_M05                  = 237  #  M05 functional of Minnesota
# const XC_MGGA_C_M05_2X               = 238  #  M05-2X functional of Minnesota
# const XC_MGGA_C_PKZB                 = 239  #  Perdew, Kurth, Zupan, and Blaha
# const XC_MGGA_C_BC95                 = 240  #  Becke correlation 95
# const XC_MGGA_C_REVTPSS              = 241  #  revised TPSS correlation
# const XC_MGGA_XC_TPSSLYP1W           = 242  #  Functionals fitted for water
# const XC_MGGA_X_MK00B                = 243  #  Exchange for accurate virtual orbital energies (v. B
# const XC_MGGA_X_BLOC                 = 244  #  functional with balanced localization
# const XC_MGGA_X_MODTPSS              = 245  #  Modified Perdew, Tao, Staroverov & Scuseria exchange
# const XC_HYB_MGGA_X_M11              = 225  #  M11 functional of Minnesota
# const XC_HYB_MGGA_XC_M05             = 438  #  M05 functional of Minnesota
# const XC_HYB_MGGA_XC_M05_2X          = 439  #  M05-2X functional of Minnesota
# const XC_HYB_MGGA_XC_B88B95          = 440  #  Mixture of B88 with BC95 (B1B95
# const XC_HYB_MGGA_XC_B86B95          = 441  #  Mixture of B86 with BC95
# const XC_HYB_MGGA_XC_PW86B95         = 442  #  Mixture of PW86 with BC95
# const XC_HYB_MGGA_XC_BB1K            = 443  #  Mixture of B88 with BC95 from Zhao and Truhlar
# const XC_HYB_MGGA_XC_M06_HF          = 444  #  M06-HF functional of Minnesota
# const XC_HYB_MGGA_XC_MPW1B95         = 445  #  Mixture of mPW91 with BC95 from Zhao and Truhlar
# const XC_HYB_MGGA_XC_MPWB1K          = 446  #  Mixture of mPW91 with BC95 for kinetics
# const XC_HYB_MGGA_XC_X1B95           = 447  #  Mixture of X with BC95
# const XC_HYB_MGGA_XC_XB1K            = 448  #  Mixture of X with BC95 for kinetics
# const XC_HYB_MGGA_XC_M06             = 449  #  M06 functional of Minnesota
# const XC_HYB_MGGA_XC_M06_2X          = 450  #  M06-2X functional of Minnesota
# const XC_HYB_MGGA_XC_PW6B95          = 451  #  Mixture of PW91 with BC95 from Zhao and Truhlar
# const XC_HYB_MGGA_XC_PWB6K           = 452  #  Mixture of PW91 with BC95 from Zhao and Truhlar for kinetics
# const XC_HYB_MGGA_XC_TPSSH           = 457  #     TPSS hybrid
# const XC_HYB_MGGA_XC_REVTPSSH        = 458  #  revTPSS hybrid

max_stars_repo_name: johncsnyder/Element


text:
using SymPy, Plots
gr()

#TODO consider looking into how Michal Kaveslov did path following -> discrete points (states), but how did he know which ones to choose??

# track data
type Track
  R
  x0
  y0
  Xn
  Yn
end

# EX4 test data -> https://www.varsitytutors.com/hotmath/hotmath_help/topics/shortest-distance-between-a-point-and-a-circle
function Track()
  Track(7,
        4,
       -5,
       -4,
       -11);
 end

function initTrack(;R=R,x0=x0,y0=y0,Xn=Xn,Yn=Yn)
  t=Track();
  t.R=R;
  t.x0=x0;
  t.y0=y0;
  t.Xn=Xn;
  t.Yn=Yn;
  return t
end

function minD(t::Track)
  sign(sqrt((t.Xn-t.x0)^2 + (t.Yn-t.y0)^2)-t.R)*(sqrt((t.Xn-t.x0)^2 + (t.Yn-t.y0)^2)-t.R);
end

function plotTrack(t::Track)
  pp=plot(0,leg=:false);
  title!(@sprintf("min D (m) = %0.2f",minD(t)))
  pts = Plots.partialcircle(0,2π,100,t.R);
  x, y = Plots.unzip(pts)
  x += t.x0;
  y += t.y0;
  pts = collect(zip(x, y));
  plot!(pts,line=(:path,5,:solid,:green),label="Track");
  XX=[t.Xn;t.x0];YY=[t.Yn;t.y0];
  plot!(XX,YY,line=(:path,5,:solid),label="D # ");

  pts = Plots.partialcircle(0,2π,100,0.25);
  x, y = Plots.unzip(pts);
  x += t.Xn;  y += t.Yn;
  pts = collect(zip(x, y));
  plot!(pts,c=:black,fill=:true,label="Vehicle")
  #xlabel!("X (m)");  ylabel!("Y (m)");
  #xlims!(-10,10); ylims!(-15,10); #TODO --> use this as a position plot
end

format=:png;

# ex4
t4=Track();
plotTrack(Track(t4));
savefig(string("t4.",format));

# ex3
t3=initTrack(;R=5,x0=-3,y0=3,Xn=-2,Yn=0);
plotTrack(t3);
savefig(string("t3.",format));

# real track
t1=initTrack(;R=100,x0=100,y0=100,Xn=0,Yn=0);
plotTrack(t1);
savefig(string("t1.",format));

# this script is designed to calcualte the analytical distance from a point to an elipse (if there is one...)
# not able to find the analytical solution

function elllipseD()
  a,b,x,y=symbols("a,b,x,y",positive=true);
  x0,y0,Xn,Yn=symbols("x0,y0,Xn,Yn",positive=true);
  D=symbols("D",positive=true);
  epsilon=symbols("epsilon",positive=true)
  # equation for distance
  exp1=(x-Xn)^2 + (y-Yn)^2 - D^2;

  # derivative of D w.r.t x == slope of D
  exp3=diff(exp1,x) - (y-Yn)/(x-Xn+sign(x-Xn)*epsilon);  # only Chuck Noris can divide by zero

  # equation for ellipse
  exp2=(((x-x0)/a)^2 + ((y-y0)/b)^2)-1;

  # both equations
  exp = [exp1, exp2, exp3];

  # solve the system
  solve(exp,[D,x])

  ####################
  # test the solution
  ####################
  # test data - > all variables and parameters must be positive
  A_=5; B_=10;
  x0_=5; y0_=20;
  Xn_=10; Yn_=0;

  # find the nearest point
  solve(exp1,)# solve for x

  pp=plot(0,leg=:false);

  A=5; B=10;
  x0=5; y0=20;
  Xn=10; Yn=0;

  pts = Plots.partialcircle(0,2π,100,A)
  x, y = Plots.unzip(pts)
  x += x0;
  y = B/A*y + y0;
  #pts = collect(zip(x, y))
  plot!(Shape(x,y), c=:red,leg=true,label="Obstacle",leg=:bottomleft)
  scatter!((Xn,Yn), markershape = :square, markercolor = :black, markersize = 4,fill = (0, 1, :black),leg=true, grid=true,label="Vehicle")
end