text:
@testset "ME tools" begin
    @test MatrixExponential <: AbstractMatrixExponential
    @test_throws DimensionMismatch MatrixExponential(ones(1,2),-ones(1,1),ones(1))
    @test_throws DomainError MatrixExponential(ones(1,1),-[1.0 1.0],ones(1))
    @test_throws DomainError MatrixExponential(ones(1,1),-[1.0 1.0],ones(2))
    exp_rv = MatrixExponential(ones(1,1),-ones(1,1),ones(1))
    # pdf
    @test DiscretisedFluidQueues._order(exp_rv)==1
    exp_pdf = pdf(exp_rv)
    @test typeof(exp_pdf(1.0))<:Float64
    r = -log.(rand(10))
    @test exp_pdf.(r)≈exp.(-r)
    @test_throws DomainError pdf([1.0 1.0],exp_rv)
    @test pdf(exp_rv,r)≈pdf(exp_rv).(r)
    @test pdf(exp(-1).*ones(1,1),exp_rv,r)≈pdf(exp_rv).(r.+1.0)
    # ccdf
    exp_ccdf = ccdf(exp_rv)
    @test typeof(exp_ccdf(1.0))<:Float64
    r = -log.(rand(10))
    @test exp_ccdf.(r)≈exp.(-r)
    @test_throws DomainError ccdf([1.0 1.0],exp_rv)
    @test ccdf(exp_rv,r)≈ccdf(exp_rv).(r)
    @test ccdf(exp(-1).*ones(1,1),exp_rv,r)≈ccdf(exp_rv).(r.+1.0)
    # cdf
    exp_cdf = cdf(exp_rv)
    @test typeof(exp_cdf(1.0))<:Float64
    r = -log.(rand(10))
    @test exp_cdf.(r)≈1.0.-exp.(-r)
    @test_throws DomainError cdf([1.0 1.0],exp_rv)
    @test cdf(exp_rv,r)≈cdf(exp_rv).(r)
    # build_me
    me_1 = build_me(cme_params[1])
    @test cdf(me_1,r)≈cdf(exp_rv,r)
    #
    me_3 = build_me(cme_params[3])
    @test typeof(me_3)<:AbstractMatrixExponential
    @test typeof(me_3)<:ConcentratedMatrixExponential
    me_3_cdf = cdf(me_3)
    @test typeof(me_3_cdf(1.0))<:Float64
    for f in fieldnames(ConcentratedMatrixExponential)
        @eval begin
            me_3 = build_me(cme_params[3])
            @test ConcentratedMatrixExponential(3).$f==me_3.$f
        end
    end
    @test -sum(me_3.a/me_3.S)≈1.0
    # build_erlang
    erl_1 = DiscretisedFluidQueues.build_erlang(1)
    @test typeof(erl_1)<:AbstractMatrixExponential
    @test typeof(erl_1)<:MatrixExponential
    @test ccdf(erl_1,r)≈ccdf(exp_rv,r)
    erl_3 = DiscretisedFluidQueues.build_erlang(3)
    @test -sum(erl_3.a/erl_3.S)≈1.0
    e1 = zeros(1,3)
    e1[1]=1.0
    @test erl_3.a==e1
    @test erl_3.S≈[-3.0 3.0 0.0; 0.0 -3.0 3.0; 0.0 0.0 -3.0]
    @test erl_3.s≈3*e1[end:-1:1]
    # orbits
    exp_rv_orbit_fun = DiscretisedFluidQueues.orbit(exp_rv)
    @test exp_rv_orbit_fun(2.0)≈[1.0]

    me_3_orbit_fun = DiscretisedFluidQueues.orbit(me_3)
    me_3_mean_2 = ConcentratedMatrixExponential(3; mean=2.0)
    me_3_mean_2_orbit_fun = DiscretisedFluidQueues.orbit(me_3_mean_2)
    @test me_3_mean_2_orbit_fun.(r)≈me_3_orbit_fun.(0.5*r)
    me_3_not_CME_Type = MatrixExponential(me_3.a,Matrix(me_3.S),me_3.s,me_3.D)
    me_3_not_CME_Type_orbit_fun = DiscretisedFluidQueues.orbit(me_3_not_CME_Type)
    @test me_3_orbit_fun.(r)≈me_3_not_CME_Type_orbit_fun.(r)
    me_3_mean_2_not_CME_Type =
        MatrixExponential(me_3_mean_2.a,Matrix(me_3_mean_2.S),me_3_mean_2.s,me_3_mean_2.D)
    me_3_mean_2_not_CME_Type_orbit_fun = DiscretisedFluidQueues.orbit(me_3_mean_2_not_CME_Type)
    @test me_3_mean_2_not_CME_Type_orbit_fun.(r)≈me_3_not_CME_Type_orbit_fun.(0.5*r)
    # Expected orbits
    ω = abs(me_3.S[3,2])
    period = 2*pi/ω
    exp_orbit =
        DiscretisedFluidQueues.expected_orbit_from_pdf(
            x->ccdf(me_3,x)/-sum(me_3.a/me_3.S),
            me_3,0.0,period,10000)/(1-exp(me_3.S[1,1]*period))
    stationary = (me_3.a/me_3.S)/sum(me_3.a/me_3.S)
    @test exp_orbit≈stationary atol=1e-6
    exp_orbit_from_cdf =
        DiscretisedFluidQueues.expected_orbit_from_cdf(
            x->cdf(stationary,me_3,x)/-sum(me_3.a/me_3.S),
            me_3,0.0,period,10000)/(1-exp(me_3.S[1,1]*period))
    @test exp_orbit_from_cdf≈stationary atol=1e-6
end


text:
using MeasureTheory
using Base.Iterators
using Statistics
using Soss
using Random

# 5 7 10

rng = Random.Xoshiro(12)

m = @model begin
    latent ~ Chain(Normal(μ=1)) do x Normal(μ=x, σ=0.2) end
    observed ~ For(latent) do x Poisson(logλ=x) end
end

truth = rand(rng, m())

xv = take(truth.latent, 100) |> collect
yv = take(truth.observed, 100) |> collect

using Plots

plt = scatter(yv, label="observations")
plot!(exp.(xv), lw=3, label="latent process")

using Statistics

scatter(exp.(xv), poiscdf.(exp.(xv), yv), label=false)

x = Chain(Normal(μ=1)) do xj Normal(μ=xj, σ=0.2) end

xvals = rand(rng, x)
y = For(xobs) do xj Poisson(logλ= xj) end
yvals = rand(rng, y)
xv = take(xobs, 100) |> collect
yv = take(yobs, 100) |> collect

take(xobs.parent, 10) |> collect
take(yobs.parent, 10) |> collect

exp.(xv)
yv

# plt = scatter(normcdf.(xv, 1, yv), label=false)

# for j in 1:10
#     xobs = rand(rng, x)
#     yobs = rand(rng, y)
#     xv = take(xobs, 100) |> collect;
#     yv = take(yobs, 100) |> collect;
#     plt = scatter!(plt, normcdf.(xv, 1, yv), label=false)
# end
# plt

using Soss

m = @model begin
    x ~ Chain(Normal()) do xj Normal(μ=xj) end
    y ~ For(xobs) do xj Poisson(logλ=xj) end
end

truth = rand(rng, m())

xobs = take(truth.x, 10) |> collect
yobs = take(truth.y, 10) |> collect

logdensity_def(m(), (x=xobs, y=yobs))

ext: jl


text:
using ConcurrentCalls
using Base.Test

@testset "RemoteCall" begin
    @test_throws ArgumentError serialize(IOBuffer(), ConcurrentCalls.RemoteCall(@task 1))
end

@testset begin
    scheduler = ConcurrentCalls.Scheduler()
    a = ConcurrentCalls.call(scheduler, +, 1, 1)
    ConcurrentCalls.teardown(scheduler)
    @test fetch(a) == 2
    @test sprint(show, scheduler.stats) isa String
end

@testset begin
    scheduler = ConcurrentCalls.Scheduler()
    a = ConcurrentCalls.call(scheduler, +, 1, 1)
    b = ConcurrentCalls.call(scheduler, +, 1, 2)
    c = ConcurrentCalls.call(scheduler, +, a, b)
    ConcurrentCalls.teardown(scheduler)
    @test fetch(c) == 5
    @test sprint(show, scheduler.stats) isa String
end

@testset begin
    scheduler = ConcurrentCalls.Scheduler()
    a = ConcurrentCalls.call(scheduler, +, 0, 0)
    for _ in 1:5
        a = ConcurrentCalls.call(scheduler, +, a, 1)
    end
    ConcurrentCalls.teardown(scheduler)
    @test fetch(a) == 5
    @test sprint(show, scheduler.stats) isa String
end

@testset "pmap" begin
    scheduler = ConcurrentCalls.Scheduler()
    a = ConcurrentCalls.pmap(scheduler, x->2x, collect(1:100))
    ConcurrentCalls.teardown(scheduler)
    @test fetch.(a) == map(x->2x, collect(1:100))

    scheduler = ConcurrentCalls.Scheduler()
    a = ConcurrentCalls.pmap(scheduler, x->2x, collect(1:100))
    ConcurrentCalls.teardown(scheduler)
    @test fetch.(a) == map(x->2x, collect(1:100))

    scheduler = ConcurrentCalls.Scheduler()
    a = ConcurrentCalls.pmap(scheduler, x->2x, (i for i in 1:100 if rand() ≥ 0))
    ConcurrentCalls.teardown(scheduler)
    @test fetch.(a) == map(x->2x, collect(1:100))
end

@testset "ChunkingGenerator" begin
    scheduler = ConcurrentCalls.Scheduler()
    gen = ConcurrentCalls.ChunkingGenerator(scheduler, 1:100, x->x/2)
    @test collect(gen) == map(x->x/2, 1:100)
    ConcurrentCalls.teardown(scheduler)

    scheduler = ConcurrentCalls.Scheduler()
    gen = ConcurrentCalls.ChunkingGenerator(scheduler, 1:100, identity)
    result = mapreduce(x->2x, +, gen)
    ConcurrentCalls.teardown(scheduler)
    @test result === mapreduce(x->2x, +, 1:100)
end

@testset begin
    # fetch before teardown
    scheduler = ConcurrentCalls.Scheduler()
    ret = ConcurrentCalls.call(scheduler, +, 1, 1)
    val = fetch(ret)
    ConcurrentCalls.teardown(scheduler)
    @test val == 2
end

@testset "@cc" begin
    val = @cc 1
    @test val == 1

    val = @cc 1 + 1
    @test val == 2

    val = @cc begin
        1 + 2
    end
    @test val == 3

    val = @cc begin
        a = 1 + 1
        a * 2
    end
    @test val == 4

    x = 4
    val = @cc begin
        x + 1
    end
    @test val == 5

    x = 2
    val = @cc begin
        (y -> x + y)(4)
    end
    @test val == 6

    val = @cc begin
        a = pmap(1:3) do x
            x
        end
        a[1] + a[2] + a[3] + 1
    end
    @test val == 7

    val = @cc pmap(x->2x, 1:3)
    @test val == [2,4,6]

    #=
    val = @cc begin
        a = [i * 2 for i in 1:10]
        sum(a)
    end
    @test val == 110
    =#

    val = @cc (x + 1 for x in 1:3)
    @test val == [2,3,4]

    val = @cc [x + 1 for x in 1:3]
    @test val == [2,3,4]
end

#=
@testset begin
    foo(x) = x + 1
    bar(x) = 2x
    baz(x) = 3x
    xs = [1,2,3,4]
    val = @cc begin
        xs = [foo(x) for x in xs]
        xs = [bar(x) for x in xs]
        xs = [baz(x) for x in xs]
        sum(xs)
    end
    @test val == 84
end
=#


text:
@testset "Base Extensions" begin
    # Define test model and pref
    m = InfiniteModel()
    m2 = InfiniteModel()
    pref = ParameterRef(m, 1)
    # test variable compare
    @testset "(==)" begin
        @test pref == pref
        @test pref == ParameterRef(m, 1)
        @test !(pref == ParameterRef(m, 2))
        @test !(pref == ParameterRef(m2, 1))
        @test !(pref != ParameterRef(m, 1))
    end
    # test copy(v)
    @testset "copy(v)" begin
        @test copy(pref) == pref
    end
    # test copy(v, m)
    @testset "copy(v, m)" begin
        @test copy(pref, m2) == ParameterRef(m2, 1)
    end
    # test broadcastable
    @testset "broadcastable" begin
        @test isa(Base.broadcastable(pref), Base.RefValue{ParameterRef})
    end
end

# Test macro methods
@testset "Macro Helpers" begin
    @testset "Symbol Methods" begin
        @test Parameter == :Parameter
        @test InfiniteOpt._is_set_keyword(:(lower_bound = 0))
    end
    # test ParameterInfoExpr datatype
    @testset "_ParameterInfoExpr" begin
        @test InfiniteOpt._ParameterInfoExpr isa DataType
        @test InfiniteOpt._ParameterInfoExpr(ones(Bool, 8)...).has_lb
        @test !InfiniteOpt._ParameterInfoExpr().has_lb
    end
    # test JuMP._set_lower_bound_or_error
    @testset "JuMP._set_lower_bound_or_error" begin
        info = InfiniteOpt._ParameterInfoExpr()
        # test normal operation
        @test isa(JuMP._set_lower_bound_or_error(error, info, 0), Nothing)
        @test info.has_lb && info.lower_bound == 0
        # test double/lack of input errors
        @test_throws ErrorException JuMP._set_lower_bound_or_error(error,
                                                                   info, 0)
        info.has_lb = false; info.has_dist = true
        @test_throws ErrorException JuMP._set_lower_bound_or_error(error,
                                                                   info, 0)
        info.has_dist = false; info.has_set = true
        @test_throws ErrorException JuMP._set_lower_bound_or_error(error,
                                                                   info, 0)
    end
    # test JuMP._set_upper_bound_or_error
    @testset "JuMP._set_upper_bound_or_error" begin
        info = InfiniteOpt._ParameterInfoExpr()
        # test normal operation
        @test isa(JuMP._set_upper_bound_or_error(error, info, 0), Nothing)
        @test info.has_ub && info.upper_bound == 0
        # test double/lack of input errors
        @test_throws ErrorException JuMP._set_upper_bound_or_error(error,
                                                                   info, 0)
        info.has_ub = false; info.has_dist = true
        @test_throws ErrorException JuMP._set_upper_bound_or_error(error,
                                                                   info, 0)
        info.has_dist = false; info.has_set = true
        @test_throws ErrorException JuMP._set_upper_bound_or_error(error,
                                                                   info, 0)
    end
    # _dist_or_error
    @testset "_dist_or_error" begin
        info = InfiniteOpt._ParameterInfoExpr()
        # test normal operation
        @test isa(InfiniteOpt._dist_or_error(error, info, 0), Nothing)
        @test info.has_dist && info.distribution == 0
        # test double/lack of input errors
        @test_throws ErrorException InfiniteOpt._dist_or_error(error, info, 0)
        info.has_dist = false; info.has_lb = true
        @test_throws ErrorException InfiniteOpt._dist_or_error(error, info, 0)
        info.has_lb = false; info.has_set = true
        @test_throws ErrorException InfiniteOpt._dist_or_error(error, info, 0)
    end
    # _set_or_error
    @testset "_set_or_error" begin
        info = InfiniteOpt._ParameterInfoExpr()
        # test normal operation
        @test isa(InfiniteOpt._set_or_error(error, info, 0), Nothing)
        @test info.has_set && info.set == 0
        # test double/lack of input errors
        @test_throws ErrorException InfiniteOpt._set_or_error(error, info, 0)
        info.has_set = false; info.has_lb = true
        @test_throws ErrorException InfiniteOpt._set_or_error(error, info, 0)
        info.has_lb = false; info.has_dist = true
        @test_throws ErrorException InfiniteOpt._set_or_error(error, info, 0)
    end
    # _constructor_set
    @testset "InfiniteOpt._constructor_set" begin
        info = InfiniteOpt._ParameterInfoExpr()
        @test_throws ErrorException InfiniteOpt._constructor_set(error, info)
        info.has_lb = true; info.lower_bound = 0
        @test_throws ErrorException InfiniteOpt._constructor_set(error, info)
        info.has_ub = true; info.upper_bound = 1
        check = :(isa($(info.lower_bound), Number))
        expected = :($(check) ? IntervalSet($(info.lower_bound), $(info.upper_bound)) : error("Bounds must be a number."))
        @test InfiniteOpt._constructor_set(error, info) == expected
        info = InfiniteOpt._ParameterInfoExpr(distribution = Normal())
        check = :(isa($(info.distribution), Distributions.NonMatrixDistribution))
        expected = :($(check) ? DistributionSet($(info.distribution)) : error("Distribution must be a subtype of Distributions.NonMatrixDistribution."))
        @test InfiniteOpt._constructor_set(error, info) == expected
        info = InfiniteOpt._ParameterInfoExpr(set = IntervalSet(0, 1))
        check1 = :(isa($(info.set), AbstractInfiniteSet))
        check2 = :(isa($(info.set), Distributions.NonMatrixDistribution))
        expected = :($(check1) ? $(info.set) : ($(check2) ? DistributionSet($(info.set)) : error("Set must be a subtype of AbstractInfiniteSet.")))
        @test InfiniteOpt._constructor_set(error, info) == expected
    end
    # _parse_one_operator_parameter
    @testset "_parse_one_operator_parameter" begin
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_one_operator_parameter(error, info,
                                                          Val(:<=), 0), Nothing)
        @test info.has_ub && info.upper_bound == 0
        @test isa(InfiniteOpt._parse_one_operator_parameter(error, info,
                                                          Val(:>=), 0), Nothing)
        @test info.has_lb && info.lower_bound == 0
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_one_operator_parameter(error, info,
                                                          Val(:in), esc(0)), Nothing)
        @test info.has_set && info.set == esc(0)
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_one_operator_parameter(error, info,
                                                          Val(:in), esc(:([0, 1]))), Nothing)
        @test info.has_lb && info.has_ub
        @test_throws ErrorException InfiniteOpt._parse_one_operator_parameter(error, info,
                                                                              Val(:d), 0)
    end
    # _parse_ternary_parameter
    @testset "_parse_ternary_parameter" begin
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_ternary_parameter(error, info, Val(:<=), 0,
                                                       Val(:<=), 0), Nothing)
        @test info.has_ub && info.upper_bound == 0
        @test info.has_lb && info.lower_bound == 0
        info = InfiniteOpt._ParameterInfoExpr()
        @test isa(InfiniteOpt._parse_ternary_parameter(error, info, Val(:>=), 0,
                                                       Val(:>=), 0), Nothing)
        @test info.has_ub && info.upper_bound == 0
        @test info.has_lb && info.lower_bound == 0
        @test_throws ErrorException InfiniteOpt._parse_ternary_parameter(error,
                                                 info, Val(:<=), 0, Val(:>=), 0)
    end
    # _parse_parameter
    @testset "_parse_parameter" begin
        info = InfiniteOpt._ParameterInfoExpr()
        @test InfiniteOpt._parse_parameter(error, info, :<=, :x, 0) == :x
        @test info.has_ub && info.upper_bound == 0
        @test InfiniteOpt._parse_parameter(error, info, :<=, 0, :x) == :x
        @test info.has_lb && info.lower_bound == 0
        info = InfiniteOpt._ParameterInfoExpr()
        @test InfiniteOpt._parse_parameter(error, info, 0, :<=, :x, :<=, 0) == :x
        @test info.has_ub && info.upper_bound == 0
        @test info.has_lb && info.lower_bound == 0
    end
end

# Test precursor functions needed for add_parameter
@testset "Basic Reference Queries" begin
    m = InfiniteModel()
    pref = ParameterRef(m, 1)
    # JuMP.index
    @testset "JuMP.index" begin
        @test JuMP.index(pref) == 1
    end
    # JuMP.owner_model
    @testset "JuMP.owner_model" begin
        @test owner_model(pref) == m
    end
end

# Test name methods
@testset "Name" begin
    m = InfiniteModel()
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    m.params[1] = param
    m.param_to_name[1] = "test"
    m.name_to_param = nothing
    pref = ParameterRef(m, 1)
    # JuMP.name
    @testset "JuMP.name" begin
        @test name(pref) == "test"
    end
    # JuMP.set_name
    @testset "JuMP.set_name" begin
        @test isa(set_name(pref, "new"), Nothing)
        @test name(pref) == "new"
    end
    # parameter_by_name
    @testset "parameter_by_name" begin
        @test parameter_by_name(m, "new") == pref
        @test isa(parameter_by_name(m, "test2"), Nothing)
        m.params[2] = param
        m.param_to_name[2] = "new"
        m.name_to_param = nothing
        @test_throws ErrorException parameter_by_name(m, "new")
    end
end

# Test parameter definition methods
@testset "Definition" begin
    # _check_supports_in_bounds
    @testset "_check_supports_in_bounds" begin
        set = IntervalSet(0, 1)
        @test isa(InfiniteOpt._check_supports_in_bounds(error, 0, set), Nothing)
        @test_throws ErrorException InfiniteOpt._check_supports_in_bounds(error,
                                                                        -1, set)
        @test_throws ErrorException InfiniteOpt._check_supports_in_bounds(error,
                                                                         2, set)
        set = DistributionSet(Uniform())
        @test isa(InfiniteOpt._check_supports_in_bounds(error, 0, set), Nothing)
        @test_throws ErrorException InfiniteOpt._check_supports_in_bounds(error,
                                                                        -1, set)
        @test_throws ErrorException InfiniteOpt._check_supports_in_bounds(error,
                                                                         2, set)
    end
    # build_parameter
    @testset "build_parameter" begin
        set = DistributionSet(Multinomial(3, [1/2, 1/2]))
        supps = Vector(0:1)
        expected = InfOptParameter(set, supps, false)
        @test build_parameter(error, set, 2, supports = supps) == expected
        expected = InfOptParameter(set, Number[], true)
        @test build_parameter(error, set, 2, independent = true) == expected
        @test_throws ErrorException build_parameter(error, set, 2, bob = 42)
        @test_throws ErrorException build_parameter(error, set, 1)
        warn = "Support points are not unique, eliminating redundant points."
        @test_logs (:warn, warn) build_parameter(error, set, 2,
                                                 supports = ones(3),
                                                 independent = true)
        warn = "Ignoring num_supports since supports is not empty."
        @test_logs (:warn, warn) build_parameter(error, set, 2, supports = [1, 2],
                                                 num_supports = 1)
        @test_throws ErrorException build_parameter(error, set, 2,
                                                    num_supports = 3)
        repeated_supps = [1, 1]
        expected = InfOptParameter(set, [1, 1], false)
        @test build_parameter(error, set, 2, supports = repeated_supps) == expected
        expected = InfOptParameter(set, [1], true)
        warn = "Support points are not unique, eliminating redundant points."
        @test_logs (:warn, warn) build_parameter(error, set, 2,
                       supports = repeated_supps,independent = true) == expected
        set = IntervalSet(0, 1)
        expected = InfOptParameter(set, [0., 0.5, 1.], false)
        @test build_parameter(error, set, num_supports = 3) == expected
    end
    # _check_supports_dimensions
    @testset "_check_supports_dimensions" begin
        model = InfiniteModel()
        param1 = InfOptParameter(IntervalSet(0, 1), Number[], false)
        param2 = InfOptParameter(IntervalSet(0, 1), [0], false)
        model.param_to_group_id[1] = 1
        model.param_to_group_id[2] = 1
        model.params[1] = param1
        model.params[2] = param2
        model.next_param_id = 1
        @test InfiniteOpt._check_supports_dimensions(model, param1, 1) isa Nothing
        @test_throws ErrorException InfiniteOpt._check_supports_dimensions(model, param2, 2)
        model.param_to_group_id[3] = 2
        model.params[3] = param1
        model.next_param_id = 2
        @test InfiniteOpt._check_supports_dimensions(model, param1, 3) isa Nothing
        model.params[1] = param2
        @test InfiniteOpt._check_supports_dimensions(model, param2, 2) isa Nothing
    end
    # add_parameter
    @testset "add_parameter" begin
        m = InfiniteModel()
        param = InfOptParameter(IntervalSet(0, 1), Number[], false)
        expected = ParameterRef(m, 1)
        @test add_parameter(m, param) == expected
        @test m.params[1] isa InfOptParameter
        @test m.param_to_group_id[1] == 1
    end
end

# Test the parameter macro
@testset "Macro" begin
    m = InfiniteModel()
    # single parameter
    @testset "Single" begin
        pref = ParameterRef(m, 1)
        @test @infinite_parameter(m, 0 <= a <= 1) == pref
        @test m.params[1].set == IntervalSet(0, 1)
        @test name(pref) == "a"
        pref = ParameterRef(m, 2)
        @test @infinite_parameter(m, b in Normal(), supports = [1; 2]) == pref
        @test m.params[2].set == DistributionSet(Normal())
        @test m.params[2].supports == [1, 2]
        pref = ParameterRef(m, 3)
        @test @infinite_parameter(m, c in IntervalSet(0, 1)) == pref
        @test m.params[3].set == IntervalSet(0, 1)
        pref = ParameterRef(m, 4)
        @test @infinite_parameter(m, set = IntervalSet(0, 1),
                                  base_name = "d") == pref
        @test name(pref) == "d"
        pref = ParameterRef(m, 5)
        @test @infinite_parameter(m, set = IntervalSet(0, 1)) == pref
        @test name(pref) == ""
        pref = ParameterRef(m, 6)
        @test @infinite_parameter(m, z in [0, 1]) == pref
        @test m.params[6].set == IntervalSet(0, 1)
        @test name(pref) == "z"
    end
    # multiple parameters
    @testset "Array" begin
        prefs = [ParameterRef(m, 7), ParameterRef(m, 8)]
        @test @infinite_parameter(m, 0 <= e[1:2] <= 1) == prefs
        @test m.params[7].set == IntervalSet(0, 1)
        @test m.params[8].set == IntervalSet(0, 1)
        prefs = [ParameterRef(m, 9), ParameterRef(m, 10)]
        @test @infinite_parameter(m, [1:2], set = IntervalSet(0, 1)) == prefs
        @test m.params[8].set == IntervalSet(0, 1)
        @test m.params[10].set == IntervalSet(0, 1)
        prefs = [ParameterRef(m, 11), ParameterRef(m, 12)]
        sets = [IntervalSet(0, 1), IntervalSet(-1, 2)]
        @test @infinite_parameter(m, f[i = 1:2], set = sets[i]) == prefs
        @test m.params[11].set == IntervalSet(0, 1)
        @test m.params[12].set == IntervalSet(-1, 2)
        prefs = [ParameterRef(m, 13), ParameterRef(m, 14)]
        @test @infinite_parameter(m, [i = 1:2], set = sets[i]) == prefs
        @test m.params[13].set == IntervalSet(0, 1)
        @test m.params[14].set == IntervalSet(-1, 2)
        prefs = [ParameterRef(m, 15), ParameterRef(m, 16)]
        @test @infinite_parameter(m, [0, -1][i] <= g[i = 1:2] <= [1, 2][i]) == prefs
        @test m.params[15].set == IntervalSet(0, 1)
        @test m.params[16].set == IntervalSet(-1, 2)
        prefs = [ParameterRef(m, 17), ParameterRef(m, 18)]
        @test @infinite_parameter(m, 0 <= h[1:2] <= 1,
                                  independent = true) == prefs
        @test m.params[17].independent
        @test m.params[18].independent
        prefs = [ParameterRef(m, 19), ParameterRef(m, 20)]
        prefs = convert(JuMP.Containers.SparseAxisArray, prefs)
        @test @infinite_parameter(m, 0 <= i[1:2] <= 1,
                                  container = SparseAxisArray) == prefs
        @test m.params[19].set == IntervalSet(0, 1)
        @test m.params[20].set == IntervalSet(0, 1)
    end
    # test for errors
    @testset "Errors" begin
        @test_macro_throws ErrorException @infinite_parameter(m, 0 <= [1:2] <= 1)
        @test_macro_throws ErrorException @infinite_parameter(m, 0 <= "bob" <= 1)
        @test_macro_throws ErrorException @infinite_parameter(m, 0 <= a <= 1)
        @test_macro_throws ErrorException @infinite_parameter(m, 0 <= j)
        @test_macro_throws ErrorException @infinite_parameter(m, j)
        @test_macro_throws ErrorException @infinite_parameter(m, j, foo = 42)
        @test_macro_throws ErrorException @infinite_parameter(m, j in Multinomial(3, [1/3, 1/3]))
        @test_macro_throws ErrorException @infinite_parameter(m, 0 <= k <= 1, Int)
    end
end

# Test if used
@testset "Used" begin
    m = InfiniteModel()
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    pref = add_parameter(m, param, "test")
    # used_by_constraint
    @testset "used_by_constraint" begin
        @test !used_by_constraint(pref)
        m.param_to_constrs[JuMP.index(pref)] = [1]
        @test used_by_constraint(pref)
        delete!(m.param_to_constrs, JuMP.index(pref))
    end
    # used_by_measure
    @testset "used_by_measure" begin
        @test !used_by_measure(pref)
        m.param_to_meas[JuMP.index(pref)] = [1]
        @test used_by_measure(pref)
        delete!(m.param_to_meas, JuMP.index(pref))
    end
    # used_by_measure
    @testset "used_by_variable" begin
        @test !used_by_variable(pref)
        m.param_to_vars[JuMP.index(pref)] = [1]
        @test used_by_variable(pref)
        delete!(m.param_to_vars, JuMP.index(pref))
    end
    # is_used
    @testset "is_used" begin
        @test !is_used(pref)
        m.param_to_constrs[JuMP.index(pref)] = [1]
        @test is_used(pref)
        delete!(m.param_to_constrs, JuMP.index(pref))
        m.param_to_meas[JuMP.index(pref)] = [1]
        @test is_used(pref)
        delete!(m.param_to_meas, JuMP.index(pref))
        m.param_to_vars[JuMP.index(pref)] = [1]
        @test is_used(pref)
        delete!(m.param_to_vars, JuMP.index(pref))
    end
end

# Test parameter set methods
@testset "Infinite Set" begin
    m = InfiniteModel()
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    pref = add_parameter(m, param, "test")
    m.param_to_vars[JuMP.index(pref)] = [1]
    # _parameter_set
    @testset "_parameter_set" begin
        @test InfiniteOpt._parameter_set(pref) == IntervalSet(0, 1)
    end
    # _update_parameter_set
    @testset "_update_parameter_set " begin
        @test isa(InfiniteOpt._update_parameter_set(pref,
                                                    IntervalSet(1, 2)), Nothing)
        @test InfiniteOpt._parameter_set(pref) == IntervalSet(1, 2)
    end
    # infinite_set
    @testset "infinite_set" begin
        @test infinite_set(pref) == IntervalSet(1, 2)
    end
    # set_infinite_set
    @testset "set_infinite_set" begin
        @test isa(set_infinite_set(pref, IntervalSet(1, 3)), Nothing)
        @test infinite_set(pref) == IntervalSet(1, 3)
    end
end

# Test parameter support methods
@testset "Supports" begin
    m = InfiniteModel()
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    pref = add_parameter(m, param, "test")
    prefs = @infinite_parameter(m, x[1:3] in [0, 1])
    # _parameter_supports
    @testset "_parameter_supports" begin
        @test InfiniteOpt._parameter_supports(pref) == Number[]
    end
    # _update_parameter_supports
    @testset "_update_parameter_supports " begin
        @test isa(InfiniteOpt._update_parameter_supports(pref, [1]), Nothing)
        @test InfiniteOpt._parameter_supports(pref) == [1]
    end
    # num_supports
    @testset "num_supports" begin
        @test num_supports(pref) == 1
    end
    # has_supports
    @testset "has_supports" begin
        @test has_supports(pref)
        InfiniteOpt._update_parameter_supports(pref, Number[])
        @test !has_supports(pref)
    end
    # supports
    @testset "supports" begin
        @test_throws ErrorException supports(pref)
        InfiniteOpt._update_parameter_supports(pref, [1])
        @test supports(pref) == [1]
    end
    # set_supports
    @testset "set_supports" begin
        @test isa(set_supports(pref, [0, 1], force = true), Nothing)
        @test supports(pref) == [0, 1]
        @test_throws ErrorException set_supports(pref, [2, 3])
        warn = "Support points are not unique, eliminating redundant points."
        @test_logs (:warn, warn) set_supports(pref, [1, 1], force = true)
        @test_throws ErrorException set_supports(pref, [0.5])
        @test set_supports(prefs[1], [1, 1]) isa Nothing
        @test supports(prefs[1]) == [1, 1]
    end
    # add_supports
    @testset "add_supports" begin
        @test isa(add_supports(pref, 0.5), Nothing)
        @test supports(pref) == [1, 0.5]
        @test isa(add_supports(pref, [0, 0.25, 1]), Nothing)
        @test supports(pref) == [1, 0.5, 0, 0.25]
        @test add_supports(prefs[2], [1, 1]) isa Nothing
        @test supports(prefs[2]) == [1, 1]
    end
    # delete_supports
    @testset "delete_supports" begin
        @test isa(delete_supports(pref), Nothing)
        @test_throws ErrorException supports(pref)
    end
end

# Test lower bound functions
@testset "Lower Bound" begin
    m = InfiniteModel()
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    pref = add_parameter(m, param, "test")
    # JuMP.has_lower_bound
    @testset "JuMP.has_lower_bound" begin
        @test has_lower_bound(pref)
        set_infinite_set(pref, DistributionSet(Normal()))
        @test has_lower_bound(pref)
        set_infinite_set(pref, DistributionSet(Multinomial(3, [1/2, 1/2])))
        @test !has_lower_bound(pref)
        set_infinite_set(pref, BadSet())
        @test !has_lower_bound(pref)
    end
    # JuMP.lower_bound
    @testset "JuMP.lower_bound" begin
        set_infinite_set(pref, DistributionSet(Multinomial(3, [1/2, 1/2])))
        @test_throws ErrorException lower_bound(pref)
        set_infinite_set(pref, IntervalSet(0, 1))
        @test lower_bound(pref) == 0
        set_infinite_set(pref, DistributionSet(Normal()))
        @test lower_bound(pref) == -Inf
    end
    # JuMP.set_lower_bound
    @testset "JuMP.set_lower_bound" begin
        @test_throws ErrorException set_lower_bound(pref, 2)
        set_infinite_set(pref, BadSet())
        @test_throws ErrorException set_lower_bound(pref, 2)
        set_infinite_set(pref, IntervalSet(0, 1))
        @test isa(set_lower_bound(pref, -1), Nothing)
        @test lower_bound(pref) == -1
    end
end

# Test upper bound functions
@testset "Upper Bound" begin
    m = InfiniteModel()
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    pref = add_parameter(m, param, "test")
    # JuMP.has_upper_bound
    @testset "JuMP.has_upper_bound" begin
        @test has_upper_bound(pref)
        set_infinite_set(pref, DistributionSet(Normal()))
        @test has_upper_bound(pref)
        set_infinite_set(pref, DistributionSet(Multinomial(3, [1/2, 1/2])))
        @test !has_upper_bound(pref)
        set_infinite_set(pref, BadSet())
        @test !has_upper_bound(pref)
    end
    # JuMP.upper_bound
    @testset "JuMP.upper_bound" begin
        set_infinite_set(pref, DistributionSet(Multinomial(3, [1/2, 1/2])))
        @test_throws ErrorException upper_bound(pref)
        set_infinite_set(pref, IntervalSet(0, 1))
        @test upper_bound(pref) == 1
        set_infinite_set(pref, DistributionSet(Normal()))
        @test upper_bound(pref) == Inf
    end
    # JuMP.set_upper_bound
    @testset "JuMP.set_upper_bound" begin
        @test_throws ErrorException set_upper_bound(pref, 2)
        set_infinite_set(pref, BadSet())
        @test_throws ErrorException set_upper_bound(pref, 2)
        set_infinite_set(pref, IntervalSet(0, 1))
        @test isa(set_upper_bound(pref, 2), Nothing)
        @test upper_bound(pref) == 2
    end
end

# Test the independent manipulation functions
@testset "Independent" begin
    m = InfiniteModel()
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    pref = add_parameter(m, param, "test")
    # is_independent
    @testset "is_indepentent" begin
        @test !is_independent(pref)
    end
    # set_independent
    @testset "set_indepentent" begin
        @test isa(set_independent(pref), Nothing)
        @test is_independent(pref)
    end
    # unset_independent
    @testset "unset_indepentent" begin
        @test isa(unset_independent(pref), Nothing)
        @test !is_independent(pref)
    end
end

# Test the internal helper functions
@testset "Internal Helpers" begin
    m = InfiniteModel()
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    pref = add_parameter(m, param, "test")
    pref2 = add_parameter(m, param, "θ")
    prefs = @infinite_parameter(m, x[1:2], set = IntervalSet(0, 1),
                                container = SparseAxisArray)
    # group_id
    @testset "group_id" begin
        @test group_id(pref) == 1
    end
    # _group (Parameter)
    @testset "_group (Parameter)" begin
        @test InfiniteOpt._group(pref) == 1
    end
    # _group (Array)
    @testset "_group (Array)" begin
        @test InfiniteOpt._group(prefs) == 3
    end
    # group_id (array)
    @testset "group_id (array)" begin
        @test group_id(prefs) == 3
        @test_throws ErrorException group_id([pref; pref2])
    end
    # _group (over Tuple)
    @testset "_group (Tuple)" begin
        @test InfiniteOpt._group.((prefs, pref)) == (3, 1)
    end
    # _only_one_group
    @testset "_only_one_group" begin
        @test InfiniteOpt._only_one_group(pref)
        @test InfiniteOpt._only_one_group(prefs)
        @test !InfiniteOpt._only_one_group(convert(JuMP.Containers.SparseAxisArray,
                                                                 [pref2; pref]))
    end
    # _root_name
    @testset "_root_name" begin
        @test InfiniteOpt._root_name(pref) == "test"
        @test InfiniteOpt._root_name(prefs[1]) == "x"
        @test InfiniteOpt._root_name(pref2) == "θ"
    end
    # _root_names
    @testset "_root_names" begin
        @test InfiniteOpt._root_names((pref, prefs, pref2)) == ("test", "x", "θ")
    end
    # _list_parameter_refs
    @testset "_list_parameter_refs" begin
        @test InfiniteOpt._list_parameter_refs((pref, prefs)) isa Vector
        result = InfiniteOpt._list_parameter_refs((pref, prefs))
        @test length(result) == 3
        @test result[1] == pref
        @test result[2] == prefs[1] || result[2] == prefs[2]
        @test result[3] == prefs[1] || result[3] == prefs[2]
    end
end

# Test everything else
@testset "Other Queries" begin
    m = InfiniteModel()
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    pref = add_parameter(m, param, "test")
    prefs = @infinite_parameter(m, [1:2], set = IntervalSet(0, 1),
                                container = SparseAxisArray)
    # JuMP.is_valid
    @testset "JuMP.is_valid" begin
        @test is_valid(m, pref)
        pref2 = ParameterRef(InfiniteModel(), 1)
        @test !is_valid(m, pref2)
        pref3 = ParameterRef(m, 5)
        @test !is_valid(m, pref3)
    end
    # supports (array)
    @testset "supports (array)" begin
        @test_throws ErrorException supports(prefs)
        @test_throws ErrorException supports([prefs[1]; pref])
        set_supports(prefs[1], [0])
        set_supports(prefs[2], [0, 1])
        @test_throws ErrorException supports(prefs)
        add_supports(prefs[1], [1])
        expected = JuMP.Containers.SparseAxisArray[]
        push!(expected, convert(JuMP.Containers.SparseAxisArray, [0; 0]))
        push!(expected, convert(JuMP.Containers.SparseAxisArray, [1; 1]))
        @test supports(prefs) == expected
        expected = JuMP.Containers.SparseAxisArray[]
        push!(expected, convert(JuMP.Containers.SparseAxisArray, [0; 0]))
        push!(expected, convert(JuMP.Containers.SparseAxisArray, [0; 1]))
        push!(expected, convert(JuMP.Containers.SparseAxisArray, [1; 0]))
        push!(expected, convert(JuMP.Containers.SparseAxisArray, [1; 1]))
        set_independent(prefs[1])
        set_independent(prefs[2])
        @test sort(supports(prefs)) == sort(expected)
    end
    # num_parameters
    @testset "num_parameters" begin
        @test num_parameters(m) == 3
        delete!(m.params, JuMP.index(pref))
        @test num_parameters(m) == 2
        m.params[JuMP.index(pref)] = param
    end
    # all_parameters
    @testset "all_parameters" begin
        @test length(all_parameters(m)) == 3
        @test all_parameters(m)[1] == pref
    end
end

# Test methods for finite parameters
@testset "Finite Parameters" begin
    # initialize the model
    m = InfiniteModel()
    # test @finite_parameter
    @testset "@finite_parameter" begin
        m2 = Model()
        # test errors
        @test_macro_throws ErrorException @finite_parameter(m)
        @test_macro_throws ErrorException @finite_parameter(m, a, 2, 3)
        @test_macro_throws ErrorException @finite_parameter(m, (2, 3, 4), 2)
        @test_macro_throws ErrorException @finite_parameter(m2, 2)
        @test_macro_throws ErrorException @finite_parameter(m, "bob")
        @test_macro_throws ErrorException @finite_parameter(m2, test, 2)
        @test_macro_throws ErrorException @finite_parameter(m, test, 2, bob = 2)
        # test anonymous definition
        pref = ParameterRef(m, 1)
        @test @finite_parameter(m, 42) == pref
        @test supports(pref) == [42]
        @test infinite_set(pref) == IntervalSet(42, 42)
        # test vector anonymous definition
        prefs = [ParameterRef(m, 2), ParameterRef(m, 3)]
        @test @finite_parameter(m, [1:2], 42, base_name = "a") == prefs
        @test supports(prefs[1]) == [42]
        @test infinite_set(prefs[2]) == IntervalSet(42, 42)
        @test name.(prefs) == ["a[1]", "a[2]"]
        # test named definition
        pref = ParameterRef(m, 4)
        @test @finite_parameter(m, b, 42) == pref
        @test supports(pref) == [42]
        @test infinite_set(pref) == IntervalSet(42, 42)
        @test name(pref) == "b"
        # test named vector definition
        prefs = [ParameterRef(m, 5), ParameterRef(m, 6)]
        prefs = convert(JuMPC.SparseAxisArray, prefs)
        @test @finite_parameter(m, c[i = 1:2], [3, 7][i],
                                container = SparseAxisArray) == prefs
        @test supports(prefs[2]) == [7]
        @test infinite_set(prefs[1]) == IntervalSet(3, 3)
        @test name(prefs[2]) == "c[2]"
    end
    # test is_finite_parameter
    @testset "is_finite_parameter" begin
        @infinite_parameter(m, 0 <= d <= 1)
        @finite_parameter(m, e, 1)
        @test !is_finite_parameter(d)
        @test is_finite_parameter(e)
    end
    # test JuMP.value
    @testset "JuMP.value" begin
        @infinite_parameter(m, 0 <= f <= 1)
        @finite_parameter(m, g, 1)
        @test value(g) == 1
        @test_throws ErrorException value(f)
    end
    # test JuMP.set_value
    @testset "JuMP.set_value" begin
        @infinite_parameter(m, 0 <= h <= 1)
        @finite_parameter(m, i, 1)
        @test isa(set_value(i, 42), Nothing)
        @test value(i) == 42
        @test_throws ErrorException set_value(h, 42)
    end
end

# Test support flll-in and geneartion functions
@testset "Support Fill-in and Generation" begin
    @testset "generate_and_add_supports! (AbstractInfiniteSet)" begin
        m = InfiniteModel()
        pref1 = @infinite_parameter(m, 0 <= a <= 1)
        set1 = JuMP.owner_model(pref1).params[JuMP.index(pref1)].set
        dist = Normal(0., 1.)
        pref2 = @infinite_parameter(m, c in dist)
        set2 = JuMP.owner_model(pref2).params[JuMP.index(pref2)].set
        @test generate_and_add_supports!(pref1, set1, num_supports = 10) isa Nothing
        @test generate_and_add_supports!(pref2, set2, num_supports = 10) isa Nothing
        @test length(supports(pref1)) == 10
        @test length(supports(pref2)) == 10
    end
    @testset "generate_and_add_supports! (Multivariate DistributionSet)" begin
        m = InfiniteModel()
        dist = MvNormal([0.; 0.], [1. 0.; 0. 2.])
        prefs = @infinite_parameter(m, x[1:2] in dist)
        set = JuMP.owner_model(prefs[1]).params[JuMP.index(prefs[1])].set
        @test generate_and_add_supports!(prefs[1], set, num_supports = 10, sig_fig = 7) isa Nothing
        @test length(supports(prefs[1])) == 10
        @test length(supports(prefs[2])) == 10
    end
    # fill_in_supports! (ParameterRef)
    @testset "fill_in_supports! (ParameterRef)" begin
        m = InfiniteModel()
        pref1 = @infinite_parameter(m, 0 <= a <= 1)
        pref2 = @infinite_parameter(m, 0 <= b[1:2] <= 1)
        dist = Normal(0., 1.)
        pref3 = @infinite_parameter(m, c in dist, supports = [-0.5, 0.5])
        @test fill_in_supports!(pref1, num_supports = 11, sig_fig = 3) isa Nothing
        @test fill_in_supports!.(pref2, num_supports = 11, sig_fig = 3) isa Array{Nothing}
        @test fill_in_supports!(pref3, num_supports = 11, sig_fig = 3) isa Nothing
        @test length(supports(pref1)) == 11
        @test length(supports(pref2[1])) == 11
        @test length(supports(pref2[2])) == 11
        @test supports(pref3) == [-0.5, 0.5]
        @test fill_in_supports!(pref1, num_supports = 20) isa Nothing
        @test length(supports(pref1)) == 11
    end
    # fill_in_supports! (InfiniteModel)
    @testset "fill_in_supports! (InfiniteModel)" begin
        m = InfiniteModel()
        dist1 = Normal(0., 1.)
        dist2 = MvNormal([0.; 0.], [1. 0.;0. 2.])
        pref_a = @infinite_parameter(m, 0 <= a <= 1)
        prefs_b = @infinite_parameter(m, 1 <= b[1:2] <= 2)
        pref_c = @infinite_parameter(m, c in dist1)
        prefs_d = @infinite_parameter(m, d[1:2] in dist2)
        pref_e = @infinite_parameter(m, 2 <= e <= 3, supports = [2.3, 2.7])
        @test fill_in_supports!(m, num_supports = 10) isa Nothing
        @test supports(pref_a) == round.(collect(range(0., stop = 1., step = 1/9)), sigdigits = 5)
        @test supports(prefs_b[1]) == round.(collect(range(1., stop = 2., step = 1/9)), sigdigits = 5)
        @test supports(prefs_b[2]) == round.(collect(range(1., stop = 2., step = 1/9)), sigdigits = 5)
        @test length(supports(pref_c)) == 10
        @test length(supports(prefs_d[1])) == 10
        @test length(supports(prefs_d[2])) == 10
        @test supports(pref_e) == [2.3, 2.7]
    end
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI

pathToFMU = joinpath(dirname(@__FILE__), "../model/SpringPendulum1D.fmu")

myFMU = fmiLoad(pathToFMU)

#create an instance and simulate it
comp1 = fmiInstantiate!(myFMU; loggingOn=true)
fmiSetupExperiment(comp1, 0.0)
fmiEnterInitializationMode(comp1)
fmiExitInitializationMode(comp1)

dt = 0.01
t_start = 0.0
t_stop = 8.0

data1 = fmiSimulateCS(comp1, dt, t_start, t_stop, ["mass.s"])
fmiPlot(data1)

#create another instance, change the spring stiffness and simulate it
comp2 = fmiInstantiate!(myFMU; loggingOn=true)
fmiSetupExperiment(comp2, 0.0)
fmiEnterInitializationMode(comp2)
springConstant = fmiGetReal(comp2, "spring.c") * 0.1
fmiSetReal(comp2, "spring.c", springConstant)
fmiExitInitializationMode(comp2)
data2 = fmiSimulateCS(comp2, dt, t_start, t_stop, ["mass.s"])
fmiPlot(data2)

fmiUnload(myFMU)


text:
using LibClang
using Test

@testset "LibClang.jl" begin
    cxstr = clang_getClangVersion()
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    @test match(r"[0-9]+.[0-9]+.[0-9]+", s).match == string(Base.libllvm_version)
end


text:
# Autogenerated wrapper script for AmazonCACerts_jll for any
export amazon_ca_root_certificate, combined_ca_root_certificates, rds_ca_root_certificate

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `amazon_ca_root_certificate`
const amazon_ca_root_certificate_splitpath = ["share", "AmazonRootCA1.pem"]

# This will be filled out by __init__() for all products, as it must be done at runtime
amazon_ca_root_certificate_path = ""

# amazon_ca_root_certificate-specific global declaration
# This will be filled out by __init__()
amazon_ca_root_certificate = ""

# Relative path to `combined_ca_root_certificates`
const combined_ca_root_certificates_splitpath = ["share", "combined-amazon-root-ca.pem"]

# This will be filled out by __init__() for all products, as it must be done at runtime
combined_ca_root_certificates_path = ""

# combined_ca_root_certificates-specific global declaration
# This will be filled out by __init__()
combined_ca_root_certificates = ""

# Relative path to `rds_ca_root_certificate`
const rds_ca_root_certificate_splitpath = ["share", "rds-ca-2019-root.pem"]

# This will be filled out by __init__() for all products, as it must be done at runtime
rds_ca_root_certificate_path = ""

# rds_ca_root_certificate-specific global declaration
# This will be filled out by __init__()
rds_ca_root_certificate = ""

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"AmazonCACerts")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    global amazon_ca_root_certificate_path = normpath(joinpath(artifact_dir, amazon_ca_root_certificate_splitpath...))

    global amazon_ca_root_certificate = amazon_ca_root_certificate_path
    global combined_ca_root_certificates_path = normpath(joinpath(artifact_dir, combined_ca_root_certificates_splitpath...))

    global combined_ca_root_certificates = combined_ca_root_certificates_path
    global rds_ca_root_certificate_path = normpath(joinpath(artifact_dir, rds_ca_root_certificate_splitpath...))

    global rds_ca_root_certificate = rds_ca_root_certificate_path
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
function nutrients_init(arch, g)
    fields = (Field(arch, g), Field(arch, g),
              Field(arch, g), Field(arch, g),
              Field(arch, g), Field(arch, g),
              Field(arch, g), Field(arch, g),
              Field(arch, g), Field(arch, g))

    nut = NamedTuple{nut_names}(fields)
    return nut
end

"""
    default_nut_init()
Generate defalut nutrient initial conditions.
"""
function default_nut_init()
    init = (DIC=20.0, NH4=0.2, NO3=0.5, PO4=0.03, DOC=1.0, DON=0.1, DOP=0.05, POC=0.0, PON=0.0, POP=0.0)
    rand_noise = (DIC=0.1, NH4=0.1, NO3=0.1, PO4=0.1, DOC=0.1, DON=0.1, DOP=0.1, POC=0.1, PON=0.1, POP=0.1)
    return (initial_condition = init, rand_noise = rand_noise)
end

"""
    generate_nutrients(arch, grid, source; mask = nothing)
Set up initial nutrient fields according to `grid`.

Keyword Arguments
=================
- `arch`: `CPU()` or `GPU()`. The computer architecture used to time-step `model`.
- `grid`: The resolution and discrete geometry on which nutrient fields are solved.
- `source`: A `NamedTuple` containing 10 numbers each of which is the uniform initial condition of one tracer,
            or a `Dict` containing the file paths pointing to the files of nutrient initial conditions.
- `mask` (optional): Mask out the tracers generated out of the domain, a 3D array with size `(Nx, Ny, Nz)`.
"""
function generate_nutrients(arch, g, source::Union{Dict,NamedTuple}; mask=nothing)
    total_size = (g.Nx+g.Hx*2, g.Ny+g.Hy*2, g.Nz+g.Hz*2)
    nut = nutrients_init(arch, g)
    pathkeys = collect(keys(source))

    if typeof(source) <: NamedTuple
        pathkeys = collect(keys(source.initial_condition))
    end

    for name in nut_names
        if length(findall(x->x==name, pathkeys)) == 0
            throw(ArgumentError("NUT_INIT: nutrient not found $(name)"))
        else
            if typeof(source) <: Dict # file paths
                tmp = deserialize(source[name]) |> array_type(arch)
                if sum(tmp .< 0.0) > 0
                    throw(ArgumentError("NUT_INIT:  The initial condition should be none-negetive."))
                end

                if size(tmp) == (g.Nx, g.Ny, g.Nz)
                    @views @. nut[name].data[g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz] = tmp[:,:,:]
                else
                    throw(ArgumentError("NUT_INIT:  grid mismatch"))
                end
            elseif typeof(source) <: NamedTuple # value
                if source.initial_condition[name] < 0.0
                    throw(ArgumentError("NUT_INIT:  The initial condition should be none-negetive."))
                end
                lower = 1.0 - source.rand_noise[name]
                upper = 1.0 + source.rand_noise[name]
                nut[name].data .= fill(source.initial_condition[name],total_size) .* rand(lower:1e-4:upper, total_size) |> array_type(arch)
            end
        end

        if mask ≠ nothing
            if size(mask) == (g.Nx, g.Ny, g.Nz)
                @views @. nut[name].data[g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz] *= mask
            else
                throw(ArgumentError("nut_mask: grid mismatch, size(mask) must equal to (grid.Nx, grid.Ny, grid.Nz)."))
            end
        end
    end

    fill_halo_nut!(nut,g)

    return nut
end


text:
#__precompile__(true) #needed?

module ForestBiometrics

#dependencies
using DataFrames
using Plots
using StatPlots

#end dependencies

export emc
export sdi, qmd
export gingrich_chart
export HeightDiameter, calculate_height

export Curtis,
Michailoff,
Meyer,
Micment,
Micment2,
Naslund,
Naslund2,
Naslund3,
Naslund4,
Power,
Wyckoff

#3 parameter equations, mainly from LMFOR R package
#sorted alphabetically
export Chapman,
Gompertz,
HossfeldIV,
Korf,
Logistic,
Monserud,
Prodan,
Ratkowsky,
Sibbesen,
Weibull

export K, KMETRIC,VolumeEquation, MerchSpecs, Sawtimber, Pulp, Fiber

export limiting_distance

export sdi_chart

export  Log, LogSegment,
        Shape, Cylinder,
        Paraboloid,
        Cone,Neiloid,
        ParaboloidFrustrum,
        NeiloidFrustrum,
        ConeFrustrum,
        area,volume,
        scribner_volume,
        international_volume,
        doyle_volume

#Alphabetic order

include("EquilibriumMoistureContent.jl")
include("ForestStocking.jl")
include("GingrichStocking_chart.jl")
include("HeightDub.jl")
include("LimitingDistance.jl")
include("SDI_chart.jl")
include("VolumeEquations.jl")

end
#end module


text:
######################################################################################################
# Copy from https://github.com/xuedong/LinBAI.jl  by Xuedong Shang                                   #
######################################################################################################
using JLD2;
using Printf;
using StatsPlots;
using LaTeXStrings;
include("runit_bai.jl"); # for types
include("../utilities/experiment_helpers.jl");
include("../utilities/thresholds.jl");

name = ARGS[1];

@load "$name.dat" dist μ pep srs data δs βs repeats seed

dump_stats(pep, μ, δs, βs, srs, data, repeats);

for i in 1:length(δs)
    plot(_boxes(pep, μ, δs[i], βs[i], srs, getindex.(data, i), repeats));
    savefig("$(name)_$(δs[i])_no_outlier.pdf");
end


text:
## HELPER

"""to display abbreviated versions of integers"""
abbreviated(n) = "..."*string(n)[end-2:end]

## EXPOSED SHOW METHODS

# short version of showing a `BaseType` object:
function Base.show(stream::IO, object::BaseType)
    str = string(typeof(object).name.name,
                        "@", abbreviated(objectid(object)))
    if !isempty(fieldnames(typeof(object)))
        printstyled(IOContext(stream, :color=> true), str,
                    color=:blue)
    else
        print(stream, str)
    end
end

# longer version of showing a `BaseType` object (used at REPL):
function Base.show(stream::IO, ::MIME"text/plain", object::BaseType)
    _recursive_show(stream, object, 1, DEFAULT_SHOW_DEPTH)
end

# version showing a `BaseType` object to arbitrary depth:
function Base.show(stream::IO, object::BaseType, depth::Int)
    _recursive_show(stream, object, 1, depth)
end

# for convenience:
Base.show(object::BaseType, depth::Int) = show(stdout, object, depth)
macro more()
    esc(quote
        show(Main.ans, 100)
    end)
end

## METHODS TO SUPRESS THE DISPLAY OF LARGE NON-BASETYPE OBJECTS

istoobig(::Any) = true
istoobig(::Number) = false
istoobig(x::AbstractArray) = maximum(size(x)) > 5
istoobig(d::Dict{T,Any}) where T <: Union{Number,Symbol} = length(keys(d)) > 5

## THE `_show` METHOD

# Note: The `_show` method controls how field values are displayed in
# the table generated by `_recursive_show`. See top of file.

# _show fallback:
function _show(stream::IO, object)
    if !istoobig(object)
        show(stream, MIME("text/plain"), object)
        println(stream)
    else
        println(stream, "omitted ", typeof(object))
    end
end

# _show for BaseType:
_show(stream::IO, object::BaseType) = println(stream, object)

# _show for other types:
function _show(stream::IO, df::AbstractDataFrame)
    println(stream, "omitted $(typeof(df)) of size $(size(df))")
end

function _show(stream::IO, A::Array{T}) where T
    if !istoobig(A)
        show(stream, MIME("text/plain"), A)
        println(stream)
    else
        println(stream, "omitted Array{$T} of size $(size(A))")
    end
end

function _show(stream::IO, d::Dict{T, Any}) where T <: Union{Number,Symbol}
    if isempty(d)
        println(stream, "empty $(typeof(d))")
    elseif !istoobig(d)
        println(stream, "omitted $(typeof(d)) with keys: ")
        show(stream, MIME("text/plain"), collect(keys(d)))
        println(stream)
    else
        println(stream, "omitted $(typeof(d))")
    end
end

function _show(stream::IO, v::Array{T, 1}) where T
    if !istoobig(v)
        show(stream, MIME("text/plain"), v)
        println(stream)
    else
        println(stream, "omitted Vector{$T} of length $(length(v))")
    end
end

_show(stream::IO, ::Nothing) = println(stream, "nothing")

"""
    _recursive_show(stream, object, current_depth, depth)

Generate a table of the field values of the `BaseType` object,
dislaying each value by calling the method `_show` on it. The
behaviour of `_show(stream, f)` is as follows:

1. If `f` is itself a `BaseType` object, then its short form is shown
and `_recursive_show` generates as separate table for each of its
field values (and so on, up to a depth of argument `depth`).

2. Otherwise `f` is displayed as "omitted T" where `T = typeof(f)`,
unless `istoobig(f)` is false (the `istoobig` fall-back for arbitrary
types being `true`). In the latter case, the long (ie,
MIME"plain/text") form of `f` is shown. To override this behaviour,
overload the `_show` method for the type in question.

"""
function _recursive_show(stream::IO, object::BaseType, current_depth, depth)
    if depth == 0 || isempty(fieldnames(typeof(object)))
        println(stream, object)
    elseif current_depth <= depth
        names = fieldnames(typeof(object))
        print(stream, "#"^current_depth, " ")
        show(stream, object)
        println(stream, ": ")
        if isempty(names)
            println(stream)
            return
        end
        for fld in names
            fld_string = string(fld)*" "^(max(0,COLUMN_WIDTH - length(string(fld))))*"=>   "
            print(stream, fld_string)
            if isdefined(object, fld)
                _show(stream, getfield(object, fld))
                #           println(stream)
            else
                println(stream, "(undefined)")
                #           println(stream)
            end
        end
        println(stream)
        for fld in names
            if isdefined(object, fld)
                subobject = getfield(object, fld)
                if isa(subobject, BaseType) && !isempty(fieldnames(typeof(subobject)))
                    _recursive_show(stream, getfield(object, fld),
                                    current_depth + 1, depth)
                end
            end
        end
    end
end


text:
using TypeClasses

# Monoid instances
# ================

# this is standard Applicative combine implementation
# there is no other sensible definition for combine and hence if this does fail, it fails correctly
TypeClasses.combine(a::State, b::State) = mapn(combine, a, b)

# we don't do the same for orelse, as orelse lives on the container level, but there is  no default definition of orelse for State

# FunctorApplicativeMonad
# =======================

# there is no implementation of foreach for State, as we cannot look into the state function

TypeClasses.map(f, a::State) = State() do state
  value, newstate = a(state)
  f(value), newstate
end

TypeClasses.pure(::Type{<:State}, a) = State() do state
  a, state
end

TypeClasses.ap(f::State, a::State) = State() do state0
  func, state1 = f(state0)
  value, state2 = a(state1)
  func(value), state2
end

# we cannot overload this generically, because `Base.map(f, ::Vector...)` would get overwritten as well (even without warning surprisingly)
TypeClasses.map(f, a::State, b::State, more::State...) = mapn(f, a, b, more...)

TypeClasses.flatmap(f, a::State) = State() do state0
  value, state1 = a(state0)
  convert(State, f(value))(state1)
end

# FlipTypes
# =========

# flip_types does not makes sense for State, as we cannot evalute the State without running the statefunctions on top
# of a given initial state


text:
## Verify that the data is self-consistent and consistent with the likelihood ##
function check_data!(X::AbstractArray{T₁,N₁},y::AbstractArray{T₂,N₂},likelihood::Union{Distribution,Likelihood}) where {T₁<:Real,T₂,N₁,N₂}
    @assert (size(y,1)==size(X,1)) "There is not the same number of samples in X and y";
    @assert N₁ <= 2 "The input matrix X can only be a vector or a matrix"
    y,nLatent,likelihood = treat_labels!(y,likelihood)
    if N₁ == 1 #Reshape a X vector as a matrix
        return reshape(X,:,1),y,nLatent,likelihood
    else
        return X,y,nLatent,likelihood
    end
end

## Verify that the likelihood and inference are compatible (are implemented) ##
function check_implementation(model::Symbol,likelihood::L,inference::I) where {I<:Inference,L<:Likelihood}
    if isa(likelihood,GaussianLikelihood)
        if model == :GP && inference isa Analytic
            return true
        elseif (model == :SVGP || model == :VStP) && inference isa AnalyticVI
            return true
        else
            return false
        end
    elseif likelihood isa StudentTLikelihood
        if inference isa AnalyticVI || inference isa QuadratureVI
            return true
        elseif model == :VGP && inference isa GibbsSampling
            return true
        else
            return false
        end
    elseif likelihood isa LaplaceLikelihood
        if inference isa AnalyticVI || inference isa QuadratureVI
            return true
        elseif model == :MCGP && inference isa GibbsSampling
            return true
        else
            return false
        end
    elseif likelihood isa HeteroscedasticLikelihood
        if inference isa AnalyticVI
            return true
        else
            return false
        end
    elseif likelihood isa HeteroscedasticLikelihood
        if inference isa AnalyticVI || inference isa QuadratureVI
            return true
        elseif model == :VGP && inference isa GibbsSampling
            return true
        else
            return false
        end
    elseif likelihood isa LogisticLikelihood
        if inference isa AnalyticVI || inference isa QuadratureVI
            return true
        elseif model == :MCGP && inference isa GibbsSampling
            return true
        else
            return false
        end
    elseif likelihood isa BayesianSVM
        if inference isa AnalyticVI
            return true
        else
            return false
        end
    elseif likelihood isa SoftMaxLikelihood
        if inference isa MCIntegrationVI
            return true
        else
            return false
        end
    elseif likelihood isa LogisticSoftMaxLikelihood
        if inference isa AnalyticVI
            return true
        elseif model == :MCGP && inference isa GibbsSampling
            return true
        elseif inference isa MCIntegrationVI
            return true
        else
            return false
        end
    elseif likelihood isa PoissonLikelihood
        if inference isa AnalyticVI
            return true
        else
            return false
        end
    elseif likelihood isa NegBinomialLikelihood
        if inference isa AnalyticVI
            return true
        else
            return false
        end
    elseif likelihood isa CustomLikelihood
        if inference isa NumericalVI || inference isa HMCSampling
            return true
        else
            return false
        end
    else
        return true
    end
end


text:
include("onda/preamp/AH2020_1238_0dB.jl")
include("onda/preamp/AH2020_1238_20dB.jl")
include("onda/hydrophone/HGL0200_2322.jl")
include("onda/hydrophone/HGL0200_2512.jl")
include("precision_acoustics/PA_3197.jl")


text:
module MRIgeneralizedBloch

using QuadGK
using DifferentialEquations
using Interpolations
using ApproxFun
import Cubature
using SpecialFunctions
using StaticArrays
using LinearAlgebra
using NLsolve
using ExponentialUtilities
using Polyester

export apply_hamiltonian_gbloch!
export apply_hamiltonian_graham_superlorentzian!
export apply_hamiltonian_sled!
export hamiltonian_linear
export d_hamiltonian_linear_dω1

export greens_lorentzian
export greens_gaussian
export greens_superlorentzian
export dG_o_dT2s_x_T2s_lorentzian
export dG_o_dT2s_x_T2s_gaussian
export dG_o_dT2s_x_T2s_superlorentzian
export interpolate_greens_function

export calculatesignal_gbloch_ide
export calculatesignal_graham_ode

export precompute_R2sl
export evaluate_R2sl_vector
export calculatesignal_linearapprox

export grad_m0s
export grad_R1a
export grad_R1f
export grad_R1s
export grad_R2f
export grad_Rx
export grad_T2s
export grad_ω0
export grad_B1

abstract type grad_param end
struct grad_m0s <: grad_param end
struct grad_R1a <: grad_param end
struct grad_R1f <: grad_param end
struct grad_R1s <: grad_param end
struct grad_R2f <: grad_param end
struct grad_Rx  <: grad_param end
struct grad_T2s <: grad_param end
struct grad_ω0  <: grad_param end
struct grad_B1  <: grad_param end

include("DiffEq_Hamiltonians.jl")
include("Linearized_R2s.jl")
include("MatrixExp_Solvers.jl")
include("DiffEq_Sovlers.jl")
include("Greens_functions.jl")
include("MatrixExp_Hamiltonians.jl")
include("OptimalControl.jl")
include("Gradient_Hamiltonians.jl")

end


text:
using Distributed

"""
    read_trajectory(fname::AbstractString) -> Trajectory

Reads trajectory from given filename.
Type of file determined based on sufix.

## List of suported file formats
- XYZ
- PDB
"""
function read_trajectory(fname::AbstractString)
    sufix=split(fname, ".")[end]
    if sufix in ("pdb", "PDB")
        return read_pdb(fname)
    elseif sufix in ("xyz", "XYZ")
        return read_xyz(fname)
    elseif sufix in ("dcd", "DCD")
        return read_dcd(fname)
    else
        error("Trajectory file format not recognized")
    end
end

"""
    read_trajectory(fnames::AbstractString...) -> Vector{Trajectory}

Reads given files using threading
"""
function read_trajectory(fnames::AbstractString...)
    return fetch.( [Threads.@spawn(read_trajectory(x)) for x in fnames] )
end

"""
    read_dcd(fname::AbstractString) -> Trajectory

Reads DCD trajectory file format file.

Format definition is based on CP2K.
[Reference](https://github.com/cp2k/cp2k/blob/master/src/motion/dumpdcd.F)
"""
function read_dcd(fname::AbstractString)
    function read_cell(io::IO)
        a = read(io, Float64)
        γ = read(io, Float64)
        b = read(io, Float64)
        β = read(io, Float64)
        α = read(io, Float64)
        c = read(io, Float64)
        return parsecell(a,b,c,α,β,γ)
    end
    x = Float32[]
    y = Float32[]
    z = Float32[]
    natom_dcd = zero(Int32)
    cell = Vector{AbstractUnitCell}()
    frame = 0
    open(fname,"r") do io
        read(io,4) # Skip empty space

        id_dcd = String(read(io,4))
        istep = read(io, Int)
        iskip = read(io, Int32)

        read(io,24) # Skip empty space

        dt = read(io, Float32)
        have_unit_cell = read(io, Int32)

        read(io, 44) # Skip empty space

        nremark = read(io, Int32)
        remark_dcd = String(read(io, 80))
        remark_dcd2 = String(read(io, 80))

        read(io, 8) # Skip empty space

        natom_dcd = read(io, Int32)

        @debug "1st phase of DCD-file" natom_dcd, dt, id_dcd, istep, iskip, have_unit_cell, nremark

        while true
            try
                read(io, 8)
                push!(cell, read_cell(io))
                read(io, 8)
                append!(x , [read(io, Float32) for _ in 1:natom_dcd])
                read(io, 8)
                append!(y, [read(io, Float32) for _ in 1:natom_dcd])
                read(io, 8)
                append!(z, [read(io, Float32) for _ in 1:natom_dcd])
                frame += 1
            catch e
                if isa(e, EOFError)
                    break
                else
                    throw(e)
                end
            end
        end
    end
    @debug "Frame = $frame, natoms = $natom_dcd, length = $(length(x))"
    @assert length(x) == length(y) == length(z) "x-, y- or z coordinate lenghts dont match"
    @debug "array" Array(hcat(x,y,z)')
    return Trajectory( reshape(Array(hcat(x,y,z)'), 3, natom_dcd, frame), cell)
end

read_xyz(fname::AbstractString) = open(read_xyz, fname, "r")

function read_xyz(io::IO)
    lines = Vector{String}()
    atoms = Vector{String}()
    xyz = Vector{Float64}()

    na = parse(Int,readline(io))  # Number of atoms
    line = readline(io)   # Comment line
    for i in 1:na
        line = readline(io)
        cont = split(line)
        push!(atoms,cont[1])
        append!(xyz, parse.(Float64, cont[2:4]))
    end
    i = 1
    for line in eachline(io)
        cont = split(line)
        if i == na+2
            i = 1
            continue
        elseif  i > 2
            if length(cont) >= 4
                append!(xyz, parse.(Float64, cont[2:4]))
            else
                error("File is not in xyz format")
            end
        end
        i += 1
    end
    @debug "help " na, length(xyz), length(xyz)/(3*na)
    return Trajectory( reshape(xyz, 3, na, Int(length(xyz)/(3*na))), atoms)
end

function write_xyz(io::IO, traj::Trajectory{AtomNames, TC}) where {TC<:AbstractUnitCell}
    for x in traj
        println(io, "  ", length(traj.names),"\n")
        for (n, r) in zip(traj.names, eachcol(x))
            println(io, n, "   ", r[1], "  ", r[2], "  ", r[3])
        end
    end
    return nothing
end

function write_xyz(fname::AbstractString, traj::Trajectory)
    open(fname, "w") do file
        write_xyz(file, traj)
    end
    return nothing
end

function read_pdb(fname::AbstractString)
    xyz = Vector{Float64}()
    atoms = Vector{String}()
    crystal = Vector{AbstractUnitCell}()
    open(fname, "r") do file
        lineiterator = eachline(file)

        # Read data that need to read only once, like atom names
        for line in lineiterator
            if occursin(r"^CRYST1", line)
                terms = split(line)
                push!(crystal, parsecell(parse.(Float64, terms[2:7])...))
            elseif occursin(r"^ATOM", line) || occursin(r"^HETATM", line)
                push!(atoms, line[77:78])
                push!(xyz, parse(Float64, line[31:38]))
                push!(xyz, parse(Float64, line[39:46]))
                push!(xyz, parse(Float64, line[47:54]))
            elseif occursin(r"^END", line)
                break
            end
        end

        # Read bulk data
        for line in lineiterator
            if occursin(r"^ATOM", line) || occursin(r"^HETATM", line)
                push!(xyz, parse(Float64, line[31:38]))
                push!(xyz, parse(Float64, line[39:46]))
                push!(xyz, parse(Float64, line[47:54]))
            elseif occursin(r"^CRYST1", line)
                terms = split(line)
                push!(crystal, parsecell(parse.(Float64, terms[2:7])...))
            end
        end
    end

    @assert length(atoms)*3*length(crystal) == length(xyz) "Number of atoms and parsed coordinates do not match."

    return Trajectory(reshape(xyz ,3, length(atoms), length(crystal)), atoms, crystal)
end


text:
"""
OptimizeMe is a demonstration module used in illustrating how to improve code and generate effective `precompile` directives.
It has deliberate weaknesses in its design, and the analysis of these weaknesses via `@snoopi_deep` is discussed
in the documentation.
"""
module OptimizeMeFixed

struct Container{T}
    value::T
end

function lotsa_containers()
    list = Any[1, 0x01, 0xffff, 2.0f0, 'a', [0], ("key", 42)]
    cs = Container{Any}.(list)
    println("lotsa containers:")
    display(cs)
end

concat_string(c1::Container, c2::Container) = string(c1.value) * ' ' * string(c2.value)

function contain_concrete(item1, item2)
    c1 = Container(item1)
    c2 = Container(item2)
    return concat_string(c1, c2)
end

function contain_list(list)
    length(list) == 2 || throw(DimensionMismatch("list must have length 2"))
    item1 = convert(Float64, list[1])::Float64
    item2 = list[2]::String
    return contain_concrete(item1, item2)
end

struct Object
    x::Int
end
Base.show(io::IO, o::Object) = print(io, "Object x: ", o.x)

function makeobjects()
    xs = [1:5; 7:7]
    return Object.(xs)
end

# "Stub" callers for precompilability
function warmup()
    mime = MIME("text/plain")
    io = Base.stdout::Base.TTY
    v = [Container{Any}(0)]
    show(io, mime, v)
    show(IOContext(io), mime, v)
    v = [Object(0)]
    show(io, mime, v)
    show(IOContext(io), mime, v)
    return nothing
end

function main()
    lotsa_containers()
    println(contain_concrete(3.14, "is great"))
    list = [2.718, "is jealous"]
    println(contain_list(list))
    display(makeobjects())
end

if Base.VERSION >= v"1.4.2"
    precompile(Tuple{typeof(main)})   # time: 0.4204474
    precompile(Tuple{typeof(warmup)})
end

end


text:
"""
Helper functions

RT_to_scans: Converts retention time to indexes in spectrum
filter_XIC: Returns spectrum filtered at chosen mz values
batch_import: Imports spectra present in given folder
retrieve_sample_name: Retrieves sample names for the spectra from the agilent .D folder

"""

"""Returns index range of given retention time range"""
function RT_to_scans(spectrum, RT_range)
	# Only possible if Rt is sorted
	if !issorted(spectrum["Rt"])
		error("Error: Retention time array is not in order")
	end

	# Convert Integer to tuple
	if RT_range isa(Number)
		RT_range = (RT_range, RT_range)
	end

	return searchsortedfirst(spectrum["Rt"], RT_range[1]), searchsortedlast(spectrum["Rt"], RT_range[2])
end

"""Returns filtered spectrum based on given mass (mz) values"""
function filter_XIC(spectrum, mass_values)
	# Return complete spectrum
	if mass_values == 0
		return reduce(+, spectrum["Mz_intensity"], dims=2)
	# Convert to tuple
	elseif typeof(mass_values) == Int
		mass_values = (mass_values,)
	end

	indices = Array{CartesianIndex{2}, 1}()

	# Add indices that correspond to given mz values
	for mass in mass_values
		mass_range = [mass - MAX_MZ_DEVIATION, mass + MAX_MZ_DEVIATION]

		# Store which indices are between mass range
		append!(indices, findall(mz -> (mz .>= mass_range[1]) .& (mz .<= mass_range[2]), spectrum["Mz_values"]))
	end

	# Create XIC spectrum
	spectrum_XIC = zeros(Int, size(spectrum["Mz_values"], 1))

	# Filter spectrum using stored indices
	for I in indices
		spectrum_XIC[I[1]] += spectrum["Mz_intensity"][I]
	end

	return spectrum_XIC
end

"""
	parse_data(row)
Parses mz values and RT overlap values from string to array
return mz_values, RT_overlap_values
"""
function parse_data(row)
	# Parse mz values
	mz_values_string = split(string(row.mz), ";")
	mz_values = Array{Any,1}(undef, length(mz_values_string))
	for (i, mz) in enumerate(mz_values_string)

		# Tuple of mz values
		if startswith(mz, "(")
			mz_tuple_string = split(mz[2:end - 1], ",")
			mz_values[i] = Tuple(([parse(Float32, sub_mass) for sub_mass in mz_tuple_string]))

		# Separate mz value
		else
			mz_values[i] = parse(Float32, mz)
		end
	end

	# Parse overlap values
	RT_overlap_values_string = split(string(row.overlap), ",")
	RT_overlap_values = zeros(Float16, 2)
	for RT_overlap in RT_overlap_values_string
		RT_overlap = RT_overlap == "missing" ? 0 : parse(Float16, RT_overlap)
		index = RT_overlap < row.RT ? 1 : 2
		RT_overlap_values[index] = RT_overlap
	end

	return mz_values, RT_overlap_values

end

"""
	batch_import(pathin, settings_json; use_filters=true)
Returns list containing all imported spectra from folder and metadata headers
Set use_filters to false to ignore filters and import all spectra,
"""
function batch_import(pathin, settings_json; use_filters=true)

	files = readdir(pathin, sort=false)
	files_supported = zeros(Bool, length(files),1)

	# Determine supported files
	for i=1:length(files)
		if length(files[i]) > length(".mzxml") && lowercase(files[i][end-4:end]) == "mzxml"
			files_supported[i] = true
		end
	end

	# Save indices of supported files
	supported_indices=findall(x -> x == true, files_supported)
	num_of_spectra = length(supported_indices)

	# Retrieve which metadata to save
	info_containing = settings_json[Symbol("output_metadata")]
	sample_info_metadata = Dict(i["element_name"] => i["header"] for i in info_containing if i["type"] == "sample_info_xml")
	metadata_headers = append!([], values(sample_info_metadata), ["Filename", "Folder Name"])
	# Retrieve filters
	filters = settings_json[Symbol("filters")]
	filters_must_contain = Dict(i["element_name"] => i["value"] for i in filters if i["filter"] == "must_contain")
	empty_field_allowed = Dict(i["element_name"] => i["empty_field_allowed"] for i in filters)

	spectra = []
	# Load files into spectra array
	for (i, file_index) in enumerate(supported_indices)
		filename = files[file_index]

		# Retrieve sample info and check filter requirements
		element_names = append!([], keys(sample_info_metadata), keys(filters_must_contain))
		sample_info = retrieve_sample_info(pathin, filename, element_names)
		filters_accepted = true
		for (element_name, value) in filters_must_contain
			if !(contains(sample_info[element_name], value) || (empty_field_allowed[element_name] && isempty(sample_info[element_name])))
				filters_accepted = false
				break
			end
		end

		if filters_accepted == false && use_filters == true
			@info "Spectrum $(i[1]) of $num_of_spectra did not meet filter requirements"
			continue
		end

		# Read spectrum
		spectrum = import_files_light(pathin,[filename])
		push!(spectra, spectrum)

		# Add metadata
		spectrum["MS1"]["Filename"] = filename
		merge!(spectrum["MS1"], Dict(header => sample_info[element_name] for (element_name, header) in sample_info_metadata))
		spectrum["MS1"]["Folder Name"] = split(pathin, "\\")[end]

		@info "Read spectra $(i[1]) of $num_of_spectra"
	end

	return spectra, metadata_headers
end

function retrieve_sample_info(pathin, filename, element_names)
	# Look for file in same folder (created by .bat file)
	filename_without_ext = filename[begin:end-6]
	sample_info_file = joinpath(pathin, filename_without_ext*"_sample_info.xml")

	# If it does not exist, look for file in .D folder
	if !isfile(sample_info_file)
		folder_loc = joinpath(pathin, filename_without_ext*".D")
		sample_info_file = joinpath(folder_loc, "AcqData\\sample_info.xml")
	end

	# File does not exist at all
	if !isfile(sample_info_file)
		@warn "Sample info file for $filename_without_ext does not exist, metadata will not be saved"
		return Dict(element_name => "" for element_name in element_names)
	end

	# Load xml file into array
	xdoc = parse_file(sample_info_file)
	xroot = root(xdoc)
	xarray = xroot["Field"]

	# Find index where element is stored and retrieve its value
	elements = Dict{String, String}()
	for element_name in element_names
		element_index = findfirst(i -> content(find_element(i, "Name")) == element_name, xarray)
		element_value = content(find_element(xarray[element_index], "Value"))
		elements[element_name] = element_value
	end

	# Free memory
	free(xdoc)

	return elements
end


text:
# JuMP example for the Julia Working Group

# Reproduce the example HS071.
# HS071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
# Start at (1,5,5,1)
# End at (1.000..., 4.743..., 3.821..., 1.379...)

# First install the package, and a solver
#Pkg.add("JuMP")
#Pkg.add("Ipopt") #solver

# Load these packages
using JuMP
using Ipopt
#using NLopt

tic()  #start clock

#Create the model object
m = Model()
#m = Model(solver = NLoptSolver(algorithm=:LD_SLSQP))
#m=Model(solver=IpoptSolver(tol=1e-8, max_iter=5000)) # can use the argument solver to pass solver options

#Defining Variables - These are defined using the @defVar macro. The first argument of
# @defVar(m,x) will always be the model to associate the variable with. You can also
# impose the bounds on the variable here; and start values (using the start option)
@defVar(m, 1<=x[1:4]<=5)
# Set start value
setValue(x[1],1)
setValue(x[2],5)
setValue(x[3],5)
setValue(x[4],1)

# One can also create variables manually one by one as:
# x = Variable(m::Model, lower::Number, upper::Number, category::Symbol, name::String)
# Bounds can also be set using setLower(x::Variable, lower) and setUpper(x::Variable, lower)

# Objective function - I will show the syntax for a nonlinear objective function. For the
# syntax for linear objective functions, see the documentation. Essentially one needs to
# use setObjective() or @setObjective() rather than @setNLObjective()
# For nonlinear objectives, these are defined using the macro @setNLObjective()

@setNLObjective(m::Model, :Min, x[1]*x[4]*(x[1]+x[2]+x[3])+x[3] )

# Next, set constraints.
# Linear constraints are added through the macro @addConstraint()
# Nonlinear constraints are added through the @addNLConstraint() macro
@addNLConstraint(m::Model, x[1]*x[2]*x[3]*x[4]>=25 )
@addNLConstraint(m::Model, x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 == 40)

# Finally, solve the model
solve(m)

# Let's look at the solutions
println("x=", getValue(x))

toc()

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
module TestDoctest

using Documenter
using Test
using Maybe

@testset "/docs" begin
    doctest(Maybe; manual = true)
end

end  # module


text:
using CSV, DataFrames, Printf
using Random

function calc_pi(n::Integer)
   x = rand(n)
   y = rand(n)
   num_in_unit_circle = sum(x.*x.+y.*y .< 1)
   num_in_unit_circle * 4/n
end

job_array_data = CSV.read("ex1_job_array_in.csv", DataFrame)
job_arrayid = parse(Int64,ARGS[1])

idx = findfirst(x-> x==job_arrayid, job_array_data[!,:array_id] )
@assert idx != nothing
@assert 1 <= idx  <= size(job_array_data, 1)

n = job_array_data[idx,:n]
s = job_array_data[idx,:seed]

Random.seed!(s)

calc_pi(1)
@time pi_estimate = calc_pi(n)
println("# After ", n, " itterations, array_id= ", job_arrayid, " estimated pi to be...")
println(pi_estimate)

output_filename = @sprintf("ex1_out_%02d.csv",job_arrayid)
df_out = DataFrame(result = [pi_estimate])
CSV.write(output_filename, df_out)


text:
module nGraph

# stdlib
using Random

# external dependcies
import Cassette
#using Distributions
import Flux
import NNlib
import ProgressMeter

using Dates

# Turn on CPU code generation by default.
function __init__()
    enable_codegen()
end

const SRCDIR = @__DIR__
const PKGDIR = dirname(SRCDIR)
const DEPSDIR = joinpath(PKGDIR, "deps")
const MODELDIR = joinpath(PKGDIR, "models")

# For convenient overloading in ops.jl
import Base: broadcasted

include("env.jl")
include("lib.jl")
import .Lib

ngraph_convert(x) = x

# Expected Transformation
#
# f(a,b,c) -> Lib.f(ngraph_convert(a), ngraph_convert(b), ngraph_convert(c))
#
# Pretty straightforward.
#
# If it becomes relevant, I've left code in that will do the transformation
#
# f(a,b,c) ->
#    #genysym1 = ngraph_convert(a)
#    #genysym2 = ngraph_convert(b)
#    #genysym2 = ngraph_convert(c)
#    GC.@preserve #gensym1 #gensym2 @gensym3 begin
#       Lib.f(#gensym1, #gensym2, @gensym3)
#    end
macro ngraphcall(expr)
    if expr.head != :call
        error("Only call `@ngraphcall` on function calls")
    end

    # Prefix "Lib." in front of the function call
    fname = expr.args[1]
    if !isa(fname, Symbol)
        error("Don't put the modlue `Lib` in front of the function call")
    end
    fname = :(Lib.$fname)

    args = expr.args[2:end]
    for i in eachindex(args)
        args[i] = :(ngraph_convert($(esc(args[i]))))
    end

    return :($fname($(args...)))

    # vars = [gensym() for i in 1:length(args)]

    # # Wrap the arguments in ngraph_convert
    # newargs = []
    # for (var, arg) in zip(vars, args)
    #     newarg = :($(esc(var)) = ngraph_convert($(esc(arg))))
    #     push!(newargs, newarg)
    # end

    # newcall = :($fname($(vars...)))

    # return quote
    #     $(newargs...)
    #     GC.@preserve $(vars...) $newcall
    # end
end

include("types.jl")
include("ops.jl")
include("tracer/tracer.jl")

end # module


text:
using GeometricIntegrators.Integrators
using GeometricIntegrators.Tableaus
using GeometricIntegrators.TestProblems.Oscillator
using GeometricIntegrators.Utils
using Test

using GeometricIntegrators.TestProblems.Oscillator: Δt, nt, refx

sode = oscillator_sode()

sint = Integrator(sode, getTableauLieA(), Δt)
ssol = integrate(sint, nt)

@test rel_err(ssol.q, refx) < 5E-2

sint = Integrator(sode, getTableauLieB(), Δt)
ssol = integrate(sint, nt)

@test rel_err(ssol.q, refx) < 5E-2

sint = Integrator(sode, getTableauStrang(), Δt)
ssol = integrate(sint, nt)

@test rel_err(ssol.q, refx) < 1E-3

sint = Integrator(sode, getTableauMcLachlan2(), Δt)
ssol = integrate(sint, nt)

@test rel_err(ssol.q, refx) < 1E-4

sint = Integrator(sode, getTableauMcLachlan4(), Δt)
ssol = integrate(sint, nt)

@test rel_err(ssol.q, refx) < 5E-4

sint = Integrator(sode, getTableauTripleJump(), Δt)
ssol = integrate(sint, nt)

@test rel_err(ssol.q, refx) < 5E-6

sint = Integrator(sode, getTableauSuzukiFractal(), Δt)
ssol = integrate(sint, nt)

@test rel_err(ssol.q, refx) < 5E-7


text:
using ReverseDiff: SpecialInstruction, istracked, increment_deriv!, deriv, TrackedArray
using ReverseDiff: InstructionTape, TrackedReal, seed!, unseed!, reverse_pass!, record!
import ReverseDiff: special_reverse_exec!, track

# if a value can't be tracked, return the untracked value instead, silently
track(value, tape::InstructionTape) = value

mutable struct GFBackpropParamsState
    trace::GFTrace
    score::TrackedReal
    tape::InstructionTape
    visitor::AddressVisitor
    tracked_params::Dict{Symbol,Any}
end

function GFBackpropParamsState(trace::GFTrace, tape, gf::GenFunction)
    tracked_params = Dict{Symbol,Any}()
    for (name, value) in gf.params
        tracked_params[name] = track(value, tape)
    end
    score = track(0., tape)
    GFBackpropParamsState(trace, score, tape, AddressVisitor(), tracked_params)
end

get_args_change(state::GFBackpropParamsState) = nothing
get_addr_change(state::GFBackpropParamsState, addr) = nothing
set_ret_change!(state::GFBackpropParamsState, value) = begin end

function read_param(state::GFBackpropParamsState, name::Symbol)
    value = state.tracked_params[name]
    value
end

function addr(state::GFBackpropParamsState, dist::Distribution{T}, args, addr) where {T}
    visit!(state.visitor, addr)
    call::CallRecord = get_primitive_call(state.trace, addr)
    retval::T = call.retval
    state.score += logpdf(dist, retval, args...)
    retval
end

struct BackpropParamsRecord
    generator::Generator
    subtrace::Any
end

function addr(state::GFBackpropParamsState, gen::Generator{T}, args, addr) where {T}
    visit!(state.visitor, addr)
    subtrace = get_subtrace(state.trace, addr)
    call::CallRecord = get_call_record(subtrace) # use the return value recorded in the trace
    retval::T = call.retval
    retval_maybe_tracked = track(retval, state.tape)
    if accepts_output_grad(gen) && !istracked(retval_maybe_tracked)
        error("Could not track return value at address $addr on AD tape.")
    end
    # some of the args may be tracked (see special_reverse_exec!)
    # note: we still need to run backprop_params on gen, even if it does not
    # accept an output gradient, because it may make random choices.
    record!(state.tape, SpecialInstruction,
        BackpropParamsRecord(gen, subtrace), (args...,), retval_maybe_tracked)
    retval_maybe_tracked
end

function splice(state::GFBackpropParamsState, gf::GenFunction, args::Tuple)
    exec(gf, state, args)
end

function maybe_track(arg, has_argument_grad::Bool, tape)
    has_argument_grad ? track(arg, tape) : arg
end

function backprop_params(gf::GenFunction, trace::GFTrace, retval_grad)
    tape = InstructionTape()
    state = GFBackpropParamsState(trace, tape, gf)
    call = get_call_record(trace)
    args = call.args
    args_maybe_tracked = (map(maybe_track, args, gf.has_argument_grads, fill(tape, length(args)))...,)
    retval_maybe_tracked = exec(gf, state, args_maybe_tracked)
    if gf.accepts_output_grad
        if !istracked(retval_maybe_tracked)
            error("Output of $gf accepts gradient but return value was not tracked on AD tape.")
        end
        if retval_grad != nothing
            # it is not an error if retval_grad is nothing, it means the gradient is zero
            # this is a convenience that
            deriv!(retval_maybe_tracked, retval_grad)
        end
    end
    seed!(state.score)
    reverse_pass!(tape)

    # increment the gradient accumulators for static parameters
    for (name, tracked) in state.tracked_params
        gf.params_grad[name] += deriv(tracked)
    end

    # return gradients with respect to inputs
    # NOTE: if a value isn't tracked the gradient is nothing
    input_grads::Tuple = (map((arg, has_grad) -> has_grad ? deriv(arg) : nothing,
                             args_maybe_tracked, gf.has_argument_grads)...,)
    input_grads
end

@noinline function special_reverse_exec!(instruction::SpecialInstruction{BackpropParamsRecord})
    record = instruction.func
    gen = record.generator
    args_maybe_tracked = instruction.input
    retval_maybe_tracked = instruction.output
    if accepts_output_grad(gen)
        @assert istracked(retval_maybe_tracked)
        retval_grad = deriv(retval_maybe_tracked)
    else
        retval_grad = nothing
    end
    arg_grads = backprop_params(gen, record.subtrace, retval_grad)
    for (arg, grad, has_grad) in zip(args_maybe_tracked, arg_grads, has_argument_grads(gen))
        if has_grad && istracked(arg)
            increment_deriv!(arg, grad)
        end
    end
    nothing
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
"""
    time_change_null(events, maxT)

Calculate the residuals for some events from a constant Poisson process.

# Arguments
* `events::Array{<:Number}` the events from the process.
* `maxT::Number` the maximum window time.

"""
function time_change_null(events::Array{<:Number}, maxT::Number)
    rate = length(events) / maxT
    Lambda = rate .* events
    _time_change(Lambda)
end

"""
    time_change_hawkes(events, bg, kappa, kern)

Calculate the residuals using the time change theorem for a Hawkes process.

# Arguments
* `events::Array{<:Number}` the events from the process.
* `bg` the background rate of the Hawkes process.
* `kappa` the kappa parameter of the Hawkes process.
* `kern` the kernel distriubtion of the Hawkes process.

"""
function time_change_hawkes(events::Array{<:Number}, bg::Number, kappa::Number, kern::Distributions.Distribution)

    Lambda = Array{Float64, 1}(undef, length(events))

    for (i, event) in enumerate(events)
        Lambda[i] = compensator(event, events, bg, kappa, kern)
    end

    _time_change(Lambda)
end

function _time_change(Lambda::Array{<:Number})

    tau = diff(Lambda)

    u = 1 .- exp.(-tau)
    b = (collect(1:length(u)) .- 0.5) / length(u)

    return (b, sort(u))
end


text:
using Juxta
using Test

ja = JuxtArray(randn(5,10), ["x","y"],
                     Dict("x"=>collect(1:5),"y"=>collect(1:10)))
ja1 = JuxtArray(randn(5,10), ["x","y"],
                     Dict("x"=>collect(1:5),"y"=>collect(1:10) .* 2))
ja2 = JuxtArray(randn(5,10), ["x","y"],
                     Dict("x"=>collect(1:5),"y"=>collect(1:10) .* 2))
ja3 = JuxtArray(randn(5,10), ["x","y"],
                     Dict("x"=>collect(1:5),"y"=>collect(1:10) .* 2))
ja4 = JuxtArray(randn(5,10,1), ["x","y","za"],
                     Dict("x"=>collect(1:5),"y"=>collect(1:10) .* 2,"za"=>[2]))
ja5 = JuxtArray(randn(5,10,1), ["x","y","za"],
                     Dict("x"=>collect(1:5),"y"=>collect(1:10) .* 2,"za"=>[2]))
@testset "juxta.jl" begin
    @test typeof(ja) == JuxtArray
    @test ja.indices["x"] == 1:5
    @test ja.indices["y"] == 1:10
    @test ja1.indices["x"] == 1:5
    @test ja1.indices["y"] == 1:10
    @test isel!(ja, x=2:4).indices["x"] == 1:3
    @test isel!(ja, y=2).indices["y"] == 1:1
    @test ja.coords["y"][1] == 2
    @test size(ja.array,1) == 3
    @test isel!(ja1, x=2:4, y=3:2:8).indices["x"] == 1:3
    @test ja1.indices["y"] == 1:3
    @test ja1.coords["y"] == [3,5,7] .* 2
    @test ja1.coords["x"] == collect(2:4)
    @test sel!(ja2, x=1.1:0.2:4.1).indices["x"] == 1:3
    @test ja2.coords["x"] == collect(2:4)
    @test sel!(ja3, "nearest", x=1.1:0.2:4.1).indices["x"] == 1:4
    @test ja3.coords["x"] == collect(1:4)
    @test sel!(ja3, "nearest", x=0.8).indices["x"] == 1:1
    @test ja3.coords["x"] == Number[1]
    @test sel!(ja3, y=10.3).indices["y"] == 1:2
    @test ja3.coords["y"] == Number[10,12]
    @test size(ja3) == (1,2)
    @test size(ja3, "x") == 1
    @test size(ja3, "y") == 2
    @test size(dropdims(ja3, ["x"])) == (2,)
    @test ja3.dims == ["y"]
    @test (ja4
           |> j->isel!(j,x=2)
           |> j->dropdims(j,["x","za"])
           |> j->size(j)) == (10,)
    @test ja4.dims == ["y"]
    @test typeof(show(ja4)) == Nothing
    @test ja5(x=1.1:0.2:4.1).indices["x"] == 1:3
end