text:
## all_values
## all_pairs
## all_triples
## all_n_tuples
## values_excursion
## pairs_excursion
## n_tuples_excursion
## full_factorial
## driven_factorial

"""
In-parameter-order General (IPOG).

This algorithm generates test cases quickly. It will generate
tuples of any order. It always generates the same set of test
cases for the same set of input values.

Lei, Yu, Raghu Kacker, D. Richard Kuhn, Vadim Okun, and James Lawrence. 2008. “IPOG/IPOG-D: Efficient Test Generation for Multi-Way Combinatorial Testing.” Software Testing, Verification & Reliability 18 (3): 125–48.
"""
struct IPOG
end

"""
Greedy Non-deterministic (GND).

This algorithm searches for test cases. It will generate tuples
of any order. It generates a different set every time it is invoked.

# Arguments

- `rng::Random.AbstractRNG`: This option is a random number
  generator. Set this if you want to generate the same test cases
  twice in a row.
- `M::Int`: The number of times it should create a candidate test case
  each time it creates a candidate. The default is 50. Raising this
  number could improve test cases and slow generation.

# Extended

This algorithm starts with the seeded test cases and then
adds test cases, one at a time. It chooses those parameters
that are least used, so far, and then chooses values of those
parameters that are least covered by previous tuples. At each
step, there is some probability of choosing among nearly-equal
next values.

It's not complicated, but it can be slow because every next choice
checks against all possible tuples. For large numbers of parameters
or large numbers of possible values of each parameter, this generator
can be slow, so test it for fewer values first and gradually increase
the number of parameters or parameter values.
"""
struct GND
    rng::Random.AbstractRNG
    M::Int

    function GND(; rng = nothing, M = 50)
        if rng === nothing
            new(Random.MersenneTwister(), M)
        else
            new(rng, M)
        end
    end
end

"""
This class requests tests that are excursions from a base case.
"""
struct Excursion
end

"""
The user specifies forbidden test cases in terms of parameter values,
but inner code thinks of parameter values as integers, so this function
wraps the user's function in order to do the translation.
"""
function wrap_disallow(disallow, parameters)
    if disallow !== nothing
        inner_filter = let filter = disallow, params = parameters
            choices -> begin
                as_values = ((c != 0 ? p[c] : nothing) for (p, c) in zip(params, choices))
                filter(as_values...)
            end
        end
    else
        inner_filter = nothing
    end
end

function seeds_to_integers(seed, parameters, Counter = Int)
    param_cnt = length(parameters)
    if length(seed) > 0
        seed_array = zeros(Counter, param_cnt, length(seed))
        for seed_idx in 1:length(seed)
            seed_array[:, seed_idx] = [indexin([c], p)[1] for (p, c) in zip(parameters, seed[seed_idx])]
        end
    else
        seed_array = []
    end
    seed_array
end

function generate_tuples(engine::IPOG, n_way, parameters, disallow, seeds, wayness, Counter)
    # We convert from parameters to integers here so that different generators
    # can use different internal representations.
    arity = Counter[length(p) for p in parameters]
    param_cnt = length(arity)
    if disallow !== nothing
        disallow_integer = wrap_disallow(disallow, parameters)
    else
        disallow_integer = x -> false
    end
    if seeds !== nothing && length(seeds) > 0
        seeds_int = seeds_to_integers(seeds, parameters, Counter)
    else
        seeds_int = zeros(Counter, param_cnt, 0)
    end
    if wayness !== nothing
        result = ipog_multi_way(arity, n_way, wayness, disallow_integer, seeds_int)
    elseif disallow !== nothing || seeds !== nothing
        result = ipog_multi(arity, n_way, disallow_integer, seeds_int)
    else
        result = ipog(arity, n_way)
    end
    [[p[c] for (p, c) in zip(parameters, result[:, i])] for i in 1:size(result, 2)]
end

function generate_tuples(engine::GND, n_way, parameters, disallow, seeds, wayness, Counter)
    arity = Counter[length(p) for p in parameters]
    if disallow !== nothing
        disallow_integer = wrap_disallow(disallow, parameters)
    else
        disallow_integer = x -> false
    end
    if seeds !== nothing && length(seeds) > 0
        seeds_int = seeds_to_integers(seeds, parameters, Counter)
    else
        seeds_int = zeros(Int, length(arity), 0)
    end
    if wayness !== nothing
        mwc = multi_way_coverage(arity, wayness, n_way)
        mc = UnitTestDesign.MatrixCoverage(mwc, size(mwc, 1), arity)
        result = n_way_coverage_multi(mc, disallow_integer, seeds_int, engine.M, engine.rng)
    else
        result = n_way_coverage_filter(arity, n_way, disallow_integer, seeds_int, engine.M, engine.rng)
    end
    [[p[c] for (p, c) in zip(parameters, answer)] for answer in result]
end

function generate_tuples(engine::Excursion, n_way, parameters, disallow, seeds, wayness, Counter)
    # We convert from parameters to integers here so that different generators
    # can use different internal representations.
    arity = Counter[length(p) for p in parameters]
    param_cnt = length(arity)
    if disallow !== nothing
        disallow_integer = wrap_disallow(disallow, parameters)
    else
        disallow_integer = x -> false
    end
    if seeds !== nothing && length(seeds) > 0
        seeds_int = seeds_to_integers(seeds, parameters, Counter)
    else
        seeds_int = zeros(Counter, param_cnt, 0)
    end
    if wayness !== nothing
        result = build_excursion_multi(arity, n_way, wayness, disallow_integer, seeds_int)
    else
        result = build_excursion(arity, n_way, disallow_integer, seeds_int)
    end
    [[p[c] for (p, c) in zip(parameters, result[:, i])] for i in 1:size(result, 2)]
end

"""
    all_tuples(parameters...; n_way, engine, disallow, seeds, wayness, Counter)

Given a tuple of parameters, generate all test cases that cover all `n_way` combinations
of those parameters.

# Arguments

- `engine=IPOG()`: The `engine` is `IPOG()`, `GND()` or `Excursion()`.
- `disallow=nothing`: The disallow function is a function of the parameters that
  returns `true` when that combination should be forbidden.
- `seeds=[]`: is a list of test cases that must be included among those
generated.
- `wayness` is a dictionary that specifies subsets of parameters for
  which to increase the `n_way` combinations. For instance, if the combinations are
  two-way, and you want the third-sixth parameters to be three-way covered, use,
  `wayness = Dict(3 => [[3, 4, 5, 6]])`.
- `Counter=Int`The Counter is an integer type to use for
the computation. It must be large enough to hold the integer number of the parameters.

# Examples
```julia
all_tuples([1, 2, 3], ["a", "b", "c"], [true, false]; n_way = 2)
parameters = fill(collect(1:3), 10)
all_tuples(parameters...; n_way = 4, engine = Excursion())
```
"""
function all_tuples(
    parameters...;
    n_way::Integer = 2, engine = IPOG(), disallow = nothing, seeds = nothing, wayness = nothing, Counter = Int
    )
    if length(parameters) < 2
        throw(DomainError(
            parameters,
            "Arguments should be lists of parameter values for more than one parameter"))
    end
    for param in parameters
        if length(param) < 2
            throw(DomainError(param, "Each argument should be a list of parameter values"))
        end
    end
    if disallow !== nothing && !isa(disallow, Function)
        throw(DomainError(disallow, "The disallow argument should be a function that
        returns true or false to exclude or include a test case."))
    end
    if seeds !== nothing && length(seeds) > 0
        if !all(length(s) .== length(parameters) for s in seeds)
            throw(DomainError(seeds, "The seeded test cases should have the same number
            of parameter values as the passed parameter values."))
        end
    end
    if wayness !== nothing
        keys_are_ints = all(isa(k, Int) for k in keys(wayness))
        if !keys_are_ints
            throw(DomainError(wayness, "The wayness argument is a dictionary with integer keys."))
        end
        values_are_lists = all(isa(v, AbstractArray) for v in values(wayness))
        lists_of_lists = all(isa(first(v), AbstractArray) for v in values(wayness))
        if !values_are_lists || !lists_of_lists
            throw(DomainError(wayness, "The wayness argument is a dictionary with values
            that are lists of lists of parameters."))
        end
    end
    generate_tuples(engine, n_way, parameters, disallow, seeds, wayness, Counter)
end

"""
    all_values(parameters...; kwargs...)

Ensure that the test cases include every value of every parameter
at least once.

See also: [`all_tuples`](@ref)
"""
function all_values(parameters...; kwargs...)
    all_tuples(parameters...; n_way = 1, kwargs...)
end

"""
    all_pairs(parameters...; kwargs...)

Ensure that the returned test cases include every pair of
parameters at least once.

# Examples
```julia
all_pairs([1, 2, 3], ["a", "b", "c"], [true, false])
```

See also: [`all_tuples`](@ref)
"""
function all_pairs(parameters...; kwargs...)
    all_tuples(parameters...; n_way = 2, kwargs...)
end

"""
    all_triples(parameters...; kwargs...)

Ensure that the returned test cases include every combination
of three parameters at least once.

See also: [`all_tuples`](@ref)
"""
function all_triples(parameters...; kwargs...)
    all_tuples(parameters...; n_way = 3, kwargs...)
end

"""
    values_excursion(parameters...; kwargs...)

This starts with the first choice for each of the parameters.
It creates test cases by varying each parameter, one at a time,
through its possible values.

# Examples
```julia
values_excursion([:a, :b, :c], [1, 2, 3])
```

See also: [`all_tuples`](@ref)
"""
function values_excursion(parameters...; kwargs...)
    all_tuples(parameters...; n_way = 1, engine = Excursion(), kwargs...)
end

"""
    pairs_excursion(parameters...; kwargs...)

This starts with the first choice for each of the parameters.
It creates test cases by varying each parameter, one at a time,
through its possible values. Then it walks pairs of parameters
away from the base case.

# Examples
```julia
pairs_excursion([:a, :b], [1, 2], [1, 2], ["a", "b"])
```

See also: [`all_tuples`](@ref)
"""
function pairs_excursion(parameters...; kwargs...)
    all_tuples(parameters...; n_way = 2, engine = Excursion(), kwargs...)
end

"""
    pairs_excursion(parameters...; kwargs...)

This starts with the first choice for each of the parameters.
It creates test cases by varying each parameter, one at a time,
through its possible values. Then it walks pairs of parameters
away from the base case, and finally triples.

See also: [`all_tuples`](@ref)
"""
function triples_excursion(parameters...; kwargs...)
    all_tuples(parameters...; n_way = 3, engine = Excursion(), kwargs...)
end

"""
    full_factorial(parameters...)
    full_factorial(parameters...; disallow = filter_function)

Generates a test case for every combination of the parameters.

# Examples
```julia
full_factorial([0.1, 0.2, 0.3], ["low", "high"], [false, true])
```

If you specify a filter function, it will remove combinations
that are disallowed.
"""
function full_factorial(parameters...; disallow = nothing)
    arity = [length(p) for p in parameters]
    param_cnt = length(arity)
    if disallow !== nothing
        disallow = wrap_disallow(disallow, parameters)
    else
        disallow = x -> false
    end
    result = full_factorial(arity, disallow)
    [[p[c] for (p, c) in zip(parameters, result[:, i])] for i in 1:size(result, 2)]
end


text:
using JuMP: REPLMode, IJuliaMode

# Test string creation
@testset "String Creators" begin
    # initialize model and attributes
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par1 <= 1)
    @infinite_parameter(m, pars[1:2] in MvNormal([1, 1], 1))
    @infinite_variable(m, x(par1))
    @infinite_variable(m, z(pars))
    @hold_variable(m, y)
    @objective(m, Min, 2 + y)
    @constraint(m, c1, x + y -2 <= 0)
    @constraint(m, c2, y^2 - 3 == 0)
    @BDconstraint(m, c3(par1 in [0, 0.5]), x == 5)
    # test in_set_string (IntervalSet)
    @testset "JuMP.in_set_string (Interval)" begin
        # test simple case
        set = IntervalSet(0, 1)
        str = JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        @test in_set_string(IJuliaMode, set) == str
        # test rounding case
        set = IntervalSet(-0, 1)
        str = JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        @test in_set_string(IJuliaMode, set) == str
        # test decimal case
        set = IntervalSet(0.1, 1.3)
        str = JuMP._math_symbol(REPLMode, :in) * " [0.1, 1.3]"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :in) * " [0.1, 1.3]"
        @test in_set_string(IJuliaMode, set) == str
        # test finite case
        set = IntervalSet(0.1, 0.1)
        str = JuMP._math_symbol(REPLMode, :eq) * " 0.1"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :eq) * " 0.1"
        @test in_set_string(IJuliaMode, set) == str
    end
    # test in_set_string (Distribution)
    @testset "JuMP.in_set_string (Distribution)" begin
        # test univariate set
        set = DistributionSet(Uniform())
        str = JuMP._math_symbol(REPLMode, :in) * " Uniform(a=0.0, b=1.0)"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :in) * " Uniform(a=0.0, b=1.0)"
        @test in_set_string(IJuliaMode, set) == str
        # test mulivariate set
        set = DistributionSet(MvNormal([1], 1))
        str = JuMP._math_symbol(REPLMode, :in) * " IsoNormal(dim: 1, μ: " *
                 "[1.0], Σ: [1.0])"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :in) * " IsoNormal(dim: 1, μ: " *
                 "[1.0], Σ: [1.0])"
        @test in_set_string(IJuliaMode, set) == str
    end
    # test in_set_string (Fallback)
    @testset "JuMP.in_set_string (Fallback)" begin
        set = BadSet()
        in1 = JuMP._math_symbol(REPLMode, :in)
        in2 = JuMP._math_symbol(IJuliaMode, :in)
        @test in_set_string(REPLMode, set) == in1 * " BadSet()"
        @test in_set_string(IJuliaMode, set) == in2 * " BadSet()"
    end
    # test bound_string
    @testset "bound_string" begin
        # test with single bound
        bounds = ParameterBounds(Dict(par1 => IntervalSet(0.5, 0.7)))
        str = "par1 " * JuMP._math_symbol(REPLMode, :in) * " [0.5, 0.7]"
        @test InfiniteOpt.bound_string(REPLMode, bounds) == str
        str = "par1 " *  JuMP._math_symbol(IJuliaMode, :in) * " [0.5, 0.7]"
        @test InfiniteOpt.bound_string(IJuliaMode, bounds) == str
    end
    # test constraint_string (bounded constraints)
    @testset "JuMP.constraint_string (Bounded)" begin
        str = "x(par1) " * JuMP._math_symbol(REPLMode, :eq) * " 5.0, " *
              JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 0.5]"
        @test constraint_string(REPLMode, constraint_object(c3)) == str
        str =  "x(par1) " * JuMP._math_symbol(IJuliaMode, :eq) * " 5.0, " *
               JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
               JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5]"
        @test constraint_string(IJuliaMode, constraint_object(c3)) == str
    end
    # test constraint_string of GeneralConstraintRefs
    @testset "JuMP.constraint_string (Reference)" begin
        # test normal constraint
        str = "c1 : x(par1) + y " * JuMP._math_symbol(REPLMode, :leq) *
              " 2.0"
        @test constraint_string(REPLMode, c1) == str
        str = "c1 : \$ x(par1) + y " * JuMP._math_symbol(IJuliaMode, :leq) *
              " 2.0 \$"
        @test constraint_string(IJuliaMode, c1) == str
        # test bounded constraint
        str = "c3 : x(par1) " * JuMP._math_symbol(REPLMode, :eq) * " 5.0, " *
              JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 0.5]"
        @test constraint_string(REPLMode, c3) == str
        str =  "c3 : \$ x(par1) " * JuMP._math_symbol(IJuliaMode, :eq) * " 5.0, " *
               JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
               JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5] \$"
        @test constraint_string(IJuliaMode, c3) == str
    end
    # test constraints_string
    @testset "JuMP.constraints_string" begin
        # test REPLMode
        strings = Vector{String}(undef, 5)
        strings[1] = "x(par1) + y " * JuMP._math_symbol(REPLMode, :leq) * " 2.0"
        strings[2] = "y² " * JuMP._math_symbol(REPLMode, :eq) * " 3.0"
        strings[3] = "x(par1) " * JuMP._math_symbol(REPLMode, :eq) * " 5.0, " *
                     JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
                     JuMP._math_symbol(REPLMode, :in) * " [0, 0.5]"
        strings[4] = "par1 " * JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        strings[5] = "pars " * JuMP._math_symbol(REPLMode, :in) *
                     " IsoNormal(dim: 2, μ: [1.0, 1.0], Σ: [1.0 0.0; 0.0 1.0])"
        @test constraints_string(REPLMode, m) == strings
        # test IJuliaMode
        strings = Vector{String}(undef, 5)
        strings[1] = "x(par1) + y " * JuMP._math_symbol(IJuliaMode, :leq) * " 2.0"
        strings[2] = "y^2 " * JuMP._math_symbol(IJuliaMode, :eq) * " 3.0"
        strings[3] = "x(par1) " * JuMP._math_symbol(IJuliaMode, :eq) * " 5.0, " *
                     JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
                     JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5]"
        strings[4] = "par1 " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        strings[5] = "pars " * JuMP._math_symbol(IJuliaMode, :in) *
                     " IsoNormal(dim: 2, μ: [1.0, 1.0], Σ: [1.0 0.0; 0.0 1.0])"
        @test constraints_string(IJuliaMode, m) == strings
    end
    # test objective_function_string
    @testset "JuMP.objective_function_string" begin
        @test objective_function_string(REPLMode, m) == "y + 2"
        @test objective_function_string(IJuliaMode, m) == "y + 2"
    end
end

# Helper function to test IO methods work correctly
function show_test(mode, obj, exp_str; repl=:both)
    if mode == REPLMode
        repl != :show  && @test sprint(print, obj) == exp_str
        repl != :print && @test sprint(show,  obj) == exp_str
    else
        @test sprint(show, "text/latex", obj) == exp_str
    end
end

# Make another helper function for other io methods
io_test(f::Function, exp_str::String, args...) = begin
    io = IOBuffer()
    f(io, args...)
    @test String(take!(io)) == exp_str
end

# test infinite model show
@testset "Show InfiniteModel" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par1 <= 1)
    @infinite_parameter(m, pars[1:2] in MvNormal([1, 1], 1))
    @infinite_variable(m, x(par1))
    @infinite_variable(m, z(pars))
    @hold_variable(m, y)
    @objective(m, Min, 2 + y)
    @constraint(m, c1, x + y -2 <= 0)
    @constraint(m, c2, y^2 - 3 == 0)
    bounds = ParameterBounds(Dict(par1 => IntervalSet(0.1, 1)))
    @BDconstraint(m, c3(par1 in [0, 0.5]), x == 5)
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    # test Base.show (constraint in REPL)
    @testset "Base.show (REPL Constraint)" begin
        # test normal
        str = "c1 : x(par1) + y " * JuMP._math_symbol(REPLMode, :leq) *
              " 2.0"
        show_test(REPLMode, c1, str)
        # test bounded
        str = "c3 : x(par1) " * JuMP._math_symbol(REPLMode, :eq) * " 5.0, " *
              JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 0.5]"
        show_test(REPLMode, c3, str)
    end
    # test Base.show (constraint in IJulia)
    @testset "Base.show (IJulia Constraint)" begin
        # test normal
        str = "c1 : \$ x(par1) + y " * JuMP._math_symbol(IJuliaMode, :leq) *
              " 2.0 \$"
        show_test(IJuliaMode, c1, str)
        # test bounded
        str =  "c3 : \$ x(par1) " * JuMP._math_symbol(IJuliaMode, :eq) * " 5.0, " *
               JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
               JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5] \$"
        show_test(IJuliaMode, c3, str)
    end
    # test Base.show (ParameterBounds in REPL)
    @testset "Base.show (REPL ParameterBounds)" begin
        # test normal
        str = "Subdomain bounds (1): par1 " * JuMP._math_symbol(REPLMode, :in) *
              " [0.1, 1]"
        show_test(REPLMode, bounds, str)
    end
    # test Base.show (ParameterBounds in IJulia)
    @testset "Base.show (IJulia ParameterBounds)" begin
        # test normal
        str = "Subdomain bounds (1): par1 " * JuMP._math_symbol(IJuliaMode, :in) *
              " [0.1, 1]"
        show_test(IJuliaMode, bounds, str)
    end
    # test Base.show (IntervalSet in REPL)
    @testset "Base.show (REPL IntervalSet)" begin
        # test normal
        show_test(REPLMode, IntervalSet(0, 1), "[0, 1]")
    end
    # test Base.show (IntervalSet in IJulia)
    @testset "Base.show (IJulia IntervalSet)" begin
        # test normal
        show_test(IJuliaMode, IntervalSet(0, 1), "[0, 1]")
    end
    # test Base.show (DistributionSet in REPL)
    @testset "Base.show (REPL IntervalSet)" begin
        # test normal
        show_test(REPLMode, DistributionSet(Uniform()),
                  "Uniform{Float64}(a=0.0, b=1.0)")
    end
    # test Base.show (DistributionSet in IJulia)
    @testset "Base.show (IJulia IntervalSet)" begin
        # test normal
        show_test(IJuliaMode, DistributionSet(Uniform()),
                  "Uniform{Float64}(a=0.0, b=1.0)")
    end
    # test show_backend_summary
    @testset "JuMP.show_backend_summary" begin
        # test without optimizer
        str = "Optimizer model backend information: \nModel mode: AUTOMATIC\n" *
              "CachingOptimizer state: NO_OPTIMIZER\nSolver name: No optimizer" *
              " attached."
        io_test(show_backend_summary, str, m)
        # test with optimizer
        set_optimizer(optimizer_model(m), mockoptimizer)
        str = "Optimizer model backend information: \nModel mode: AUTOMATIC\n" *
              "CachingOptimizer state: EMPTY_OPTIMIZER\nSolver name: Mock"
        io_test(show_backend_summary, str, m)
    end
    # test show_objective_function_summary
    @testset "JuMP.show_objective_function_summary" begin
        str = "Objective function type: GenericAffExpr{Float64,HoldVariableRef}\n"
        io_test(show_objective_function_summary, str, m)
    end
    # test show_constraints_summary
    @testset "JuMP.show_constraints_summary" begin
        # test _plural
        @test InfiniteOpt._plural(1) == ""
        @test InfiniteOpt._plural(2) == "s"
        # test the main function
        str = "`GenericAffExpr{Float64,GeneralVariableRef}`-in-`MathOptInter" *
              "face.LessThan{Float64}`: 1 constraint\n`GenericQuadExpr{" *
              "Float64,HoldVariableRef}`-in-`MathOptInterface.EqualTo{" *
              "Float64}`: 1 constraint\n`GenericAffExpr{Float64,Infinite" *
              "VariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 1 " *
              "constraint\n"
        io_test(show_constraints_summary, str, m)
    end
    # test show_objective_function_summary
    @testset "Base.show (InfiniteModel)" begin
        # test minimization
        str = "An InfiniteOpt Model\nMinimization problem with:\nVariables: " *
              "3\nObjective function type: GenericAffExpr{Float64,Hold" *
              "VariableRef}\n`GenericAffExpr{Float64,GeneralVariableRef}`-in-" *
              "`MathOptInterface.LessThan{Float64}`: 1 constraint\n`Generic" *
              "QuadExpr{Float64,HoldVariableRef}`-in-`MathOptInterface.Equal" *
              "To{Float64}`: 1 constraint\n`GenericAffExpr{Float64,Infinite" *
              "VariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 1 " *
              "constraint\nNames registered in the model: c1, c2, c3, par1, " *
              "pars, x, y, z\nOptimizer model backend information: \nModel " *
              "mode: AUTOMATIC\nCachingOptimizer state: EMPTY_OPTIMIZER\n" *
              "Solver name: Mock"
        show_test(REPLMode, m, str, repl=:show)
        # test maximization
        set_objective_sense(m, MOI.MAX_SENSE)
        str = "An InfiniteOpt Model\nMaximization problem with:\nVariables: " *
              "3\nObjective function type: GenericAffExpr{Float64,Hold" *
              "VariableRef}\n`GenericAffExpr{Float64,GeneralVariableRef}`-in-" *
              "`MathOptInterface.LessThan{Float64}`: 1 constraint\n`Generic" *
              "QuadExpr{Float64,HoldVariableRef}`-in-`MathOptInterface.Equal" *
              "To{Float64}`: 1 constraint\n`GenericAffExpr{Float64,Infinite" *
              "VariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 1 " *
              "constraint\nNames registered in the model: c1, c2, c3, par1, " *
              "pars, x, y, z\nOptimizer model backend information: \nModel " *
              "mode: AUTOMATIC\nCachingOptimizer state: EMPTY_OPTIMIZER\n" *
              "Solver name: Mock"
        show_test(REPLMode, m, str, repl=:show)
        # test feasibility
        set_objective_sense(m, MOI.FEASIBILITY_SENSE)
        str = "An InfiniteOpt Model\nFeasibility problem with:\nVariables: 3" *
              "\n`GenericAffExpr{Float64,GeneralVariableRef}`-in-`MathOpt" *
              "Interface.LessThan{Float64}`: 1 constraint\n`GenericQuadExpr{" *
              "Float64,HoldVariableRef}`-in-`MathOptInterface.EqualTo{" *
              "Float64}`: 1 constraint\n`GenericAffExpr{Float64,Infinite" *
              "VariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 1 " *
              "constraint\nNames registered in the model: c1, c2, c3, par1, " *
              "pars, x, y, z\nOptimizer model backend information: \nModel " *
              "mode: AUTOMATIC\nCachingOptimizer state: EMPTY_OPTIMIZER\n" *
              "Solver name: Mock"
        show_test(REPLMode, m, str, repl=:show)
    end
end

ext: jl


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	co2_cap(EP::Model, inputs::Dict, setup::Dict)

This policy constraints mimics the CO$_2$ emissions cap and permit trading systems, allowing for emissions trading across each zone for which the cap applies. The constraint $p \in \mathcal{P}^{CO_2}$ can be flexibly defined for mass-based or rate-based emission limits for one or more model zones, where zones can trade CO$_2$ emissions permits and earn revenue based on their CO$_2$ allowance. Note that if the model is fully linear (e.g. no unit commitment or linearized unit commitment), the dual variable of the emissions constraints can be interpreted as the marginal CO$_2$ price per tonne associated with the emissions target. Alternatively, for integer model formulations, the marginal CO$_2$ price can be obtained after solving the model with fixed integer/binary variables.

The CO$_2$ emissions limit can be defined in one of the following ways: a) a mass-based limit defined in terms of annual CO$_2$ emissions budget (in million tonnes of CO2), b) a load-side rate-based limit defined in terms of tonnes CO$_2$ per MWh of demand and c) a generation-side rate-based limit defined in terms of tonnes CO$_2$ per MWh of generation.

**Mass-based emissions constraint**

Mass-based emission limits are implemented in the following expression. For each constraint, $p \in \mathcal{P}^{CO_2}_{mass}$, we define a set of zones $z \in \mathcal{Z}^{CO_2}_{p,mass}$ that can trade CO$_2$ allowance. Input data for each constraint  $p \in \mathcal{P}^{CO_2}_{mass}$ requires the CO$_2$ allowance/ budget for each model zone, $\epsilon^{CO_{2}}_{z,p, mass}$, to be provided in terms of million metric tonnes. For every generator $y$, the parameter $\epsilon_{y,z}^{CO_2}$ reflects the specific $CO_2$ emission intensity in tCO$_2$/MWh associated with its operation.  The resulting constraint is given as:

```math
\begin{aligned}
    \sum_{z \in \mathcal{Z}^{CO_2}_{p,mass}} \sum_{y \in \mathcal{G}} \sum_{t \in \mathcal{T}} \left(\epsilon_{y,z}^{CO_2} \times \omega_{t} \times \Theta_{y,z,t} \right)
   & \leq \sum_{z \in \mathcal{Z}^{CO_2}_{p,mass}} \epsilon^{CO_{2}}_{z,p, mass} \hspace{1 cm}  \forall p \in \mathcal{P}^{CO_2}_{mass}
\end{aligned}
```

In the above constraint, we include both power discharge and charge term for each resource to account for the potential for CO$_2$ emissions (or removal when considering negative emissions technologies) associated with each step. Note that if a limit is applied to each zone separately, then the set $\mathcal{Z}^{CO_2}_{p,mass}$ will contain only one zone with no possibility of trading. If a system-wide emission limit constraint is applied, then $\mathcal{Z}^{CO_2}_{p,mass}$ will be equivalent to a set of all zones.

**Load-side rate-based emissions constraint**

We modify the right hand side of the above mass-based constraint, $p \in \mathcal{P}^{CO_2}_{load}$, to set emissions target based on a CO$_2$ emission rate limit in tCO$_2$/MWh $\times$ the total demand served in each zone. In the following constraint, total demand served takes into account non-served energy and storage related losses. Here, $\epsilon_{z,p,load}^{maxCO_2}$ denotes the emission limit in terms on tCO$_2$/MWh.

```math
\begin{aligned}
    \sum_{z \in \mathcal{Z}^{CO_2}_{p,load}} \sum_{y \in \mathcal{G}} \sum_{t \in \mathcal{T}} \left(\epsilon_{y,z}^{CO_2} \times \omega_{t} \times \Theta_{y,t,z} \right)
    \leq & \sum_{z \in \mathcal{Z}^{CO_2}_{p,load}} \sum_{t \in \mathcal{T}}  \left(\epsilon_{z,p,load}^{CO_2} \times  \omega_{t} \times D_{z,t} \right) \\  + & \sum_{z \in \mathcal{Z}^{CO_2}_{p,load}} \sum_{y \in \mathcal{O}}  \sum_{t \in \mathcal{T}} \left(\epsilon_{z,p,load}^{CO_2} \times \omega_{t} \times \left(\Pi_{y,t,z} - \Theta_{y,t,z} \right) \right) \\  - & \sum_{z \in \mathcal{Z}^{CO_2}_{p,load}} \sum_{s \in \mathcal{S} } \sum_{t \in \mathcal{T}}  \left(\epsilon_{z,p,load}^{CO_2} \times \omega_{t} \times \Lambda_{s,z,t}\right) \hspace{1 cm}  \forall p \in \mathcal{P}^{CO_2}_{load}
\end{aligned}
```

**Generator-side emissions rate-based constraint**

Similarly, a generation based emission constraint is defined by setting the emission limit based on the total generation times the carbon emission rate limit in tCO$_2$/MWh of the region. The resulting constraint is given as:

```math
\begin{aligned}
\sum_{z \in \mathcal{Z}^{CO_2}_{p,gen}} \sum_{y \in \mathcal{G}} \sum_{t \in \mathcal{T}} & \left(\epsilon_{y,z}^{CO_2} \times \omega_{t} \times \Theta_{y,t,z} \right) \\
    \leq \sum_{z \in \mathcal{Z}^{CO_2}_{p,gen}} \sum_{y \in \mathcal{G}} \sum_{t \in \mathcal{T}} & \left(\epsilon_{z,p,gen}^{CO_2} \times  \omega_{t} \times \Theta_{y,t,z} \right)  \hspace{1 cm}  \forall p \in \mathcal{P}^{CO_2}_{gen}
\end{aligned}
```

Note that the generator-side rate-based constraint can be used to represent a fee-rebate (``feebate'') system: the dirty generators that emit above the bar ($\epsilon_{z,p,gen}^{maxCO_2}$) have to buy emission allowances from the emission regulator in the region $z$ where they are located; in the same vein, the clean generators get rebates from the emission regulator at an emission allowance price being the dual variable of the emissions rate constraint.
"""
function co2_cap(EP::Model, inputs::Dict, setup::Dict)

	println("C02 Policies Module")

	dfGen = inputs["dfGen"]
	SEG = inputs["SEG"]  # Number of lines
	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones

	### Constraints ###

	## Mass-based: Emissions constraint in absolute emissions limit (tons)
	if setup["CO2Cap"] == 1
		@constraint(EP, cCO2Emissions_systemwide[cap=1:inputs["NCO2Cap"]],
			sum(inputs["omega"][t] * EP[:eEmissionsByZone][z,t] for z=findall(x->x==1, inputs["dfCO2CapZones"][:,cap]), t=1:T) <=
			sum(inputs["dfMaxCO2"][z,cap] for z=findall(x->x==1, inputs["dfCO2CapZones"][:,cap]))
		)

	## Load + Rate-based: Emissions constraint in terms of rate (tons/MWh)
	elseif setup["CO2Cap"] == 2 ##This part moved to non_served_energy.jl
		@constraint(EP, cCO2Emissions_systemwide[cap=1:inputs["NCO2Cap"]],
			sum(inputs["omega"][t] * EP[:eEmissionsByZone][z,t] for z=findall(x->x==1, inputs["dfCO2CapZones"][:,cap]), t=1:T) <=
			sum(inputs["dfMaxCO2Rate"][z,cap] * sum(inputs["omega"][t] * (inputs["pD"][t,z] - sum(EP[:vNSE][s,t,z] for s in 1:SEG)) for t=1:T) for z = findall(x->x==1, inputs["dfCO2CapZones"][:,cap])) +
			sum(inputs["dfMaxCO2Rate"][z,cap] * setup["StorageLosses"] *  EP[:eELOSSByZone][z] for z=findall(x->x==1, inputs["dfCO2CapZones"][:,cap]))
		)

	## Generation + Rate-based: Emissions constraint in terms of rate (tons/MWh)
	elseif (setup["CO2Cap"]==3)
		@constraint(EP, cCO2Emissions_systemwide[cap=1:inputs["NCO2Cap"]],
			sum(inputs["omega"][t] * EP[:eEmissionsByZone][z,t] for z=findall(x->x==1, inputs["dfCO2CapZones"][:,cap]), t=1:T) <=
			sum(inputs["dfMaxCO2Rate"][z,cap] * inputs["omega"][t] * EP[:eGenerationByZone][z,t] for t=1:T, z=findall(x->x==1, inputs["dfCO2CapZones"][:,cap]))
		)
	end

	return EP

end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
include("model_search/mountain_car_model_search_utils.jl")

struct MountainCarModelSearchAgent
    mountaincar::MountainCar
    model::MountainCar
    data::Array{MountainCarContTransition}
    optimization_params::MountainCarOptimizationParameters
    horizon::Int64
end

function return_based_model_search(
    mountaincar::MountainCar,
    data::Array{MountainCarContTransition},
    optimization_params::MountainCarOptimizationParameters,
    horizon::Int64;
    hardcoded::Bool = false,
    num_episodes_eval::Int64 = 1,
    debug::Bool = false,
    eval_distance::Bool = false,
    optimistic::Bool = false,
    model_evaluation::Bool = false,
)
    params = true_params
    least_squares_params = get_least_squares_fit(mountaincar, params, data)
    x_matrices_array, _, x_next_array, _, cost_array = preprocess_data(data)
    function eval_fn(p)
        policy, values, converged = value_iteration(mountaincar, p)
        if converged
            if optimistic
                return mfmc_optimistic_evaluation(
                    mountaincar,
                    policy,
                    p,
                    horizon,
                    x_matrices_array,
                    x_next_array,
                    cost_array,
                    num_episodes_eval,
                    0.015;
                    debug = debug,
                )
            elseif model_evaluation
                return mfmc_model_evaluation(
                    mountaincar,
                    policy,
                    p,
                    horizon,
                    x_matrices_array,
                    x_next_array,
                    cost_array,
                    num_episodes_eval,
                    0.01;
                    debug = debug,
                )
            else
                return mfmc_evaluation(
                    mountaincar,
                    policy,
                    values,
                    horizon,
                    x_matrices_array,
                    x_next_array,
                    cost_array,
                    num_episodes_eval;
                    hardcoded = hardcoded,
                    debug = debug,
                )
            end
        else
            if debug
                println("Value Iteration did not converge. Skipping parameters ", p)
            end
            return Inf
        end
    end

    params = hill_climb(
        eval_fn,
        optimization_params;
        least_squares_params = least_squares_params,
        debug = debug,
    )
    if eval_distance
        policy, values, _ = value_iteration(mountaincar, params)
        estimated_return, eval_distances = mfmc_evaluation(
            mountaincar, policy, values, horizon,
            x_matrices_array, x_next_array, cost_array,
            num_episodes_eval;
            eval_distance = eval_distance,
        )
        println("Mean distance ", mean(eval_distances),
                " Std distance ", std(eval_distances),
                " Max distance ", maximum(eval_distances),
                " Min distance ", minimum(eval_distances))
        println("Estimated return is ", estimated_return)
    end
    params
end

function bellman_based_model_search(
    mountaincar::MountainCar,
    data::Array{MountainCarContTransition},
    optimization_params::MountainCarOptimizationParameters,
    horizon::Int64;
    num_episodes_eval::Int64 = 1,
    debug::Bool = false,
)
    params = true_params
    least_squares_params = get_least_squares_fit(mountaincar, params, data)
    x_matrices_array, _, xnext_array, _, _ = preprocess_data(data)
    function eval_fn(p)
        policy, values, converged = value_iteration(mountaincar, p)
        if converged
            return bellman_evaluation(
                mountaincar,
                p,
                policy,
                values,
                horizon,
                x_matrices_array,
                xnext_array,
                num_episodes_eval;
                debug = debug,
            )
        else
            if debug
                println("Value Iteration did not converge. Skipping parameters ", p)
            end
            return Inf
        end
    end
    params = hill_climb(
        eval_fn,
        optimization_params;
        least_squares_params = least_squares_params,
        debug = debug,
    )
    params
end

function MountainCarModelSearchAgent(
    mountaincar::MountainCar,
    model::MountainCar,
    horizon::Int64,
    num_episodes_offline::Int64,
)
    data = generate_batch_data(mountaincar, true_params, num_episodes_offline, horizon)
    # println("Generated ", length(data), " transitions offline")
    MountainCarModelSearchAgent(mountaincar, model, horizon, data)
end

function MountainCarModelSearchAgent(
    mountaincar::MountainCar,
    model::MountainCar,
    horizon::Int64,
    data::Array{MountainCarContTransition},
)
    maximum_num_evaluations = 90
    optimization_params =
        MountainCarOptimizationParameters([0.0012, 0.5], [-0.0025, 3], maximum_num_evaluations, false)
    MountainCarModelSearchAgent(mountaincar, model, data, optimization_params, horizon)
end

function run_return_based_model_search(
    agent::MountainCarModelSearchAgent;
    ensemble = false,
    hardcoded = false,
    max_steps = 500,
    debug = false,
    eval_distance = false,
    optimistic = false,
    model_evaluation = false,
)
    params = return_based_model_search(
        agent.model,
        agent.data,
        agent.optimization_params,
        agent.horizon;
        hardcoded = hardcoded,
        debug = debug,
        eval_distance = eval_distance,
        optimistic = optimistic,
        model_evaluation = model_evaluation,
    )
    println("Params found is ", params)
    run(agent, params; max_steps = max_steps)
end

function run_bellman_based_model_search(
    agent::MountainCarModelSearchAgent;
    max_steps = 500,
    debug = false,
)
    params = bellman_based_model_search(
        agent.model,
        agent.data,
        agent.optimization_params,
        agent.horizon;
        debug = debug,
    )
    run(agent, params, max_steps = max_steps)
end

function run(
    agent::MountainCarModelSearchAgent,
    params::Array{Float64};
    max_steps = 500,
    debug = false,
    num_eval = 10,
)
    policy, _, _ = value_iteration(agent.model, params)
    rng = MersenneTwister(0)
    num_steps = 0.0
    for _ in 1:num_eval
        num_steps += length(
            simulate_episode(
                agent.mountaincar,
                true_params,
                agent.horizon;
                policy = policy,
                rng = rng,
            )
        )
    end
    num_steps/num_eval
end


text:
"""
Constructs SeqModel object from config dict.
"""
function construct_model(config::Dict,
                         max_time::Float64,
                         num_neurons::Int64)

    # Prior on sequence type proportions / relative frequencies.
    seq_type_proportions = SymmetricDirichlet(
        config[:seq_type_conc_param],
        config[:num_sequence_types]
    )

    # Prior on expected number of spikes induces by a sequence events.
    seq_event_amplitude = specify_gamma(
        config[:mean_event_amplitude],    # mean of gamma; α / β
        config[:var_event_amplitude]      # variance of gamma; α / β²
    )

    # Prior on relative response amplitudes per neuron to each sequence type.
    neuron_response_proportions = SymmetricDirichlet(
        config[:neuron_response_conc_param],
        num_neurons
    )

    # Prior on the response offsets and widths for each neuron.
    neuron_response_profile = NormalInvChisq(
        config[:neuron_offset_pseudo_obs],
        0.0, # prior mean
        config[:neuron_width_pseudo_obs],
        config[:neuron_width_prior],
    )

    # Prior on expected number of background spikes in a unit time interval.
    bkgd_amplitude = specify_gamma(
        config[:mean_bkgd_spike_rate],    # mean of gamma; α / β
        config[:var_bkgd_spike_rate]      # variance of gamma; α / β²
    )

    # Prior on relative background firing rates across neurons.
    bkgd_proportions = SymmetricDirichlet(
        config[:bkgd_spikes_conc_param],
        num_neurons
    )

    model = SeqModel(
        # constants
        max_time,
        config[:max_sequence_length],

        # warp parameters
        config[:num_warp_values],
        config[:max_warp],
        config[:warp_variance],

        # priors
        config[:seq_event_rate],
        seq_type_proportions,
        seq_event_amplitude,
        neuron_response_proportions,
        neuron_response_profile,
        bkgd_amplitude,
        bkgd_proportions
    )

    if (:num_threads in keys(config)) && config[:num_threads] > 0
        return DistributedSeqModel(model, config[:num_threads])
    else
        return model
    end
end


text:
module FBD

export example, MEPS, hypergraphsize,
    DCHSBM_sampler, Kronecker_sampler, Typing_sampler, ER_sampler,
    hyper_pa, er,
    AliasTable

include("alias_table.jl")
include("io.jl")

include("dchsbm/intensity_functions.jl")
include("dchsbm/model.jl")

include("kronecker/util.jl")
include("kronecker/model.jl")

include("hyper_pa.jl")

include("typing.jl")

include("ER.jl")
include("simple_er.jl")

include("examples.jl")

end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
nums = rand(100000)

for i = 1:length(nums)
    r = rand(rng)
    nums[i] *= exp10(r)
end

open("testfloats", "w") do f
    for num in nums
        println(f, num)
    end
end

function strtod()
    for line in eachline("testfloats")
        parse(Float64, line)
    end
end

using Mmap
function xparse()
    io = IOBuffer(Mmap.mmap("testfloats"))
    for i = 1:100_000
        Main.Parsers.xparse(io, Float64)
        Main.Parsers.readbyte(io)
    end
end

function bench(n=100_000)
    open("results", "w") do f
        println(f, "exp,float,strtod,xparse")
        rng = -322:307
        for i = 1:n
            exp = rand(rng)
            r = string(exp10(exp) * round(rand(), digits=rand(1:9)))
            strtod = @elapsed begin
                a1 = parse(Float64, r)
            end
            io = IOBuffer(r)
            xparse = @elapsed begin
                a2 = Main.Parsers.xparse(io, Float64)
            end
            if a1 != a2.result
                @warn "a1: $a1, a2: $a2"
            end
            println(f, exp, ',', r, ',', strtod, ',', xparse)
        end
    end
end

function prof(str, n)
    io = IOBuffer(str)
    res = Parsers.Result(Float64)
    for i = 1:n
        seekstart(io)
        r = Parsers.defaultparser(io, res)
    end
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, HolidayCalendars

makedocs(
    modules = [HolidayCalendars],
    clean = false,
    format = Documenter.HTML(),
    sitename = "HolidayCalendars.jl",
    authors = "Iain Skett",
    pages = [
        "Holiday Calendars" => "index.md",
    ],
)

deploydocs(
    repo = "github.com/InfiniteChai/HolidayCalendars.jl.git",
    devbranch = "main",
    devurl = "latest"
)


text:
using OctaveREPL
using Test

@testset "OctaveREPL.jl" begin
    # Write your tests here.
end


text:
using BioSequences, StatsBase, DataFrames, FASTX
import ..enzyme_list
using ..wgregseq: promoter_finder.Promoter_Calculator

function join_seqs(up::String, down::String)
    return up * down
end

function join_seqs(up::BioSequences.LongDNA, down::BioSequences.LongDNA)
    return up * down
end

function add_sites_oligo(seq::String, up::String, down::String)
    return
end

function add_sites_oligo(seq::BioSequences.LongDNA, up::BioSequences.LongDNA, down::BioSequences.LongDNA)
    return join_seqs(join_seqs(up, seq), down)
end

function add_re_sites(oligo::BioSequences.LongDNA, re1::String, re2::String)
    re1_site = LongDNA{4}(enzyme_list[enzyme_list.enzyme.==re1, "site"][1])
    re2_site = LongDNA{4}(enzyme_list[enzyme_list.enzyme.==re2, "site"][1])
    return add_sites_oligo(oligo, re1_site, re2_site)
end

function import_primer(index::Int, direction::String)
    dir = @__DIR__
    path = joinpath(split(dir, '/')[1:end-1])
    if direction == "fwd"
        record = open(FASTA.Reader, "/$path/data/forward_finalprimers.fasta") do r
            collect(r)[index]
        end
        return sequence(record)
    elseif direction == "rev"
        record = open(FASTA.Reader, "/$path/data/reverse_finalprimers.fasta") do r
            collect(r)[index]
        end
        return sequence(record) |> reverse_complement
    else
        throw(ArgumentError("dir has to be either \"fwd\" or \"rev\""))
    end

end

function add_primer(oligos::Vector{LongSequence{DNAAlphabet{4}}}, index::Int, direction::String="both")
    if direction == "fwd"
        fwd_primer = import_primer(index, "fwd")
        rev_primer = LongDNA{4}("")
    elseif direction == "rev"
        rev_primer = import_primer(index, "rev")
        fwd_primer = LongDNA{4}("")
    elseif direction == "both"
        fwd_primer = import_primer(index, "fwd")
        rev_primer = import_primer(index, "rev")
    else
        throw(ArgumentError("direction has to be either \"fwd\", \"rev\" or \"both\"."))
    end
    return add_sites_oligo.(oligos, [fwd_primer], [rev_primer])
end

"""
    find_seq(TSS, strand, up, dn, genome)

Find sequence upstream and downstream of TSS. If gene is reverse
transcribed, return reverse complement of sequence.

Parameters
----------
TSS : int
    Transcription start site. Refers to position in the genome.
strand : string
    Either '+' for fwd strand or '-' for reverse strand.
up : int
    Number of basepairs upstream of TSS to grab. Relative to TSS according to strand / direction.
dn : int
    Number of basepairs downstream of TSS to grab. Relative to TSS according to strand / direction.
genome : string
    Reference genome as a string

Returns
-------
output : string
    Genetic sequence upstream and downstream of TSS
"""
function find_seq(TSS::Int, strand::String, up::Int, dn::Int, genome::BioSequences.LongDNA)

    if strand == "-"
        gene = genome[TSS-dn:TSS+up-1]

        outgene = LongDNA{4}(gene) |> reverse_complement

        left_pos = TSS-dn
        right_pos = TSS+up

    elseif strand == "+"
        outgene = genome[TSS-up+1:TSS+dn]
        left_pos = TSS-up
        right_pos = TSS+dn
    else
        throw(ArgumentError("`strand` has to be either \"+\" or \"-\""))
    end
    return (outgene, left_pos, right_pos)
end

function _random_mutation_generator(sequence, rate)::Array{Tuple{Int, Int}, 1}
    num_mutations = Int(floor(length(sequence) * rate))
    positions = sample(1:length(sequence), num_mutations, replace=false)
    mutants = sample(1:3, num_mutations)
    return [(x, y) for (x, y) in zip(positions, mutants)]
end

function random_mutation_generator(sequence, rate, num_mutants)
    mutant_list = Vector{Vector{Tuple{Int, Int}}}([])
    for i in 1:num_mutants
        push!(mutant_list, _random_mutation_generator(sequence, rate))
    end
    return mutant_list
end

"""
    function mutate_from_index(sequence, index, alph)

Generate mutated sequence from wild type sequence given a an index and
an alphabet.

Parameters
----------
- sequence : string
    Wild type seqence
- index : array of tuples
    Each tuple is a position and the mutation at that position.
- alph : array
    Alphabet to pick mutation from

Returns
-------
- str
    Sequence including mutation
"""
function mutate_from_index(sequence, index; alphabet=[DNA_A, DNA_C, DNA_G, DNA_T])
    mut_seq = deepcopy(sequence)
    for (locus, mutation) in index
        mut_seq[locus] = filter(x-> x != mut_seq[locus], alphabet)[mutation]
    end
    return mut_seq
end

"""

Creates single or double mutants.

#Parameters
----------
- sequence : string
    DNA sequence that is going to be mutated.
- num_mutants : `Int`, default None
    Number of mutant sequences. If None, all possible mutatants are created.
- mut_per_seq : `Int`, default 1
    Number of mutations per sequence.
- site_start : `int``, default 0
    Beginning of the site that is about to be mutated.
- site_end : int, default -1
    End of the site that is about to be mutated.
- alph_type : string, default "DNA"
    Can either be "DNA" for letter sequences, or "Numeric" for integer sequences.
- number_fixed : bool
    If True, the number of mutations is fixed as the rate times length of the sequence.
- keep_wildtype : bool, default False
    If True, adds wild type sequence as first sequence in the list.

#Returns
-------
- mutants : list
    List of mutant sequences. Each element is a string.
"""
function mutations_rand(
    sequence::BioSequences.LongDNA{4},
    rate::Float64,
    num_mutants::Int;
    site_start=1,
    site_end=length(sequence)
)

    mutation_window = sequence[site_start:site_end]

    # Create list
    mutants = BioSequences.LongDNA{4}[]
    push!(mutants, sequence)

    mutant_indices = random_mutation_generator(mutation_window, rate, num_mutants)

    for x in mutant_indices
        push!(mutants, sequence[1:site_start-1] * mutate_from_index(mutation_window, x) * sequence[site_end+1:end])
    end
    return mutants
end

"""
    find_restriction_sites(enzyme, sequence_list)

Searches for restriction sites of a specific enzyme in a list of sequences

#Parameters
----------
enzyme : string
    Name of the enzyme.
sequence_list: array-type
    list of suquences

"""
function find_restriction_sites(enzyme::String, sequence_list::Vector{LongSequence{DNAAlphabet{4}}})

    ind = findfirst(x -> x == enzyme, enzyme_list.enzyme)
    site = LongDNA{4}(enzyme_list[ind, :site])
    return sum(occursin.((ExactSearchQuery(site),), sequence_list))
end

"""
    find_restriction_sites(enzyme, sequence_list)

Searches for restriction sites of a specific enzyme in a list of sequences

#Parameters
----------
enzyme : string
    Name of the enzyme.
sequence_list: array-type
    list of suquences

"""
function find_restriction_sites(enzyme::Vector{String}, sequence_list::Vector{LongSequence{DNAAlphabet{4}}})
    sites = find_restriction_sites.(enzyme, (sequence_list,))
    return DataFrame(enzyme=enzyme, sites=sites)

end

"""
    function find_best_promoter(df)

Find the promoter with the lowest predicited free energy, using the model from
La Fleur et al, 2022.

#Parameters
-----------
df : DataFrame
    DataFrame containing the transcription start sites for promoters.
wt_sequence : BioSequences.LongDNA
    Sequence of wild type genome.

#Returns
df[ind, :] : Dataframe Row
    Row of DataFrame containing the strongest predicted site.

"""
function find_best_promoter(df::DataFrames.DataFrame, wt_sequence::BioSequences.LongDNA)
    if "tss" ∉ names(df)
        throw(ArgumentError("DataFrame has no column \"tss\"."))
    end
    if "direction" ∉ names(df)
        throw(ArgumentError("DataFrame has no column \"tss\"."))
    end
    # Find range of possible start sites
    min_tss = df.tss |> minimum
    max_tss = df.tss |> maximum

    # Get sequence around start sites
    sequence = wt_sequence[Int(min_tss)-115:Int(max_tss)+115]

    # Run model on sequence
    p = Promoter_Calculator()
    r = p(sequence)
    if df.direction[1] == "-"
        sites = [r["Reverse_Predictions_per_TSS"][x - (Int(min_tss)-115)] for x in df.tss]
    elseif df.direction[1] == "+"
        sites = [r["Forward_Predictions_per_TSS"][x - (Int(min_tss)-115)] for x in df.tss]
    else
        throw(ArgumentError("Direction has to be either \"+\" or \"-\"."))
    end
    # Find strongest site
    ind = argmin([site["dG_total"] for site in sites])
    return df[ind, :]
end

"""
    check_primers_re_sites(enz1, enz2, primer, direction)

Check if restriction enzymes occur in primer site.
"""
function check_primers_re_sites(enz1, enz2, primer, direction)
    site1 = enzyme_list[enzyme_list.enzyme .== enz1, "site"][1] |> LongDNA{4}
    site2 = enzyme_list[enzyme_list.enzyme .== enz2, "site"][1] |> LongDNA{4}

    clear = false
    if direction == "both"
        while clear == false
            fwd_primer = import_primer(primer, "fwd")
            rev_primer = import_primer(primer, "rev")
            if occursin(ExactSearchQuery(site1), fwd_primer)
                primer += 1
            elseif occursin(ExactSearchQuery(site1), rev_primer)
                primer += 1
            elseif occursin(ExactSearchQuery(site2), fwd_primer)
                primer += 1
            elseif occursin(ExactSearchQuery(site2), rev_primer)
                primer += 1
            else
                clear = true
            end
        end
    elseif direction in ["fwd", "rev"]
        while clear == false
            _primer = import_primer(primer, direction)
            if occursin(ExactSearchQuery(site1), _primer)
                primer += 1
            elseif occursin(ExactSearchQuery(site2), _primer)
                primer += 1
            else
                clear = true
            end
        end
    else
        throw(ArgumentError("direction has to be in [\"both\", \"fwd\", \"rev\"]."))
    end
    return primer
end


text:
## Umbrella functions to read and write graph file formats ##

"""
    readgraph(filename[; format, gzip, ...])

Load a network from file.

Currently supported formats are:

* `:auto` — try to gess the format from the file's extenstion,
* `:delim` — delimited plain-text file containing one edge per line
  (typically, `*.csv`, `*.tsv`, and `*.txt` extension),
* `:gml` — Graph Modeling Language file (`*.gml`),
* `:pajek` — Pajek file (`*.paj`, `*.net`),
* `:ngraph` — ngraph binary format (`*.bin` + `*.json`).

In the works:

* `:ucinet`  — UCINet files (`*.ucn`),
* `:graphml` — GraphML files (`*.graphml`),
* `:gdf` — GUESS file (https://gephi.org/users/supported-graph-formats/gdf-format/)

See the excellent [Gephi manual](https://gephi.org/users/supported-graph-formats/)
for the overview of different formats.
"""
function readgraph(filename; format=:auto, gzip=:auto, kvargs...)
    extensions = split(lowercase(basename(filename)), '.')[2:end]
    matchext(exts...) = any(i -> i in extensions, exts)
    if format == :auto
        if matchext("txt", "tsv", "csv", "edgelist")
            format = :delim
        elseif matchext("net", "paj")
            format = :pajek
        elseif matchext("gml")
            format = :gml
        elseif matchext("dl")
            format = :ucinet
        else
            error("could not determine file format automatically")
        end
    end
    if gzip == :auto
        gzip = matchext("gz")
    end
    f = gzip ? GZip.open(filename) : open(filename)
    if format == :delim
        g = read_edgelist(f; kvargs...)
    elseif format == :pajek
        g = read_pajek(f; kvargs...)
    elseif format == :gml
        g = read_gml(f; kvargs...)
    elseif format == :ucinet
        g = read_dl(f)
    end
    close(f)
    return g
end


text:
"""
    CircularBuffer{T}(n)

The CircularBuffer type implements a circular buffer of fixed capacity
where new items are pushed to the back of the list, overwriting values
in a circular fashion.

Allocate a buffer of elements of type `T` with maximum capacity `n`.
"""
mutable struct CircularBuffer{T} <: AbstractVector{T}
    capacity::Int
    first::Int
    length::Int
    buffer::Vector{T}

    CircularBuffer{T}(capacity::Int) where {T} = new{T}(capacity, 1, 0, Vector{T}(undef, capacity))
end

CircularBuffer(capacity) = CircularBuffer{Any}(capacity)

Base.IndexStyle(::Type{<:CircularBuffer}) = IndexLinear()

"""
    empty!(cb::CircularBuffer)

Reset the buffer.
"""
function Base.empty!(cb::CircularBuffer)
    cb.length = 0
    return cb
end

@inline function _buffer_index(cb::CircularBuffer, i::Int)
    @boundscheck checkbounds(cb, i)
    n = capacity(cb)
    idx = cb.first + i - 1
    return idx > n ? idx - n : idx
end

"""
    cb[i]

Get the i-th element of CircularBuffer.

* `cb[1]` to get the element at the front
* `cb[end]` to get the element at the back
"""
Base.@propagate_inbounds function Base.getindex(cb::CircularBuffer, i::Int)
    j = _buffer_index(cb, i)
    @inbounds return cb.buffer[j]
end

"""
    cb[i] = data

Store data to the `i`-th element of `CircularBuffer`.
"""
Base.@propagate_inbounds function Base.setindex!(cb::CircularBuffer, data, i::Int)
    j = _buffer_index(cb, i)
    @inbounds cb.buffer[j] = data
    return cb
end

"""
    pop!(cb::CircularBuffer)

Remove the element at the back.
"""
function Base.pop!(cb::CircularBuffer)
    # To be consistent with `Base.pop!`, always check for emptiness and throw an
    # `ArgumentError` if empty. We use `checkbounds` here (and in `popfirst!`) instead of
    # `isempty` because benchmarking shows substantially (~30%) improved performance on some
    # machines.
    checkbounds(Bool, cb, 1) || throw(ArgumentError("circular buffer must be non-empty"))
    @inbounds out = last(cb)
    cb.length -= 1
    return out
end

"""
    push!(cb::CircularBuffer, data)

Add an element to the back and overwrite front if full.
"""
function Base.push!(cb::CircularBuffer, data)
    # if full, increment and overwrite, otherwise push
    if isfull(cb)
        cb.first = (cb.first == capacity(cb) ? 1 : cb.first + 1)
    else
        cb.length += 1
    end
    @inbounds cb[length(cb)] = data
    return cb
end

"""
    popfirst!(cb::CircularBuffer)

Remove the element from the front of the `CircularBuffer`.
"""
function popfirst!(cb::CircularBuffer)
    checkbounds(Bool, cb, 1) || throw(ArgumentError("circular buffer must be non-empty"))
    @inbounds out = first(cb)
    cb.first = (cb.first == capacity(cb) ? 1 : cb.first + 1)
    cb.length -= 1
    return out
end

"""
    pushfirst!(cb::CircularBuffer, data)

Insert one or more items at the beginning of CircularBuffer
and overwrite back if full.
"""
function pushfirst!(cb::CircularBuffer, data)
    # if full, decrement and overwrite, otherwise pushfirst
    cb.first = (cb.first == 1 ? capacity(cb) : cb.first - 1)
    if !isfull(cb)
        cb.length += 1
    end
    @inbounds cb.buffer[cb.first] = data
    return cb
end

"""
    append!(cb::CircularBuffer, datavec::AbstractVector)

Push at most last `capacity` items.
"""
function Base.append!(cb::CircularBuffer, datavec::AbstractVector)
    # push at most last `capacity` items
    n = length(datavec)
    for i in max(1, n-capacity(cb)+1):n
        push!(cb, datavec[i])
    end
    return cb
end

"""
    fill!(cb::CircularBuffer, data)

Grows the buffer up-to capacity, and fills it entirely.
It doesn't overwrite existing elements.
"""
function Base.fill!(cb::CircularBuffer, data)
    for i in 1:capacity(cb)-length(cb)
        push!(cb, data)
    end
    return cb
end

"""
    length(cb::CircularBuffer)

Return the number of elements currently in the buffer.
"""
Base.length

"""
    size(cb::CircularBuffer)

Return a tuple with the size of the buffer.
"""
Base.size(cb::CircularBuffer) = (cb.length,)

"""
    isempty(cb::CircularBuffer)

Test whether the buffer is empty.
"""
Base.isempty

""""
    capacity(cb::CircularBuffer)

Return capacity of CircularBuffer.
"""
capacity(cb::CircularBuffer) = cb.capacity

"""
    isfull(cb::CircularBuffer)

Test whether the buffer is full.
"""
isfull(cb::CircularBuffer) = length(cb) == capacity(cb)

"""
    first(cb::CircularBuffer)

Get the first element of CircularBuffer.
"""
@inline function Base.first(cb::CircularBuffer)
    @boundscheck checkbounds(cb, 1)
    @inbounds return cb.buffer[cb.first]
end

"""
    last(cb::CircularBuffer)

Get the last element of CircularBuffer.
"""
Base.@propagate_inbounds Base.last(cb::CircularBuffer) = cb[end]


text:
# Import needed libraries
using Plots
using XLSX
using DataFrames
using Statistics

# Read the relevant excel workbook
df = DataFrames.DataFrame(XLSX.readtable("Folds5x2_pp.xlsx", "Sheet5")...)

# Split into design matrix and target vector
X = df[:, 1:4]
design_matrix = convert(Matrix, X)

y = df[:, 5]  # Target vector

# Split into training and test
train_size = 0.80
data_size = size(design_matrix)[1]

train_index = trunc(Int, train_size * data_size)

# Split using the desired train size
X_train = design_matrix[1:train_index, :]
X_test = design_matrix[train_index+1:end, :]

y_train = y[1:train_index]
y_test = y[train_index+1:end]

"""
    scale_features(X)

This function attempts to normalise the design matrix (X) user pass.
The input data X is standardised and return along with other learned metrics.

A tuple with 3 elements is returned representing (Standardised data, mean, std deviation).

"""
function scale_features(X)

    μ = mean(X, dims=1)
    σ = std(X, dims=1)

    X_norm = (X .- μ) ./ σ

    return (X_norm, μ, σ)
end

"""
    transform_features(X, μ, σ)

This functions uses the mean and standard deviation values users pass to
normalise a new design matrix.
"""
function transform_features(X, μ, σ)
    X_norm = (X .- μ) ./ σ
    return X_norm
end

# Scale training features and get artificats for future use
X_train_scaled, μ, σ = scale_features(X_train)

# Transform the testing features by using the learned artifacts
X_test_scaled = transform_features(X_test, μ, σ)

"""
    mean_squared_cost(X, y, θ)

This function computes the batch cost based on the values of the design matrix (X),
target vector (y), and the weights (θ) passed to it.
"""
function mean_squared_cost(X, y, θ)
    # Sample size
    m = size(X)[1]

    # Vectorised Prediction loss
    preds = X * θ
    loss = preds - y

    # Half mean squared loss
    cost =  (1/(2m)) * (loss' * loss)

    return cost
end

"""
    lin_reg_grad_descent(X, y, α, fit_intercept=true, n_iter=1000)

This function uses gradient descent algorithm to find the best weights (θ)
that minimises the mean squared loss between the predictions that the model
generates and the target vector (y).

A tuple of 1D vectors representing the weights (θ)
and a history of loss at each iteration (𝐉) is returned.
"""
function lin_reg_grad_descent(X, y, α, fit_intercept=true, n_iter=2000)
    # Initialize some useful values
    m = length(y) # number of training examples

    if fit_intercept
        # Add a constant of 1s if fit_intercept is specified
        constant = ones(m, 1)
        X = hcat(constant, X)
    else
        X # Assume user added constants
    end

    # Use the number of features to initialise the theta θ vector
    n = size(X)[2]
    θ = zeros(n)

    # Initialise the cost vector based on the number of iterations
    𝐉 = zeros(n_iter)

    for iter in range(1, stop=n_iter)
        pred = X * θ

        # Calcaluate the cost for each iter
        𝐉[iter] = mean_squared_cost(X, y, θ)

        # Update the theta θ at each iter
        θ = θ - ((α/m) * X') * (pred - y);
    end
    return (θ, 𝐉)
end

θ, 𝐉 = lin_reg_grad_descent(X_train_scaled, y_train, 0.05, true, 3000)

plot(𝐉,
     label="Cost per iter",
     ylabel="Cost",
     xlabel="Number of Iteration",
     title="Cost Per Iteration")

savefig("cost_plot")

"""
    predict(X, θ, fit_intercept=true)

This function uses the learned weights (θ) to make predictions based on the
design matrix passed as (X).

The 1D vector representing these predictions is finally returned.
"""
function predict(X, θ, fit_intercept=true)
    m = size(X)[1]

    if fit_intercept
        constant = ones(m)
        X = hcat(constant, X)
    else
        X
    end

    predictions = X * θ

    return predictions
end

# Make predictions for both training and testing datasets
train_preds = predict(X_train_scaled, θ)
test_preds = predict(X_test_scaled, θ)

"""
    rmse_score(y_true, y_pred)

This function calculates the RMSE based on the values of ground truth (y_true)
and the predicted values generated by a model (y_pred)
"""
function rmse_score(y_true, y_pred)

    errors = y_pred - y_true
    errors² = errors .^ 2
    mse = mean(errors²)
    rmse = sqrt(mse)

    return rmse

end

println("RMSE for Training Set: ", rmse_score(y_train, train_preds))
println("RMSE for Testing Set: ", rmse_score(y_test, test_preds))

"""
    r_squared_score(y_pred, y_true)

This function returns the R squared value based on the predictions (y_pred)
and the ground truth values (y_true)passed to it.
"""
function r_squared_score(y_pred, y_true)
    # Just a convinient way of using notations
    ∑ = sum
    μ = mean

    # Compute sum of explained variance (SST) and sum of squares of residuals
    sst = ∑(((y_true .- μ(y_true)) .^ 2))
    ssr = ∑(((y_pred .- y_true) .^ 2))

    r² = 1 - (ssr / sst)

    return r²
end

# Get the r-squared score for training and test datasets
train_r² = r_squared_score(train_preds, y_train)
test_r² = r_squared_score(test_preds, y_test)

println("Training R² score for test sets: ", train_r²)
println("Testing R² score for test sets: ", test_r²)

max_stars_repo_name: stratosmacker/Python-for-Data-Science


text:
using Base.Test, JSON

const JSON_FILE = tempname()
println("Path to JSON_FILE: $JSON_FILE")

Lumberjack.add_truck(JsonTruck(open(JSON_FILE, "w")))

# test without dates
log("debug", "some-msg")
log("info", "some-msg")
log("warn", "some-msg")
log("error", "some-msg")
log("crazy", "some-msg")

# test with msec_date_saw
add_saw(Lumberjack.msec_date_saw)

log("debug", "some-msg")
log("info", "some-msg")
log("warn", "some-msg")
log("error", "some-msg")
log("crazy", "some-msg")

remove_saws()

# test with fn_call_saw
add_saw(Lumberjack.fn_call_saw)

let
    @noinline caller(mode, msg) = log(mode, msg)
    caller("debug", "some-msg")
    caller("info", "some-msg")

    remove_saws()

    # test with stacktrace_saw
    add_saw(Lumberjack.stacktrace_saw)

    @noinline child_caller(mode, msg) = log(mode, msg)
    @noinline parent_caller(mode, msg) = child_caller(mode, msg)
    parent_caller("warn", "some-msg")
    parent_caller("error", "some-msg")

    remove_saws()

    # test with extra params
    log("debug", "some-msg", Dict{Any,Any}( :thing1 => "thing1" ))
    log("info", "some-msg", Dict{Any,Any}( :thing1 => "thing1", :thing2 => 69 ))
    log("warn", "some-msg", Dict{Any,Any}( :thing1 => "thing1", :thing2 => 69, :thing3 => [1, 2, 3] ))
    log("error", "some-msg", Dict{Any,Any}( :thing1 => "thing1", :thing2 => 69, :thing3 => [1, 2, 3], :thing4 => Dict{Any,Any}( "a" => "apple" )))
    log("crazy", "some-msg", Dict{Any,Any}( :thing1 => "thing1", :thing2 => 69, :thing3 => [1, 2, 3], :thing4 => Dict{Any,Any}( "a" => "apple" ), :thing5 => :some_symbol ))
end

js = JSON.parse("[$(join(readlines(open(JSON_FILE)), ','))]")

for j in js
    @test j["msg"] == "some-msg"
end

for i = 6:10
    @test haskey(js[i], "date")
end

for i = 11:12
    @test js[i]["lookup"]["name"] == "caller"
end

for i = 13:14
    @test js[i]["stacktrace"][1]["name"] == "child_caller"
    @test js[i]["stacktrace"][1]["file"] == basename(string(@__FILE__))
    @test js[i]["stacktrace"][2]["name"] == "parent_caller"
    @test js[i]["stacktrace"][2]["file"] == basename(string(@__FILE__))
end

for i = 15:19
    @test js[i]["thing1"] == "thing1"
end

for i = 16:19
    @test js[i]["thing2"] == 69
end

for i = 17:19
    @test js[i]["thing3"] == [1, 2, 3]
end

for i = 18:19
    @test js[i]["thing4"]["a"] == "apple"
end

@test js[19]["thing5"] == "some_symbol"

# clean up
@test success(`rm $JSON_FILE`)


text:
using GraphPlot
using Plots
pyplot()
using PyPlot
pygui(false)
include("loadbuild.jl")
include("graphing.jl")
include("beliefprop.jl")
include("Psis.jl")

EDGEFILE = "data/save.feather"
BIASFILE = "data/bias.csv"
plotdir = "results/plots/"
figext = "eps"

problem_types = ["bias", "fake"]
epsilon = 0.44
MAXITER = 2

Plots.PyPlot.rc("lines", lw=10)
# Plots.PyPlot.rc("font", family="monospace")
# Plots.PyPlot.rc("font", weight="bold")
Plots.PyPlot.rc("font", size=22)

if EDGEFILE == "data/save.feather"
  by_domain_flag = true
elseif EDGEFILE == "data/save2.feather"
  by_domain_flag = false
end

save = loadSave(EDGEFILE)

# size(unique(save[:sdom]))

#println(size(unique(vcat(save[:rsrc],save[:rdest]))))

bias = loadBias(BIASFILE)

cats = buildWeb(save, by_domain_flag) #only mutual links for "a"

edgeunion = countunion(cats["a"], cats["img"])
@time for x in ["link", "script"]
    countunion!(edgeunion, cats[x])
end

sets = [("union", 1)] #[("union",10)], [("mut", 1000), ("mut", 50), ("mut", 2), ("mut", 1), ("a", 1000), ("a", 500),("a", 100), ("img", 200), ("img", 50)]

#s2p = Dict{String, Array}()

roc_x = Vector{Float64}()
roc_y = Vector{Float64}()

k = 3 #number of folds

for (sn, set) in enumerate(sets)
   edgetype = set[1]
    if edgetype == "union"
      edgeset = edgeunion
    else
      edgeset = cats[set[1]]
    end

    # if ~by_domain_flag
    #   domainmap = Dict{String, String}()
    #   for j in edgeset
    #     domainmap[j.first] = split(j.first, '/')[1]
    #     for i in j.second
    #       domainmap[i.first] = split(i.first, '/')[1]
    #     end
    #   end
    # end

    G, s2i, i2s = constructGraph(edgeset, minsamp=set[2], giant=false, by_domain=by_domain_flag)
    #G2, s2i2, i2s2 = constructGraph(edgeset, minsamp=set[2], giant=true, by_domain=by_domain_flag)
    mindeg = minimum([length([x for x in all_neighbors(G, v) if x != v]) for v in 1:nv(G)])
    println("Minimum Degree: $mindeg")
    maxdeg = maximum(degree(G))
    println("Maximum Degree: $maxdeg")
    # lx, ly = spring_layout(G)
    #lx2, ly2 = spring_layout(G2)

    for key in problem_types
        p = Plots.plot()
        plot!(p, [0, 1], [0, 1], line=:dash, label="")
        println(key*string(set))
        try mkdir(plotdir*string(set[1])) catch x end
        try mkdir(plotdir*string(set[1])*"/"*string(set[2])) catch x end
        mydir = plotdir*string(set[1])*"/"*string(set[2])*"/"

        # plotGraph(bias, G, s2i, i2s, lx, ly, color=key,
        #             path=mydir*key*"train.png")
        # ϕ = makeBeliefs(bias, G, s2i, i2s, key)
        # pr, b, lodds = beliefprop(G, ϕ, Psis(0.44), 1);
        # plotGraph(bias, G, s2i, i2s, lx, ly, color=b[:,1],
                    # path=mydir*key*"prop.png")

        # plotGraph(bias, G2, s2i2, i2s2, lx2, ly2, color=key,
        #             path=mydir*key*"train_main.png")
        # ϕ = makeBeliefs(bias, G2, s2i2, i2s2, key)
        # pr, b, lodds = beliefprop(G2, ϕ, Psis(0.44), 1);
        # plotGraph(bias, G2, s2i2, i2s2, lx2, ly2, color=b[:,1],
        #             path=mydir*key*"prop_main.png")

        folds = makeFolds(G, k) # makeFolds(G, k, i2s, domainmap)
        aucs = []
        for f in 1:k
          println("Fold " *string(f))
            ϕ = makeBeliefs(bias, G, s2i, i2s, key, folds=folds, fold=f)
            pr, b, lodds = beliefprop(G, ϕ, Psis(epsilon), MAXITER);
            # for v in 1:nv(G)
            #     if folds[v] == f
            #         site = i2s[v]
            #         if ~haskey(s2p, site)
            #             s2p[site] = ones(1, size(sets, 1)+2)./2.0
            #             s2p[site][size(sets, 1)+1] = reality(bias, site, "bias")
            #             s2p[site][size(sets, 1)+2] = folds[v]
            #         end
            #         s2p[site][sn] = b[v]
            #     end
            # end

            roc_x, roc_y, acc, lvls, cms = getROC(bias, G, s2i, i2s, key, b, folds=folds, fold=f)
            auc = AUC(roc_x, roc_y)
            #println("Max "*string(maximum(acc))*" at "*string(lvls[find(x->x==maximum(acc), acc)[1]]))
            println("AUC : "*string(auc))
            best_cutoff = collect(keys(acc))[indmax(collect(values(acc)))]
            println("best cutoff threshold : "*string(best_cutoff))
            println("confusion matrix (best threshold) : ")
            show(cms[best_cutoff])
            println("")
            println("confusion matrix (0.5 threshold) : ")
            show(cms[0.5])
            println("")
            #println(typeof(roc_x))
            #println(typeof(roc_y))
            #println(roc_x)
            #println(roc_y)
            plot!(p, roc_x, roc_y,label="Fold $f (AUC=$(@sprintf("%.3f", auc)))")
            push!(aucs, auc)
        end
        avg_auc = sum(aucs)/k
        title!("Receiver Operating Characteristic (Average AUC=$(@sprintf("%.3f", avg_auc)))")
        println("Avg AUC : "*string(avg_auc))
        figpath = mydir*key*"_ROC_curve.$figext"
        println("Saving figure to: $figpath")
        savefig(figpath)

        open(mydir*"score.txt", "w") do f
            write(f, string(sum(aucs)/k))
        end
    end
end

# for site in keys(s2p)
#     s2p[site][10] = reality(bias, site, "bias")
# end

# s2p

# s2p["cnn.com"]
# s2p["bluebirdbanter.com"]

# using JSON
# open("metadata.json", "w") do f
#     write(f, JSON.json(s2p))
# end

max_stars_repo_name: jpfairbanks/NewsGraphing

max_issues_repo_name: jpfairbanks/NewsGraphing

max_forks_repo_name: jpfairbanks/NewsGraphing


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct CVPipe <: Pipe
        name::String
        arc::Arc
        available::Bool
        diameter::Float64
        length::Float64
        roughness::Float64
        flowlimits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        flow::Union{Nothing,Float64}
        headloss::Union{Nothing,Float64}
        internal::InfrastructureSystemsInternal
    end

A unidirectional-flow pipe that contains a check valve (CV).

# Arguments
- `name::String`
- `arc::Arc`
- `available::Bool`
- `diameter::Float64`
- `length::Float64`
- `roughness::Float64`
- `flowlimits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `flow::Union{Nothing,Float64}`
- `headloss::Union{Nothing,Float64}`
- `internal::InfrastructureSystemsInternal`: internal reference, do not modify
"""
mutable struct CVPipe <: Pipe
    name::String
    arc::Arc
    available::Bool
    diameter::Float64
    length::Float64
    roughness::Float64
    flowlimits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    flow::Union{Nothing,Float64}
    headloss::Union{Nothing,Float64}
    "internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function CVPipe(name, arc, available, diameter, length, roughness, flowlimits, flow, headloss, )
    CVPipe(name, arc, available, diameter, length, roughness, flowlimits, flow, headloss, InfrastructureSystemsInternal(), )
end

function CVPipe(; name, arc, available, diameter, length, roughness, flowlimits, flow, headloss, )
    CVPipe(name, arc, available, diameter, length, roughness, flowlimits, flow, headloss, )
end

# Constructor for demo purposes; non-functional.
function CVPipe(::Nothing)
    CVPipe(;
        name="init",
        arc=Arc(nothing),
        available=true,
        diameter=0.0,
        length=0.0,
        roughness=0.0,
        flowlimits=(min=0.0, max=10.0),
        flow=nothing,
        headloss=nothing,
    )
end

"""Get CVPipe name."""
get_name(value::CVPipe) = value.name
"""Get CVPipe arc."""
get_arc(value::CVPipe) = value.arc
"""Get CVPipe available."""
get_available(value::CVPipe) = value.available
"""Get CVPipe diameter."""
get_diameter(value::CVPipe) = value.diameter
"""Get CVPipe length."""
get_length(value::CVPipe) = value.length
"""Get CVPipe roughness."""
get_roughness(value::CVPipe) = value.roughness
"""Get CVPipe flowlimits."""
get_flowlimits(value::CVPipe) = value.flowlimits
"""Get CVPipe flow."""
get_flow(value::CVPipe) = value.flow
"""Get CVPipe headloss."""
get_headloss(value::CVPipe) = value.headloss
"""Get CVPipe internal."""
get_internal(value::CVPipe) = value.internal


text:
using JLD2, Distributions, Plots, StatPlots, AverageShiftedHistograms, DynNets
 plotly()

N_est = 900
Tvals =  [50 500 5000]
N_ind = length(Tvals)[1]

N = 10 ;
GBAdgp =1;
GBAest =1;   useStartVal = false

useStartVal? startValFlag = "T" : startValFlag = "F";# false #false #
# save
save_fold = "./data/estimatesTest/"
    file_name = "EstTest_N_$(N)_T_$(Tvals)_$(N_est)_GBAdgp_$(GBAdgp)_GBAest_$(GBAest)_$(useStartVal)_target.jld"
    save_path = save_fold*file_name#
estDataT = @load(save_path,groupsIndsEst,ObsAllT,ParTvDgpAllT,simPar_1,simPar_2,simPar_3,
        estPar_1,estPar_2,estPar_3,est_conv_flag,est_times)

##

RavgMSE(A,B) = sqrt.(mean((A .- B).^2,1))
AvgBias(A,B) = mean((A .- B),1)
tmpModelEst = DynNets.SdErgm1(zeros(10,10),[ones(N),ones(GBAest),ones(GBAest)],groupsIndsEst)
tmpModelDgp = DynNets.SdErgm1(zeros(10,10),[ones(N),ones(GBAdgp),ones(GBAdgp)],groupsIndsEst)

rmseTotGas = zeros(N_est,N,N_ind)
biasTotGas = zeros(N_est,N,N_ind)
rmseTotSnap = zeros(N_est,N,N_ind)
biasTotSnap = zeros(N_est,N,N_ind)
dgpTvParT = Array{Array{Float64,3},1}(N_ind)
filTvParT = Array{Array{Float64,3},1}(N_ind)
ParSnapT =  Array{Array{Float64,3},1}(N_ind)

#Remove Nans
if sum(isnan(estPar_3))/(prod(size(estPar_2)))<0.01
    println("NaN Fraction = $(sum(isnan(estPar_3))/(prod(size(estPar_2))))")
    for indT = 1:length(Tvals)
        #kind of bootstrap
        indNan  = find(isnan , estPar_3[indT,:,1])
        indReal = find(isfinite , estPar_3[indT,:,1])
        randInd = rand(indReal,length(indNan))
        estPar_1[indT,indNan,:] = estPar_1[indT,randInd,:]
        estPar_2[indT,indNan,1] = estPar_2[indT,randInd,1]
        estPar_3[indT,indNan,1] = estPar_3[indT,randInd,1]
    end
else
    #too many nans
    error()
end
# Compute MSE Bias and Snap estimates
@time for indT = 1:N_ind
#    indT=1
            T =Tvals[indT]
            tmpParDgp = zeros(T,N,N_est)
            tmpParFil = zeros(T,N,N_est)
            tmpParSnap = zeros(T,N,N_est)
        for n=1:N_est
#    n=100
    obsT = ObsAllT[indT][:,:,n]
    vecParEst =  [estPar_1[indT,n,:]; estPar_2[indT,n];  estPar_3[indT,n]]
    vecParDgp =  [simPar_1[indT,n,:]; simPar_2[indT,n];  simPar_3[indT,n]]

    GroupIndsEst  = groupsIndsEst ;
    tmpParFil[:,:,n] =   DynNets.score_driven_filter_or_dgpAndLikeliood(tmpModelEst,vecParEst,
                                                degsT = obsT,
                                                groupsInds = GroupIndsEst )[1]

    fitSnap = DynNets.estSnapSeq(DynNets.SdErgm1(obsT)) # sequence of snapshots estimate
    # remove the infinites
    for i=1:5
        infInd = find(.!isfinite.(fitSnap))

        replInd = infInd[infInd .> 2]
        fitSnap[replInd] = fitSnap[replInd.-1]
        infInd = find(.!isfinite.(fitSnap))
        replInd = infInd[infInd .< (Int(N*T) -2) ]
        fitSnap[replInd] = fitSnap[replInd.+1]
    end
    tmpParSnap[:,:,n] = fitSnap

    GroupIndsDgp =  groupsIndsEst
    tmpParDgp[:,:,n] = DynNets.score_driven_filter_or_dgpAndLikeliood(tmpModelDgp,vecParDgp,
                                                degsT = obsT,
                                                groupsInds = GroupIndsDgp )[1]

            rmseTotGas[n,:,indT] = RavgMSE(tmpParFil[:,:,n],tmpParDgp[:,:,n] )
            biasTotGas[n,:,indT] = AvgBias(tmpParFil[:,:,n],tmpParDgp[:,:,n] )
            rmseTotSnap[n,:,indT] = RavgMSE(fitSnap,tmpParDgp[:,:,n] )
            biasTotSnap[n,:,indT] = AvgBias(fitSnap,tmpParDgp[:,:,n] )
        end
    dgpTvParT[indT]  = tmpParDgp
    filTvParT[indT]  = tmpParFil
    ParSnapT[indT] = tmpParSnap
end

## Example single realization plots
Nlines=2
indPl = round(Int,linspace(1,N,Nlines))
C(g::ColorGradient) = RGB[g[z] for z=linspace(0,1,Nlines)]
lineColors = C(cgrad(:rainbow))

indT = 3
n=1
T = Tvals[indT]

fit_T = dgpTvParT[indT][:,:,n]
fil_fit_T =   filTvParT[indT][:,:,n]
fitSnap =   ParSnapT[indT][:,:,n]
 p=plot()
for i = 1:Nlines p = plot!(fit_T[:,indPl[i]],  color = lineColors[i]  ) end #lineColors
for i = 1:Nlines p = plot!(fil_fit_T[:,indPl[i]],linestyle=:dash,  color = lineColors[i]  ) end #lineColors
for i = 1:Nlines p= plot!(fitSnap[:,indPl[i]],linestyle=:dot,markershape=:circle,
    markersize=1,linewidth = 0, markerstrokecolor=  lineColors[i], color = lineColors[i]  ) end

p =plot(p,titlefont = font(16),xlims = (0,T) ,titleloc=:center,legend=:none,
    title = "Gas DGP,N $(N) T = $(T) B = $(simPar_2[indT,n,1]), A = $(simPar_3[indT,n,1])",size=(1350,600))

## Box Plots
yLim = (0,1)
    pRT1 = boxplot( rmseTotGas[:,:,1]  , ylims = yLim  , legend = false, title = "T = $(Tvals[1]), N = $(N)")
    pRT2 = boxplot( rmseTotGas[:,:,2]  , ylims = yLim  , legend = false, title = "T = $(Tvals[2])")
     pRT3 = boxplot( rmseTotGas[:,:,3]  , ylims = yLim , legend = false, title = "T = $(Tvals[3])")
     yLim = (-0.8,1.1)
    pBT1 = boxplot( biasTotGas[:,:,1] , ylims = yLim , legend = false, title = "T = $(Tvals[1])")
    pBT2 = boxplot( biasTotGas[:,:,2] , ylims = yLim , legend = false, title = "T = $(Tvals[2])")
    pBT3 = boxplot( biasTotGas[:,:,3] , ylims = yLim ,legend = false, title = "T = $(Tvals[3])")
    ptot = plot(pRT1,pRT2,pRT3,pBT1,pBT2,pBT3,titlefont = font(12) ,titleloc=:center,size=(1350,600) )
#
# yLim = (0,1.5)
#     pRT1 = boxplot( rmseTotSnap[:,:,1]  , ylims = yLim  , legend = false, title = "T = $(Tvals[1]), N = $(N), GWdgp =  $(GWdgp), Single Snapshot ")
#     pRT2 = boxplot( rmseTotSnap[:,:,2]  , ylims = yLim  , legend = false, title = "T = $(Tvals[2])")
#      pRT3 = boxplot( rmseTotSnap[:,:,3]  , ylims = yLim , legend = false, title = "T = $(Tvals[3])")
#      yLim = (-1.4,0.8)
#     pBT1 = boxplot( biasTotSnap[:,:,1] , ylims = yLim , legend = false, title = "T = $(Tvals[1])")
#     pBT2 = boxplot( biasTotSnap[:,:,2] , ylims = yLim , legend = false, title = "T = $(Tvals[2])")
#     pBT3 = boxplot( biasTotSnap[:,:,3] , ylims = yLim ,legend = false, title = "T = $(Tvals[3])")
#     ptot = plot(pRT1,pRT2,pRT3,pBT1,pBT2,pBT3,titlefont = font(12) ,titleloc=:center,size=(1350,600) )

#    title = ["($i)" for j = 1:1, i=1:11], titleloc = :right, titlefont = font(8)

## histograms of static parameters

function Bias(A::Array{<:Real,2},B::Array{<:Real,2})
    return  Bias =    (A.-B)./std(A .-B)
end
Bias(A::SharedArray{<:Real,2},B::SharedArray{<:Real,2}) = BiasRel( convert(Array{Real},A),convert(Array{Real},B))

for ind = 1:length(Tvals)

#ind=3
errA = Bias(simPar_3[ind,:,:],estPar_3[ind,:,:])
errB = Bias(simPar_2[ind,:,:],estPar_2[ind,:,:])
errW =  (simPar_1[ind,:,:].- estPar_1[ind,:,:])./(repmat(std(estPar_1[ind,:,:],1),N_est,1))
plot_tot =Vector{Plots.Plot{Plots.PlotlyBackend}}
rng_space = -5:.1:5

err= errW
tmpAsh = ash(err;rng =rng_space )
pW1 = plot(tmpAsh,
    title ="  T= $(Tvals[ind]),
    W    , mean(std) =  $(round(mean(std(err,1)),3)) ")
plot!(Normal(0,1),linewidth = 5)

err= errB
tmpAsh = ash(err;rng = rng_space)
pB1 = plot(tmpAsh,
title =    " $(N_est) Est   B   mu =  $(round(mean(err[:]),3)) ")
plot!(Normal(0,1),linewidth = 5)

err= errA

tmpAsh = ash(err;rng = rng_space)
pA1 = plot(tmpAsh,
title = startValFlag * "    $(GBAest) / $(N) A  mu =  $(round(mean(err[:]),3)), std =  $(round(std(err[:]),3)) ")
plot!(Normal(0,1),linewidth = 5)

exp = "pl_T_$(ind)=  plot(pW1,pB1,pA1,layout=(3,1) )"
eval(parse(exp))

end
std(err[:])

exp2 = " " ; for i = 1:ind exp2 *=  "pl_T_$(i) , " end

exp3 = "plot( " * exp2 * " layout = (1,ind),size=(1350,600),legend =:none)"
eval(parse(exp3))
#plot(pl_T_1 , pl_T_2,pl_T_3 , pl_T_3, layout = (1,4),size=(1350,600),legend =:none)

##

plot(errW)

##


text:
push!(LOAD_PATH,joinpath(@__DIR__, "../src/"))
using LinearAlgebra
using StableApproxEPW

# ## Target of the approximation problem

# We consider the Helmholtz solution in the unit disk, with wavenumber
k = 5;
# We need to define the maximum Fourier mode number `P` in the approximation
# target:
P = 25;
# The target approximation will be of the form
# ```math
# u = \mathbf{x} ↦ \sum_{|p| \leq P} u_p b_{p}(\mathbf{x}).
# ```
# Next we construct the vector of coefficients in the basis `b_p`
# for `p` in `[-P,P]`:
U = zeros(ComplexF64, 2P+1);
U[P+1]     = 0.5;  # This is the constant mode (`p=0`)
U[P+1 + P] = 1im;  # This is mode `p = P`
# Here the only non-zero coefficients are ``u_0 = \frac{1}{2}`` and
# ``u_P = \imath`` so that
# ```math
# u = \frac{1}{2} b_{0} + \imath b_{P}.
# ```
# The target of the approximation problem can then be constructed as:
u = solution_surrogate(U; k=k);
# ``u`` can be evaluated at any `(r,θ)` point.
# Alternatively, the target ``u`` could have been a single mode.
# For instance, to get the circular wave with mode number `p=15`, simply set
# ``u = bp(15; k=k)``.
# Of course, any other function (defined in polar coordinates) can be defined
# by the user as target of the approximation problem.

# ## Reconstruction method

# The approximation will be reconstructed by sampling the target on the
# boundary of the unit disk.
# To do so, we need to know now the dimension `N` of the approximation sets
# that we are going to use:
N = 100;
# We can determine the number of sampling nodes necessary for a successful
# reconstruction, based on `N`, `P` (to avoid aliasing) and the oversampling
# ratio `η`:
S = number_of_boundary_sampling_nodes(N; η=2, P=P);
# The (equispaced) boundary nodes are then constructed as:
X = boundary_sampling_nodes(S);
# The right-hand-side of the linear system can then be readily constructed:
b = samples_from_nodes(u, X);

# ## Approximation with **propagative** plane waves

# We construct the approximation set of PPW:
Φppw = approximation_set(N; k=k);
# The matrix and its (SVD) factorization
Appw = samples_from_nodes(Φppw, X);
iAppw = RegularizedSVDPseudoInverse(Appw; ϵ=1e-14);
# The coefficients of the approximation are computed:
ξppw = solve_via_regularizedSVD(iAppw, b);
ũppw = (r,θ) -> sum([ξi * ϕi(r,θ) for (ξi, ϕi) in zip(ξppw, Φppw)]);
# Absolute error function
eppw = (r,θ) -> ũppw(r,θ) - u(r,θ); eppw(1,π/2)
# Let's compute the relative residual:
resppw = norm(Appw * ξppw - b) / norm(b)
# We did not obtained any accuracy whatsoever!
# The reason is that the coefficients are too large:
nrmppw = norm(ξppw) / norm(U)

# ## Approximation with **evanescent** plane waves

# We construct the approximation set of EPW:
Φepw = approximation_set(N, P, sobol_sampling; k=k);
# It is possible to choose other types of sampling methods, instead of
# `sobol_sampling`, for instance `uniform_sampling` and `random_sampling`.
# The matrix and its (SVD) factorization
Aepw = samples_from_nodes(Φepw, X);
iAepw = RegularizedSVDPseudoInverse(Aepw; ϵ=1e-14);
# The coefficients of the approximation are computed:
ξepw = solve_via_regularizedSVD(iAepw, b);
ũepw = (r,θ) -> sum([ξi * ϕi(r,θ) for (ξi, ϕi) in zip(ξepw, Φepw)]);
# Absolute error function
eepw = (r,θ) -> ũepw(r,θ) - u(r,θ); eepw(1,π/2)
# Let's compute the relative residual:
resepw = norm(Aepw * ξepw - b) / norm(b)
# We get almost 8 digits of accuracy!
# The size of the coefficients remains quite high:
nrmepw = norm(ξepw) / norm(U)

# ### Down to machine precision

# Let's double the number of EPW
N = 200
# The above approximation process can be obtained much more rapidly with the
# following convenience function
_, ξepw, ũepw, resepw, nrmepw, eepw = Dirichlet_sampling(k, U, N; smpl_type=sobol_sampling);
# One can check that we obtain now a relative residual very close to machine
# precision (13 digits of accuracy):
resepw
# The size of the coefficients is also greatly reduced:
nrmepw

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
const DEFAULT_REGISTRIES = Pkg.RegistrySpec[Pkg.RegistrySpec(;
    name="General",
    uuid="23338594-aafe-5451-b93e-139f81909106",
    url="https://github.com/JuliaRegistries/General.git",
)]

"""
    main(
        env::AbstractDict=ENV,
        ci_cfg::CIService=auto_detect_ci_service(; env=env);
        entry_type::EntryType=KeepEntry(),
        registries::Vector{Pkg.RegistrySpec}=DEFAULT_REGISTRIES,
        subdirs::AbstractVector{<:AbstractString}=[""],
        master_branch::Union{DefaultBranch,AbstractString}=DefaultBranch(),
        bump_compat_containing_equality_specifier=true,
        pr_title_prefix::String="",
        include_jll::Bool=false,
        unsub_from_prs=false,
        cc_user=false,
        bump_version=false,
    )

Main entry point for the package.

# Arguments
- `env::AbstractDict=ENV`: Optional dictionary of environment variables, see README for overview
- `ci_cfg::CIService=auto_detect_ci_service(; env=env)`: CI Configuration, default to what is auto-detected

# Keywords
- `entry_type::EntryType=KeepEntry()`: How to handle bumps for entry types
- `registries::Vector{Pkg.RegistrySpec}=DEFAULT_REGISTRIES`: RegistrySpec of all registries to use
- `subdirs::AbstractVector{<:AbstractString}=[""]`: Subdirectories for nested packages
- `master_branch::Union{DefaultBranch,AbstractString}=DefaultBranch()`: Name of the master branch
- `bump_compat_containing_equality_specifier=true`: Bump compat entries with equality specifiers
- `pr_title_prefix::String=""`: Prefix for pull request titles
- `include_jll::Bool=false`: Include JLL packages to bump
- `unsub_from_prs=false`: Unsubscribe the user from the pull requests
- `cc_user=false`: CC the user on the pull requests
- `bump_version=false`: When set to true, the version in Project.toml will be bumped if a pull request is made. Minor bump if >= 1.0, or patch bump if < 1.0
"""
function main(
    env::AbstractDict=ENV,
    ci_cfg::CIService=auto_detect_ci_service(; env=env);
    entry_type::EntryType=KeepEntry(),
    registries::Vector{Pkg.RegistrySpec}=DEFAULT_REGISTRIES,
    subdirs::AbstractVector{<:AbstractString}=[""],
    master_branch::Union{DefaultBranch,AbstractString}=DefaultBranch(),
    bump_compat_containing_equality_specifier=true,
    pr_title_prefix::String="",
    include_jll::Bool=false,
    unsub_from_prs=false,
    cc_user=false,
    bump_version=false,
)
    generated_prs = Vector{Union{GitHub.PullRequest,GitLab.MergeRequest}}()

    api, repo = get_api_and_repo(ci_cfg)

    for subdir in subdirs
        deps = get_project_deps(
            api,
            ci_cfg,
            repo;
            subdir=subdir,
            include_jll=include_jll,
            master_branch=master_branch,
        )
        get_latest_version_from_registries!(deps, registries)

        for dep in deps
            pr = @mock make_pr_for_new_version(
                api,
                repo,
                dep,
                entry_type,
                ci_cfg;
                subdir=subdir,
                master_branch=master_branch,
                env=env,
                bump_compat_containing_equality_specifier=bump_compat_containing_equality_specifier,
                pr_title_prefix=pr_title_prefix,
                unsub_from_prs=unsub_from_prs,
                cc_user=cc_user,
                bump_version=bump_version,
            )

            if !isnothing(pr)
                push!(generated_prs, pr)
            end
        end
    end

    return generated_prs
end