text:
module FitTools

using LinearAlgebra
using JSON
using DelimitedFiles
using Glob
using Printf

using JuLIP, ACE, ACEtb
using ACEtb.Bonds: get_env
using ACEtb.Utils: read_json, write_json, get_data
using ACEtb.Predictions: predict, train_and_predict

struct FitContext
    base_dir::String
    dftb_path::String
    output_dir::String
    param_basename::String
end

function make_params(mode::Symbol, training_datasets, test_datasets; rcut=12,
    renv=3, zenv=1, deg=5, order=2,
    envdeg=4, pcut=1, unitcell=false,
    method=:QR, rtol=1e-15)

    if mode == :fit
        predict = 1
        fit = 1
    elseif mode == :predict
        predict = 1
        fit = 0
    else
        error("mode should be either :fit or :predict")
    end

    param_dict = Dict(
        "training_datasets" => training_datasets,
        "test_datasets" => test_datasets,
        "model" => Dict(
                "predict" => predict,
                "fit" => fit,
                "cutoff_params" =>
                        Dict("rcut" => rcut,
                            "renv" => renv,
                            "zenv" => zenv,
                            "pcut" => pcut),
                "fit_params" => Dict(
                    "degree" => deg,
                    "order" => order,
                    "env_deg" => envdeg,
                    "method" => method,
                    "rtol" => rtol)),
        "onsite-terms" => Dict(
        "Al" => [
            -5.5708730114137865e+01,
            -4.0432979354580754e+00,
            -4.9611320942331338e-01,
            -2.6243905672312331e+00,
            -2.6243905672300967e+00,
            -2.6241786637179012e+00,
            -4.6415677514380915e-01,
            -4.6415677508608638e-01,
            -4.6415677490900636e-01,
            -2.8478483849360947e-01,
            -2.8478483839241397e-01,
            -2.8478483828616641e-01,
            -1.8885840374342641e-01,
            -1.8885840355491695e-01
        ]))
    return param_dict, params_to_string(param_dict)
end

function params_to_string(params::Dict)
    input_dict = Dict{String}{Any}()
    for (k, v) in params
        if k in ["training_datasets", "test_datasets", "model", "onsite-terms"]
            input_dict[k] = copy(v)
        end
    end
    delete!(input_dict["model"], "fit")
    delete!(input_dict["model"], "predict")
    return @sprintf "%x" hash(input_dict)
end

function do_fit(pot_dict::Dict, elm_names=["Al"])
    trdata = pot_dict["training_datasets"]
    cutoff_params = pot_dict["model"]["cutoff_params"]
    fit_params = pot_dict["model"]["fit_params"]
    @info "│    Fitting..."
    Bint_table, cutf_func, train_dict = train_and_predict(trdata, elm_names, cutoff_params,
                                                          fit_params; MPIproc=1)

    return train_dict, Bint_table, cutf_func
end

function load_model(fc::FitContext, param_hash::String)
    pot_file = joinpath(fc.output_dir, "pot_" * param_hash, "pot.json")
    pot_dict = read_json(pot_file)

    cutoff_params = pot_dict["model"]["cutoff_params"]
    fit_params = pot_dict["model"]["fit_params"]
    @info "│    Loading model..."
    Bint_table, cutf_func = predict(pot_dict, cutoff_params, fit_params; MPIproc=1)
    @info "│    Model is loaded."
    return pot_dict, param_hash, Bint_table, cutf_func
end

function test_error(test, Bint_table)
    errors = []
    for (i, (R0, Renv, V)) in enumerate(test)
        push!(errors, Bint_table(R0, Renv) - V)
    end
    return errors
end

function fit_model(fc::FitContext, params::Dict, training_datasets, test_datasets; band_errors=true)
    pot_dict, param_hash = make_params(:fit, training_datasets, test_datasets; params...)
    pot_dir = joinpath(fc.output_dir, "pot_" * param_hash)
    pot_file = joinpath(pot_dir, "pot.json")

    if isfile(pot_file)
        println("Skipping fit as $pot_file already present")
        pot_dict, param_hash, Bint_table, cutf_func = load_model(fc, param_hash)
    else
        isdir(pot_dir) || mkdir(pot_dir)
        println("Fitting model $param_hash in directory $pot_dir")

        train_dict, Bint_table, cutf_func = do_fit(pot_dict, ["Al"])
        merge!(pot_dict, train_dict)
        @info "Train error $(train_dict["error_train"])"
    end

    pot_dict["model"]["fit"] = 0 # disable fitting for subsequent runs

    if "test_error" ∉ keys(pot_dict) || "test_error_all" ∉ keys(pot_dict)
        # compute test error on independent set
        test_dataset = pot_dict["test_datasets"]
        test_data = ACEtb.Utils.get_data(test_dataset, cutf_func, get_env)
        pot_dict["test_error_all"] = test_error(test_data, Bint_table)
        pot_dict["test_error"] = norm(test_error_all)
        @info "Test error $(pot_dict["test_error"])"
    end

    # save so we can compute Fermi level; will be overwritten later
    open(joinpath(pot_dir, "pot.json"), "w") do json_file
        JSON.print(json_file, pot_dict, 4)
    end

    # compute and save Fermi level of the model for FCC and BCC
    if "fermi_level" ∉ keys(pot_dict)
        fcc_fermi_level = fermi_level(fc, :FCC, param_hash)
        bcc_fermi_level = fermi_level(fc, :BCC, param_hash)
        pot_dict["fermi_level"] = Dict("FCC" => fcc_fermi_level,
                                       "BCC" => bcc_fermi_level)
    end

    # include calculated Fermi level in potential file
    open(joinpath(pot_dir, "pot.json"), "w") do json_file
        JSON.print(json_file, pot_dict, 4)
    end

    # use the new model to compute FCC and BCC band errors
    band_errors && compute_band_errors(fc, param_hash)

    return pot_dict, param_hash, Bint_table, cutf_func
end

function fermi_level(fc::FitContext, phase::Symbol, param_hash::String)
    @assert phase in [:FCC, :BCC]
    println("Calculating Fermi level for $phase with params $param_hash")
    pot_dir = joinpath(fc.output_dir, "pot_" * param_hash)

    cp("dftb_in.hsd-$phase", joinpath(pot_dir, "dftb_in.hsd"), force=true)
    pot_file = joinpath(pot_dir,  "pot.json")
    cp(pot_file, joinpath(pot_dir, fc.param_basename); force=true)
    open(joinpath(pot_dir, "dftb_fermi_$phase.log"), "w") do logfile
        run(pipeline(setenv(`$(fc.dftb_path)/dftb+`, dir=pot_dir), stdout=logfile, stderr=logfile))
    end
    return parse(Float64, split(readchomp(`grep "Fermi level" $pot_dir/detailed.out`))[5])
end

function predict_bands(fc::FitContext, phase::Symbol, param_hash::String)
    @assert phase in [:FCC, :BCC]
    pot_dir = joinpath(fc.output_dir, "pot_" * param_hash)
    band_file = joinpath(pot_dir, "BANDS_$(phase).dat")
    if isfile(band_file)
        println("Skipping band calculation for $phase as $band_file already present")
        return read_bands(fc, phase, param_hash)
    end
    println("Calculating bands for $phase with params $param_hash")
    phase = String(phase)
    pot_file = joinpath(pot_dir, "pot.json")
    cp(pot_file, joinpath(pot_dir, fc.param_basename); force=true)

    # read template with correct band path definition for this phase
    lines = readlines("dftb_in.hsd-$phase-BANDS")

    # set the Fermi level to match value from model
    pot_dict = JSON.parsefile(pot_file)
    fermi_level = pot_dict["fermi_level"][phase]
    new_lines = []
    for line in lines
        m = match(r"(^\s*)FixedFermiLevel.*$", line)
        if m !== nothing
            line = m.captures[1] * "FixedFermiLevel = $(fermi_level)"
        end
        push!(new_lines, line)
    end
    open(joinpath(pot_dir, "dftb_in.hsd"), "w") do dftb_in
        write(dftb_in, join(new_lines, "\n"))
    end

    open(joinpath(pot_dir, "dftb_bands_$phase.log"), "w") do logfile
        run(pipeline(setenv(`$(fc.dftb_path)/dftb+`, dir=pot_dir), stdout=logfile, stderr=logfile))
        run(pipeline(setenv(`$(fc.dftb_path)/dp_bands band.out band`, dir=pot_dir),
                     stdout=logfile, stderr=logfile))
    end
    cp(joinpath(pot_dir, "band_tot.dat"), band_file)

    bands = readdlm(band_file)[:, 2:end]
    return fermi_level, bands
end

function read_bands(fc::FitContext, phase::Symbol, param_hash::String)
    @assert phase in [:FCC, :BCC]
    phase = String(phase)

    pot_dir = joinpath(fc.output_dir, "pot_" * param_hash)
    pot_file = joinpath(pot_dir, "pot.json")
    pot_dict = JSON.parsefile(pot_file)
    fermi_level = get(pot_dict, "fermi_level", missing)
    if fermi_level !== missing
        fermi_level = get(fermi_level, phase, missing)
    end

    band_file = joinpath(pot_dir, "BANDS_$(phase).dat")
    bands = readdlm(band_file)[:, 2:end]

    return fermi_level, bands
end

function read_FHI_bands(fc::FitContext, phase::Symbol)
    @assert phase in [:FCC, :BCC]
    index = phase == :FCC ? 1 : 2
    phase = lowercase(String(phase))
    dict = JSON.parsefile(joinpath(fc.base_dir, "FCC_BCC_BANDS.json"))
    fermi_level = dict["Bands"]["fermi"][phase]
    bands = hcat(dict["Bands"][phase]["TB_FHIaims"][index]...)'
    # undo shift to put back in absolute terms
    bands = bands .+ fermi_level
    return fermi_level, bands
end

kB = 8.6173303e-5
T = 100.0

fermi_dirac(x, μ, T) = 1. / (1 + exp((x - μ) / (kB * T)))

function band_error(fermi_level, bands, fermi_level_ref, bands_ref)
    error_k =  [sum(abs.((bands[n, :] .* fermi_dirac.(bands[n, :], fermi_level_ref, T) -
                          bands_ref[n, :] .* fermi_dirac.(bands_ref[n, :], fermi_level_ref, T))))
                          for n = 1:size(bands, 1)]
    RMSE = sqrt(sum(error_k .^ 2)/length(error_k))
    return RMSE
end

function compute_band_errors(fc::FitContext, param_hash::String)
    pot_file = joinpath(fc.output_dir, "pot_" * param_hash, "pot.json")
    pot_dict = JSON.parsefile(pot_file)

    fcc_fermi_level, fcc_bands = predict_bands(fc, :FCC, param_hash)
    bcc_fermi_level, bcc_bands = predict_bands(fc, :BCC, param_hash)

    fcc_fermi_level_FHI, fcc_bands_FHI = read_FHI_bands(fc, :FCC)
    bcc_fermi_level_FHI, bcc_bands_FHI = read_FHI_bands(fc, :BCC)

    RMSE_FCC = band_error(fcc_fermi_level, fcc_bands, fcc_fermi_level_FHI, fcc_bands_FHI)
    RMSE_BCC = band_error(bcc_fermi_level, bcc_bands, bcc_fermi_level_FHI, bcc_bands_FHI)

    band_error_dict = Dict("FCC" => RMSE_FCC,
                           "BCC" => RMSE_BCC)

    # save results in potential dict
    pot_dict["band_error"] = band_error_dict
    open(pot_file, "w") do json_file
        JSON.print(json_file, pot_dict, 4)
    end

    return band_error_dict
end

export fit_model, load_model, FitContext, read_bands, params_to_string, read_FHI_bands, test_error

end


text:
using NonlinearEigenproblems: nep_gallery, compute_Mder
using FEASTSolver

function info(Λ, X, residuals, c, r)
    inside(x) = in_contour(x, c, r)
    in_eig = Λ[inside.(Λ)]
    in_res = residuals[inside.(Λ)]
    print("number inside : ")
    println(size(Λ[inside.(Λ)])[1])
    if sum(inside.(Λ)) > 0
        in_res_conv = in_res[in_res .<= 1e-5]
        in_eig_conv = in_eig[in_res .<= 1e-5]
        print("number inside converged : ")
        println(size(in_eig_conv)[1])
        print("max res inside: ")
        println(maximum(residuals[inside.(Λ)]))
        if size(in_res_conv, 1) > 0
            print("max res inside non spurious: ")
            println(maximum(in_res_conv))
        end
    end
    # display(in_eig_conv)
end

nep = nep_gallery("nlevp_native_fiber")
T(x) = compute_Mder(nep, x)

n = size(nep, 1)
C = complex(99.5, 0.0)
R = 5.0

e, v, res = nlfeast_moments!(T, ComplexF64.(rand(Float64,n,10)), 2^4, 10, c=C, moments=10, r=R, debug=true, ϵ=10e-16, store=true)

info(e,v,res,C,R)
println(e)

ext: jl


text:
#=
Force generation and troponin model by Rice et al. Use in Li et al (2015)
=#
using Parameters
include("common.jl")
# Force generation, troponin and crossbridge (XB) parameters
@with_kw struct XBParams
    TOT_LTRPN = 0.07  # Total pool of low affinity troponin ca binding sites (mM)
    TOT_HTRPN = 0.14  # Total pool of high affinity troponin ca binding sites (mM)
    K_P_LTRPN = 100.0  # Forward constant of low-affinity troponin sites
    K_M_LTRPN = 4e-2  # Backward constant of low-affinity troponin sites
    K_P_HTRPN = 100.  # Forward constant of high-affinity troponin sites
    K_M_HTRPN = 3.3e-4  # Backward constant of high-affinity troponin sites
    K_PN_TROP = 0.04  # Troponin rate constant (permissive -> nonpermissive)
    SL = 2.15  # Sarcomere length (μm)
    F_XB = 0.05
    G_MIN = 0.1
    G_OFF = 0.01
    V_MAX_AM = 7.2E-3
    KM_ATP_AM = 0.03
    KI_ADP_AM = 0.26
    ζ = 0.1  # Force constant (N/mm²)
    G01 = G_MIN
    G12 = 2 * G_MIN
    G23 = 3 * G_MIN
    F01 = 3 * F_XB
    F12 = 10 * F_XB
    F23 = 7 * F_XB
    ΣPATHS = G01 * G12 * G23 + F01 * G12 * G23 + F01 * F12 * G23 + F01 * F12 * F23
    P1_MAX = F01 * G12 * G23 / ΣPATHS
    P2_MAX = F01 * F12 * G23 / ΣPATHS
    P3_MAX = F01 * F12 * F23 / ΣPATHS
    Φ = 1 + (2.3 - SL) / ((2.3 - 1.7) ^ 1.6)
    G01_SL = Φ * G01
    G12_SL = Φ * G12
    G23_SL = Φ * G23
    G01_OFF = Φ * G_OFF
    N_TROP = 3.5 * SL - 2.0
    K_CA_TRPN = K_M_LTRPN / K_P_LTRPN  # Activation factor for calcium of low affinity troponin sites
    K½_TRPN = _mm(1.7e-3 - 0.8e-3 * (SL - 1.7) / 0.6, K_CA_TRPN)
    FORCE_NORM_DENOM = P1_MAX + P2_MAX + P3_MAX
    FORCE_DENOM = P1_MAX + 2 * P2_MAX + 3 * P3_MAX
    V_MAX_AM_SCALE = V_MAX_AM / (F01 + F12 + F23)
end

# Contractile force (N/mm^2)
force(p1, p2, p3, n1, ζ, FORCE_DENOM) = ζ / FORCE_DENOM * (p1 + n1 + 2 * p2 + 3 * p3)

# ATP hydrolysis rate by myofibrils (mM/ms)
function vam(p0, p1, p2, atp_i, adp_i, V_MAX_AM_SCALE, F01, F12, F23, KM_ATP_AM, KI_ADP_AM)
    vAm = V_MAX_AM_SCALE * (F01 * p0 + F12 * p1 + F23 * p2) * _mm(atp_i * _mm(KI_ADP_AM, adp_i), KM_ATP_AM)
end
function vam(p0, p1, p2, atp_i, adp_i, pXB::XBParams)
    @unpack V_MAX_AM_SCALE, F01, F12, F23, KM_ATP_AM, KI_ADP_AM = pXB
    vam(p0, p1, p2, atp_i, adp_i, V_MAX_AM_SCALE, F01, F12, F23, KM_ATP_AM, KI_ADP_AM)
end

function _force_denom(p1, p2, p3, n1, DENOM)
    return (p1 + p2 + p3 + n1) / DENOM
end

# Insider rates of crossbridge(xb) system
function _d_ltr_ca(ltr_ca, ca_i, p1, p2, p3, n1, K_P_LTRPN, K_M_LTRPN, TOT_LTRPN, FORCE_NORM_DENOM)
    K_P_LTRPN * ca_i * (TOT_LTRPN - ltr_ca) - K_M_LTRPN * ltr_ca * (1 - 2 / 3 * _force_denom(p1, p2, p3, n1, FORCE_NORM_DENOM))
end
function _d_htr_ca(htr_ca, ca_i, K_P_HTRPN, K_M_HTRPN, TOT_HTRPN)
    K_P_HTRPN * ca_i * (TOT_HTRPN - htr_ca) - K_M_HTRPN * htr_ca
end

# Calculates the rates of the troponin and crossbridge(xb) system (scalar version)
function xb_system(p0, p1, p2, p3, n1, ltr_ca, htr_ca, ca_i, TOT_LTRPN, TOT_HTRPN,
                  K_PN_TROP, K½_TRPN, N_TROP, F01, F12, F23, G01_SL, G12_SL,
                  G23_SL, G01_OFF, FORCE_NORM_DENOM, K_P_LTRPN, K_P_HTRPN, K_M_HTRPN, K_M_LTRPN)
    k_np_trop = K_PN_TROP * (ltr_ca / (TOT_LTRPN * K½_TRPN)) ^ N_TROP
    n0 = one(p0) - (p0 + p1 + p2 + p3 + n1)
    v01 = F01 * p0
    v10 = G01_SL * p1
    v12 = F12 * p1
    v21 = G12_SL * p2
    v23 = F23 * p2
    v32 = G23_SL * p3
    v04 = K_PN_TROP * p0
    v40 = k_np_trop * n0
    v15 = K_PN_TROP * p1
    v51 = k_np_trop * n1
    v54 = G01_OFF * n1
    d_p0 = v10 + v40 - v01 - v04
    d_p1 = v01 + v21 + v51 - v10 - v12 - v15
    d_p2 = v12 + v32 - v21 - v23
    d_p3 = v23 - v32
    d_n1 = v15 - v51 - v54
    d_ltr_ca = _d_ltr_ca(ltr_ca, ca_i, p1, p2, p3, n1, K_P_LTRPN, K_M_LTRPN, TOT_LTRPN, FORCE_NORM_DENOM)
    d_htr_ca = _d_htr_ca(htr_ca, ca_i, K_P_HTRPN, K_M_HTRPN, TOT_HTRPN)
    return d_p0, d_p1, d_p2, d_p3, d_n1, d_ltr_ca, d_htr_ca
end

function xb_system(p0, p1, p2, p3, n1, ltr_ca, htr_ca, ca_i, pXB::XBParams)
    @unpack (TOT_LTRPN, TOT_HTRPN, K_PN_TROP, K½_TRPN, N_TROP, F01, F12, F23, G01_SL, G12_SL,
             G23_SL, G01_OFF, FORCE_NORM_DENOM, K_P_LTRPN, K_P_HTRPN, K_M_HTRPN, K_M_LTRPN) = pXB
    xb_system(p0, p1, p2, p3, n1, ltr_ca, htr_ca, ca_i, TOT_LTRPN, TOT_HTRPN,
             K_PN_TROP, K½_TRPN, N_TROP, F01, F12, F23, G01_SL, G12_SL,
             G23_SL, G01_OFF, FORCE_NORM_DENOM, K_P_LTRPN, K_P_HTRPN, K_M_HTRPN, K_M_LTRPN)
end

# Calcium flux to troponin (mM/ms)
jtrpn(d_ltr_ca, d_htr_ca) = d_ltr_ca + d_htr_ca


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

using Markdown
using InteractiveUtils

# ╔═╡ bbfa7498-30bf-4264-b4d9-ee0c238c9984
using Plots, Combinatorics, PlutoTest

# ╔═╡ f0d31c21-f866-4d51-86ec-d473147f1677
using LightGraphs, SimpleWeightedGraphs, GraphPlot

# ╔═╡ f20ff74b-16f0-4a9d-9969-2d240ba18177
plotly()

# ╔═╡ 56cedfe0-b02c-11eb-1f56-3557d6c170bc
md"""
# [Day 9: All in a Single Night](https://adventofcode.com/2015/day/9)

Every year, Santa manages to deliver all of his presents in a single night.

This year, however, he has some new locations to visit; his elves have provided him the distances between every pair of locations. He can start and end at any two (different) locations he wants, but he must visit each location exactly once. What is the shortest distance he can travel to achieve this?

For example, given the following distances:

```
London to Dublin = 464
London to Belfast = 518
Dublin to Belfast = 141

```

The possible routes are therefore:

```
Dublin -> London -> Belfast = 982
London -> Dublin -> Belfast = 605
London -> Belfast -> Dublin = 659
Dublin -> Belfast -> London = 659
Belfast -> Dublin -> London = 605
Belfast -> London -> Dublin = 982

```

The shortest of these is `London -> Dublin -> Belfast = 605`, and so the answer is `605` in this example.

**What is the distance of the shortest route?**
"""

# ╔═╡ 1609ee1f-1d90-4889-b39d-5fe0fb6de8c4
puzzle_input = "AlphaCentauri to Snowdin = 66
AlphaCentauri to Tambi = 28
AlphaCentauri to Faerun = 60
AlphaCentauri to Norrath = 34
AlphaCentauri to Straylight = 34
AlphaCentauri to Tristram = 3
AlphaCentauri to Arbre = 108
Snowdin to Tambi = 22
Snowdin to Faerun = 12
Snowdin to Norrath = 91
Snowdin to Straylight = 121
Snowdin to Tristram = 111
Snowdin to Arbre = 71
Tambi to Faerun = 39
Tambi to Norrath = 113
Tambi to Straylight = 130
Tambi to Tristram = 35
Tambi to Arbre = 40
Faerun to Norrath = 63
Faerun to Straylight = 21
Faerun to Tristram = 57
Faerun to Arbre = 83
Norrath to Straylight = 9
Norrath to Tristram = 50
Norrath to Arbre = 60
Straylight to Tristram = 27
Straylight to Arbre = 81
Tristram to Arbre = 90"

# ╔═╡ 3ffc6481-2267-4d97-8065-d2fd78c637e9
function build_graph(paths)
	edges = []
	for line in split(paths, "\n")
		m = match(r"(\w+) to (\w+) = (\d+)", line)
		push!(edges, (m[1], m[2], parse(Int, m[3])))
	end

	nodes = Dict()
	next_id = 1
	for (from, to, weight) in edges
		if !haskey(nodes, from)
			nodes[from] = next_id
			next_id += 1
		end
		if !haskey(nodes, to)
			nodes[to] = next_id
			next_id += 1
		end
	end

	g = SimpleWeightedGraph(length(nodes));

	node_labels = []
	edge_labels = []

	for (from, to, weight) in edges
		push!(edge_labels, string(weight))

		if from ∉ node_labels
			push!(node_labels, from)
		end
		if to ∉ node_labels
			push!(node_labels, to)
		end

		add_edge!(g, nodes[from], nodes[to], weight);
	end

	g, node_labels, edge_labels
end

# ╔═╡ 78d4f47f-96b6-44f9-b9ab-2cb84bbb19b0
example_graph, example_nodes, example_edges = build_graph("""London to Dublin = 464
London to Belfast = 518
Dublin to Belfast = 141""")

# ╔═╡ 7e17c6f5-3972-4777-9611-f4c321108144
gplot(example_graph, nodelabel=example_nodes, edgelabel=example_edges)

# ╔═╡ 8e78e6fe-cb4a-46f7-839b-42aa35f7b71a
puzzle_graph, puzzle_nodes, puzzle_edges = build_graph(puzzle_input)

# ╔═╡ 8dde0e83-121f-4720-80b6-00020bd33a6b
gplot(puzzle_graph, nodelabel=puzzle_nodes, edgelabel=puzzle_edges)

# ╔═╡ 76183faa-4f76-4da6-bcff-255e7d70b730
function read(paths)
	nodes = []
	distances = Dict()

	for line in split(paths, "\n")
		m = match(r"(\w+) to (\w+) = (\d+)", line)
		distances[m[1] * m[2]] = distances[m[2] * m[1]] = parse(Int, m[3])

		if m[1] ∉ nodes
			push!(nodes, m[1])
		end
		if m[2] ∉ nodes
			push!(nodes, m[2])
		end
	end

	nodes, distances
end

# ╔═╡ 8abefd65-5a57-4717-8bf4-4e92ebe9d562
function smallest((nodes, distances),)
	smallest_distance = Inf

	for c in permutations(nodes)
		d = 0
		for idx in 1:length(c)-1
			m = c[idx]
			n = c[idx+1]
			name = m * n
			d += distances[name]
		end

		if d < smallest_distance
			smallest_distance = d
		end
	end

	smallest_distance
end

# ╔═╡ 684af9bd-4282-40ea-85d8-2d2d89c258c0
@assert smallest(read("""London to Dublin = 464
London to Belfast = 518
Dublin to Belfast = 141""")) == 605

# ╔═╡ be574368-3f0d-4863-b420-ad14106b5175
part1 = smallest(read(puzzle_input))

# ╔═╡ a1dafb48-4984-431a-82de-d5280f30af21
@test part1 == 141

# ╔═╡ 1d1f08c8-7092-4ea9-8884-43e28062c036
md"Your puzzle answer was `141`."

# ╔═╡ 64028478-8bc7-412a-82be-de768aaa1fe6
md"""
# Part Two

The next year, just to show off, Santa decides to take the route with the longest distance instead.

He can still start and end at any two (different) locations he wants, and he still must visit each location exactly once.

For example, given the distances above, the longest route would be `982` via (for example) `Dublin -> London -> Belfast`.

**What is the distance of the longest route?**
"""

# ╔═╡ ccd1c5e7-1c13-4608-8293-e4e0a76cca37
function biggest((nodes, distances),)
	biggest_distance = 0

	for c in permutations(nodes)
		d = 0
		for idx in 1:length(c)-1
			m = c[idx]
			n = c[idx+1]
			name = m * n
			d += distances[name]
		end

		if d > biggest_distance
			biggest_distance = d
		end
	end

	biggest_distance
end

# ╔═╡ a548a1da-187a-4aa7-a338-c6c0b4114177
@assert biggest(read("""London to Dublin = 464
London to Belfast = 518
Dublin to Belfast = 141""")) == 982

# ╔═╡ dbdedef1-f657-460f-93a7-e27a2563fcae
part2 = biggest(read(puzzle_input))

# ╔═╡ 7d51e278-0442-4967-9bba-0cbc1518c693
@test part2 == 736

# ╔═╡ 0d5b2494-5454-41ec-9a3d-e627e4c27f51
md"Your puzzle answer was `736`."

# ╔═╡ Cell order:
# ╠═bbfa7498-30bf-4264-b4d9-ee0c238c9984
# ╠═f0d31c21-f866-4d51-86ec-d473147f1677
# ╠═f20ff74b-16f0-4a9d-9969-2d240ba18177
# ╟─56cedfe0-b02c-11eb-1f56-3557d6c170bc
# ╟─1609ee1f-1d90-4889-b39d-5fe0fb6de8c4
# ╠═3ffc6481-2267-4d97-8065-d2fd78c637e9
# ╠═78d4f47f-96b6-44f9-b9ab-2cb84bbb19b0
# ╠═7e17c6f5-3972-4777-9611-f4c321108144
# ╠═8e78e6fe-cb4a-46f7-839b-42aa35f7b71a
# ╠═8dde0e83-121f-4720-80b6-00020bd33a6b
# ╠═76183faa-4f76-4da6-bcff-255e7d70b730
# ╠═8abefd65-5a57-4717-8bf4-4e92ebe9d562
# ╠═684af9bd-4282-40ea-85d8-2d2d89c258c0
# ╠═be574368-3f0d-4863-b420-ad14106b5175
# ╠═a1dafb48-4984-431a-82de-d5280f30af21
# ╟─1d1f08c8-7092-4ea9-8884-43e28062c036
# ╟─64028478-8bc7-412a-82be-de768aaa1fe6
# ╠═ccd1c5e7-1c13-4608-8293-e4e0a76cca37
# ╠═a548a1da-187a-4aa7-a338-c6c0b4114177
# ╠═dbdedef1-f657-460f-93a7-e27a2563fcae
# ╠═7d51e278-0442-4967-9bba-0cbc1518c693
# ╟─0d5b2494-5454-41ec-9a3d-e627e4c27f51


text:
# benchmarking several functions on views
using BenchmarkTools
using SparseArrays

import SparseArrays:   blockdiag, nnz, findnz, ftranspose
import Base:    isone

if !hasmethod(blockdiag, (AbstractMatrix,))
    blockdiag(A::AbstractMatrix...) = blockdiag(sparse.(A)...)
end
if !hasmethod(nnz, (AbstractMatrix,))
    nnz(A::AbstractMatrix) = nnz(sparse(A))
end
if !hasmethod(findnz, (AbstractMatrix,))
    findnz(A::AbstractMatrix) = findnz(sparse(A))
end
if !hasmethod(ftranspose, (AbstractMatrix,Any))
    ftranspose(A::AbstractMatrix, f) = ftranspose(sparse(A), f)
end
isone(A::AbstractMatrix) = isone(sparse(A))

function view_benchmarks(n::Integer=1000)

    x_findall(x) = findall(Base.Fix2(isless, 0.0), x)
    x_ftranspos(x) = ftranspose(x, transpose)
    fname(f) = string(nameof(f))

    suite = BenchmarkGroup()

    rng = MersenneTwister(20190131)
    res = []
    n1, n2 = (n ÷ 10) + 1, n - n ÷ 10
    B = []
    A = sprand(rng, n, n, 0.01)
    push!(B, ("CSC", A))
    push!(B, ("CSCView-slice", view(A, :, n1:n2)))
    push!(B, ("CSCInterface-1", view(A, 1:n, n1:n2)))
    push!(B, ("CSCInterface-2", view(A, n1:n2, n2:-1:n1)))
    push!(B, ("CSCView-rev", view(A, n2:-1:n1, n1:n2)))
    N = axes(B, 1)

    for f in (nnz, findnz, x_ftranspos, Matrix )
        suite[fname(f)] = BenchmarkGroup(["views", string(n)])
        for i in N
            name, X = B[i]
            suite[fname(f)][name] = @benchmarkable $f($X)
        end
    end

    for f in (x_findall, iszero, isone)
        suite[fname(f)] = BenchmarkGroup(["views", string(n)])
        for i in N
            name, X = B[i]
            suite[fname(f)][name] = @benchmarkable $f($X)
        end
    end

    for f in (mapreduce,)
        suite[fname(f)] = BenchmarkGroup(["views", string(n)])
        for i in N
            name, X = B[i]
            suite[fname(f)][name*"dims1"] = @benchmarkable $f(floor, +, $X, dims=1)
            suite[fname(f)][name*"dims2"] = @benchmarkable $f(floor, +, $X, dims=2)
        end
    end

    for f in (hcat, vcat, blockdiag)
        suite[fname(f)] = BenchmarkGroup(["views", string(n)])
        for i in N
            name, X = B[i]
            suite[fname(f)][name] = @benchmarkable $f($X, $X)
        end
    end

    AS = A + A'
    AS = [AS AS]
    B = []
    push!(B, ("Symm", AS))
    push!(B, ("SymmView-slice", view(AS, :, 1:n)))
    push!(B, ("SymmInterface-2", view(AS, n1:n2, n1:n2)))
    push!(B, ("SymmView-rev", view(AS, n2:-1:n1, n2:-1:n1)))
    N = axes(B, 1)
    for f in (issymmetric, ishermitian)
        suite[fname(f)] = BenchmarkGroup(["views", string(n)])
        for i in N
            name, X = B[i]
            suite[fname(f)][name] = @benchmarkable $f($X)
        end
    end

    ASU = triu(A)
    ASL = tril(A)
    B = []
    push!(B, ("Triu", ASU))
    push!(B, ("Triu-slice", view(ASU, :, 1:n)))
    push!(B, ("TriuInterface", view(ASU, n1:n2, n1:n2)))
    push!(B, ("Triu-rev", view(ASL, n2:-1:n1, n2:-1:n1)))
    N = axes(B, 1)
    for f in (istriu,)
        suite[fname(f)] = BenchmarkGroup(["views", string(n)])
        for i in N
            name, X = B[i]
            suite[fname(f)][name] = @benchmarkable $f($X)
        end
    end

    B = []
    push!(B, ("Tril", ASL))
    push!(B, ("Tril-slice", view(ASL, :, 1:n)))
    push!(B, ("TrilInterface", view(ASL, n1:n2, n1:n2)))
    push!(B, ("Tril-rev", view(ASU, n2:-1:n1, n2:-1:n1)))
    N = axes(B, 1)
    for f in (istril,)
        suite[fname(f)] = BenchmarkGroup(["views", string(n)])
        for i in N
            name, X = B[i]
            suite[fname(f)][name] = @benchmarkable $f($X)
        end
    end

    suite
end


text:
using Fontconfig, Luxor

function buildfontlist()
    fonts = []
    for font in Fontconfig.list()
        families = Fontconfig.format(font, "%{family}")
        for family in split(families, ",")
            push!(fonts, family)
        end
    end
    filter!(f -> !ismatch(r".LastResort|Agenda|Topaz|Bodoni Ornaments|System",
			           f), fonts)
    return sort(unique(fonts))
end

function tabulatepi()
    fonts = buildfontlist()
    ncols = 25
    nrows = convert(Int, ceil(length(fonts))) ÷ ncols
    @svg begin
        background("ivory")
        setopacity(1)
        t = Table(nrows, ncols, 30, 25)
        sethue("black")
        cellnumber = 1
        for n in 1:length(fonts)
            fontface(fonts[n])
            te = textextents("π")
            if te[3] > 0.0
                fontsize(18)
                text("π", t[cellnumber], halign=:center)
                setfont("Lucida-Sans", 3)
                settext(fonts[n], t[cellnumber] + (0, isodd(cellnumber) ? 6 : 10), halign="center")
                cellnumber += 1
            end
        end
    end 800 1200
end

tabulatepi()

max_stars_repo_name: MarsBighead/mustang

max_issues_repo_name: MarsBighead/mustang

max_forks_repo_name: MarsBighead/mustang


text:
using ArchGDAL: AbstractRasterBand, getdataset, getgeotransform, blocksize, pixeltype,
        readblock!, writeblock!, getband, getgeotransform, width, height, getnodatavalue
using ArchGDAL
using OffsetArrays

"""
Find the origin and linear transform from pixel to x-y.
"""
geo_to_transform(geo) = ([geo[1], geo[4]], [geo[2] geo[3]; geo[5] geo[6]])

"""
    geo_to_xy_corner(origin, transform, i, j)

i and j are 1-based.

In case of north up images, the GT(2) and GT(4) coefficients are zero, and the
GT(1) is pixel width, and GT(5) is pixel height. The (GT(0),GT(3)) position is
the top left corner of the top left pixel of the raster.

Note that the pixel/line coordinates in the above are from (0.0,0.0) at the top
left corner of the top left pixel to (width_in_pixels,height_in_pixels) at the
bottom right corner of the bottom right pixel. The pixel/line location of the
center of the top left pixel would therefore be (0.5,0.5).
"""
geo_to_xy_corner(origin, transform, i, j) = transform * [i - 1, j - 1] + origin

"""
    geo_to_xy_center(origin, transform, i, j)

The center of the pixel. See also: [`geo_to_xy_corner`](@ref).
"""
geo_to_xy_center(origin, transform, i, j) = (
        transform * [i - 0.5, j - 0.5] + origin
        )

"""
Which pixel contains this xy point.
"""
function pixel_containing(origin, transform, xy)
    Int.(floor.(transform \ (xy - origin))) .+ 1  # +1 for one-based indexing.
end

"""
Which pixel is closest to this xy point.
"""
function nearest_pixel(origin, transform, xy)
    Int.(round.(transform / (xy - origin))) .+ 1  # +1 for one-based indexing.
end

"""
The rect is a tuple of (ul, lr).
This tends to pull in the next pixel when
the grids are exactly aligned.
"""
function ij_cover_rect(origin, transform, rect)
    ul_ij = pixel_containing(origin, transform, rect[1])
    lr_ij = pixel_containing(origin, transform, rect[2])
    ul_ij, lr_ij
end

function pointinrect(ij, rect)
    # For the y, the up and down are switched, so it's rect[2][2].
    (rect[1][1] <= ij[1] < rect[2][1]) && (rect[2][2] <= ij[2] < rect[1][2])
end

function block_cover_ij_rect(ij_rect, blocksize)
    ((1 + (ij_rect[1][1] - 1) ÷ blocksize[1], 1 + (ij_rect[1][2] - 1) ÷ blocksize[2]),
     (1 + (ij_rect[2][1] - 1) ÷ blocksize[1], 1 + (ij_rect[2][2] - 1) ÷ blocksize[2]))
end

function corner_pixel(block_idx::CartesianIndex, blocksize)
    CartesianIndex((block_idx[1] - 1) * blocksize[1] + 1, (block_idx[2] - 1) * blocksize[2] + 1)
end

### PixelGrid methods.

"""
Represents a grid extent and its geometry, without data.

For GDAL, coordinates are (x-offset, y-offset), and it's stored as
(row, col) in the array. We think of the x-offset, visually, as the column,
but forget that. This world is (x-offset, y-offset).

The grid extent is an upper-lefthand `(i0, j0)` and a lower right-hand
`(i1, j1)`, where `i0 < i1` and `j0 < j1`. This orientation comes from
the geo tradition to flip the y-axis, so increasing y goes negative.
The `origin` and `transform` are from the GDAL geo array and represent
the origin of the grid and how each point offsets that origin.
From the GDAL documentation, the geo is used this way.
```
Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
```
"""
struct PixelGrid
    ul::Tuple{Int64,Int64}  # As x, y
    lr::Tuple{Int64,Int64}  # As x, y
    origin::Array{Float64,1}  # x, y, too.
    transform::Array{Float64,2}
end

"""
Make a function to get the geotransform so that we can stub the Band during testing.
"""
ArchGDAL.getgeotransform(band::AbstractRasterBand) = getgeotransform(getdataset(band))

"""
Make a PixelGrid that represents all data available in a Band.
"""
function pixelgrid(band::AbstractRasterBand)
    geo = getgeotransform(band)
    w = width(band)
    h = height(band)
    origin, transform = geo_to_transform(geo)
    PixelGrid((1, 1), (w, h), origin, transform)
end

"""
Compute the area of the xy intersection of two pixelgrids.

If the returned value is negative, they don't intersect.
"""
function intersection_area(a::PixelGrid, b::PixelGrid)
    # The +1 is because this gives the same upper-left corner in both cases.
    ula = geo_to_xy_corner(a.origin, a.transform, a.ul[1], a.ul[2])
    lra = geo_to_xy_corner(a.origin, a.transform, a.lr[1] + 1, a.lr[2] + 1)
    ulb = geo_to_xy_corner(b.origin, b.transform, b.ul[1], b.ul[2])
    lrb = geo_to_xy_corner(b.origin, b.transform, b.lr[1] + 1, b.lr[2] + 1)
    ul = [max(ula[1], ulb[1]), min(ula[2], ulb[2])]
    lr = [min(lra[1], lrb[1]), max(lra[2], lrb[2])]
    (lr[1] - ul[1]) * (ul[2] - lr[2])
end

function xy_bounds(pg::PixelGrid)
    ul = geo_to_xy_corner(pg.origin, pg.transform, pg.ul[1], pg.ul[2])
    lr = geo_to_xy_corner(pg.origin, pg.transform, pg.lr[1] + 1, pg.lr[2] + 1)
    ul, lr
end

function crop_to(large::PixelGrid, xybounds)
    ul_xy, lr_xy = xybounds
    ul_ij, lr_ij = ij_cover_rect(large.origin, large.transform, (ul_xy, lr_xy))
    uln = (max(large.ul[1], ul_ij[1]), max(large.ul[2], ul_ij[2]))
    lrn = (min(large.lr[1], lr_ij[1]), min(large.lr[2], lr_ij[2]))
    PixelGrid(uln, lrn, large.origin, large.transform)
end

invalid(pg::PixelGrid) = (pg.ul[1] > pg.lr[1]) || (pg.ul[2] > pg.lr[2])
valid(pg::PixelGrid) = !invalid(pg)

"""
Given a pixel grid, pull out a pixel grid for a single pixel.
"""
function onepixelgrid(pg::PixelGrid, pij)
    PixelGrid((pij[1], pij[2]), (pij[1], pij[2]), pg.origin, pg.transform)
end

"""
Iterate over all blocks that cover a given pixel rectangle.
"""
function iter_block_indices(pg::PixelGrid, blocksize)
    coarse_block_bounds = block_cover_ij_rect((pg.ul, pg.lr), blocksize)
    i, j = 1, 2
    CartesianIndices((coarse_block_bounds[1][i]:coarse_block_bounds[2][i],
            coarse_block_bounds[1][j]:coarse_block_bounds[2][j]))
end

### Stub of GDAL's AbstractRasterBand, for testing.

struct StubBand{T} <: AbstractRasterBand
    A::Array{T, 2}  # data
    bs::Array{Int32,1}  # blocksize
    o::Array{Float64,1}  # origin in xy.
    t::Array{Float64,2}  # linear transform from pixel to xy.
end
function Base.show(io::IO, b::StubBand{T}) where {T}
    print(io, "StubBand{$(T)} A=$(size(b.A)) blocksize=$(b.bs):\n\t$(b.o)+$(b.t)")
end
ArchGDAL.width(b::StubBand) = size(b.A, 1)
ArchGDAL.height(b::StubBand) = size(b.A, 2)
ArchGDAL.blocksize(b::StubBand) = b.bs
ArchGDAL.pixeltype(b::StubBand{T}) where {T} = T
function ArchGDAL.readblock!(b::StubBand, i::Integer, j::Integer, buffer)
    @assert 0 <= b.bs[1] * i < width(b)  # beginning of block is less than width.
    @assert 0 <= b.bs[2] * j < height(b)
    # The last block will likely extend past the end of the actual data.
    xaxis = (b.bs[1]*i + 1):min(b.bs[1]*(i+1), size(b.A, 1))
    yaxis = (b.bs[2]*j + 1):min(b.bs[2]*(j+1), size(b.A, 2))
    buffer[1:length(xaxis), 1:length(yaxis)] = b.A[xaxis, yaxis]
    if length(xaxis) < b.bs[1]
        buffer[(length(xaxis) + 1):b.bs[1], :] .= getnodatavalue(b)
    end
    if length(yaxis) < b.bs[2]
        buffer[:, (length(yaxis) + 1):b.bs[2]] .= getnodatavalue(b)
    end
end
ArchGDAL.getnodatavalue(b::StubBand{T}) where {T <: Real} = -floatmax(T)
ArchGDAL.getnodatavalue(b::StubBand{T}) where {T <: Integer} = typemin(T)
function ArchGDAL.writeblock!(b::StubBand, i::Integer, j::Integer, buffer)
    @assert 0 <= b.bs[1] * i < width(b)  # beginning of block is less than width.
    @assert 0 <= b.bs[2] * j < height(b)
    # The last block will likely extend past the end of the actual data.
    xaxis = (b.bs[1]*i + 1):min(b.bs[1]*(i+1), size(b.A, 1))
    yaxis = (b.bs[2]*j + 1):min(b.bs[2]*(j+1), size(b.A, 2))
    b.A[xaxis, yaxis] = buffer[1:length(xaxis), 1:length(yaxis)]
end
ArchGDAL.getgeotransform(b::StubBand) = [b.o[1], b.t[1,1], b.t[1,2], b.o[2], b.t[2,1], b.t[2,2]]

### DataGrid methods.
"""
A DataGrid represents a matrix of raster values.

The data is stored as an array of data blocks, but its extent is usually smaller than
the extent of the blocks. We can always assume that the data is complete blocks.
We keep the data, A, indexed in absolute ids, so the axis range starts with the actual
values, not usually from 1.
"""
struct DataGrid{T}
    A::AbstractArray{T,2}
    nodatavalue::T
    blocksize::Array{Int32,1}
    ulb::Tuple{Int32,Int32}  # upper left block
    lrb::Tuple{Int32,Int32}  # lower right block
    pg::PixelGrid
end

function load_data_grid(band::AbstractRasterBand, pg::PixelGrid)
    bs = blocksize(band)
    ulb, lrb = block_cover_ij_rect((pg.ul, pg.lr), bs)
    dtype = pixeltype(band)
    ul_ij = ((ulb[1] - 1) * bs[1] + 1, (ulb[2] - 1) * bs[2] + 1)
    lr_ij = (lrb[1] * bs[1], lrb[2] * bs[2])
    A = zeros(dtype, lr_ij[1] - ul_ij[1] + 1, lr_ij[2] - ul_ij[2] + 1)
    buffer = zeros(dtype, bs...)
    for bidx in CartesianIndices(((lrb[1] - ulb[1] + 1), (lrb[2] - ulb[2] + 1)))
        # readblock! is zero-based and is col, row b/c it's xaxis, yaxis.
        readblock!(band, bidx[1] - 1, bidx[2] - 1, buffer)
        # Blocksize is still x, y, so it's flipped here.
        A[((bidx[1] - 1) * bs[1] + 1):(bidx[1] * bs[1]),
            ((bidx[2] - 1) * bs[2] + 1):(bidx[2] * bs[2])] = buffer
    end
    offset = OffsetArray(A, ul_ij[1]:lr_ij[1], ul_ij[2]:lr_ij[2])
    nodata = getnodatavalue(band)
    DataGrid{dtype}(offset, nodata, bs, ulb, lrb, pg)
end

"""
Count nonzero pixels that are inside the boundsrect.
"""
function count_nonzero_pixels(dg::DataGrid, boundsrect)::Int64
    cnt = zero(Int64)
    ul = dg.pg.ul
    lr = dg.pg.lr
    # i, j are x, y for GDAL but flipped for the array.
    for j in ul[2]:lr[2]
        for i in ul[1]:lr[1]
            if dg.A[i, j] != dg.nodatavalue
                pt = geo_to_xy_center(dg.pg.origin, dg.pg.transform, i, j)
                if pointinrect(pt, boundsrect)
                    cnt += 1
                end
            end
        end
    end
    return cnt
end

function fill_nonzero_pixels!(dg::DataGrid, value, boundsrect)
    ul = dg.pg.ul
    lr = dg.pg.lr
    for j in ul[2]:lr[2]
        for i in ul[1]:lr[1]
            if dg.A[i, j] != dg.nodatavalue
                pt = geo_to_xy_center(dg.pg.origin, dg.pg.transform, i, j)
                if pointinrect(pt, boundsrect)
                    dg.A[i, j] = value
                end
            end
        end
    end
    return nothing
end

function write_datagrid!(dg::DataGrid, band)
    bs = dg.blocksize
    for j in dg.ulb[2]:dg.lrb[2]
        for i in dg.ulb[1]:dg.lrb[1]
            xrange = ((i - 1) * bs[1] + 1):(i * bs[1])
            yrange = ((j - 1) * bs[2] + 1):(j * bs[2])
            writeblock!(band, i - 1, j - 1, dg.A[xrange, yrange])
        end
    end
    return nothing
end

function load_single_block(band::AbstractRasterBand, pg::PixelGrid, block)
    bs = blocksize(band)
    dtype = pixeltype(band)
    A = zeros(dtype, bs[1], bs[2])
    readblock!(band, block[1] - 1, block[2] - 1, A)
    ul = ((block[1] - 1) * bs[1] + 1, (block[2] - 1) * bs[2] + 1)
    lr = (block[1] * bs[1], block[2] * bs[2])
    offset = OffsetArray(A, ul[1]:lr[1], ul[2]:lr[2])
    npg = PixelGrid(ul, lr, pg.origin, pg.transform)
    nodata = getnodatavalue(band)
    DataGrid{dtype}(offset, nodata, bs, Tuple(block), Tuple(block), npg)
end

function tune_band(fine_band, coarse_band, max_pixels)
    fine_pg = pixelgrid(fine_band)
    coarse_pg = pixelgrid(coarse_band)
    coarse_crop_pg = crop_to(coarse_pg, xy_bounds(fine_pg))
    coarse_blocksize = blocksize(coarse_band)

    outside_cnt = 0
    pixel_limit = 0
    coarse_buffer = zeros(Int32, coarse_blocksize...)
    for block_idx in iter_block_indices(coarse_crop_pg, coarse_blocksize)
        block = load_single_block(coarse_band, coarse_crop_pg, block_idx)
        cartesian = CartesianIndices(block.A)
        for linear_index in eachindex(block.A)
            pij = cartesian[linear_index]
            single_coarse_pg = onepixelgrid(coarse_crop_pg, pij)
            coarse_rect = xy_bounds(single_coarse_pg)
            fine_cover_pg = crop_to(fine_pg, coarse_rect)
            area = intersection_area(single_coarse_pg, fine_cover_pg)
            if valid(fine_cover_pg) && area > 1e-9
                pixel_limit += 1
                fine_grid = load_data_grid(fine_band, fine_cover_pg)
                @show fine_grid.A .> 0
                @show fine_grid.pg.ul
                @show fine_grid.pg.lr
                landscan_value = block.A[
                    single_coarse_pg.ul[1],
                    single_coarse_pg.ul[2]
                ]
                # This should be restricted to pixels under the coarse.
                @assert fine_grid.nodatavalue == -floatmax(Float64)
                pixel_cnt = count_nonzero_pixels(fine_grid, coarse_rect)
                if pixel_cnt > 0
                    rate = min(landscan_value / pixel_cnt, 100)
                    @show block_idx
                    @show pij
                    @show coarse_rect
                    @show pixel_cnt
                    @show sum(fine_grid.A .> 0)
                    fill_nonzero_pixels!(fine_grid, rate, coarse_rect)
                    @show sum(fine_grid.A .> 0)
                    before = sum(fine_band.A .> 0)
                    write_datagrid!(fine_grid, fine_band)
                    after = sum(fine_band.A .> 0)
                    @show (before, after)
                end
            else
                outside_cnt += 1
            end
            if pixel_limit >= max_pixels
                return nothing
            end
        end
    end
    @show outside_cnt
end

function whole_band(fine_band, coarse_band, max_pixels)
    fine_pg = pixelgrid(fine_band)
    coarse_pg = pixelgrid(coarse_band)
    coarse_crop_pg = crop_to(coarse_pg, xy_bounds(fine_pg))
    ccpg = coarse_crop_pg
    base = (ccpg.ul[1] - 1, ccpg.ul[2] - 1)
    coarse = ArchGDAL.read(coarse_band, ccpg.ul[1]:ccpg.lr[1], ccpg.ul[1]:ccpg.lr[2])
    fine = ArchGDAL.read(fine_band)
    fine_no_data = ArchGDAL.getnodatavalue(fine_band)

    coarse_count = zeros(Int32, size(coarse)...)
    for jc in size(coarse, 2)
        for ic in size(coarse, 1)
            single_pg = onepixelgrid(coarse_crop_pg, base + [ic, jc])
            rect = xy_bounds(single_pg)
            ulf, lrf = ij_cover_rect(single_pg.origin, single_pg.transform, rect)
            for jff in ulf[2]:lrf[2]
                for iff in ulf[1]:lrf[1]
                    center = geo_to_xy_center(fine_pg.origin, fine_pg.transform, iff, jff)
                    if pointinrect(center, rect) && fine[iff, jff] != fine_no_data
                        coarse_count[ic, jc] += 1
                    end
                end
            end
        end
    end

    for jff in 1:size(fine, 2)
        for iff in 1:size(fine, 1)
            if fine[iff, jff] != fine_no_data
                center = geo_to_xy_center(fine_pg.origin, fine_pg.transform, iff, jff)
                ls_pixel = pixel_containing(coarse_pg.origin, coarse_pg.transform, center)
                ls_shift = ls_pixel - base
                count = coarse_count[ls_shift...]
                pop = coarse[ls_shift...]
                rate = min(pop / count, 100)
                fine[iff, jff] = rate
            end
        end
    end
    ArchGDAL.write!(fine_band, fine)
end

"""
Given HRSL on 30m grid and LandScan (LS) on 900m grid,
transfer landscan density to HRSL settlement pixels that are nonzero.

1. Find bounds on points b/c that's a country's HRSL.
2. Find ls rect that covers those points.
3. Find block rect that covers the ls rect.
4. For each block in ls.
   1. Load block.
   2. For each pixel inside ls rect that covers points.
      1. Get points rect within pixel. (some pixels only partially covered)
      2. Load points blocks that cover that rect.
      3. Sum number of nonzero points within rect.
      4. Rescale by pixel value.
      5. Write modified points blocks.
"""
function assignweights(
        points_path, weight_path, reweighted_path;
        max_pixels = typemax(Int64)
        )
    # ds = dataset
    fine_copy_ds = ArchGDAL.read(expanduser(points_path)) do fine_ds
        reweighted_path = expanduser(reweighted_path)
        ArchGDAL.copy(fine_ds, filename = reweighted_path)
    end
    ArchGDAL.read(expanduser(weight_path)) do coarse_ds
        fine_band = getband(fine_copy_ds, 1)
        coarse_band = getband(coarse_ds, 1)
        whole_band(fine_band, coarse_band, max_pixels)
    end
end


text: Dict{String, Any}("tbl_multiline" => Dict{String, Any}("c" => Dict{String, Any}("value" => "and yet\nanother line", "type" => "string"), "b" => Dict{String, Any}("value" => "multiline\n", "type" => "string"), "a" => Dict{String, Any}("value" => "1", "type" => "integer"), "d" => Dict{String, Any}("value" => "4", "type" => "integer")))


text:
using Unitful

export MDFv2InMemoryPart, MDFv2Variables, MDFv2InMemoryPart, MDFv2Root, MDFv2Study,
       MDFv2Experiment, MDFv2Tracer, MDFv2Scanner, MDFv2Drivefield, MDFv2Receiver,
       MDFv2Acquisition, MDFv2Measurement, MDFv2Calibration, MDFv2Reconstruction, MDFv2InMemory

export defaultMDFv2Root, defaultMDFv2Study, defaultMDFv2Experiment,
       defaultMDFv2Tracer, defaultMDFv2Scanner, defaultMDFv2Drivefield,
       defaultMDFv2Receiver, defaultMDFv2Acquisition, defaultMDFv2Measurement,
       defaultMDFv2Calibration, defaultMDFv2Reconstruction, defaultMDFv2InMemory

export checkConsistency, inMemoryMDFToDict, inMemoryMDFFromDict

abstract type MDFv2InMemoryPart end

mutable struct MDFv2Variables
  "tracer materials/injections for multi-color MPI"
  A::Union{Int64, Nothing}
  "acquired frames (N = O + E), same as a spatial position for calibration"
  N::Union{Int64, Nothing}
  "acquired background frames (E = N − O)"
  E::Union{Int64, Nothing}
  "acquired foreground frames (O = N − E)"
  O::Union{Int64, Nothing}
  "coefficients stored after sparsity transformation (B \\le O)"
  B::Union{Int64, Nothing}
  "periods within one frame"
  J::Union{Int64, Nothing}
  "partitions of each patch position"
  Y::Union{Int64, Nothing}
  "receive channels"
  C::Union{Int64, Nothing}
  "drive-field channels"
  D::Union{Int64, Nothing}
  "frequencies describing the drive-field waveform"
  F::Union{Int64, Nothing}
  "points sampled at receiver during one drive-field period"
  V::Union{Int64, Nothing}
  "sampling points containing processed data (W = V if no frequency selection or bandwidth reduction has been applied)"
  W::Union{Int64, Nothing}
  "frequencies describing the processed data (K = V/2 + 1 if no frequency selection or bandwidth reduction has been applied)"
  K::Union{Int64, Nothing}
  "frames in the reconstructed MPI data set"
  Q::Union{Int64, Nothing}
  "voxels in the reconstructed MPI data set"
  P::Union{Int64, Nothing}
  "channels in the reconstructed MPI data set"
  S::Union{Int64, Nothing}

  function MDFv2Variables()
    return new(nothing, nothing, nothing, nothing, nothing, nothing, nothing,
               nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing)
  end
end

"Root group of an in-memory MDF"
mutable struct MDFv2Root <: MDFv2InMemoryPart
  "UTC creation time of MDF data set"
  time::Union{DateTime, Missing}
  "Universally Unique Identifier (RFC 4122) of MDF file"
  uuid::Union{UUID, Missing}
  "Version of the file format"
  version::Union{VersionNumber, Missing}

  function MDFv2Root(;
    time = missing,
    uuid = missing,
    version = missing)

    return new(
      time,
      uuid,
      version
    )
  end
end

defaultMDFv2Root() = MDFv2Root(time=Dates.now(), uuid=UUIDs.uuid4(), version=VersionNumber("2.1.0"))

"Study group of an in-memory MDF"
mutable struct MDFv2Study <: MDFv2InMemoryPart
  "Short description of the study"
  description::Union{String, Missing}
  "Name of the study"
  name::Union{String, Missing}
  "Number of the study"
  number::Union{Int64, Missing}
  "UTC creation time of study; optional"
  time::Union{DateTime, Nothing}
  "Universally Unique Identifier (RFC 4122) of study"
  uuid::Union{UUID, Missing}

  function MDFv2Study(;
    description = missing,
    name = missing,
    number = missing,
    time = nothing,
    uuid = missing)

    return new(
      description,
      name,
      number,
      time,
      uuid
    )
  end
end

defaultMDFv2Study() = MDFv2Study()

"Experiment group of an in-memory MDF"
mutable struct MDFv2Experiment <: MDFv2InMemoryPart
  "Short description of the experiment"
  description::Union{String, Missing}
  "Flag indicating if the data in this file is simulated rather than measured"
  isSimulation::Union{Bool, Missing}
  "Experiment name"
  name::Union{String, Missing}
  "Experiment number within study"
  number::Union{Int64, Missing}
  "Name of the subject that was imaged"
  subject::Union{String, Missing}
  "Universally Unique Identifier (RFC 4122) of experiment"
  uuid::Union{UUID, Missing}

  function MDFv2Experiment(;
    description = missing,
    isSimulation = missing,
    name = missing,
    number = missing,
    subject = missing,
    uuid = missing)

    return new(
      description,
      isSimulation,
      name,
      number,
      subject,
      uuid
    )
  end
end

defaultMDFv2Experiment() = MDFv2Experiment() # Should we create the UUID automatically?

"Tracer group of an in-memory MDF; optional"
mutable struct MDFv2Tracer <: MDFv2InMemoryPart
  "Batch of tracer"
  batch::Union{Vector{String}, Missing}
  "A mol(solute)/L no Molar concentration of solute per litre"
  concentration::Union{Vector{Float64}, Missing}
  "UTC time at which tracer injection started; optional"
  injectionTime::Union{Vector{DateTime}, Nothing}
  "Name of tracer used in experiment"
  name::Union{Vector{String}, Missing}
  "Solute, e.g. Fe"
  solute::Union{Vector{String}, Missing}
  "Name of tracer supplier"
  vendor::Union{Vector{String}, Missing}
  "Total volume of applied tracer"
  volume::Union{Vector{Float64}, Missing}

  function MDFv2Tracer(;
    batch = missing,
    concentration = missing,
    injectionTime = nothing,
    name = missing,
    solute = missing,
    vendor = missing,
    volume = missing)

    return new(
      batch,
      concentration,
      injectionTime,
      name,
      solute,
      vendor,
      volume
    )
  end
end

defaultMDFv2Tracer() = MDFv2Tracer()

"Scanner group of an in-memory MDF"
mutable struct MDFv2Scanner <: MDFv2InMemoryPart
  "Diameter of the bore; optional"
  boreSize::Union{Float64, Nothing}
  "Facility where the MPI scanner is installed"
  facility::Union{String, Missing}
  "Scanner manufacturer"
  manufacturer::Union{String, Missing}
  "Scanner name"
  name::Union{String, Missing}
  "User who operates the MPI scanner"
  operator::Union{String, Missing}
  "Scanner topology (e.g. FFP, FFL, MPS)"
  topology::Union{String, Missing}

  function MDFv2Scanner(;
    boreSize = nothing,
    facility = missing,
    manufacturer = missing,
    name = missing,
    operator = missing,
    topology = missing)

    return new(
      boreSize,
      facility,
      manufacturer,
      name,
      operator,
      topology
    )
  end
end

defaultMDFv2Scanner() = MDFv2Scanner()

"Drivefield subgroup of acquisition group of an in-memory MDF"
mutable struct MDFv2Drivefield <: MDFv2InMemoryPart
  "Base frequency to derive drive field frequencies"
  baseFrequency::Union{Float64, Missing}
  "Trajectory cycle is determined by lcm(divider)/baseFrequency. It will not change
  when averaging was applied. The duration for measuring the V data points (i.e. the
  drive-field period) is given by the product of period and numAverages"
  cycle::Union{Float64, Missing}
  "Divider of the baseFrequency to determine the drive field frequencies"
  divider::Union{Array{Int64, 2}, Missing}
  "Number of drive field channels, denoted by D"
  numChannels::Union{Int64, Missing}
  "Applied drive field phase"
  phase::Union{Array{Float64, 3}, Missing}
  "Applied drive field strength"
  strength::Union{Array{Float64, 3}, Missing}
  "Waveform type: sine, triangle or custom"
  waveform::Union{Array{String, 2}, Missing}

  function MDFv2Drivefield(;
    baseFrequency = missing,
    cycle = missing,
    divider = missing,
    numChannels = missing,
    phase = missing,
    strength = missing,
    waveform = missing)

    return new(
      baseFrequency,
      cycle,
      divider,
      numChannels,
      phase,
      strength,
      waveform
    )
  end
end

defaultMDFv2Drivefield() = MDFv2Drivefield()

"Receiver subgroup of acquisition group of an in-memory MDF"
mutable struct MDFv2Receiver <: MDFv2InMemoryPart
  "Bandwidth of the receiver unit"
  bandwidth::Union{Float64, Missing}
  "Dimension less scaling factor and offset (a_c, b_c) to convert raw data into a
  physical quantity with corresponding unit of measurement unit; optional"
  dataConversionFactor::Union{Array{Float64, 2}, Nothing}
  "Induction factor mapping the projection of the magnetic moment to the voltage in the receive coil; optional"
  inductionFactor::Union{Vector{Float64}, Nothing}
  "Number of receive channels C"
  numChannels::Union{Int64, Missing}
  "Number of sampling points during one period, denoted by V"
  numSamplingPoints::Union{Int64, Missing}
  "Transfer function of the receive channels in Fourier domain. unit is the field
  from the /measurement group; optional"
  transferFunction::Union{Array{ComplexF64, 2}, Nothing}
  "SI unit of the measured quantity, usually Voltage V"
  unit::Union{String, Missing}

  function MDFv2Receiver(;
    bandwidth = missing,
    dataConversionFactor = nothing,
    inductionFactor = nothing,
    numChannels = missing,
    numSamplingPoints = missing,
    transferFunction = nothing,
    unit = missing)

    return new(
      bandwidth,
      dataConversionFactor,
      inductionFactor,
      numChannels,
      numSamplingPoints,
      transferFunction,
      unit
    )
  end
end

defaultMDFv2Receiver() = MDFv2Receiver(unit = "V")

"Acquisition group of an in-memory MDF"
mutable struct MDFv2Acquisition <: MDFv2InMemoryPart
  "Gradient strength of the selection field in x, y, and z directions; optional"
  gradient::Union{Array{Float64, 4}, Nothing}
  "Number of block averages per drive-field period."
  numAverages::Union{Int64, Missing}
  "Number of available measurement frames"
  numFrames::Union{Int64, Missing}
  "Number of drive-field periods within a frame denoted by J"
  numPeriodsPerFrame::Union{Int64, Missing}
  "Offset field applied; optional"
  offsetField::Union{Array{Float64, 3}, Nothing}
  "UTC start time of MPI measurement"
  startTime::Union{DateTime, Missing}

  drivefield::Union{MDFv2Drivefield, Missing}
  receiver::Union{MDFv2Receiver, Missing}

  function MDFv2Acquisition(;
    gradient = nothing,
    numAverages = missing,
    numFrames = missing,
    numPeriodsPerFrame = missing,
    offsetField = nothing,
    startTime = missing,
    drivefield = missing,
    receiver = missing)

    return new(
      gradient,
      numAverages,
      numFrames,
      numPeriodsPerFrame,
      offsetField,
      startTime,
      drivefield,
      receiver
    )
  end
end

defaultMDFv2Acquisition() = MDFv2Acquisition(drivefield=defaultMDFv2Drivefield(), receiver=defaultMDFv2Receiver())

"Measurement group of an in-memory MDF"
mutable struct MDFv2Measurement <: MDFv2InMemoryPart
  "Measured data at a specific processing stage"
  data::Union{AbstractArray{<:Number, 4}, Missing}
  "Indices of original frame order; optional if !isFramePermutation"
  framePermutation::Union{Vector{Int64}, Nothing}
  "Indices of selected frequency components; optional if !isFrequencySelection"
  frequencySelection::Union{Vector{Int64}, Nothing}
  "Flag, if the background has been subtracted"
  isBackgroundCorrected::Union{Bool, Missing}
  "Mask indicating for each of the N frames if it is a background measurement (true) or not"
  isBackgroundFrame::Union{Vector{Bool}, Missing}
  "Flag, if the frame dimension N has been moved to the last dimension"
  isFastFrameAxis::Union{Bool, Missing}
  "Flag, if the data is stored in frequency space"
  isFourierTransformed::Union{Bool, Missing}
  "Flag, if the order of frames has been changed, see framePermutation"
  isFramePermutation::Union{Bool, Missing}
  "Flag, if only a subset of frequencies has been selected and stored, see frequencySelection"
  isFrequencySelection::Union{Bool, Missing}
  "Flag, if the foreground frames are compressed along the frame dimension"
  isSparsityTransformed::Union{Bool, Missing}
  "Flag, if spectral leakage correction has been applied"
  isSpectralLeakageCorrected::Union{Bool, Missing}
  "Flag, if the data has been corrected by the transferFunction"
  isTransferFunctionCorrected::Union{Bool, Missing}
  "Name of the applied sparsity transformation; optional if !isSparsityTransformed"
  sparsityTransformation::Union{String, Nothing}
  "Subsampling indices \\beta{j,c,k,b}; optional if !isSparsityTransformed"
  subsamplingIndices::Union{Array{Integer, 4}, Nothing}

  function MDFv2Measurement(;
    data = missing,
    framePermutation = nothing,
    frequencySelection = nothing,
    isBackgroundCorrected = missing,
    isBackgroundFrame = missing,
    isFastFrameAxis = missing,
    isFourierTransformed = missing,
    isFramePermutation = missing,
    isFrequencySelection = missing,
    isSparsityTransformed = missing,
    isSpectralLeakageCorrected = missing,
    isTransferFunctionCorrected = missing,
    sparsityTransformation = nothing,
    subsamplingIndices = nothing)

    return new(
      data,
      framePermutation,
      frequencySelection,
      isBackgroundCorrected,
      isBackgroundFrame,
      isFastFrameAxis,
      isFourierTransformed,
      isFramePermutation,
      isFrequencySelection,
      isSparsityTransformed,
      isSpectralLeakageCorrected,
      isTransferFunctionCorrected,
      sparsityTransformation,
      subsamplingIndices
    )
  end
end

defaultMDFv2Measurement() = MDFv2Measurement()

"Calibration group of an in-memory MDF"
mutable struct MDFv2Calibration <: MDFv2InMemoryPart
  "Size of the delta sample used for calibration scan; optional"
  deltaSampleSize::Union{Vector{Float64}, Nothing}
  "Field of view of the system matrix; optional"
  fieldOfView::Union{Vector{Float64}, Nothing}
  "Center of the system matrix (relative to origin/center); optional"
  fieldOfViewCenter::Union{Vector{Float64}, Nothing}
  "Method used to obtain calibration data. Can for instance be robot, hybrid or simulation"
  method::Union{String, Missing}
  "Applied offset field strength to emulate a spatial position (x, y, z); optional"
  offsetFields::Union{Array{Float64, 2}, Nothing}
  "Ordering of the dimensions, default is xyz; optional"
  order::Union{String, Nothing}
  "Position of each of the grid points, stored as (x, y, z) triples; optional"
  positions::Union{Array{Float64, 2}, Nothing}
  "Number of voxels in each dimension, inner product is O; optional"
  size::Union{Vector{Int64}, Nothing}
  "Signal-to-noise estimate for recorded frequency components; optional"
  snr::Union{Array{Float64, 3}, Nothing}

  function MDFv2Calibration(;
    deltaSampleSize = nothing,
    fieldOfView = nothing,
    fieldOfViewCenter = nothing,
    method = missing,
    offsetFields = nothing,
    order = nothing,
    positions = nothing,
    size = nothing,
    snr = nothing)

    return new(
      deltaSampleSize,
      fieldOfView,
      fieldOfViewCenter,
      method,
      offsetFields,
      order,
      positions,
      size,
      snr
    )
  end
end

defaultMDFv2Calibration() = MDFv2Calibration(order="xyz")

"Reconstruction group of an in-memory MDF"
mutable struct MDFv2Reconstruction <: MDFv2InMemoryPart
  "Reconstructed data"
  data::Union{Array{Number, 3}, Missing}
  "Field of view of reconstructed data; optional"
  fieldOfView::Union{Vector{Float64}, Nothing}
  "Center of the reconstructed data (relative to scanner origin/center); optional"
  fieldOfViewCenter::Union{Vector{Float64}, Nothing}
  "Mask indicating for each of the P voxels if it is part of the overscan region (true) or not; optional"
  isOverscanRegion::Union{Vector{Bool}, Nothing}
  "Ordering of the dimensions, default is xyz; optional"
  order::Union{String, Nothing}
  "Position of each of the grid points, stored as (x, y, z) tripels; optional"
  positions::Union{Array{Float64, 2}, Nothing}
  "Number of voxels in each dimension, inner product is P; optional"
  size::Union{Vector{Int64}, Nothing}

  function MDFv2Reconstruction(;
    data = missing,
    fieldOfView = nothing,
    fieldOfViewCenter = nothing,
    isOverscanRegion = nothing,
    order = nothing,
    positions = nothing,
    size = nothing)

    return new(
      data,
      fieldOfView,
      fieldOfViewCenter,
      isOverscanRegion,
      order,
      positions,
      size
    )
  end
end

defaultMDFv2Reconstruction() = MDFv2Reconstruction(order="xyz")

mutable struct MDFv2InMemory <: MPIFile # TODO: Not sure, if MPIFile is a good fit
  root::Union{MDFv2Root, Missing}
  study::Union{MDFv2Study, Missing}
  experiment::Union{MDFv2Experiment, Missing}
  tracer::Union{MDFv2Tracer, Nothing}
  scanner::Union{MDFv2Scanner, Missing}
  acquisition::Union{MDFv2Acquisition, Missing}
  measurement::Union{MDFv2Measurement, Nothing}
  calibration::Union{MDFv2Calibration, Nothing}
  reconstruction::Union{MDFv2Reconstruction, Nothing}
  custom::Dict{String, Any}
  variables::MDFv2Variables

  function MDFv2InMemory(;
    root = missing,
    study = missing,
    experiment = missing,
    tracer = nothing,
    scanner = missing,
    acquisition = missing,
    measurement = nothing,
    calibration = nothing,
    reconstruction = nothing,
    custom = Dict{String, Any}(),
    variables = MDFv2Variables())

    return new(
      root,
      study,
      experiment,
      tracer,
      scanner,
      acquisition,
      measurement,
      calibration,
      reconstruction,
      custom,
      variables
    )
  end
end

MDFv2InMemory(mdfFile::MDFFileV2) = inMemoryMDFFromMDFFileV2(mdfFile)
MDFv2InMemory(dict::Dict) = inMemoryMDFFromDict(dict)

function defaultMDFv2InMemory()
  return MDFv2InMemory(
    root = defaultMDFv2Root(),
    study = defaultMDFv2Study(),
    experiment = defaultMDFv2Experiment(),
    scanner = defaultMDFv2Scanner(),
    acquisition = defaultMDFv2Acquisition()
  )
end

function check(part::MDFv2Root, variables::MDFv2Variables)
  @assert part.version > VersionNumber("2") "The version should be at least version 2."
end

function check(part::MDFv2Study, variables::MDFv2Variables)
  # nothing to check yet
end

function check(part::MDFv2Experiment, variables::MDFv2Variables)
  # nothing to check yet
end

function check(part::MDFv2Tracer, variables::MDFv2Variables)
  # Init comparison
  if isnothing(variables.A)
    variables.A = length(part.vendor)
  end

  # Check if all sizes match
  for fieldname in fieldnames(MDFv2Tracer)
    field = getproperty(part, fieldname)
    @assert length(field) == variables.A "Inconsistent dimensions in `$fieldname` in `tracer`."
  end
end

function check(part::MDFv2Scanner, variables::MDFv2Variables)
  # nothing to check yet
end

function check(part::MDFv2Acquisition, variables::MDFv2Variables)
  # Pick variables first
  if isnothing(variables.N)
    variables.N = part.numFrames
  end
  if isnothing(variables.J)
    variables.J = part.numPeriodsPerFrame
  end

  # Check dimensions of `gradient` field
  if !isnothing(part.gradient)
    if isnothing(variables.J)
      variables.J = size(part.gradient, 1)
    else
      @assert variables.J == size(part.gradient, 1) "Inconsistent dimension J in `gradient` in `acquisition`."
    end

    if isnothing(variables.Y)
      variables.Y = size(part.gradient, 2)
    else
      @assert variables.Y == size(part.gradient, 1) "Inconsistent dimension Y in `gradient` in `acquisition`."
    end

    @assert size(part.gradient, 3) == 3 "Inconsistent third dimension in `gradient` in `acquisition`."
    @assert size(part.gradient, 4) == 3 "Inconsistent fourth dimension in `gradient` in `acquisition`."
  end

  # Check dimensions of `offsetField` field
  if !isnothing(part.offsetField)
    if isnothing(variables.J)
      variables.J = size(part.offsetField, 1)
    else
      @assert variables.J == size(part.offsetField, 1) "Inconsistent dimension J in `offsetField` in `acquisition`."
    end

    if isnothing(variables.Y)
      variables.Y = size(part.offsetField, 2)
    else
      @assert variables.Y == size(part.offsetField, 2) "Inconsistent dimension Y in `offsetField` in `acquisition`."
    end

    @assert size(part.offsetField, 3) == 3 "Inconsistent third dimension in `offsetField` in `acquisition`."
  end
end

function check(part::MDFv2Drivefield, variables::MDFv2Variables)
  # Pick variables first
  if isnothing(variables.D)
    variables.D = part.numChannels
  end
  if isnothing(variables.F)
    variables.F = size(part.divider, 2)
  end
  if isnothing(variables.J)
    variables.J = size(part.phase, 1)
  end

  # Then check dimensions for multidimensional fields
  @assert variables.D == size(part.divider, 1) "Inconsistent dimension D in `divider` in `drivefield`."
  @assert variables.F == size(part.divider, 2) "Inconsistent dimension F in `divider` in `drivefield`."

  @assert variables.J == size(part.phase, 1) "Inconsistent dimension J in `phase` in `drivefield`."
  @assert variables.D == size(part.phase, 2) "Inconsistent dimension D in `phase` in `drivefield`."
  @assert variables.F == size(part.phase, 3) "Inconsistent dimension F in `phase` in `drivefield`."

  @assert variables.J == size(part.strength, 1) "Inconsistent dimension J in `strength` in `drivefield`."
  @assert variables.D == size(part.strength, 2) "Inconsistent dimension D in `strength` in `drivefield`."
  @assert variables.F == size(part.strength, 3) "Inconsistent dimension F in `strength` in `drivefield`."

  @assert variables.D == size(part.waveform, 1) "Inconsistent dimension D in `waveform` in `drivefield`."
  @assert variables.F == size(part.waveform, 2) "Inconsistent dimension F in `waveform` in `drivefield`."
end

function check(part::MDFv2Receiver, variables::MDFv2Variables)
  # C is defined by numChannels
  if isnothing(variables.C)
    variables.C = part.numChannels
  end
  # V is defined by numSamplingPoints
  if isnothing(variables.V)
    variables.V = part.numSamplingPoints
  end

  if !isnothing(part.dataConversionFactor)
    @assert variables.C == size(part.dataConversionFactor, 1) "Inconsistent dimension C in `dataConversionFactor` in `receiver`."
    @assert size(part.dataConversionFactor, 2) == 2 "Inconsistent second dimension in `dataConversionFactor` in `receiver`."
  end

  if !isnothing(part.inductionFactor)
    @assert variables.C == size(part.inductionFactor, 1) "Inconsistent dimension C in `inductionFactor` in `receiver`."
  end

  if !isnothing(part.transferFunction)
    @assert variables.C == size(part.transferFunction, 1) "Inconsistent dimension C in `transferFunction` in `receiver`."

    if isnothing(variables.K)
      variables.K = size(part.transferFunction, 2)
    else
      @assert variables.K == size(part.transferFunction, 2) "Inconsistent dimension K in `transferFunction` in `receiver`."
    end
  end
end

function check(part::MDFv2Measurement, variables::MDFv2Variables)
  # TODO: check data consistency; I don't know how to do that since we can't know all variables
  # N × J × C × K or
  # J × C × K × N or
  # N × J × C × W or
  # J × C × W × N or
  # J×C×K×(B+E)

  # N is defined by the length of isBackgroundFrame here; should have been retrieved earlier
  if isnothing(variables.N)
    variables.N = length(part.isBackgroundFrame)
  end
  # E is defined by the number of `true` values in isBackgroundFrame here; should have been retrieved earlier
  if isnothing(variables.E)
    variables.E = length([x for x in part.isBackgroundFrame if x == true])
  end
  # E is defined by the number of `false` values in isBackgroundFrame here; should have been retrieved earlier
  if isnothing(variables.O)
    variables.O = length([x for x in part.isBackgroundFrame if x == false])
  end
  if isnothing(variables.W)
    if isnothing(variables.V)
      @warn "Can't determine variable W since variable V is not defined. Should happen in receiver."
    else
      variables.W = variables.V
    end
  end

  if isnothing(variables.J)
    if !part.isSparsityTransformed
      @warn "Can't determine variable J for measurement from `subsamplingIndices` since the dataset is not sparsity transformed."
    else
      variables.J = size(part.subsamplingIndices, 1)
    end
  end
  if isnothing(variables.C)
    if !part.isSparsityTransformed
      @warn "Can't determine variable C for measurement from `subsamplingIndices` since the dataset is not sparsity transformed. Should happen in receiver."
    else
      variables.C = size(part.subsamplingIndices, 2)
    end
  end
  if isnothing(variables.K)
    if !part.isSparsityTransformed
      @warn "Can't determine variable K for measurement from `subsamplingIndices` since the dataset is not sparsity transformed. Should happen in receiver."
    else
      variables.K = size(part.subsamplingIndices, 3)
    end
  end
  if isnothing(variables.B)
    if !part.isSparsityTransformed
      @warn "Can't determine variable B for measurement from `subsamplingIndices` since the dataset is not sparsity transformed. Should happen in receiver."
    else
      variables.B = size(part.subsamplingIndices, 4)
    end
  end

  if part.isFramePermutation
    @assert variables.N == length(part.framePermutation) "Inconsistent dimension N in `framePermutation` in `measurement`."
  end

  if part.isFrequencySelection
    if isnothing(variables.K)
      variables.K = length(part.frequencySelection)
    else
      @assert variables.K == length(part.frequencySelection) "Inconsistent dimension K in `frequencySelection` in `measurement`."
    end
  end

  @assert variables.N == length(part.isBackgroundFrame) "Inconsistent dimension N in `isBackgroundFrame` in `measurement`."

  if part.isSparsityTransformed
    @assert !isnothing(part.sparsityTransformation) "Field `sparsityTransformation` must be set when `isSparsityTransformed` is set in in `measurement`."
  end

  if part.isSparsityTransformed
    # J, C, K and B should be defined by now
    @assert variables.J == size(part.subsamplingIndices, 1) "Inconsistent dimension J in `subsamplingIndices` in `measurement`."
    @assert variables.C == size(part.subsamplingIndices, 2) "Inconsistent dimension C in `subsamplingIndices` in `measurement`."
    @assert variables.K == size(part.subsamplingIndices, 3) "Inconsistent dimension K in `subsamplingIndices` in `measurement`."
    @assert variables.B == size(part.subsamplingIndices, 4) "Inconsistent dimension B in `subsamplingIndices` in `measurement`."
  end
end

function check(part::MDFv2Calibration, variables::MDFv2Variables)
  if !isnothing(part.deltaSampleSize)
    @assert length(part.deltaSampleSize) == 3 "Inconsistent length in `deltaSampleSize` in `calibration`."
  end

  if !isnothing(part.fieldOfView)
    @assert length(part.fieldOfView) == 3 "Inconsistent length in `fieldOfView` in `calibration`."
  end

  if !isnothing(part.fieldOfViewCenter)
    @assert length(part.fieldOfViewCenter) == 3 "Inconsistent length in `fieldOfViewCenter` in `calibration`."
  end

  if !isnothing(part.offsetFields)
    if isnothing(variables.O)
      variables.O = size(part.offsetFields, 1)
    else
      @assert variables.O == size(part.offsetFields, 1) "Inconsistent dimension O in `offsetFields` in `calibration`."
    end
    @assert size(part.offsetFields, 2) == 3 "Inconsistent second dimension in `offsetFields` in `calibration`."
  end

  if !isnothing(part.order)
    @assert part.order in ["xyz", "xzy", "yxz", "yzx", "zyx", "zxy"] "Wrong `order` of `$(part.order)` in `calibration`."
  end

  if !isnothing(part.positions)
    # O must be defined by now
    @assert variables.O == size(part.positions, 1) "Inconsistent dimension O in `positions` in `calibration`."
    @assert size(part.positions, 2) == 3 "Inconsistent second dimension in `positions` in `calibration`."
  end

  if !isnothing(part.size)
    @assert length(part.size) == 3 "Inconsistent length in `size` in `calibration`."
    # O must be defined by now
    @assert variables.O == prod(part.size) "The product of `size` with `$(part.size)` must equal O."
  end

  if !isnothing(part.snr)
    if isnothing(variables.J)
      variables.J = size(part.snr, 1)
    else
      @assert variables.J == size(part.snr, 1) "Inconsistent dimension J in `snr` in `calibration`."
    end

    if isnothing(variables.C)
      variables.C = size(part.snr, 2)
    else
      @assert variables.C == size(part.snr, 2) "Inconsistent dimension C in `snr` in `calibration`."
    end

    if isnothing(variables.K)
      variables.K = size(part.snr, 3)
    else
      @assert variables.K == size(part.snr, 3) "Inconsistent dimension K in `snr` in `calibration`."
    end
  end
end

function check(part::MDFv2Reconstruction, variables::MDFv2Variables)
  # Pick variables first
  if isnothing(variables.Q)
    variables.Q = size(part.data, 1)
  end
  if isnothing(variables.P)
    variables.P = size(part.data, 2)
  end
  if isnothing(variables.S)
    variables.S = size(part.data, 1)
  end

  if !isnothing(part.fieldOfView)
    @assert length(part.fieldOfView) == 3 "Inconsistent length in `fieldOfView` in `reconstruction`."
  end

  if !isnothing(part.fieldOfViewCenter)
    @assert length(part.fieldOfViewCenter) == 3 "Inconsistent length in `fieldOfViewCenter` in `reconstruction`."
  end

  if !isnothing(part.isOverscanRegion)
    @assert variables.P == length(part.isOverscanRegion) "Inconsistent length in `isOverscanRegion` in `reconstruction`."
  end

  if !isnothing(part.order)
    @assert part.order in ["xyz", "xzy", "yxz", "yzx", "zyx", "zxy"] "Wrong `order` of `$(part.order)` in `reconstruction`."
  end

  if !isnothing(part.positions)
    @assert variables.P == size(part.positions, 1) "Inconsistent dimension P in `positions` in `reconstruction`."
    @assert size(part.positions, 2) == 3 "Inconsistent second dimension in `positions` in `reconstruction`."
  end

  if !isnothing(part.size)
    @assert length(part.size) == 3 "Inconsistent length in `size` in `reconstruction`."
    # P must be defined by now
    @assert variables.P == prod(part.size) "The product of `size` with `$(part.size)` must equal P."
  end
end

"Check for missing fields in MDF parts"
function checkMissing(part::T) where T <: MDFv2InMemoryPart
  for (fieldname, fieldtype) in zip(fieldnames(T), fieldtypes(T))
    if !(fieldtype == MDFv2Variables)
      field = getproperty(part, fieldname)
      @assert !ismissing(field) "The field `$fieldname` is missing in the given in-memory MDF."
    end
  end
end

"Check, whether all non-optional fields have been
set and if the dimensions of the fields match"
function checkConsistency(mdf::MDFv2InMemory)
  for (fieldname, fieldtype) in zip(fieldnames(MDFv2InMemory), fieldtypes(MDFv2InMemory))
    if !(fieldtype == MDFv2Variables || fieldname == :custom)
      # At the moment, this should be a Union
      fieldtype = (fieldtype.b <: MDFv2InMemoryPart) ? fieldtype.b : fieldtype.a

      @debug "Checking consistency of `$fieldname`."

      field = getproperty(mdf, fieldname)
      @assert !ismissing(field) "The field `$fieldname` is missing in the given in-memory MDF."

      if !isnothing(field) # Only optional fields can be `nothing`
        checkMissing(field)
        check(field, mdf.variables)
      end

      # Check subgroups of acquisition
      if fieldtype == MDFv2Acquisition
        checkMissing(field.drivefield)
        check(field.drivefield, mdf.variables)

        checkMissing(field.receiver)
        check(field.receiver, mdf.variables)
      end
    end
  end

  return true # If no assertions failed, return true for checking in tests
end

### Create getters and setters

customSymbols = Dict{Symbol, String}(
  :dfCustomWaveform => "/acquisition/drivefield/customWaveform",
  :calibIsMeanderingGrid => "/calibration/isMeanderingGrid",
  :calibTemperatures => "/calibration/_temperatures",
  :rxTransferFunctionFileName => "/acquisition/receiver/transferFunctionFileName",
  :recoParameters => "/reconstruction/_parameters",
  :auxiliaryData => "/custom/auxiliaryData",
)

specificationSymbols = Dict{Symbol, String}()

aliases = Dict{Symbol, Symbol}(
  :measIsTransferFunctionCorrected => :measIsTFCorrected,
  :measIsBackgroundCorrected => :measIsBGCorrected,
  :measIsBackgroundFrame => :measIsBGFrame,
  :calibSnr => :calibSNR,
  :calibFieldOfView => :calibFov,
  :calibFieldOfViewCenter => :calibFovCenter,
  :recoFieldOfView => :recoFov,
  :recoFieldOfViewCenter => :recoFovCenter
)

prefixes = Dict{String, String}(
  "MDFv2Root" => "",
  "MDFv2Study" => "study",
  "MDFv2Experiment" => "experiment",
  "MDFv2Tracer" => "tracer",
  "MDFv2Scanner" => "scanner",
  "MDFv2Acquisition" => "acq",
  "MDFv2Drivefield" => "df",
  "MDFv2Receiver" => "rx",
  "MDFv2Measurement" => "meas",
  "MDFv2Calibration" => "calib",
  "MDFv2Reconstruction" => "reco"
)
for (fieldname, fieldtype) in zip(fieldnames(MDFv2InMemory), fieldtypes(MDFv2InMemory))
  fieldnameStr = string(fieldname)
  if !(fieldnameStr == "variables" || fieldnameStr == "custom")
    # At the moment, this should be a Union
    missingOrNothing = (fieldtype.b <: MDFv2InMemoryPart) ? fieldtype.a : fieldtype.b
    fieldtype = (fieldtype.b <: MDFv2InMemoryPart) ? fieldtype.b : fieldtype.a

    # Create getter and for the whole group
    @eval begin
      # TODO: Add docstring from struct; I did not yet find a way to retrieve it
      function $(fieldname)(mdf::MDFv2InMemory)::Union{$fieldtype, $missingOrNothing}
        return mdf.$fieldname
      end

      # TODO: Add docstring from struct; I did not yet find a way to retrieve it
      function $(fieldname)(mdf::MDFv2InMemory, value::Union{$fieldtype, $missingOrNothing})
        mdf.$fieldname = value
      end
    end

    for (partFieldname, partFieldtype) in zip(fieldnames(fieldtype), fieldtypes(fieldtype))
      partFieldnameStr = string(partFieldname)

      # The acquisition group has subgroups, so we need to go deeper there
      if !(partFieldnameStr == "drivefield" || partFieldnameStr == "receiver")
        if fieldtype != MDFv2Root
          capitalizedPartFieldname = uppercase(partFieldnameStr[1])*partFieldnameStr[2:end]
        else
          capitalizedPartFieldname = partFieldnameStr
        end
        functionSymbol = Symbol(prefixes[replace(string(fieldtype), "MPIFiles." => "")]*capitalizedPartFieldname)

        # Save symbols for later use in conversion
        if fieldnameStr != "root"
          specificationSymbols[functionSymbol] = "/"*fieldnameStr*"/"*partFieldnameStr
        else
          specificationSymbols[functionSymbol] = "/"*partFieldnameStr
        end

        @eval begin
          # TODO: Add docstring from struct; I did not yet find a way to retrieve it
          function $(functionSymbol)(mdf::MDFv2InMemory)::Union{$partFieldtype, $missingOrNothing}
            if !(isnothing(mdf.$fieldname) || ismissing(mdf.$fieldname))
              return mdf.$fieldname.$partFieldname
            else
              return mdf.$fieldname
            end
          end

          # TODO: Add docstring from struct; I did not yet find a way to retrieve it
          function $(functionSymbol)(mdf::MDFv2InMemory, value::Union{$partFieldtype, $missingOrNothing})
            # Automatically create fields if they do not exist
            if isnothing(mdf.$fieldname) || ismissing(mdf.$fieldname)
              @debug "Creating field $($fieldnameStr)"
              mdf.$fieldname = $fieldtype()
            end

            mdf.$fieldname.$partFieldname = value
          end
        end

        # If needed, create aliases
        if haskey(aliases, functionSymbol)
          alias = aliases[functionSymbol]
          @eval begin
            $(alias)(mdf::MDFv2InMemory)::Union{$partFieldtype, $missingOrNothing} = $(functionSymbol)(mdf)
            $(alias)(mdf::MDFv2InMemory, value::$partFieldtype) = $(functionSymbol)(mdf, value)
            $(functionSymbol)(f::MDFFileV2)::Union{$partFieldtype, $missingOrNothing} = $(alias)(f) # Should this be here or in MDF.jl?
          end
        end
      else
        # At the moment, this should be a Union
        missingOrNothing = (partFieldtype.b <: MDFv2InMemoryPart) ? partFieldtype.a : partFieldtype.b
        partFieldtype = (partFieldtype.b <: MDFv2InMemoryPart) ? partFieldtype.b : partFieldtype.a

        # Create getter and for the whole group
        @eval begin
          # TODO: Add docstring from struct; I did not yet find a way to retrieve it
          function $(partFieldname)(mdf::MDFv2InMemory)::Union{$partFieldtype, $missingOrNothing}
            if !(isnothing(mdf.$fieldname) || ismissing(mdf.$fieldname))
              return mdf.$fieldname.$partFieldname
            else
              return mdf.$fieldname
            end
          end

          # TODO: Add docstring from struct; I did not yet find a way to retrieve it
          function $(partFieldname)(mdf::MDFv2InMemory, value::Union{$partFieldtype, $missingOrNothing})
            # Automatically create fields if they do not exist
            if isnothing(mdf.$fieldname) || ismissing(mdf.$fieldname)
              @debug "Creating field $($fieldnameStr)"
              mdf.$fieldname = $fieldtype()
            end

            mdf.$fieldname.$partFieldname = value
          end
        end

        for (subPartFieldname, subPartFieldtype) in zip(fieldnames(partFieldtype), fieldtypes(partFieldtype))
          subPartFieldnameStr = string(subPartFieldname)
          capitalizedSubPartFieldname = uppercase(subPartFieldnameStr[1])*subPartFieldnameStr[2:end]
          functionSymbol = Symbol(prefixes[replace(string(partFieldtype), "MPIFiles." => "")]*capitalizedSubPartFieldname)

          # Save symbols for later use in conversion
          specificationSymbols[functionSymbol] = "/"*fieldnameStr*"/"*partFieldnameStr*"/"*subPartFieldnameStr

          # Create getter and setter for the respective field  within the naming scheme
          @eval begin
            # TODO: Add docstring from struct; I did not yet find a way to retrieve it
            function $(functionSymbol)(mdf::MDFv2InMemory)::Union{$subPartFieldtype, $missingOrNothing}
              if !(isnothing(mdf.$fieldname) || ismissing(mdf.$fieldname))
                if !(isnothing(mdf.$fieldname.$partFieldname) || ismissing(mdf.$fieldname.$partFieldname))
                  return mdf.$fieldname.$partFieldname.$subPartFieldname
                else
                  return mdf.$fieldname.$partFieldname
                end
              else
                return mdf.$fieldname
              end
            end

            # TODO: Add docstring from struct; I did not yet find a way to retrieve it
            function $(functionSymbol)(mdf::MDFv2InMemory, value::Union{$subPartFieldtype, $missingOrNothing})
              # Automatically create fields if they do not exist
              if isnothing(mdf.$fieldname) || ismissing(mdf.$fieldname)
                @debug "Creating field $($fieldnameStr)"
                mdf.$fieldname = $fieldtype()
              end

              if isnothing(mdf.$fieldname.$partFieldname) || ismissing(mdf.$fieldname.$partFieldname)
                @debug "Creating field $($fieldnameStr).$($partFieldnameStr)"
                mdf.$fieldname.$partFieldname = $partFieldtype()
              end

              mdf.$fieldname.$partFieldname.$subPartFieldname = value
            end
          end

          # If needed, create aliases
          if haskey(aliases, functionSymbol)
            alias = aliases[functionSymbol]
            @eval begin
              $(alias)(mdf::MDFv2InMemory)::Union{$partFieldtype, $missingOrNothing} = $(functionSymbol)(mdf)
              $(alias)(mdf::MDFv2InMemory, value::Union{$partFieldtype, $missingOrNothing}) = $(functionSymbol)(mdf, value)
              $(functionSymbol)(f::MDFFileV2)::Union{$partFieldtype, $missingOrNothing} = $(alias)(f) # Should this be here or in MDF.jl?
            end
          end
        end
      end
    end
  end
end

# The MDF specification uses the plural
calibOffsetFields(f::MDFFile) = calibOffsetField(f)

# And some utility functions
measIsCalibProcessed(mdf::MDFv2InMemory)::Union{Bool, Missing} = measIsFramePermutation(mdf) &&
                                                                 measIsFourierTransformed(mdf) &&
                                                                 measIsFastFrameAxis(mdf)

experimentHasReconstruction(mdf::MDFv2InMemory)::Bool = !isnothing(mdf.reconstruction)
experimentHasMeasurement(mdf::MDFv2InMemory)::Bool = !isnothing(mdf.measurement)
rxHasTransferFunction(mdf::MDFv2InMemory) = !isnothing(mdf.acquisition.receiver.transferFunction)

# Last part is a little workaround for the testcases, since the file somehow contains a `deltaSampleSize` field
experimentIsCalibration(mdf::MDFv2InMemory)::Bool = !isnothing(mdf.calibration) && !ismissing(mdf.calibration.method)

# Creation and conversion

"Create an in-memory MDF from a dict matching the respective function names."
function inMemoryMDFFromDict(dict::Dict{Symbol, Any})::MDFv2InMemory
  mdf = MDFv2InMemory()

  for (functionSymbol, value) in dict
    # Conversion uses aliases => replace these symbols with their alias
    if haskey(aliases, functionSymbol)
      functionSymbol = aliases[functionSymbol]
    end

    try
      f = getfield(MPIFiles, functionSymbol)
      f(mdf, value)
    catch e
      if e isa UndefVarError
        @warn "The function `$(string(functionSymbol))` corresponding to a key in the given dict could not be found."
      else
        rethrow()
      end
    end
  end

  return mdf
end
inMemoryMDFFromDict(dict::Dict{String, Any})::MDFv2InMemory = inMemoryMDFFromDict(Dict(map(x -> Symbol(x.first)=>x.second, collect(dict))))

"Create a dict from an in-memory MDF by matching the respective function names."
function inMemoryMDFToDict(mdf::MDFv2InMemory)::Dict{String, Any}
  resultDict = Dict{String, Any}()

  # Add standard-compliant and non-standard fields
  for functionSymbol in vcat(collect(keys(specificationSymbols)), collect(keys(customSymbols)))
    f = getfield(MPIFiles, functionSymbol)
    result = f(mdf)
    if !(isnothing(result) || ismissing(result))
      # Conversion uses aliases => replace these symbols with their alias
      if haskey(aliases, functionSymbol)
        functionSymbol = aliases[functionSymbol]
      end

      resultDict[string(functionSymbol)] = result
    end
  end

  # Add measurements data
  resultDict["measData"] = measDataRaw(mdf)

  return resultDict
end

"Alias to inMemoryMDFToDict()"
toDict(mdf::MDFv2InMemory)::Dict{String, Any} = inMemoryMDFToDict(mdf)

"Create an in-memory MDF from an MDFFile by calling the corresponding functions."
function inMemoryMDFFromMDFFileV2(mdfFile::MDFFileV2)::MDFv2InMemory
  mdf = MDFv2InMemory()

  # Add standard-compliant fields
  for functionSymbol in keys(specificationSymbols)
    f = getfield(MPIFiles, functionSymbol)
    result = f(mdfFile)
    if !(isnothing(result) || ismissing(result))
      f(mdf, result) # Call the setter of an MDFv2InMemory with a getter from an MDFFileV2
    end
  end

  # Add non-standard fields
  for functionSymbol in keys(customSymbols)
    f = getfield(MPIFiles, functionSymbol)
    result = f(mdfFile)
    if !(isnothing(result) || ismissing(result))
      mdf.custom[string(functionSymbol)] = result
    end
  end

  # Add measurements data
  measDataRaw(mdf, mdfFile.mmap_measData)

  return mdf
end

"Create an MDFFile from an in-memory MDF; alias to `saveasMDF`."
function inMemoryMDFToMDFFileV2(filename::String, mdf::MDFv2InMemory)
  saveasMDF(filename, mdf)
end

function saveasMDF(filename::String, mdf::MDFv2InMemory; failOnInconsistent::Bool=false)
  # file has to be removed if exists. Otherwise h5create fails.
  isfile(filename) && rm(filename)
  h5open(filename, "w") do file
    saveasMDF(file, mdf)
  end
end

"Create an MDFFile from an in-memory MDF."
function saveasMDF(file::HDF5.File, mdf::MDFv2InMemory; failOnInconsistent::Bool=false)
  if failOnInconsistent
    checkConsistency(mdf)
  else
    try
      checkConsistency(mdf)
    catch e
      if e isa AssertionError
        @warn "There is an inconsistency in the given in-memory MDF. The message is: `$(e.msg)`."
      else
        rethrow()
      end
    end
  end

  for (functionSymbol, key) in merge(customSymbols, specificationSymbols)
    f = getfield(MPIFiles, functionSymbol)
    result = f(mdf)
    if !(isnothing(result) || ismissing(result))
      # Convert datatypes to something compatible with HDF5
      if eltype(result) == Bool
        result = Int8.(result)
      end
      if result isa VersionNumber || result isa DateTime || result isa UUID
        result = string(result)
      end
      if eltype(result) == DateTime
        result = string.(result)
      end

      file[key] = result
    end
  end
end

# This is non-standard (add new non-standard functions to `customSymbols` in order to have the custom fields set!)
dfCustomWaveform(mdf::MDFv2InMemory)::Union{String, Nothing} = @keyoptional mdf.custom["dfCustomWaveform"] # TODO: Should this be a 2D array?
dfCustomWaveform(mdf::MDFv2InMemory, customWaveform::String) = mdf.custom["dfCustomWaveform"] = customWaveform

calibIsMeanderingGrid(mdf::MDFv2InMemory)::Union{Bool, Nothing} = @keyoptional Bool(mdf.custom["calibIsMeanderingGrid"])
calibIsMeanderingGrid(mdf::MDFv2InMemory, meandering::Bool) = mdf.custom["calibIsMeanderingGrid"] = Int(meandering)

rxTransferFunctionFileName(mdf::MDFv2InMemory)::Union{String, Nothing} = @keyoptional mdf.custom["rxTransferFunctionFileName"]
rxTransferFunctionFileName(mdf::MDFv2InMemory, filename::String) = mdf.custom["rxTransferFunctionFileName"] = filename

recoParameters(mdf::MDFv2InMemory) = @keyoptional mdf.custom["recoParameters"]
recoParameters(mdf::MDFv2InMemory, parameters) = mdf.custom["recoParameters"] = parameters

calibTemperatures(mdf::MDFv2InMemory) = @keyoptional mdf.custom["calibTemperatures"]
calibTemperatures(mdf::MDFv2InMemory, calibTemperatures) = mdf.custom["calibTemperatures"] = calibTemperatures

auxiliaryData(mdf::MDFv2InMemory) = @keyoptional mdf.custom["auxiliaryData"]
auxiliaryData(mdf::MDFv2InMemory, auxiliaryData) = mdf.custom["auxiliaryData"] = auxiliaryData

function measDataRaw(mdf::MDFv2InMemory)
  if !(isnothing(mdf.measurement))
    return mdf.measurement.data
  else
    return mdf.measurement
  end
end

function measDataRaw(mdf::MDFv2InMemory, value)
  # Automatically create fields if they do not exist
  if isnothing(mdf.measurement)
    @debug "Creating field measurement"
    mdf.measurement = MDFv2Measurement()
  end

  mdf.measurement.data = value
end

filepath(mdf::MDFv2InMemory) = nothing # Has to be implemented...


text:
"""
Implement some ball matrix arithmetic, i.e., verified matrix operations
for (abstract) interval matrices stored in midpoint-radius format.

The 'radius' for a matrix can be either a matrix of the same dimension, or
a scalar which represents a norm.

A type parameter specifies the "norm type" for the radius. It can be `:Abs` for
entrywise radii, or `OpnormL1`, `OpnormLinf` for normwise radii (specifying which norm to use)
"""

using FastRounding

abstract type NormKind end
struct OpnormL1 <: NormKind end
struct OpnormLinf <: NormKind end
struct Abs <: NormKind end

"""
Certified upper bound to ||A||, abs(A), or whatever is used in the specified radius
"""
function normbound(A::AbstractVecOrMat{T}, ::OpnormL1) where T
    # partly taken from JuliaLang's LinearAlgebra/src/generic.jl
    Tnorm = typeof(float(real(zero(T))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for j = 1:size(A, 2)
            nrmj::Tsum = 0
            for i = 1:size(A, 1)
                nrmj = nrmj ⊕₊ abs(A[i,j])
            end
            nrm = max(nrm,nrmj)
        end
    end
    return convert(Tnorm, nrm)
end

function normbound(A::AbstractVecOrMat{T}, ::OpnormLinf) where T
    # partly taken from JuliaLang's LinearAlgebra/src/generic.jl
    Tnorm = typeof(float(real(zero(T))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for i = 1:size(A, 1)
            nrmi::Tsum = 0
            for j = 1:size(A, 2)
                nrmi = nrmi ⊕₊ abs(A[i,j])
            end
            nrm = max(nrm,nrmi)
        end
    end
    return convert(Tnorm, nrm)
end

function normbound(A::AbstractArray{T}, ::Abs) where T
    return abs.(A)
end

radiustype(MatrixType, ::OpnormLinf) = Float64
radiustype(MatrixType, ::OpnormL1) = Float64
radiustype(MatrixType, ::Abs) = MatrixType

"""
Type for ball arrays. Comes in two variants, one with cached normbound(midpoint)
(handy if you have to compute many products with it) and one without.
"""
abstract type BallArray{N<:NormKind, MatrixType<:Union{AbstractArray,Number}, RadiusType<:Union{AbstractArray, Number}} end

struct BallArrayWithoutCachedNorm{N, MatrixType, RadiusType} <: BallArray{N, MatrixType, RadiusType}
    midpoint::MatrixType
    radius::RadiusType
    function BallArrayWithoutCachedNorm{N, MatrixType, RadiusType}(midpoint, radius=zero(RadiusType)) where {N, MatrixType, RadiusType}
        if !isa(radius, radiustype(MatrixType, N()))
            error("Wrong radius type specified")
        end
        return new(midpoint, radius)
    end
end

struct BallArrayWithCachedNorm{N, MatrixType, RadiusType} <: BallArray{N, MatrixType, RadiusType}
    midpoint::MatrixType
    radius::RadiusType
    norm::RadiusType
    function BallArrayWithCachedNorm{N, MatrixType, RadiusType}(midpoint, radius=zero(RadiusType)) where {N, MatrixType, RadiusType}
        if !isa(radius, radiustype(MatrixType, N()))
            error("Wrong radius type specified")
        end
        return new{N, MatrixType, RadiusType}(midpoint, radius, normbound(midpoint, N()))
    end
end

function BallArray{N}(midpoint, radius) where {N}
    return BallArrayWithoutCachedNorm{N, typeof(midpoint), typeof(radius)}(midpoint, radius)
end

function BallArrayWithCachedNorm{N}(midpoint, radius) where {N}
    return BallArrayWithCachedNorm{N, typeof(midpoint), typeof(radius)}(midpoint, radius)
end

function midpoint_norm(A::BallArrayWithCachedNorm)
    return A.norm
end

function midpoint_norm(A::BallArrayWithoutCachedNorm{N, MT, RT}) where {N, MT, RT}
    return normbound(A, N())
end

"""
Product between ball matrices / vectors.
"""
function Base.:*(A::BallArray, B::BallArray)
    # A*maximum(B, dims=2) TODO
end

# TODO: using this file for tests

function matvec(A, v)
    w = zero(v)
    rows = rowvals(A)
    vals = nonzeros(A)
    m, n = size(A)
    for j = 1:n
       for k in nzrange(A, j)
          @inbounds i = rows[k]
          @inbounds Aij = vals[k]
          @inbounds w[i] +=  Aij * v[j]
       end
    end
end

function matvec_roundup(A, v)
    w = zero(v)
    rows = rowvals(A)
    vals = nonzeros(A)
    m, n = size(A)
    for j = 1:n
       for k in nzrange(A, j)
          @inbounds i = rows[k]
          @inbounds Aij = vals[k]
          @inbounds w[i] = w[i] ⊕₊ Aij ⊗₊ v[j]
       end
    end
end


text:
"""
    COFDetector(k = 5,
                metric = Euclidean(),
                algorithm = :kdtree,
                leafsize = 10,
                reorder = true,
                parallel = false)

Local outlier density based on chaining distance between graphs of neighbors, as described in [1].

Parameters
----------
$K_PARAM

$KNN_PARAMS

Examples
--------
$(SCORE_UNSUPERVISED("COFDetector"))

References
----------
[1] Tang, Jian; Chen, Zhixiang; Fu, Ada Wai-Chee; Cheung, David Wai-Lok (2002): Enhancing Effectiveness of Outlier
Detections for Low Density Patterns.
"""
OD.@detector mutable struct COFDetector <: UnsupervisedDetector
    k::Integer = 5::(_ > 0)
    metric::DI.Metric = DI.Euclidean()
    algorithm::Symbol = :kdtree::(_ in (:kdtree, :balltree))
    static::Union{Bool,Symbol} = :auto::(_ in (true, false, :auto))
    leafsize::Integer = 10::(_ ≥ 0)
    reorder::Bool = true
    parallel::Bool = false
end

struct COFModel <: DetectorModel
    # An efficient COF prediction requires us to store the full pairwise distance matrix of the training examples in
    # addition to the learned tree as well as the ACDs of the training examples.
    tree::NN.NNTree
    pdists::AbstractArray
    acds::AbstractVector
end

function OD.fit(detector::COFDetector, X::Data; verbosity)::Fit
    X_prep = prepare_data(X, detector.static)

    # calculate pairwise distances in addition to building the tree;
    # TODO: we could remove this once NearestNeighbors.jl exports something like `allpairs`
    pdists = DI.pairwise(detector.metric, X, dims=2)

    # use tree to calculate distances
    tree = @tree detector X_prep

    # We need k + 1 neighbors to calculate the chaining distance and have to make sure the indices are sorted
    idxs, _ = knn_others(tree, X_prep, detector.k + 1)
    acds = _calc_acds(idxs, pdists, detector.k)
    scores = _cof(idxs, acds, detector.k)
    return COFModel(tree, pdists, acds), scores
end

function OD.transform(detector::COFDetector, model::COFModel, X::Data)::Scores
    X = prepare_data(X, detector.static)

    if detector.parallel
        idxs, _ = knn_parallel(model.tree, X, detector.k + 1, true)
        return _cof(idxs, model.pdists, model.acds, detector.k)
    else
        idxs, _ = NN.knn(model.tree, X, detector.k + 1, true)
        return _cof(idxs, model.pdists, model.acds, detector.k)
    end
end

function _cof(idxs::AbstractVector{<:AbstractVector}, acds::AbstractVector, k:: Int)::Scores
    # Calculate the connectivity-based outlier factor from given acds
    cof = Vector{Float64}(undef, length(idxs))
    for (i, idx) in enumerate(idxs)
        @inbounds cof[i] = acds[i] * k / sum(acds[@view idx[2:end]])
    end
    cof
end

function _cof(idxs::AbstractVector{<:AbstractVector}, pdists::AbstractMatrix, acds:: AbstractVector, k::Int)::Scores
    # Calculate the connectivity-based outlier factor for test examples with given training distances and acds.
    cof = Vector{Float64}(undef, length(idxs))
    acdsTest = _calc_acds(idxs, pdists, k)
    for (i, idx) in enumerate(idxs)
        @inbounds cof[i] = acdsTest[i] * k / sum(acds[@view idx[2:end]])
    end
    cof
end

function _calc_acds(idxs::AbstractVector{<:AbstractVector}, pdists::AbstractMatrix, k::Int)::AbstractVector
    kplus1 = k + 1
    acds = zeros(length(idxs))
    for (i, idx) in enumerate(idxs)
        for j in 1:k
            # calculate the minimum distance (from all reachable points). That is, we sort the distances of a specific
            # point (given by idx[j+1]) according to the order of the current idx, where idx[1] specifies the idx of the
            # nearest neighbors and idx[k] specifies the idx of the k-th neighbor. We then restrict this so-called
            # set-based nearest path (SBN) to the points that are reachable with [begin:j]
            cost = minimum(pdists[idx, idx[j + 1]][begin:j])
            @inbounds acds[i] += ((2 * (kplus1 - j)) / (k * kplus1)) * cost
        end
    end
    acds
end

# """
# _cof_unused(distances, k)
# Unused version that does not rely on a pre-calculated tree, but is slower.
# """
# function _cof_unused(distances::AbstractMatrix, k:: Int)::AbstractVector
#     # pre allocate arrays
#     samples = size(distances, 2)
#     acds = zeros(samples)
#     cof = Vector{Float64}(undef, samples)
#     sbns = Array{Int, 2}(undef, (k, samples))
#     kplus1 = k + 1
#     # calculate sbns and acds
#     for i in 1:samples
#         # calculate the set based nearest path (nearest neighbors of current)
#         sbn_path = sortperm(distances[:, i], alg=PartialQuickSort(kplus1))[begin:kplus1]
#         sbns[:, i] = sbn_path[2:end]
#         for j in 1:k
#             # calculate the minimum distance (from all reachable points)
#             # e.g. point sbn_path[3] has two possible reachable points, sbn_path[1] and sbn_path[2]
#             cost = minimum(distances[sbn_path, sbn_path[j + 1]][begin:j])
#             acd = ((2 * (kplus1 - j)) / (k * kplus1)) * cost
#             acds[i] += acd
#         end
#     end
#     # calculate cof
#     for i in 1:samples
#         cof[i] = acds[i] * k / sum(acds[sbns[:, i]])
#     end
#     cof
# end