text:
using Test
using Dorothy

@DorothyAll()

@testset "Secondary structure" begin

	@testset "Inference" begin
		m1 = readf("$(Dorothy.datapath)/1BTL.gro")
		inferss!(m1)
		@test m1.SS == ["C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "C", "C", "C", "C", "C", "C", "C", "C",
				"C", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "C", "C", "C", "C", "C", "C", "C", "C", "C", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "C",
				"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "C", "C", "C", "C", "C", "C", "C", "C", "G", "G", "G",
				"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
				"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
				"G", "G", "G", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "C", "C", "C", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "C", "C", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "G", "G", "G",
				"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
				"G", "G", "G", "G", "G", "G", "G", "G", "G", "C", "C", "C", "C",
				"C", "C", "C", "C", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
				"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C",
				"C", "C", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "C", "C", "C", "C", "C", "C", "C", "C", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "C", "C", "C",
				"C", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "B", "B", "B", "B", "B", "B",
				"B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "G", "G", "G",
				"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
				"G", "G", "G", "G", "G", "G", "G", "G", "G", "C", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "C", "C", "C", "C", "C", "C", "C", "C",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "C", "C", "C", "C", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
				"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
				"G", "G", "G", "G", "G", "G", "G", "C", "C", "C", "C", "C", "C",
				"C", "C", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T",
				"T", "T", "T", "T", "T", "T", "T", "T", "C", "C", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "C", "C", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
				"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C",
				"C", "C", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
				"H", "H", "H", "H", "C", "C", "C", "C", "C", "C", "C", "C", "C",
				"C", "C", "C", "C", "C", "C", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",
				"", "", "", ""]
	end

end # @testset


text:
function get_blocks(cliques)
    number_elements = length(collect(flatten(cliques)))
    blocks = Array{String,2}(undef, number_elements * 2, 2)
    blocks_list = Dict("B$i" => clique for (i, clique) in enumerate(cliques))
    i = 1
    for block in blocks_list
        for vertex in block.second
            blocks[i, 1] = block.first
            blocks[i, 2] = "VOLT_$(vertex)_Re"
            blocks[i + 1, 1] = block.first
            blocks[i + 1, 2] = "VOLT_$(vertex)_Im"
            i += 2
        end
    end
    return blocks
end

function get_cliquetree_mat(cliquetree)
    cliquetree_mat = Array{String,2}(undef, length(cliquetree), 2)
    for (index, neighbor) in enumerate(cliquetree)
        cliquetree_mat[index, 1] = "B$(neighbor[1])"
        cliquetree_mat[index, 2] =  "B$(neighbor[2])"
    end
    return cliquetree_mat
end

function get_opf_mat(path)
    return readdlm(path, skipstart=1)
end

function get_opf_ctr(path)
    return readdlm(path, skipstart=1)
end

function read_data(cliques, cliquetree, path_opf_ctr, path_opf_mat)
    return get_blocks(cliques), get_cliquetree_mat(cliquetree),
           get_opf_ctr(path_opf_ctr), get_opf_mat(path_opf_mat)
end

function read_mosek_log(log_file)
    search = ["nb_coupling_constraints primal :",
              "  Constraints            :",
              "Julia resolution time :",
              "Objective value: ",
              "memory : "]
    count_iter = false
    values = zeros(length(search) + 1)
    index = 1
    for line in readlines(log_file)
        if occursin(search[index], line)
            values[index] = parse(Float64, split(line, ':')[2])
            index < length(search) - 1 && (index += 1)
        elseif occursin("ITE", line)
            count_iter = true
            continue
        elseif occursin("Optimizer terminated.", line)
            count_iter = false
        end
        count_iter && (values[end] += 1)
    end
    return values
    #lines = readlines(log_file)
    #i = 1
    #while split(lines[i], ':')[1]!= "nb_coupling_constraints primal " && split(lines[i], ':')[1]!="nb_coupling_constraints "
    #    i+=1
    #end
    #nlc = parse(Int, split(lines[i], ':')[2])
    #while split(lines[i], ':')[1] != "  Constraints            "
    #    i+=1
    #end
    #n_total_constraints = parse(Int,split(lines[i], ':')[2])
    #m = n_total_constraints - nlc
    #while split(lines[i], ':')[1] != "Optimizer terminated. Time" && i < length(lines)
    #    i+=1
    #end
    #time = split(lines[i], ':')[2][2:end]
    #nb_iter = split(lines[i-1], ' ')[1]
    #objective = lines[i+2][18:end]
    #return time, nb_iter, objective, m, nlc
    return values
end

ext: jl

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using Enrich
using Test

max_stars_repo_name: Fredrikp-ume/Enrich

max_issues_repo_name: Fredrikp-ume/Enrich

max_forks_repo_name: Fredrikp-ume/Enrich


text:
struct UnionVector{
    T,
    ETS,
    TD <: AbstractVector,
    TM <: Vector{UInt8},
    TV <: Tuple
} <: Abstract.UnionVector{T}

    data::TD
    typeid::TM
    views::TV

    function UnionVector(
        ::Type{ETS},
        data::TD,
        typeid::TM,
        views::TV,
    ) where {
        ETS <: Tuple,
        TD <: AbstractVector,
        TM <: Vector{UInt8},
        TV <: Tuple,
    }
        return new{asunion(ETS), ETS, TD, TM, TV}(data, typeid, views)
    end
end

function UnionVector(ETS::Type, data::AbstractVector, typeid::Vector)
    views = foldltupletype((), ETS) do views, ET
        return (views..., reinterpret(ofsamesize(eltype(data), ET), data))
    end
    return UnionVector(ETS, data, typeid, views)
end

const ElTypeSpec = Union{Type{<:Tuple}, TypeTuple}
aseltypetuple(::Type{ETS}) where {ETS <: Tuple} = ETS
aseltypetuple(ETS::TypeTuple) = Tuple{ETS...}

eltypetuple(A::UnionVector{<:Any, ETS}) where ETS = ETS

Base.size(A::UnionVector) = size(A.data)

function UnionVector(::UndefInitializer, ETS::ElTypeSpec, n::Integer)
    elsize = max(sizeof.(astupleoftypes(ETS))...)
    typeid = zeros(UInt8, n)
    BT = Tuple{ntuple(_ -> UInt8, elsize)...}
    data = Vector{BT}(undef, n)
    return UnionVector(aseltypetuple(ETS), data, typeid)
end

function UnionVector(ETS::Tuple, items::AbstractVector)
    A = UnionVector(undef, ETS, length(items))
    for i in axes(A, 1)
        # TODO: find the best way to support arrays with offset
        @inbounds A[i] = items[i - firstindex(A) + firstindex(items)]
    end
    return A
end

function UnionVector(data::AbstractVector)
    ETS = foldl(data, init=()) do ETS, d
	T = typeof(d)
	if T ∉ ETS
	    ETS = (ETS..., T)
	end
	return ETS
    end
    return UnionVector(ETS, data)
end

Base.@propagate_inbounds typeat(A::UnionVector{<:Any, ETS}, i) where ETS =
    fieldtype(ETS, Int(A.typeid[i]))

Base.@propagate_inbounds Base.getindex(A::UnionVector, i::Int) =
    A.views[A.typeid[i]][i] |> unpad

# TODO: handle conversion
view_and_id(A::UnionVector, T::Type) =
    foldlargs(1, A.views...) do id, xs
        if paddedtype(eltype(xs)) === T
            reduced((xs, id))
        else
            id + 1
        end
    end |> ifunreduced() do _
        error("unsupported element type (type conversion not implemented yet)")
    end

# Base.checkbounds(::Type{Bool}, A::UnionVector, i) =
#     checkbounds(Bool, A.data, i) && checkbounds(Bool, A.typeid, i)

Base.@propagate_inbounds function Base.setindex!(A::UnionVector, v, i::Int)
    xs, id = view_and_id(A, typeof(v))
    A.typeid[i] = id
    xs[i] = v
    return
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    jdcnv(date) -> julian_days

### Purpose ###

Convert proleptic Gregorian Calendar date in UTC standard to number of Julian
days.

### Explanation ###

Takes the given proleptic Gregorian date in UTC standard and returns the number
of Julian calendar days since epoch `-4713-11-24T12:00:00`.

### Argument ###

* `date`: date in proleptic Gregorian Calendar.  Each element can be either a `DateTime` or
  anything that can be converted directly to `DateTime`.

### Output ###

Number of Julian days, as a floating point.

### Example ###

Find the Julian days number at 2016 August 23, 03:39:06.

```jldoctest
julia> using AstroLib, Dates

julia> jdcnv(DateTime(2016, 08, 23, 03, 39, 06))
2.4576236521527776e6

julia> jdcnv(2016, 08, 23, 03, 39, 06)
2.4576236521527776e6

julia> jdcnv("2016-08-23T03:39:06")
2.4576236521527776e6
```

### Notes ###

This is the inverse of `daycnv`.

`get_juldate` returns the number of Julian days for current time.  It is
equivalent to `jdcnv(now(Dates.UTC))`.

For the conversion of Julian date to number of Julian days, use `juldate`.
"""
const jdcnv = Dates.datetime2julian
jdcnv(date...) = jdcnv(DateTime(date...))


text:
using LinearAlgebra # for ⋅
using MathOptInterface
using JuMP
using SCIP  # for solver-specific callbacks
using SparseArrays

const MOI = MathOptInterface

include("common.jl")

# iterative algorithms (decomposition with two separate models)
include("itergbd.jl")
include("itertopo.jl")

# all-in-one integrated model
include("minlp.jl")

# constraint handlers
include("treetopohdlr.jl")
include("semisubhdlr.jl")
include("gbdsubhdlr.jl")

# callback-based algorithms: keep one master problem
include("cbtopo.jl")  # needs itertopo.jl
include("cbgbd.jl")   # needs itergbd.jl


text:
using Measurements: Measurement, measurement
using SymEngine: Basic, symbols
using Test: @test, @testset, @test_throws
using Unitful: Quantity, DimensionlessQuantity, @u_str

using EquationsOfStateOfSolids:
    Murnaghan,
    Murnaghan1st,
    BirchMurnaghan,
    BirchMurnaghan2nd,
    BirchMurnaghan3rd,
    BirchMurnaghan4th,
    PoirierTarantola,
    PoirierTarantola2nd,
    PoirierTarantola3rd,
    PoirierTarantola4th,
    Vinet,
    AntonSchmidt

@testset "Promoting `eltype`" begin
    @testset "Promoting to floating-point numbers" begin
        @test eltype(Murnaghan1st(1, 2, 3.0, 0)) === Float64
        @test eltype(BirchMurnaghan2nd(1, 2.0, 0)) === Float64
        @test eltype(BirchMurnaghan3rd(1, 2, 3.0, 0)) === Float64
        @test eltype(BirchMurnaghan4th(1, 2.0, 3, 4, 0)) === Float64
        @test eltype(PoirierTarantola2nd(1, 2.0, 0)) === Float64
        @test eltype(PoirierTarantola3rd(1, 2, 3.0, 0)) === Float64
        @test eltype(Vinet(1, 2, 3.0, 0)) === Float64
        @test eltype(AntonSchmidt(1, 2, 3.0, 0)) === Float64
        # @test eltype(BreenanStacey(1, 2, 3.0, 0)) === Float64
        @test eltype(Murnaghan1st{Float32}(Int(1), 2 // 1, Int8(3), Float64(4))) === Float32
        @test eltype(Murnaghan1st{BigFloat}(Int(1), 2 // 1, Int8(3), Float64(4))) ===
              BigFloat
        @test eltype(PoirierTarantola4th{Float16}(Int8(1), 2 // 1, 4, Int16(5), 6)) ===
              Float16
        @test eltype(BirchMurnaghan4th(Int8(1), 2 // 1, big(4.0), Int16(5), 6)) === BigFloat
        @test eltype(BirchMurnaghan4th(Int8(1), 2 // 1, big(4), Int16(5), 6.0)) === BigFloat
        @test_throws InexactError Vinet{Float64}(41 + 0.1im, 1.2, 4)
    end

    @testset "Promoting to rationals" begin
        @test eltype(PoirierTarantola4th(1, 2, 3, 4, 0)) === Int
        @test eltype(Murnaghan1st(Int32(1), Int16(2), Int8(3), 0)) === Int
        @test eltype(Murnaghan1st(1, 2 // 1, Int8(3), 0)) === Rational{Int}
        @test eltype(BirchMurnaghan2nd(1, Int8(2), 0)) === Int
        @test eltype(BirchMurnaghan2nd(1 // 1, Int32(2))) === Rational{Int}
        @test eltype(BirchMurnaghan3rd(Int8(1), 2, 4, 0)) === Int
        @test eltype(BirchMurnaghan3rd(Int8(1), 2 // 1, 4, 0)) === Rational{Int}
        @test eltype(BirchMurnaghan4th(Int8(1), 2, 4, Int16(5), 6)) === Int
        @test eltype(BirchMurnaghan4th(Int8(1), 2 // 1, 4, Int16(5), 6)) === Rational{Int}
        @test eltype(BirchMurnaghan4th(Int8(1), 2 // 1, big(4), Int16(5), 6)) ===
              Rational{BigInt}
        @test eltype(PoirierTarantola2nd(Int8(1), 2, 3)) === Int
        @test eltype(PoirierTarantola2nd(Int8(1), 2 // 1, 3)) === Rational{Int}
        @test eltype(PoirierTarantola3rd(Int8(1), 2, 3, Int16(4))) === Int
        @test eltype(PoirierTarantola3rd(Int8(1), 2 // 1, 3 // 1, Int16(4))) ===
              Rational{Int}
        @test eltype(PoirierTarantola4th(Int8(1), 2, 3, Int16(4), 5)) === Int
        @test eltype(PoirierTarantola4th(Int8(1), 2 // 1, 3, Int16(4), 5)) === Rational{Int}
        @test eltype(Vinet(Int8(1), 2, 3, Int16(4))) === Int
        @test eltype(Vinet(Int8(1), 2 // 1, 3, Int16(4))) === Rational{Int}
        @test eltype(AntonSchmidt(Int8(1), 2, 3, 0)) === Int
        @test eltype(AntonSchmidt(Int8(1), 2 // 1, 3, 0)) === Rational{Int}
        @test_throws InexactError PoirierTarantola3rd{BigInt}(41, 1 // 2, 4)
        @test_throws InexactError PoirierTarantola3rd{BigInt}(41, 1, 2.2)
    end

    @testset "Promoting with units" begin
        @test Murnaghan1st(1u"angstrom^3", 2u"eV/angstrom^3", 3.0, 4u"eV") ===
              Murnaghan1st(1.0u"angstrom^3", 2.0u"eV/angstrom^3", 3.0, 4.0u"eV")
        @test Murnaghan1st(1u"angstrom^3", 2u"eV/nm^3", 3 // 2, 4u"eV") ===
              Murnaghan1st((1 // 1)u"angstrom^3", (2 // 1)u"eV/nm^3", 3 // 2, (4 // 1)u"eV")
        @test BirchMurnaghan2nd(1u"angstrom^3", 2u"eV/angstrom^3", 3.0u"J") ===
              BirchMurnaghan2nd(1.0u"angstrom^3", 2.0u"eV/angstrom^3", 3.0u"J")
        BirchMurnaghan2nd((1 // 1)u"pm^3", (2 // 1)u"eV/angstrom^3", (3 // 1)u"eV")
        @test BirchMurnaghan3rd(1u"angstrom^3", 2u"GPa", 4.0, 3u"eV") ===
              BirchMurnaghan3rd(1.0u"angstrom^3", 2.0u"GPa", 4.0, 3.0u"eV")
        @test BirchMurnaghan3rd(1u"angstrom^3", 2u"GPa", 4 // 1, 3u"eV") ===
              BirchMurnaghan3rd(
            (1 // 1)u"angstrom^3",
            (2 // 1)u"GPa",
            4 // 1,
            (3 // 1)u"eV",
        )
        @test BirchMurnaghan4th(1u"nm^3", 2u"GPa", 3.0, 4u"GPa^-1", 5u"eV") ===
              BirchMurnaghan4th(1.0u"nm^3", 2.0u"GPa", 3.0, 4.0u"GPa^-1", 5.0u"eV")
        @test BirchMurnaghan4th(1u"nm^3", 2u"GPa", 3 // 1, 4u"1/GPa", 5u"J") ===
              BirchMurnaghan4th(
            (1 // 1)u"nm^3",
            (2 // 1)u"GPa",
            3 // 1,
            (4 // 1)u"1/GPa",
            (5 // 1)u"J",
        )
        @test PoirierTarantola2nd(1u"pm^3", 2u"GPa", 3.0u"eV") ===
              PoirierTarantola2nd(1.0u"pm^3", 2.0u"GPa", 3.0u"eV")
        @test PoirierTarantola2nd(1u"nm^3", 2u"GPa", (3 // 1)u"eV") ===
              PoirierTarantola2nd((1 // 1)u"nm^3", (2 // 1)u"GPa", (3 // 1)u"eV")
        @test PoirierTarantola3rd(1u"nm^3", 2u"GPa", 3, 4.0u"eV") ===
              PoirierTarantola3rd(1.0u"nm^3", 2.0u"GPa", 3, 4.0u"eV")
        @test PoirierTarantola3rd(1u"nm^3", 2u"GPa", 3, (4 // 1)u"eV") ===
              PoirierTarantola3rd((1 // 1)u"nm^3", (2 // 1)u"GPa", 3 // 1, (4 // 1)u"eV")
        @test PoirierTarantola4th(1u"nm^3", 2u"GPa", 3, 1 / 0.25u"Pa", 5.0u"eV") ===
              PoirierTarantola4th(1.0u"nm^3", 2.0u"GPa", 3.0, 4.0u"1/Pa", 5.0u"eV")
        @test PoirierTarantola4th(1u"nm^3", 2u"GPa", 3 // 1, 4u"GPa^(-1)", 5u"eV") ===
              PoirierTarantola4th(
            (1 // 1)u"nm^3",
            (2 // 1)u"GPa",
            3 // 1,
            (4 // 1)u"GPa^(-1)",
            (5 // 1)u"eV",
        )
        @test Vinet(1u"nm^3", 2u"GPa", 3, 4.0u"eV") ===
              Vinet(1.0u"nm^3", 2.0u"GPa", 3.0, 4.0u"eV")
        @test Vinet(1u"nm^3", 2u"GPa", 3, (4 // 1)u"eV") ===
              Vinet((1 // 1)u"nm^3", (2 // 1)u"GPa", 3 // 1, (4 // 1)u"eV")
        @test AntonSchmidt(1u"nm^3", 2u"GPa", 3.0, 4u"eV") ===
              AntonSchmidt(1.0u"nm^3", 2.0u"GPa", 3.0, 4.0u"eV")
        @test AntonSchmidt(1u"nm^3", 2u"GPa", 3 // 1, 4u"eV") ===
              AntonSchmidt((1 // 1)u"nm^3", (2 // 1)u"GPa", 3 // 1, (4 // 1)u"eV")
        # @test BreenanStacey(1u"nm^3", 2u"GPa", 3.0, 0u"eV") ===
        #       BreenanStacey{Quantity{Float64}}
        # @test BreenanStacey(1u"nm^3", 2u"GPa", 3 // 1, 0u"eV") ===
        #   BreenanStacey{Quantity{Rational{Int}}}
        @test BirchMurnaghan3rd(1u"angstrom^3", 2u"GPa", 4 // 1, 3u"eV").b′0 isa
              DimensionlessQuantity
        @test AntonSchmidt(1u"nm^3", 2u"GPa", 3.0, 4u"eV").b′0 isa DimensionlessQuantity
    end
end

@testset "Parameter `e0` and promotion" begin
    @test Murnaghan1st(1, 2, 3.0) === Murnaghan1st(1.0, 2.0, 3.0, 0.0)
    @test BirchMurnaghan2nd(1, 2.0) === BirchMurnaghan2nd(1.0, 2.0, 0.0)
    @test BirchMurnaghan3rd(1, 2, 3.0) === BirchMurnaghan3rd(1.0, 2.0, 3.0, 0.0)
    @test BirchMurnaghan4th(1, 2.0, 3, 4) === BirchMurnaghan4th(1.0, 2.0, 3.0, 4.0, 0.0)
    @test Vinet(1, 2, 3.0) === Vinet(1.0, 2.0, 3.0, 0.0)
    @test PoirierTarantola2nd(1, 2.0) === PoirierTarantola2nd(1.0, 2.0, 0.0)
    @test PoirierTarantola3rd(1, 2, 3.0) === PoirierTarantola3rd(1.0, 2.0, 3.0, 0.0)
    @test PoirierTarantola4th(1, 2, 3, 4) === PoirierTarantola4th(1, 2, 3, 4, 0)
    @test AntonSchmidt(1, 2, 3.0) === AntonSchmidt(1.0, 2.0, 3.0, 0.0)
    # @test BreenanStacey(1, 2, 3.0) === BreenanStacey(1.0, 2.0, 3.0, 0.0)
    @test eltype(Murnaghan1st(1u"angstrom^3", 2u"eV/angstrom^3", 3)) === Quantity{Int}
    @test eltype(Murnaghan1st(1u"angstrom^3", 2u"eV/angstrom^3", 3.0)) === Quantity{Float64}
    @test eltype(BirchMurnaghan2nd(1u"nm^3", 2u"GPa")) === Quantity{Int}
    @test eltype(BirchMurnaghan2nd(1u"nm^3", 2.0u"GPa")) === Quantity{Float64}
    @test eltype(BirchMurnaghan3rd(1u"nm^3", 2u"GPa", 4)) === Quantity{Int}
    @test eltype(BirchMurnaghan3rd(1u"nm^3", 2u"GPa", 4.0)) === Quantity{Float64}
    @test eltype(BirchMurnaghan4th(1u"nm^3", 2u"GPa", 3, 4u"1/GPa")) === Quantity{Int}
    @test eltype(BirchMurnaghan4th(1u"nm^3", 2u"GPa", 3, 4.0u"1/GPa")) === Quantity{Float64}
    @test eltype(PoirierTarantola2nd(1u"nm^3", 2u"GPa")) === Quantity{Int}
    @test eltype(PoirierTarantola2nd(1u"nm^3", 2.0u"GPa")) === Quantity{Float64}
    @test eltype(PoirierTarantola3rd(1u"nm^3", 2u"GPa", 3)) === Quantity{Int}
    @test eltype(PoirierTarantola3rd(1u"nm^3", 2u"GPa", 3.0)) === Quantity{Float64}
    @test eltype(PoirierTarantola4th(1u"nm^3", 2u"GPa", 3, 4u"1/GPa")) === Quantity{Int}
    @test eltype(PoirierTarantola4th(1u"nm^3", 2u"GPa", 3, 4.0u"1/GPa")) ===
          Quantity{Float64}
    @test eltype(Vinet(1u"nm^3", 2u"GPa", 3)) === Quantity{Int}
    @test eltype(Vinet(1u"nm^3", 2u"GPa", 3.0)) === Quantity{Float64}
    @test eltype(AntonSchmidt(1u"nm^3", 2u"GPa", 3)) === Quantity{Int}
    @test eltype(AntonSchmidt(1u"nm^3", 2u"GPa", 3.0)) === Quantity{Float64}
    #@test eltype(BreenanStacey(1u"nm^3", 2u"GPa", 3))
    #@test eltype(BreenanStacey(1u"nm^3", 2u"GPa", 3.0)
end

@testset "Constructors" begin
    @test_throws ArgumentError Murnaghan(1, 2, 3.0)
    @test_throws ArgumentError Murnaghan(1, 2, 3.0, 4, 5 // 1, Int32(6))
    @test Murnaghan(1, 2, 3.0, 4) === Murnaghan1st(1.0, 2.0, 3.0, 4.0)
    @test Murnaghan(1u"angstrom^3", 2u"eV/angstrom^3", 3.0, 4u"eV") ===
          Murnaghan1st(1.0u"angstrom^3", 2.0u"eV/angstrom^3", 3.0, 4.0u"eV")
    @test_throws ArgumentError BirchMurnaghan(1, 2)
    @test BirchMurnaghan(1, 2 // 1, 3.0) === BirchMurnaghan2nd(1.0, 2.0, 3.0)
    @test BirchMurnaghan(1, 2, 3.0, 4) === BirchMurnaghan3rd(1.0, 2.0, 3.0, 4.0)
    @test BirchMurnaghan(1, 2, 3.0, 4, Int32(5)) ===
          BirchMurnaghan4th(1.0, 2.0, 3.0, 4.0, 5.0)
    @test BirchMurnaghan(1u"angstrom^3", 2u"eV/angstrom^3", 3.0, 4u"eV") ===
          BirchMurnaghan3rd(1.0u"angstrom^3", 2.0u"eV/angstrom^3", 3.0, 4.0u"eV")
    @test_throws ArgumentError PoirierTarantola(1, 2)
    @test PoirierTarantola(1, 2 // 1, 3.0) === PoirierTarantola2nd(1.0, 2.0, 3.0)
    @test PoirierTarantola(1, 2, 3.0, 4) === PoirierTarantola3rd(1.0, 2.0, 3.0, 4.0)
    @test PoirierTarantola(1, 2, 3.0, 4, Int32(5)) ===
          PoirierTarantola4th(1.0, 2.0, 3.0, 4.0, 5.0)
    @test PoirierTarantola(1u"angstrom^3", 2u"eV/angstrom^3", 3.0, 4u"eV") ===
          PoirierTarantola3rd(1.0u"angstrom^3", 2.0u"eV/angstrom^3", 3.0, 4.0u"eV")
end

@testset "`float` on an EOS" begin
    @test eltype(float(Vinet(1, 2, 3))) === Float64
    @test eltype(float(Vinet(1u"nm^3", 2u"GPa", 3))) === Quantity{Float64}
    @test eltype(float(Murnaghan1st(big(2), 3, 4, 5))) === BigFloat
end

@testset "Other element types" begin
    @test eltype(
        BirchMurnaghan4th(measurement("1 +- 0.1"), 3 // 1, 2, measurement("-123.4(56)")),
    ) === Measurement{Float64}
    @testset "`SymEngine.Basic`" begin
        v0, b0, b′0, b″0, e0 = symbols("v0, b0, b′0, b″0, e0")
        @test eltype(BirchMurnaghan4th(v0, b0, b′0, b″0, e0)) === Basic
        @test eltype(BirchMurnaghan4th(v0, b0, b′0, b″0)) === Basic
        @test iszero(BirchMurnaghan4th(v0, b0, b′0, b″0).e0)
    end
end

@testset "Test equality" begin
    @test BirchMurnaghan3rd(1, 2, 3, 4) == BirchMurnaghan3rd(1.0, 2.0, 3.0, 4.0)
    @test BirchMurnaghan3rd(1 // 1, 2 // 1, 6 // 2, 12 // 3) ==
          BirchMurnaghan3rd(1.0, 2.0, 3.0, 4.0)
    @test BirchMurnaghan3rd(1u"angstrom^3", 2u"GPa", 3) ==
          BirchMurnaghan3rd(1.0u"angstrom^3", 2.0u"GPa", 3.0)
end


text:
###############################################
## Reading and writing Julia data .jld files ##
###############################################

module JLD00
using Printf
using HDF5
# Add methods to...
import HDF5: close, dump, exists, file, getindex, setindex!, g_create, g_open, o_delete, name, names, read, size, write,
             HDF5ReferenceObj, HDF5BitsKind, ismmappable, readmmap
import Base: length, show, delete!, iterate
import ..JLD

# See julia issue #8907
replacements = Any[]
push!(replacements, :(s = replace(s, r"Uint(?=\d{1,3})" => "UInt")))
push!(replacements, :(s = replace(s, r"ASCIIString|UTF8String|ByteString" => "String")))
ex = Expr(:block, replacements...)
@eval function julia_type(s::AbstractString)
    $ex
    _julia_type(s)
end

const magic_base = "Julia data file (HDF5), version "
const version_current = "0.0.2"
const pathrefs = "/_refs"
const pathtypes = "/_types"
const pathrequire = "/_require"
const name_type_attr = "julia type"

### Dummy types used for converting attribute strings to Julia types
mutable struct UnsupportedType; end
mutable struct UnconvertedType; end
mutable struct CompositeKind; end   # here this means "a type with fields"

# The Julia Data file type
# Purpose of the nrefs field:
# length(group) only returns the number of _completed_ items in a group. Since
# we'll write recursively, we need to keep track of the number of reference
# objects _started_.
mutable struct JldFile <: HDF5.DataFile
    plain::HDF5File
    version::String
    toclose::Bool
    writeheader::Bool
    mmaparrays::Bool

    function JldFile(plain::HDF5File, version::AbstractString=version_current, toclose::Bool=true,
                     writeheader::Bool=false, mmaparrays::Bool=false)
        f = new(plain, version, toclose, writeheader, mmaparrays)
        if toclose
            finalizer(close, f)
        end
        f
    end
end

struct JldGroup
    plain::HDF5Group
    file::JldFile
end

struct JldDataset
    plain::HDF5Dataset
    file::JldFile
end

file(x::JldFile) = x
file(x::Union{JldGroup, JldDataset}) = x.file

function close(f::JldFile)
    if f.toclose
        close(f.plain)
        if f.writeheader
            magic = zeros(UInt8, 512)
            tmp = string(magic_base, f.version)
            magic[1:length(tmp)] = Vector{UInt8}(codeunits(tmp))
            rawfid = open(f.plain.filename, "r+")
            write(rawfid, magic)
            close(rawfid)
        end
        f.toclose = false
    end
    nothing
end
close(g::Union{JldGroup, JldDataset}) = close(g.plain)
show(io::IO, fid::JldFile) = isvalid(fid.plain) ? print(io, "Julia data file version ", fid.version, ": ", fid.plain.filename) : print(io, "Julia data file (closed): ", fid.plain.filename)

function jldopen(filename::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Bool, ff::Bool; mmaparrays::Bool=false)
    local fj
    if ff && !wr
        error("Cannot append to a write-only file")
    end
    if !cr && !isfile(filename)
        error("File ", filename, " cannot be found")
    end
    version = version_current
    pa = p_create(HDF5.H5P_FILE_ACCESS)
    try
        pa["fclose_degree"] = HDF5.H5F_CLOSE_STRONG
        if cr && (tr || !isfile(filename))
            # We're truncating, so we don't have to check the format of an existing file
            # Set the user block to 512 bytes, to save room for the header
            p = p_create(HDF5.H5P_FILE_CREATE)
            local f
            try
                p["userblock"] = 512
                f = HDF5.h5f_create(filename, HDF5.H5F_ACC_TRUNC, p.id, pa.id)
            finally
                close(p)
            end
            fj = JldFile(HDF5File(f, filename), version, true, true, mmaparrays)
            # initialize empty require list
            write(fj, pathrequire, String[])
        else
            # Test whether this is a jld file
            sz = filesize(filename)
            if sz < 512
                error("File size indicates $filename cannot be a Julia data file")
            end
            magic = Vector{UInt8}(undef, 512)
            rawfid = open(filename, "r")
            try
                magic = read!(rawfid, magic)
            finally
                close(rawfid)
            end
            if length(magic) ≥ ncodeunits(magic_base) && view(magic, 1:ncodeunits(magic_base)) == Vector{UInt8}(codeunits(magic_base))
                f = HDF5.h5f_open(filename, wr ? HDF5.H5F_ACC_RDWR : HDF5.H5F_ACC_RDONLY, pa.id)
                version = unsafe_string(pointer(magic) + length(magic_base))
                fj = JldFile(HDF5File(f, filename), version, true, true, mmaparrays)
                # Load any required files/packages
                if exists(fj, pathrequire)
                    r = read(fj, pathrequire)
                    for fn in r
                        require(fn)
                    end
                end
            else
                if ishdf5(filename)
                    println("$filename is an HDF5 file, but it is not a recognized Julia data file. Opening anyway.")
                    fj = JldFile(h5open(filename, rd, wr, cr, tr, ff), version_current, true, false, mmaparrays)
                else
                    error("$filename does not seem to be a Julia data or HDF5 file")
                end
            end
        end
    finally
        close(pa)
    end
    return fj
end

function jldopen(fname::AbstractString, mode::AbstractString="r"; mmaparrays::Bool=false)
    mode == "r"  ? jldopen(fname, true , false, false, false, false, mmaparrays=mmaparrays) :
    mode == "r+" ? jldopen(fname, true , true , false, false, false, mmaparrays=mmaparrays) :
    mode == "w"  ? jldopen(fname, false, true , true , true , false, mmaparrays=mmaparrays) :
#     mode == "w+" ? jldopen(fname, true , true , true , true , false) :
#     mode == "a"  ? jldopen(fname, false, true , true , false, true ) :
#     mode == "a+" ? jldopen(fname, true , true , true , false, true ) :
    error("invalid open mode: ", mode)
end

function jldopen(f::Function, args...)
    jld = jldopen(args...)
    try
        f(jld)
    finally
        close(jld)
    end
end

function jldobject(obj_id::HDF5.Hid, parent)
    obj_type = HDF5.h5i_get_type(obj_id)
    obj_type == HDF5.H5I_GROUP ? JldGroup(HDF5Group(obj_id, file(parent.plain)), file(parent)) :
    obj_type == HDF5.H5I_DATATYPE ? HDF5Datatype(obj_id) :
    obj_type == HDF5.H5I_DATASET ? JldDataset(HDF5Dataset(obj_id, file(parent.plain)), file(parent)) :
    error("Invalid object type for path ", path)
end

getindex(parent::Union{JldFile, JldGroup}, path::String) =
    jldobject(HDF5.h5o_open(parent.plain.id, path), parent)

function getindex(parent::Union{JldFile, JldGroup, JldDataset}, r::HDF5ReferenceObj)
    if r == HDF5.HDF5ReferenceObj_NULL; error("Reference is null"); end
    obj_id = HDF5.h5r_dereference(parent.plain.id, HDF5.H5R_OBJECT, r)
    jldobject(obj_id, parent)
end

### "Inherited" behaviors
g_create(parent::Union{JldFile, JldGroup}, args...) = JldGroup(g_create(parent.plain, args...), file(parent))
function g_create(f::Function, parent::Union{JldFile, JldGroup}, args...)
    g = JldGroup(g_create(parent.plain, args...), file(parent))
    try
        f(g)
    finally
        close(g)
    end
end
g_open(parent::Union{JldFile, JldGroup}, args...) = JldGroup(g_open(parent.plain, args...), file(parent))
name(p::Union{JldFile, JldGroup, JldDataset}) = name(p.plain)
exists(p::Union{JldFile, JldGroup, JldDataset}, path::String) = exists(p.plain, path)
root(p::Union{JldFile, JldGroup, JldDataset}) = g_open(file(p), "/")
o_delete(parent::Union{JldFile, JldGroup}, args...) = o_delete(parent.plain, args...)
function ensurepathsafe(path::String)
    if any([startswith(path, s) for s in (pathrefs,pathtypes,pathrequire)])
        error("$name is internal to the JLD format, use o_delete if you really want to delete it")
    end
end
function delete!(o::JldDataset)
    fullpath = name(o)
    ensurepathsafe(fullpath)
    o_delete(o.file, fullpath)
    refspath = joinpath(pathrefs, fullpath[2:end])
    exists(o.file, refspath) && o_delete(o.file, refspath)
end
function delete!(g::JldGroup)
    fullpath = name(g)
    ensurepathsafe(fullpath)
    for o in g typeof(o) == JldDataset && delete!(o) end
    o_delete(g.file,name(g))
end
function delete!(parent::Union{JldFile, JldGroup}, path::String)
    exists(parent, path) || error("$path does not exist in $parent")
    delete!(parent[path])
end
delete!(parent::Union{JldFile, JldGroup}, args::Tuple{Vararg{String}}) = for a in args delete!(parent,a) end
ismmappable(obj::JldDataset) = ismmappable(obj.plain)
readmmap(obj::JldDataset, args...) = readmmap(obj.plain, args...)
setindex!(parent::Union{JldFile, JldGroup}, val, path::String) = write(parent, path, val)

Base.iterate(parent::Union{JldFile, JldGroup}, state=(names(parent), 1)) = state[2] > length(state[1]) ? nothing :
                                                     (parent[state[1][state[2]]], (state[1], state[2]+1))

### Julia data file format implementation ###

### Read ###

function read(parent::Union{JldFile, JldGroup}, name::String)
    local val
    obj = parent[name]
    try
        val = read(obj)
    finally
        close(obj)
    end
    val
end
read(parent::Union{JldFile,JldGroup}, name::Symbol) = read(parent, string(name))

# read and readsafely differ only in how they handle CompositeKind
function read(obj::Union{JldFile, JldDataset})
    if !exists(attrs(obj.plain), name_type_attr)
        # Fallback to plain read
        return read(obj.plain)
    end
    # Read the type
    typename = a_read(obj.plain, name_type_attr)
    if typename == "Tuple"
        return read_tuple(obj)
    end
    # Convert to Julia type
    T = julia_type(typename)
    if T == CompositeKind
        # Use type information in the file to ensure we find the right module
        typename = a_read(obj.plain, "CompositeKind")
        try
            gtypes = root(obj)[pathtypes]
            try
                objtype = gtypes[typename]
                try
                    modnames = a_read(objtype.plain, "Module")
                    mod = Main
                    for mname in modnames
                        mod = Core.eval(mod, Symbol(mname))
                    end
                    T = Core.eval(mod, Symbol(typename))
                finally
                    close(objtype)
                end
            finally
                close(gtypes)
            end
        catch
            error("Type ", typename, " is not recognized. As a fallback, you can load ", name(obj), " with readsafely().")
        end
    end
    if T == UnsupportedType
        error("Type expression ", typename, ", is not recognized.")
    end
    read(obj, T)
end
function readsafely(obj::Union{JldFile, JldDataset})
    if !exists(attrs(obj.plain), name_type_attr)
        # Fallback to plain read
        return read(obj.plain)
    end
    # Read the type
    typename = a_read(obj.plain, name_type_attr)
    if typename == "Tuple"
        return read_tuple(obj)
    end
    # Convert to Julia type
    T = julia_type(typename)
    local ret
    if T == CompositeKind
        # Read as a dict
        typename = a_read(obj.plain, "CompositeKind")
        gtypes = root(obj)[pathtypes]
        try
            objtype = gtypes[typename]
            try
                n = read(objtype)
                v = getrefs(obj, Any)
                ret = Dict(zip(n[1,:], v))
            finally
                close(objtype)
            end
        finally
            close(gtypes)
        end
    else
        ret = read(obj, T)
    end
    ret
end
function readsafely(parent::Union{JldFile, JldGroup}, name::String)
    local ret
    obj = parent[name]
    try
        ret = readsafely(JldDataset(obj.plain, obj.file))
    finally
        close(obj)
    end
    return ret
end
readsafely(parent::Union{JldFile,JldGroup}, name::Symbol) = readsafely(parent, string(symbol))

# Basic types
const BitsKindOrString = Union{HDF5BitsKind, String}
read(obj::JldDataset, ::Type{T}) where {T<:BitsKindOrString} = read(obj.plain, T)
function read(obj::JldDataset, ::Type{Array{T}}) where T<:HDF5BitsKind
    A = obj.file.mmaparrays && HDF5.iscontiguous(obj.plain) ? readmmap(obj.plain, Array{T}) : read(obj.plain, Array{T})
    if isempty(A) && exists(obj, "dims")
        dims = a_read(obj.plain, "dims")
        A = reshape(A, dims...)
    end
    A
end
read(obj::JldDataset, ::Type{Array{T}}) where {T<:String} = read(obj.plain, Array{T})
read(obj::JldDataset, ::Type{Array{T,N}}) where {T<:BitsKindOrString,N} = read(obj, Array{T})

# Arrays-of-arrays of basic types
function read(obj::JldDataset, ::Type{Array{Array{T,N},M}}) where {T<:HDF5BitsKind,M,N}
    # fallback for backwards compatibility with pre-v0.2.27 format
    HDF5.hdf5_to_julia_eltype(datatype(obj.plain)) == HDF5ReferenceObj &&
        return getrefs(obj, Array{T,N})
    A = read(obj.plain, HDF5.HDF5Vlen{T})
    if isempty(A) && exists(obj, "dims")
        dims = a_read(obj.plain, "dims")
        A = reshape(A, dims...)
    end
    convert(Array{Array{T,N},M}, A)
end

# Nothing
read(obj::JldDataset, ::Type{Nothing}) = nothing
read(obj::JldDataset, ::Type{Bool}) = read(obj, UInt8) != 0

# Types
read(obj::JldDataset, ::Type{Type{T}}) where {T} = T

# Bool
function read(obj::JldDataset, ::Type{Array{Bool,N}}) where N
    format = a_read(obj.plain, "julia_format")
    if format == "EachUint8"
        bool(read(obj.plain, Array{UInt8}))
    else
        error("bool format not recognized")
    end
end

# Complex
function read(obj::JldDataset, ::Type{Complex{T}}) where T
    a = read(obj.plain, Array{T})
    a[1]+a[2]*im
end
function read(obj::JldDataset, ::Type{Array{T,N}}) where {T<:Complex,N}
    A = read(obj, Array{realtype(T)})
    reshape(reinterpret(T, vec(A)), ntuple(i->size(A, i+1), ndims(A)-1))
end

# Symbol
read(obj::JldDataset, ::Type{Symbol}) = Symbol(read(obj.plain, String))
read(obj::JldDataset, ::Type{Array{Symbol,N}}) where {N} = map(Symbol, read(obj.plain, Array{String}))

# Char
read(obj::JldDataset, ::Type{Char}) = Char(read(obj.plain, UInt32))

# General arrays
read(obj::JldDataset, t::Type{Array{T,N}}) where {T,N} = getrefs(obj, T)

# Tuple
function read_tuple(obj::JldDataset)
    t = read(obj, Array{Any, 1})
    return tuple(t...)
end
function read_tuple(obj::JldDataset, indices::AbstractVector)
    t = read(obj, Array{Any, 1})
    return tuple(t...)
end

# Dict
function read(obj::JldDataset, ::Type{T}) where T<:AbstractDict
    kv = getrefs(obj, Any)
    ret = T()
    for (cn, c) in zip(kv[1], kv[2])
        ret[cn] = c
    end
    ret
end

# Expressions
function read(obj::JldDataset, ::Type{Expr})
    a = getrefs(obj, Any)
    Expr(a[1], a[2]...)
end

# CompositeKind
read(obj::JldDataset, T::UnionAll) = read(obj, Base.unwrap_unionall(T))
function read(obj::JldDataset, T::DataType)
    if isempty(fieldnames(T)) && T.size > 0
        return read_bitstype(obj, T)
    end
    local x
    # Add the parameters
    if exists(obj, "TypeParameters")
        params = a_read(obj.plain, "TypeParameters")
        if !isempty(params)
            p = Vector{Any}(undef, length(params))
            for i = 1:length(params)
                p[i] = Core.eval(@__MODULE__, Meta.parse(params[i]))
            end
            T = T.name.wrapper
            T = T{p...}
        end
    end
    v = getrefs(obj, Any)
    if length(v) == 0
        x = ccall(:jl_new_struct, Any, (Any,Any...), T)
    else
        n = fieldnames(T)
        if length(v) != length(n)
            error("Wrong number of fields")
        end
        if !T.mutable
            x = ccall(:jl_new_structv, Any, (Any,Ptr{Cvoid},UInt32), T, v, length(fieldnames(T)))
        else
            x = ccall(:jl_new_struct_uninit, Any, (Any,), T)
            for i = 1:length(v)
                if isassigned(v, i)
                    setfield!(x, n[i], v[i])
                end
            end
        end
    end
    x
end

function read_bitstype(obj::JldDataset, T::DataType)
    a = read(obj.plain)
    reinterpret(T, a[1])
end

# Read an array of references
function getrefs(obj::JldDataset, ::Type{T}) where T
    refs = read(obj.plain, Array{HDF5ReferenceObj})
    out = Array{T}(undef, size(refs))
    f = file(obj)
    for i = 1:length(refs)
        if refs[i] != HDF5.HDF5ReferenceObj_NULL
            ref = f[refs[i]]
            try
                out[i] = read(ref)
            finally
                close(ref)
            end
        end
    end
    return out
end
function getrefs(obj::JldDataset, ::Type{T}, indices::Union{Integer, AbstractVector}...) where T
    refs = read(obj.plain, Array{HDF5ReferenceObj})
    refs = refs[indices...]
    f = file(obj)
    local out
    if isa(refs, HDF5ReferenceObj)
        # This is a scalar, not an array
        ref = f[refs]
        try
            out = read(ref)
        finally
            close(ref)
        end
    else
        out = Array{T}(undef, size(refs))
        for i = 1:length(refs)
            ref = f[refs[i]]
            try
                out[i] = read(ref)
            finally
                close(ref)
            end
        end
    end
    return out
end

### Writing ###

# Write "basic" types
function write(parent::Union{JldFile, JldGroup}, name::String,
                   data::Union{T, StridedArray{T}}, astype::String) where T<:Union{HDF5BitsKind, String}
    # Create the dataset
    dset, dtype = d_create(parent.plain, name, data)
    try
        # Write the attribute
        a_write(dset, name_type_attr, astype)
        isa(data, StridedArray) && isempty(data) && a_write(dset, "dims", [size(data)...])
        # Write the data
        HDF5.writearray(dset, dtype.id, data)
    finally
        close(dset)
        close(dtype)
    end
end
write(parent::Union{JldFile, JldGroup}, name::String, data::Union{T, Array{T}}) where {T<:Union{HDF5BitsKind, String}} =
    write(parent, name, data, full_typename(typeof(data)))

# Arrays-of-arrays of basic types
write(parent::Union{JldFile, JldGroup}, name::String,
            data::Array{Array{T,1}}, astype::String) where {T<:Union{HDF5BitsKind, String}} =
    write(parent, name, HDF5.HDF5Vlen(data), astype)
write(parent::Union{JldFile, JldGroup}, name::String,
            data::Array{Array{T,1}}) where {T<:Union{HDF5BitsKind, String}} =
    write(parent, name, data, full_typename(typeof(data)))
function write(parent::Union{JldFile, JldGroup}, name::String,
                  data::HDF5.HDF5Vlen{T}, astype::String) where T
    # Create the dataset
    dset, dtype = d_create(parent.plain, name, data)
    try
        # Write the attribute
        a_write(dset, name_type_attr, astype)
        isa(data, Array) && isempty(data) && a_write(dset, "dims", [size(data)...])
        # Write the data
        HDF5.writearray(dset, dtype.id, data)
    finally
        close(dset)
        close(dtype)
    end
end

# Write nothing
function write(parent::Union{JldFile, JldGroup}, name::String, n::Nothing, astype::String)
    local dspace, dset
    try
        dspace = dataspace(nothing)
        dset = HDF5Dataset(HDF5.h5d_create(HDF5.parents_create(HDF5.checkvalid(parent.plain), name, HDF5.H5T_NATIVE_UINT8, dspace.id,
                           HDF5.H5P_DEFAULT, HDF5.H5P_DEFAULT, HDF5.H5P_DEFAULT)...), file(parent.plain))
        a_write(dset, name_type_attr, astype)
    finally
        close(dspace)
        close(dset)
    end
end
write(parent::Union{JldFile, JldGroup}, name::String, n::Nothing) = write(parent, name, n, "Nothing")

# Types
# the first is needed to avoid an ambiguity warning
write(parent::Union{JldFile, JldGroup}, name::String, t::Tuple{Vararg{Type{T}}}) where {T} =
    write(parent, name, Any[t...], "Tuple")
write(parent::Union{JldFile, JldGroup}, name::String, t::Type{T}) where {T} =
    write(parent, name, nothing, string("Type{", full_typename(t), "}"))

# Bools
write(parent::Union{JldFile, JldGroup}, name::String, tf::Bool) = write(parent, name, uint8(tf), "Bool")
function write(parent::Union{JldFile, JldGroup}, name::String, tf::Array{Bool})
    write(parent, name, uint8(tf), full_typename(typeof(tf)))
    a_write(parent[name].plain, "julia_format", "EachUint8")
end

# Complex
realtype(::Type{Complex{T}}) where {T} = T
function write(parent::Union{JldFile, JldGroup}, name::String, c::Complex)
    reim = [real(c), imag(c)]
    write(parent, name, reim, full_typename(typeof(c)))
end
function write(parent::Union{JldFile, JldGroup}, name::String, C::Array{T}) where T<:Complex
    reim = reinterpret(realtype(T), reshape(C, ntuple(i->i==1 ? 1 : size(C,i-1), ndims(C)+1)))
    write(parent, name, reim, full_typename(typeof(C)))
end

# Int128/UInt128

# Symbols
write(parent::Union{JldFile, JldGroup}, name::String, sym::Symbol) = write(parent, name, string(sym), "Symbol")
write(parent::Union{JldFile, JldGroup}, name::String, syms::Array{Symbol}) = write(parent, name, map(string, syms), full_typename(typeof(syms)))

# Char
write(parent::Union{JldFile, JldGroup}, name::String, char::Char) = write(parent, name, UInt32(char), "Char")

# General array types (as arrays of references)
function write(parent::Union{JldFile, JldGroup}, path::String, data::Array{T}, astype::AbstractString) where T
    local gref  # a group, inside /_refs, for all the elements in data
    local refs
    # Determine whether parent already exists in /_refs, so we can avoid group/dataset conflict
    pname = name(parent)
    if startswith(pname, pathrefs)
        gref = g_create(parent, path*"g")
    else
        pathr = HDF5.joinpathh5(pathrefs, pname, path)
        if exists(file(parent), pathr)
            gref = g_open(file(parent), pathr)
        else
            gref = g_create(file(parent), pathr)
        end
    end
    grefname = name(gref)
    try
        # Write the items to the reference group
        refs = Array{HDF5ReferenceObj}(undef, size(data))
        # pad with zeros to keep in order
        nd = ndigits(length(data))
        z = "0"
        z = z[ones(Int, nd-1)]
        nd = 1
        for i = 1:length(data)
            if isassigned(data, i)
                if ndigits(i) > nd
                    nd = ndigits(i)
                    z = z[1:end-1]
                end
                itemname = z*string(i)
                write(gref, itemname, data[i])
                # Extract references
                tmp = gref[itemname]
                refs[i] = HDF5ReferenceObj(tmp.plain, grefname*"/"*itemname)
                close(tmp)
            else
                refs[i] = HDF5.HDF5ReferenceObj_NULL
            end
        end
    finally
        close(gref)
    end
    # Write the references as the chosen variable
    cset, ctype = d_create(parent.plain, path, refs)
    try
        HDF5.writearray(cset, ctype.id, refs)
        a_write(cset, name_type_attr, astype)
    finally
        close(ctype)
        close(cset)
    end
end
write(parent::Union{JldFile, JldGroup}, path::String, data::Array{T}) where {T} =
    write(parent, path, data, full_typename(typeof(data)))

# Tuple
write(parent::Union{JldFile, JldGroup}, name::String, t::Tuple) = write(parent, name, Any[t...], "Tuple")

# AbstractDict
function write(parent::Union{JldFile, JldGroup}, name::String, d::AbstractDict)
    tn = full_typename(typeof(d))
    if tn == "DataFrame"
        return write_composite(parent, name, d)
    end
    n = length(d)
    ks = Vector{keytype(d)}(undef, n)
    vs = Vector{valtype(d)}(undef, n)
    i = 0
    for (k,v) in d
        ks[i+=1] = k
        vs[i] = v
    end
    da = Any[ks, vs]
    write(parent, name, da, tn)
end

# Expressions
function write(parent::Union{JldFile, JldGroup}, name::String, ex::Expr)
    args = ex.args
    # Discard "line" expressions
    keep = trues(length(args))
    for i = 1:length(args)
        if isa(args[i], Expr) && args[i].head == :line
            keep[i] = false
        end
    end
    args = args[keep]
    a = Any[ex.head, args]
    write(parent, name, a, "Expr")
end

# CompositeKind
write(parent::Union{JldFile, JldGroup}, name::String, s; rootmodule="") = write_composite(parent, name, s; rootmodule=rootmodule)

function write_composite(parent::Union{JldFile, JldGroup}, name::String, s; rootmodule="")
    T = typeof(s)
    if isempty(fieldnames(T))
        if T.size > 0
            return write_bitstype(parent, name, s)
        end
        isdefined(T, :instance) || error("This is the write function for CompositeKind, but the input is of type ", T)
    end
    if has_pointer_field(s, name)
        return
    end
    Tname = string(T.name.name)
    n = fieldnames(T)
    local gtypes
    if !exists(file(parent), pathtypes)
        gtypes = g_create(file(parent), pathtypes)
    else
        gtypes = parent[pathtypes]
    end
    try
        if !exists(gtypes, Tname)
            # Write names to a dataset, so that other languages reading this file can
            # at least create a sensible dict
            nametype = Matrix{String}(2, length(n))
            t = T.types
            for i = 1:length(n)
                nametype[1, i] = string(n[i])
                nametype[2, i] = string(t[i])
            end
            write(gtypes.plain, Tname, nametype)
            obj = gtypes[Tname]
            # Write the module name as an attribute
            mod = Base.fullname(T.name.module)
            modnames = [map(string, mod)...]
            indx = findfirst(x->x==rootmodule, modnames)
            if indx > 0
                modnames = modnames[indx+1:end]
            end
            a_write(obj.plain, "Module", modnames)
            close(obj)
        end
    finally
        close(gtypes)
    end
    # Write the data
    v = Vector{Any}(undef, length(n))
    for i = 1:length(v)
        if isdefined(s, n[i])
            v[i] = getfield(s, n[i])
        end
    end
    write(parent, name, v, "CompositeKind")
    obj = parent[name]
    a_write(obj.plain, "CompositeKind", Tname)
    params = [map(full_typename, T.parameters)...]
    a_write(obj.plain, "TypeParameters", params)
    close(obj)
end

function write_bitstype(parent::Union{JldFile, JldGroup}, name::String, s)
    T = typeof(s)
    if T.size == 1
        ub = reinterpret(UInt8, s)
    elseif T.size == 2
        ub = reinterpret(UInt16, s)
    elseif T.size == 4
        ub = reinterpret(UInt32, s)
    elseif T.size == 8
        ub = reinterpret(UInt64, s)
    else
        error("Unsupported bitstype $T of size $(T.size)")
    end
    write(parent, name, [ub], "$(full_typename(T))")
end

function has_pointer_field(obj::Tuple, name)
    for o in obj
        if has_pointer_field(o, name)
            return true
        end
    end
    false
end

function has_pointer_field(obj, name)
    names = fieldnames(typeof(obj))
    for fieldname in names
        if isdefined(obj, fieldname)
            x = getfield(obj, fieldname)
            if isa(x, Ptr)
                @warn("Skipping $name because field \"$fieldname\" is a pointer")
                return true
            end
            if !isa(x, AbstractDict) && has_pointer_field(x, name)
                return true
            end
        end
    end
    false
end

### Size, length, etc ###
function size(dset::JldDataset)
    if !exists(attrs(dset.plain), name_type_attr)
        return size(dset.plain)
    end
    # Read the type
    typename = a_read(dset.plain, name_type_attr)
    if typename == "Tuple"
        return size(dset.plain)
    end
    # Convert to Julia type
    T = julia_type(typename)
    if T == CompositeKind || T <: AbstractDict || T == Expr
        return ()
    elseif T <: Complex
        return ()
    elseif isarraycomplex(T)
        sz = size(dset.plain)
        return sz[2:end]
    end
    size(dset.plain)
end
length(dset::JldDataset) = prod(size(dset))
lastindex(dset::JldDataset) = length(dset)

isarraycomplex(::Type{Array{T, N}}) where {T<:Complex, N} = true
isarraycomplex(t) = false

### Read/write via getindex/setindex! ###
function getindex(dset::JldDataset, indices::Union{Integer, Base.RangeIndex}...)
    if !exists(attrs(dset.plain), name_type_attr)
        # Fallback to plain read
        return getindex(dset.plain, indices...)
    end
    # Read the type
    typename = a_read(dset.plain, name_type_attr)
    if typename == "Tuple"
        return read_tuple(dset, indices...)
    end
    # Convert to Julia type
    T = julia_type(typename)
    _getindex(dset, T, indices...)
end

_getindex(dset::JldDataset, ::Type{Array{T,N}}, indices::Base.RangeIndex...) where {T<:HDF5BitsKind,N} =
    HDF5._getindex(dset.plain, T, indices...)
function _getindex(dset::JldDataset, ::Type{Array{T,N}}, indices::Base.RangeIndex...) where {T<:Complex,N}
    reinterpret(T, HDF5._getindex(dset.plain, realtype(T), 1:2, indices...), ntuple(i->length(indices[i]), length(indices)))
end
function _getindex(dset::JldDataset, ::Type{Array{Bool,N}}, indices::Base.RangeIndex...) where N
    tf = HDF5._getindex(dset.plain, UInt8, indices...)
    bool(tf)
end
_getindex(dset::JldDataset, ::Type{Array{T,N}}, indices::Union{Integer, Base.RangeIndex}...) where {T,N} =
    getrefs(dset, T, indices...)
function setindex!(dset::JldDataset, X::Array, indices::Base.RangeIndex...)
    if !exists(attrs(dset.plain), name_type_attr)
        # Fallback to plain read
        return setindex!(dset.plain, X, indices...)
    end
    # Read the type
    typename = a_read(dset.plain, name_type_attr)
    if typename == "Tuple"
        return read_tuple(dset, indices...)
    end
    # Convert to Julia type
    T = julia_type(typename)
    HDF5._setindex!(dset, T, X, indices...)
end

length(x::Union{JldFile, JldGroup}) = length(names(x))

### Dump ###
function dump(io::IO, parent::Union{JldFile, JldGroup}, n::Int, indent)
    nms = names(parent)
    println(io, typeof(parent), " len ", length(nms))
    if n > 0
        i = 1
        for k in nms
            print(io, indent, "  ", k, ": ")
            v = parent[k]
            if isa(v, HDF5Group)
                dump(io, v, n-1, string(indent, "  "))
            else
                if exists(attrs(v.plain), name_type_attr)
                    typename = a_read(v.plain, name_type_attr)
                    if length(typename) >= 5 && (typename[1:5] == "Array" || typename[1:5] == "Tuple")
                        println(io, typename, " ", size(v))
                    else
                        println(io, typename)
                    end
                else
                    dump(io, v, 1, indent)
                end
            end
            close(v)
            if i > n
                println(io, indent, "  ...")
                break
            end
            i += 1
        end
    end
end

### Converting attribute strings to Julia types

const _typedict = Dict{String,Type}()
_typedict["CompositeKind"] = CompositeKind
function _julia_type(s::AbstractString)
    typ = get(_typedict, s, UnconvertedType)
    if typ == UnconvertedType
        e = Meta.parse(s)
        e = JLD.fixtypes(e)
        typ = UnsupportedType
        if JLD.is_valid_type_ex(e)
            try     # try needed to catch undefined symbols
                typ = eval(e)
                if !isa(typ, Type)
                    typ = UnsupportedType
                end
            catch
                try
                    typ = Core.eval(Main, e)
                catch
                    typ = UnsupportedType
                    if !isa(typ, Type)
                        typ = UnsupportedType
                    end
                end
            end
        else
            typ = UnsupportedType
        end
        if typ != UnsupportedType
            _typedict[s] = typ
        end
    end
    typ
end

### Converting Julia types to fully qualified names
full_typename(jltype::Union) = @sprintf "Union(%s)" join(map(full_typename, jltype.types), ",")
function full_typename(tv::TypeVar)
    if is(tv.lb, Union{}) && is(tv.ub, Any)
        "TypeVar(:$(tv.name))"
    elseif is(tv.lb, Union{})
        "TypeVar(:$(tv.name),$(full_typename(tv.ub)))"
    else
        "TypeVar(:$(tv.name),$(full_typename(tv.lb)),$(full_typename(tv.ub)))"
    end
end
full_typename(jltype::Tuple{Vararg{Type}}) =
    length(jltype) == 1 ? @sprintf("(%s,)", full_typename(jltype[1])) :
                          @sprintf("(%s)", join(map(full_typename, jltype), ","))
full_typename(x) = string(x)
function full_typename(jltype::DataType)
    #tname = "$(jltype.name.module).$(jltype.name)"
    tname = string(jltype.name.module, ".", jltype.name.name)  # NOTE: performance bottleneck
    if isempty(jltype.parameters)
        tname
    else
        @sprintf "%s{%s}" tname join([full_typename(x) for x in jltype.parameters], ",")
    end
end

### Version number utilities
versionnum(v::AbstractString) = map(int, split(v, '.'))
versionstring(v::Array{Int}) = join(v, '.')
function isversionless(l::Array{Int}, r::Array{Int})
    len = min(length(l), length(r))
    for i = 1:len
        if l[i] < r[i]
            return true
        end
    end
    if length(r) > len
        for i = len+1:length(r)
            if r[i] > 0
                return true
            end
        end
    end
    false
end

function names(parent::Union{JldFile, JldGroup})
    n = names(parent.plain)
    keep = trues(length(n))
    reserved = [pathrefs[2:end], pathtypes[2:end], pathrequire[2:end]]
    for i = 1:length(n)
        if in(n[i], reserved)
            keep[i] = false
        end
    end
    n[keep]
end

macro save(filename, vars...)
    if isempty(vars)
        # Save all variables in the current module
        writeexprs = Vector{Expr}(undef, 0)
        m = @__MODULE__
        for vname in names(m)
            s = string(vname)
            if !ismatch(r"^_+[0-9]*$", s) # skip IJulia history vars
                v = Core.eval(m, vname)
                if !isa(v, Module)
                    push!(writeexprs, :(if !isa($(esc(vname)), Function) write(f, $s, $(esc(vname))) end))
                end
            end
        end
    else
        writeexprs = Vector{Expr}(undef, length(vars))
        for i = 1:length(vars)
            writeexprs[i] = :(write(f, $(string(vars[i])), $(esc(vars[i]))))
        end
    end
    Expr(:block,
         :(local f = jldopen($(esc(filename)), "w")),
         Expr(:try, Expr(:block, writeexprs...), false, false,
              :(close(f))))
end

macro load(filename, vars...)
    if isempty(vars)
        if isa(filename, Expr)
            filename = Core.eval(@__MODULE__, filename)
        end
        # Load all variables in the top level of the file
        readexprs = Vector{Expr}(undef, 0)
        vars = Vector{Expr}(undef, 0)
        f = jldopen(filename)
        nms = names(f)
        for n in nms
            obj = f[n]
            if isa(obj, JldDataset)
                sym = esc(Symbol(n))
                push!(readexprs, :($sym = read($f, $n)))
                push!(vars, sym)
            end
        end
        return Expr(:block,
                    Expr(:global, vars...),
                    Expr(:try,  Expr(:block, readexprs...), false, false,
                         :(close($f))),
                    Symbol[v.args[1] for v in vars]) # "unescape" vars
    else
        readexprs = Vector{Expr}(undef, length(vars))
        for i = 1:length(vars)
            readexprs[i] = :($(esc(vars[i])) = read(f, $(string(vars[i]))))
        end
        return Expr(:block,
                    Expr(:global, map(esc, vars)...),
                    :(local f = jldopen($(esc(filename)))),
                    Expr(:try,  Expr(:block, readexprs...), false, false,
                         :(close(f))),
                    Symbol[v for v in vars]) # vars is a tuple
    end
end

# Save all the key-value pairs in the dict as top-level variables of the JLD
function save(filename::AbstractString, dict::AbstractDict)
    jldopen(filename, "w") do file
        for (k,v) in dict
            write(file, String(k), v)
        end
    end
end
# Or the names and values may be specified as alternating pairs
function save(filename::AbstractString, name::AbstractString, value, pairs...)
    if isodd(length(pairs)) || !isa(pairs[1:2:end], (AbstractString...))
        throw(ArgumentError("arguments must be in name-value pairs"))
    end
    jldopen(filename, "w") do file
        write(file, String(name), value)
        for i = 1:2:length(pairs)
            write(file, String(pairs[i]), pairs[i+1])
        end
    end
end

# load with just a filename returns a dictionary containing all the variables
function load(filename::AbstractString)
    jldopen(filename, "r") do file
        Dict{String,Any}([(var, read(file, var)) for var in names(file)])
    end
end
# When called with explicitly requested variable names, return each one
function load(filename::AbstractString, varname::AbstractString)
    jldopen(filename, "r") do file
        read(file, varname)
    end
end
load(filename::AbstractString, varnames::AbstractString...) = load(filename, varnames)
function load(filename::AbstractString, varnames::Tuple{Vararg{AbstractString}})
    jldopen(filename, "r") do file
        map((var)->read(file, var), varnames)
    end
end

function addrequire(file::JldFile, filename::AbstractString)
    files = read(file, pathrequire)
    push!(files, filename)
    o_delete(file, pathrequire)
    write(file, pathrequire, files)
end

export
    addrequire,
    ismmappable,
    jldopen,
    o_delete,
    readmmap,
    readsafely,
    @load,
    @save,
    load,
    save
end


text:
# vect: Not exported
# oldstyle_vcat_warning: Not exported
# size: interfaces.jl
# eltype: interfaces.jl
# elsize: Not exported
# ndims: interfaces.jl
# length: interfaces.jl
# endof: interfaces.jl
# first: compound using start/next
# last: compound using getindex/endof
# stride: interfaces.jl
# strides: interfaces.jl
# isassigned: interfaces.jl
# trailingsize: Not exported
# linearindexing: Not exported
# checkbounds: interfaces.jl
# throw_boundserror: Not exported
# _internal_checkbounds: Not exported
# similar: interfaces.jl
# reshape
@primitive reshape(x,i...),dy  reshape(dy,size(x))
addtest(reshape,rand(2,2),(4,1))
# copy!: Overwriting operation
# copy: interfaces.jl
# copy_transpose!: Not exported
# zero: interfaces.jl
# start: interfaces.jl
# next: interfaces.jl
# done: interfaces.jl
# eachindex: interfaces.jl
# _maxlength: Not exported
# isempty: interfaces.jl
# convert: Cannot support.
# full
@primitive full(x),dy dy
# map: Cannot support.
# pointer: interfaces.jl
# getindex: interfaces.jl
# unsafe_getindex: Not exported
# _getindex: Not exported
# _unsafe_getindex: Not exported
# setindex!: interfaces.jl, not supported
# unsafe_setindex!: Not exported
# _setindex!: Not exported
# _unsafe_setindex!: Not exported
# get (getindex with a default value)
# This can be left as a composite function, it will get its gradient from getindex if necessary.
get{T<:AbstractArray}(A::Value{T}, i::Integer, default) = (if checkbounds(Bool, length(A), i); A[i]; else; default; end)
get{T<:AbstractArray}(A::Value{T}, I::Tuple{}, default) = similar(A, typeof(default), 0)
get{T<:AbstractArray}(A::Value{T}, I::Dims, default)    = (if checkbounds(Bool, size(A), I...); A[I...]; else; default; end)
# get!: Overwriting function
# promote_eltype: Not exported

# cat(dims,A...): Concatenate the input arrays along the specified
# dimensions in the iterable dims. For dimensions not in dims, all
# input arrays should have the same size, which will also be the size
# of the output array along that dimension. For dimensions in dims,
# the size of the output array is the sum of the sizes of the input
# arrays along that dimension. If dims is a single number, the
# different arrays are tightly stacked along that dimension. If dims
# is an iterable containing several dimensions, this allows one to
# construct block diagonal matrices and their higher-dimensional
# analogues by simultaneously increasing several dimensions for every
# new input array and putting zero blocks elsewhere. For example,
# cat([1,2], matrices...) builds a block diagonal matrix, i.e. a block
# matrix with matrices[1], matrices[2], ... as diagonal blocks and
# matching zero blocks away from the diagonal.

# After catdims, cat can take 0 or more arguments of any type.  We
# only catch the cases where one of the first two args is a Value.

typealias CatDims Union{Int,Tuple,Array}
cat_r = recorder(cat)
cat(i::CatDims,a::Value,b::Value,c...)=cat_r(i,a,b,c...)
cat(i::CatDims,a,b::Value,c...)=cat_r(i,a,b,c...)
cat(i::CatDims,a::Value,b...)=cat_r(i,a,b...)
cat{N}(::Type{Grad{N}},dy,y,i,x...)=uncat(dy,N-1,i,x...)   # N-1 because first arg is catdims

# We need to extract the n'th block from dy which has the same shape
# as y=cat(catdims,x...).
# TODO: make uncat a primitive and define its gradient for higher order gradients.
function uncat(dy,n,catdims,x...)
    idx = []
    for d=1:ndims(dy)
        if in(d,catdims)
            pos = 0
            for j=1:n-1
                pos += size(x[j],d)
            end
            push!(idx,(pos+1):(pos+size(x[n],d)))
        else
            push!(idx,1:size(dy,d))
        end
    end
    dx = dy[idx...]
    if isa(getval(x[n]),AbstractArray)
        dx = reshape(dx, size(x[n]))
    else
        length(dx)==1 || error("uncat mismatch")
        dx = dx[1]
    end
    return dx
end

# TODO: add some test cases for cat.

# Same deal with vcat and hcat, catch if one of the first two args is
# a Value.  TODO: these only work for vectors and matrices, do general arrays.

# vcat: vcat(X...) = cat(1, X...)
vcat_r = recorder(vcat)
vcat(a::Value,b::Value,c...)=vcat_r(a,b,c...)
vcat(a,b::Value,c...)=vcat_r(a,b,c...)
vcat(a::Value,b...)=vcat_r(a,b...)
vcat{N}(::Type{Grad{N}},dy::Value,y,x...)=dy[xrange(x,1,N),:] # ambiguity fix
vcat{N}(::Type{Grad{N}},dy,y,x...)=dy[xrange(x,1,N),:]
addtest(vcat, rand(1,2), rand(2,2))

# hcat: hcat(X...) = cat(2, X...)
hcat_r = recorder(hcat)
hcat(a::Value,b::Value,c...)=hcat_r(a,b,c...)
hcat(a,b::Value,c...)=hcat_r(a,b,c...)
hcat(a::Value,b...)=hcat_r(a,b...)
hcat{N}(::Type{Grad{N}},dy::Value,y,x...)=dy[:,xrange(x,2,N)] # ambiguity fix
hcat{N}(::Type{Grad{N}},dy,y,x...)=dy[:,xrange(x,2,N)]
addtest(hcat, rand(2,2), rand(2,1))

function xrange(x, d, n)
    s = 0
    for i=1:n-1; s += size(x[i],d); end
    return (s+1):(s+size(x[n],d))
end

# typed_vcat: Not exported
# typed_hcat: Not exported
# cat_t: Not exported
# hvcat: TODO
# hvcat_fill: Not exported
# typed_hvcat: Not exported
# isequal: interfaces.jl
# lexcmp: interfaces.jl
# ==: interfaces.jl
# sub2ind: Not an array op.
# _sub2ind: Not exported
# ind2sub: Not an array op.
# ind2sub!: Not an array op.
# mapslices: Cannot support.
# promote_to!: Not exported
# map_promote: Not exported
# map!: Cannot support.
# map_to!: Not exported
# ith_all: Not exported
# map_n!: Not exported
# map_to_n!: Not exported
# push!: Overwriting function
# unshift!: Overwriting function
# hash: Works for Value.


text:
# This file is a part of KCenters.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using SimilaritySearch
using CategoricalArrays, StatsBase, MLDataUtils
export enet, dnet, kcenters, associate_centroids, ClusteringData

"""
    struct ClusteringData{DataType<:AbstractVector}
        # n elements in the dataset, m centers
        centers::DataType # centers, m entries
        freqs::Vector{Int32} # number of elements associated to each center, m entries
        dmax::Vector{Float32} # stores the distant element associated to each center, m entries
        codes::Vector{Int32} # id of the associated center, n entries
        distances::Vector{Float32} # from each element to its nearest center (label), n entries
        err::Vector{Float32} # dynamic of the error function, at least one entry
    end

The datastructure output of our clustering procedures
"""
struct ClusteringData{DataType<:AbstractVector}
    # n elements in the dataset, m centers
    centers::DataType # centers, m entries
    freqs::Vector{Int32} # number of elements associated to each center, m entries
    dmax::Vector{Float32} # stores the distant element associated to each center, m entries
    codes::Vector{Int32} # id of the associated center, n entries
    distances::Vector{Float32} # from each element to its nearest center (label), n entries
    err::Vector{Float32} # dynamic of the error function, at least one entry
end

"""
    kcenters(dist::PreMetric, X::AbstractVector{T}, y::CategoricalArray, sel::AbstractCenterSelection=CentroidSelection()) where T

Computes a center per region (each region is defined by the set of items having the same label in `y`).
The output is compatible with `kcenters` function when `eltype(y)` is Int
"""
function kcenters(dist::PreMetric, X::AbstractVector{T}, y::CategoricalArray, sel::AbstractCenterSelection=CentroidSelection()) where T
    m = length(levels(y))
    centers = Vector{T}(undef, m)
    freqs = zeros(Int32, m)
    invindex = labelmap(y.refs)

    for i in 1:m
        elements = @view X[invindex[i]]
        centers[i] = center(sel, elements)
        freqs[i] = length(elements)
    end

    distances = Float32[evaluate(dist, X[i], centers[y.refs[i]]) for i in eachindex(X)]
    codes = Int32.(y.refs)
    ClusteringData(centers, freqs, compute_dmax(m, codes, distances), codes, distances, Float32[sum(distances)])
end

function compute_dmax(m, codes, distances)
    dmax = zeros(Float32, m)
    for i in eachindex(codes)
        code = codes[i]
        d = distances[i]
        dmax[code] = max(dmax[code], d)
    end
    dmax
end

"""
    kcenters(dist::PreMetric, X::AbstractVector{T}, k::Integer; sel::AbstractCenterSelection=CentroidSelection(), initial=:fft, maxiters=0, tol=0.001, recall=1.0) where T
    kcenters(dist::PreMetric, X::AbstractVector{T}, C::AbstractzVector{T}; sel::AbstractCenterSelection=CentroidSelection(), maxiters=30, tol=0.001, recall=1.0) where T

Performs a kcenters clustering of `X` using `dist` as distance function and `sel` to compute center objects.
It is based on the Lloyd's algorithm yet using different algorithms as initial clusters:

- `:fft` the _farthest first traversal_ selects a set of farthest points among them to serve as cluster seeds.
- `:dnet` the _density net_ algorithm selects a set of points following the same distribution of the datasets; in contrast with a random selection, `:dnet` ensures that the selected points are not ``\\lfloor n/k \\rfloor`` nearest neighbors.
- `:sfft` the `:fft` over a ``k + \\log n`` random sample
- `:sdnet` the `:dnet` over a ``k + \\log n`` random sample
- `:rand` selects the set of random points along the dataset.

If recall is 1.0 then an exhaustive search is made to find associations of each item to its nearest cluster; if ``0 < recall < 0`` then an approximate index
(`SearchGraph` from `SimilaritySearch.jl`) will be used for the same purpose; the `recall` controls the expected search quality (trade with search time).
"""
function kcenters(dist::PreMetric, X::AbstractVector{T}, k::Integer; sel::AbstractCenterSelection=CentroidSelection(), initial=:fft, maxiters=10, tol=0.001, recall=1.0, verbose=false) where T
    if initial == :fft
        m = 0
        irefs = enet(dist, X, k+m).irefs
        if m > 0
            irefs = irefs[1+m:end]
        end
        initial = X[irefs]
    elseif initial == :dnet
        irefs = dnet(dist, X, k).irefs
        resize!(irefs, k)
        initial = X[irefs]
    elseif initial == :sfft
        n = length(X)
        m = min(n, ceil(Int, sqrt(n)) + k)
        X_ = X[unique(rand(1:n, m))]
        C = enet(dist, X_, k, verbose=verbose)
        initial = X_[C.irefs]
    elseif initial == :sdnet
        n = length(X)
        m = min(n, ceil(Int, sqrt(n)) + k)
        X_ = X[unique(rand(1:n, m))]
        irefs = dnet(dist, X_, k, verbose=verbose).irefs
        resize!(irefs, k)
        initial = X_[irefs]
    elseif initial == :fftdensity
        n = length(X)
        m = min(n, ceil(Int, log(n)) + 2 * k)

        E = enet(dist, X, m)

        C = Dict{Int, Int}()
        # D = Dict{Int, Float64}()
        for p in E.seq
            i = argmin(p)
            C[i] = get(C, i, 0) + 1
            # d = maximum(p)
            # D[i] = max(get(D, i, 0.0), d)
        end

        irefs = [s[1] for s in sort!(collect(C), by=p->p[2], rev=true)[1:2*k]]
        XX = X[irefs]
        C = enet(dist, XX, k)
        initial = XX[C.irefs]
    elseif initial == :rand
        initial = rand(X, k)
    elseif initial isa Symbol
        error("Unknown kind of initial value $initial")
    else
        initial = initial::AbstractVector{T}
    end

    kcenters(dist, X, initial, sel=sel, maxiters=maxiters, tol=tol, recall=recall, verbose=verbose)
end

function kcenters(dist::PreMetric, X::AbstractVector{T}, C::AbstractVector{T}; sel::AbstractCenterSelection=CentroidSelection(), maxiters=-1, tol=0.001, recall=1.0, verbose=true) where T
    # Lloyd's algoritm
    n = length(X)
    numcenters = length(C)

    if maxiters == -1
        maxiters = ceil(Int, log2(n))
    end

    function create_index(CC)
        if recall >= 1.0
            ExhaustiveSearch(dist, CC)
        else
            idx = SearchGraph(; dist)
            append!(idx, CC)
        end
    end

    freqs = zeros(Int32, numcenters)
    codes = Vector{Int32}(undef, n)
    distances = zeros(Float32, n)
    err = Float32[typemax(Float32), associate_centroids_and_compute_error!(X, create_index(C), codes, distances, freqs)]
    iter = 0

    while iter < maxiters && err[end-1] - err[end] >= tol
        iter += 1
        verbose && println(stderr, "*** starting iteration: $iter; err: $err ***")
        clusters = [Int[] for i in 1:numcenters]
        for (objID, plist) in enumerate(codes)
            for refID in plist
                push!(clusters[refID], objID)
            end
        end

        verbose && println(stderr, "*** computing centroids ***")
        Threads.@threads for i in 1:length(clusters)
            plist = clusters[i]
            # C[i] can be empty because we could be using approximate search
            if length(plist) > 0
                C[i] = center(sel, X[plist])
            end
        end

        verbose && println(stderr, "*** computing $(numcenters) nearest references ***")
        s = associate_centroids_and_compute_error!(X, create_index(C), codes, distances, freqs)
        push!(err, s)
        isnan(err[end]) && error("ERROR invalid score $err")
        verbose && println(stderr, "*** new score with $(numcenters) references: $err ***")
    end

    verbose && println(stderr, "*** finished computation of $(numcenters) references, err: $err ***")
    ClusteringData(C, freqs, compute_dmax(numcenters, codes, distances), codes, distances, err)
end

function associate_centroids_and_compute_error!(X, index::AbstractSearchContext, codes, distances, counters)
    Threads.@threads for objID in 1:length(X)
        #for objID in 1:length(X)
        res = KnnResult(1)
        search(index, X[objID], res)
        codes[objID] = argmin(res)
        distances[objID] = maximum(res)
    end

    for refID in codes
        counters[refID] += 1
    end

    mean(distances)
end

"""
    associate_centroids(dist::PreMetric, X, centers)

Returns the named tuple `(codes=codes, freqs=freqs, distances=distances, err=s)` where codes contains the nearest center
index for each item in `X` under the context of the `dist` distance function. `C` is the set of centroids and
`X` the dataset of objects. `C` also can be provided as a SimilaritySearch's Index.
"""
function associate_centroids(dist::PreMetric, X, centers)
    n = length(X)
    freqs = zeros(Int, length(centers))
    codes = Vector{Int}(undef, n)
    distances = Vector{Float64}(undef, n)
    if centers isa AbstractVector
        centers = ExhaustiveSearch(dist, centers)
    end
    s = associate_centroids_and_compute_error!(X, centers, codes, distances, freqs)
    (codes=codes, distances=distances, err=s)
end


text:
using SubsetSelection, SubsetSelectionCIO
using DataFrames
include("cv.jl")

##############################################
##SOLVER FOR CROSS-VALIDATING k/λ
##############################################
function oa_cv_constraint(ℓ::LossFunction, kRange, Γmax, Y, X;
    holdout = .3,
    anticycling = false, averaging = true, gradUp = 10, maxIter = 100, δ = 1e-3,
    TimeLimit = 60.)

  Path = DataFrame(k=Float64[], indices=Array[], w=Array[],
             cuts=Integer[], γ=Float64[], Δt = Float64[], error=Float64[])

  n, p = size(X)

  train, val = split_data(n, holdout) #Split data into train and validation sets

  indInit = SubsetSelection.ind_init(Constraint(kRange[1]), p)
  αInit=SubsetSelection.alpha_init(ℓ, Y[train]) #Initialization

  for c in 1:length(kRange)
    # indInit = partial_min(Constraint(Range[c]), X[train, :], αInit, γ)
    k = kRange[c]
    γ0 = 1. *p / k / (maximum(sum(X[train,:].^2,2))*n)
    factor = 1.

    stop = false
    for inner_epoch in 1:Γmax
      if !stop
        γ = factor * γ0

        result = subsetSelection(ℓ, Constraint(k), Y[train], X[train, :],
            indInit = indInit, αInit=αInit,
            anticycling = anticycling, averaging = averaging, gradUp = gradUp, γ = γ, maxIter = maxIter, δ = δ)

        indices_cio, w_cio, Δt_cio, status, Gap, cutCount = oa_formulation(ℓ, Y[train], X[train, :],
                                  k, γ, indices0=result.indices, ΔT_max=TimeLimit)

        valError = error(ℓ, Y[val], X[val, :], indices_cio, w_cio)
        push!(Path, [k, indices_cio, w_cio, cutCount, γ, Δt_cio, valError])

        αInit = result.α[:]
        indInit = indices_cio[:]

        g_old = norm(grad_primal(ℓ, Y[train], X[train, indices_cio], 0*w_cio, 2*γ))
        g_new = norm(grad_primal(ℓ, Y[train], X[train, indices_cio], w_cio, 2*γ))
        stop = (g_new / g_old < 1e-2)

        factor *= 2
      end
    end
  end

  return Path
end


text:
# using JuMP

"""
The `def` macro is used to build other macros that can insert the same block of
Julia code into different parts of a program.
This is macro is used to generate a standard set of fields inside a model type
hierarchy. (Thanks: InfrastructureModels.jl; see http://www.stochasticlifestyle.com/type-dispatch-design-post-object-oriented-programming-julia/)
"""
macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

abstract type AbstractDieterModel end

"""
A macro for adding the basic fields to an AbstractDieterModel type definition.
"""
@def dm_fields begin
    model::JuMP.AbstractModel
    data::Dict{String,<:Any}
    sets::Dict{Symbol,Array{<:Any,1}}
    parameters::Dict{Symbol,<:Any}
    settings::Dict{Symbol,<:Any}
    results::Dict{Symbol,<:Any}

    # ref::Dict{Symbol,<:Any} # reference data
    # var::Dict{Symbol,<:Any} # JuMP variables
    # con::Dict{Symbol,<:Any} # JuMP constraint references
end

mutable struct DieterModel <: AbstractDieterModel @dm_fields end

#=
 Note: We avoid defining a custom inner constructor method for subtypes of AbstractDieterModel
 while defining the basic fields using a custom macro @dm_fields.
 If we want to enforce constraints or invariants on the fields data, we may need an inner constructor.
 This gives the option to define outer constructors as needed for AbstractDieterModel subtype.
=#

const current_path = @__DIR__

"""
Default constructor for a Dieter model, initialised with the Model DataType and
a user-specified `data` Dict associated with the model in the constructor.
"""
function InitialiseDieterModel(ModelType::Type, data::Dict{String,T} where T;
            datapath::String="",
            jump_model::JuMP.AbstractModel=JuMP.Model(),
            settings=Dict{Symbol,Any}()
         )

        @assert ModelType <: AbstractDieterModel
        # initialise data structures for the model

        sets = Dict{Symbol, Array{Any,1}}()
        parameters = Dict{Symbol, Dict}()

        # initialise settings

        default_settings = Dict{Symbol, Union{String,Number}}(
            :scen => "Default",
            :datapath => datapath,
            :interest => 0.04,
            :co2 => 0,     # Carbon Price
            :min_res => 0,  # Minimum Renewable Share (%)
            :cu_cost => 0,   # Curtailment Cost
            :ev => 0,
            :heat => 0,
            :h2 => 0
        )

        initial_settings = merge(default_settings, settings)

        # initialise results
        # results = Dict{String,Any}()
        results = Dict{Symbol,DataFrame}()

        dtr = ModelType(
            jump_model,
            data,
            sets,
            parameters,
            initial_settings,
            results
        )

    return dtr
end

# %%
"""
Initialise a basic file structure for accessing the model data.
The default directory structure is:
<datapath>/
├── base
│   ├── availability.csv
│   ├── load.csv
│   ├── storages.csv
│   └── technologies.csv
├── ev
│   ├── ev.csv
│   ├── ev__demand.csv
│   └── ev__power.csv
├── h2
│   └── h2__technologies.csv
└── heat
    ├── buildings.csv
    ├── dhw__demand.csv
    ├── heat.csv
    ├── heat__demand.csv
    ├── heat__technologies.csv
    └── temperature.csv
"""
function initialise_data_dir_structure(datapath;remove=false)
    if !(isdir(datapath))
        @info "The given input directory \"$datapath\" does not exist. It will be created."
        # return false
    end

    datadirs = ["base","ev","h2","heat"]

    for dir in datadirs
        thispath=joinpath(datapath,dir)
        if isdir(thispath)
            @info "Directory \"$dir\" in \"$datapath\"/ already exists, skipping..."
            continue
        else
            @info "Creating directory \"$dir\" in \"$datapath\"/ "
            mkpath(thispath)
        end
    end

    return true
end

"""
Initialise the strings giving the file-path references to the model data.
The argument `sfx` is the file suffix of each model data file.
Currently, `sfx` can be either "csv" or "sql"
"""
function initialise_data_file_dict!(dtr::AbstractDieterModel,sfx::String)

    datapath = dtr.settings[:datapath]

    dtr.data["files"] = Dict{String,String}()

    fileDict = dtr.data["files"]

    # %% Base
    fileDict["nodes"] = joinpath(datapath,"base","nodes."*sfx)
    fileDict["arcs"] = joinpath(datapath,"base","arcs."*sfx)
    fileDict["tech"] = joinpath(datapath,"base","technologies."*sfx)
    fileDict["storage"] = joinpath(datapath,"base","storages."*sfx)
    fileDict["load"] = joinpath(datapath,"base","load."*sfx)
    fileDict["avail"] = joinpath(datapath,"base","availability."*sfx)
    fileDict["fuel_prices"] = joinpath(datapath,"base","fuel_prices."*sfx)

    fileDict["map_node_demand"] = joinpath(datapath,"base","map_node_demand."*sfx)
    fileDict["map_node_tech"] = joinpath(datapath,"base","map_node_tech."*sfx)
    fileDict["map_node_storages"] = joinpath(datapath,"base","map_node_storages."*sfx)

    # %% EV
    if !(dtr.settings[:ev] |> ismissing)
        fileDict["ev"] = joinpath(datapath,"ev","ev."*sfx)
        fileDict["ev_demand"] = joinpath(datapath,"ev","ev_demand."*sfx)
        fileDict["ev_power"] = joinpath(datapath,"ev","ev_power."*sfx)
    end

    # %% Heat
    if !(dtr.settings[:heat] |> ismissing)
        fileDict["heat"] = joinpath(datapath,"heat","heat."*sfx)
        fileDict["heat_technologies"] = joinpath(datapath,"heat","heat_technologies."*sfx)
        fileDict["buildings"] = joinpath(datapath,"heat","buildings."*sfx)
        fileDict["temperature"] = joinpath(datapath,"heat","temperature."*sfx)
        fileDict["heat_demand"] = joinpath(datapath,"heat","heat_demand."*sfx)
        fileDict["dhw_demand"] = joinpath(datapath,"heat","dhw_demand."*sfx)
    end

    # %% Hydrogen (h2)
    if !(dtr.settings[:h2] |> ismissing)
        fileDict["h2_technologies"] = joinpath(datapath,"h2","h2_technologies."*sfx)
    end

    return dtr
end


text:
using GeneralizedDynamicsFromData
using OrderedCollections
using HDF5
using FileIO

experiment_name = "selkov_iv_study_MSE"
repetitions = 50

weight_decay = 1e-4
η = 1e-2
η_decay_rate = 0.5
η_decay_step = 800
η_limit = 1e-4
noise = 1e-4

construct_optimiser() = Flux.Optimiser(WeightDecay(weight_decay),
                                       ExpDecay(η,η_decay_rate,η_decay_step,η_limit),
                                       ADAM())
construct_loss(θ, y, predict) = mse_loss(θ, y, predict)

net_config = OrderedDict([
    :inputs => 2,
    :outputs => 1,
    :neurons => 16,
    :layers => 1,
    :non_lin => tanh,
    :initialization => Flux.glorot_normal
])

initial_values = [ Float64[1.0 ,1.0],
                   Float64[0.5 ,1.0],
                   Float64[0.0 ,1.0],
                   Float64[0.25,0.5],
                   Float64[0.25,3.0],
                   Float64[0.5 ,6.0] ]

params = [ [0.1, 0.6],    # oscillation
           [0.1, 0.15] ]   # steadystate

for (id, u0) in enumerate(initial_values)

  problem = Dict([:equation => selkov,
                  :parameters => params[1],
                  :u0 => u0,
                  :tspan => (0.0f0, 15.0f0),
                  :ts => 0.1,
                  :solver => Tsit5,
                  :optimizer => construct_optimiser,
                  :max_iter => 3000,
                  :loss => construct_loss
                 ])

  summary, _ = repeat_experiment(problem,
                                 net_config,
                                 repetitions;
                                 longterm = 60.0,
                                 ε = noise,
                                 progress=false)

  min_losses = convert(Array{Float64,1}, summary["losses"])
  longterm_predictions = cat([sol for sol in summary["longterm_predictions"]]...; dims=3)
  longterm_solution = cat([sol for sol in summary["longterm_solution"]]...; dims=3)

  fn = joinpath("./data/", experiment_name*"_oscillation_"*"$(id).h5")

  h5open(fn, "w") do file
      file["parameters"] = params[1]
      file["u0"] = u0
      file["losses"] = min_losses
      file["longterm_predictions"] = longterm_predictions
      file["longterm_solution"] = longterm_solution
  end

end

for (id, u0) in enumerate(initial_values)

  problem = Dict([:equation => selkov,
                  :parameters => params[2],
                  :u0 => u0,
                  :tspan => (0.0f0, 15.0f0),
                  :ts => 0.1,
                  :solver => Tsit5,
                  :optimizer => construct_optimiser,
                  :max_iter => 3000,
                  :loss => construct_loss
                 ])

  summary, _ = repeat_experiment(problem,
                                 net_config,
                                 repetitions;
                                 longterm = 60.0,
                                 ε = noise,
                                 progress=false)

  min_losses = convert(Array{Float64,1}, summary["losses"])
  longterm_predictions = cat([sol for sol in summary["longterm_predictions"]]...; dims=3)
  longterm_solution = cat([sol for sol in summary["longterm_solution"]]...; dims=3)

  fn = joinpath("./data/", experiment_name*"_steadystate_"*"$(id).h5")

  h5open(fn, "w") do file
      file["parameters"] = params[2]
      file["u0"] = u0
      file["losses"] = min_losses
      file["longterm_predictions"] = longterm_predictions
      file["longterm_solution"] = longterm_solution
  end

end

max_stars_repo_name: pnieters/GeneralizedDynamicsFromData

max_issues_repo_name: pnieters/GeneralizedDynamicsFromData

max_forks_repo_name: pnieters/GeneralizedDynamicsFromData


text:
# The following code is in the public domain.

using Base64

showimg(mime, fn; tag="img") = open(fn) do f
    base64 = base64encode(f)
    display("text/html", """<$tag src="data:$mime;base64,$base64" />""")
end


text:
@testset "Numerics" begin
    @test_throws ErrorException assert_approx_equal(1, 1 + 1e-5, 1e-10, 1e-6, "assertion")
    @test_throws ErrorException assert_approx_equal(1, 1 + 1e-5, 1e-10, 1e-6)
    @test assert_approx_equal(1, 1 + 1e-7, 1e-10, 1e-6, "assertion")
    @test assert_approx_equal(1, 1 + 1e-7, 1e-10, 1e-6)
    @test_throws ErrorException assert_approx_equal(0, 1e-9, 1e-10, 1e-6, "assertion")
    @test_throws ErrorException assert_approx_equal(0, 1e-9, 1e-10, 1e-6)
    @test assert_approx_equal(0, 1e-11, 1e-10, 1e-6, "assertion")
    @test assert_approx_equal(0, 1e-11, 1e-10, 1e-6)
end


text:
module menu_

using IUP

export
	menu

include("../src/libiup_h.jl")

#------------------------example html/examples/C/menu.c ----------------------------------
function menu()
	IupOpen()	#Initializes IUP

	item_open = IupItem ("Open", "");
	IupSetAttribute(item_open, "KEY", "O");

	item_save = IupItem ("Save", "");
	IupSetAttribute(item_save, "KEY", "S");

	item_undo = IupItem ("Undo", "");
	IupSetAttribute(item_undo, "KEY", "U");
	IupSetAttribute(item_undo, "ACTIVE", "NO");
	item_exit = IupItem ("Exit", "");
	IupSetAttribute(item_exit, "KEY", "x");
	IupSetCallback(item_exit, "ACTION", cfunction(exit_cb, Cint, (Ptr{Ihandle},)))

	file_menu = IupMenu(item_open,
	                  item_save,
	                  IupSeparator(),
	                  item_undo,
	                  item_exit)

	sub1_menu = IupSubmenu("File", file_menu);

	thismenu = IupMenu(sub1_menu)
	IupSetHandle("mymenu", thismenu)
	dlg = IupDialog(IupCanvas(""))
	IupSetAttribute(dlg, "MENU", "mymenu")
	IupSetAttribute(dlg, "TITLE", "IupMenu")

	IupShow(dlg)
	IupMainLoop()                               # Initializes IUP main loop
	IupClose()                                  # And close it when ready
end

function exit_cb(h::Ptr{Ihandle})
	return convert(Cint, IUP_CLOSE)
end

end 	# module