text:
@testset "desiralizers" begin
    @test LightOSM.string_deserializer(:xml) == LightXML.parse_string
    @test LightOSM.string_deserializer(:osm) == LightXML.parse_string
    @test LightOSM.string_deserializer(:json) == JSON.parse
    @test_throws ArgumentError LightOSM.string_deserializer(:other)
end

@testset "string utils" begin
    #tryparse_string_to_number
    @test LightOSM.tryparse_string_to_number("1") === 1
    @test LightOSM.tryparse_string_to_number("1.0") === 1.0
    @test LightOSM.tryparse_string_to_number("a") === "a"

    # remove_non_numeric
    @test LightOSM.remove_non_numeric("1.0") === 1.0
    @test LightOSM.remove_non_numeric("1") === 1
    @test LightOSM.remove_non_numeric("a") === 0

    # remove_sub_string_after
    @test LightOSM.remove_sub_string_after("after_this_remove", "after_this") === ""
    @test LightOSM.remove_sub_string_after("keep_after_this_remove", "after_this") === "keep_"
    @test LightOSM.remove_sub_string_after("keep_after_this_remove_after_this", "after_this") === "keep_"

    # remove_numeric
    @test LightOSM.remove_numeric("a1a") === "aa"
    @test LightOSM.remove_numeric("a1.0a") === "aa"
    @test LightOSM.remove_numeric("aa1bb2") === "aabb"
end

@testset "array tools" begin
    # trailing_elements
    @test LightOSM.trailing_elements([1]) == [1,1]
    @test LightOSM.trailing_elements([1,3]) == [1,3]
    @test LightOSM.trailing_elements([1,2,3]) == [1,3]

    # first_common_trailing_element
    @test LightOSM.first_common_trailing_element([1,2],[2,1]) == 1
    @test LightOSM.first_common_trailing_element([1,2],[2,2]) == 2
    @test LightOSM.first_common_trailing_element([1,3, 2],[2,3, 2]) == 2
    @test_throws ArgumentError LightOSM.first_common_trailing_element([1,2],[3,4]) == 3
    @test_throws ArgumentError LightOSM.first_common_trailing_element([1,5, 2],[3,5, 4]) == 3

    # join_two_arrays_on_common_trailing_elements
    @test LightOSM.join_two_arrays_on_common_trailing_elements([1,2,3],[3,4,5]) == [1,2,3,4,5]
    @test LightOSM.join_two_arrays_on_common_trailing_elements([1,2,3],[5,4,3]) == [1,2,3,4,5]
    @test LightOSM.join_two_arrays_on_common_trailing_elements([3,2,1],[5,4,3]) == [1,2,3,4,5]
    @test LightOSM.join_two_arrays_on_common_trailing_elements([3,2,1],[3,4,5]) == [1,2,3,4,5]
    @test LightOSM.join_two_arrays_on_common_trailing_elements([3,2,1],[1,4,3]) == [1,2,3,4,1]

    # join_arrays_on_common_trailing_elements
    @test LightOSM.join_arrays_on_common_trailing_elements([1,2],[2,3],[3,4]) == [1,2,3,4]
    @test LightOSM.join_arrays_on_common_trailing_elements([1,2],[4,3],[3,2]) == [1,2,3,4]
    @test LightOSM.join_arrays_on_common_trailing_elements([1,2]) == [1,2]
    @test_throws ErrorException LightOSM.join_arrays_on_common_trailing_elements([1,2],[4,3],[5,6]) == [1,2,3,4]

    # flatten
    @test LightOSM.flatten([1]) == [1]
    @test LightOSM.flatten([[1,2],[[3], [4]],1,2]) == [1,2,3,4,1,2]
end

@testset "file_deserializer" begin
    # check_valid_filename
    @test LightOSM.check_valid_filename("map.osm")
    @test LightOSM.check_valid_filename("map.json")
    @test LightOSM.check_valid_filename("map.xml")
    @test_throws ArgumentError LightOSM.check_valid_filename("map.osm.doc")
    @test_throws ArgumentError LightOSM.check_valid_filename("map.json.doc")
    @test_throws ArgumentError LightOSM.check_valid_filename("map.xml.doc")
    @test_throws ArgumentError LightOSM.check_valid_filename("map")

    # file_deserializer
    files = ["data.osm", "data.xml", "data.json", "data.doc"]
    foreach(touch, files)

    @test LightOSM.file_deserializer("data.osm") == LightXML.parse_file
    @test LightOSM.file_deserializer("data.xml") == LightXML.parse_file
    @test LightOSM.file_deserializer("data.json") == JSON.parsefile
    @test_throws ArgumentError LightOSM.file_deserializer("data.doc")

    foreach(rm, files) # tidy up
end

@testset "validate_save_location" begin
    osmfilename = "folder/filename.osm"
    @test LightOSM.validate_save_location(osmfilename, :osm) == osmfilename
    @test_throws ArgumentError LightOSM.validate_save_location(osmfilename, :xml)
    @test_throws ArgumentError LightOSM.validate_save_location(osmfilename, :json)

    jsonfilename = "folder/filename.json"
    @test LightOSM.validate_save_location(jsonfilename, :json) == jsonfilename
    @test_throws ArgumentError LightOSM.validate_save_location(jsonfilename, :osm)
    @test_throws ArgumentError LightOSM.validate_save_location(jsonfilename, :xml)

    xmlfilename = "folder/filename.xml"
    @test LightOSM.validate_save_location(xmlfilename, :xml) == xmlfilename
    @test_throws ArgumentError LightOSM.validate_save_location(xmlfilename, :osm)
    @test_throws ArgumentError LightOSM.validate_save_location(xmlfilename, :json)

    noextension = "filename"
    @test LightOSM.validate_save_location(noextension, :xml) == noextension * ".xml"
    @test LightOSM.validate_save_location(noextension, :json) == noextension * ".json"
    @test LightOSM.validate_save_location(noextension, :osm) == noextension * ".osm"
end


text:
abstract type NaiveModel end

function isfitted(model::NaiveModel)
    return model.fitted
end

isunivariate(model::NaiveModel) = true
length_observations(model::NaiveModel) = length(model.y)
observations(model::NaiveModel) = model.y
get_standard_residuals(model::NaiveModel) = model.residuals ./ sqrt(model.sigma2)
typeof_model_elements(model::NaiveModel) = eltype(model.y)

function assert_zero_missing_values(model::NaiveModel)
    for i in 1:length(model.y)
        if isnan(model.y[i])
            return error("model $(typeof(model)) does not support missing values.)")
        end
    end
    return nothing
end

@doc raw"""
    Naive(y::Vector{<:Real})

A naive model where the h step ahead forecast is
```math
y_{T+h|T} = y_T
```

# References
 * Hyndman, Rob J., Athanasopoulos, George. "Forecasting: Principles and Practice"
"""
mutable struct Naive <: NaiveModel
    y::Vector{<:Real}
    residuals::Vector{<:Real}
    sigma2::Real
    fitted::Bool

    function Naive(y::Vector{<:Real})
        Fl = eltype(y)
        n = length(y)
        return new(y, zeros(Fl, n), zero(Fl), false)
    end
end

function fit!(model::Naive)
    assert_zero_missing_values(model)
    residuals = model.y[2:end] - model.y[1:end-1]
    model.residuals = residuals
    model.sigma2 = var(residuals)
    model.fitted = true
    return model
end

function forecast(model::Naive, steps_ahead::Int; bootstrap::Bool = false)
    @assert isfitted(model)
    Fl = eltype(model.y)
    expected_value = Vector{Vector{Fl}}(undef, steps_ahead)
    covariance = Vector{Matrix{Fl}}(undef, steps_ahead)
    if bootstrap
        scenarios = simulate_scenarios(model, steps_ahead, 10_000)
        for i in 1:steps_ahead
            expected_value[i] = [mean(scenarios[i, 1, :])]
            covariance[i] = [var(scenarios[i, 1, :])][:, :]
        end
    else
        for i in 1:steps_ahead
            expected_value[i] = [model.y[end]]
            covariance[i] = [model.sigma2 * i][:, :]
        end
    end
    return Forecast{Fl}(expected_value, covariance)
end

function simulate_scenarios(model::Naive, steps_ahead::Int, n_scenarios::Int)
    Fl = typeof_model_elements(model)
    scenarios = Array{Fl, 3}(undef, steps_ahead, 1, n_scenarios)
    for s in 1:n_scenarios, i in 1:steps_ahead
        scenarios[i, 1, s] = model.y[end] + rand(model.residuals)
    end
    return scenarios
end

function reinstantiate(::Naive, y::Vector{<:Real})
    return Naive(y)
end

@doc raw"""
    SeasonalNaive(y::Vector{<:Real}, seasoanl::Int)

A seasonal naive model where the h step ahead forecast is
```math
y_{T+h|T} = y_{T + h - m(k+1)}
```
where m is the seasonal period and k is the integer part of (h-1)/m.

# References
 * Hyndman, Rob J., Athanasopoulos, George. "Forecasting: Principles and Practice"
"""
mutable struct SeasonalNaive <: NaiveModel
    y::Vector{<:Real}
    residuals::Vector{<:Real}
    seasonal::Int
    sigma2::Real
    fitted::Bool

    function SeasonalNaive(y::Vector{<:Real}, seasonal::Int)
        Fl = eltype(y)
        n = length(y)
        return new(y, zeros(Fl, n), seasonal, zero(Fl), false)
    end
end

function fit!(model::SeasonalNaive)
    assert_zero_missing_values(model)
    residuals = model.y[model.seasonal+1:end] - model.y[1:end-model.seasonal]
    model.residuals = residuals
    model.sigma2 = var(residuals)
    model.fitted = true
    return model
end

function forecast(model::SeasonalNaive, steps_ahead::Int; bootstrap::Bool = false)
    @assert isfitted(model)
    Fl = eltype(model.y)
    expected_value = Vector{Vector{Fl}}(undef, steps_ahead)
    covariance = Vector{Matrix{Fl}}(undef, steps_ahead)
    if bootstrap
        scenarios = simulate_scenarios(model, steps_ahead, 10_000)
        for i in 1:steps_ahead
            expected_value[i] = [mean(scenarios[i, 1, :])]
            covariance[i] = [var(scenarios[i, 1, :])][:, :]
        end
    else
        for i in 1:steps_ahead
            if i in 1:model.seasonal
                expected_value[i] = [model.y[end - model.seasonal + i]]
            else
                expected_value[i] = copy(expected_value[i - model.seasonal])
            end
            covariance[i] = [model.sigma2 * (fld(i - 1, model.seasonal) + 1)][:, :]
        end
    end
    return Forecast{Fl}(expected_value, covariance)
end

function simulate_scenarios(model::SeasonalNaive, steps_ahead::Int, n_scenarios::Int)
    Fl = typeof_model_elements(model)
    scenarios = Array{Fl, 3}(undef, steps_ahead, 1, n_scenarios)
    for s in 1:n_scenarios, i in 1:steps_ahead
        scenarios[i, 1, s] = model.y[end - model.seasonal + mod(i, model.seasonal)] + rand(model.residuals)
    end
    return scenarios
end

function reinstantiate(model::SeasonalNaive, y::Vector{<:Real})
    return SeasonalNaive(y, model.seasonal)
end

@doc raw"""
    ExperimentalSeasonalNaive(y::Vector{<:Real}, seasonal::Int; S::Int = 10_000)

A seasonal naive model where the h step ahead forecast is the mean of the simulation
of S scenarios
```math
y_{T+h|T} = y_{T + h - m(k+1)} + \varepsilon_t
```
where m is the seasonal period, k is the integer part of (h-1)/m and ``\varepsilon_t`` is
a sampled error.

We call it experimental because so far we could not find a good reference and implementation.
If you know something please post it as an issue.
"""
mutable struct ExperimentalSeasonalNaive <: NaiveModel
    y::Vector{<:Real}
    residuals::Vector{<:Real}
    seasonal::Int
    sigma2::Real
    fitted::Bool
    S::Int

    function ExperimentalSeasonalNaive(y::Vector{<:Real}, seasonal::Int; S::Int = 10_000)
        Fl = eltype(y)
        n = length(y)
        return new(y, zeros(Fl, n), seasonal, zero(Fl), false, S)
    end
end

function fit!(model::ExperimentalSeasonalNaive)
    assert_zero_missing_values(model)
    residuals = model.y[model.seasonal+1:end] - model.y[1:end-model.seasonal]
    model.residuals = residuals
    model.sigma2 = var(residuals)
    model.fitted = true
    return model
end

function forecast(model::ExperimentalSeasonalNaive, steps_ahead::Int)
    @assert isfitted(model)
    Fl = eltype(model.y)
    expected_value = Vector{Vector{Fl}}(undef, steps_ahead)
    covariance = Vector{Matrix{Fl}}(undef, steps_ahead)
    scenarios = simulate_scenarios(model, steps_ahead, model.S)
    for i in 1:steps_ahead
        expected_value[i] = [mean(scenarios[i, 1, :])]
        covariance[i] = [var(scenarios[i, 1, :])][:, :]
    end
    return Forecast{Fl}(expected_value, covariance)
end

function simulate_scenarios(model::ExperimentalSeasonalNaive, steps_ahead::Int, n_scenarios::Int)
    Fl = typeof_model_elements(model)
    scenarios = Array{Fl, 3}(undef, steps_ahead, 1, n_scenarios)
    for i in 1:steps_ahead
        if i <= model.seasonal
            for s in 1:n_scenarios
                scenarios[i, 1, s] = model.y[end - model.seasonal + i] + rand(model.residuals)
            end
        else
            for s in 1:n_scenarios
                scenarios[i, 1, s] = scenarios[i - model.seasonal, 1, s] + rand(model.residuals)
            end
        end
    end
    return scenarios
end

function reinstantiate(model::ExperimentalSeasonalNaive, y::Vector{<:Real})
    return ExperimentalSeasonalNaive(y, model.seasonal; S = model.S)
end

function cross_validation(model::NaiveModel, steps_ahead::Int, start_idx::Int;
                  n_scenarios::Int = 10_000)
    Fl = typeof_model_elements(model)
    num_fits = length(model.y) - start_idx - steps_ahead
    b = CrossValidation{Fl}(num_fits, steps_ahead)
    for i in 1:num_fits
        println("CrossValidation: step $i of $num_fits")
        y_to_fit = model.y[1:start_idx - 1 + i]
        y_to_verify = model.y[start_idx + i:start_idx - 1 + i + steps_ahead]
        model_to_fit = reinstantiate(model, y_to_fit)
        fit!(model_to_fit)
        forec = forecast(model_to_fit, steps_ahead)
        scenarios = simulate_scenarios(model_to_fit, steps_ahead, n_scenarios)
        expected_value_vector = forecast_expected_value(forec)[:]
        abs_errors = evaluate_abs_error(y_to_verify, expected_value_vector)
        crps_scores = evaluate_crps(y_to_verify, scenarios[:, 1, :])
        b.abs_errors[:, i] = abs_errors
        b.crps_scores[:, i] = crps_scores
    end
    for i in 1:steps_ahead
        b.mae[i] = mean(b.abs_errors[i, :])
        b.mean_crps[i] = mean(b.crps_scores[i, :])
    end
    return b
end

ext: jl


text:
# from Posit to Float64,32,16
Float64(x::Posit8) = ccall((:convertP8ToDouble, SoftPositPath), Float64, (Posit8,),x)
Float64(x::Posit32) = ccall((:convertP32ToDouble, SoftPositPath), Float64, (Posit32,),x)

# conversion from PX2 to Float64,32,16
Float64(x::Posit8_2) = ccall((:convertPX2ToDouble, SoftPositPath), Float64, (Posit8_2,),x)
Float64(x::Posit16_2) = ccall((:convertPX2ToDouble, SoftPositPath), Float64, (Posit16_2,),x)
Float64(x::Posit24_2) = ccall((:convertPX2ToDouble, SoftPositPath), Float64, (Posit24_2,),x)

# conversion from PX1 to Float64,32,16
Float64(x::Posit8_1) = ccall((:convertPX1ToDouble, SoftPositPath), Float64, (Posit8_1,),x)
Float64(x::Posit16_1) = ccall((:convertPX1ToDouble, SoftPositPath), Float64, (Posit16_1,),x)
Float64(x::Posit24_1) = ccall((:convertPX1ToDouble, SoftPositPath), Float64, (Posit24_1,),x)

# conversion to Float32,16
Float16(x::AbstractPosit) = Float16(Float64(x))
Float32(x::AbstractPosit) = Float32(Float64(x))
Float16(x::Posit16) = Float16(Float32(x))
Float64(x::Posit16) = Float64(Float32(x))

function Float32(x::Posit16)
    ui = reinterpret(UInt16,x)

    signbitx = signbit(x)                   # sign of number
    abs_p16 = (signbitx ? -ui : ui) << 1    # two's complement for negative, remove sign bit

    # determine exponent sign from 2nd bit: 1=x>
    sign_exponent = reinterpret(Bool,((abs_p16 & 0x8000) >> 15) % UInt8)
    n_regimebits = sign_exponent ? leading_ones(abs_p16) : leading_zeros(abs_p16)

    # extract exponent and mantissa bits
    exponent_bit = ((abs_p16 << (n_regimebits+1)) & 0x8000) >> 15
    mantissa = ((abs_p16 << (n_regimebits + 2)) % UInt32) << 7

    # assemble float exponent from posit # of regime bits and exponent bit
    exponent = ((sign_exponent ? n_regimebits : -n_regimebits+1) << 1) + exponent_bit + 125

    # set exponent (and 1st mantissa bit) to NaN for NaR inputs
    # set exponent to 0 for zero(Posit16) input
    exponent = n_regimebits == 16 ? signbitx*0x7fc00000 : (exponent % UInt32) << 23

    # assemble sign, exponent and mantissa bits
    sign = (signbitx % UInt32) << 31		# isolate sign bit
    f32 = sign | exponent | mantissa		# concatenate sign, exponent and mantissa
    return reinterpret(Float32,f32)
end

# legacy
Float32_old(x::AbstractPosit) = Float32(Float64_old(x))
Float64_old(x::Posit16) = ccall((:convertP16ToDouble, SoftPositPath), Float64, (Posit16,),x)


text:
module Book

using ..Chess, ..Chess.PGN, Dates, Printf, StatsBase

export BookEntry

export createbook,
    findbookentries, pickbookmove, printbookentries, purgebook, writebooktofile

"""
    BookEntry

A struct representing an opening book entry.

Book entries contain the following slots:

- `key`: The hash key of the board position this book entry represents.
- `move`: The move played, encoded as an `Int32`. In order to get the actual
  `Move` value representing the move stored in a book entry `e`, you should do
  `Move(e.move)`.
- `elo`: The highest Elo rating of a player who played this move.
- `oppelo`: The highest Elo of the opponent in a game where this move was
  played.
- `wins`: The number of times the player who played this move won the game.
- `draws`: The number of times the player who played this move drew the game.
- `losses`: The number of times the player who played this move lost the game.
- `firstyear`: The year this move was first played.
- `lastyear`: The year this move was last played.
- `score`: The score of this move, used to obtain a probability distribution
  when picking a book move for a position. The score is computed based on the
  W/L/D stats for the move, the ratings of the players who have played it, and
  on its popularity in more recent games.
"""
mutable struct BookEntry
    key::UInt64
    move::Int32
    elo::Int16
    oppelo::Int16
    wins::Int32
    draws::Int32
    losses::Int32
    firstyear::Int16
    lastyear::Int16
    score::Float32
end

function BookEntry()
    BookEntry(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
end

function Base.show(io::IO, e::BookEntry)
    println(io, "BookEntry:")
    println(io, "  key: $(e.key)")
    println(io, "  move: $(Move(e.move))")
    println(io, "  top elo: $(e.elo)")
    println(io, "  top opponent elo: $(e.oppelo)")
    println(io, "  (w, d, l): ($(e.wins), $(e.draws), $(e.losses))")
    println(io, "  first played: $(e.firstyear)")
    println(io, "  last played: $(e.lastyear)")
    println(io, "  score: $(e.score)")
end

const ENTRY_SIZE = 36
const COMPACT_ENTRY_SIZE = 16

function entrytobytes(entry::BookEntry, compact::Bool)::Vector{UInt8}
    io = IOBuffer(
        UInt8[],
        read = true,
        write = true,
        maxsize = compact ? COMPACT_ENTRY_SIZE : ENTRY_SIZE,
    )
    write(io, entry.key)
    write(io, entry.move)
    write(io, entry.score)
    if !compact
        write(io, entry.elo)
        write(io, entry.oppelo)
        write(io, entry.wins)
        write(io, entry.draws)
        write(io, entry.losses)
        write(io, entry.firstyear)
        write(io, entry.lastyear)
    end
    take!(io)
end

function entryfrombytes(bytes::Vector{UInt8}, compact::Bool)::BookEntry
    io = IOBuffer(bytes)
    result = BookEntry()
    result.key = read(io, UInt64)
    result.move = read(io, Int32)
    result.score = read(io, Float32)
    if !compact
        result.elo = read(io, Int16)
        result.oppelo = read(io, Int16)
        result.wins = read(io, Int32)
        result.draws = read(io, Int32)
        result.losses = read(io, Int32)
        result.firstyear = read(io, Int16)
        result.lastyear = read(io, Int16)
    end
    result
end

function computescore(
    result,
    color,
    elo,
    date,
    scorewhitewin,
    scorewhitedraw,
    scorewhiteloss,
    scoreblackwin,
    scoreblackdraw,
    scoreblackloss,
    scoreunknown,
    highelofactor,
    yearlydecay,
)::Float32
    base = if result == "1-0" && color == WHITE
        scorewhitewin
    elseif result == "1/2-1/2" && color == WHITE
        scorewhitedraw
    elseif result == "0-1" && color == WHITE
        scorewhiteloss
    elseif result == "0-1" && color == BLACK
        scoreblackwin
    elseif result == "1/2-1/2" && color == BLACK
        scoreblackdraw
    elseif result == "1-0" && color == BLACK
        scoreblackloss
    else
        scoreunknown
    end
    base *
    max(1.0, 0.01 * highelofactor * (2300 - elo)) *
    exp(log(yearlydecay) * (Dates.value(today() - date) / 365.25))
end

function mergeable(e1::BookEntry, e2::BookEntry)::Bool
    e1.key == e2.key && e1.move == e2.move
end

function merge(e1::BookEntry, e2::BookEntry)::BookEntry
    BookEntry(
        e1.key,
        e1.move,
        max(e1.elo, e2.elo),
        max(e1.oppelo, e2.oppelo),
        e1.wins + e2.wins,
        e1.draws + e2.draws,
        e1.losses + e2.losses,
        min(e1.firstyear, e2.firstyear),
        max(e1.lastyear, e2.lastyear),
        e1.score + e2.score,
    )
end

function compress!(entries::Vector{BookEntry})
    i = 1
    j = 1
    n = length(entries)
    iterations = 0
    while j + 1 < n
        for k ∈ j+1:n
            if !mergeable(entries[i], entries[k]) || k == n
                j = k
                i += 1
                entries[i] = entries[k]
                break
            end
            entries[i] = merge(entries[i], entries[k])
        end
    end
    entries[1:i-1]
end

function compareentries(e1::BookEntry, e2::BookEntry)::Bool
    e1.key < e2.key || (e1.key == e2.key && e1.move < e2.move)
end

function sortentries!(entries::Vector{BookEntry})
    sort!(entries, lt = compareentries)
end

function addgame!(
    entries::Vector{BookEntry},
    g::SimpleGame,
    scorewhitewin,
    scorewhitedraw,
    scorewhiteloss,
    scoreblackwin,
    scoreblackdraw,
    scoreblackloss,
    scoreunknown,
    highelofactor,
    yearlydecay,
    maxply,
    minelo,
)
    result = headervalue(g, "Result")
    if result ≠ "*"
        w = result == "1-0" ? 1 : 0
        d = result == "1/2-1/2" ? 1 : 0
        l = result == "0-1" ? 1 : 0
        welo = !isnothing(whiteelo(g)) ? whiteelo(g) : 0
        belo = !isnothing(blackelo(g)) ? blackelo(g) : 0
        date = !isnothing(dateplayed(g)) ? dateplayed(g) : Date(1900, 1, 1)
        year = Dates.year(date)
        wscore = computescore(
            result,
            WHITE,
            welo,
            date,
            scorewhitewin,
            scorewhitedraw,
            scorewhiteloss,
            scoreblackwin,
            scoreblackdraw,
            scoreblackloss,
            scoreunknown,
            highelofactor,
            yearlydecay,
        )
        bscore = computescore(
            result,
            BLACK,
            belo,
            date,
            scorewhitewin,
            scorewhitedraw,
            scorewhiteloss,
            scoreblackwin,
            scoreblackdraw,
            scoreblackloss,
            scoreunknown,
            highelofactor,
            yearlydecay,
        )

        tobeginning!(g)
        while !isatend(g) && g.ply <= maxply
            b = board(g)
            wtm = sidetomove(b) == WHITE
            m = nextmove(g)
            if wtm && welo >= minelo || !wtm && belo >= minelo
                push!(
                    entries,
                    BookEntry(
                        b.key,
                        m.val,
                        wtm ? welo : belo,
                        wtm ? belo : welo,
                        wtm ? w : l,
                        d,
                        wtm ? l : w,
                        year,
                        year,
                        wtm ? wscore : bscore,
                    ),
                )
            end
            forward!(g)
        end
    end
end

function addgamefile!(
    entries::Vector{BookEntry},
    filename::String,
    scorewhitewin,
    scorewhitedraw,
    scorewhiteloss,
    scoreblackwin,
    scoreblackdraw,
    scoreblackloss,
    scoreunknown,
    highelofactor,
    yearlydecay,
    maxply,
    minelo,
    count,
)
    for g ∈ PGN.gamesinfile(filename)
        addgame!(
            entries,
            g,
            scorewhitewin,
            scorewhitedraw,
            scorewhiteloss,
            scoreblackwin,
            scoreblackdraw,
            scoreblackloss,
            scoreunknown,
            highelofactor,
            yearlydecay,
            maxply,
            minelo,
        )
        count += 1
        if count % 1000 == 0
            println("$count games added, $(length(entries)) entries.")
        end
    end
    count
end

"""
    createbook(pgnfiles::Vararg{String};
               scorewhitewin = 8.0,
               scorewhitedraw = 4.0,
               scorewhiteloss = 1.0,
               scoreblackwin = 8.0,
               scoreblackdraw = 5.0,
               scoreblackloss = 1.0,
               scoreunknown = 0.0,
               highelofactor = 6.0,
               yearlydecay = 0.85,
               maxply = 60,
               minelo = 0)

Creates an opening book tree from one or more PGN files.

The opening tree is stored in RAM. You will probably want to save it to disk
using `writebooktofile` afterwards, for instance like this:

```julia-repl
julia> bk = createbook("my-pgn-file.pgn");

julia> writebooktofile(bk, "my-book.obk")
```

The createbook function takes a number of optional named parameters that
can be used to control what moves are included in the opening tree, and the
scoring of the moves (which is used to produce move probabilities when picking
a move using `pickbookmove`). These are:

- `scorewhitewin`: The base score for all white moves in a game won by white.
- `scorewhitedraw`: The base score for all white moves in a drawn game.
- `scorewhiteloss`: The base score for all black moves in a game won by black.
- `scoreblackwin`: The base score for all black moves in a game won by black.
- `scoreblackdraw`: The base score for all black moves in a drawn game.
- `scoreblackloss`: The base score for all black moves in a game won by white.
- `scoreunknown`: The base score for all moves in a game with an unknown result.
- `highelofactor`: Score multiplier for moves played by a player with high
  rating. The base score is multiplied by
  `max(1.0 0.01 * highelofactor * (2300 - elo))`
- `yearlydecay`: Controls exponential yearly reduction of scores. If a game was
  played `n` years ago, all scores are multiplied by `n^yearlydecay`.
- `maxply`: Maximum depth of the opening tree. If `maxply` equals 60 (the
  default), no moves after move 30 are included in the opening tree.
- `minelo`: Minimum Elo for book moves. Moves played by players below this
  number are not included in the opening tree.
"""
function createbook(
    pgnfiles::Vararg{String};
    scorewhitewin = 8.0,
    scorewhitedraw = 4.0,
    scorewhiteloss = 1.0,
    scoreblackwin = 8.0,
    scoreblackdraw = 5.0,
    scoreblackloss = 1.0,
    scoreunknown = 0.0,
    highelofactor = 6.0,
    yearlydecay = 0.85,
    maxply = 60,
    minelo = 0,
)
    result = Vector{BookEntry}()
    count = 0
    for filename ∈ pgnfiles
        count = addgamefile!(
            result,
            filename,
            scorewhitewin,
            scorewhitedraw,
            scorewhiteloss,
            scoreblackwin,
            scoreblackdraw,
            scoreblackloss,
            scoreunknown,
            highelofactor,
            yearlydecay,
            maxply,
            minelo,
            count,
        )
    end
    result = compress!(sortentries!(result))
    result
end

"""
    writebooktofile(entries::Vector{BookEntry}, filename::String,
                    compact = false)

Writes a book (as created by `createbookfile`) to a binary file.

If the optional parameter `compact` is `true`, the book is written in a more
compact format that does not include W/L/D counts, Elo numbers and years.
"""
function writebooktofile(entries::Vector{BookEntry}, filename::String, compact = false)
    open(filename, "w") do f
        write(f, UInt8(compact ? 1 : 0))
        for e ∈ entries
            write(f, entrytobytes(e, compact))
        end
    end
end

"""
    purgebook(infilename::String, outfilename::String;
              minscore = 0, mingamecount = 5, compact = false)

Creates a smaller version of an opening book file by removing unimportant lines.

Book moves with score lower than `minscore` or which have been played in fewer
than `mingamecount` games are not included in the output file.

If the optional parameter `compact` is `true`, the output file is written in a
more compact format that does not include W/L/D counts, Elo numbers and years.
"""
function purgebook(
    infilename::String,
    outfilename::String;
    minscore = 0,
    mingamecount = 5,
    compact = false,
)
    open(infilename, "r") do inf
        open(outfilename, "w") do outf
            write(outf, UInt8(compact ? 1 : 0))
            incompact = read(inf, UInt8) == 1
            entrysize = incompact ? COMPACT_ENTRY_SIZE : ENTRY_SIZE
            entrycount = div(filesize(infilename) - 1, entrysize)
            for i = 1:entrycount
                e = readentry(inf, i - 1, incompact)
                if e.wins + e.draws + e.losses ≥ mingamecount && e.score > minscore
                    write(outf, entrytobytes(e, compact))
                end
                if i % 100000 == 0
                    @printf(
                        "%d/%d (%.1f%%) entries processed.\n",
                        i,
                        entrycount,
                        100 * i / entrycount
                    )
                end
            end
        end
    end
end

function readkey(f::IO, index::Int, entrysize::Int)::UInt64
    seek(f, 1 + entrysize * index)
    read(f, UInt64)
end

function findkey(
    f::IO,
    key::UInt64,
    left::Int,
    right::Int,
    entrysize::Int,
)::Union{Int,Nothing}
    while left <= right
        middle = div(left + right, 2)
        midkey = readkey(f, middle, entrysize)
        if key == midkey && (middle == 0 || key ≠ readkey(f, middle - 1, entrysize))
            return middle
        elseif midkey < key
            left = middle + 1
        else
            right = middle - 1
        end
    end
end

function readentry(f::IO, index::Int, compact = false)::Union{BookEntry,Nothing}
    entrysize = compact ? COMPACT_ENTRY_SIZE : ENTRY_SIZE
    seek(f, 1 + entrysize * index)
    e = entryfrombytes(read(f, entrysize), compact)
end

"""
    findbookentries(b::Board, bookfilename::String)
    findbookentries(key::UInt64, bookfilename::String)

Returns all book entries for the given board or key.

The return value is a (possibly empty) `Vector{BookEntry}`, sorted by
descending scores.
"""
function findbookentries(key::UInt64, bookfilename::String)::Vector{BookEntry}
    result = Vector{BookEntry}()
    open(bookfilename, "r") do f
        compact = read(f, UInt8) == 1
        entrysize = compact ? COMPACT_ENTRY_SIZE : ENTRY_SIZE
        entrycount = div(filesize(bookfilename) - 1, entrysize)
        i = findkey(f, key, 0, entrycount - 1, entrysize)
        if !isnothing(i)
            for j = i:entrycount
                e = readentry(f, j, compact)
                if e.key == key
                    push!(result, e)
                else
                    break
                end
            end
        end
    end
    sort(result, by = e -> -e.score)
end

function findbookentries(b::Board, bookfilename::String)::Vector{BookEntry}
    findbookentries(b.key, bookfilename)
end

"""
    printbookentries(b::Board, bookfilename::String)

Pretty-print the move entries for the provided board.
"""
function printbookentries(b::Board, bookfilename::String)
    entries = findbookentries(b, bookfilename)
    scoresum = sum(map(e -> e.score, entries))
    for e ∈ entries
        @printf(
            "%s %.3f %.1f%% (+%d, =%d, -%d) %d %d %d %d\n",
            movetosan(b, Move(e.move)),
            e.score / scoresum,
            100 * ((e.wins + 0.5 * e.draws) / (e.wins + e.draws + e.losses)),
            e.wins,
            e.draws,
            e.losses,
            e.elo,
            e.oppelo,
            e.firstyear,
            e.lastyear
        )
    end
end

"""
    pickbookmove(b::Board, bookfilename::String;
                 minscore = 0, mingamecount = 1)

Picks a book move for the board `b`, returning `nothing` when out of book.

The move is selected with probabilities given by the `score` slots in the
`BookEntry` objects. The `minscore` and `mingamecount` parameters can be used
to exclude moves with low score or low play counts.
"""
function pickbookmove(
    b::Board,
    bookfilename::String;
    minscore = 0,
    mingamecount = 1,
)::Union{Move,Nothing}
    entries = filter(
        e -> e.wins + e.draws + e.losses >= mingamecount && e.score >= minscore,
        findbookentries(b, bookfilename),
    )
    if length(entries) == 0
        nothing
    else
        w = Weights(map(e -> e.score, entries))
        Move(sample(entries, w).move)
    end
end

end # module


text:
@test coord_distance(1.0, -1.0, 2.0, -2.0) > 0
@test coord_distance(1.0, -1.0, 2.0, -2.0, calcfunc=vincenty, radius=EARTH_RADIUS_EQUATORIAL) > 0
@test coord_distance(1.0, -1.0, 2.0, -2.0, calcfunc=haversine) > 0
@test coord_distance((1., -1.),(2.,-2.), calcfunc=vincenty) > 0
@test haversine(1., -1., 2., -2., EARTH_RADIUS_POLAR) > 0
@test flatpythagorean(1., -1., 2., -2., EARTH_RADIUS_EQUATORIAL) > 0
@test vincenty(1., -1., 2., -2., EARTH_RADIUS_POLAR) > 0
@test coord_distance(1.0, -1.0, 2.0, -2.0, units=u"mi") > 0
@test coord_distance(1.0, -1.0, 2.0, -2.0, units=u"mi", returntype=Rational{Int64}) > 0


text:
# tunelogistic uses the logistic function to scale the acceptance rate by a factor ranging from 0 to 2
# k gives the curve's steepness. For larger k, the curve becomes more steep
function tunelogistic(x::Real,k::Real)
    r = zero(x)
    if x < 0
        r = one(x)/(one(x)+exp(-k*x)) + one(x)/2
    else
        r = 2one(x)/(one(x)+exp(-k*x))
    end
    sqrt(r)
end

# tuneerf uses the error function (erf) to scale the acceptance rate by a factor ranging from 0 to 2
# k gives the curve's steepness. For larger k, the curve becomes more steep
function tuneerf(x::Real,k::Real)
    r = zero(x)
    if x < 0
        r = erf(k*x)/2 + one(x)
    else
        r = erf(k*x) + one(x)
    end
    sqrt(r)
end


text:
# TODO: add a type for kinetic energy
# TODO: add cache for gradients by letting ∂logπ∂θ return both log density and gradient

struct Hamiltonian{M<:AbstractMetric, Tlogπ, T∂logπ∂θ}
    metric::M
    ℓπ::Tlogπ
    ∂ℓπ∂θ::T∂logπ∂θ
end
Base.show(io::IO, h::Hamiltonian) = print(io, "Hamiltonian(metric=$(h.metric))")

struct DualValue{Tv<:AbstractFloat, Tg<:AbstractVector{Tv}}
    value::Tv    # Cached value, e.g. logπ(θ).
    gradient::Tg # Cached gradient, e.g. ∇logπ(θ).
end

# `∂H∂θ` now returns `(logprob, -∂ℓπ∂θ)`
function ∂H∂θ(h::Hamiltonian, θ::AbstractVector)
    res = h.∂ℓπ∂θ(θ)
    return DualValue(res[1], -res[2])
end

∂H∂r(h::Hamiltonian{<:UnitEuclideanMetric}, r::AbstractVector) = copy(r)
∂H∂r(h::Hamiltonian{<:DiagEuclideanMetric}, r::AbstractVector) = h.metric.M⁻¹ .* r
∂H∂r(h::Hamiltonian{<:DenseEuclideanMetric}, r::AbstractVector) = h.metric.M⁻¹ * r

struct PhasePoint{T<:AbstractVector, V<:DualValue}
    θ::T  # Position variables / model parameters.
    r::T  # Momentum variables
    ℓπ::V # Cached neg potential energy for the current θ.
    ℓκ::V # Cached neg kinect energy for the current r.
    function PhasePoint(θ::T, r::T, ℓπ::V, ℓκ::V) where {T,V}
        @argcheck length(θ) == length(r) == length(ℓπ.gradient) == length(ℓπ.gradient)
        if any(isfinite.((θ, r, ℓπ, ℓκ)) .== false)
            @warn "The current proposal will be rejected due to numerical error(s)." isfinite.((θ, r, ℓπ, ℓκ))
            ℓπ = DualValue(-Inf, ℓπ.gradient)
            ℓκ = DualValue(-Inf, ℓκ.gradient)
        end
        new{T,V}(θ, r, ℓπ, ℓκ)
    end
end

phasepoint(
    h::Hamiltonian,
    θ::T,
    r::T;
    ℓπ=∂H∂θ(h, θ),
    ℓκ=DualValue(neg_energy(h, r, θ), ∂H∂r(h, r))
) where {T<:AbstractVector} = PhasePoint(θ, r, ℓπ, ℓκ)

# If position variable and momentum variable are in different containers,
# move the momentum variable to that of the position variable.
# This is needed for AHMC to work with CuArrays (without depending on it).
phasepoint(
    h::Hamiltonian,
    θ::T1,
    _r::T2;
    r=T1(_r),
    ℓπ=∂H∂θ(h, θ),
    ℓκ=DualValue(neg_energy(h, r, θ), ∂H∂r(h, r))
) where {T1<:AbstractVector,T2<:AbstractVector} = PhasePoint(θ, r, ℓπ, ℓκ)

Base.isfinite(v::DualValue) = all(isfinite, v.value) && all(isfinite, v.gradient)
Base.isfinite(v::AbstractVector) = all(isfinite, v)
Base.isfinite(z::PhasePoint) = isfinite(z.ℓπ) && isfinite(z.ℓκ)

###
### Negative energy (or log probability) functions.
### NOTE: the general form (i.e. non-Euclidean) of K depends on both θ and r.
###

neg_energy(z::PhasePoint) = z.ℓπ.value + z.ℓκ.value

neg_energy(h::Hamiltonian, θ::AbstractVector) = h.ℓπ(θ)

neg_energy(
    h::Hamiltonian{<:UnitEuclideanMetric},
    r::T,
    θ::T
) where {T<:AbstractVector} = -sum(abs2, r) / 2

function neg_energy(
    h::Hamiltonian{<:DiagEuclideanMetric},
    r::T,
    θ::T
) where {T<:AbstractVector}
    _r = [abs2(r[i]) * h.metric.M⁻¹[i] for i in 1:length(r)]
    return -sum(_r) / 2
end

function neg_energy(
    h::Hamiltonian{<:DenseEuclideanMetric},
    r::T,
    θ::T
) where {T<:AbstractVector}
    mul!(h.metric._temp, h.metric.M⁻¹, r)
    return -dot(r, h.metric._temp) / 2
end

energy(args...) = -neg_energy(args...)

####
#### Momentum refreshment
####

phasepoint(
    rng::AbstractRNG,
    θ::AbstractVector{T},
    h::Hamiltonian
) where {T<:Real} = phasepoint(h, θ, rand(rng, h.metric))

refresh(
    rng::AbstractRNG,
    z::PhasePoint,
    h::Hamiltonian
) = phasepoint(h, z.θ, rand(rng, h.metric))


text:
"""
    Lennard-Jones Potential
"""
mutable struct LennardJones <: Potential
    ε::Float64
    σ::Float64
end

"""
    LennardJones(params::Dict)

Creates the LennardJones potential.
"""
function LennardJones(params::Dict)
    # Creates the LennardJones model
    ε = params["ε"]
    σ = params["σ"]
    return LennardJones(ε, σ)
end

"""
    potential_energy(p::LennardJones, r::Position, args...)

Calculates LennardJones potential energy.
"""
function potential_energy(p::LennardJones, r::Position, args...)
    return 4.0 * p.ε * ((p.σ / norm(r))^12 - (p.σ / norm(r))^6)
end


text:
using Test, LinearMaps, LinearAlgebra, Quaternions

# type piracy because Quaternions.jl doesn't have it right
Base.:(*)(z::Complex, q::Quaternion) = quat(z) * q
Base.:(*)(q::Quaternion, z::Complex) = q * quat(z)
Base.:(+)(q::Quaternion, z::Complex) = q + quat(z)

@testset "noncommutative number type" begin
    x = Quaternion.(rand(10), rand(10), rand(10), rand(10))
    v = rand(10)
    A = Quaternion.(rand(10,10), rand(10,10), rand(10,10), rand(10,10))
    B = rand(ComplexF64, 10, 10)
    C = similar(A)
    γ = Quaternion.(rand(4)...) # "Number"
    α = UniformScaling(γ)
    β = UniformScaling(Quaternion.(rand(4)...))
    λ = rand(ComplexF64)
    L = LinearMap(A)
    F = LinearMap{eltype(A)}(x -> A*x, y -> A'y, 10)
    @test Array(F) == A
    @test Array(F') == A'
    @test Array(transpose(F)) == transpose(A)
    @test Array(α * F) == α * A
    @test Array(F * α) == A * α
    @test Array(α * F) == α * A
    @test Array(F * α ) == A * α
    @test Array(α * F') == α * A'
    for M in (L, F)
        @test mul!(C, transpose(A), M) ≈ transpose(A)*A
        @test mul!(C, A', M) ≈ A'A
        @test mul!(C, A, M) ≈ A*A
        @test mul!(copy(C), M, A, γ, λ) ≈ A*A*γ + C*λ
        @test mul!(copy(C), A, M, γ, λ) ≈ A*A*γ + C*λ
        @test mul!(copy(C), A, M, γ, 0) ≈ A*A*γ
        @test mul!(copy(C), transpose(A), M, γ, λ) ≈ transpose(A)*A*γ + C*λ
        @test mul!(copy(C), adjoint(A), M, γ, λ) ≈ A'*A*γ + C*λ
    end
    @test Array((α * F')') ≈ (γ * A')' ≈ A * conj(γ)
    @test L * x ≈ A * x
    @test L' * x ≈ A' * x
    @test α * (L * x) ≈ γ * (A * x)
    @test α * L * x ≈ γ * A * x
    @test L * α * x ≈ A * γ * x
    @test 3L * x ≈ 3A * x
    @test 3L' * x ≈ 3A' * x
    @test λ*L isa LinearMaps.CompositeMap
    @test γ * (λ * LinearMap(B)) isa LinearMaps.CompositeMap
    @test (λ * LinearMap(B)) * γ isa LinearMaps.CompositeMap
    @test λ*L * x ≈ λ*A * x
    @test λ*L' * x ≈ λ*A' * x
    @test α * (3L * x) ≈ γ * (3A * x)
    @test (@inferred α * 3L) * x ≈ γ * 3A * x
    @test (@inferred 3L * α) * x ≈ 3A * γ * x
    @test (α * L') * x ≈ (γ * A') * x
    @test (α * L')' * x ≈ (γ * A')' * x
    @test (α * L')' * v ≈ (γ * A')' * v
    @test Array(@inferred adjoint(α * L * β)) ≈ conj(β) * A' * conj(γ)
    @test Array(@inferred transpose(α * L * β)) ≈ β * transpose(A) * γ
    J = LinearMap(α, 10)
    @test (β * J) * x ≈ LinearMap(β*α, 10) * x ≈ β*γ*x
    @test (J * β) * x ≈ LinearMap(α*β, 10) * x ≈ γ*β*x
    M = β.λ * (γ * L * L)
    @test M == β * (γ * L * L)
    @test length(M.maps) == 3
    @test M.maps[end].λ == β.λ * γ
    @test γ * (β * L * L) == γ * (β.λ * L * L) == α * (β.λ * L * L) == α * (β * L * L)
    @test length(M.maps) == 3
    @test M.maps[end].λ == β.λ * γ
    M = (L * L * γ) * β.λ
    @test M == (L * L * γ) * β == (L * L * α) * β == (L * L * α) * β.λ
    @test length(M.maps) == 3
    @test M.maps[1].λ == γ*β.λ
    @test γ*FillMap(γ, (3, 4)) == FillMap(γ^2, (3, 4)) == FillMap(γ, (3, 4))*γ
    U = LinearMap(quat(1.0)*I, 10)
    for β in (0, 1, rand())
        @test mul!(copy(x), J, x, γ, β) == γ * x * γ + x * β
        @test mul!(copy(x), U, x, γ, β) == x * γ + x * β
    end
    # exercise non-RealOrComplex scalar operations
    @test Array(γ * (L'*L)) ≈ γ * (A'*A) # CompositeMap
    @test Array((L'*L) * γ) ≈ (A'*A) * γ
    @test Array(-L) == -A
    @test Array(γ \ L) ≈ γ \ A
    @test Array(L / γ) ≈ A / γ
    M = rand(ComplexF64, 10, 10); α = rand(ComplexF64);
    y = α * M * x; Y = α * M * A
    @test (α * LinearMap(M)) * x ≈ (quat(α) * LinearMap(M)) * x ≈ y
    @test mul!(copy(y), α * LinearMap(M), x, α, false) ≈ α * M * x * α
    @test mul!(copy(y), α * LinearMap(M), x, quat(α), false) ≈ α * M * x * α
    @test mul!(copy(Y), α * LinearMap(M), A) ≈ α * M * A
    @test mul!(copy(Y), α * LinearMap(M), A, α, false) ≈ α * M * A * α
end

@testset "nonassociative number type" begin
    x = Octonion.(rand(10), rand(10), rand(10), rand(10),rand(10), rand(10), rand(10), rand(10))
    v = rand(10)
    A = Octonion.(rand(10,10), rand(10,10), rand(10,10), rand(10,10),rand(10,10), rand(10,10), rand(10,10), rand(10,10))
    α = UniformScaling(Octonion.(rand(8)...))
    β = UniformScaling(Octonion.(rand(8)...))
    L = LinearMap(A)
    @test Array(L) == A
    @test Array(α * L) == α * A
    @test Array(L * α) == A * α
    @test Array(α * L) == α * A
    @test Array(L * α) == A * α
    @test (α * L')' * v ≈ (α * A')' * v
end


text:
# using Ipopt
include(joinpath(pwd(),"src/Ipopt.jl"))
using Test, ForwardDiff, LinearAlgebra

nc = 1
nf = 2
nq = 3
nu = 2
nβ = nc*nf
ns = 2

nx = nq+nu+nc+nβ+nc+2
np = nq+nβ+4nc

idx_s = nq+nu+nc+nβ+nc .+ (1:ns)

dt = 0.1

M(q) = 1.0*Matrix(I,nq,nq)
B(q) = [1. 0. 0.;0. 1. 0.]
P(q) = [1. 0. 0.;0. 1. 0.]

G(q) = [0; 0; 9.8]

N(q) = [0; 0; 1]

qpp = [0.,0.,0.1]
v0 = [10.,-7.0, 0.]
v1 = v0 - G(qpp)*dt
qp = qpp + 0.5*dt*(v0 + v1)

v2 = v1 - G(qp)*dt
q1 = qp + 0.5*dt*(v1 + v2)

qf = [0.; 0.; 0.]
uf = [0.; 0.]

W = 1.0*Matrix(I,nq,nq)
w = -W*qf
R = 1.0e-1*Matrix(I,nu,nu)
r = -R*uf
obj_c = 0.5*qf'*W*qf + 0.5*uf'*R*uf

function unpack(x)
    q = x[1:nq]
    u = x[nq .+ (1:nu)]
    y = x[nq+nu+nc]
    β = x[nq+nu+nc .+ (1:nβ)]
    ψ = x[nq+nu+nc+nβ+nc]
	s = x[nq+nu+nc+nβ+nc .+ (1:ns)]

    return q,u,y,β,ψ,s
end

function obj(x)
    q,u,y,β,ψ,s = unpack(x)
    return 0.5*q'*W*q + w'*q + 0.5*u'*R*u + r'*u + obj_c + λ_global'*s + 0.5*mu_prev*s'*s
end

function con!(c,x)
    q,u,y,β,ψ,s = unpack(x)

    c[1] = (N(q)'*q)
    c[2] = (((0.5*y)^2 - β'*β))
    c[3:5] = (M(q)*(2*qp - qpp - q)/dt - G(q)*dt + B(q)'*u + P(q)'*β + N(q)*y)
    c[6:7] = (P(q)*(q-qp)/dt + 2.0*β*ψ)
    c[8] = s[1] - y*(N(q)'*q)
    c[9] = s[2] - ψ*((0.5*y)^2 - β'*β)

	return c
end

n = nx
m = np

x_L = zeros(nx)
x_L[1:(nq+nu)] .= -1.0e8
x_L[nq+nu+nc .+ (1:nβ)] .= -1.0e8
x_U = 1.0e8*ones(nx)
cI_idx = zeros(Bool,m)
cI_idx[1:nc+nc] .= 1

g_L = zeros(m)
g_U = zeros(m)
g_U[1:2nc] .= 1.0e8

q0 = q1
u0 = 1.0e-8*ones(nu)
y0 = 1.0e-8*ones(1)[1]
β0 = 1.0e-8*ones(nβ)
ψ0 = 1.0e-8*ones(1)[1]
s0 = 1.0e-8*ones(2)
x0 = [q0;u0;y0;β0;ψ0;s0]

function row_col!(row,col,r,c)
    for cc in c
        for rr in r
            push!(row,convert(Int,rr))
            push!(col,convert(Int,cc))
        end
    end
    return row, col
end

function eval_f(x::Vector{Float64})
  x_global .= x
  return obj(x)
end

function eval_g(x::Vector{Float64}, g::Vector{Float64})
	con!(g,x)

	return g
end

function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
	grad_f .= ForwardDiff.gradient(obj,x)
end

function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
  if mode == :Structure
		row = []
		col = []
    row_col!(row,col,(1:m),(1:n))
		rows .= vec(row)
		cols .= vec(col)
  else

  	values .= vec(ForwardDiff.jacobian(con!,zeros(m),x))
  end
end

function convergence(x,tol)
	q,u,y,β,ψ,s = unpack(x)

	cond1 = norm(s,Inf) < tol
    c_tmp = zeros(m)
	con!(c_tmp,x_global)
	c_tmp[1:2nc] .= min.(0.0,c_tmp[1:2:nc])
	cond2 = norm(c_tmp,Inf) < tol

	return cond1 && cond2
end

# This tests callbacks.
function intermediate(alg_mod::Int, iter_count::Int,
  obj_value::Float64, inf_pr::Float64, inf_du::Float64, mu::Float64,
  d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64,
  ls_trials::Int)

	println("CALLBACK")
  println("   mu_prev: $(mu_prev)")
  println("   mu: $(mu)")
  println("   x: $(x_global)")

  if mu != mu_prev
    println("-mu difference-")
    println("-AL update-")
    λ_global .+= 1/mu_prev*x_global[idx_s]
    global mu_prev = mu
  end

  return true#!convergence(x_global,1.0e-4) #iter_count < 1  # Interrupts after one iteration.
end

function eval_h(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
  if mode == :Structure

  else

  end
end

prob = createProblem(n, x_L, x_U, m, g_L, g_U, m*n, n*n,
                     eval_f, eval_g, eval_grad_f, eval_jac_g, )
addOption(prob,"hessian_approximation","limited-memory")
addOption(prob,"mu_strategy","monotone")

addOption(prob,"constr_viol_tol",1.0e-3)
addOption(prob,"tol",1.0e-3)
# addOption(prob,"")
setIntermediateCallback(prob, intermediate)

global x_global = zeros(n)
global λ_global = zeros(ns)
global mu_prev = 0.1

prob.x .= 1.0e-5*rand(n)
prob.x[1:nq] = copy(q1)
x_global .= prob.x
@time solvestat = solveProblem(prob)

λ_global
mu_prev
# @test Ipopt.ApplicationReturnStatus[solvestat] == :User_Requested_Stop

prob.x[idx_s]


text:
"""
Minimal example for WhithamGreenNaghdi
"""
using WaterWaves1D

"""
	solveandplot()

Solves an initial value problem with WhithamGreenNaghdi model.
See parameters below.

"""
function solveandplot()
	# Set up
	param = ( μ  = 1, ϵ  = 1,
				N  = 2^10,   # number of collocation points
	            L  = 10*π,   # mesh of size 2L
	            T  = 1,      # final time
	            dt = 1/10^4, # timestep
				ns=1000      # stores data every ns time steps
				)

	mesh=Mesh(param); # construct mesh of collocation points, Fourier modes, etc.
	η= exp.(-(mesh.x).^2);v= 0*η; # Initial data

	init     = Init(mesh,η,v)
	model = WhithamGreenNaghdi(param;
				iterate=false) # iterate = true for GMRES iterate solver
	problem = Problem(model, init, param)

	# Compute
	@time solve!( problem )

	# Plot
	fftηfin=last(problem.data.U)[:,1]
	plt = plot(layout=(1,2))
	plot!(plt[1,1],fftshift(mesh.k),fftshift(log10.(abs.(fftηfin)));
			title = "frequencies")
	plot!(plt[1,2],mesh.x,real.(ifft(fftηfin));
			title = "physical space")
	display(plt)
end

solveandplot()


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

doc"""
    struct ExecContext
        use_threads::Bool
        onprocs::Vector{Int64}
    end

Functions that take an `ExecContext` argument must limit their use of
threads and processes accordingly. Depending on `use_threads`, the function
may use all (or only a single) thread(s) on each process in `onprocs` (in
addition to the current thread on the current process).

The caller may choose to change the `ExecContext` from call to call, based on
execution time and latency measurements, etc.

Functions can announce their [`BAT.ExecCapabilities`](@ref) via [`exec_capabilities`](@ref).
"""
struct ExecContext
    use_threads::Bool
    onprocs::Vector{Int64}
end

export ExecContext

ExecContext() = ExecContext(true, workers())

doc"""
    struct ExecCapabilities
        nthreads::Int
        threadsafe::Bool
        nprocs::Int
        remotesafe::Bool
    end

Specifies the execution capabilities of functions that support an
`ExecContext` argument.

`nthreads` specifies the maximum number of threads the function can utilize
efficiently, internally. If `nthreads <= 1`, the function implementation is
single-threaded. `nthreads == 0` indicates that the function is cheap and
that when used in combination with other functions, their capabilities should
dominate.

`threadsafe` specifies whether the function is thread-safe, and can be can be
run on multiple threads in parallel by the caller.

`nprocs` specifies the maximum number of worker processes the function can
utilize efficiently, internally. If `procs <= 1`, the function cannot
use worker processes. `nthreads == 0` carries equivalent meaning to
`nthreads == 0`.

`remotesafe` specifies that the function can be run on a remote thread,
it implies that the function arguments can be (de-)serialized safely.

Functions with an `ExecContext` argument should announce their capabilities
via methods of `exec_capabilities`. Functions should, ideally, either support
internal multithreading (`nthreads > 1`) or be thread-safe
(`threadsafe` == true). Likewise, functions should either utilize worker
processes (`nprocs` > 1) internally or support remote execution
(`remotesafe` == true) by the caller.
"""
struct ExecCapabilities
    nthreads::Int
    threadsafe::Bool
    nprocs::Int
    remotesafe::Bool
end

export ExecCapabilities

doc"""
    intersect(a:ExecCapabilities, b:ExecCapabilities)

Get the intersection of execution capabilities of a and b, i.e. the
`ExecCapabilities` that should be used when to functions are used in
combination (e.g. in sequence).
"""
Base.intersect(a::ExecCapabilities, b::ExecCapabilities) = ExecCapabilities(
    if a.nthreads == 0
        b.nthreads
    elseif b.nthreads == 0
        a.nthreads
    else
        # ToDo/Decision: Better use maximum?
        min(a.nthreads, b.nthreads)
    end,
    a.threadsafe && b.threadsafe,
    if a.nprocs == 0
        b.nprocs
    elseif b.nprocs == 0
        a.nprocs
    else
        # ToDo/Decision: Better use maximum?
        min(a.nprocs, b.nprocs)
    end,
    a.remotesafe && a.remotesafe
)

doc"""
    exec_capabilities(f, args...)::ExecCapabilities

Determines the execution capabilities of a function `f` that supports an
`ExecContext` argument, when called with arguments `args...`. The
`ExecContext` argument itself is excluded from `args...` for
`exec_capabilities`.

Before calling `f`, the caller must use

    exec_capabilities(f, args...)

to determine the execution capabilities of `f` with the intended arguments,
and take the resulting `ExecCapabilities` into account. If `f` is not
thread-safe (but remote-safe), and the caller needs to run it on multiple
threads, the caller may deep-copy the function arguments.
"""
function exec_capabilities end

function negotiate_exec_context(context::ExecContext, target_caps::AbstractVector{ExecCapabilities})
    target_caps_threadsafe = all(x -> x.threadsafe, target_caps)
    target_caps_nthreads = minimum(x -> x.nthreads, target_caps)

    if context.use_threads
        target_use_threads = target_caps_nthreads > 1
        self_use_threads = target_caps_threadsafe && !target_use_threads
    else
        target_use_threads = false
        self_use_threads = false
    end

    target_caps_remotesafe = all(x -> x.remotesafe, target_caps)

    # ToDo: Take different nprocs into account
    target_caps_nprocs = minimum(x -> x.nprocs, target_caps)

    nprocs_avail = length(context.onprocs)
    if isempty(context.onprocs)
        target_onprocs = [myid()]
        self_onprocs = [myid()]
    elseif nprocs_avail == 1 && context.onprocs[1] == myid()
        target_onprocs = context.onprocs
        self_onprocs = context.onprocs
    else
        if target_caps_remotesafe
            # ToDo: Distribute available processes for targets
            error("Not implemented yet")
            target_onprocs = error("Not implemented yet")
            self_onprocs = error("Not implemented yet")
        else
            # ToDo: Distribute available processes over targets
            target_onprocs = error("Not implemented yet")
            self_onprocs = [myid()]
        end
    end

    self_context = ExecContext(self_use_threads, target_onprocs)
    # ToDo: Improve naive implementation:
    target_context = map(x -> ExecContext(target_use_threads, self_onprocs), target_caps)

    (self_context, target_context)
end


text:
export Abstract_Model2Fit_Peak

"""
A specialization of [`Abstract_Model2Fit`](@ref) dedicated to peak types.

TODO: in the future some extra methods are to be defined like
- peak_position
- peak_area...
"""
abstract type Abstract_Model2Fit_Peak <: Abstract_Model2Fit end


text:
# that's still a very lame function, but might be useful once more sampling
# algorithms are implemented (or not...)
function sample(model::Model, algorithm)
    run(algorithm, model)
end


text:
module reflectometryNN
    using ForwardDiff
    using Flux
    using BSON: @save
    include("groupDelay.jl")
    mutable struct trainingData
        data
        density
        freqs
        measurement_max
        normalizationRad
        XMode
        B0
        R0
    end

    struct NNOutput
        NN
        normalizationRad::Float64
    end

"""
Function to create neural network training data.

Inputs:
        freq                         -- Vector of frequencies of the reflectometry measurements (GHz)
        dif_fits::Int64              -- The number of different profiles in a single batch (batch size)
        coeffs_fit::Array{Float64}   -- The 6 fit coefficients for the 5th order polynomial for plasma density in units of 1e19 (per cubic meter)
                                            n(r) = 1e19*(coeffs_fit[1]*r^5+coeffs_fit[2]*r^4+coeffs_fit[3]*r^3+coeffs_fit[4]*r^2+coeffs_fit[5]*r^1+coeffs_fit[6])
                                            where r=radius into the plasma from the plasma facing wall.
        normalizationRad::Number     -- The output radius normalization factor, that helps with training. Make output order 1-10.
        calibration::Array{Float64}  -- Vector of calibrated group delay length (non-plasma length) at each frequency, for a given reflectometry measurement. Length in METERS.

Optional inputs:
        trainingSpread::Float64      -- The variation in the profiles, fraction of each coefficient. Default value is 0.5.
        XMode::Bool                  -- A boolean if XMode is the desired profile (true or false). Default is false (OMode).
        Xcutoff::String              -- The type of XMode cutoff ("LEFT" or "RIGHT").
        Bmag::Number                 -- The magnetic field magnitude measured at major-radius R0 (Tesla)
        R0::Number                   -- The major-radial location of the magnetic field measurement (meters)

    e.g. for the O-Mode HFS reflectometer on DIII-D
    makeTrainingData(range(6, 27, length = 200), 1000, ,[2.181e5, -2.901e4, 1252, -16.14, 2.381, 0.03], 5e7; trainingSpread = 0.5)

Output is a trainingData struct that includes:
        trainingData.data             -- the simulated training profiles
        trainingData.freqs            -- the frequency vector used for the measurements (GHz)
        trainingData.measurement_max  -- the normalization used for each group delay in the batch
        trainingData.normalizationRad -- the output radius normalization inputted into the function
        trainingData.Xmode            -- Boolean that says if XMode was enabled
        trainingData.B0               -- Magnitude of magnetic field for XMode (Tesla)
        trainingData.R0               -- Major-radial location of magnetic field measurement for XMode (meters)
"""
    function makeTrainingData(freq,dif_fits::Int64, coeffs_fit::Array{Float64}, normalizationRad::Number, calibration::Array{Float64}; trainingSpread::Float64 = 0.5, XMode::Bool = false, Xcutoff::String="LEFT", Bmag::Number=0, R0::Number=0)
        if XMode
            data = (Array{Float64}(undef,length(freq)+1,dif_fits),Array{Float64}(undef,length(freq),dif_fits))
        else
            data = (Array{Float64}(undef,length(freq),dif_fits),Array{Float64}(undef,length(freq),dif_fits))
        end
        measurement_max = Array{Float64}(undef,(dif_fits))
        radius = Array{Float64}(undef,length(freq))
        density = Array{Float64}(undef,length(freq))
        for j in 1:1:dif_fits
            monocheck=1
            local coeffs
            monocheck=true
            while monocheck
                coeffs = [coeffs_fit[1]+trainingSpread*rand()*(-1)^rand(1:2), coeffs_fit[2]+trainingSpread*rand()*(-1)^rand(1:2), coeffs_fit[3]+trainingSpread*rand()*(-1)^rand(1:2),
                                coeffs_fit[4]+trainingSpread*rand()*(-1)^rand(1:2), coeffs_fit[5]+trainingSpread*rand()*(-1)^rand(1:2), coeffs_fit[6]+trainingSpread*rand()*(-1)^rand(1:2)]

                    density(x) = 1e19*(coeffs[1]*(x)^5+coeffs[2]*(x)^4+coeffs[3]*(x)^3+coeffs[4]*(x)^2+coeffs[5]*(x)+coeffs[6])
                    xtest=0
                    while  density(xtest) < (freq[end]*2*pi*1e9)^2*8.854e-12*9.109e-31/(1.6022e-19)^2
                        xtest += 0.0001
                        if ForwardDiff.derivative(density,xtest) < 0
                            monocheck=true
                            break
                        end
                    end
                    if density(xtest) > (freq[end]*2*pi*1e9)^2*8.854e-12*9.109e-31/(1.6022e-19)^2
                        monocheck=false
                    end
            end
            if XMode
                if Xcutoff=="LEFT"
                    radius,density = radCalc.XModeL(freq,coeffs,Bmag,R0)
                elseif Xcutoff == "RIGHT"
                    radius,density = radCalc.XModeR(freq,coeffs,Bmag,R0)
                end
            else
                radius,density = radCalc.OMode(freq,coeffs)
            end
            if XMode
                measurement = dphidw_xmode(freq,coeffs,calibration,Bmag,R0,radius)
                data[1][1,j] = Bmag
                data[1][2:end,j] = measurement/maximum(measurement)
                measurement_max[j] = maximum(measurement)
                data[2][:,j] = radius/maximum(measurement)/normalizationRad
            else
                measurement = dphidw_omode(freq,coeffs,calibration,radius)
                data[1][:,j] = measurement/maximum(measurement)
                measurement_max[j] = maximum(measurement)
                data[2][:,j] = radius/maximum(measurement)/normalizationRad
            end
        end

        if XMode
            data_formated = Flux.Data.DataLoader(data;batchsize=dif_fits,shuffle=false,partial=true)
            output = trainingData(data_formated, density, freq, measurement_max, normalizationRad, XMode, Bmag, R0)
        else
            data_formated = Flux.Data.DataLoader(data;batchsize=dif_fits,shuffle=false,partial=true)
            output = trainingData(data_formated, density, freq, measurement_max, normalizationRad, XMode, NaN, NaN)
        end

        output
    end

"""
Function to create the neural network object

Inputs:
        dataNN:trainingData  -- a trainingData struct from the makeTrainingData() function
        epochs_num::Int64    -- the number of epochs you wish to train the neural network for

Optional inputs:
        actication           -- the activation function used in the neural network. Default is sigmoid.
        optimizer            -- the optimizer (descent algorithm) used in the neural network. Default is NADAM
        learningrate         -- the initial optimizer learning rate. Default is 1e-5
        learningdecay        -- the exponential factor of decay of the learning rate. Decay is 1e-4 per epoch.
        neurons              -- the number of neurons. Default is 200.
        layers               -- the number of layers in the neural network (2 or 3). Default is 2.

Output is a NNOutput struct, with form:
        NNOutput.NN               -- the neural network object.
        NNoutput.normalizationRad -- the output radius normalization used for this neural network object

"""
    function NeuralNet(dataNN::trainingData, epochs_num::Int64; activation = swish, optimizer = RMSProp, learningrate::Float64 = 1e-4, learningdecay::Float64 = 0.001, neurons::Int64=200, layers::Int64=2, min_learningrate::Float64=1e-6,print_epoch::Number = 1000)
        local loss_check = 1
        local loss_count = 0
        if layers==2
            if dataNN.XMode
                NN = Chain(Dense(length(dataNN.freqs)+1,neurons, activation ),
                        Dense(neurons,length(dataNN.freqs)))
            else
                NN = Chain(Dense(length(dataNN.freqs),neurons, activation ),
                        Dense(neurons,length(dataNN.freqs)))
            end

        else
            if dataNN.XMode
                NN = Chain(Dense(length(dataNN.freqs)+1,neurons, activation ),
                        Dense(neurons, neurons, activation),
                        Dense(neurons,length(dataNN.freqs)))
            else
                NN = Chain(Dense(length(dataNN.freqs),neurons, activation ),
                        Dense(neurons, neurons, activation),
                        Dense(neurons,length(dataNN.freqs)))
            end
        end
        function losss(x,y)
            loss_check = sum(abs2,NN(x)-(y))
        end
        for j in 1:1:epochs_num
                if j % print_epoch == 0
                    loss_count += 1
                    print(string(loss_count), string(" / "), string(round(epochs_num/print_epoch)), string(":  "), string(loss_check), string("\n"))
                end
                descent_param = max(min_learningrate,learningrate*exp(-1*j*learningdecay))
                Flux.train!(losss,Flux.params(NN),dataNN.data,optimizer(descent_param))
        end
       NNOutput(NN, dataNN.normalizationRad)
    end

"""
A function that simply saves an object (typically used for the NNOutput object to be used later).

Inputs:
        filename::String  -- the filename of the object, saved in current directory
        object            -- the object to be saved, typically the NNOutput object for this use case
"""
    function save(filename::String, object)
        @save filename object
    end

end


text:
# Verify that we handle correctly the case where the starting point is a zero
@testset "already converged" begin
function f!(fvec, x)
    fvec[1] = (x[1]+3)*(x[2]^3-7)+18
    fvec[2] = sin(x[2]*exp(x[1])-1)
end

df = OnceDifferentiable(f!, rand(2), rand(2))

r = nlsolve(df, [ 0.; 1.], method = :trust_region)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-8
@test r.iterations == 0

r = nlsolve(df, [ 0.; 1.], method = :newton)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-8
@test r.iterations == 0
end


text:
"""
    Layer

Defines a new layer within the video.

# Fields
- `frames::Frames`: A range of frames for which the `Layer` exists
- `width::Int`: Width of the layer
- `height::Int`: hegiht of the layer
- `position::Point`: initial positon of the center of the layer on the main canvas
- `layer_objects::Vector{AbstractObject}`: Objects defined under the layer
- `actions::Vector{AbstractAction}`: a list of actions applied to the entire layer
- `current_setting::LayerSetting`: The current state of the layer see [`LayerSetting`](@ref)
- `opts::Dict{Symbol,Any}`: can hold any options defined by the user
- `image_matrix::Vector`: Hold the Drwaing of the layer as a Luxor image matrix
"""
mutable struct Layer <: AbstractObject
    frames::Frames
    width::Int
    height::Int
    position::Point
    layer_objects::Vector{AbstractObject}
    actions::Vector{AbstractAction}
    current_setting::LayerSetting
    opts::Dict{Symbol,Any}
    image_matrix::Union{Array,Nothing}
    layer_cache::LayerCache
end

"""
    CURRENT_LAYER

holds the current layer in an array to be declared as a constant
The current layer can be accessed using CURRENT_LAYER[1]
"""
const CURRENT_LAYER = Array{Layer,1}()

# for width, height and position defaults are defined in the to_layer_m function
function Layer(
    frames,
    width,
    height,
    position::Point;
    layer_objects::Vector{AbstractObject} = AbstractObject[],
    actions::Vector{AbstractAction} = AbstractAction[],
    setting::LayerSetting = LayerSetting(),
    misc::Dict{Symbol,Any} = Dict{Symbol,Any}(),
    mat = nothing,
    layer_cache::LayerCache = LayerCache(),
)
    if width === nothing
        width = CURRENT_VIDEO[1].width
    end
    if height === nothing
        height = CURRENT_VIDEO[1].height
    end

    layer = Layer(
        frames,
        width,
        height,
        position,
        layer_objects,
        actions,
        setting,
        misc,
        mat,
        layer_cache,
    )

    if isempty(CURRENT_LAYER)
        push!(CURRENT_LAYER, layer)
    else
        CURRENT_LAYER[1] = layer
    end
    push!(CURRENT_VIDEO[1].layers, layer)

    return layer
end


text:
module VTK

using ...TopOptProblems: TopOptProblems, StiffnessTopOptProblem, Ferrite
using WriteVTK

export save_mesh

function save_mesh(filename, problem::StiffnessTopOptProblem)
    topology = ones(getncells(TopOptProblems.getdh(problem).grid))
    vtkfile = WriteVTK.vtk_grid(filename, problem, topology)
    return outfiles = WriteVTK.vtk_save(vtkfile)
end
function save_mesh(filename, problem, solver)
    return save_mesh(filename, problem, solver.vars)
end
function save_mesh(filename, alg)
    problem = alg.obj.problem
    vars = alg.optimizer.obj.solver.vars
    return save_mesh(filename, problem, vars)
end
function save_mesh(filename, problem, vars::AbstractVector)
    vtkfile = WriteVTK.vtk_grid(filename, problem, vars)
    return outfiles = WriteVTK.vtk_save(vtkfile)
end
function WriteVTK.vtk_grid(
    filename::AbstractString,
    problem::StiffnessTopOptProblem{dim,T},
    vars::AbstractVector{T},
) where {dim,T}
    varind = problem.varind
    black = problem.black
    white = problem.white
    grid = problem.ch.dh.grid
    full_top = length(vars) == length(TopOptProblems.getdh(problem).grid.cells)

    celltype = Ferrite.cell_to_vtkcell(Ferrite.getcelltype(grid))
    cls = Ferrite.MeshCell[]
    for (i, cell) in enumerate(Ferrite.CellIterator(grid))
        if full_top
            if vars[i] >= 0.5
                push!(cls, Ferrite.MeshCell(celltype, copy(Ferrite.getnodes(cell))))
            end
        else
            if black[i]
                push!(cls, Ferrite.MeshCell(celltype, copy(Ferrite.getnodes(cell))))
            elseif !white[i]
                if vars[varind[i]] >= 0.5
                    push!(cls, Ferrite.MeshCell(celltype, copy(Ferrite.getnodes(cell))))
                end
            end
        end
    end
    coords = reshape(reinterpret(T, Ferrite.getnodes(grid)), (dim, Ferrite.getnnodes(grid)))
    return vtk_grid(filename, coords, cls)
end

end


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Random
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "statevector" begin

    @testset "test statevector int constructor" begin
        N = rand(1:10)
        s = Statevector(N)
        @test s.N == N
        @test s.state == zeros(ComplexF64, 2^N)
    end

    @testset "test statevector vector constructor" begin
        N = rand(1:10)
        ψ = randn(ComplexF64, 2^N)
        s = Statevector(ψ)
        @test s.N == N
        @test s.state == ψ
    end

    @testset "test statevector indexing" begin
        N = rand(1:10)
        ψ = randn(ComplexF64, 2^N)
        s = Statevector(ψ)
        r = rand(1:2^N)
        @test s[r] == ψ[r]
    end

    @testset "test statevector attributes" begin
        N = rand(1:10)
        ψ = randn(ComplexF64, 2^N)
        s = Statevector(ψ)
        @test size(s) == size(ψ)
        @test length(s) == length(ψ)
        @test IndexStyle(typeof(s)) == IndexLinear()
    end
end


text:
"""
    Second experimental scenario.

    In this experimental scenario, we fixed each node threshold
    to a random value between 1 and its degree, varying it at each run of the experiment.
    We used a majority policy for the hyperedge thresholds, instead.
"""

using Pkg, Distributed
#addprocs(4)
Pkg.activate(".")
@everywhere using Distributed, Pkg
@everywhere Pkg.activate(".")
using LTMSim, SimpleHypergraphs, Random, Serialization
@everywhere using LTMSim, SimpleHypergraphs, Random, Serialization

#
#
#
project_path = dirname(pathof(LTMSim))
data_path = joinpath(project_path, "..", "data", "hgs")
res_path = joinpath(project_path, "..", "res", "journal", "randV_propE05_w_wo_opt.data")

hg_files = readdir(data_path)
hgs = [prune_hypergraph!(hg_load(joinpath(data_path, hg_file))) for hg_file in hg_files]

runs = 50
data = Dict{String, Vector{Vector{Int}}}()

push!(data, "BinarySearch(H)" => Vector{Int}[])
push!(data, "Greedy(H)" => Vector{Int}[])
push!(data, "Greedy([H]₂)" => Vector{Int}[])
push!(data, "SubTSS(H)" => Vector{Int}[])

push!(data, "BinarySearch(H)-noOpt" => Vector{Int}[])
push!(data, "Greedy(H)-noOpt" => Vector{Int}[])
push!(data, "Greedy([H]₂)-noOpt" => Vector{Int}[])
push!(data, "SubTSS(H)-noOpt" => Vector{Int}[])

for index in 1:length(hgs)
    println("Index=$index, $(hg_files[index])")

    h = hgs[index]

    results = @distributed (append!) for run=1:runs
		println("run=$run at proc $(myid())")
        Random.seed!(run)

        metaV = randMetaV(h)
        metaE = proportionalMetaE(h,0.5)

        # with optimization
        r1 = greedy_tss_2section(h,metaV,metaE; opt=true)
        r2 = bisect(h,metaV,metaE; opt=true)
        r3 = greedy_tss(h,metaV,metaE; opt=true)
        r4 = sub_tss(h,metaV,metaE; opt=true)

        # without optimization
        r5 = greedy_tss_2section(h,metaV,metaE; opt=false)
        r6 = bisect(h,metaV,metaE; opt=false)
        r7 = greedy_tss(h,metaV,metaE; opt=false)
        r8 = sub_tss(h,metaV,metaE; opt=false)

        [(r1,r2,r3,r4,r5,r6,r7,r8)]
    end

    push!(data["Greedy([H]₂)"], [r[1] for r in results])
    push!(data["BinarySearch(H)"], [r[2] for r in results])
    push!(data["Greedy(H)"], [r[3] for r in results])
    push!(data["SubTSS(H)"], [r[4][1] for r in results])

    push!(data["Greedy([H]₂)-noOpt"], [r[5] for r in results])
    push!(data["BinarySearch(H)-noOpt"], [r[6] for r in results])
    push!(data["Greedy(H)-noOpt"], [r[7] for r in results])
    push!(data["SubTSS(H)-noOpt"], [r[8][1] for r in results])
end

serialize(res_path, data)