text:
module TestDirichlet
using Test
using EvidentialDeepLearning
const EDL = EvidentialDeepLearning
import Distributions; const DI = Distributions
using Distributions: logpdf, entropy, Uniform, var
using SpecialFunctions: loggamma, digamma

@testset "Against Distributions" begin
    for k in 1:10
        D111 = Dirichlet(Tuple(1 for _ in 1:k))
        @test EDL.nclasses(D111) === k
        p = rand(k); p = p / sum(p)
        @test logpdf(D111, p) ≈ loggamma(k)
        @test entropy(D111) ≈ -loggamma(k)

        @test EDL.kl_uniform(D111) ≈ 0 atol=100eps(Float64)
    end

    for i in 1:10
        k = rand(1:10)
        d_edl  = Dirichlet(Tuple(10*rand(k)))
        d_dist = convert(DI.Distribution, d_edl)
        @test d_dist isa DI.Dirichlet
        p = rand(k)
        p = p / sum(p)
        @test logpdf(d_dist, p) ≈ logpdf(d_edl, p)
        @test entropy(d_dist) ≈ entropy(d_edl)
        @test var(d_dist) ≈ var(d_edl)
        if EDL.nclasses(d_edl) == 1
            @test EDL.kl_uniform(d_edl) == 0
        else
            @test EDL.kl_uniform(d_edl) > 0
        end
    end
end

@testset "against evidential-deep-learning" begin
    function _Dirichlet_SOS_KL_devec(d::EDL.Dirichlet)
        # Line by line port of
        # https://github.com/aamini/evidential-deep-learning/blob/d1d8e395fb083308d14fa92c5ce766e97b2a066a/evidential_deep_learning/losses/discrete.py#L6
        #
        α = d.α
        # beta=tf.constant(np.ones((1,alpha.shape[1])),dtype=tf.float32)
        β   = one.(α)
        # S_alpha = tf.reduce_sum(alpha,axis=1,keepdims=True)
        S_α = sum(α)
        # S_beta = tf.reduce_sum(beta,axis=1,keepdims=True)
        S_β = sum(β)
        # lnB = tf.math.lgamma(S_alpha) - tf.reduce_sum(tf.math.lgamma(alpha),axis=1,keepdims=True)
        lnB = loggamma(S_α) - sum(loggamma, α)
        # lnB_uni = tf.reduce_sum(tf.math.lgamma(beta),axis=1,keepdims=True) - tf.math.lgamma(S_beta)
        # lnB_uni = tf.reduce_sum(tf.math.lgamma(beta),axis=1,keepdims=True) - tf.math.lgamma(S_beta)
        lnB_uni = sum(loggamma, β) - loggamma(S_β)
        # dg0 = tf.math.digamma(S_alpha)
        dg0 = digamma(S_α)
        # dg1 = tf.math.digamma(alpha)
        dg1 = digamma.(α)
        # kl = tf.reduce_sum((alpha - beta)*(dg1-dg0),axis=1,keepdims=True) + lnB + lnB_uni
        kl = sum((α .- β).*(dg1 .- dg0)) + lnB + lnB_uni
        # return kl
        return kl
    end
    for _ in 1:100
        k = rand(1:10)
        α = Tuple(rand(Uniform(0,19), k))
        d = Dirichlet(α)
        @test _Dirichlet_SOS_KL_devec(d) ≈ EDL.kl_uniform(d)
    end

end

@testset "remove_non_misleading_evidence" begin
    EDL.remove_non_misleading_evidence(Dirichlet(1,2), 2) === Dirichlet(1,1)
    EDL.remove_non_misleading_evidence(Dirichlet(1f0,0.1f0, 3f0), 2) === Dirichlet(1f0,1f0,3f0)
end

end#module


text:
export FewShotDataLoader, MetaDataSample

using Serialization
using Random: shuffle!
import StatsBase
import Base

struct FewShotDataLoader
    labels::AbstractArray{L,1} where {L}
    samples::AbstractArray
    label2Idices::Dict

    function FewShotDataLoader(path::String)
        data = deserialize(path)
        @assert typeof(data) <: Dict "$path should containing a Dict, get $(typeof(data))"
        labels = data["labels"]
        samples = data["data"]
        @assert ndims(samples) >= 2 "ndims of samples should be at least 2, get $(ndims(samples))"
        N = size(labels)[end]
        label2Idices = Dict{eltype(labels),Array{Integer,1}}()
        for i = 1:N
            label = labels[i]
            if !haskey(label2Idices, label)
                label2Idices[label] = []
            end
            push!(label2Idices[label], i)
        end
        new(labels, samples, label2Idices)
    end
end

function Base.show(io::IO, dloader::FewShotDataLoader)
    print(
        io,
        "FewShotDataLoader(sample:$(size(dloader.samples)), labels:$(typeof(dloader.labels)))",
    )
end

"""
Fields
======

`support_samples`: array of support samples with size (`feature_size`..., `num_support`, `num_tasks`)
`support_labels`: array of labels of support samples with size (`num_support`, `num_tasks`)
`query_samples`: array of query samples with size (`feature_size`..., `num_query`, `num_tasks`)
`query_labels`: array of labels of query samples with size (`num_query`, `num_tasks`)
"""
struct MetaDataSample{T,L}
    support_n_ways::Integer
    support_k_shots::Integer
    query_n_ways::Integer
    query_k_shots::Integer
    support_samples::AbstractArray{T}
    support_labels::AbstractArray{L,2}
    query_samples::AbstractArray{T}
    query_labels::AbstractArray{L,2}
end

function Base.getproperty(m::MetaDataSample, name::Symbol)
    if name == :n_support
        return m.support_k_shots * m.support_n_ways
    elseif name == :n_query
        return m.query_k_shots * m.query_n_ways
    end
    getfield(m, name)
end

function Base.show(io::IO, meta_sample::MetaDataSample)
    print(
        io,
        "MetaDataSample(",
        "support_n_ways: $(meta_sample.support_n_ways), ",
        "support_k_shots: $(meta_sample.support_k_shots), ",
        "query_n_ways: $(meta_sample.query_n_ways), ",
        "query_k_shots: $(meta_sample.query_k_shots), ",
        "support_samples: [$(eltype(meta_sample.support_samples))]$(size(meta_sample.support_samples)), ",
        "support_labels: [$(eltype(meta_sample.support_labels))]$(size(meta_sample.support_labels)), ",
        "query_samples: [$(eltype(meta_sample.query_samples))]$(size(meta_sample.query_samples)), ",
        "query_labels: [$(eltype(meta_sample.query_labels))]$(size(meta_sample.query_labels)))",
    )
end

function MetaDataSample(;
    support_n_ways,
    support_k_shots,
    query_n_ways,
    query_k_shots,
    support_samples,
    support_labels,
    query_samples,
    query_labels,
)
    @assert support_n_ways > 0 "support_n_ways must be positive: $support_n_ways"
    @assert support_k_shots > 0 "support_k_shots must be positive: $support_k_shots"
    @assert query_n_ways > 0 "query_n_ways must be positive: $query_n_ways"
    @assert query_k_shots > 0 "query_k_shots must be positive: $query_k_shots"
    @assert ndims(support_samples) >= 3 "support_samples must be at least 3-dim"
    @assert ndims(query_samples) >= 3 "query_samples must be at least 3-dim"
    @assert support_n_ways * support_k_shots == size(support_samples)[end-1]
    @assert query_n_ways * query_k_shots == size(query_samples)[end-1]
    @assert size(support_samples)[1:end-2] == size(query_samples)[1:end-2] "size of support samples and query samples should be the same"
    MetaDataSample(
        support_n_ways,
        support_k_shots,
        query_n_ways,
        query_k_shots,
        support_samples,
        support_labels,
        query_samples,
        query_labels,
    )
end

function Base.size(meta_data_sample::MetaDataSample)
    return size(meta_data_sample.support_samples)[end]
end

function _sample(
    dloader::FewShotDataLoader;
    support_n_ways = 2,
    support_k_shots = 5,
    query_k_shots = support_k_shots,
)
    # find exclusive train/test samples
    uniq_labels = unique(dloader.labels)
    target_labels = StatsBase.sample(uniq_labels, support_n_ways, replace = false)
    support_idxs = []
    query_idxs = []
    for label in target_labels
        idx_map = Dict([(v, i) for (i, v) in enumerate(dloader.label2Idices[label])])
        selected_idxs = StatsBase.sample(dloader.label2Idices[label], support_k_shots)
        support_idxs = vcat(
            support_idxs,
            selected_idxs,
        )
        weights = ones(length(dloader.label2Idices[label]))
        for idx in selected_idxs
            weights[idx_map[idx]] = 0
        end
        weights = StatsBase.Weights(weights)
        selected_idxs =
            StatsBase.sample(dloader.label2Idices[label], weights, query_k_shots)
        query_idxs = vcat(
            query_idxs,
            selected_idxs
        )
    end
    @assert intersect(support_idxs, query_idxs) |> isempty
    shuffle!(support_idxs)
    shuffle!(query_idxs)
    support_samples = Utils.add_dim(
        dloader.samples[repeat([:], ndims(dloader.samples)-1)..., support_idxs]
    )
    support_labels = Utils.add_dim(dloader.labels[support_idxs])
    query_samples = Utils.add_dim(
        dloader.samples[repeat([:], ndims(dloader.samples)-1)..., query_idxs]
    )
    query_labels = Utils.add_dim(dloader.labels[query_idxs])
    return support_samples, support_labels, query_samples, query_labels
end

function StatsBase.sample(
    dloader::FewShotDataLoader;
    support_n_ways = 2,
    support_k_shots = 5,
    query_k_shots = support_k_shots,
)
    support_samples, support_labels, query_samples, query_labels = _sample(
        dloader,
        support_n_ways = support_n_ways,
        support_k_shots = support_k_shots,
        query_k_shots = query_k_shots,
    )
    MetaDataSample(
        support_n_ways = support_n_ways,
        support_k_shots = support_k_shots,
        query_n_ways = support_n_ways,
        query_k_shots = query_k_shots,
        support_samples = support_samples,
        support_labels = support_labels,
        query_samples = query_samples,
        query_labels = query_labels,
    )
end

function StatsBase.sample(
    dloader::FewShotDataLoader,
    n_tasks::Integer;
    support_n_ways = 2,
    support_k_shots = 5,
    query_k_shots = support_k_shots,
)
    support_sets = []
    support_labels = []
    query_sets = []
    query_labels = []
    for i = 1:n_tasks
        support, support_label, query, query_label = _sample(
            dloader,
            support_n_ways = support_n_ways,
            support_k_shots = support_k_shots,
            query_k_shots = query_k_shots,
        )
        push!(support_sets, support)
        push!(support_labels, support_label)
        push!(query_sets, query)
        push!(query_labels, query_label)
    end
    ndim_sample = ndims(support_sets[1])
    MetaDataSample(
        support_n_ways = support_n_ways,
        support_k_shots = support_k_shots,
        query_n_ways = support_n_ways,
        query_k_shots = query_k_shots,
        support_samples = cat(support_sets..., dims = Val(ndim_sample)),
        support_labels = cat(support_labels..., dims = Val(2)),
        query_samples = cat(query_sets..., dims = Val(ndim_sample)),
        query_labels = cat(query_labels..., dims = Val(2)),
    )
end

ext: jl


text:
@testset "Test Reserves from Thermal Dispatch" begin
    template = get_thermal_dispatch_template_network(CopperPlatePowerModel)
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve1"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve11"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveDown}, RangeReserve, "Reserve2"),
    )
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve, "ORDC1"),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves=true)
    model = DecisionModel(template, c_sys5_uc)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 648, 0, 120, 216, 72, false)
    reserve_variables = [
        :ActivePowerReserveVariable__VariableReserve_ReserveUp_Reserve1
        :ActivePowerReserveVariable__ReserveDemandCurve_ReserveUp_ORDC1
        :ActivePowerReserveVariable__VariableReserve_ReserveDown_Reserve2
        :ActivePowerReserveVariable__VariableReserve_ReserveUp_Reserve11
    ]
    for (k, var_array) in model.internal.container.variables
        if PSI.encode_key(k) in reserve_variables
            for var in var_array
                @test JuMP.has_lower_bound(var)
                @test JuMP.lower_bound(var) == 0.0
            end
        end
    end
end

@testset "Test Ramp Reserves from Thermal Dispatch" begin
    template = get_thermal_dispatch_template_network(CopperPlatePowerModel)
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RampReserve, "Reserve1"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RampReserve, "Reserve11"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveDown}, RampReserve, "Reserve2"),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves=true)
    model = DecisionModel(template, c_sys5_uc)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 384, 0, 336, 192, 24, false)
    reserve_variables = [
        :ActivePowerReserveVariable__VariableReserve_ReserveDown_Reserve2,
        :ActivePowerReserveVariable__VariableReserve_ReserveUp_Reserve1,
        :ActivePowerReserveVariable__VariableReserve_ReserveUp_Reserve11,
    ]
    for (k, var_array) in model.internal.container.variables
        if PSI.encode_key(k) in reserve_variables
            for var in var_array
                @test JuMP.has_lower_bound(var)
                @test JuMP.lower_bound(var) == 0.0
            end
        end
    end
end

@testset "Test Reserves from Thermal Standard UC" begin
    template = get_thermal_standard_uc_template()
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve1"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve11"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveDown}, RangeReserve, "Reserve2"),
    )
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve, "ORDC1"),
    )
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves=true)

    model = DecisionModel(template, c_sys5_uc; optimizer=cbc_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 1008, 0, 480, 216, 192, true)
end

@testset "Test Reserves from Thermal Standard UC with NonSpinningReserve" begin
    template = get_thermal_standard_uc_template()
    set_device_model!(
        template,
        DeviceModel(ThermalMultiStart, ThermalStandardUnitCommitment),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserveNonSpinning, NonSpinningReserve, "NonSpinningReserve"),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc_non_spin"; add_reserves=true)
    model = DecisionModel(template, c_sys5_uc; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 1032, 0, 888, 192, 288, true)
end

@testset "Test Upwards Reserves from Renewable Dispatch" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve3"),
    )
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve, "ORDC1"),
    )

    c_sys5_re = PSB.build_system(PSITestSystems, "c_sys5_re"; add_reserves=true)
    model = DecisionModel(template, c_sys5_re)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 360, 0, 72, 48, 72, false)
end

@testset "Test Reserves from Storage" begin
    template = get_thermal_dispatch_template_network(CopperPlatePowerModel)
    set_device_model!(template, DeviceModel(GenericBattery, BatteryAncillaryServices))
    set_device_model!(template, RenewableDispatch, FixedOutput)
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve3"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveDown}, RangeReserve, "Reserve4"),
    )
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve, "ORDC1"),
    )

    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat"; add_reserves=true)
    model = DecisionModel(template, c_sys5_bat)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 432, 0, 288, 264, 96, true)
end

@testset "Test Reserves from Hydro" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchRunOfRiver)
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve5"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveDown}, RangeReserve, "Reserve6"),
    )
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve, "ORDC1"),
    )

    c_sys5_hyd = PSB.build_system(PSITestSystems, "c_sys5_hyd"; add_reserves=true)
    model = DecisionModel(template, c_sys5_hyd)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 240, 0, 48, 96, 72, false)
end

@testset "Test Reserves from with slack variables" begin
    template = get_thermal_dispatch_template_network(
        NetworkModel(CopperPlatePowerModel; use_slacks=true),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve1"; use_slacks=true),
    )
    set_service_model!(
        template,
        ServiceModel(
            VariableReserve{ReserveUp},
            RangeReserve,
            "Reserve11";
            use_slacks=true,
        ),
    )
    set_service_model!(
        template,
        ServiceModel(
            VariableReserve{ReserveDown},
            RangeReserve,
            "Reserve2";
            use_slacks=true,
        ),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves=true)
    model = DecisionModel(template, c_sys5_uc;)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 504, 0, 120, 192, 24, false)
end

@testset "Test AGC" begin
    c_sys5_reg = PSB.build_system(PSITestSystems, "c_sys5_reg")
    @test_throws ArgumentError template_agc_reserve_deployment(; dummy_arg=0.0)

    template_agc = template_agc_reserve_deployment()
    set_service_model!(template_agc, ServiceModel(PSY.AGC, PIDSmoothACE, "AGC_Area1"))
    agc_problem = DecisionModel(AGCReserveDeployment, template_agc, c_sys5_reg)
    @test build!(agc_problem; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    # These values might change as the AGC model is refined
    moi_tests(agc_problem, false, 696, 0, 480, 0, 384, false)
end

@testset "Test GroupReserve from Thermal Dispatch" begin
    template = get_thermal_dispatch_template_network()
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve1"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve11"),
    )
    set_service_model!(
        template,
        ServiceModel(VariableReserve{ReserveDown}, RangeReserve, "Reserve2"),
    )
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve, "ORDC1"),
    )
    set_service_model!(
        template,
        ServiceModel(StaticReserveGroup{ReserveDown}, GroupReserve, "init"),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves=true)
    services = get_components(Service, c_sys5_uc)
    contributing_services = Vector{Service}()
    for service in services
        if !(typeof(service) <: PSY.ReserveDemandCurve)
            push!(contributing_services, service)
        end
    end
    groupservice = StaticReserveGroup{ReserveDown}(;
        name="init",
        available=true,
        requirement=0.0,
        ext=Dict{String, Any}(),
    )
    add_service!(c_sys5_uc, groupservice, contributing_services)

    model = DecisionModel(template, c_sys5_uc)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 648, 0, 120, 240, 72, false)
end

@testset "Test GroupReserve Errors" begin
    template = get_thermal_dispatch_template_network()
    set_service_model!(template, ServiceModel(VariableReserve{ReserveUp}, RangeReserve))
    set_service_model!(template, ServiceModel(VariableReserve{ReserveDown}, RangeReserve))
    set_service_model!(
        template,
        ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )
    set_service_model!(
        template,
        ServiceModel(StaticReserveGroup{ReserveDown}, GroupReserve),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves=true)
    services = get_components(Service, c_sys5_uc)
    contributing_services = Vector{Service}()
    for service in services
        if !(typeof(service) <: PSY.ReserveDemandCurve)
            push!(contributing_services, service)
        end
    end
    groupservice = StaticReserveGroup{ReserveDown}(;
        name="init",
        available=true,
        requirement=0.0,
        ext=Dict{String, Any}(),
    )
    add_service!(c_sys5_uc, groupservice, contributing_services)

    off_service = VariableReserve{ReserveUp}("Reserveoff", true, 0.6, 10)
    push!(groupservice.contributing_services, off_service)

    model = DecisionModel(template, c_sys5_uc)
    @test build!(
        model;
        output_dir=mktempdir(cleanup=true),
        console_level=Logging.AboveMaxLevel,
    ) == BuildStatus.FAILED
end

@testset "Test StaticReserve" begin
    template = get_thermal_dispatch_template_network()
    set_service_model!(
        template,
        ServiceModel(StaticReserve{ReserveUp}, RangeReserve, "Reserve3"),
    )

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    static_reserve = StaticReserve{ReserveUp}("Reserve3", true, 30, 100)
    add_service!(c_sys5_uc, static_reserve, get_components(ThermalGen, c_sys5_uc))
    model = DecisionModel(template, c_sys5_uc)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    @test typeof(model) <: DecisionModel{<:PSI.DecisionProblem}
end

@testset "Test Reserves with Feedforwards" begin
    template = get_thermal_dispatch_template_network()
    service_model = ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve1")
    ff_lb = LowerBoundFeedforward(
        component_type=VariableReserve{ReserveUp},
        source=ActivePowerReserveVariable,
        affected_values=[ActivePowerReserveVariable],
        meta="Reserve1",
    )
    PSI.attach_feedforward!(service_model, ff_lb)

    set_service_model!(template, service_model)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves=true)
    model = DecisionModel(template, c_sys5_uc; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 240, 0, 120, 264, 24, false)
end


text:
# Given test reports from unit test runs, select a subset of tests.
using EzXML

"""
    test_outcomes_from_xml(xmlpath)

Reads the XML from TestReports, which is a standard JUnit XML format.
The xmlpath is a document object model (DOM) from EzXML's `readxml()`.
Count each test suite as a unit test, not individual tests within
the test suite. Return a dictionary from test name to number of failures.
"""
function test_outcomes_from_xml(xmlpath)
    outcomes = Dict{String, Int}()
    for testsuite in findall("//testsuite", xmlpath)
        name = nothing
        failures = nothing
        for attrib in attributes(testsuite)
            if attrib.name == "name"
                name = attrib.content
            elseif attrib.name == "failures"
                failures = parse(Int, attrib.content)
            end
        end
        if !isnothing(name) && !isnothing(failures)
            outcomes[name] = failures
        end
    end
    outcomes
end

"""
    read_reports(dir)

Read unit test reports in a directory as a matrix of which passed and
failed for each run. This assumes every XML file in the directory is
a unit test log report.
"""
function read_reports(dir)
    reports = filter(x -> startswith(x, "log") && endswith(x, ".xml"), readdir(dir))
    sample = test_outcomes_from_xml(EzXML.readxml(joinpath(dir, first(reports))))
    println("sample $(sample)")
    outcomes = zeros(Bool, length(sample), length(reports))
    key = Dict{String, Int}((b, a) for (a, b) in enumerate(keys(sample)))

    for (report_idx, report) in enumerate(reports)
        log = EzXML.readxml(joinpath(dir, report))
        outcome = test_outcomes_from_xml(log)
        for (test_name, test_fails) in outcome
            if !haskey(key, test_name)
                # In case a test case was added mid-way through.
                key[test_name] = length(key) + 1
                next_outcomes = zeros(Bool, length(key), length(reports))
                next_outcomes[1:(length(key) - 1), :] .= outcomes
                outcomes = next_outcomes
            end
            outcomes[key[test_name], report_idx] = test_fails > 0
        end
    end
    (outcomes, key)
end

"""
    select_tests(outcomes)

Given a Bool array of pass/fail outcomes that is (test cases) x (test runs),
choose a subset of test cases that would have found every failure.

This is a greedy algorithm. It begins by choosing the test case that
caught the most failed tests. Then it forgets about those tests and
chooses the next test case that eliminated the most failed tests among the
remaining.
"""
function select_tests(outcomes)
    detected = outcomes[:, vec(sum(outcomes, dims = 1) .> 0)]
    choices = zeros(Int, size(outcomes, 1))
    choice_cnt = 0
    while sum(detected) > 0
        choice_cnt += 1
        choice = argmax(vec(sum(detected, dims = 2)))
        choices[choice_cnt] = choice
        detected = detected[:, .!detected[choice, :]]
    end
    choices[1:choice_cnt]
end


text:
#
# Xlib.jl --
#
# Julia wrapper of X11 library.
#
#------------------------------------------------------------------------------
#
# Copyright (C) 2017, Éric Thiébaut.
#
# This file is part of Xlib.jl which is licensed under the MIT "Expat" License.
#

__precompile__(true)

module Xlib

using Compat

const _XLIB = "libX11.so" # FIXME: should be set in ../deps

include("exports.jl")
include("constants.jl")
include("types.jl")
include("methods.jl")
include("utils.jl")

end


text:
# An instruction builder represents a point within a basic block and is the exclusive means
# of building instructions using the C interface.

export Builder,
       position!,
       debuglocation, debuglocation!

@checked struct Builder
    ref::API.LLVMBuilderRef
end

Base.unsafe_convert(::Type{API.LLVMBuilderRef}, builder::Builder) = builder.ref

Builder(ctx::Context) = Builder(API.LLVMCreateBuilderInContext(ctx))

dispose(builder::Builder) = API.LLVMDisposeBuilder(builder)

function Builder(f::Core.Function, args...; kwargs...)
    builder = Builder(args...; kwargs...)
    try
        f(builder)
    finally
        dispose(builder)
    end
end

Base.position(builder::Builder) = BasicBlock(API.LLVMGetInsertBlock(builder))
position!(builder::Builder, inst::Instruction) =
    API.LLVMPositionBuilderBefore(builder, inst)
position!(builder::Builder, bb::BasicBlock) =
    API.LLVMPositionBuilderAtEnd(builder, bb)
position!(builder::Builder) = API.LLVMClearInsertionPosition(builder)

Base.insert!(builder::Builder, inst::Instruction) =
    API.LLVMInsertIntoBuilder(builder, inst)
Base.insert!(builder::Builder, inst::Instruction, name::String) =
    API.LLVMInsertIntoBuilderWithName(builder, inst, name)

function debuglocation(builder::Builder)
    ref = API.LLVMGetCurrentDebugLocation2(builder)
    ref == C_NULL ? nothing : Metadata(ref)
end
debuglocation!(builder::Builder) =
    API.LLVMSetCurrentDebugLocation2(builder, C_NULL)
debuglocation!(builder::Builder, loc::Metadata) =
    API.LLVMSetCurrentDebugLocation2(builder, loc)
debuglocation!(builder::Builder, loc::MetadataAsValue) =
    API.LLVMSetCurrentDebugLocation2(builder, Metadata(loc))
debuglocation!(builder::Builder, inst::Instruction) =
    API.LLVMSetInstDebugLocation(builder, inst)

## build methods

# TODO/IDEAS:
# - dynamic dispatch based on `llvmtype` (eg. disambiguating `add!` and `fadd!`)

# NOTE: the return values for these operations are, according to the C API, always a Value.
#       however, the C++ API learns us that we can be more strict.

export ret!, br!, switch!, indirectbr!, invoke!, resume!, unreachable!,

       binop!, add!, nswadd!, nuwadd!, fadd!, sub!, nswsub!, nuwsub!, fsub!, mul!, nswmul!,
       nuwmul!, fmul!, udiv!, sdiv!, exactsdiv!, fdiv!, urem!, srem!, frem!, neg!, nswneg!,
       nuwneg!, fneg!,

       shl!, lshr!, ashr!, and!, or!, xor!, not!,

       extract_element!, insert_element!, shuffle_vector!,

       extract_value!, insert_value!,

       alloca!, array_alloca!, malloc!, array_malloc!, free!, load!, store!, fence!,
       atomic_rmw!, atomic_cmpxchg!, gep!, inbounds_gep!, struct_gep!,

       trunc!, zext!, sext!, fptoui!, fptosi!, uitofp!, sitofp!, fptrunc!, fpext!,
       ptrtoint!, inttoptr!, bitcast!, addrspacecast!, zextorbitcast!, sextorbitcast!,
       truncorbitcast!, cast!, pointercast!, intcast!, fpcast!,

       icmp!, fcmp!, phi!, select!, call!, va_arg!, landingpad!,

       globalstring!, globalstring_ptr!, isnull!, isnotnull!, ptrdiff!

# terminator instructions

ret!(builder::Builder) =
    Instruction(API.LLVMBuildRetVoid(builder))

ret!(builder::Builder, V::Value) =
    Instruction(API.LLVMBuildRet(builder, V))

ret!(builder::Builder, RetVals::Vector{<:Value}) =
    Instruction(API.LLVMBuildAggregateRet(builder, RetVals, length(RetVals)))

br!(builder::Builder, Dest::BasicBlock) =
    Instruction(API.LLVMBuildBr(builder, Dest))

br!(builder::Builder, If::Value, Then::BasicBlock, Else::BasicBlock) =
    Instruction(API.LLVMBuildCondBr(builder, If, Then, Else))

switch!(builder::Builder, V::Value, Else::BasicBlock, NumCases::Integer=10) =
    Instruction(API.LLVMBuildSwitch(builder, V, Else, NumCases))

indirectbr!(builder::Builder, Addr::Value, NumDests::Integer=10) =
    Instruction(API.LLVMBuildIndirectBr(builder, Addr, NumDests))

invoke!(builder::Builder, Fn::Value, Args::Vector{<:Value}, Then::BasicBlock, Catch::BasicBlock, Name::String="") =
    Instruction(API.LLVMBuildInvoke(builder, Fn, Args, length(Args), Then, Catch, Name))

resume!(builder::Builder, Exn::Value) =
    Instruction(API.LLVMBuildResume(builder, Exn))

unreachable!(builder::Builder) =
    Instruction(API.LLVMBuildUnreachable(builder))

# binary operations

binop!(builder::Builder, Op::API.LLVMOpcode, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildBinOp(builder, Op, LHS, RHS, Name))

add!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildAdd(builder, LHS, RHS, Name))

nswadd!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNSWAdd(builder, LHS, RHS, Name))

nuwadd!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNUWAdd(builder, LHS, RHS, Name))

fadd!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFAdd(builder, LHS, RHS, Name))

sub!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildSub(builder, LHS, RHS, Name))

nswsub!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNSWSub(builder, LHS, RHS, Name))

nuwsub!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNUWSub(builder, LHS, RHS, Name))

fsub!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFSub(builder, LHS, RHS, Name))

mul!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildMul(builder, LHS, RHS, Name))

nswmul!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNSWMul(builder, LHS, RHS, Name))

nuwmul!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNUWMul(builder, LHS, RHS, Name))

fmul!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFMul(builder, LHS, RHS, Name))

udiv!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildUDiv(builder, LHS, RHS, Name))

sdiv!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildSDiv(builder, LHS, RHS, Name))

exactsdiv!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildExactSDiv(builder, LHS, RHS, Name))

fdiv!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFDiv(builder, LHS, RHS, Name))

urem!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildURem(builder, LHS, RHS, Name))

srem!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildSRem(builder, LHS, RHS, Name))

frem!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFRem(builder, LHS, RHS, Name))

# bitwise binary operations

shl!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildShl(builder, LHS, RHS, Name))

lshr!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildLShr(builder, LHS, RHS, Name))

ashr!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildAShr(builder, LHS, RHS, Name))

and!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildAnd(builder, LHS, RHS, Name))

or!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildOr(builder, LHS, RHS, Name))

xor!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildXor(builder, LHS, RHS, Name))

# vector operations

extract_element!(builder::Builder, VecVal::Value, Index::Value, Name::String="") =
    Value(API.LLVMBuildExtractElement(builder, VecVal, Index, Name))

insert_element!(builder::Builder, VecVal::Value, EltVal::Value, Index::Value, Name::String="") =
    Value(API.LLVMBuildInsertElement(builder, VecVal, EltVal, Index, Name))

shuffle_vector!(builder::Builder, V1::Value, V2::Value, Mask::Value, Name::String="") =
    Value(API.LLVMBuildShuffleVector(builder, V1, V2, Mask, Name))

# aggregate operations

extract_value!(builder::Builder, AggVal::Value, Index, Name::String="") =
    Value(API.LLVMBuildExtractValue(builder, AggVal, Index, Name))

insert_value!(builder::Builder, AggVal::Value, EltVal::Value, Index, Name::String="") =
    Value(API.LLVMBuildInsertValue(builder, AggVal, EltVal, Index, Name))

# memory access and addressing operations

alloca!(builder::Builder, Ty::LLVMType, Name::String="") =
    Instruction(API.LLVMBuildAlloca(builder, Ty, Name))

array_alloca!(builder::Builder, Ty::LLVMType, Val::Value, Name::String="") =
    Instruction(API.LLVMBuildArrayAlloca(builder, Ty, Val, Name))

malloc!(builder::Builder, Ty::LLVMType, Name::String="") =
    Instruction(API.LLVMBuildMalloc(builder, Ty, Name))

array_malloc!(builder::Builder, Ty::LLVMType, Val::Value, Name::String="") =
    Instruction(API.LLVMBuildArrayMalloc(builder, Ty, Val, Name))

free!(builder::Builder, PointerVal::Value) =
    Instruction(API.LLVMBuildFree(builder, PointerVal))

load!(builder::Builder, PointerVal::Value, Name::String="") =
    Instruction(API.LLVMBuildLoad(builder, PointerVal, Name))

store!(builder::Builder, Val::Value, Ptr::Value) =
    Instruction(API.LLVMBuildStore(builder, Val, Ptr))

fence!(builder::Builder, ordering::API.LLVMAtomicOrdering, singleThread::Core.Bool=false, Name::String="") =
    Instruction(API.LLVMBuildFence(builder, ordering, convert(Bool, singleThread), Name))
Ptr
atomic_rmw!(builder::Builder, op::API.LLVMAtomicRMWBinOp, Ptr::Value, Val::Value, ordering::API.LLVMAtomicOrdering, singleThread::Core.Bool) =
    Instruction(API.LLVMBuildAtomicRMW(builder, op, Ptr, Val, ordering, convert(Bool, singleThread)))

atomic_cmpxchg!(builder::Builder, Ptr::Value, Cmp::Value, New::Value, SuccessOrdering::API.LLVMAtomicOrdering, FailureOrdering::API.LLVMAtomicOrdering, SingleThread::Core.Bool) =
    Instruction(API.LLVMBuildAtomicCmpXchg(builder, Ptr, Cmp, New, SuccessOrdering,FailureOrdering, convert(Bool, SingleThread)))

gep!(builder::Builder, Pointer::Value, Indices::Vector{<:Value}, Name::String="") =
    Value(API.LLVMBuildGEP(builder, Pointer, Indices, length(Indices), Name))

inbounds_gep!(builder::Builder, Pointer::Value, Indices::Vector{<:Value}, Name::String="") =
    Value(API.LLVMBuildInBoundsGEP(builder, Pointer, Indices, length(Indices), Name))

struct_gep!(builder::Builder, Pointer::Value, Idx, Name::String="") =
    Value(API.LLVMBuildStructGEP(builder, Pointer, Idx, Name))

# conversion operations

trunc!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildTrunc(builder, Val, DestTy, Name))

zext!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildZExt(builder, Val, DestTy, Name))

sext!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildSExt(builder, Val, DestTy, Name))

fptoui!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPToUI(builder, Val, DestTy, Name))

fptosi!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPToSI(builder, Val, DestTy, Name))

uitofp!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildUIToFP(builder, Val, DestTy, Name))

sitofp!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildSIToFP(builder, Val, DestTy, Name))

fptrunc!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPTrunc(builder, Val, DestTy, Name))

fpext!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPExt(builder, Val, DestTy, Name))

ptrtoint!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildPtrToInt(builder, Val, DestTy, Name))

inttoptr!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildIntToPtr(builder, Val, DestTy, Name))

bitcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildBitCast(builder, Val, DestTy, Name))

addrspacecast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildAddrSpaceCast(builder, Val, DestTy, Name))

zextorbitcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildZExtOrBitCast(builder, Val, DestTy, Name))

sextorbitcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildSExtOrBitCast(builder, Val, DestTy, Name))

truncorbitcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildTruncOrBitCast(builder, Val, DestTy, Name))

cast!(builder::Builder, Op::API.LLVMOpcode, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildCast(builder, Op, Val, DestTy, Name))

pointercast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildPointerCast(builder, Val, DestTy, Name))

intcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildIntCast(builder, Val, DestTy, Name))

fpcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPCast(builder, Val, DestTy, Name))

# other operations

icmp!(builder::Builder, Op::API.LLVMIntPredicate, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildICmp(builder, Op, LHS, RHS, Name))

fcmp!(builder::Builder, Op::API.LLVMRealPredicate, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFCmp(builder, Op, LHS, RHS, Name))

phi!(builder::Builder, Ty::LLVMType, Name::String="") =
    Instruction(API.LLVMBuildPhi(builder, Ty, Name))

select!(builder::Builder, If::Value, Then::Value, Else::Value, Name::String="") =
    Value(API.LLVMBuildSelect(builder, If, Then, Else, Name))

call!(builder::Builder, Fn::Value, Args::Vector{<:Value}=Value[], Name::String="") =
    Instruction(API.LLVMBuildCall(builder, Fn, Args, length(Args), Name))

va_arg!(builder::Builder, List::Value, Ty::LLVMType, Name::String="") =
    Instruction(API.LLVMBuildVAArg(builder, List, Ty, Name))

landingpad!(builder::Builder, Ty::LLVMType, PersFn::Value, NumClauses::Integer, Name::String="") =
    Instruction(API.LLVMBuildLandingPad(builder, Ty, PersFn, NumClauses, Name))

neg!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildNeg(builder, V, Name))

nswneg!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildNSWNeg(builder, V, Name))

nuwneg!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildNUWNeg(builder, V, Name))

fneg!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildFNeg(builder, V, Name))

not!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildNot(builder, V, Name))

# other build methods

globalstring!(builder::Builder, Str::String, Name::String="") =
    Value(API.LLVMBuildGlobalString(builder, Str, Name))

globalstring_ptr!(builder::Builder, Str::String, Name::String="") =
    Value(API.LLVMBuildGlobalStringPtr(builder, Str, Name))

isnull!(builder::Builder, Val::Value, Name::String="") =
    Value(API.LLVMBuildIsNull(builder, Val, Name))

isnotnull!(builder::Builder, Val::Value, Name::String="") =
    Value(API.LLVMBuildIsNotNull(builder, Val, Name))

ptrdiff!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildPtrDiff(builder, LHS, RHS, Name))


text:
# add some basic wrappers for Lapack, because Julia Base was too lazy to
# expose an interface that allows reusing all the arrays
import Base.LinAlg.LAPACK.liblapack, Base.LinAlg.LAPACK.getrf!, Base.LinAlg.LAPACK.getrs!, Base.LinAlg.BlasInt
import Base.LinAlg.BLAS.@blasfunc

@eval begin
"""
  A better wrapper for Lapack DGETRF.  The Base function getrs! is adequate
  (although it does error checking every time) for solving the system using
  this factorization.  This function accepts Contiguous ArrayViews

  **Inputs/Outputs**

   * A: matrix to factorize (inplace)
   * ipiv: array of pivot indices

  **Outputs**

   * info: the Lapack info field
"""
function getrf!(A::ContiguousArrays{Float64}, ipiv::AbstractArray{BlasInt})

  # TODO: put this in a debug block
  @ifsafeview begin
    @assert stride(A, 1) == 1
    @assert stride(A, 2) == size(A, 1)
    @assert stride(ipiv, 1) == 1
    @assert length(ipiv) == min( size(A, 1), size(A, 2) )
  end

  info = Ref{BlasInt}()
  m, n = size(A)
  lda  = size(A, 1)
  ccall( (@blasfunc($(Symbol("dgetrf_"))), liblapack), Void,
#  ccall(( :dgetrf_64_, liblapack), Void,
        (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{Float64},
         Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
        &m, &n, A, &lda, ipiv, info)
  return info[]
end

"""
  A wrapper for DGETRS that does less error checking.  This function accepts
  Contiguous ArrayViews

  **Inputs**

   * trans: char
   * A: the matrix, previously factored by DGETRF
   * ipiv: ipiv from DGETRF

"""
function getrs2!(trans::Char, A::ContiguousArrays{Float64}, ipiv::AbstractVector{BlasInt}, B::ContiguousArrays{Float64})

  @ifsafeview begin
    @assert stride(B, 1) == 1
    @assert stride(B, 2) == size(B, 1)
    @assert size(A, 1) == size(A, 2)
    if size(A, 1) != size(B, 1)
        throw(DimensionMismatch("B has leading dimension $(size(B,1)), but needs $n"))
    end
  end

  m, n = size(A)
  nrhs = size(B, 2)
  info = Ref{BlasInt}()
  ccall((@blasfunc($(Symbol("dgetrs_"))), liblapack), Void,
        (Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{Float64}, Ptr{BlasInt},
         Ptr{BlasInt}, Ptr{Float64}, Ptr{BlasInt}, Ptr{BlasInt}),
        &trans, &n, &size(B,2), A, &max(1,stride(A,2)), ipiv, B, &max(1,stride(B,2)), info)
  return info[]
end

"""
  Wrapper for DLASWP
"""
function laswp!(A::ContiguousArrays{Float64, 1}, k1::BlasInt, k2::BlasInt,
                ipiv::AbstractVector{BlasInt}, incx::Integer=1)

  @assert length(ipiv) == k2*abs(incx)

  n = 1
  lda = 1  # length(A) ?
#  lda = max(1, stride(A, 2))
  ccall(( @blasfunc($(Symbol("dlaswp_"))), liblapack), Void, (Ptr{BlasInt}, Ptr{Float64},
            Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
            Ptr{BlasInt}), &n, A, &lda, &k1, &k2, ipiv, &incx)

  return nothing

end

end  # end @eval

"""
  This function applies the permutation matrix specified by `ipiv` from
  `getrf`.  Use `laswp!` for the inverse permutation matrix.
"""
function applyIpiv!(ipiv::AbstractVector{I}, x::AbstractVector) where {I <: Integer}

  @ifsafeview @assert length(x) == length(ipiv)

  # ipiv specifies the swaps to perform for the *inverse* permutation
  # matrix.  Loop backwards to get piv, the non-inverse permutations.
  for i=length(x):-1:1
    destval = x[ipiv[i]]
    x[ipiv[i]] = x[i]
    x[i] = destval
  end

  return nothing
end

import Base.SparseArrays.UMFPACK: UmfpackLU, umf_ctrl, umf_info,
                                  UmfpackIndexTypes, umf_nm, UMFPACK_A,
                                  UMFPACK_At, UMFPACK_Aat, umferror

for itype in UmfpackIndexTypes
  sol_r = umf_nm("solve", :Float64, itype)

  @eval begin
    function solve_suitesparse(U::UmfpackLU{Float64, $itype}, b::AbstractVector{Float64},
                              typ::Integer, x::AbstractVector{Float64})

      if length(b) != U.m
        throw(DimensionMismatch())
      end

      ret = ccall( ($sol_r, :libumfpack), $itype, ($itype, Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}), typ, U.colptr, U.rowval, U.nzval, x, b, U.numeric, umf_ctrl, umf_info)

      umferror(ret)

      return nothing
    end

  end  # end eval begin
end # loop over itype


text:
using Libtask
using Test

@testset "TArray Copying Test" begin
    DATA = Dict{Task, Array}()

    function f_ct()
        ta = tzeros(UInt64, 4);
        for i in 1:4
            ta[i] = hash(Libtask._current_task())
            DATA[Libtask._current_task()] = convert(Array, ta)
            produce(ta[i])
        end
    end

    t = CTask(f_ct)

    @test consume(t) == hash(t) # index = 1
    @test consume(t) == hash(t) # index = 2

    a = copy(t);

    @test consume(a) == hash(a) # index = 3
    @test consume(a) == hash(a) # index = 4

    @test consume(t) == hash(t) # index = 3

    @test DATA[t] == [hash(t), hash(t), hash(t), 0]
    @test DATA[a] == [hash(t), hash(t), hash(a), hash(a)]

end


text:
import LinearAlgebra, Distributions, Random, Statistics, DataFrames

"""
    simulate_coefs_correlation(coefs_mean::Number=0.1; coefs_sd::Number=0.1, n::Int=10)

Generate a vector of random correlation coefficients from a normal distribution.

# Arguments
- `coefs_mean::Number`: Mean of the normal distribution from which to get the coefs.
- `coefs_sd::Number`: SD of the normal distribution.
- `n::Int`: Number of coefficients.

# Examples
```julia
simulate_coefs_correlation(0.5)
```
"""
function simulate_coefs_correlation(coefs_mean::Number=0.1; coefs_sd::Number=0.1, n::Int=10)
  # Generate outcome
  coefs = Random.rand(Distributions.Normal(coefs_mean, coefs_sd), n)
  coefs = [x = [x] for x in coefs]
end

"""
    simulate_data_correlation(coefs; n::Int=100, noise::Number=0.0, groupnames=:random)

Generate a DataFrame of correlated variables.

# Multiple Variables / Groups
- If `coefs` is a vector (*e.g., `[0.1, 0.2]`), the DataFrame will contain `length(coefs)` variables (`Var1, Var2, ...`). Altough uncorrelated between them, they are  correlated to the outcome (`y`) by the specified coefs.
- If `coefs` is a vector of vectors (*e.g., `[[0.1], [0.2]]`), it will create `length(coefs)` groups, *i.e.*, stacked DataFrames with a correlation between the variables and the outcome varying between groups. It is possible to specify the `groupnames`.

# Arguments
- `coefs`: Correlation coefficients. Can be a number, a vector of numbers or a vector of vectors.
- `n::Int`: Number of observations.
- `noise::Number`: The SD of the random gaussian noise.
- `groupnames::Vector`: Vector of group names (default to `:random`).
- `kwargs...`: Arguments to pass to other functions.

!!! note

    **Ideas / help required:**
    - Different group sizes (See [#9](https://github.com/neuropsychology/Psycho.jl/issues/9))
    - Bug in some cases (*e.g.*, `simulate_data_correlation([0.2, 0.9, 0.5])`) related to failure in Cholesky factorization (See [#11](https://github.com/neuropsychology/Psycho.jl/issues/11))

# Examples
```julia
simulate_data_correlation(0.2)
```
"""
function simulate_data_correlation end

function simulate_data_correlation(coefs::Vector{<:Number}; n::Int=100, noise::Number=0.0)

  # Generate outcome
  y = standardize(Random.rand(Distributions.Normal(0, 1), n))

  n_var = length(coefs)
  X = standardize(Random.rand(Distributions.Normal(0, 1), n, n_var))
  X = hcat(y, X)

  # find the current correlation matrix
  cor_structure = Statistics.cov(X)
  chol = LinearAlgebra.cholesky(cor_structure).U

  # cholesky decomposition to get independence
  X = X / chol

  # create new correlation structure (zeros can be replaced with other r vals)
  coefs_structure = Matrix{Float64}(LinearAlgebra.I, n_var+1, n_var+1)
  coefs_structure[1, 2:end] = coefs
  coefs_structure[2:end, 1] = coefs

  X = X * LinearAlgebra.cholesky(coefs_structure).U
  X = X * Statistics.std(y) .+ Statistics.mean(y)

  X = X[:, 2:end]

  # Add noise
  if noise != 0
    y = y .+ Random.rand(Distributions.Normal(0, noise), n)
  end

  data = DataFrames.DataFrame(hcat(y, X),
    Symbol.(vcat(["y"], "Var" .* string.(1:n_var))))
  return data
end

function simulate_data_correlation(coefs::Number; n::Int=100, noise::Number=0.0)
  simulate_data_correlation([coefs], n=n, noise=noise)
end

function simulate_data_correlation(coefs::Vector{<:Vector}; n::Int=100, noise::Number=0.0, groupnames=:random, kwargs...)
  data = simulate_data_groups(simulate_data_correlation, coefs=coefs, n=n, noise=noise, groupnames=groupnames, kwargs...)
  return data
end


text:
using Documenter, AMDGPU

function main()
    ci = get(ENV, "CI", "") == "true"

    @info "Generating Documenter site"
    makedocs(
        sitename="AMDGPU.jl",
        pages = [
            "Home" => "index.md",
            "Quick Start" => "quickstart.md",
            "Devices/Agents" => "devices.md",
            "Queues and Signals" => "queues_signals.md",
            "Kernel Dependencies" => "kernel_deps.md",
            "Global Variables" => "globals.md",
            "Exceptions" => "exceptions.md",
            "Printing" => "printing.md",
            "Memory" => "memory.md",
            "Intrinsics" => [
                "Execution Control" => "execution_control.md",
                "Wavefront Operations" => "wavefront_ops.md",
            ],
            "API Reference" => "api.md"
        ]
    )
    if ci
        @info "Deploying to GitHub"
        deploydocs(
            repo = "github.com/JuliaGPU/AMDGPU.jl.git",
            push_preview = true
        )
    end
end

isinteractive() || main()


text:
##### weighted least squares
 struct SemWLS{Vt <: Union{AbstractArray, UniformScaling{Bool}},
        St <: AbstractArray} <: LossFunction
    V::Vt
    s::St
end

### Constructor
function SemWLS(observed::T, objective, grad; V = nothing) where {T <: SemObs}
    ind = CartesianIndices(observed.obs_cov)
    ind = filter(x -> (x[1] >= x[2]), ind)
    s = observed.obs_cov[ind]
    # compute V here
    if isnothing(V)
        D = duplication_matrix(observed.n_man)
        S = inv(observed.obs_cov)
        S = kron(S,S)
        V = 0.5*(D'*S*D)
    end
    return SemWLS(V, s)
end

### Loss
function (semwls::SemWLS)(par, model)
    diff = semwls.s - model.imply.imp_cov
    F = dot(diff, semwls.V, diff)
    return F
end

############## SWLS
struct SemSWLS{Vt <: Union{AbstractArray, UniformScaling{Bool}},
        St <: AbstractArray} <: LossFunction
    V::Vt
    sᵀV::St
    #s::St
end

### Constructor
function SemSWLS(observed::T, objective, grad; V = nothing) where {T <: SemObs}
    ind = CartesianIndices(observed.obs_cov)
    ind = filter(x -> (x[1] >= x[2]), ind)
    s = observed.obs_cov[ind]

    if isnothing(V)
        D = duplication_matrix(observed.n_man)
        S = inv(observed.obs_cov)
        S = kron(S,S)
        V = 0.5*(D'*S*D)
    end

    sᵀV = transpose(s)*V
    return SemSWLS(V, sᵀV)
end

### Loss
function (semswls::SemSWLS)(par, model)
    # s = semswls.s
    # V = semswls.V
    # G = model.imply.G
    outer = semswls.sᵀV*model.imply.G
    b = cholesky(Symmetric(model.imply.G'*semswls.V*model.imply.G))
    inner = b\outer'
    F = -outer*inner
    # F = transpose(s)*V*s - transpose(s)*V*G*inv(transpose(G)*V*G)*(transpose(G)*V*s)
#=     θ = inv(transpose(G)*V*G)*(transpose(G)*V*s)
    σ = G*θ
    diff = s - σ
    F = dot(diff, V, diff) =#
    return F
end


text:
export SchedulerHaltedException
export sch_handle, halt!, exec!, get_dag_ids, add_thunk!

"Identifies a thunk by its ID."
struct ThunkID
    id::Int
end

"A handle to the scheduler, used by dynamic thunks."
struct SchedulerHandle
    thunk_id::ThunkID
    out_chan::RemoteChannel
    inp_chan::RemoteChannel
end

"Gets the scheduler handle for the currently-executing thunk."
sch_handle() = task_local_storage(:sch_handle)::SchedulerHandle

"Thrown when the scheduler halts before finishing processing the DAG."
struct SchedulerHaltedException <: Exception end

function safepoint(state)
    if state.halt[]
        # Force dynamic thunks and listeners to terminate
        for (inp_chan,out_chan) in values(state.worker_chans)
            close(inp_chan)
            close(out_chan)
        end
        # Throw out of scheduler
        throw(SchedulerHaltedException())
    end
end

"Processes dynamic messages from worker-executing thunks."
function dynamic_listener!(ctx, state)
    task = current_task()
    for tid in keys(state.worker_chans)
        inp_chan, out_chan = state.worker_chans[tid]
        @async begin
            while isopen(inp_chan)
                tid, f, data = try
                    take!(inp_chan)
                catch err
                    if !(unwrap_nested_exception(err) isa Union{SchedulerHaltedException,
                                                                ProcessExitedException,
                                                                InvalidStateException})
                        @error exception=(err,catch_backtrace())
                    end
                    break
                end
                res = try
                    (false, lock(state.lock) do
                        Base.invokelatest(f, ctx, state, task, tid, data)
                    end)
                catch err
                    (true, RemoteException(CapturedException(err,catch_backtrace())))
                end
                try
                    put!(out_chan, res)
                catch err
                    if !(unwrap_nested_exception(err) isa Union{SchedulerHaltedException,
                                                                ProcessExitedException,
                                                                InvalidStateException})
                        @error exception=(err,catch_backtrace())
                    end
                end
            end
        end
    end
end

## Worker-side methods for dynamic communication

"Executes an arbitrary function within the scheduler, returning the result."
function exec!(f, h::SchedulerHandle, args...)
    put!(h.out_chan, (h.thunk_id.id, f, args))
    failed, res = take!(h.inp_chan)
    failed && throw(res)
    res
end

"Commands the scheduler to halt execution immediately."
halt!(h::SchedulerHandle) = exec!(_halt, h, nothing)
function _halt(ctx, state, task, tid, _)
    state.halt[] = true
    put!(state.chan, (1, nothing, SchedulerHaltedException(), nothing))
    Base.throwto(task, SchedulerHaltedException())
end

"Waits on a thunk to complete, and fetches its result."
function Base.fetch(h::SchedulerHandle, id::ThunkID)
    future = ThunkFuture(Future(1))
    exec!(_register_future, h, future, id.id)
    fetch(future; proc=thunk_processor())
end
"Waits on a thunk to complete, and fetches its result."
function register_future!(h::SchedulerHandle, id::ThunkID, future::ThunkFuture)
    exec!(_register_future, h, future, id.id)
end
function _register_future(ctx, state, task, tid, (future, id))
    @assert tid != id "Cannot fetch own result"
    thunk = state.thunk_dict[id]
    ownthunk = state.thunk_dict[tid]
    dominates(target, t) = (t == target) || any(_t->dominates(target, _t), filter(istask, t.inputs))
    @assert !dominates(ownthunk, thunk) "Cannot fetch result of dominated thunk"
    if thunk in state.finished || thunk in state.errored
        error = thunk in state.errored
        if haskey(state.cache, thunk)
            put!(future, state.cache[thunk]; error=error)
        else
            put!(future, nothing)
        end
    else
        futures = get!(()->ThunkFuture[], state.futures, thunk)
        push!(futures, future)
        schedule!(ctx, state)
    end
    nothing
end

# TODO: Optimize wait() to not serialize a Chunk
"Waits on a thunk to complete."
function Base.wait(h::SchedulerHandle, id::ThunkID; future=ThunkFuture(1))
    register_future!(h, id, future)
    wait(future)
end

"Returns all Thunks IDs as a Dict, mapping a Thunk to its downstream dependents."
get_dag_ids(h::SchedulerHandle) =
    exec!(_get_dag_ids, h, nothing)::Dict{ThunkID,Set{ThunkID}}
function _get_dag_ids(ctx, state, task, tid, _)
    deps = Dict{ThunkID,Set{ThunkID}}()
    for (key,val) in state.dependents
        deps[ThunkID(key.id)] = Set(map(t->ThunkID(t.id), collect(val)))
    end
    deps
end

"Adds a new Thunk to the DAG."
add_thunk!(f, h::SchedulerHandle, args...; kwargs...) =
    ThunkID(exec!(_add_thunk!, h, f, args, kwargs)::Int)
function _add_thunk!(ctx, state, task, tid, (f, args, kwargs))
    _args = map(arg->arg isa ThunkID ? state.thunk_dict[arg.id] : arg, args)
    thunk = Thunk(f, _args...; kwargs...)
    state.thunk_dict[thunk.id] = thunk
    state.dependents[thunk] = Set{Thunk}()
    @assert reschedule_inputs!(state, thunk) || (thunk in state.errored)
    schedule!(ctx, state)
    return thunk.id
end


text:
const ldastanmodel = "
data {
  int<lower=2> K;               // num topics
  int<lower=2> V;               // num words
  int<lower=1> M;               // num docs
  int<lower=1> N;               // total word instances
  int<lower=1,upper=V> w[N];    // word n
  int<lower=1,upper=M> doc[N];  // doc ID for word n
  vector<lower=0>[K] alpha;     // topic prior
  vector<lower=0>[V] beta;      // word prior
}
parameters {
  simplex[K] theta[M];   // topic dist for doc m
  simplex[V] phi[K];     // word dist for topic k
}
model {
  for (m in 1:M)
    theta[m] ~ dirichlet(alpha);  // prior
  for (k in 1:K)
    phi[k] ~ dirichlet(beta);     // prior
  for (n in 1:N) {
    real gamma[K];
    for (k in 1:K)
      gamma[k] <- log(theta[doc[n],k]) + log(phi[k,w[n]]);
    increment_log_prob(log_sum_exp(gamma));  // likelihood
  }
}
"


text:
export Grid, δx, δy, δz, δ, x, y, z, x⃗, xc, yc, zc, x⃗c, N, g⃗, _fftaxes, Nranges, ∫

struct Grid{ND,T}
	Δx::T
    Δy::T
    Δz::T
	Nx::Int
    Ny::Int
    Nz::Int
end

Grid(Δx::T, Δy::T, Δz::T, Nx::Int, Ny::Int, Nz::Int) where {T<:Real} = Grid{3,T}(
	Δx,
    Δy,
    Δz,
	Nx,
    Ny,
    Nz,
)

Grid(Δx::T, Δy::T, Nx::Int, Ny::Int) where {T<:Real} = Grid{2,T}(
	Δx,
    Δy,
    1.,
	Nx,
    Ny,
    1,
)

function Base.size(gr::Grid{3})::NTuple{3,Int}
	(gr.Nx, gr.Ny, gr.Nz)
end
function Base.size(gr::Grid{2})::NTuple{2,Int}
	(gr.Nx, gr.Ny)
end
Base.size(gr::Grid,d::Int) = Base.size(gr)[d]

δx(g::Grid) = g.Δx / g.Nx
δy(g::Grid) = g.Δy / g.Ny
δz(g::Grid) = g.Δz / g.Nz

δ(g::Grid{2}) = g.Δx * g.Δy / ( g.Nx * g.Ny )
δ(g::Grid{3}) = g.Δx * g.Δy * g.Δz / ( g.Nx * g.Ny * g.Nz )

∫(intgnd;g::Grid) = sum(intgnd)*δ(g)	# integrate a scalar field over the grid

# voxel/pixel center positions
function x(g::Grid{ND,T})::Vector{T}  where {ND,T<:Real}
 	# ( ( g.Δx / g.Nx ) .* (0:(g.Nx-1))) .- g.Δx/2.
	LinRange(-g.Δx/2, g.Δx/2 - g.Δx/g.Nx, g.Nx)
end
function y(g::Grid{ND,T})::Vector{T}  where {ND,T<:Real}
 	# ( ( g.Δy / g.Ny ) .* (0:(g.Ny-1))) .- g.Δy/2.
	LinRange(-g.Δy/2, g.Δy/2 - g.Δy/g.Ny, g.Ny)
end
function z(g::Grid{ND,T})::Vector{T}  where {ND,T<:Real}
 	# ( ( g.Δz / g.Nz ) .* (0:(g.Nz-1))) .- g.Δz/2.
	LinRange(-g.Δz/2, g.Δz/2 - g.Δz/g.Nz, g.Nz)
end
function x⃗(g::Grid{2,T})::Array{SVector{3,T},2} where T<:Real
	( (xx,yy) = (x(g),y(g)); [SVector{3,T}(xx[ix],yy[iy],0.) for ix=1:g.Nx,iy=1:g.Ny] ) # (Nx × Ny ) 2D-Array of (x,y,z) vectors at pixel/voxel centers
end
function x⃗(g::Grid{3,T})::Array{SVector{3,T},3} where T<:Real
	( (xx,yy,zz) = (x(g),y(g),z(g)); [SVector{3,T}(xx[ix],yy[iy],zz[iz]) for ix=1:g.Nx,iy=1:g.Ny,iz=1:g.Nz] ) # (Nx × Ny × Nz) 3D-Array of (x,y,z) vectors at pixel/voxel centers
end

# voxel/pixel corner positions
function xc(g::Grid{ND,T})::Vector{T} where {ND,T<:Real}
	collect( ( ( g.Δx / g.Nx ) .* (0:g.Nx) ) .- ( g.Δx/2. * ( 1 + 1. / g.Nx ) ) )
	# collect(range(-g.Δx/2.0, g.Δx/2.0, length=g.Nx+1))
end
function yc(g::Grid{ND,T})::Vector{T} where {ND,T<:Real}
	collect( ( ( g.Δy / g.Ny ) .* (0:g.Ny) ) .- ( g.Δy/2. * ( 1 + 1. / g.Ny ) ) )
	# collect(range(-g.Δy/2.0, g.Δy/2.0, length=g.Ny+1))
end
function zc(g::Grid{3,T})::Vector{T} where T<:Real
	collect( ( ( g.Δz / g.Nz ) .* (0:g.Nz) ) .- ( g.Δz/2. * ( 1 + 1. / g.Nz ) ) )
	# collect(range(-g.Δz/2.0, g.Δz/2.0, length=g.Nz+1))
end
function x⃗c(g::Grid{2,T})::Array{SVector{3,T},2}  where T<:Real
	( (xx,yy) = (xc(g),yc(g)); [SVector{3}(xx[ix],yy[iy],0.) for ix=1:(g.Nx+1),iy=1:(g.Ny+1)] )
end
function x⃗c(g::Grid{3,T})::Array{SVector{3,T},3}  where T<:Real
	( (xx,yy,zz) = (xc(g),yc(g),zc(g)); [SVector{3}(xx[ix],yy[iy],zz[iz]) for ix=1:(g.Nx+1),iy=1:(g.Ny+1),iz=1:(g.Nz+1)] )
end

# grid size
@inline N(g::Grid)::Int = *(size(g)...)

import Base: eachindex
@inline eachindex(g::Grid) = CartesianIndices(size(g)) #(1:NN for NN in size(g))

# reciprocal lattice vectors (from fftfreqs)
function g⃗(gr::Grid{3,T})::Array{SVector{3, T}, 3} where T<:Real
	[	SVector(gx,gy,gz) for   gx in fftfreq(gr.Nx,gr.Nx/gr.Δx),
								gy in fftfreq(gr.Ny,gr.Ny/gr.Δy),
							    gz in fftfreq(gr.Nz,gr.Nz/gr.Δz)		]
end
_fftaxes(gr::Grid{3}) = (2:4)

function g⃗(gr::Grid{2,T})::Array{SVector{3, T}, 2} where T
	 [	SVector(gx,gy,0.) for   gx in fftfreq(gr.Nx,gr.Nx/gr.Δx),
								gy in fftfreq(gr.Ny,gr.Ny/gr.Δy)		]
end
_fftaxes(gr::Grid{2}) = (2:3)

"""
################################################################################
#																			   #
#							   Plotting methods					   			   #
#																			   #
################################################################################
"""


text:
@testset "Basic thermo" begin
    C₁C₅ = CubicEoS.load(BrusilovskyEoSMixture; names=("methane", "n-pentane"))
    RT = 300 * CubicEoS.GAS_CONSTANT_SI
    @testset "Calls w/wo buffers" begin
        nc = ncomponents(C₁C₅)
        buffer_st = thermo_buffer(C₁C₅)
        buffer_nt = (ai=zeros(nc), aij=zeros((nc, nc)))
        buffer_dict = Dict(:ai => zeros(nc), :aij => zeros((nc, nc)))
        # the following tests are grouped by
        # 1. no buffer vs buffer struct
        # 2. no buffer vs buffer namedtuple
        # 3. no buffer vs buffer dict
        N = 2 * [0.8, 1-0.8]
        V = 2.5
        for buffer in (buffer_st, buffer_nt, buffer_dict)
            @test CubicEoS.eos_parameters(C₁C₅, N, RT) == CubicEoS.eos_parameters(C₁C₅, N, RT; buf=buffer)
            @test pressure(C₁C₅, N, V, RT) == pressure(C₁C₅, N, V, RT; buf=buffer)
        end
    end
    @testset "Pressure" begin
        @testset "Homogeneity" begin
            V = 1
            χ = [0.8, 1-0.8]
            for ∑N in 0.1:0.2:1.9
                N = χ * ∑N
                υ = V / ∑N
                @test pressure(C₁C₅, N, V, RT) ≈ pressure(C₁C₅, χ, υ, RT)
            end
        end
    end
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
using Test
using Revise
using Pkg; Pkg.activate(".")

@testset "loadtest" begin
    @test (using dlgo; true) == true
end


text:
#=
BasicPrimTest:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/BasicPrim.jl")
matrix =  GBMatrix([[0.0, 7.0, 5.0, 5.0] [7.0, 0.0, 0.0, 0.0] [5.0, 0.0, 0.0, 0.0] [5.0, 0.0, 0.0, 0.0]])
m = GBVector([Inf, 0.0, 0.0, 0.0]) #source matrix
@test b_prim(matrix, 4, m) == 17.0


text:
function get_decomp_filename(m_new::M, m_old::M, input_type::Symbol,
                             cond_new::Symbol, cond_old::Symbol,
                             product::Symbol, class::Symbol;
                             pathfcn::Function = rawpath,
                             fileformat::Symbol = :jld2,
                             forecast_string_new = "", forecast_string_old = "") where M<:AbstractDSGEModel
    output_var = Symbol(product, class)

    fn_new = get_forecast_filename(m_new, input_type, cond_new, output_var,
                                   pathfcn = pathfcn, fileformat = fileformat, forecast_string = forecast_string_new)
    fn_old = get_forecast_filename(m_old, input_type, cond_old, output_var, forecast_string = forecast_string_old)

    dir = dirname(fn_new)
    base_new, ext = splitext(basename(fn_new))
    base_old, _   = splitext(basename(fn_old))
    base_old = replace(base_old, string(output_var) => "")
    base_old = spec(m_old) * "_" * subspec(m_old) * base_old

    return joinpath(dir, base_new * "__" * base_old * ext)
end

function get_decomp_output_files(m_new::M, m_old::M, input_type::Symbol,
                                 cond_new::Symbol, cond_old::Symbol,
                                 classes::Vector{Symbol}; forecast_string_new = "", forecast_string_old = "") where M<:AbstractDSGEModel
    output_files = Dict{Symbol, String}()
    for comp in [:data, :news, :shockdec, :dettrend, :para, :total]
        for class in classes
            product = Symbol(:decomp, comp)
            output_var = Symbol(product, class)
            output_files[output_var] =
                get_decomp_filename(m_new, m_old, input_type, cond_new, cond_old,
                                    product, class, pathfcn = rawpath, fileformat = :jld2,
                                    forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old)
        end
    end
    return output_files
end

function get_decomp_mean_file(m_new::M, m_old::M, input_type::Symbol,
                              cond_new::Symbol, cond_old::Symbol,
                              class::Symbol; forecast_string_new = "", forecast_string_old = "") where M<:AbstractDSGEModel
    get_decomp_filename(m_new, m_old, input_type, cond_new, cond_old, :decomp, class,
                        pathfcn = workpath, fileformat = :jld2,
                        forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old)
end

function write_forecast_decomposition(m_new::M, m_old::M, input_type::Symbol,
                                      classes::Vector{Symbol},
                                      decomp_output_files::Dict{Symbol, String},
                                      decomps::Dict{Symbol, Array{Float64}};
                                      block_number::Nullable{Int} = Nullable{Int}(),
                                      block_inds::AbstractRange{Int} = 1:0,
                                      verbose::Symbol = :low,
                                      forecast_string_new = "", forecast_string_old = "") where M<:AbstractDSGEModel
    for comp in [:data, :news, :shockdec, :dettrend, :para, :total]
        for class in classes
            prod = Symbol(:decomp, comp)
            var = Symbol(prod, class)
            filepath = decomp_output_files[var]

            if isnull(block_number) || get(block_number) == 1
                # Write forecast metadata to a jld2 and the raw forecast output
                # to an h5. The data in the HDF5 will be transferred to the jld2
                # and the h5 file will be deleted when combine_raw_forecast_output_and_metadata
                # is executed.
                JLD2.jldopen(filepath, "w") do file
                    # Write metadata
                    # Pass in m_old because its historical and forecast dates are used
                    write_forecast_metadata(m_old, file, var)
                end

                h5open(replace(filepath, "jld2" => "h5"), "w") do file
                    # Pre-allocate HDF5 dataset which will contain all draws
                    if !isnull(block_number) && get(block_number) == 1
                        # Determine forecast output size
                        jstep = get_jstep(m_new, n_forecast_draws(m_new, input_type))
                        dims = collect(size(decomps[var]))
                        dims[1] = convert(Int, n_forecast_draws(m_new, input_type) / jstep)

                        # Determine chunk (block) size
                        block_size = convert(Int, forecast_block_size(m_new) / jstep)
                        chunk_dims = copy(dims)
                        chunk_dims[1] = block_size

                        # Initialize dataset
                        # pfile = file.plain
                        HDF5.d_create(file, "arr", datatype(Float64), dataspace(dims...), "chunk", chunk_dims)
                    end
                end
            end

            h5open(replace(filepath, "jld2" => "h5"), "r+") do file
                if isnull(block_number)
                    write(file, "arr", decomps[var])
                else
                    write_forecast_block(file, decomps[var], block_inds)
                end
            end

            println(verbose, :high, " * Wrote $(basename(filepath))")
        end # of loop over classes
    end # of loop over comps
end


text:
using CSV, DelimitedFiles, DataFrames

"""
    flat_square_2d_grid(n, a, b)

constructs a regular grid in interval [a,b]^2

# Arguments
- `n:Integer`: total number of points sqrtn^2 with sqrtn = floor(sqrt(n))^2)
- `a:Float`: start point of interval [a,b].
- `b:Float`: end point of interval [a,b].

# Output
- `X::Array{Float64,2}`:  nx2 array with coordinates of n grid nodes with [a,b]
position (i,j) -> row = j + sqrtn (i-1) for i,j = 1, ..., sqrtn

"""
function flat_square_2d_grid(n, a, b)

    sqrtn = Int64(floor(sqrt(n)));
    X = zeros(sqrtn*sqrtn, 2);
    counter = 0;
    for i in 1:sqrtn
        for j in 1:sqrtn
            counter += 1;
            X[counter,1] =  a + (b-a)*(i-1)/(sqrtn-1);
            X[counter,2] =  a + (b-a)*(j-1)/(sqrtn-1);
        end
    end

    return X;
end

"""
    add_DPP_samples_to_Fredholm_samples(s,Fredholm_sample,intensity)

Loads DPP samples and concatenate them with Fredholm samples

# Arguments
- `s:Integer`: number of DPP samples.
- `Fredholm_sample:Array`: array containing samples to approximate Fredholm determinant.
- `intensity:Integer`: intensity of the DPP producing the samples (50 or 100).

# Outputs

- `all_sples`: array containing all the concatenated samples.
- `indices_Fredholm_sple`: array containing the indices of the Fredholm samples.
- `indices_DPP_sples`: array of arrays. The l-th entry of the array is an array containing the indices of the l-th DPP sample.

"""
function add_DPP_samples_to_Fredholm_samples(s,Fredholm_sple,intensity)

    # intensity of generating process
    strg_intensity = string(intensity)

    # create an array containing s arrays of DPP samples
    indices_DPP_sples = Array{Int64,1}[];
    n = size(Fredholm_sple,1)
    indices_Fredholm_sple = collect(1:n); # the set I
    all_sples = Fredholm_sple;

    print("Loading DPP samples from: \n")
    id_last_sple = n;
    for i = 0:(s-1)
        # read files and specify that they start from first row
        file_name = "data/statspats/samples/GaussDPPsample_alpha0_00p5_rho0_"*strg_intensity*"_nb_"*string(i+1)*".csv"
        print(file_name, " \n")
        sple = CSV.File(file_name; header=true) |> Tables.matrix
        sple = sple[:,2:3]; # remove first column
        id_sple = collect((id_last_sple+1):(id_last_sple + size(sple,1)));
        id_last_sple = id_last_sple + size(sple,1);
        # add i-th dpp sample indices
        push!(indices_DPP_sples,id_sple);
        # array with all samples
        all_sples = [all_sples; sple];
    end
    print("Done \n")

    return all_sples, indices_Fredholm_sple, indices_DPP_sples
end

"""
    plot_objectives(i_start,i_stop, obj)

Plots the objective functions (obj) from time i_start to time i_stop. No output.

# Arguments
- `i_start::Integer`: index of the starting time
- `i_stop::Integer`: index of the stopping time
- `obj:Array`: array containing the objective values
"""
function plot_objectives(i_start,i_stop, obj)

    range = i_start:i_stop;
    # plotting
    plt_objectives = plot(range,obj[range],title="objective values",xlabel="number of iterations",ylabel="objective value",legend=false);
    display(plt_objectives)

end

"""
    plot_insample_likelihood(B, R, indices_Fredholm_sple, total_sples)

Plots the in-sample likelihood.

# Arguments
- `B::Array{Float64,2}`: positive definite matrix B solution of the discrete optim problem.
- `R::Array{Float64,2}`: Cholesky factor such that K = R'R.
- `indices_Fredholm_sple:Array{Int64,1}`: indices of samples for approximating Fredholm determinant.
- `total_sples`: array containing all the samples used for representing the solution.

"""
function plot_insample_likelihood(B, R, indices_Fredholm_sple, total_sples)

    diagL = diag(R'*B*R);

    # identify uniform samples
    Fredholm_sples = total_sples[indices_Fredholm_sple,:];
    color_Fredholm_sples = diagL[indices_Fredholm_sple];

    # identify DPP samples
    indices_DPP_sples = setdiff(collect(1:size(total_sples,1)),indices_Fredholm_sple);
    DPP_sples = total_sples[indices_DPP_sples,:];
    color_DPP_sples = diagL[indices_DPP_sples];

    # plot in-sample likelihood
    plt_in_sple_likelihood = scatter(Fredholm_sples[:,1],Fredholm_sples[:,2],zcolor=color_Fredholm_sples,marker = :cross,markersize = 3, title="in-sample likelihood",label="unif",legend=true)
    scatter!(DPP_sples[:,1],DPP_sples[:,2],zcolor=color_DPP_sples,marker = :circle,markersize = 3,legend = false,colorbar = true,framestyle=:box,xtickfont = font(10),ytickfont = font(10),label="DPP")
    display(plt_in_sple_likelihood)
end

"""
    heatmap_intensity_correlation_kernel(GramK,n_side)

Plots a heatmap of the intensity (diagonal of Gram matrix) on a grid.

# Arguments
- `GramK::Array{Float64,2}`: Gram matrix of a kernel on a grid.
- `n_side::Integer`: number of points along each dim of the grid [0:(1/(n_side - 1)):1]^2.

"""
function heatmap_intensity_correlation_kernel(GramK,n_side)

    x_tics = 0:(1/(n_side - 1)):1;
    y_tics = x_tics;

    # plotting intensity
    IntensityGramK = reshape(diag(GramK),(n_side,n_side));
    plt_intensity_K = heatmap(x_tics,y_tics,IntensityGramK,colorbar = true,xtickfont = font(10),ytickfont = font(10),title= "Intensity of learned DPP ");
    display(plt_intensity_K)

end

"""
    heatmap_slice_Gram_matrix_at_grid_center(GramK,n_side)

Plots a heatmap of a slice of a Gram matrix on a grid.
The slice corresponds to a choice of point in the center of the grid (approximately).

# Arguments
- `GramK::Array{Float64,2}`: Gram matrix of a kernel on a grid.
- `n_side::Integer`: number of points along each dim of the grid [0:(1/(n_side - 1)):1]^2.

"""
function heatmap_slice_Gram_matrix_at_grid_center(GramK,n_side)

        x_tics = 0:(1/(n_side -1)):1;
        y_tics = x_tics;

        ## plotting slice of Gram matrix
        i = floor(n_side/2);
        j = floor(n_side/2);
        # center of the grid
        id_slice = Int64(j + n_side*(i-1));
        GramK_slice = GramK[:,id_slice];

        # reshaping
        GramK_slice_reshaped = reshape(GramK_slice,(n_side,n_side));
        plt = heatmap(x_tics,y_tics,GramK_slice_reshaped,xtickfont = font(10),ytickfont = font(10),title= "correlation kernel for one fixed argument");
        display(plt)
end

"""
    folder_name(s,n,p,sigma,lambda,tol)

outputs a string with a folder name identifying the parameters.

# Arguments
- `s::Integer`: number of DPP samples.
- `n::Integer`: number of Fredholm samples.
- `sigma::Float64': kernel bandwidth.
- `lambda::Float64': regularization parameter.
- `tol::Float64': stopping criterion for optimization algo.

# Outputs
- `foldername::String`: a string with the foldername.

"""
function folder_name(s,n,p,sigma,lambda,tol)

    foldername = "results/result_s="*string(Int64(s))*"_n="*string(Int64(n))*"_p="*string(Int64(p))*"_divideby1e3sigma="*string(Int64(1e3*sigma))*"_divideBy1e6lambda="*string(Int64(1e6*lambda))*"divideBy1e6_tol="*string(Int64(1e6*tol));

    return foldername
end

"""
    save_results_to_jld(B, R, n, total_sples, GramK, GramA, GramK_0, obj, i_stop, foldername)

Save variables into jld file. No output.
"""
function save_results_to_jld(B, R, n, total_sples, GramK, GramA, GramK_0, obj, i_stop, foldername)

    filename = foldername*"/outputs.jld";

    save(filename, "B",B, "R",R, "n",n, "total_samples",total_sples, "GramK",GramK, "GramA",GramA, "GramK0",GramK_0, "obj",obj, "i_stop",i_stop)

end

"""
    load_results_from_jld(foldername)

Loads variables from jld file.
"""
function load_results_from_jld(foldername)

    ## load results
    filename = foldername*"/outputs.jld";
    dict = load(filename);

    B =  dict["B"];
    R = dict["R"];
    n = dict["n"];
    total_sples = dict["total_samples"];
    GramK = ["GramK"];
    GramA = dict["GramA"];
    GramK_0 = dict["GramK0"];
    obj =  dict["obj"];
    i_stop = dict["i_stop"];

    return B, R, n, total_sples, GramK, GramA, GramK_0, obj, i_stop
end

"""
    create_repo(s,n,p,sigma,lambda,tol)

creates a repo for storing output data.

# Arguments
- `s::Integer`: number of DPP samples.
- `n::Integer`: number of Fredholm samples.
- `sigma::Float64': kernel bandwidth.
- `lambda::Float64': regularization parameter.
- `tol::Float64': stopping criterion for optimization algo.

# Outputs
- `foldername::String`: a string with the foldername.
"""
function create_repo(s,n,p,sigma,lambda,tol)
    foldername = folder_name(s,n,p,sigma,lambda,tol);
    if isdir(foldername)==false
        mkdir(foldername);
    end
    return foldername
end


text:
module SignalBase
using Unitful

include("api.jl")
include("units.jl")

end # module


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    Cartesian3D, Segment3D, Point3D,
    rotation

struct Cartesian3D{T<:AbstractMatrix{Float64}} <: Coordinates
    coords::T
end

struct Segment3D{T<:AbstractMatrix{Float64}}
    coords::T
end

struct Point3D{T<:AbstractVector{Float64}}
    pos::T
end

cartesian3d(data::AbstractMatrix{Float64}) = Cartesian3D(data)
cartesian3d(c3d::Cartesian3D, vertices::Vector{Int}
    ) = Cartesian3D(@view c3d.coords[vertices, :])

"""
    cartesian2d(coords::Cartesian3D) -> Cartesian2D

Embed `Cartesian3D` into `Cartesian2D` by just removing z-axis component.
"""
cartesian2d(c3d::Cartesian3D) = cartesian2d(@view c3d.coords[:, 1:2])

segment(c3d::Cartesian3D, u, v) = Segment3D(@view c3d.coords[[u, v], :])
segment(u::Point3D, v::Point3D) = Segment3D(transpose(hcat(u.pos, v.pos)))

point(c3d::Cartesian3D, i::Int) = Point3D(@view c3d.coords[i, :])

x(point::Point3D) = point.pos[1]
y(point::Point3D) = point.pos[2]
z(point::Point3D) = point.pos[3]

u(segment::Segment3D) = Point3D(@view segment.coords[1, :])
v(segment::Segment3D) = Point3D(@view segment.coords[2, :])
ux(segment::Segment3D) = segment.coords[1]
uy(segment::Segment3D) = segment.coords[3]
uz(segment::Segment3D) = segment.coords[5]
vx(segment::Segment3D) = segment.coords[2]
vy(segment::Segment3D) = segment.coords[4]
vz(segment::Segment3D) = segment.coords[6]

x_components(c3d::Cartesian3D) = c3d.coords[:, 1]
y_components(c3d::Cartesian3D) = c3d.coords[:, 2]
z_components(c3d::Cartesian3D) = c3d.coords[:, 3]

vector(segment::Segment3D) = v(segment) - u(segment)

function rotation(axis::Point3D, angle)
    (x1, y1, z1) = (x(axis), y(axis), z(axis))
    c = cos(angle)
    s = sin(angle)
    a12 = x1 * y1 * (1 - c)
    a13 = x1 * z1 * (1 - c)
    a23 = y1 * z1 * (1 - c)
    return Matrix{Float64}([
        (c + x1^2 * (1 - c)) (a12 - z1 * s) (a13 + y1 * s);
        (a12 + z1 * s) (c + y1^2 * (1 - c)) (a23 - x1 * s);
        (a13 - y1 * s) (a23 + x1 * s) (c + z1^2 * (1 - c))
    ])
end


text:
using Dates
using Printf

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("Proj4 not installed properly, run Pkg.build(\"Proj4\"), restart Julia and try again")
end
include(depsjl_path)
check_deps()

const FILENAMES = (epsg_path,
                   esri_path)

parse_projection(line::String) = match(r"<(\d+)>\s+(.*?)\s+<>", line)

println("# Contents of this file is generated on $(now()). Do not edit by hand!")
println("# \$ julia generate_projection_codes.jl > projection_codes.jl\n")

for filename in FILENAMES
    lines = open(filename) do proj_file
        readlines(proj_file)
    end

    println("$(basename(filename)) = Dict(")
    for i=1:length(lines)
        m = parse_projection(lines[i])
        if m === nothing
            continue
        else
            proj_code, proj_string = m.captures
            @printf("  %s => \"%s\", %s\n", proj_code, proj_string, lines[i-1])
        end
    end
    println(")\n")
end


text:
#= # for local testing, need this:
using Test
using PhyloNetworks
using PhyloPlots
=#

@testset "unconstrained NNI moves" begin

str_level1 = "(((8,9),(((((1,2,3),4),(5)#H1),(#H1,(6,7))))#H2),(#H2,10));"
net_level1 = readTopology(str_level1);
# same topology as: rootatnode!(net_level1, -3). edges 1:22
str_nontreechild = "((((Ag,E))#H3,(#H1:7.159::0.056,((M:0.0)#H2:::0.996,(Ak,(#H3:0.08,#H2:0.0::0.004):0.023):0.078):2.49):2.214):0.026,((Az:2.13,As:2.027):1.697)#H1:0.0::0.944,Ap);"
net_nontreechild = readTopology(str_nontreechild);
# problem: the plot has an extra vertical segment, for a clade that's not in the major tree
# --> fix that in PhyloPlots (fixit)
str_hybridladder = "(#H2:::0.2,((C,((B)#H1)#H2:::0.8),(#H1,(A1,A2))),O);"
net_hybridladder = readTopology(str_hybridladder);

#=
plot(net_level1, :R, showNodeNumber=true, showEdgeNumber=true)
plot(net_nontreechild, :R, showNodeNumber=true, showEdgeNumber=true)
plot(net_hybridladder, :R, showNodeNumber=true, showEdgeNumber=true)
=#

@test isnothing(PhyloNetworks.nni!(net_level1, net_level1.edge[1], 0x01)) # external edge

@testset "level1 edge 3: BB undirected move $move" for move in 0x01:0x08
    undoinfo = PhyloNetworks.nni!(net_level1, net_level1.edge[3], move);
    #location of v node (number -4)
    nodes = [n.number for n in net_level1.edge[20].node] #α's connections
    if move in [1, 2, 6, 8] #check that edge α connected to v
        @test -4 in nodes
    else
        @test !(-4 in nodes)
    end
    #location of u node (number -3)
    nodes = [n.number for n in net_level1.edge[1].node]
    if move in [2, 4, 5, 6] #check that edge γ connected to u
        @test -3 in nodes
    else
        @test !(-3 in nodes)
    end
    # check directionality
    if move in [3, 4, 5, 7]
        # keeping α, β or flipping uv keeps node -4 as child of edge 3
        @test PhyloNetworks.getChild(net_level1.edge[3]).number == -4
    else
        # switching α, β AND flipping uv or doing neither makes node -3 child of edge 3
        @test PhyloNetworks.getChild(net_level1.edge[3]).number == -3
    end
    # undo move
    PhyloNetworks.nni!(undoinfo...);
    # confirm we're back to original topology
    @test writeTopology(net_level1) == str_level1
end # of level1 edge 3: BB undirected

@testset "level1 test edge 13: BR directed move $move" for move in 0x01:0x03
    # e.hybrid and tree parent:  BR case, 3 moves because e cannot contain the root
    # 3cycle test: α connected to γ so move 1 will create 3 cycles.
    if move == 0x01
        @test isnothing(PhyloNetworks.nni!(net_level1, net_level1.edge[13], move)) # would create a 3cycle
    # move 1 would create a 3 cycle, but should work if we don't forbid 3cycles
        undoinfo = PhyloNetworks.nni!(net_level1, net_level1.edge[13], move, false); #no3cycle=false
        nodes = [n.number for n in net_level1.edge[17].node]
        @test 8 in nodes # check that edge α connected to v
        nodes = [n.number for n in net_level1.edge[11].node]
        @test !(-11 in nodes)
        @test PhyloNetworks.getChild(net_level1.edge[13]).number == -11
        PhyloNetworks.nni!(undoinfo...); # undo move
        @test writeTopology(net_level1) == str_level1 # confirm we're back to original topology
    else
        undoinfo = PhyloNetworks.nni!(net_level1, net_level1.edge[13], move);
        # test that move was made
        nodes = [n.number for n in net_level1.edge[17].node]
        if move == 0x03
            @test 8 in nodes # check that edge α connected to v
        else
            @test !(8 in nodes)
        end
        #location of u node (number -11)
        nodes = [n.number for n in net_level1.edge[11].node]
        if move == 0x03 # check that edge γ connected to u
            @test -11 in nodes
        else # check that edge δ connected to u
            @test !(-11 in nodes)
        end
        # check directionality
        @test PhyloNetworks.getChild(net_level1.edge[13]).number == 8
        PhyloNetworks.nni!(undoinfo...); # undo move
        @test writeTopology(net_level1) == str_level1 # confirm we're back to original topology
    end
end # of level1 edge 13: BR directed

@testset "level1 edge 16: BB directed move $move" for move in 0x01:0x02
    # e not hybrid, tree parent:  BB case, 2 NNIs if directed, 8 if undirected
    undoinfo = PhyloNetworks.nni!(net_level1, net_level1.edge[16], move);
    nodes = [n.number for n in net_level1.edge[13].node] #β's connections
    @test -11 in nodes #check β connected to u node (number -11)
    nodes = [n.number for n in net_level1.edge[14].node] #γ's connections
    if move == 1 #check that edge γ connected to v
        @test -12 in nodes
    else #check that edge γ connected to u
        @test -11 in nodes
    end
    # check directionality node -11 child of edge 16 in both cases
    @test PhyloNetworks.getChild(net_level1.edge[16]).number == -11
    PhyloNetworks.nni!(undoinfo...); # undo move
    @test writeTopology(net_level1) == str_level1 #confirm we're back to original topology
end #of level1 edge 16: BB directed

# BB directed has two allowed moves, so 0x03 should throw an exception
@test_throws Exception PhyloNetworks.nni!(net_level1, net_level1.edge[16], 0x03);

@testset "level1 edge 17: BB directed with 4cycle move $move" for move in 0x01:0x02
    # no3cycle test
        # move 2 should fail because of 3cycle problems
        # β is connected in a 4 cycle with γ so a move that makes γ and β a pair
        # would create a 3 cycle
    if move == 0x01
        undoinfo = PhyloNetworks.nni!(net_level1, net_level1.edge[17], move);
        nodes = [n.number for n in net_level1.edge[12].node] #β's connections
        @test -6 in nodes #check β connected to u node (number -6)
        nodes = [n.number for n in net_level1.edge[13].node] #γ's connections
        #check that edge γ connected to v
        @test -11 in nodes
        # check directionality node -6 child of edge 17 in both cases
        @test PhyloNetworks.getChild(net_level1.edge[17]).number == -6
        PhyloNetworks.nni!(undoinfo...); # undo move
        @test writeTopology(net_level1) == str_level1 #confirm we're back to original topology
    else
        @test isnothing(PhyloNetworks.nni!(net_level1, net_level1.edge[17], move))
        undoinfo = PhyloNetworks.nni!(net_level1, net_level1.edge[17], move, false) #should work if we dont check for 3cycles
        nodes = [n.number for n in net_level1.edge[12].node] #β's connections
        @test -6 in nodes #check β connected to u node (number -6)
        nodes = [n.number for n in net_level1.edge[13].node] #γ's connections
        #check that edge γ connected to u
        @test -6 in nodes
        # check directionality node -11 child of edge 17 in both cases
        @test PhyloNetworks.getChild(net_level1.edge[17]).number == -6
        PhyloNetworks.nni!(undoinfo...); # undo move
        @test writeTopology(net_level1) == str_level1 #confirm we're back to original topology
    end
end # of level1 edge 17: BB directed

@testset "level1 edge 18: RR (directed) move $move" for move in 0x01:0x04
    # RB case, 4 moves. uv edge cannot contain the root (always directed)
    undoinfo = PhyloNetworks.nni!(net_level1, net_level1.edge[18], move);
    #test that move was made
    nodes = [n.number for n in net_level1.edge[19].node] #α's connections
    if move in [0x01, 0x03]
        @test -6 in nodes #check that edge α connected to v
    elseif move in [0x02, 0x04]
        @test !(-6 in nodes)
    end
    nodes = [n.number for n in net_level1.edge[12].node] #γ's connections
    if move in [0x03, 0x04] #check that edge γ connected to u (11)
        @test 11 in nodes
    elseif move in [0x01, 0x02]  #check that edge δ connected to u, not γ
        @test !(11 in nodes)
    end
    #check directionality (should point toward u, node 11)
    @test PhyloNetworks.getChild(net_level1.edge[18]).number == 11
    PhyloNetworks.nni!(undoinfo...);
    @test writeTopology(net_level1) == str_level1
end # of level1 edge 18: RR (directed)

@testset "non tree child net edge 3: RB (directed) move $move" for move in 0x01:0x04
    # RB case, 4 moves. uv edge cannot contain the root
    undoinfo = PhyloNetworks.nni!(net_nontreechild, net_nontreechild.edge[3], move);
    #test that move was made
    #location of v node (node -5)
    nodes = [n.number for n in net_nontreechild.edge[4].node] #α's connections
    if move in [0x01, 0x03]
        #check that edge α connected to v
        @test -5 in nodes
    else
        #check that edge β connected to v, not u
        @test !(-5 in nodes)
    end
    nodes = [n.number for n in net_nontreechild.edge[2].node] #δ's connections
    if move in [0x01, 0x02] #check that edge δ connected to u
        @test 3 in nodes
    else
        @test !(3 in nodes)
    end
    #check directionality (should point toward u, node 3)
    @test PhyloNetworks.getChild(net_nontreechild.edge[3]).number == 3
    #undo move
    PhyloNetworks.nni!(undoinfo...);
    #confirm we're back to original topology
    @test writeTopology(net_nontreechild) == str_nontreechild
end #of non tree child net edge 5: RB (directed)

@testset "hybrid ladder net edge 1: BR undirected (u at root, potential nonDAG, 3cycle) move $move" for move in 0x01:0x06
    # BR case, 6 moves. uv edge can contain the root. u at root
    # DAG check: u is root, so if α -> γ, moves 1 and 3 will create a nonDAG
        # α -> γ so moves 1 and 3 will fail.
        # moves 1 and 3 in the notes correspond to 1, 1' and 3, 3' (1, 4 and 3, 6)
    # 3cycle check: could create 3 cycle (α connected to γ) so moves 1 and 5 forbidden
    if move in [0x01, 0x03, 0x05, 0x06] # DAG check
        @test isnothing(PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[1], move, false))
    end
    if move in [0x01, 0x05] # 3cycle check
        @test isnothing(PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[1], move))
    end
    if move == 0x04
        undoinfo = PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[1], move, false);
        nodes = [n.number for n in net_hybridladder.edge[12].node] # α
        @test !(1 in nodes) # check that edge α is connected to v
        nodes = [n.number for n in net_hybridladder.edge[4].node] # δ's connections
        @test -2 in nodes # check that edge δ is connected to u
        # check directionality (edge should point toward u, node -2)
        @test PhyloNetworks.getChild(net_hybridladder.edge[1]).number == 1
        PhyloNetworks.nni!(undoinfo...);
        @test writeTopology(net_hybridladder) == str_hybridladder
    elseif move == 0x02
        undoinfo = PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[1], move);
        nodes = [n.number for n in net_hybridladder.edge[12].node] # α
        @test !(1 in nodes) # check that edge α not connected to v
        nodes = [n.number for n in net_hybridladder.edge[4].node] # δ's connections
        @test -2 in nodes # δ connected to u
        # check directionality (edge should point toward u, node -2)
        @test PhyloNetworks.getChild(net_hybridladder.edge[1]).number == 1
        PhyloNetworks.nni!(undoinfo...);
        @test writeTopology(net_hybridladder) == str_hybridladder
    end
end # of hybrid ladder net edge 1: BR undirected

@testset "hybrid ladder net edge 4: RR (directed) move $move" for move in 0x01:0x02
    # RR case, 2 moves. uv edge cannot contain the root (always directed)
    undoinfo = PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[4], move);
    #test that move was made
    nodes = [n.number for n in net_hybridladder.edge[5].node] #α's connections
    if move == 0x01 #check that edge α connected to v
        @test 4 in nodes
    else move == 0x02 #check that edge β connected to v, not α
        @test !(4 in nodes)
    end
    #check that edge δ connected to u, not γ
    nodes = [n.number for n in net_hybridladder.edge[3].node]
    @test 1 in nodes
    #check directionality (should point toward u, node 1)
    @test PhyloNetworks.getChild(net_hybridladder.edge[4]).number == 1
    #undo move
    PhyloNetworks.nni!(undoinfo...);
    #confirm we're back to original topology
    @test writeTopology(net_hybridladder) == str_hybridladder
end #of hybrid ladder net edge 4: RR (directed)

@testset "hybrid ladder net edge 5: BR undirected move $move" for move in 0x01:0x06
    # BR case, 6 moves. uv edge can contain the root
    # 3 cycle test: α connected to γ, α -> u
    #   moves 1, 5 forbidden
    # DAG test:
    #   no path from α -> γ or β -> γ so all moves should work
    if move in [0x01, 0x05]
        @test isnothing(PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[5], move)) #3cycle forbidden
        undoinfo = PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[5], move, false) #3cycles allowed
    else
        undoinfo = PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[5], move);
    end
    #tests for all moves:
    nodes = [n.number for n in net_hybridladder.edge[6].node] #α
    if move in [0x01, 0x03, 0x05, 0x06]
        @test 1 in nodes #check that edge α connected to v
    else
        @test !(1 in nodes) #check that edge α not connected to v
    end
    nodes = [n.number for n in net_hybridladder.edge[4].node] #δ's connections
    if move in [0x01, 0x02, 0x04, 0x05]
        @test -4 in nodes #δ connected to u
    else
        @test !(-4 in nodes)
    end
    #check directionality
    if move in [0x01, 0x05]
        #(edge should point toward u, node -4)
        @test PhyloNetworks.getChild(net_hybridladder.edge[5]).number == -4
    else
        @test PhyloNetworks.getChild(net_hybridladder.edge[5]).number == 1
    end
    #undo move
    PhyloNetworks.nni!(undoinfo...);
    #confirm we're back to original topology
    @test writeTopology(net_hybridladder) == str_hybridladder
end #of hybrid ladder net edge 5: BR undirected

@testset "hybrid ladder net edge 12: BB undirected (edge below root) move $move" for move in 0x01:0x08
    # BB case, 8 moves. uv edge can contain the root. no flip.
    # no3cycle: moves 1, 4, 5, 8 would create a 3 cycle because α is connected to γ
    if move in [0x01, 0x04, 0x05, 0x08]
        @test isnothing(PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[12], move))
        undoinfo = PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[12], move, false);
    else
        undoinfo = PhyloNetworks.nni!(net_hybridladder, net_hybridladder.edge[12], move);
    end
    nodes = [n.number for n in net_hybridladder.edge[1].node] #α
    if move in [0x01, 0x02, 0x06, 0x08]
        @test -3 in nodes #check that edge α connected to v
    else
        @test !(-3 in nodes) #check that edge α not connected to v
    end
    nodes = [n.number for n in net_hybridladder.edge[11].node] #δ's connections
    if move in [0x01, 0x03, 0x07, 0x08]
        @test -2 in nodes #δ connected to u
    else
        @test !(-2 in nodes)
    end
    #check directionality
    @test PhyloNetworks.getChild(net_hybridladder.edge[12]).number == -3
    #undo move
    PhyloNetworks.nni!(undoinfo...);
    # confirm we're back to original topology
    @test writeTopology(net_hybridladder) == str_hybridladder
end # of hybrid ladder net edge 12: BB undirected (edge below root)

@testset "test isdescendant and isconnected functions" begin
    net_level1 = readTopology(str_level1);
    @test PhyloNetworks.isdescendant(net_level1.node[7], net_level1.node[17]) == true # nodes -9, -6
    @test PhyloNetworks.isdescendant(net_level1.node[7], net_level1.node[3]) == false # nodes -9, -4
    @test PhyloNetworks.isdescendant(net_level1.node[15], net_level1.node[17]) == true # node -12, -6
    @test PhyloNetworks.isconnected(net_level1.node[12], net_level1.node[17]) == true #node -7, -6
    @test PhyloNetworks.isconnected(net_level1.node[12], net_level1.node[19]) == false # node -7, -3
end

#TODO add these after adding contraint functionality
# myconstraint = PhyloNetworks.TopologyConstraint(0x03,Set("8"),net_level1)
# undoinfo = PhyloNetworks.nni!(net_level1, net_level1.edge[1], myconstraint)

end # of testset on unconstrained NNIs


text:
using SeisModels, Test

@testset "IO" begin
    @testset "Mineos" begin
        mktempdir() do dir
            file = joinpath(dir, "model")
            write_mineos(file, AK135)
            ak135 = read_mineos(file)
            @test AK135 ≈ ak135
            @test read_mineos(file) == read_mineos(IOBuffer(String(read(file))))
            let io = IOBuffer()
                write_mineos(io, AK135)
                seekstart(io)
                @test read_mineos(io) ≈ AK135
            end
        end
    end

    @testset "Tvel" begin
        mktempdir() do dir
            file = joinpath(dir, "model.tvel")
            write_tvel(file, AK135)
            @test read_tvel(file) ≈ AK135
            @test read_tvel(file) == read_tvel(IOBuffer(String(read(file))))
            let io = IOBuffer()
                write_tvel(io, AK135)
                seekstart(io)
                @test read_tvel(io) ≈ AK135
            end
            write_tvel(file, AK135, comment1="Line 1", comment2="Line 2")
            lines = readlines(file)
            @test lines[1] == "Line 1"
            @test lines[2] == "Line 2"
        end
    end
end