text:
struct TeXDocument
    contents::String
end

"""
    TeXDocument(contents::AbstractString, add_defaults::Bool; requires, preamble, class, classoptions)

This constructor function creates a TeX document which can be passed to `teximg`.
All arguments are to be passed as strings.

If `add_defaults` is `false`, then we will *not* automatically add document structure.
Note that in this case, keyword arguments will be disregarded and `contents` must be
a complete LaTeX document.

Available keyword arguments are:
- `requires`: code which comes before `documentclass` in the preamble.  Default: `raw"\\RequirePackage{luatex85}"`.
- `class`: the document class.  Default (and what you should use): `"standalone"`.
- `classoptions`: the options you should pass to the class, i.e., `\\documentclass[\$classoptions]{\$class}`.  Default: `"preview, tightpage, 12pt"`.
- `preamble`: arbitrary code for the preamble (between `\\documentclass` and `\\begin{document}`).  Default: `raw"\\usepackage{amsmath, xcolor} \\pagestyle{empty}"`.
"""
function TeXDocument(
            contents::AbstractString,
            add_defaults::Bool;
            requires::AbstractString = raw"\RequirePackage{luatex85}",
            class::AbstractString = "standalone",
            classoptions::AbstractString = "preview, tightpage, 12pt",
            preamble::AbstractString = raw"""
                        \usepackage{amsmath, xcolor}
                        \pagestyle{empty}
                        """,
        )
        if add_defaults
            return TeXDocument(
                """
                $(requires)

                \\documentclass[$(classoptions)]{$(class)}

                $(preamble)

                \\begin{document}

                $(contents)

                \\end{document}
                """
            )
        else
            return TeXDocument(contents)
        end
end

"""
    texdoc(contents::AbstractString; kwargs...)

A shorthand for `TeXDocument(contents, add_defaults=true; kwargs...)`.

Available keyword arguments are:

- `requires`: code which comes before `documentclass` in the preamble.  Default: `raw"\\RequirePackage{luatex85}"`.
- `class`: the document class.  Default (and what you should use): `"standalone"`.
- `classoptions`: the options you should pass to the class, i.e., `\\documentclass[\$classoptions]{\$class}`.  Default: `"preview, tightpage, 12pt"`.
- `preamble`: arbitrary code for the preamble (between `\\documentclass` and `\\begin{document}`).  Default: `raw"\\usepackage{amsmath, xcolor} \\pagestyle{empty}"`.

"""
texdoc(contents; kwargs...) = TeXDocument(contents, true; kwargs...)

function Base.convert(::Type{String}, doc::TeXDocument)
    return Base.convert(String, doc.contents)
end

struct CachedTeX
    doc::TeXDocument
    handle::Rsvg.RsvgHandle
    raw_dims::Rsvg.RsvgDimensionData
    svg::String
end

"""
    CachedTeX(doc::TeXDocument, dpi = 72.0; kwargs...)

Compile a `TeXDocument`, compile it and return the cached TeX object.

A `CachedTeX` struct stores the document, the Rsvg handle, the Rsvg dimensions
and the SVG as a string for error checking.

In `kwargs`, one can pass anything which goes to the internal function `compile_latex`.
These are primarily:
- `engine = \`lualatex\`/\`xelatex\`/...`: the LaTeX engine to use when rendering
- `options=\`-file-line-error\``: the options to pass to `latexmk`.
"""
function CachedTeX(doc::TeXDocument, dpi = 72.0; method = :pdf, kwargs...)
    svg = if method == :dvi
        dvi2svg(latex2dvi(convert(String, doc); kwargs...))
    elseif method == :pdf
        pdf2svg(latex2pdf(convert(String, doc); kwargs...))
    else
        @error("$method not recognized!  Must be one of (:dvi, :pdf).")
    end

    handle = svg2rsvg(svg, dpi)

    dims = Rsvg.handle_get_dimensions(handle)
    return CachedTeX(
        doc,
        handle,
        dims,
        svg
    )
end

function CachedTeX(str::String, dpi = 72.0; kwargs...)
    return CachedTeX(implant_text(str), dpi; kwargs...)
end

function CachedTeX(x::LaTeXString, dpi = 72.0; kwargs...)
    x = convert(String, x)
    return if first(x) == "\$" && last(x) == "\$"
        CachedTeX(implant_math(x[2:end-1]), dpi; kwargs...)
    else
        CachedTeX(implant_text(x), dpi; kwargs...)
    end
end

# do not rerun the pipeline on CachedTeX
CachedTeX(ct::CachedTeX, dpi=72.0) = ct

# function new_cachedtex(doc::TeXDocument)
#     svg = compile_latex(convert(String, doc), format="pdf", read_format="svg")
#     handle = svg2rsvg(String(svg), 72.0)
#     dims = Rsvg.handle_get_dimensions(handle)
#     return CachedTeX(
#         doc, handle, dims, svg
#     )
# end

function Base.show(io::IO, ct::CachedTeX)
    if length(ct.doc.contents) > 1000
        println(io, "CachedTeX(TexDocument(...), $(ct.handle), $(ct.raw_dims))")
    else
        println(io, "CachedTeX($(ct.doc), $(ct.handle), $(ct.raw_dims))")
    end
end

function implant_math(str)
    return TeXDocument(
        """\\(\\displaystyle $str\\)""", true;
        requires = "\\RequirePackage{luatex85}",
        preamble = """
        \\usepackage{amsmath, amsfonts, xcolor}
        \\pagestyle{empty}
        """,
        class = "standalone",
        classoptions = "preview, tightpage, 12pt",
    )
end

function implant_text(str)
    return TeXDocument(
        str, true;
        requires = "\\RequirePackage{luatex85}",
        preamble = """
        \\usepackage{amsmath, amsfonts, xcolor}
        \\pagestyle{empty}
        """,
        class = "standalone",
        classoptions = "preview, tightpage, 12pt"
    )
end


text:
# included from test_sesam3, sp, sol, p, u0 need to be defined

@testset "non-negative pools" begin
    #st = first(states(sp))
    for st in states(sp)
        #@show st
        #@test all(sol[st] .>= 0.0)
        @test all(sol[st] .>= -eps(Float64)*100)
    end
    # minimum(sol[s.α_LP])
    # minimum(sol[s.α_RP])
end;

@testset "microbial dC balance" begin
    #Plots.plot(sol,  vars=[s.u_C, s.r_B+s.syn_B+s.syn_Enz])
    @test all(isapprox.(
        sol[s.u_C],
        sol[s.r_B+s.syn_B+s.syn_Enz],
        rtol = 1e-6))
end;

@testset "system dC balance" begin
    change = sol[s.dB + s.dL + s.dR]
    output = sol[s.r_tot]
    input = sol[s.i_L]
    #Plots.plot(sol.t, [input, output, change, change .+ output])
    @test all(isapprox.(input, change .+ output, rtol = 1e-6))
end;

@testset "microbial dN balance" begin
    uptake = p[s.ν_N]*sol[s.u_NOM] + sol[s.resorp_N]
    usage = sol[s.syn_Enz]/p[s.β_NEnz] + sol[s.syn_B]/p[s.β_NB] + sol[s.Φ_NB]
    #Plots.plot(sol.t, [uptake, usage])
    @test all(isapprox.(uptake, usage, rtol = 1e-6))
end;

@testset "system dN balance" begin
    change = sol[s.dB]/p[s.β_NB] + sol[s.dL_N + s.dR_N + s.dI_N]
    output = sol[s.u_PlantN + s.leach_N]
    input = sol[s.i_L/s.β_Ni + s.i_IN]
    #sol[s.leach_N]
    #sol[s.i_IN]
    bo = 1:10
    bo = 1:length(input)
    #Plots.plot(sol.t[bo], [input[bo], output[bo], change[bo], (change .+ output)[bo]], label = ["input" "output" "change" "change+output"])
    @test all(isapprox.(input, change .+ output, rtol = 1e-6))
end;

@testset "balanced allocation sums to one" begin
    #Plots.plot(sol, vars=[s.α_LT,s.α_RT, s.α_LT+s.α_RT])
    @test all(isapprox.(sol[s.α_LT + s.α_RT + s.α_LPT + s.α_RPT], 1.0, rtol = 1e-8))
    @test all(isapprox.(sol[s.α_L + s.α_R + s.α_LP + s.α_RP], 1.0, rtol = 1e-8))
end;

@testset "microbial dP balance" begin
    uptake = p[s.ν_P]*sol[s.u_POM] + sol[s.resorp_P];
    usage = sol[s.syn_Enz]/p[s.β_PEnz] + sol[s.syn_B]/p[s.β_PB] + sol[s.Φ_PB];
    #Plots.plot(sol.t, [uptake, usage])
    @test all(isapprox.(uptake, usage, rtol = 1e-6))
    @test all(sol[s.resorp_P] .>= 0.0)
end;

@testset "system dP balance" begin
    change = sol[s.dB]/p[s.β_PB] + sol[s.dL_P + s.dR_P + s.dI_P]
    output = sol[s.u_PlantP + s.leach_P]
    input = sol[s.i_L/s.β_Pi + s.i_IP]
    #sol[s.leach_N]
    #sol[s.i_IN]
    bo = 1:10
    bo = 1:length(input)
    bo = length(input) .- (10:-1:0)
    #Plots.plot(sol.t[bo], [input[bo], output[bo], change[bo], (change .+ output)[bo]], label = ["input" "output" "change" "change+output"])
    #Plots.plot(sol.t[bo], [input[bo], (change .+ output)[bo]], label = ["input" "change+output"])
    @test all(isapprox.(input, change .+ output, rtol = 1e-6))
end;

ext: jl


text:
# Command: Get Window Handle
"""
	window_handle(session::Session)::String

The window handle associated with the current top-level browsing context.
"""
function window_handle(session::Session)::String
    @unpack addr, id = session
    response = HTTP.get(
        "$addr/session/$id/window_handle",
        [("Content-Type" => "application/json")],
    )
    @assert response.status == 200
    convert(String, JSON3.read(response.body).value)
end


text:
"""
    abstract type AbstractClangType <: Any
Supertype for `clang::Type`s.

1. `clang::Type` has no corresponding handle type in Julia.
2. `CXType_` is an opaque pointer for `clang::Type *`.
3. `get_type_ptr`/`getTypePtr` is for converting a `QualType` to a `CXType_`.
4. `get_qual_type`/`getCanonicalTypeInternal` is for converting a `CXType_` to a `QualType`.
"""
abstract type AbstractClangType end

struct UnexposedType{T<:AbstractClangType} <: AbstractClangType
    ty::T
end

"""
    abstract type AbstractQualType <: AbstractClangType
Supertype for `clang::QualType`s.

As we use `QualType` everywhere, it's useful to create a type hierarchy of different
`QualType`s on the Julia side for multiple dispatch.
"""
abstract type AbstractQualType <: AbstractClangType end

Base.unsafe_convert(::Type{CXQualType}, x::T) where {T<:AbstractQualType} = x.ptr
Base.cconvert(::Type{CXQualType}, x::T) where {T<:AbstractQualType} = x

"""
    QualType <: AbstractQualType
Represent a qualified type.

Note that, the underlying pointer is NOT a *pointer* to a `clang::QualType` object. Instead,
it's the opaque pointer representation of the `clang::QualType` itself.
"""
struct QualType <: AbstractQualType
    ptr::CXQualType
end
QualType(x::T) where {T<:AbstractQualType} = QualType(x.ptr)

# FIXME: find a use case
# """
#     CanQualType <: AbstractClangType
# Represent a canonical, qualified type.
# """
# struct CanQualType <: AbstractClangType
#     ty::QualType
# end
# CanQualType(x::AbstractQualType) = CanQualType(QualType(x))

"""
    abstract type AbstractBuiltinType <: AbstractQualType
Supertype for Clang builtin types.
"""
abstract type AbstractBuiltinType <: AbstractQualType end

"""
    BuiltinType <: AbstractBuiltinType
A builtin `QualType`.

Note that, we treat `BuiltinType` as a `QualType`, not a `CXType_`.
"""
struct BuiltinType <: AbstractBuiltinType
    ptr::CXQualType
end

struct VoidTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct BoolTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct CharTy <: AbstractBuiltinType
    ptr::CXQualType
end

# [C++ 3.9.1p5].
struct WCharTy <: AbstractBuiltinType
    ptr::CXQualType
end

# Same as WCharTy in C++, integer type in C99.
struct WideCharTy <: AbstractBuiltinType
    ptr::CXQualType
end

# [C99 7.24.1], integer type unchanged by default promotions.
struct WIntTy <: AbstractBuiltinType
    ptr::CXQualType
end

# [C++20 proposal]
struct Char8Ty <: AbstractBuiltinType
    ptr::CXQualType
end

# [C++0x 3.9.1p5], integer type in C99.
struct Char16Ty <: AbstractBuiltinType
    ptr::CXQualType
end

# [C++0x 3.9.1p5], integer type in C99.
struct Char32Ty <: AbstractBuiltinType
    ptr::CXQualType
end

struct SignedCharTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct ShortTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct IntTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct LongTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct LongLongTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct Int128Ty <: AbstractBuiltinType
    ptr::CXQualType
end

struct UnsignedCharTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct UnsignedShortTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct UnsignedIntTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct UnsignedLongTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct UnsignedLongLongTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct UnsignedInt128Ty <: AbstractBuiltinType
    ptr::CXQualType
end

struct FloatTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct DoubleTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct LongDoubleTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct Float128Ty <: AbstractBuiltinType
    ptr::CXQualType
end

# [OpenCL 6.1.1.1], ARM NEON
struct HalfTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct BFloat16Ty <: AbstractBuiltinType
    ptr::CXQualType
end

# C11 extension ISO/IEC TS 18661-3
struct Float16Ty <: AbstractBuiltinType
    ptr::CXQualType
end

struct FloatComplexTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct DoubleComplexTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct LongDoubleComplexTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct Float128ComplexTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct VoidPtrTy <: AbstractBuiltinType
    ptr::CXQualType
end

struct NullPtrTy <: AbstractBuiltinType
    ptr::CXQualType
end

"""
    abstract type AbstractComplexType <: AbstractQualType
Supertype for `ComplexType`s.
"""
abstract type AbstractComplexType <: AbstractQualType end

"""
    struct ComplexType <: AbstractComplexType
A `ComplexType`.

Note that, we treat `ComplexType` as a `QualType`, not a `CXType_`.
"""
struct ComplexType <: AbstractComplexType
    ptr::CXQualType
end

"""
    struct PointerType <: AbstractQualType
A `PointerType`.

Note that, we treat `PointerType` as a `QualType`, not a `CXType_`.
"""
struct PointerType <: AbstractQualType
    ptr::CXQualType
end

"""
    abstract type AbstractReferenceType <: AbstractQualType
Supertype for `ReferenceType`s.
"""
abstract type AbstractReferenceType <: AbstractQualType end

"""
    struct ReferenceType <: AbstractReferenceType
A `ReferenceType`.

Note that, we treat `ReferenceType` as a `QualType`, not a `CXType_`.
"""
struct ReferenceType <: AbstractReferenceType
    ptr::CXQualType
end

"""
    struct LValueReferenceType <: AbstractReferenceType
A `LValueReferenceType`.

Note that, we treat `LValueReferenceType` as a `QualType`, not a `CXType_`.
"""
struct LValueReferenceType <: AbstractReferenceType
    ptr::CXQualType
end

"""
    struct RValueReferenceType <: AbstractReferenceType
A `RValueReferenceType`.

Note that, we treat `RValueReferenceType` as a `QualType`, not a `CXType_`.
"""
struct RValueReferenceType <: AbstractReferenceType
    ptr::CXQualType
end

"""
    struct MemberPointerType <: AbstractQualType
A `MemberPointerType`.

Note that, we treat `MemberPointerType` as a `QualType`, not a `CXType_`.
"""
struct MemberPointerType <: AbstractQualType
    ptr::CXQualType
end

"""
    abstract type AbstractArrayType <: AbstractQualType
Supertype for `ArrayType`s.
"""
abstract type AbstractArrayType <: AbstractQualType end

"""
    struct ArrayType <: AbstractArrayType
A `ArrayType`.

Note that, we treat `ArrayType` as a `QualType`, not a `CXType_`.
"""
struct ArrayType <: AbstractArrayType
    ptr::CXQualType
end

"""
    struct ConstantArrayType <: AbstractArrayType
A `ConstantArrayType`.

Note that, we treat `ConstantArrayType` as a `QualType`, not a `CXType_`.
"""
struct ConstantArrayType <: AbstractArrayType
    ptr::CXQualType
end

"""
    struct IncompleteArrayType <: AbstractArrayType
A `IncompleteArrayType`.

Note that, we treat `IncompleteArrayType` as a `QualType`, not a `CXType_`.
"""
struct IncompleteArrayType <: AbstractArrayType
    ptr::CXQualType
end

"""
    struct VariableArrayType <: AbstractArrayType
A `VariableArrayType`.

Note that, we treat `VariableArrayType` as a `QualType`, not a `CXType_`.
"""
struct VariableArrayType <: AbstractArrayType
    ptr::CXQualType
end

"""
    struct DependentSizedArrayType <: AbstractArrayType
A `DependentSizedArrayType`.

Note that, we treat `DependentSizedArrayType` as a `QualType`, not a `CXType_`.
"""
struct DependentSizedArrayType <: AbstractArrayType
    ptr::CXQualType
end

"""
    abstract type AbstractFunctionType <: AbstractQualType
Supertype for `FunctionType`s.
"""
abstract type AbstractFunctionType <: AbstractQualType end

"""
    struct FunctionType <: AbstractFunctionType
A `FunctionType`.

Note that, we treat `FunctionType` as a `QualType`, not a `CXType_`.
"""
struct FunctionType <: AbstractFunctionType
    ptr::CXQualType
end

"""
    struct FunctionNoProtoType <: AbstractFunctionType
A `FunctionNoProtoType`.

Note that, we treat `FunctionNoProtoType` as a `QualType`, not a `CXType_`.
"""
struct FunctionNoProtoType <: AbstractFunctionType
    ptr::CXQualType
end

"""
    struct FunctionProtoType <: AbstractFunctionType
A `FunctionProtoType`.

Note that, we treat `FunctionProtoType` as a `QualType`, not a `CXType_`.
"""
struct FunctionProtoType <: AbstractFunctionType
    ptr::CXQualType
end

"""
    struct TypedefType <: AbstractQualType
A `TypedefType`.

Note that, we treat `TypedefType` as a `QualType`, not a `CXType_`.
"""
struct TypedefType <: AbstractQualType
    ptr::CXQualType
end

"""
    abstract type AbstractTagType <: AbstractQualType
Supertype for `TagType`s.
"""
abstract type AbstractTagType <: AbstractQualType end

"""
    struct TagType <: AbstractTagType
A `TagType`.

Note that, we treat `TagType` as a `QualType`, not a `CXType_`.
"""
struct TagType <: AbstractTagType
    ptr::CXQualType
end

"""
    struct RecordType <: AbstractTagType
A `RecordType`.

Note that, we treat `RecordType` as a `QualType`, not a `CXType_`.
"""
struct RecordType <: AbstractTagType
    ptr::CXQualType
end

"""
    struct EnumType <: AbstractTagType
A `EnumType`.

Note that, we treat `EnumType` as a `QualType`, not a `CXType_`.
"""
struct EnumType <: AbstractTagType
    ptr::CXQualType
end

"""
    struct TemplateTypeParmType <: AbstractQualType
A `TemplateTypeParmType`.

Note that, we treat `TemplateTypeParmType` as a `QualType`, not a `CXType_`.
"""
struct TemplateTypeParmType <: AbstractQualType
    ptr::CXQualType
end

"""
    struct SubstTemplateTypeParmType <: AbstractQualType
A `SubstTemplateTypeParmType`.

Note that, we treat `SubstTemplateTypeParmType` as a `QualType`, not a `CXType_`.
"""
struct SubstTemplateTypeParmType <: AbstractQualType
    ptr::CXQualType
end

"""
    struct SubstTemplateTypeParmPackType <: AbstractQualType
A `SubstTemplateTypeParmPackType`.

Note that, we treat `SubstTemplateTypeParmPackType` as a `QualType`, not a `CXType_`.
"""
struct SubstTemplateTypeParmPackType <: AbstractQualType
    ptr::CXQualType
end

"""
    struct TemplateSpecializationType <: AbstractQualType
A `TemplateSpecializationType`.

Note that, we treat `TemplateSpecializationType` as a `QualType`, not a `CXType_`.
"""
struct TemplateSpecializationType <: AbstractQualType
    ptr::CXQualType
end

"""
    abstract type AbstractTypeWithKeyword <: AbstractQualType
Supertype for `TypeWithKeyword`s.
"""
abstract type AbstractTypeWithKeyword <: AbstractQualType end

"""
    struct TypeWithKeyword <: AbstractTypeWithKeyword
A `TypeWithKeyword`.

Note that, we treat `TypeWithKeyword` as a `QualType`, not a `CXType_`.
"""
struct TypeWithKeyword <: AbstractTypeWithKeyword
    ptr::CXQualType
end

"""
    struct ElaboratedType <: AbstractTypeWithKeyword
A `ElaboratedType`.

Note that, we treat `ElaboratedType` as a `QualType`, not a `CXType_`.
"""
struct ElaboratedType <: AbstractTypeWithKeyword
    ptr::CXQualType
end

"""
    struct DependentNameType <: AbstractTypeWithKeyword
A `DependentNameType`.

Note that, we treat `DependentNameType` as a `QualType`, not a `CXType_`.
"""
struct DependentNameType <: AbstractTypeWithKeyword
    ptr::CXQualType
end

"""
    struct DependentTemplateSpecializationType <: AbstractTypeWithKeyword
A `DependentTemplateSpecializationType`.

Note that, we treat `DependentTemplateSpecializationType` as a `QualType`, not a `CXType_`.
"""
struct DependentTemplateSpecializationType <: AbstractTypeWithKeyword
    ptr::CXQualType
end


text:
"""
    Som

Structure to hold all data of a trained SOM.

# Fields:
- `codes::Array{Float64,2}`: 2D-array of codebook vectors. One vector per row
- `xdim::Int`: number of neurons in x-direction
- `ydim::Int`: number of neurons in y-direction
- `numCodes::Int`: total number of neurons
- `grid::Array{Float64,2}`: 2D-array of coordinates of neurons on the map
          (2 columns (x,y)] for rectangular and hexagonal maps
           3 columns (x,y,z) for spherical maps)
"""
mutable struct Som
    codes::Matrix{Float64}
    xdim::Int
    ydim::Int
    numCodes::Int
    grid::Matrix{Float64}

    Som(;
        codes::Matrix{Float64},
        xdim::Int,
        ydim::Int,
        numCodes::Int = xdim * ydim,
        grid::Matrix{Float64},
    ) = new(codes, xdim, ydim, numCodes, grid)
end

Base.copy(som::Som) = Som(
    codes = som.codes,
    xdim = som.xdim,
    ydim = som.ydim,
    numCodes = som.numCodes,
    grid = som.grid,
)


text:
# -----------------------------------------------------------------------
# Radiative Forcing From Methane (based on Myhre et al. equations).
# -----------------------------------------------------------------------

# Create function to account for overlap in CH₄ and N₂O absoprtion bands.
function overlap(M, N)
    d = 1.0 + 2.01e-5*(M*N)^0.75 + 5.31e-15*M*(M*N)^1.52
    return 0.47 * log(d)
end

@defcomp rf_ch4_myhre_fair begin

    N₂O_0           = Parameter()             # Initial (pre-industrial) nitrous oxide concentration (ppb).
    CH₄_0           = Parameter()             # Initial (pre-industrial) methane concentration (ppb).
    scale_CH₄       = Parameter()             # Scaling factor for uncertainty in methane radiative forcing.
    H₂O_share       = Parameter()             # Share of direct methane forcing used to estimate stratospheric water vapor forcing due to methane oxidation.
    CH₄             = Parameter(index=[time]) # Atmospheric methane concentration (ppb).

    forcing_CH₄     = Variable(index=[time])  # Radiative forcing from methane (Wm⁻²).
    forcing_CH₄_H₂O = Variable(index=[time])  # Radiative forcing from stratospheric water vapor due to methane oxidiation (Wm⁻²).

    function run_timestep(p, v, d, t)

        if is_first(t)
            # Set initial radiative forcing to 0.
            v.forcing_CH₄[t] = 0.0
            v.forcing_CH₄_H₂O[t] = 0.0

        else

            # Calculate direct CH₄ radiative forcing.
            direct_rf = 0.036 * (sqrt(p.CH₄[t]) - sqrt(p.CH₄_0)) - (overlap(p.CH₄[t], p.N₂O_0) - overlap(p.CH₄_0, p.N₂O_0))

            # Calculate direct radiative forcing from CH4, accounting for forcing uncertainty.
            v.forcing_CH₄[t] = direct_rf * p.scale_CH₄

            #Calculate indirect CH4 forcing due to production of stratospheric H20 from CH4 oxidation.
            v.forcing_CH₄_H₂O[t] = direct_rf * p.H₂O_share
        end
    end
end


text:
using QuPhys
using Documenter

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

makedocs(;
    modules=[QuPhys],
    authors="Alberto Mercurio",
    repo="https://github.com/albertomercurio/QuPhys.jl/blob/{commit}{path}#{line}",
    sitename="QuPhys.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://albertomercurio.github.io/QuPhys.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/albertomercurio/QuPhys.jl",
    devbranch="main",
)


text:
using JSON

function xyz_to_geometry(xyzfile)
  input_file = []
  #println(joinpath(@__DIR__, xyzfile))
  open(xyzfile) do file
    input_file = readlines(file)
  end

  natoms = parse(Int, input_file[1])
  #println(natoms)
  geometry = input_file[3:end]
  #display(geometry); println()

  coords_vector = Vector{Float64}([])
  symbols = Vector{String}([])
  for geometry_line in geometry
    coords_match = eachmatch(r"([-]?[\d]{1,}\.[\d]{1,})",geometry_line)
    atom_coords = collect(coords_match)
    atom_coords_match = map(x -> x.match, atom_coords)
    append!(coords_vector, parse.(Float64,atom_coords_match[1:3]))

    symbol_match = match(r"[A-Za-z]{1,2}",geometry_line)
    push!(symbols, symbol_match.match)
  end
  #coords = reshape(coords_vector, (3,length(geometry)))
  #coords = transpose(coords)
  coords = coords_vector

  #@assert natoms == size(coords)[1]
  #@assert natoms == length(symbols)

  #display(coords); println()
  #display(symbols); println()

  return coords, symbols
end

function xyz_to_molecule(input, charge = 0)
  geom_json = "{\n"

  #== write in molecule geometry ==#
  coords, symbols = xyz_to_geometry(input)
  geom_json *= "    \"geometry\" : [\n"
  for icoord in 1:length(coords)
    geom_json *= "      "

    value = coords[icoord]
    #println(value)
    if icoord == length(coords)
      geom_json *= "$value\n"
    elseif icoord%3 == 0
      geom_json *= "$value,\n"
    else
      geom_json *= "$value,"
    end
  end
  geom_json *= "    ],\n"

  #== write molecule symbols ==#
  geom_json *= "    \"symbols\" : [\n"
  for iatom in 1:length(symbols)
    geom_json *= "      "
    symbol = symbols[iatom]
    if iatom == length(symbols)
      geom_json *= "\"$symbol\"\n"
    elseif iatom%5 == 0
      geom_json *= "\"$symbol\",\n"
    else
      geom_json *= "\"$symbol\","
    end
  end
  geom_json *= "    ],\n"

  #== write molecule charge ==#
  geom_json *= "    \"molecular_charge\": $charge\n"
  geom_json *= "}\n"

  molecule = JSON.parse(geom_json)
  return molecule
end
export xyz_to_molecule


text:
using Each
using Test

@testset "Each.jl" begin
    @testset "each()" begin
        @test sum.(each(Tuple{Vararg{Number}},((1,1),(1,2),(1,3)))) == (2,3,4)
        @test sum.(each(Tuple{Vararg{Number}},(1,1))) == (2,)
        @test_throws ErrorException each(String,((1,2),(3,4)))
        @test eltype(each(String,"something")) == String
    end
    @testset "@each" begin
        x = [5, 10, 15, 20]
        @test @each x > 0 == true
        @test @each x < 0 == false
        @test @each x % 5 == true
        @test @each x == 5 == false
    end

end


text:
using Test
using NCDatasets

# reference value from Meeus, Jean (1998)
# launch of Sputnik 1

@test CFTime.datetuple_standard(2_436_116 - 2_400_001) == (1957, 10, 4)
@test CFTime.datenum_gregjulian(1957,10,4,true) == 36115

@test CFTime.datenum_gregjulian(333,1,27,false) == -557288

function datenum_datetuple_all_calendars()
    #dayincrement = 11
    dayincrement = 11000

    for (tonum,totuple) in [
        (CFTime.datenum_standard,CFTime.datetuple_standard),
        (CFTime.datenum_julian,CFTime.datetuple_julian),
        (CFTime.datenum_prolepticgregorian,CFTime.datetuple_prolepticgregorian),
        (CFTime.datenum_allleap,CFTime.datetuple_allleap),
        (CFTime.datenum_noleap,CFTime.datetuple_noleap),
        (CFTime.datenum_360day,CFTime.datetuple_360day),
    ]
        for Z = -2_400_000 + CFTime.DATENUM_OFFSET : dayincrement : 600_000 + CFTime.DATENUM_OFFSET
            y,m,d = totuple(Z)
            @test tonum(y,m,d) == Z
        end
    end
end

datenum_datetuple_all_calendars()

# test of DateTime structures

dt = DateTimeNoLeap(1959,12,31,23,39,59,123)
@test dt + Dates.Millisecond(7) == DateTimeNoLeap(1959,12,31,23,39,59,130)
@test dt + Dates.Second(7)      == DateTimeNoLeap(1959,12,31,23,40,6,123)
@test dt + Dates.Minute(7)      == DateTimeNoLeap(1959,12,31,23,46,59,123)
@test dt + Dates.Hour(7)        == DateTimeNoLeap(1960,1,1,6,39,59,123)
@test dt + Dates.Day(7)         == DateTimeNoLeap(1960,1,7,23,39,59,123)
@test dt + Dates.Month(7)       == DateTimeNoLeap(1960,7,31,23,39,59,123)
@test dt + Dates.Year(7)        == DateTimeNoLeap(1966,12,31,23,39,59,123)
@test dt + Dates.Month(24)      == DateTimeNoLeap(1961,12,31,23,39,59,123)

@test dt - Dates.Month(0)       == DateTimeNoLeap(1959,12,31,23,39,59,123)
@test dt - Dates.Month(24)      == DateTimeNoLeap(1957,12,31,23,39,59,123)
@test dt - Dates.Year(7)        == DateTimeNoLeap(1952,12,31,23,39,59,123)

# leap day
@test DateTimeAllLeap(2001,2,28) + Dates.Day(1) == DateTimeAllLeap(2001,2,29)
@test DateTimeNoLeap(2001,2,28)  + Dates.Day(1) == DateTimeNoLeap(2001,3,1)
@test DateTimeJulian(2001,2,28)  + Dates.Day(1) == DateTimeJulian(2001,3,1)
@test DateTimeJulian(1900,2,28)  + Dates.Day(1) == DateTimeJulian(1900,2,29)
@test DateTime360Day(2001,2,28)     + Dates.Day(1) == DateTime360Day(2001,2,29)
@test DateTime360Day(2001,2,29)     + Dates.Day(1) == DateTime360Day(2001,2,30)

@test DateTimeAllLeap(2001,2,29) - DateTimeAllLeap(2001,2,28) == Dates.Day(1)
@test DateTimeNoLeap(2001,3,1)   - DateTimeNoLeap(2001,2,28)  == Dates.Day(1)
@test DateTimeJulian(2001,3,1)   - DateTimeJulian(2001,2,28)  == Dates.Day(1)
@test DateTimeJulian(1900,2,29)  - DateTimeJulian(1900,2,28)  == Dates.Day(1)
@test DateTime360Day(2001,2,29)     - DateTime360Day(2001,2,28)     == Dates.Day(1)
@test DateTime360Day(2001,2,30)     - DateTime360Day(2001,2,29)     == Dates.Day(1)

# reference values from python's cftime
@test DateTimeJulian(2000,1,1) + Dates.Day(1) == DateTimeJulian(2000,01,02)
@test DateTimeJulian(2000,1,1) + Dates.Day(12) == DateTimeJulian(2000,01,13)
@test DateTimeJulian(2000,1,1) + Dates.Day(123) == DateTimeJulian(2000,05,03)
@test DateTimeJulian(2000,1,1) + Dates.Day(1234) == DateTimeJulian(2003,05,19)
@test DateTimeJulian(2000,1,1) + Dates.Day(12345) == DateTimeJulian(2033,10,19)
@test DateTimeJulian(2000,1,1) + Dates.Day(12346) == DateTimeJulian(2033,10,20)
@test DateTimeJulian(1,1,1) + Dates.Day(1234678) == DateTimeJulian(3381,05,14)

# generic tests
function stresstest_DateTime(::Type{DT}) where DT
    t0 = DT(1,1,1)
    for n = -800000:11:800000
        #@show n
        t = t0 + Dates.Day(n)
        y, m, d, h, mi, s, ms = CFTime.datetuple(t)
        @test DT(y, m, d, h, mi, s, ms) == t
    end
end

for DT in [
    DateTimeStandard,
    DateTimeJulian,
    DateTimeProlepticGregorian,
    DateTimeAllLeap,
    DateTimeNoLeap,
    DateTime360Day
]

    dtime = DT(1959,12,30, 23,39,59,123)
    @test Dates.year(dtime) == 1959
    @test Dates.month(dtime) == 12
    @test Dates.day(dtime) == 30
    @test Dates.hour(dtime) == 23
    @test Dates.minute(dtime) == 39
    @test Dates.second(dtime) == 59
    @test Dates.millisecond(dtime) == 123

    @test string(DT(2001,2,20)) == "2001-02-20T00:00:00"
    @test CFTime.datetuple(DT(1959,12,30,23,39,59,123)) == (1959,12,30,23,39,59,123)

    stresstest_DateTime(DT)
end

@test_throws ErrorException DateTime360Day(2010,0,1)
@test_throws ErrorException DateTime360Day(2010,1,0)

# test show
io = IOBuffer()
show(io,DateTimeJulian(-1000,1,1))
isempty(findfirst("Julian",String(take!(io)))) == false

# time

t0,plength = CFTime.timeunits("days since 1950-01-02T03:04:05Z")
@test t0 == DateTimeStandard(1950,1,2, 3,4,5)
@test plength == 86400000

t0,plength = CFTime.timeunits("days since -4713-01-01T00:00:00Z")
@test t0 == DateTimeStandard(-4713,1,1)
@test plength == 86400000

t0,plength = CFTime.timeunits("days since -4713-01-01")
@test t0 == DateTimeStandard(-4713,1,1)
@test plength == 86400000

t0,plength = CFTime.timeunits("days since 2000-01-01 0:0:0")
@test t0 == DateTimeStandard(2000,1,1)
@test plength == 86400000

# issue 24
t0,plength = CFTime.timeunits("hours since 1900-01-01 00:00:00.0")
@test t0 == DateTimeStandard(1900,1,1)
@test plength == 86400000 ÷ 24

t0,plength = CFTime.timeunits("seconds since 1992-10-8 15:15:42.5")
@test t0 == DateTimeStandard(1992,10,8,15,15,42,500)
@test plength == 1000

for (calendar,DT) in [
    ("standard",DateTimeStandard),
    ("gregorian",DateTimeStandard),
    ("proleptic_gregorian",DateTimeProlepticGregorian),
    ("julian",DateTimeJulian),
    ("noleap",DateTimeNoLeap),
    ("365_day",DateTimeNoLeap),
    ("all_leap",DateTimeAllLeap),
    ("366_day",DateTimeAllLeap),
    ("360_day",DateTime360Day)]

    calendart0,calendarplength = CFTime.timeunits("days since 2000-1-1 0:0:0",calendar)
    @test calendart0 == DT(2000,1,1)
    @test calendarplength == 86400000
end

@test_throws ErrorException CFTime.timeunits("fortnights since 2000-01-01")
@test_throws ErrorException CFTime.timeunits("days since 2000-1-1 0:0:0","foo")

# value from python's cftime
# print(cftime.DatetimeJulian(-4713,1,1) + datetime.timedelta(2455512,.375 * 24*60*60))
# 2010-10-29 09:00:00

#@test timedecode([2455512.375],"days since -4713-01-01T00:00:00","julian")
#   == DateTimeJulian(2010,10,29,09)

# values from
# https://web.archive.org/web/20180212214229/https://en.wikipedia.org/wiki/Julian_day

# Modified JD
@test CFTime.timedecode([58160.6875],"days since 1858-11-17","standard") ==
    [DateTime(2018,2,11,16,30,0)]

# CNES JD
@test CFTime.timedecode([24878.6875],"days since 1950-01-01","standard") ==
    [DateTime(2018,2,11,16,30,0)]

# Unix time
# wikipedia pages reports 1518366603 but it should be 1518366600
@test CFTime.timedecode([1518366600],"seconds since 1970-01-01","standard") ==
    [DateTime(2018,2,11,16,30,0)]

# The Julian Day Number (JDN) is the integer assigned to a whole solar day in
# the Julian day count starting from noon Universal time, with Julian day
# number 0 assigned to the day starting at noon on Monday, January 1, 4713 BC,
# proleptic Julian calendar (November 24, 4714 BC, in the proleptic Gregorian
# calendar),

# Julian Day Number of 12:00 UT on January 1, 2000, is 2 451 545
# https://web.archive.org/web/20180613200023/https://en.wikipedia.org/wiki/Julian_day

@test CFTime.timedecode(DateTimeStandard,2_451_545,"days since -4713-01-01T12:00:00") ==
    DateTimeStandard(2000,01,01,12,00,00)

# Note for DateTime, 1 BC is the year 0!
# DateTime(1,1,1)-Dates.Day(1)
# 0000-12-31T00:00:00

@test CFTime.timedecode(DateTime,2_451_545,"days since -4713-11-24T12:00:00") ==
    DateTime(2000,01,01,12,00,00)

@test CFTime.timedecode(DateTimeProlepticGregorian,2_451_545,"days since -4714-11-24T12:00:00") ==
    DateTimeProlepticGregorian(2000,01,01,12,00,00)

@test CFTime.timedecode([2455512.375],"days since -4713-01-01T00:00:00","julian", prefer_datetime = false) ==
    [DateTimeJulian(2010,10,29,9,0,0)]

@test CFTime.timeencode([DateTimeJulian(2010,10,29,9,0,0)],"days since -4713-01-01T00:00:00","julian") ==
    [2455512.375]

@test CFTime.timedecode(DateTime,[22280.0f0],"days since 1950-01-01 00:00:00") == [DateTime(2011,1,1)]

@test_throws ErrorException CFTime.timeencode(
    [DateTimeJulian(2010,10,29,9,0,0)],
    "days since -4713-01-01T00:00:00","360_day")

# Transition between Julian and Gregorian Calendar

#=
In [11]: cftime.DatetimeGregorian(1582,10,4) + datetime.timedelta(1)
Out[11]: cftime.DatetimeGregorian(1582, 10, 15, 0, 0, 0, 0, -1, 1)

In [12]: cftime.DatetimeProlepticGregorian(1582,10,4) + datetime.timedelta(1)
Out[12]: cftime.DatetimeProlepticGregorian(1582, 10, 5, 0, 0, 0, 0, -1, 1)

In [13]: cftime.DatetimeJulian(1582,10,4) + datetime.timedelta(1)
Out[13]: cftime.DatetimeJulian(1582, 10, 5, 0, 0, 0, 0, -1, 1)
=#

@test DateTimeStandard(1582,10,4) + Dates.Day(1) == DateTimeStandard(1582,10,15)
@test DateTimeProlepticGregorian(1582,10,4) + Dates.Day(1) == DateTimeProlepticGregorian(1582,10,5)
@test DateTimeJulian(1582,10,4) + Dates.Day(1) == DateTimeJulian(1582,10,5)

@test CFTime.datetuple(CFTime.timedecode(0,"days since -4713-01-01T12:00:00","julian", prefer_datetime = false)) ==
    (-4713, 1, 1, 12, 0, 0, 0)

dt = NCDatasets.reinterpret(DateTimeStandard, DateTimeJulian(1900,2,28))
@test typeof(dt) == DateTimeStandard
@test CFTime.datetuple(dt) == (1900,2,28,0, 0, 0, 0)

dt = NCDatasets.reinterpret(DateTime, DateTimeNoLeap(1900,2,28))
@test typeof(dt) == DateTime
@test Dates.year(dt) == 1900
@test Dates.month(dt) == 2
@test Dates.day(dt) == 28

# check ordering

@test DateTimeStandard(2000,01,01) < DateTimeStandard(2000,01,02)
@test DateTimeStandard(2000,01,01) ≤ DateTimeStandard(2000,01,01)

@test DateTimeStandard(2000,01,03) > DateTimeStandard(2000,01,02)
@test DateTimeStandard(2000,01,03) ≥ DateTimeStandard(2000,01,01)

import NCDatasets.CFTime: datetuple
datetuple(dt::DateTime) = (Dates.year(dt),Dates.month(dt),Dates.day(dt),
                           Dates.hour(dt),Dates.minute(dt),Dates.second(dt),
                           Dates.millisecond(dt))

# check conversion

for T1 in [DateTimeProlepticGregorian,DateTimeStandard,DateTime]
    for T2 in [DateTimeProlepticGregorian,DateTimeStandard,DateTime]
        # datetuple should not change after 1582-10-15
        # for Gregorian Calendars
        dt1 = T1(2000,01,03)
        dt2 = convert(T2,dt1)

        @test CFTime.datetuple(dt1) == CFTime.datetuple(dt2)
    end
end

for T1 in [DateTimeStandard,DateTimeJulian]
    for T2 in [DateTimeStandard,DateTimeJulian]
        # datetuple should not change before 1582-10-15
        # for Julian Calendars
        dt1 = T1(200,01,03)
        dt2 = convert(T2,dt1)

        @test CFTime.datetuple(dt1) == CFTime.datetuple(dt2)
    end
end

for T1 in [DateTimeProlepticGregorian,DateTimeJulian,DateTimeStandard,DateTime]
    for T2 in [DateTimeProlepticGregorian,DateTimeJulian,DateTimeStandard,DateTime]
        # verify that durations (even accross 1582-10-15) are maintained
        # after convert
        dt1 = [T1(2000,01,03), T1(-100,2,20)]
        dt2 = convert.(T2,dt1)
        @test dt1[2]-dt1[1] == dt2[2]-dt2[1]
    end
end

# issue #12

units = "days since 1850-01-01 00:00:00"
calendar = "noleap"
data_orig = [54750.5, 54751.5, 54752.5]

# Decoding
datacal = CFTime.timedecode(data_orig, units, calendar)
# Reencoding
data_orig_back = CFTime.timeencode(datacal, units, calendar)
@test data_orig ≈ data_orig_back

# issue #17

# reference values from cftime
# for T in [cftime.DatetimeGregorian,cftime.DatetimeJulian,cftime.DatetimeProlepticGregorian,cftime.DatetimeAllLeap,cftime.DatetimeNoLeap, cftime.Datetime360Day]:
#     print(T,T(1582,11,1) - T(1582,10,1))

@test daysinmonth(DateTimeStandard(1582,10,1)) == 21
@test daysinmonth(DateTimeJulian(1582,10,1)) == 31
@test daysinmonth(DateTimeProlepticGregorian(1582,10,1)) == 31
@test daysinmonth(DateTimeAllLeap(1582,10,1)) == 31
@test daysinmonth(DateTimeNoLeap(1582,10,1)) == 31
@test daysinmonth(DateTime360Day(1582,10,1)) == 30

# import cftime
# for T in [cftime.DatetimeGregorian,cftime.DatetimeJulian,cftime.DatetimeProlepticGregorian,cftime.DatetimeAllLeap,cftime.DatetimeNoLeap, cftime.Datetime360Day]:
#    print(T,T(1583,1,1) - T(1582,1,1))

@test daysinyear(DateTimeStandard(1582,10,1)) == 355
@test daysinyear(DateTimeJulian(1582,10,1)) == 365
@test daysinyear(DateTimeProlepticGregorian(1582,10,1)) == 365
@test daysinyear(DateTimeAllLeap(1582,10,1)) == 366
@test daysinyear(DateTimeNoLeap(1582,10,1)) == 365
@test daysinyear(DateTime360Day(1582,10,1)) == 360

for T in [DateTimeStandard, DateTimeJulian, DateTimeProlepticGregorian,
          DateTimeAllLeap, DateTimeNoLeap, DateTime360Day]
    @test Dates.yearmonthday(T(2004,1,2)) == (2004, 1, 2)
    @test Dates.yearmonth(T(2004,1,2)) == (2004, 1)
    @test Dates.monthday(T(2004,1,2)) == (1, 2)

    # test constructor with argument which are not Int64
    @test T(Int16(2000),Int32(1),UInt8(1)) == T(2000,1,1)
end

# time ranges

@test length(DateTimeNoLeap(2000, 01, 01):Dates.Day(1):DateTimeNoLeap(2000, 12, 31)) == 365
@test length(DateTimeNoLeap(2000, 01, 01):Dates.Month(1):DateTimeNoLeap(2000, 12, 31)) == 12

for T in [DateTimeStandard, DateTimeJulian, DateTimeProlepticGregorian,
          DateTimeAllLeap, DateTimeNoLeap, DateTime360Day]

    # end date is inclusive
    @test length(T(2000, 01, 01):Dates.Month(1):T(2001, 1, 1)) == 13
    @test length(T(2000, 01, 01):Dates.Year(1):T(2001, 1, 1)) == 2
end

# if @isdefined DataArrays

# # DataArray
# data_orig = DataArrays.DataArray([54750.5, 54751.5, 54752.5])
# # Decoding
# datacal = CFTime.timedecode(data_orig, units, calendar)
# # Reencoding
# data_orig_back = CFTime.timeencode(datacal, units, calendar)
# @test data_orig ≈ data_orig_back
# end

# issue 21
@test parse(DateTimeNoLeap,"1999-12-05", dateformat"yyyy-mm-dd") == DateTimeNoLeap(1999,12,05)
@test DateTimeNoLeap("1999-12-05", "yyyy-mm-dd") == DateTimeNoLeap(1999,12,05)
@test DateTimeNoLeap("1999-12-05", dateformat"yyyy-mm-dd") == DateTimeNoLeap(1999,12,05)

# issue 29
@test Dates.firstdayofyear(DateTimeNoLeap(2008, 12, 31)) == DateTimeNoLeap(2008, 1, 1)
@test Dates.dayofyear(DateTimeNoLeap(2008, 12, 31)) == 365
@test Dates.dayofmonth(DateTimeAllLeap(2008, 2, 29)) == 29


text:
@testset "normal sampler" begin

	cc = dummyChoiceContext()

	@testset "getparams, numparams, paramranges, sample" begin

		@testset "using default construction" begin

		    s = DataGenerators.NormalSampler()

		    @testset "numparams and paramranges" begin
		        @test DataGenerators.numparams(s) == 2
		        prs = DataGenerators.paramranges(s)
		        @test typeof(prs) <: Vector{Tuple{Float64,Float64}}
		        @test prs == [(-realmax(Float64), realmax(Float64)), (0.0, realmax(Float64))]
			end

		    @testset "default params" begin
		        @test DataGenerators.getparams(s) == [0.0, 1.0]
		    end

		    @mtestset "default sampling" reps=Main.REPS alpha=Main.ALPHA begin
		        x, trace = DataGenerators.sample(s, (0,1), cc)
		        @test typeof(x) <: Float64
				@mtest_that_sometimes x < -1.0
				@mtest_that_sometimes -1.0 <= x < 0.0
				@mtest_that_sometimes 0.0 < x <= 1.0
				@mtest_that_sometimes 1.0 < x
		        @mtest_distributed_as Normal(0.0,1.0) x
			end

		end

		@testset "using non-default construction" begin

		    s = DataGenerators.NormalSampler([-949.88, 123.4])
		    @test DataGenerators.getparams(s) == [-949.88, 123.4]

		    @mtestset "consistent with normal" reps=Main.REPS alpha=Main.ALPHA begin
		    	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Normal(-949.88,123.4) x
		    end

		end

	end

	@testset "setparams" begin

	    s = DataGenerators.NormalSampler()
	    prs = DataGenerators.paramranges(s)
	    midparams = map(pr->robustmidpoint(pr[1],pr[2]), prs)

	    @testset "valid parameters $params" for params in [[0.8, 19.3], [87.4, 0.0],]
		# note special case of sigma = 0

	        DataGenerators.setparams!(s, params)

			@mtestset "consistent with normal" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
				mu = params[1]
				sigma = params[2] == 0.0 ? eps(params[2]) : params[2]
		        @mtest_distributed_as Normal(mu,sigma) x
			end

	    end

	    @testset "boundary values for parameter index $pidx" for pidx in 1:length(prs)

	        pr = prs[pidx]
	        params = copy(midparams)

			@testset "bound index $bidx" for bidx in 1:2

	            params[pidx] = pr[bidx]
	            DataGenerators.setparams!(s, params)

				@mtestset "consistent with normal" reps=Main.REPS alpha=Main.ALPHA begin
	            	x, trace = DataGenerators.sample(s, (0,1), cc)
					mu = params[1]
					sigma = params[2] == 0.0 ? eps(params[2]) : params[2]
			        @mtest_distributed_as Normal(mu,sigma) x
				end

				@testset "range check exception" begin
		            params[pidx] = bidx == 1 ? prevfloat(pr[bidx]) : nextfloat(pr[bidx])
		            @test_throws ErrorException DataGenerators.setparams!(s, params)
				end

			end

	    end

	    @testset "wrong number of parameters" begin

	        @test_throws ErrorException DataGenerators.setparams!(s, midparams[1:end-1])
	        @test_throws ErrorException DataGenerators.setparams!(s, [midparams; 0.5])

	    end

	end

	@testset "estimateparams" begin

	    otherparams = [155.6, 13.7]

	    @testset "from parameters $params" for params in [[-56.0, 18.9], [3.44, 0.0],]
		# note special case of sigma = 0

	        s1 = DataGenerators.NormalSampler(params)
	        s2 = DataGenerators.NormalSampler(otherparams)
	        traces = map(1:100) do i
	            x, trace = DataGenerators.sample(s1, (0,1), cc)
	            trace
	        end
	        estimateparams!(s2, traces)

			@mtestset "consistent with normal" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s2, (0,1), cc)
				if params[2] > 0.0
					@mtest_distributed_as Normal(params[1],params[2]) x
				else
					@test isapprox(x, params[1])
					# necessary since above distributed_as fails when sigma=0.0 owing to rounding errors
				end
			end

	    end

	    @testset "too few traces" begin

	        params = [21.0, 0.7]
	        s1 = DataGenerators.NormalSampler(params)
	        s2 = DataGenerators.NormalSampler(otherparams)
	        traces = map(1:1) do i
	            x, trace = DataGenerators.sample(s1, (0,1), cc)
	            trace
	        end

			@mtestset "consistent with normal" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s2, (0,1), cc)
		        @mtest_distributed_as Normal(otherparams[1],otherparams[2]) x
			end

	    end

	end

end


text:
function true_test(model::MOI.ModelLike, config::MOIT.Config)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(model, MOI.VectorOfVariables, CP.True)

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x1, x1]),
        CP.IfThenElse(CP.True(), MOI.EqualTo(1), MOI.EqualTo(0)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, c1)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
    end
end

function false_test(model::MOI.ModelLike, config::MOIT.Config)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(model, MOI.VectorOfVariables, CP.False)

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x1, x1]),
        CP.IfThenElse(CP.False(), MOI.EqualTo(1), MOI.EqualTo(0)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, c1)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 0
    end
end

const truefalsetests = Dict("true" => true_test, "false" => false_test)

MOIT.@moitestset truefalse


text:
using PipeLayout
import PipeLayout: has_geosteiner, euclidean_steiner_tree

if has_geosteiner()
    @testset "compute ESMT using GeoSteiner" begin
        nodes = [Node(0,0), Node(2,0), Node(1,1)]
        topo = euclidean_steiner_tree(nodes)
        @test length(topo.nodes) == 4
        @test topo.nodes[1] == nodes[1]
        @test topo.nodes[2] == nodes[2]
        @test topo.nodes[3] == nodes[3]
        @test topo.nodes[4].x ≈ 1.0
        @test topo.nodes[4].y ≈ 0.57735 atol=1e-5

        @test length(topo.arcs) == 3
        for arc in topo.arcs
            # head or tail is steiner node
            @test arc.tail == 4 || arc.head == 4
            # other node is terminal
            @test arc.tail in [1,2,3] || arc.head in [1,2,3]
        end
    end
end


text:
export NLS_ForwardDiff, create_NLS_problem_using_ForwardDiff

using ForwardDiff: jacobian

@doc raw"""

```julia
struct NLS_ForwardDiff <: AbstractNLS
    ...
end

```

A specialization that uses the `ForwardDiff` package to compute the Jacobian.

By comparison with [`AbstractNLS`](@ref) you only have to define these
functions:
- [`parameter_size`](@ref) : returns ``n_θ``
- [`residue_size`](@ref) : returns ``n_S``
- [`eval_r`](@ref) : computation of ``\mathbf{r}``

See: [`create_NLS_problem_using_ForwardDiff`](@ref)
"""
struct NLS_ForwardDiff <: AbstractNLS
    _eval_r_function::Function
    _residue_size::Int
    _parameter_size::Int
end

@doc raw"""
```julia
create_NLS_problem_using_ForwardDiff(r::Function;domain_image_dim::Pair{Int,Int})
```

Creates an [`AbstractNLS`](@ref) specialized instance where the
[`eval_r_J`](@ref) function is automatically defined using automatic
differentiation.

- `r` is a function that maps a parameter vector θ to its residue. The
  Jacobian matrix is computed using the `ForwardDiff` package.
- `domain_image_dim` is a pair of the form `θ length => r length` that
  defines domain and codomain dimensions.

# Usage example

An example defining the Rosenbrock function

```math
\frac{1}{2}\|r(\theta)\|^2\text{ where }r = \sqrt{2} \left( \begin{array}{c}  1-\theta_1 \\ 10(\theta_2-\theta_1^2) \end{array} \right)
```

```julia
nls = create_NLS_problem_using_ForwardDiff(2 => 2) do θ
     sqrt(2) sqrt(2)* [ 1-θ[1], 10*(θ[2]-θ[1]^2) ]* [ 1-θ[1], 10*(θ[2]-θ[1]^2) ]
end
```
"""
function create_NLS_problem_using_ForwardDiff(r::Function,domain_image_dim::Pair{Int,Int})
    NLS_ForwardDiff(r,last(domain_image_dim),first(domain_image_dim))
end

parameter_size(nls::NLS_ForwardDiff) = nls._parameter_size
residue_size(nls::NLS_ForwardDiff) = nls._residue_size

function eval_r(nls::NLS_ForwardDiff,θ::AbstractVector)
    nls._eval_r_function(θ)
end

function eval_r_J(nls::NLS_ForwardDiff, θ::AbstractVector{T}) where {T}

    r_evaluation = (r,θ)->(r.=eval_r(nls,θ))

    r=Vector{T}(undef,residue_size(nls))

    J = jacobian(r_evaluation, r, θ)

    r,J
end


text:
using PowerSystems
using JuMP
using Ipopt
using PowerSimulations
const PS = PowerSimulations

ipopt_optimizer = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, print_level=0)

base_dir = dirname(dirname(pathof(PowerSystems)))
include(joinpath(base_dir,"data/data_5bus_pu.jl"))
sys5 = PowerSystem(nodes5, generators5, loads5_DA, branches5, nothing,  100.0);

@test try
    @info "testing copper plate"
    Net = PS.CopperPlatePowerModel
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys5);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys5)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 3400, atol = 1000)
finally end

# Flow Models
@test try
    @info "testing net flow"
    Net = PS.StandardNetFlow
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys5);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys5)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 3400, atol = 1000)
finally end

# Flow Models
@test try
    @info "testing PTDF 5-bus"
    Net = PS.StandardPTDF
    m = Model(ipopt_optimizer);
    ptdf,  A = PowerSystems.buildptdf(sys5.branches, sys5.buses)
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys5);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys5, PTDF = ptdf)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 3400, atol = 1000)
finally end

@test try
    @info "testing AngleDC-OPF 5-bus"
    Net = PS.DCAngleForm
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys5);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys5)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 3400, atol = 1000)
finally end

@test try
    @info "testing ACP-OPF 5-bus"
    Net = PS.StandardAC
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys5);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys5)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 3400, atol = 1000)
true finally end

@test try
    @info "testing ACP- QCWForm 5-bus"
    Net = PM.QCWRForm
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys5);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys5)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 3400, atol = 1000)
true finally end

include(joinpath(base_dir,"data/data_14bus_pu.jl"))
sys14 = PowerSystem(nodes14, generators14, loads14, branches14, nothing,  100.0);

@test try
    @info "testing copper plate 14-bus"
    Net = PS.CopperPlatePowerModel
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys14);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys14);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys14)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 1200, atol = 1000)
finally end

# Flow Models
@test try
    @info "testing net 14-bus"
    Net = PS.StandardNetFlow
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys14);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys14);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys14)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 1200, atol = 1000)
finally end

# Flow Models
@test_skip try
    @info "testing PTDF 14-bus"
    Net = PS.StandardPTDF
    m = Model(ipopt_optimizer);
    ptdf,  A = PowerSystems.buildptdf(sys14.branches, sys14.buses)
    netinjection = PS.instantiate_network(Net, sys14);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys14);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys14, PTDF = ptdf)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 1200, atol = 1000)
finally end

@test try
    @info "testing AngleDC-OPF 14-bus"
    Net = PS.DCAngleForm
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys14);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys14);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys14)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 1200, atol = 1000)
finally end

@test try
    @info "testing ACP-OPF 14-bus"
    Net = PS.StandardAC
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys14);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys14);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys14)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 1200, atol = 1000)
true finally end

@test try
    @info "testing ACP-QCWForm 14-bus"
    Net = PM.QCWRForm
    m = Model(ipopt_optimizer);
    netinjection = PS.instantiate_network(Net, sys14);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys14);
    PS.constructnetwork!(m, [(device=Line, formulation=PS.PiLine)], netinjection, Net, sys14)
    @objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 1200, atol = 1000)
true finally end


text:
function printPassed()
  print_with_color(:green, "PASSED\n")
end


text:
#=
This module contains functions that generate NormalFormGame instances of the
2-player games studied by Fearnley, Igwe, and Savani (2015):

* Colonel Blotto Games (`blotto_game`): A non-zero sum extension of the Blotto
  game as studied by Hortala-Vallve and Llorente-Saguer (2012), where opposing
  parties have asymmetric and heterogeneous battlefield valuations.

* Ranking Games (`ranking_game`): In these games, as studied by Goldberg et al.
  (2013), each player chooses an effort level associated with a cost and a
  score. The players are ranked according to their scores, and the player with
  the higher score wins the prize. Each player's payoff is given by the value of
  the prize minus the cost of the effort.

* SGC Games (`sgc_game`): These games were introduced by Sandholm, Gilpin, and
  Conitzer (2005) as a worst case scenario for support enumeration as it has a
  unique equilibrium where each player uses half of his actions in his support.

* Tournament Games (`tournament_game`): These games are constructed by
  Anbalagan et al. (2013) as games that do not have interim epsilon-Nash
  equilibria with constant cardinality supports for epsilon smaller than a
  certain threshold.

* Unit Vector Games (`unit_vector_game`): These games are games where the payoff
  matrix of one player consists of unit (column) vectors, used by Savani and von
  Stengel (2016) to construct instances that are hard, in terms of computational
  complexity, both for the Lemke-Howson and support enumeration algorithms.

Large part of the code here is based on the C code available at
https://github.com/bimatrix-games/bimatrix-generators distributed under BSD
3-Clause License.

References
----------
* Y. Anbalagan, S. Norin, R. Savani, and A. Vetta, "Polylogarithmic Supports
  Are Required for Approximate Well-Supported Nash Equilibria below 2/3," WINE,
  2013.

* J. Fearnley, T. P. Igwe, and R. Savani, "An Empirical Study of Finding
  Approximate Equilibria in Bimatrix Games," International Symposium on
  Experimental Algorithms (SEA), 2015.

* L.A. Goldberg, P.W. Goldberg, P. Krysta, and C. Ventre, "Ranking Games that
  have Competitiveness-based Strategies", Theoretical Computer Science, 2013

* R. Hortala-Vallve and A. Llorente-Saguer, "Pure Strategy Nash Equilibria in
  Non-Zero Sum Colonel Blotto Games", International Journal of Game Theory,
  2012.

* T. Sandholm, A. Gilpin, and V. Conitzer, "Mixed-Integer Programming Methods
  for Finding Nash Equilibria," AAAI, 2005.

* R. Savani and B. von Stengel, "Unit Vector Games," International Journal of
  Economic Theory, 2016.

=#

# BSD 3-Clause License
#
# Copyright (c) 2015, John Fearnley, Tobenna P. Igwe, Rahul Savani
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# * Neither the name of the copyright holder nor the names of its
#   contributors may be used to endorse or promote products derived from
#   this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

using QuantEcon: MVNSampler, simplex_grid, next_k_array!, k_array_rank
using LightGraphs: random_tournament_digraph

# blotto_game
"""
    blotto_game([rng=GLOBAL_RNG], h, t, rho[, mu=0])

Return a NormalFormGame instance of a 2-player non-zero sum Colonel Blotto game
(Hortala-Vallve and Llorente-Saguer, 2012), where the players have an equal
number `t` of troops to assign to `h` hills (so that the number of actions for
each player is equal to (t+h-1) choose (h-1) = (T+h-1)!/(T!*(h-1)!)). Each
player has a value for each hill that he receives if he assigns strictly more
troops to the hill than his opponent (ties are broken uniformly at random),
where the values are drawn from a multivariate normal distribution with
covariance `rho`. Each player’s payoff is the sum of the values of the hills won
by that player.

# Arguments

- `rng::AbstractRNG=GLOBAL_RNG`: Random number generator used.
- `h::Integer` : Number of hills.
- `t::Integer` : Number of troops.
- `rho::Real` : Covariance of the players' values of each hill. Must be in
  [-1, 1].
- `mu::Real=0` : Mean of the players' values of each hill.

# Returns

- `g::NormalFormGame`

# Examples

```julia
julia> rng = MersenneTwister(1234);

julia> g = blotto_game(rng, 2, 3, 0.5)
4×4 NormalFormGame{2,Float64}

julia> g.players[1]
4×4 Player{2,Float64}:
 0.186434  -0.494479  -0.494479  -0.494479
 0.867347   0.186434  -0.494479  -0.494479
 0.867347   0.867347   0.186434  -0.494479
 0.867347   0.867347   0.867347   0.186434

julia> g.players[2]
4×4 Player{2,Float64}:
 -0.688223  -1.02919   -1.02919   -1.02919
 -0.347259  -0.688223  -1.02919   -1.02919
 -0.347259  -0.347259  -0.688223  -1.02919
 -0.347259  -0.347259  -0.347259  -0.688223
```
"""
function blotto_game(rng::AbstractRNG, h::Integer, t::Integer, rho::Real,
                     mu::Real=0)
    d_mean = [mu; mu]
    d_cov = [1 rho; rho 1]
    dist = MVNSampler(d_mean, d_cov)
    values = rand(rng, dist, h)

    actions = simplex_grid(h, t)
    n = size(actions)[2]

    payoff_arrays = [Array{Float64}(undef, n, n) for i in 1:2]

    payoffs = Array{Float64}(undef, 2)
    @inbounds for i = 1:n, j = 1:n
        fill!(payoffs, 0)
        for k = 1:h
            if actions[k, i] == actions[k, j]
                for p = 1:2
                    payoffs[p] += values[p, k] / 2
                end
            else
                winner = 1 + Int(actions[k, i] < actions[k, j])
                payoffs[winner] += values[winner, k]
            end
        end
        payoff_arrays[1][i, j] = payoffs[1]
        payoff_arrays[2][j, i] = payoffs[2]
    end

    g = NormalFormGame(
        [Player(payoff_array) for payoff_array::Array in payoff_arrays]
    )
    return g
end

blotto_game(h::Integer, t::Integer, rho::Real, mu::Real=0) =
    blotto_game(Random.GLOBAL_RNG, h, t, rho, mu)

# ranking_game
"""
    ranking_game([rng=GLOBAL_RNG], n[, steps=10])

Return a NormalFormGame instance of (the 2-player version of) the "ranking game"
studied by Goldberg et al. (2013), where each player chooses an effort level
associated with a score and a cost which are both increasing functions with
randomly generated step sizes. The player with the higher score wins the first
prize, whose value is 1, and the other player obtains the "second prize" of
value 0; in the case of a tie, the first prize is split and each player receives
a value of 0.5. The payoff of a player is given by the value of the prize minus
the cost of the effort.

# Arguments

- `rng::AbstractRNG=GLOBAL_RNG`: Random number generator used.
- `n::Integer` : Number of actions, i.e, number of possible effort levels.
- `steps::Integer=10` : Parameter determining the upper bound for the size of
  the random steps for the scores and costs for each player: The step sizes for
  the scores are drawn from `1`, ..., `steps`, while those for the costs are
  multiples of `1/(n*steps)`, where the cost of effort level `1` is 0, and the
  maximum possible cost of effort level `n` is less than or equal to 1.

# Returns

- `g::NormalFormGame`

# Examples

```julia
julia> rng = MersenneTwister(1234);

julia> g = ranking_game(rng, 5)
5×5 NormalFormGame{2,Float64}

julia> g.players[1]
5×5 Player{2,Float64}:
 0.5    0.0    0.0    0.0    0.0
 0.84  -0.16  -0.16  -0.16  -0.16
 0.8    0.8   -0.2   -0.2   -0.2
 0.78   0.78   0.78  -0.22  -0.22
 0.6    0.6    0.6    0.6   -0.4

julia> g.players[2]
5×5 Player{2,Float64}:
 0.5   0.0    0.0    0.0    0.0
 0.84  0.84  -0.16  -0.16  -0.16
 0.7   0.7    0.7   -0.3   -0.3
 0.5   0.5    0.5    0.5   -0.5
 0.36  0.36   0.36   0.36   0.36
```
"""
function ranking_game(rng::AbstractRNG, n::Integer, steps::Integer=10)
    payoff_arrays = [Array{Float64}(undef, n, n) for i in 1:2]

    scores = rand(rng, 1:steps, (n, 2))
    cumsum!(scores, scores, dims=1)

    costs = Array{Float64}(undef, n-1, 2)
    rand!(rng, costs, 1:steps)
    cumsum!(costs, costs, dims=1)
    costs ./= n * steps

    for (p, payoff_array) in enumerate(payoff_arrays)
        fill!(view(payoff_array, 1, :), 0)
        for j in 1:n, i in 2:n
            @inbounds payoff_array[i, j] = -costs[i-1, p]
        end
    end

    prize = 1.
    for j in 1:n, i in 1:n
        if scores[i, 1] > scores[j, 2]
            payoff_arrays[1][i, j] += prize
        elseif scores[i, 1] < scores[j, 2]
            payoff_arrays[2][j, i] += prize
        else
            payoff_arrays[1][i, j] += prize / 2
            payoff_arrays[2][j, i] += prize / 2
        end
    end

    g = NormalFormGame(
        [Player(payoff_array) for payoff_array::Array in payoff_arrays]
    )
    return g
end

ranking_game(n::Integer, steps::Integer=10) =
    ranking_game(Random.GLOBAL_RNG, n, steps)

# sgc_game
"""
    sgc_game(k)

Return a NormalFormGame instance of the 2-player game introduced by Sandholm,
Gilpin, and Conitzer (2005), which has a unique Nash equilibrium, where each
player plays half of the actions with positive probabilities. Payoffs are
normalized so that the minimum and the maximum payoffs are 0 and 1,
respectively.

# Arguments

- `k::Integer` : Positive integer determining the number of actions. The
  returned game will have `4*k-1` actions for each player.

# Returns

- `g::NormalFormGame`

# Examples

```julia
julia> g = sgc_game(2)
7×7 NormalFormGame{2,Float64}

julia> g.players[1]
7×7 Player{2,Float64}:
 0.75  0.5   1.0   0.5   0.5   0.5   0.5
 1.0   0.75  0.5   0.5   0.5   0.5   0.5
 0.5   1.0   0.75  0.5   0.5   0.5   0.5
 0.0   0.0   0.0   0.75  0.0   0.0   0.0
 0.0   0.0   0.0   0.0   0.75  0.0   0.0
 0.0   0.0   0.0   0.0   0.0   0.75  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.75

julia> g.players[2]
7×7 Player{2,Float64}:
 0.75  0.5   1.0   0.5   0.5   0.5   0.5
 1.0   0.75  0.5   0.5   0.5   0.5   0.5
 0.5   1.0   0.75  0.5   0.5   0.5   0.5
 0.0   0.0   0.0   0.0   0.75  0.0   0.0
 0.0   0.0   0.0   0.75  0.0   0.0   0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.75
 0.0   0.0   0.0   0.0   0.0   0.75  0.0
```
"""
function sgc_game(k::Integer)
    n, m = 4*k-1, 2*k-1
    payoff_arrays = [Array{Float64}(undef, n, n) for i in 1:2]

    for payoff_array in payoff_arrays
        for j in 1:m
            for i in 1:m
                payoff_array[i, j] = 0.75
            end
            for i in (m+1):n
                payoff_array[i, j] = 0
            end
        end
        for j in (m+1):n
            for i in 1:m
                payoff_array[i, j] = 0.5
            end
            for i in (m+1):n
                payoff_array[i, j] = 0
            end
        end

        payoff_array[1, m] = 1
        payoff_array[1, 2] = 0.5
        for i in 2:(m-1)
            payoff_array[i, i-1] = 1
            payoff_array[i, i+1] = 0.5
        end
        payoff_array[m, m-1] = 1
        payoff_array[m, 1] = 0.5
    end

    for h in 0:(k-1)
        i, j = m + 1 + 2*h, m + 1 + 2*h
        payoff_arrays[1][i, j] = payoff_arrays[1][i+1, j+1] = 0.75
        payoff_arrays[2][j, i+1] = payoff_arrays[2][j+1, j] = 0.75
    end

    g = NormalFormGame(
        [Player(payoff_array) for payoff_array::Array in payoff_arrays]
    )
    return g
end

# tournament_game
"""
    tournament_game(n, k; seed=-1)

Return a NormalFormGame instance of the 2-player win-lose game, whose payoffs
are either 0 or 1, introduced by Anbalagan et al. (2013). Player 1 has n
actions, which constitute the set of nodes {1, ..., n}, while player 2 has n
choose k actions, each corresponding to a subset of k elements of the set of n
nodes. Given a randomly generated tournament graph on the n nodes, the payoff
for player 1 is 1 if, in the tournament, the node chosen by player 1 dominates
all the nodes in the k-subset chosen by player 2. The payoff for player 2 is 1
if player 2's k-subset contains player 1's chosen node.

# Notes

The actions of player 2 are ordered according to the
[combinatorial number system]
(https://en.wikipedia.org/wiki/Combinatorial_number_system),
which is different from the order used in the original library in C.

# Arguments

- `n::Integer` : Number of nodes in the tournament graph.
- `k::Integer` : Size of subsets of nodes in the tournament graph.
- `seed::Integer=-1`: Seed for random number generator. If seed is negative,
  then `GLOBAL_RNG` is used.

# Returns

- `g::NormalFormGame`

# Examples

```julia
julia> seed = 1234;

julia> g = tournament_game(5, 2; seed=seed)
5×10 NormalFormGame{2,Float64}

julia> g.players[1]
5×10 Player{2,Float64}:
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  1.0  1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> g.players[2]
10×5 Player{2,Float64}:
 1.0  1.0  0.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 0.0  1.0  1.0  0.0  0.0
 1.0  0.0  0.0  1.0  0.0
 0.0  1.0  0.0  1.0  0.0
 0.0  0.0  1.0  1.0  0.0
 1.0  0.0  0.0  0.0  1.0
 0.0  1.0  0.0  0.0  1.0
 0.0  0.0  1.0  0.0  1.0
 0.0  0.0  0.0  1.0  1.0
```
"""
function tournament_game(n::Integer, k::Integer; seed::Integer=-1)
    m = zero(Csize_t)
    try
        m = binomial(Csize_t(n), Csize_t(k))
    catch InexactError
        throw(ArgumentError("Maximum allowed size exceeded"))
    end

    R = zeros(Float64, n, m)
    C = zeros(Float64, m, n)
    tourn_graph = random_tournament_digraph(n; seed=seed)
    fadjlist = tourn_graph.fadjlist

    # populate matrix C
    X = collect(1:k)
    for j = 1:m
        C[j, X] .= 1.
        next_k_array!(X)
    end

    # populate matrix R
    # continue to use array `X` to store indices
    a = Vector{Int}(undef, k)
    for i = 1:n
        d = length(fadjlist[i])
        if d >= k
            for j = 1:k
                a[j] = j
            end
            while a[end] <= d
                for j in 1:k
                    X[j] = fadjlist[i][a[j]]
                end
                R[i, k_array_rank(Csize_t, X)] = 1.
                next_k_array!(a)
            end
        end
    end

    g = NormalFormGame([Player(R), Player(C)])
    return g
end

# unit_vector_game
"""
    unit_vector_game([rng=GLOBAL_RNG], n; avoid_pure_nash=false)

Return a NormalFormGame instance of the 2-player game "unit vector game" (Savani
and von Stengel, 2016). Payoffs for player 2 are chosen randomly from the [0, 1)
range. For player 1, each column contains exactly one 1 payoff and the rest is
0.

# Arguments

- `rng::AbstractRNG=GLOBAL_RNG`: Random number generator used.
- `n::Integer` : Number of actions.
- `avoid_pure_nash::Bool=false` : If true, player 1's payoffs will be placed in
  order to avoid pure Nash equilibria. (If necessary, the payoffs for player 2
  are redrawn so as not to have a dominant action.)

# Returns

- `g::NormalFormGame`

# Examples

```julia
julia> rng = MersenneTwister(1234);

julia> g = unit_vector_game(rng, 5)
5×5 NormalFormGame{2,Float64}

julia> g.players[1]
5×5 Player{2,Float64}:
 0.0  0.0  0.0  0.0  0.0
 1.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0

julia> g.players[2]
5×5 Player{2,Float64}:
 0.590845  0.854147  0.648882   0.112486   0.950498
 0.766797  0.200586  0.0109059  0.276021   0.96467
 0.566237  0.298614  0.066423   0.651664   0.945775
 0.460085  0.246837  0.956753   0.0566425  0.789904
 0.794026  0.579672  0.646691   0.842714   0.82116

julia> pure_nash(g)
1-element Array{Tuple{Int64,Int64},1}:
 (2, 1)
```

With `avoid_pure_nash=true`:

```julia
julia> rng = MersenneTwister(1234);

julia> g = unit_vector_game(rng, 5; avoid_pure_nash=true)
5×5 NormalFormGame{2,Float64}

julia> g.players[1]
5×5 Player{2,Float64}:
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  0.0  1.0
 0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0

julia> g.players[2]
5×5 Player{2,Float64}:
 0.590845  0.854147  0.648882   0.112486   0.950498
 0.766797  0.200586  0.0109059  0.276021   0.96467
 0.566237  0.298614  0.066423   0.651664   0.945775
 0.460085  0.246837  0.956753   0.0566425  0.789904
 0.794026  0.579672  0.646691   0.842714   0.82116

julia> pure_nash(g)
0-element Array{Tuple{Int64,Int64},1}
```
"""
function unit_vector_game(rng::AbstractRNG, n::Integer;
                          avoid_pure_nash::Bool=false)
    payoff_arrays = [zeros(Float64, n, n), rand(rng, Float64, (n, n))]

    if !avoid_pure_nash
        ones_ind = rand(rng, 1:n, n)
        for j in 1:n
            payoff_arrays[1][ones_ind[j], j] = 1
        end
    else
        n == 1 && throw(ArgumentError("Cannot avoid pure Nash with n=1"))
        maxes = maximum(payoff_arrays[2], dims=1)
        is_subotimal = payoff_arrays[2] .< maxes
        nums_suboptimal = sum(is_subotimal, dims=2)

        while any(nums_suboptimal .== 0)
            rand!(rng, payoff_arrays[2])
            maximum!(maxes, payoff_arrays[2])
            for i in 1:n, j in 1:n
                is_subotimal[j, i] = payoff_arrays[2][j, i] < maxes[i]
            end
            sum!(nums_suboptimal, is_subotimal)
        end

        for j in 1:n
            one_ind = rand(rng, 1:n)
            while !is_subotimal[j, one_ind]
                one_ind = rand(rng, 1:n)
            end
            payoff_arrays[1][one_ind, j] = 1
        end
    end

    g = NormalFormGame(
        [Player(payoff_array) for payoff_array::Array in payoff_arrays]
    )
    return g
end

unit_vector_game(n::Integer; avoid_pure_nash::Bool=false) =
    unit_vector_game(Random.GLOBAL_RNG, n, avoid_pure_nash=avoid_pure_nash)


text:
module FedeAlignment
include("clustering.jl")
include("dtw.jl")
include("fede_dtw.jl")
include("fede_dist.jl")
import Clustering: cutree
import Distances: pairwise!, Cityblock, PreMetric
using DelimitedFiles: readdlm
import PyCall: PyObject

"""
Given repeated notes, chose the the best one
"""
@views function _chose_sample(matched_idx::AbstractMatrix{Int64},
                              score_reverse::Vector{Int64}, perfm_reverse::Vector{Int64},
                              dist_mat::AbstractMatrix{Float64})::Array{Bool,1}

    to_keep = ones(Bool, size(matched_idx, 1))
    nonunique_indices = FedeDTW.my_nonunique(matched_idx[:, 1])
    for (val, indices) in nonunique_indices
        clusters1 = score_reverse[val]
        clusters2 = perfm_reverse[matched_idx[indices, 2]]
        # here clusters1 is only 1 number, corresponding to the cluster of the
        # element which is being repeated
        chosen = argmin(dist_mat[clusters1, clusters2])
        to_keep[indices] .= false
        to_keep[indices[chosen]] = true
    end
    return to_keep

end

"""
Computes a unique path from a path among notes by removing repeated notes.
It chooses the couple of notes associated with the two nearest clusters. In
case of parity, the first couple is taken.
"""
@views function _get_unique_path(matched_idx::Matrix{Int64}, score_reverse::Vector{Int64},
                                 perfm_reverse::Vector{Int64},
                                 dist_mat::Array{Float64,2})::Matrix{Int64}

    to_keep = _chose_sample(matched_idx, score_reverse, perfm_reverse, dist_mat)
    matched_idx = matched_idx[to_keep, :]
    to_keep = _chose_sample(matched_idx[:, end:-1:1],
                            # view(matched_idx, :, size(matched_idx, 2):-1:1),
                            perfm_reverse, score_reverse, transpose(dist_mat))
    matched_idx = matched_idx[to_keep, :]
    return matched_idx
end

"""
Compute distance matrix using `dist_func`
"""
function _compute_dist_mat(clusters1::Vector{Vector{Int64}},
                           clusters2::Vector{Vector{Int64}}, dist_func::PreMetric,
                           win_fn::Function)

    N = length(clusters1)
    M = length(clusters2)
    out = fill(Inf::Float64, N, M)

    for j in 1:M
        @simd for i in 1:N
            # limit the computation to radius
            # see slanted band in
            # https://github.com/DynamicTimeWarping/dtw-python/blob/master/dtw/window.py
            if win_fn(i, j)
                out[i, j] = dist_func(clusters1[i], clusters2[j])
            end
        end
    end
    # initial and last note have zero distance from any note
    out[1, :] .= 0
    out[end, :] .= 0
    out[:, 1] .= 0
    out[:, end] .= 0
    return out
end

struct ClusterResult
    features::Vector{Vector{Int64}}
    clusters::Vector{Vector{Int64}}
    reverse::Vector{Int64}
end

struct DistData{T<:AbstractMatrix,K<:Function}
    matrix::T
    window::K
    th::Float64
    simcarn::Bool
    weight::Float64
end

"""
Takes a list of cluster labels in [1, Inf) and transforms it in a list of lists of
indices/features; also returns a list of labels where labels are referred to the
transformed index of each cluster; this last list is similar to the input but
with sorted labels
"""
function _transform_clusters(clusters::Vector{Int64}, mat::Array{Float64,2})::ClusterResult
    # creating the list of clusters
    N = maximum(clusters)
    L = length(clusters)
    features = [Int64[] for _ in 1:N]
    transformed_clusters = [Int64[] for _ in 1:N]
    reverse = Vector{Int64}(undef, L)
    # a map from older clusters to new (ordered) number of clusters
    new_clusters = zeros(N)
    this_cluster = -1
    counter = 0
    new_cluster = -1
    for i in 1:L
        # for each note
        if clusters[i] != this_cluster
            this_cluster = clusters[i]
            new_cluster = new_clusters[this_cluster]
            if new_cluster == 0
                # a new cluster
                counter += 1
                new_cluster = counter
                new_clusters[this_cluster] = counter
            end
        end
        push!(features[new_cluster], floor(Int64, mat[i, 1]))
        push!(transformed_clusters[new_cluster], i)
        reverse[i] = new_cluster
    end
    return ClusterResult(features, transformed_clusters, reverse)
end

"""
Here, a cluster is a list of pitch. We take the list of clusters and
compare each note inside each cluster to match the corresponding notes. The
match is done based on the features (pitch), but the indices of the
corresponding notes are returned; this is why we need both the list of
cluster features and the list of cluster note indices.
"""
function _matching_notes_clusters(score_clusters, score_features, perfm_clusters,
                                  perfm_features, th::Float64, simcarn::Bool)::Matrix{Int64}
    out = Vector{Matrix{Int64}}(undef, 0)
    for i in 1:length(score_clusters)
        matched_idx1, matched_idx2 = FedeDist.my_intersect_idx(score_features[i],
                                                               perfm_features[i], th,
                                                               simcarn)

        if length(matched_idx1) > 0
            push!(out,
                  hcat(score_clusters[i][matched_idx1], perfm_clusters[i][matched_idx2]))
        end

    end
    return vcat(out...)
end

"""
Returns clusters of notes by onsets using single-linkage and stopping
agglomeration procedure when `threshold` is reached.

Returns the features of each cluster, ordered by onsets, the index of notes
in each cluster, ordered by onset, and the list of cluster of each note.
Features is only one (pitch).

If `num_clusters` is not None, it should be a int and `maxclust` criterion
is used in place of the threshold (which is not used)
"""
@views function _clusterize(mat::Array{Float64,2};
                            threshold::Union{Float64,Nothing}=nothing,
                            num_clusters::Union{Int64,Nothing}=nothing)::ClusterResult
    dist = Matrix{Float64}(undef, size(mat, 1), size(mat, 1))
    Z = Cluster.single_linkage!(pairwise!(dist, Cityblock(), mat[:, 2]))
    if threshold === nothing
        # clusterize the onsets using default linkage method in O(n²)
        clusters = cutree(Z; k=num_clusters)
    else
        clusters = cutree(Z; h=threshold)
    end

    # creating the list of clusters
    return _transform_clusters(clusters, mat)
end

function add_dummy_notes(matscore::Matrix{Float64}, matperfm::Matrix{Float64})
    # inserting starting and ending notes
    s::Float64 = -1.0
    e::Float64 = max(maximum(matscore[:, 3]), maximum(matperfm[:, 3])) + 1
    starting = fill(s, (1, size(matscore, 2)))
    ending = fill(e, (1, size(matscore, 2)))
    matscore = vcat(starting, matscore, ending)
    matperfm = vcat(starting, matperfm, ending)

    return matscore, matperfm
end

function _clustering(matscore, matperfm,
                     score_th::Union{Float64,Nothing}=0.0)::Tuple{ClusterResult,
                                                                  ClusterResult}
    # compute clusters
    perfm_clusters = _clusterize(matperfm; threshold=0.05)
    if score_th !== nothing
        score_clusters = _clusterize(matscore; threshold=score_th)
    else
        score_clusters = _clusterize(matscore; num_clusters=length(perfm_clusters.features))
    end
    return score_clusters, perfm_clusters

end

function _get_dist_data(score::ClusterResult, perfm::ClusterResult, th::Float64,
                        simcarn::Bool, α::Int64, β::Int64)::DistData

    dist_fn = FedeDist.JaccardDist(th, simcarn)
    window_fn = FedeDTW.get_fede_win(FedeDTW.FedeWindowParam(dist_fn, 5, α, β),
                                     score.features, perfm.features)

    dist_mat = _compute_dist_mat(score.features, perfm.features, dist_fn, window_fn)
    weight = 2 / 3
    if simcarn
        weight = 1 - weight
    end
    return DistData(dist_mat, window_fn, th, simcarn, weight)
end

"""
Returns a mapping of indices between notes in `matscore` and in `matperfm`
with shape (N, 2), where N is the number of matching notes.
Also returns `matscore` and `matperfm` with initial and ending virtual
notes: you can safely discard them by slicing with `[1:-1]`, but returned
indices are referred to the returned mats.

If `score_th` is None, then the number of clusters for the score is
inferred from the performance (whose threshold is fixed in `settings.py`)
"""
@views function _get_matching_notes(score::ClusterResult, perfm::ClusterResult,
                                    distdata::DistData,
                                    step_pattern::DTW.StepPattern)::Matrix{Int64}

    path = DTW.dtw_path(distdata.matrix; window_fn=distdata.window,
                        step_pattern=step_pattern)

    if path === nothing
        # # no path available
        # return an empty matrix
        return Matrix{Int64}(undef, 0, 0)
    end
    path1 = path[:, 1]
    path2 = path[:, 2]

    # computing the matched notes
    matched_idx = _matching_notes_clusters(score.clusters[path1], score.features[path1],
                                           perfm.clusters[path2], perfm.features[path2],
                                           distdata.th, distdata.simcarn)

    # compute unique matched notes
    matched_idx = _get_unique_path(matched_idx, score.reverse, perfm.reverse,
                                   distdata.matrix)
    return matched_idx
end

"""
Reapeatly calls ``_get_matching_notes`` and merge the indices.

Returns a mapping of indices between notes in `matscore` and in `matperfm`
with shape (N, 2), where N is the number of matching notes.
Also returns `matscore` and `matperfm` with initial and ending virtual
notes: you can safely discard them by slicing with `[1:-1]`, but returned
indices are referred to the returned mats.

If `score_th` is None, then the number of clusters for the score is
inferred from the performance (whose threshold is fixed in `settings.py`)

`thresholds` is a list of thresholds used; a value of 1 causes simcarn not
being used, a value != 1 causes simcarn to be used

`dist` is the function used for thresholds != 1 (using simcarn)

`sp_weights` must be a matrix where each row is a weight for the
`DTW.weight_step_pattern` function; it should have size (13, N), where N is the
number of desired different step-patterns DTW.

"""
function get_matching_notes(matscore::Array{Float64,2}, matperfm::Array{Float64,2},
                            α::Int64, β::Int64, sp_weights::Matrix{Float64},
                            score_th::Union{Float64,Nothing} = nothing,
                            thresholds::Vector{T}=[1.0, 0.5]) where {T<:Real}
    matscore, matperfm = add_dummy_notes(matscore, matperfm)

    # computing clustering once for all
    clusters::Tuple{ClusterResult,ClusterResult} = _clustering(matscore, matperfm, score_th)

    # computing distances once for each threshold
    distdata = Vector{DistData{Matrix{Float64},FedeDTW.FedeWindowFunc}}(undef,
                                                                        length(thresholds))
    # for (i, th) in collect(enumerate(thresholds))
    Threads.@threads for (i, th) in collect(enumerate(thresholds))
        distdata[i] = _get_dist_data(clusters[1], clusters[2], th, th != 1, α, β)
    end

    # compute step_patterns
    min_weight = sum(@view(sp_weights[end, :])) * 0.01
    step_patterns = [
        DTW.weight_step_pattern(col) for col in eachcol(sp_weights) if col[end] > min_weight
    ]

    config = [(dist, step, step.weight * dist.weight)
              for dist::DistData in distdata, step::DTW.StepPattern in step_patterns]
    matched_indices = Vector{Matrix{Int64}}(undef, length(config))

    # for (i, conf) in collect(enumerate(config))
    Threads.@threads for (i, conf) in collect(enumerate(config))
        matched_indices[i] = _get_matching_notes(clusters[1], clusters[2], conf[1], conf[2])
    end

    # merging matching indices
    matched_idx = FedeDTW.merge_matching_indices(matched_indices, [conf[3] for conf in config])
    return copy(matched_idx), matscore, matperfm
end

function get_matching_notes(matscore::Array{Float64,2}, matperfm::Array{Float64,2},
                            α::PyObject, β::PyObject,
                            sp_weights::Matrix{Float64},
                            score_th=Union{Float64,Nothing} = nothing,
                            thresholds::Vector{T}=[1.0, 0.5]) where {T<:Real}
    _α::Int64 = convert(Int64, α)
    _β::Int64 = convert(Int64, β)
    # _sp_weights = convert(Matrix{Float64}, sp_weights)
    return get_matching_notes(matscore, matperfm, _α, _β, sp_weights, score_th, thresholds)
end

"""
Just a simple function that you can use to precompile the core functions
"""
function precompile(path1, path2)
    matscore = readdlm(path1, ',', Float64)
    matperfm = readdlm(path2, ',', Float64)
    clusters::Tuple{ClusterResult,ClusterResult} = _clustering(matscore, matperfm, nothing)
    clusters = _clustering(matscore, matperfm, 0.05)
    step = DTW.get_default_step_patterns()[1]
    distdata = [_get_dist_data(clusters[1], clusters[2], 0.5, true, 5, 1),
                _get_dist_data(clusters[1], clusters[2], 1.0, false, 1, 1)]
    m1 = _get_matching_notes(clusters[1], clusters[2], distdata[1], step)
    m2 = _get_matching_notes(clusters[1], clusters[2], distdata[2], step)
    return matched_idx = FedeDTW.merge_matching_indices([m1, m2], [1.0, 1.0])
end
end # module


text:
# Plotting data along two X-axes, in original space
using Plots

""" Plots the feasible and infeasible samples of BlackBoxClassifier in two axes. """
function plot_2d(bbc::BlackBoxClassifier; axes = [1, 2])
    feas_idxs = findall(x -> x .>= 0, bbc.Y)
    infeas_idxs = findall(x -> x .< 0, bbc.Y)
    scatter(Matrix(bbc.X)[feas_idxs, axes[1]], Matrix(bbc.X)[feas_idxs, axes[2]],
            legend = false, color = :green,  fmt = :png)
    scatter!(Matrix(bbc.X)[infeas_idxs, axes[1]], Matrix(bbc.X)[infeas_idxs, axes[2]],
            title = "Feasible vs. infeasible samples",
    legend = false, color = :red,  fmt = :png)
end

""" Plots whether or not BlackBoxClassifier predictions are correct, in 2D. """
function plot_2d_predictions(bbc::BlackBoxClassifier; axes = [1, 2])
    Y = IAI.predict(bbc.learners[end], bbc.X)
    trues = [i for i=1:length(Y) if Y[i] == (bbc.Y[i] >= 0)];
    falses = [i for i=1:length(Y) if Y[i] != (bbc.Y[i] >= 0)];
    scatter(Matrix(bbc.X)[trues, axes[1]], Matrix(bbc.X)[trues, axes[2]], legend = false,
            color = :purple,  fmt = :png)
    scatter!(Matrix(bbc.X)[falses, axes[1]], Matrix(bbc.X)[falses, axes[2]], legend = false,
             title = "Correct vs. incorrect predictions",
             color = :yellow,  fmt = :png)
end

""" Plots the data in BlackBoxRegressor."""
function plot_2d(bbr::BlackBoxRegressor; axes = [1, 2])
    plt3d = Plots.plot(bbr.X[:, axes[1]], bbr.X[:, axes[2]], bbr.Y,
    seriestype=:surface, markersize = 2, camera=(20,30))
    xlabel!(string(bbr.vars[axes[1]]))
    ylabel!(string(bbr.vars[axes[2]]))
    title!("Data")
    display(plt3d)
end

""" Plots whether or not BlackBoxRegressor predictions are correct, in 2D. """
function plot_2d_predictions(bbr::BlackBoxRegressor; axes = [1, 2])
    plt3d2 = Plots.plot(bbr.X[:, axes[1]], bbr.X[:, axes[2]], bbr.predictions,
    seriestype=:surface, markersize = 2, camera=(20,30))
    xlabel!(string(bbr.vars[axes[1]]))
    ylabel!(string(bbr.vars[axes[2]]))
    title!("Predictions")
    display(plt3d2)
end

""" Plots the accuracies of different BlackBoxLearners. """
function plot_accuracies(bbls::Array{BlackBoxClassifier})
        bar(1:length(bbls), [bbl.accuracies[end] for bbl in bbls], xlabel="Constraint number",
                        xticks=1:length(ineqs), title="Constraint accuracies", legend=false)
end

# function

# colors = ["yellow", "orange", "red"]
# plt = plot()
# plt = quiver!([gm.solution_history[end, "x"]], [gm.solution_history[end,"y"]],
#     quiver = (obj_gradient[:, "x"]./sqrt.((obj_gradient[:, "x"].^2 .+ obj_gradient[:, "y"].^2)),
#                 obj_gradient[:, "y"]./sqrt.((obj_gradient[:, "x"].^2 .+ obj_gradient[:, "y"].^2))),
#     markershape = :star5, color = "green", label = "+ Objective")
# for i=1:length(bbls)
#     bbl = bbls[i]
#     infeas_idxs = findall(x -> x .< 0, bbl.Y)
#     plt = scatter!(bbl.X[infeas_idxs,"x"], bbl.X[infeas_idxs, "y"], color = colors[i])
#     plt = quiver!([gm.solution_history[end, "x"]], [gm.solution_history[end,"y"]],
#     quiver = (constr_grads[[i], "x"]./sqrt.((constr_grads[[i], "x"].^2 .+ constr_grads[[i], "y"].^2)),
#                 constr_grads[[i], "y"]./sqrt.((constr_grads[[i], "x"].^2 .+ constr_grads[[i], "y"].^2))),
#     markershape = :circle, color = colors[i], label = "+ BBL$(i)")
# end
# plt = quiver!([gm.solution_history[end, "x"]], [gm.solution_history[end,"y"]],
#     quiver = ([d_vals[1]], [d_vals[2]]),
#     markershape = :square, color = "purple", label = "Descent direction")

# using Plots
# plt = scatter(gm.solution_history[:,"x"], gm.solution_history[:,"y"])
# display(plt)