text:
using Plots, Test

@testset "Shapes" begin
    @testset "Type" begin
        square = Shape([(0,0.0),(1,0.0),(1,1.0),(0,1.0)])
        @test isa(square, Shape{Int64, Float64})
    end

    @testset "Copy" begin
        square = Shape([(0,0),(1,0),(1,1),(0,1)])
        square2 = Shape(square)
        @test square2.x == square.x
        @test square2.y == square.y
    end

    @testset "Center" begin
        square = Shape([(0,0),(1,0),(1,1),(0,1)])
        @test Plots.center(square) == (0.5,0.5)
    end

    @testset "Translate" begin
        square = Shape([(0,0),(1,0),(1,1),(0,1)])
        squareUp = Shape([(0,1),(1,1),(1,2),(0,2)])
        squareUpRight = Shape([(1,1),(2,1),(2,2),(1,2)])

        @test Plots.translate(square,0,1).x == squareUp.x
        @test Plots.translate(square,0,1).y == squareUp.y

        @test Plots.center(translate!(square,1)) == (1.5,1.5)
    end

    @testset "Rotate" begin
        # 2 radians rotation matrix
        R2 = [cos(2) sin(2); -sin(2) cos(2)]
        coords = [0 0; 1 0; 1 1; 0 1]'
        coordsRotated2 = R2*(coords.-0.5).+0.5

        square = Shape([(0,0),(1,0),(1,1),(0,1)])

        # make a new, rotated square
        square2 = Plots.rotate(square, -2)
        @test square2.x ≈ coordsRotated2[1,:]
        @test square2.y ≈ coordsRotated2[2,:]

        # unrotate the new square in place
        rotate!(square2, 2)
        @test square2.x ≈ coords[1,:]
        @test square2.y ≈ coords[2,:]
    end
end

@testset "Brush" begin
    @testset "Colors" begin
        baseline = brush(1, RGB(0,0,0))
        @test brush(:black) == baseline
        @test brush("black") == baseline
    end
    @testset "Weight" begin
        @test brush(10).size == 10
        @test brush(0.1).size == 1
    end
    @testset "Alpha" begin
        @test brush(0.4).alpha == 0.4
        @test brush(20).alpha == nothing
    end
    @testset "Bad Argument" begin
        # using test_logs because test_warn seems to not work anymore
        @test_logs (:warn,"Unused brush arg: nothing (Nothing)") begin
            brush(nothing)
        end
    end
end

@testset "Fonts" begin
    @testset "Scaling" begin
        sizesToCheck = [:titlefontsize, :legendfontsize, :legendtitlefontsize,
            :xtickfontsize, :ytickfontsize, :ztickfontsize,
            :xguidefontsize, :yguidefontsize, :zguidefontsize,]
        # get inital font sizes
        initialSizes = [Plots.default(s) for s in sizesToCheck ]

        #scale up font sizes
        scalefontsizes(2)

        # get inital font sizes
        doubledSizes = [Plots.default(s) for s in sizesToCheck ]

        @test doubledSizes == initialSizes*2

        # reset font sizes
        resetfontsizes()

        finalSizes = [Plots.default(s) for s in sizesToCheck ]

        @test finalSizes == initialSizes
    end
end

@testset "Series Annotations" begin
    square = Shape([(0,0),(1,0),(1,1),(0,1)])
    @test_logs (:warn,"Unused SeriesAnnotations arg: triangle (Symbol)") begin
        p = plot([1,2,3],
            series_annotations=(["a"],
                                 2,    # pass a scale factor
                                 (1,4), # pass two scale factors (overwrites first one)
                                 square, # pass a shape
                                 font(:courier), # pass an annotation font
                                 :triangle  # pass an incorrect argument
                            ))
        sa = p.series_list[1].plotattributes[:series_annotations]
        @test sa.strs == ["a"]
        @test sa.font.family == "courier"
        @test sa.baseshape == square
        @test sa.scalefactor == (1,4)
    end
end


text:
# Common Subexpression elimination

# Original version together with Shashi Gowda:
# https://github.com/JuliaSymbolics/SymbolicUtils.jl/issues/121

## Works by keeping track of a dict of expressions that have already been seen
## Each new expression is assigned a fresh symbolic name

## Remove powers

# "Replace powers by repeated multiplies"
# function remove_powers(ex)

#     rule = @rule((~y)^2 => (~y) * (~y))
#     rule2 = @rule((~y)^(~n) => (~y) * (~y)^(~n - 1))

#     ex2 = Fixpoint(Postwalk(Chain([rule, rule2])))(ex)

#     return ex2
# end

# "Common subexpression elimination"
# function cse(expr)

#     # expr = remove_powers(expr)

#     dict = OrderedDict()

#     r = @rule ~x::(x -> x isa Term) =>
#         haskey(dict, ~x) ? dict[~x] : dict[~x] = make_symbol()

#     final = Postwalk(Chain([r]))(expr)

#     return [var => ex for (ex, var) in pairs(dict)], final
# end

cse(vars, ex::Num) = cse(vars, value(ex))
cse(vars, ex::Num, dict::OrderedDict) = cse(vars, value(ex), dict)

"Common subexpression elimination.
This version starts with the given variables"
# function cse(vars, expr)

#     dict = OrderedDict{Any, Any}(v => v for v in vars)

#     r = @rule ~x::(x -> x isa Term) =>
#         haskey(dict, ~x) ? dict[~x] : dict[~x] = make_symbol()

#     final = Postwalk(Chain([r]))(expr)
#     [var => ex for (ex, var) in pairs(dict)], final
# end
3
"CSE with pre-existing dict"
function cse(vars, expr, dict::OrderedDict=initialize_dict(vars))

    orig_length = length(dict)

    rules = [
        @rule ~x::(x -> (istree(x) || x isa Sym)) =>
            haskey(dict, ~x) ? dict[~x] : dict[~x] = make_symbol()
    ]

    final = Postwalk(Chain(rules))(expr)  # *modifies* dict
    # @show final
    # @show dict

    if length(dict) > orig_length
        new_additions = dict.keys[orig_length+1:end]
    else
        new_additions = []
    end

    [dict[var] => var for var in new_additions], final
end

initialize_dict(vars) = OrderedDict{Any, Any}(v => v for v in value.(vars))

make_Assignment(p::Pair) = Assignment(p[1], p[2])

cse_total(vars, ex::Num) = cse_total(vars, value(ex))
cse_total(vars, ex::Num, dict::OrderedDict) = cse_total(vars, value(ex), dict)

"Perform CSE and return a vector of `Assignment`s"
function cse_total(vars, ex, dict=initialize_dict(vars))
    ex2, final = cse(vars, ex, dict)
    ex3 = make_Assignment.(ex2)

    # remove redefinitions of input variables:
    ex4 = [eq for eq in ex3 if !(any(x -> isequal(x, eq.lhs), vars))]

    return ex4, final
end

function cse_total(vars, eq::Assignment, dict::OrderedDict)
    eqns, final = cse_total(vars, rhs(eq), dict)

    # replace final Assignment with the correct variable
    # from the left-hand side of eq:

    # this currently "wastes" a generated variable
    # eqns[end] = Assignment(lhs(eq), rhs(eqns[end]))

    # alternative: add a new Assignment:
    push!(eqns, Assignment(lhs(eq), final))

    return eqns
end

function cse_total(vars, eqs::Vector{Assignment})
    dict = initialize_dict(vars)

    all_eqns = []

    for eq in eqs
        push!(all_eqns, cse_total(vars, eq, dict))
        # modifies dict
        # @show dict
    end

    return reduce(vcat, all_eqns)
end

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Clingo_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Clingo")
JLLWrappers.@generate_main_file("Clingo", UUID("e1d4ad0e-dc8b-5d2d-8155-e9a55b531f08"))
end  # module Clingo_jll


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: outposts
using AWS.Compat
using AWS.UUIDs

"""
    create_order(line_items, outpost_identifier, payment_option)
    create_order(line_items, outpost_identifier, payment_option, params::Dict{String,<:Any})

Creates an order for an Outpost.

# Arguments
- `line_items`: The line items that make up the order.
- `outpost_identifier`:  The ID or the Amazon Resource Name (ARN) of the Outpost.
- `payment_option`: The payment option for the order.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"PaymentTerm"`: The payment terms for the order.
"""
function create_order(
    LineItems,
    OutpostIdentifier,
    PaymentOption;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return outposts(
        "POST",
        "/orders",
        Dict{String,Any}(
            "LineItems" => LineItems,
            "OutpostIdentifier" => OutpostIdentifier,
            "PaymentOption" => PaymentOption,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_order(
    LineItems,
    OutpostIdentifier,
    PaymentOption,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return outposts(
        "POST",
        "/orders",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "LineItems" => LineItems,
                    "OutpostIdentifier" => OutpostIdentifier,
                    "PaymentOption" => PaymentOption,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_outpost(name, site_id)
    create_outpost(name, site_id, params::Dict{String,<:Any})

Creates an Outpost. You can specify AvailabilityZone or AvailabilityZoneId.

# Arguments
- `name`:
- `site_id`:

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AvailabilityZone"`:
- `"AvailabilityZoneId"`:
- `"Description"`:
- `"Tags"`: The tags to apply to the Outpost.
"""
function create_outpost(Name, SiteId; aws_config::AbstractAWSConfig=global_aws_config())
    return outposts(
        "POST",
        "/outposts",
        Dict{String,Any}("Name" => Name, "SiteId" => SiteId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_outpost(
    Name,
    SiteId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return outposts(
        "POST",
        "/outposts",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("Name" => Name, "SiteId" => SiteId), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_outpost(outpost_id)
    delete_outpost(outpost_id, params::Dict{String,<:Any})

Deletes the Outpost.

# Arguments
- `outpost_id`:  The ID of the Outpost.

"""
function delete_outpost(OutpostId; aws_config::AbstractAWSConfig=global_aws_config())
    return outposts(
        "DELETE",
        "/outposts/$(OutpostId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_outpost(
    OutpostId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return outposts(
        "DELETE",
        "/outposts/$(OutpostId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_site(site_id)
    delete_site(site_id, params::Dict{String,<:Any})

Deletes the site.

# Arguments
- `site_id`:

"""
function delete_site(SiteId; aws_config::AbstractAWSConfig=global_aws_config())
    return outposts(
        "DELETE", "/sites/$(SiteId)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function delete_site(
    SiteId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return outposts(
        "DELETE",
        "/sites/$(SiteId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_outpost(outpost_id)
    get_outpost(outpost_id, params::Dict{String,<:Any})

Gets information about the specified Outpost.

# Arguments
- `outpost_id`:  The ID of the Outpost.

"""
function get_outpost(OutpostId; aws_config::AbstractAWSConfig=global_aws_config())
    return outposts(
        "GET",
        "/outposts/$(OutpostId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_outpost(
    OutpostId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return outposts(
        "GET",
        "/outposts/$(OutpostId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_outpost_instance_types(outpost_id)
    get_outpost_instance_types(outpost_id, params::Dict{String,<:Any})

Lists the instance types for the specified Outpost.

# Arguments
- `outpost_id`:  The ID of the Outpost.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:
- `"NextToken"`:
"""
function get_outpost_instance_types(
    OutpostId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return outposts(
        "GET",
        "/outposts/$(OutpostId)/instanceTypes";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_outpost_instance_types(
    OutpostId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return outposts(
        "GET",
        "/outposts/$(OutpostId)/instanceTypes",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_outposts()
    list_outposts(params::Dict{String,<:Any})

Create a list of the Outposts for your AWS account. Add filters to your request to return a
more specific list of results. Use filters to match an Outpost lifecycle status,
Availibility Zone (us-east-1a), and AZ ID (use1-az1).  If you specify multiple filters, the
filters are joined with an AND, and the request returns only results that match all of the
specified filters.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AvailabilityZoneFilter"`:  A filter for the Availibility Zone (us-east-1a) of the
  Outpost.   Filter values are case sensitive. If you specify multiple values for a filter,
  the values are joined with an OR, and the request returns all results that match any of the
  specified values.
- `"AvailabilityZoneIdFilter"`:  A filter for the AZ IDs (use1-az1) of the Outpost.
  Filter values are case sensitive. If you specify multiple values for a filter, the values
  are joined with an OR, and the request returns all results that match any of the specified
  values.
- `"LifeCycleStatusFilter"`:  A filter for the lifecycle status of the Outpost.   Filter
  values are case sensitive. If you specify multiple values for a filter, the values are
  joined with an OR, and the request returns all results that match any of the specified
  values.
- `"MaxResults"`:
- `"NextToken"`:
"""
function list_outposts(; aws_config::AbstractAWSConfig=global_aws_config())
    return outposts(
        "GET", "/outposts"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_outposts(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return outposts(
        "GET", "/outposts", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    list_sites()
    list_sites(params::Dict{String,<:Any})

Lists the sites for the specified AWS account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:
- `"NextToken"`:
"""
function list_sites(; aws_config::AbstractAWSConfig=global_aws_config())
    return outposts("GET", "/sites"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
end
function list_sites(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return outposts(
        "GET", "/sites", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags for the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.

"""
function list_tags_for_resource(
    ResourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return outposts(
        "GET",
        "/tags/$(ResourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    ResourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return outposts(
        "GET",
        "/tags/$(ResourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds tags to the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.
- `tags`: The tags to add to the resource.

"""
function tag_resource(ResourceArn, Tags; aws_config::AbstractAWSConfig=global_aws_config())
    return outposts(
        "POST",
        "/tags/$(ResourceArn)",
        Dict{String,Any}("Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    ResourceArn,
    Tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return outposts(
        "POST",
        "/tags/$(ResourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Tags" => Tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes tags from the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.
- `tag_keys`: The tag keys.

"""
function untag_resource(
    ResourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return outposts(
        "DELETE",
        "/tags/$(ResourceArn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    ResourceArn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return outposts(
        "DELETE",
        "/tags/$(ResourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
function add_solver_xml!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
	# Add solution display
	if !anyints(m)
		m.xmlmodel = replace(m.xmlmodel, r"(?s)<wantint><![CDATA[yes]]></wantint>", "<wantint><![CDATA[no]]></wantint>")
		# Lets hack this here. Maybe if we get fancy and offer other categories we should break this out a bit more
		m.xmlmodel = replace(m.xmlmodel, r"(?s)<category>milp</category>", "<category>lp</category>")
	end
	# Add user options
	param_string = ""
	for key in keys(m.solver.params)
		param_string *= "$(key) $(m.solver.params[key])\n"
	end
	# Add user options
	m.xmlmodel = replace(m.xmlmodel, r"(?s)<param>.*</param>", "<param><![CDATA[$(param_string)]]></param>")
end

function solution_status(m, args...)
	for s in args
		if match(Regex("SOLUTION STATUS\\W+?:\\W+?$s"), m.last_results) != nothing
			return true
		end
	end
	return false
end

function parse_status!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
	if solution_status(m, "INTEGER_OPTIMAL", "OPTIMAL")
		m.status = OPTIMAL
	elseif solution_status(m, "PRIMAL_UNBOUNDED", "UNBOUNDED")
		m.status = UNBOUNDED
	elseif solution_status(m, "DUAL_INFEASIBLE")
		m.status = UNBNDORINF
	elseif solution_status(m, "DUAL_UNBOUNDED", "PRIMAL_INFEASIBLE", "INFEASIBLE")
		m.status = INFEASIBLE
	end
end

function parse_objective!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
	sci = match(r"PRIMAL\W+?OBJECTIVE\W+?:\W+?(-?[\d\.]+e[\+\-]\d+)", m.last_results).captures
	m.objVal = parse(Float64, sci[1])
end

function parse_solution!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
	for v in matchall(r"V(\d+).+?(-?[\d\.]+e[\+\-]\d+)", m.last_results)
		regmatch = match(r"V(\d+).+?(-?[\d\.]+e[\+\-]\d+)", v)
		m.solution[parse(Int, regmatch.captures[1])] = parse(Float64, regmatch.captures[2])
	end
end

function parse_duals!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
	m.duals = zeros(m.nrow)
	for v in matchall(r"\d+\W+C(\d+).+?\n", m.last_results)
		s = split(v)
		d_lb = parse(Float64, s[7])
		d_ub = -parse(Float64, s[8])
		m.duals[parse(Int, s[2][2:end])] = abs(d_lb) > abs(d_ub)?d_lb:d_ub
	end

	m.reducedcosts = zeros(m.ncol)
	for v in matchall(r"\d+\W+V(\d+).+?\n", m.last_results)
		s = split(v)
		d_lb = parse(Float64, s[7])
		d_ub = -parse(Float64, s[8])
		m.reducedcosts[parse(Int, s[2][2:end])] = abs(d_lb) > abs(d_ub)?d_lb:d_ub
	end

end


text:
function definitions(b::Block)
  defs = [Variable(i) for i = 1:length(b.ir.defs) if b.ir.defs[i][1] == b.id]
  append!(defs, arguments(b))
end

function usages(b::Block)
  uses = Set{Variable}()
  prewalk(b) do x
    x isa Variable && push!(uses, x)
    return x
  end
  return uses
end

function dominators(ir)
  doms = Dict(b => Set(blocks(ir)) for b in blocks(ir))
  worklist = blocks(ir)
  while !isempty(worklist)
    b = popfirst!(worklist)
    # We currently special case the first block here,
    # since Julia sometimes creates blocks with no predecessors,
    # which otherwise throw off the analysis.
    ds = isempty(predecessors(b)) ? Set([b, block(ir, 1)]) :
      push!(intersect([doms[c] for c in predecessors(b)]...), b)
    if ds != doms[b]
      doms[b] = ds
      for c in successors(b)
        c in worklist || push!(worklist, c)
      end
    end
  end
  return doms
end

function domtree(ir, start = 1)
  doms = dominators(ir)
  doms = Dict(b => filter(c -> b != c && b in doms[c], blocks(ir)) for b in blocks(ir))
  children(b) = filter(c -> !(c in union(map(c -> doms[c], doms[b])...)), doms[b])
  tree(b) = Pair{Int,Any}(b.id,tree.(children(b)))
  tree(block(ir, start))
end

function domorder(ir, start = 1; full = false)
  tree = domtree(ir, start)
  flatten((b,cs)) = vcat(b, flatten.(cs)...)
  tree = flatten(tree)
  if full
    for b = 1:length(ir.blocks)
      b in tree || push!(tree, b)
    end
  end
  return tree
end

domorder!(ir::IR, start = 1) = permute!(ir, domorder(ir, start, full = true))

function verify(ir::IR)
  @assert issorted(domorder(ir)) "Blocks are not in domtree order."
  doms = dominators(ir)
  # TODO check definitions within a block
  for (b, ds) in doms
    defs = union(definitions.(ds)...)
    for x in setdiff(usages(b), defs)
      error("Variable $x in block $(b.id) is not defined.")
    end
  end
  return
end

function renumber(ir)
  p = Pipe(ir)
  for (v, st) in p
    ex = st.expr
    if isbits(ex) # Trivial expressions can be inlined
      delete!(p, v)
      substitute!(p, v, substitute(p, ex))
    end
  end
  return finish(p)
end

function merge_returns!(ir)
  bs = [b for b in blocks(ir) if isreturn(b)]
  length(bs) == 1 && bs[1] == blocks(ir)[end] && return ir
  block!(ir)
  ret = argument!(blocks(ir)[end])
  return!(ir, ret)
  for b in bs
    branches(b)[end] = branch(length(ir.blocks), arguments(branches(b)[end])[1])
  end
  return ir
end

function expand!(ir::IR)
  worklist = blocks(ir)
  spats = Dict(b => Dict() for b in blocks(ir))
  while !isempty(worklist)
    b = pop!(worklist)
    b.id == 1 && continue
    defs = definitions(b)
    uses = usages(b)
    for v in setdiff(uses, defs)
      haskey(spats[b], v) && continue
      spats[b][v] = argument!(b, v)
      for c in predecessors(b)
        c in worklist || push!(worklist, c)
      end
    end
    ir.blocks[b.id] = prewalk(x -> get(spats[b], x, x), ir.blocks[b.id])
  end
  return ir
end

function prune!(ir::IR)
  worklist = blocks(ir)
  while !isempty(worklist)
    b = popfirst!(worklist)
    isempty(arguments(b)) && continue
    brs = filter(br -> br.block == b.id, [br for a in blocks(ir) for br in branches(a)])
    isempty(brs) && continue
    inputs = [setdiff(in, (a,)) for (a, in) in zip(arguments(b), zip(arguments.(brs)...))]
    del = findall(x -> length(x) == 1, inputs)
    rename = Dict(zip(arguments(b)[del], first.(inputs[del])))
    if !isempty(rename)
      prewalk!(x -> get(rename, x, x), ir)
      deletearg!(b, del)
      for c in successors(b)
        c in worklist || push!(worklist, c)
      end
    end
  end
  return ir
end

function ssa!(ir::IR)
  current = 1
  defs = Dict(b => Dict{Slot,Any}() for b in 1:length(ir.blocks))
  todo = Dict(b => Dict{Int,Vector{Slot}}() for b in 1:length(ir.blocks))
  function reaching(b, v)
    haskey(defs[b.id], v) && return defs[b.id][v]
    b.id == 1 && return undef
    x = defs[b.id][v] = argument!(b, insert = false)
    for pred in predecessors(b)
      if pred.id < current
        for br in branches(pred, b)
          push!(br.args, reaching(pred, v))
        end
      else
        push!(get!(todo[pred.id], b.id, Slot[]), v)
      end
    end
    return x
  end
  for b in blocks(ir)
    current = b.id
    rename(ex) = prewalk(x -> x isa Slot ? reaching(b, x) : x, ex)
    for (v, st) in b
      ex = st.expr
      if isexpr(ex, :(=)) && ex.args[1] isa Slot
        defs[b.id][ex.args[1]] = rename(ex.args[2])
        delete!(ir, v)
      else
        ir[v] = rename(ex)
      end
    end
    for i = 1:length(basicblock(b).branches)
      basicblock(b).branches[i] = rename(basicblock(b).branches[i])
    end
    for (succ, ss) in todo[b.id], br in branches(b, succ)
      append!(br.args, [reaching(b, v) for v in ss])
    end
  end
  return ir
end

function inlineable!(ir)
  pushfirst!(ir, Expr(:meta, :inline))
  return ir
end

function log!(ir, msg)
  pushfirst!(ir, xcall(Core, :println, msg))
  return ir
end

function pis!(ir::IR)
  for (v, st) in ir
    ex = st.expr
    ex isa PiNode || continue
    ir[v] = xcall(Core, :typeassert, ex.val, ex.typ)
  end
  return ir
end


text:
using SymbolicUtils, Test
using SymbolicUtils: polynormalize, Term, symtype
@testset "polyform" begin
    @syms a b c d
    @test polynormalize(a * (b + -1 * c) + -1 * (b * a + -1 * c * a)) == 0
    @eqtest polynormalize(sin(a+b)+sin(c+d)) == sin(a+b) + sin(c+d)
    @eqtest simplify(polynormalize(sin((a+b)^2)^2)) == simplify(sin(a^2+2*(b*a)+b^2)^2)
    @test simplify(polynormalize(sin((a+b)^2)^2 + cos((a+b)^2)^2)) == 1
    @syms x1::Real f(::Real)::Real

    # issue 193
    @test isequal(polynormalize(f(x1 + 2.0)), f(2.0 + x1))
    @test symtype(polynormalize(f(x1 + 2.0))) == Real

    # cleanup rules
    @test polynormalize(Term{Number}(identity, 0)) == 0
    @test polynormalize(Term{Number}(one, 0)) == 1
    @test polynormalize(Term{Number}(zero, 0)) == 0
    @test polynormalize(identity(a * b) - b * a) == 0
    @test polynormalize(a * b - b * a) == 0
end


text:
"""
    fit(::Type{Fvcb}, df; Tᵣ = nothing, VcMaxRef = 0., JMaxRef = 0., RdRef = 0., TPURef = 0.)

Optimize the parameters of the [`Fvcb`](@ref) model. Also works for [`FvcbIter`](@ref).

# Arguments

- df: a DataFrame with columns A, PPFD, Tₗ and Cᵢ, where each row is an observation. The column
names should match exactly
- Tᵣ: reference temperature for the optimized parameter values. If not provided, use the average Tₗ.
- VcMaxRef, JMaxRef, RdRef: initialisation values for the parameter optimisation

# Examples

```julia
using PlantBiophysics, Plots, DataFrames

file = joinpath(dirname(dirname(pathof(PlantBiophysics))),"test","inputs","data","P1F20129.csv")
df = read_walz(file)
# Removing the Rh and light curves for the fitting because temperature varies
filter!(x -> x.curve != "Rh Curve" && x.curve != "ligth Curve", df)

# Fit the parameter values:
VcMaxRef, JMaxRef, RdRef, TPURef = fit(Fvcb, df; Tᵣ = 25.0)
# Note that Tᵣ was set to 25 °C in our response curve. You should adapt its value to what you
# had during the response curves

# Checking the results:
filter!(x -> x.curve == "CO2 Curve", df)

# Sort the DataFrame by :Cᵢ to get ordered data point
sort!(df, :Cᵢ)

# Re-simulating A using the newly fitted parameters:
leaf =
    LeafModels(
        photosynthesis = FvcbRaw(VcMaxRef = VcMaxRef, JMaxRef = JMaxRef, RdRef = RdRef, TPURef = TPURef),
        Tₗ = df.Tₗ, PPFD = df.PPFD, Cᵢ = df.Cᵢ
    )
photosynthesis!(leaf)
df_sim = DataFrame(leaf)

# Visualising the results:
ACi_struct = PlantBiophysics.ACi(VcMaxRef, JMaxRef, RdRef, df.A, df_sim.A, df[:,:Cᵢ], df_sim.Cᵢ)
plot(ACi_struct,leg=:bottomright)

# Note that we can also simulate the results using the full photosynthesis model too (Fvcb):
# Adding the windspeed to simulate the boundary-layer conductance (we put a high value):
df[!, :Wind] .= 10.0

leaf = LeafModels(
        photosynthesis = Fvcb(VcMaxRef = VcMaxRef, JMaxRef = JMaxRef, RdRef = RdRef, Tᵣ = 25.0, TPURef = TPURef),
        # stomatal_conductance = ConstantGs(0.0, df[i,:gs]),
        stomatal_conductance = Medlyn(0.03, 12.),
        Tₗ = df.Tₗ, PPFD = df.PPFD, Cₛ = df.Cₐ, Dₗ = 0.1
    )

w = Weather(select(df, :T, :P, :Rh, :Cₐ, :T => (x -> 10) => :Wind))
photosynthesis!(leaf, w)
df_sim2 = DataFrame(leaf)

# And finally we plot the results:
ACi_struct_full = PlantBiophysics.ACi(VcMaxRef, JMaxRef, RdRef, df.A, df_sim2.A, df[:,:Cᵢ], df_sim2.Cᵢ)
plot(ACi_struct_full,leg=:bottomright)
# Note that the results differ a bit because there are more variables that are re-simulated (e.g. Cᵢ)
```
"""
function fit(::T, df; Tᵣ = nothing, VcMaxRef = 0.0, JMaxRef = 0.0, RdRef = 0.0, TPURef = 0.0) where {T<:Union{Type{Fvcb},Type{FvcbIter},Type{FvcbRaw}}}
    if Tᵣ === nothing
        Tᵣ = mean(df.Tₗ)
    end

    function model(x, p)
        leaf =
            LeafModels(
                photosynthesis = FvcbRaw(VcMaxRef = p[1], JMaxRef = p[2], RdRef = p[3], TPURef = p[4]),
                Tₗ = x[:, 1], PPFD = x[:, 2], Cᵢ = x[:, 3]
            )
        photosynthesis!(leaf)
        DataFrame(leaf).A
    end

    # Fitting the A-Cᵢ curve using LsqFit.jl
    # fits = curve_fit(model, df.Cᵢ[ind], df.A[ind], [VcMaxRef, JMaxRef, RdRef, TPURef])
    fits = curve_fit(model, Array(select(df, :Tₗ, :PPFD, :Cᵢ)), df.A, [VcMaxRef, JMaxRef, RdRef, TPURef])

    return (VcMaxRef = fits.param[1], JMaxRef = fits.param[2], RdRef = fits.param[3], TPURef = fits.param[4])
end

# Plot recipes for making A/Ci curves:
mutable struct ACi
    VcMaxRef
    JMaxRef
    RdRef
    A_meas
    A_sim
    Cᵢ_meas
    Cᵢ_sim
end

ACi(VcMaxRef, JMaxRef, RdRef, A_meas, A_sim, Cᵢ_meas) = ACi(VcMaxRef, JMaxRef, RdRef, A_meas, A_sim, Cᵢ_meas, copy(Cᵢ_meas))

@recipe function f(h::ACi)
    x = h.Cᵢ_meas
    x2 = h.Cᵢ_sim
    y = h.A_meas
    y2 = h.A_sim
    # Main plot (measurement):
    xguide --> "Cᵢ (ppm)"
    yguide --> "A (μmol m⁻² s⁻¹)"

    EF_ = round(EF(y, y2), digits = 3)
    dr_ = round(dr(y, y2), digits = 3)
    RMSE_ = round(RMSE(y, y2), digits = 3)

    @series begin
        seriestype := :scatter
        label := "Measured"
        x, y
    end

    @series begin
        label := "Simulated (EF:$EF_,dr:$dr_,RMSE:$RMSE_)"
        seriestype := :line
        x2, y2
    end
end


text:
"""
    InSet(x::AbstractIntVar, s::IntSetVar, trailer::Trailer)

InSet constraint, states that x ∈ s
"""
struct InSet <: Constraint
    x       ::AbstractIntVar
    s       ::IntSetVar
    active  ::StateObject{Bool}
    function InSet(x::AbstractIntVar, s::IntSetVar, trailer)
        constraint = new(x, s, StateObject{Bool}(true, trailer))
        addOnDomainChange!(x, constraint)
        addOnDomainChange!(s, constraint)
        return constraint
    end
end

"""
    propagate!(constraint::InSet, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`InSet` propagation function.
"""
function propagate!(constraint::InSet, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    # Feasibility

    if all(!is_possible(constraint.s.domain, v) for v in constraint.x.domain)
        return false
    end

    # Filter x
    x_values = collect(constraint.x.domain)
    for v in x_values
        if !is_possible(constraint.s.domain, v)
            remove!(constraint.x.domain, v)
            addToPrunedDomains!(prunedDomains, constraint.x, [v])
            triggerDomainChange!(toPropagate, constraint.x)
        end
    end

    # Filter s
    if isbound(constraint.x)
        if !is_required(constraint.s.domain, assignedValue(constraint.x))
            require!(constraint.s.domain, assignedValue(constraint.x))
            addToPrunedDomains!(prunedDomains, constraint.s, SetModification(;required=[assignedValue(constraint.x)]))
            triggerDomainChange!(toPropagate, constraint.s)
        end
    end

    # Deactivation
    if all(is_required(constraint.s.domain, v) for v in constraint.x.domain)
        setValue!(constraint.active, false)
    end

    return true
end

variablesArray(constraint::InSet) = [constraint.x, constraint.s]

function Base.show(io::IO, ::MIME"text/plain", con::InSet)
    println(io, typeof(con), ": ", con.x.id, " ∈ ", con.s.id, ", active = ", con.active)
    println(io, "   ", con.x)
    print(io, "   ", con.s)
end

function Base.show(io::IO, con::InSet)
    print(io, typeof(con), ": ", con.x.id, " ∈ ", con.s.id)
end


text:
using EventSimulation
using Distributions
using Printf

# flow of fax states
# WAITING ---> ENTRY ---> MOVED ---> SPECIAL
#                |                      |
#                |                      v
#                +-----------------> FINISHED
# WAITING   : fax waiting in entry queue
# ENTRY     : fax processed in entry service
# MOVED     : fax waiting in special queue
# SPECIAL   : fax processed in special service
# FINISHED  : fax left the system
@enum STATE WAITING ENTRY MOVED SPECIAL FINISHED

mutable struct Fax
    tin::Float64  # time when customer entered the system
    tout::Float64 # time when customer left the system
    special::Bool # did customer require special queue?
    state::STATE  # current state of the customer
end

mutable struct ServiceCenter <: AbstractState
    ia::Dict{Float64,Exponential{Float64}} # inter arrival distribution
    entry::Uniform{Float64}   # distribution of entry processing time in hours
    special::Uniform{Float64} # distribution of specia processing time in hours
    special_prob::Float64     # probability that fax requires special clerk
    entry_queue::SimQueue{Fax}   # entry queue
    special_queue::SimQueue{Fax} # special queue
    faxes::Vector{Fax}        # archive of all processed faxes

    function ServiceCenter()
        raw_ar = Dict(8.0 => 4.37, 9.0 => 6.24, 10.0 => 5.29, 11.0 => 2.97,
                      12.0 => 2.03, 13.0 => 2.79, 14.0 => 2.36, 15.0 => 1.04)
        ia = Dict(k => Exponential(1/v/60) for (k,v) in raw_ar)
        entry = Uniform(0.5/60, 4.5/60)
        special = Uniform(2.0/60, 6.0/60)
        special_prob = 0.2
        new(ia, entry, special, special_prob,
            SimQueue{Fax}(), SimQueue{Fax}(), Vector{Fax}())
    end
end

# new fax arrivals
function arrival(s)
    s.now > 16.0 && return # stopped when we close the shop!
    fax = Fax(s.now, NaN, false, WAITING) # new fax arrived
    push!(s.state.faxes, fax) # archive it
    provide!(s, s.state.entry_queue, fax) # add to entry queue
    # and schedule arrival of next fax
    register!(s, arrival, rand(s.state.ia[floor(s.now)]))
end

# at 12:00 a new shift is planned
function change_shift(s, ec2, sc2)
    # all morning clerks leave
    # but what is in process at 12:00 is finished by morning clerks
    empty!(s.state.entry_queue.requests)
    empty!(s.state.special_queue.requests)
    # afternoon clerks are scheduled
    # they have to finish processing all faxes even after 16:00
    for i in 1:ec2
        request!(s, s.state.entry_queue, (x,y) -> serve_entry(x, y, Inf))
    end
    for i in 1:sc2
        request!(s, s.state.special_queue, (x,y) -> serve_special(x, y, Inf))
    end
end

# start service of entry clerk that ends the shift at endtime
function serve_entry(s, fax, endtime)
    fax.state = ENTRY
    register!(s, x -> finish_entry(x, fax, endtime), rand(s.state.entry))
end

# finish service of entry clerk that ends the shift at endtime
function finish_entry(s, fax, endtime)
    # was the fax hard - if yes move it to special queue
    if rand() < s.state.special_prob
        fax.special = true
        fax.state = MOVED
        provide!(s, s.state.special_queue, fax)
    else
        fax.state = FINISHED
        fax.tout = s.now
    end

    # if shift is not finished request another fax
    if s.now < endtime
        request!(s, s.state.entry_queue, (x,y) -> serve_entry(x, y, endtime))
    end
end

# start service of special clerk that ends the shift at endtime
function serve_special(s, fax, endtime)
    fax.state = SPECIAL
    register!(s, x -> finish_special(x, fax, endtime), rand(s.state.special))
end

# finish service of special clerk that ends the shift at endtime
function finish_special(s, fax, endtime)
    fax.tout = s.now
    fax.state = FINISHED
    if s.now < endtime
        request!(s, s.state.special_queue, (x,y) -> serve_special(x, y, endtime))
    end
end

# produce a snapshot of simulation state
# if final is true additionaly print summary statistics
function report(s, final=false)
    x = s.state.faxes
    println("time: ", s.now)
    for inst in instances(STATE)
        print("\t", lpad(inst,8),": ", count(y -> y.state == inst, x))
    end
    print("\n\t\t\t   idle:  ", length(s.state.entry_queue.requests))
    println("\t\t\t idle:\t  ", length(s.state.special_queue.requests))
    if final
        @printf("%% special:\t%6.4f\n",
                mean(f.special for f in x))
        @printf("avg normal:\t%6.4f\n",
                mean(f.tout-f.tin for f in x if !f.special))
        @printf("avg special:\t%6.4f\n",
                mean(f.tout-f.tin for f in x if f.special))
    end
end

# run simulation
# gets number of entry and special clerks on morning and afternoon shifts
# if tracing is true then reporting is enabled
function run(ec1, sc1, ec2, sc2, tracing::Bool=false)
    sc = ServiceCenter()
    s = Scheduler(sc)
    if tracing
        # print state every hour
        # take special care of 12:00, when there is a change of shift
        for t in [8, 9, 10, 11, 12-eps(12.0),
                  12+eps(12.0), 13, 14, 15, 16]
            register!(s, report, t)
        end
    end
    # plan to change shift at 12:00, remember that s.now is 0.0 at this point
    register!(s, x -> change_shift(x, ec2, sc2), 12.0)

    # morning clerks are scheduled; they finish job at 12:00
    for i in 1:ec1
        request!(s, sc.entry_queue, (x,y) -> serve_entry(x, y, 12.0))
    end
    for i in 1:sc1
        request!(s, sc.special_queue, (x,y) -> serve_special(x, y, 12.0))
    end

    # set time to 8:00
    s.now = 8.0
    # initiate fax arrival process
    register!(s, arrival, rand(sc.ia[8.0]))
    # and open the shop!
    # simulation will finish naturally when all faxes are processed
    go!(s, Inf)
    tracing && report(s, true)

    # return SLA statistics
    esla = mean(f.tout-f.tin<1/6 for f in sc.faxes if !f.special)
    ssla = mean(f.tout-f.tin<1/6 for f in sc.faxes if f.special)
    return esla, ssla, esla > 0.96 && ssla > 0.8
end

# example run
println("Scenario (16,6,7,3):")
run(16, 6, 7, 3, true)

# check if any other scenario with less or equal number of clerks
# ensures that SLA is met with probability higher than 99%
n = 2048
m = mean(run(17, 7, 8, 4)[3] for i in 1:n)
s = sqrt(m*(1-m)/n)
l, h = m .+ [-2s, 2s]
@printf("\nSLA prob. (95%% CI), scenario (17,7,8,4): %6.4f (%6.4f, %6.4f)", m, l, h)


text:
import Base: A_mul_B!, convert

using BandedMatrices

function Base.convert{BM<:BandedMatrix,T<:Union{SymTridiagonal,Tridiagonal}}(::Type{BM},M::T)
    m,n = size(M,1),size(M,2)
    ret=BandedMatrix(eltype(BM)==Any?eltype(M):promote_type(eltype(BM),eltype(M)),m,n,1,1)

    for i = 1:m
        ret[i,i] = M[i,i]
    end
    for i = 1:m-1
        ret[i,i+1] = M[i,i+1]
        ret[i+1,i] = M[i+1,i]
    end
    ret
end

function Base.convert{BM<:BandedMatrix}(::Type{BM},M::Diagonal)
    m,n = size(M,1),size(M,2)
    ret=BandedMatrix(eltype(BM)==Any?eltype(M):promote_type(eltype(BM),eltype(M)),m,n,0,0)

    for i = 1:m
        ret[i,i] = M[i,i]
    end
    ret
end

A_mul_B!{T<:Number}(α::T, A::BandedMatrix{T},
                    x::AbstractVector{T},
                    β::T, y::AbstractVector{T}) =
                        BLAS.gbmv!('N',A.m,A.l,A.u,α,A.data,x,β,y)


text:
using Base.Test
using Flux # Dense
using BSON # BSON.@load

BSON.@load joinpath(@__DIR__, "mymodel.bson") model

@test model.layers[1] isa Dense
@test size(model.layers[1].W) == (5, 10)


text:
## import packages
    include("load_packages.jl")

# Choose the benchmark which should be evaluated
## Benchmarks used in our article:
### 1_solution_gap_10
### 2_sensitivity_ind_100
### 2_sensitivity_md_100
### 2_sensitivity_hd_100
### 3_calculation_time_ind
### 3_calculation_time_md
### 3_calculation_time_hd
    experiment = "3_calculation_time_company"

## Alternatively one could specify new transactional data sets and capacity constellations.
## To see how the data has to be specified take a look at the "capacity_***" and "transactions_***" data.

# Set the number of cpu cores your computer has at its disposal
    cpu_cores  = 8
    ren_lock = ReentrantLock()

# Choose Optimisations and Heuristics to evaluate in the benchmark
    start = DataFrame(QMKOPT = [0], # quadratic-multiple knapsack heuristic with CPLEX as solver
                      QMK    = [1], # quadratic-multiple knapsack heuristic with SBB as solver
                      QMKLOC = [0], # quadratic-multiple knapsack heuristic with SBB as solver + local search based on the QMK objective function
                      CHI    = [1], # chi-square heuristic
                      CHILOC = [1], # chi-square heuristic + local search based on the QMK objective function
                      KLINK  = [0], # K-LINK heuristic by Zhang, W.-H. Lin, M. Huang and X. Hu (2021) https://doi.org/10.1016/j.ejor.2019.07.004
                      GP     = [1], # greedy pairs heuristic by A. Catalan and M. Fisher (2012) https://doi.org/10.2139/ssrn.2166687
                      GPLOC  = [0], # greedy pairs heuristic by A. Catalan and M. Fisher (2012) https://doi.org/10.2139/ssrn.2166687 + local search
                      GS     = [1], # greedy seeds heuristic by A. Catalan and M. Fisher (2012) https://doi.org/10.2139/ssrn.2166687
                      GSLOC  = [0], # greedy seeds heuristic by A. Catalan and M. Fisher (2012) https://doi.org/10.2139/ssrn.2166687 + local search
                      BS     = [1], # bestselling heuristic by A. Catalan and M. Fisher (2012) https://doi.org/10.2139/ssrn.2166687
                      BSLOC  = [0], # bestselling heuristic by A. Catalan and M. Fisher (2012) https://doi.org/10.2139/ssrn.2166687  + local search
                      OPT    = [0], # optimisation model to determine the optimal solution with CPLEX
                      RND    = [1]) # random allocation of SKUs (cannot be deactivated)

# Parameter for the transaction generation if no transactional data is specified under the path
# "transactions/transactions_$experiment". The benchmark will generate random and independent transactions
# while "order" specifies the number of transactions in each dataset. The parameter max_dependence
# specifies the maximal strength of the dependecies between products. max-groupsize specifies the
# maximal group size while group link specifies the ratio of outlinks of each group. Finally, ind_chance
# can be used to determine the chance that instead of group-allocations an independent product is ordered.
    orders = 500000
    min_dependence = 0.0
    max_dependence = 0.3
    max_groupsize = 100
    group_link = 0.01
    ind_chance = 0.01

# Parameters for the KLINK heuristic
## trials: number of different trials with a completly new random solution
## stagnant: number of local search iterations without improvement till termination
## strategy: select the local search strategy of the KLINK heuristic
###          -> 1: move SKU to a different warehouse if it improves the objective and the warehouse has capacity cap_left
###          -> 2: pair-wise exchange between SKUs if it improves the objective
## klinkstatus: display the current iteration during the heuristic
    trials   = 100
    stagnant = 10
    strategy = 2
    klinkstatus = 0

# Parameters for all Optimisations
## abort: number of seconds until the Optimisation is aborted
## show_opt: specify whether the status of the optimisation should be shown
## allowed_gap: specify the termination criterion in case a gap is allowed in the optimisation
## max_nodes: maximum number of nodes till termination
    abort       = 43200
    show_opt    = 0
    allowed_gap = 0.00000
    max_nodes   = 10000000

# Parameters for CHISQUARE
## sig: significance level alpha for the chi-square tests
    sig = 0.05

# Parameters for RANDOM
## iterations: number of different random allocations for the comparison
    iterations = 10

# Initialise the basic problem by loading the respective capacity constellations
## capacity_benchmark: capacity matrix with column = capacity and row = constellation
    capacity_benchmark  = readdlm("capacity/capacity_$experiment.csv", ';', Int64)
    skus_benchmark      = vec(readdlm("capacity/skus_$experiment.csv", ';', Int64))

# Run the benchmark
    parcels_benchmark,
    time_benchmark,
    cap_used,
    parcel_reduction,
    split_reduction,
    gap_optimisation = BENCHMARK(capacity_benchmark::Array{Int64,2},
                                 skus_benchmark::Vector{Int64},
                                 start::DataFrame,
                                 orders::Int64,
                                 max_groupsize::Int64,
                                 max_dependence::Float64,
                                 trials::Int64,
                                 stagnant::Int64,
                                 strategy::Int64,
                                 klinkstatus::Int64,
                                 abort::Int64,
                                 iterations::Int64,
                                 show_opt::Int64,
                                 cpu_cores::Int64,
                                 allowed_gap::Float64,
                                 max_nodes::Int64,
                                 sig::Float64)

# Export the results
CSV.write("results/parcels_sent_$experiment.csv", parcels_benchmark)
CSV.write("results/duration_$experiment.csv", time_benchmark)
CSV.write("results/capacity_used_$experiment.csv", cap_used)
CSV.write("results/parcel_reduction_$experiment.csv", parcel_reduction)
CSV.write("results/split_reduction_$experiment.csv", split_reduction)
CSV.write("results/optimisation_gap_$experiment.csv", gap_optimisation)

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
files = [
    "atlas"
    "axes"
    "chain_in_chain"
    "dice_nofriction"
    "dice_tiltedplane"
    "dice"
    "disconnected_bodies"
    "doublependulum_3d"
    "doublependulum_disconnection"
    "doublependulum_urdf"
    "elliptic_joint"
    "football"
    "fourbar_disconnection"
    "inverted_pyramid_plane"
    "joint_force"
    "joint_torque"
    "linear_pendulum"
    "n_fourbars"
    "n_pendulum"
    "nutation"
    "pendulum_forced"
    "pendulum"
    "planar_example"
    "scissor_lift"
    "slider_crank"
    "slider_crank3d"
    "wheel"
]

controlled = [
    "doublependulum_disconnection"
    "joint_force"
    "joint_torque"
    "pendulum_forced"
    "nutation"
    "football"
    "fourbar_disconnection"
    "slider_crank"
]

for i=1:length(files)
    include("../test/examples/"*files[i]*".jl")
    mech.g = 0.0

    steps = Base.OneTo(100)
    storage = Storage{Float64}(steps,length(mech.bodies))

    if files[i] in controlled
        control_function = eval(Meta.parse(files[i]*"_control!"))
        SUITE[files[i]] = @benchmarkable simulate!($mech, $steps, $storage, $control_function) samples=2
    else
        SUITE[files[i]] = @benchmarkable simulate!($mech, $steps, $storage) samples=2
    end
end


text:
module Stiffness

export StiffnessData

struct StiffnessData
    grads::Array{Float64, 3}
    gq::Array{Float64, 3}       #grads of 6 basis functions at 3 quadrature points.

    StiffnessData() = new(zeros(2, 3, 6),
                          reshape([-1.0,-1.0,1.0,1.0,-1.0,-1.0,1.0,0.0,1.0,0.0,-1.0,0.0,0.0,
                                   -1.0,0.0,1.0,0.0,1.0,0.0,-2.0,-2.0,-2.0,2.0,0.0,0.0,2.0,2.0,
                                   2.0,2.0,0.0,0.0,2.0,-2.0,-2.0,-2.0,0], 2, 3, 6)
                          )
end

function op!(S::StiffnessData, x::Array{Float64, 1}, y::Array{Float64, 1}, m::Array{Float64, 1})
    a11 = -y[1] + y[3]
    a12 =  y[1] - y[2]
    a21 =  x[1] - x[3]
    a22 = -x[1] + x[2]

    for f=1:6
        @simd for p=1:3
            S.grads[1, p, f] = a11 * S.gq[1, p, f] + a12 * S.gq[2, p, f]
            S.grads[2, p, f] = a21 * S.gq[1, p, f] + a22 * S.gq[2, p, f]
        end

        # this seems slower:
        # S.grads[1, 1:3, f] = a11 * S.gq[1, 1:3, f] + a12 * S.gq[2, 1:3, f]
        # S.grads[2, 1:3, f] = a21 * S.gq[1, 1:3, f] + a22 * S.gq[2, 1:3, f]
    end

    # This seems slower than the for loop above
    #@. @views S.grads[1, :, :] = a11 * S.gq[1, :, :] + a12 * S.gq[2, :, :]
    #@. @views S.grads[2, :, :] = a21 * S.gq[2, :, :] + a22 * S.gq[2, :, :]

    det = -(x[2] - x[3])*y[1] + (x[1] - x[3])*y[2] - (x[1] - x[2])*y[3]
    dv = 1.0 / (6.0 * det)
    ii = 1
    for i=1:6
        for j=1:i
            s = 0.0
            @simd for k=1:3
                s += S.grads[1, k, i] * S.grads[1, k, j] + S.grads[2, k, i] * S.grads[2, k, j]
            end
            m[ii] = dv * s

            # this seems slower:
            # m[ii] = dv * (dot(S.grads[1, :, i], S.grads[1, :, j]) +
            #         dot(S.grads[2, :, i], S.grads[2, :, j]))

            ii += 1
        end
    end

    return m
end

end

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
@testset "I341" begin
    @test messageset(lintstr("""
    f(; a=1) = a
    a = (:a, 1)
    f(; a)
    """)) == Set([:I341])

    @test messageset(lintstr("""
    f(; a=1) = a
    a = :b
    f(; a => 1)
    """)) == Set([:I341])

end


text:
# ---------------------------------------------------
module Sim
# ---------------------------------------------------
using ..Deuron
using ..Config
using ..Logs

# Main sim entry point
function entry()
    trace = try
        if Config.exit_state(Deuron.core.config) == "Paused"
            Logs.log_to_app(Deuron.core.logs, "Simulation resuming", true)
        else
            Logs.log_to_app(Deuron.core.logs, "Simulation starting", true)
        end

        while true
            if !Deuron.paused
                Logs.log_to_app(Deuron.core.logs, ".", true)
                sleep(1)
            else
                Logs.log_to_app(Deuron.core.logs, "--paused--", true)
                sleep(1)
            end

            if Deuron.reset
                Logs.log_to_app(Deuron.core.logs, "Sim resetting...", true)
                Deuron.resume()
                Deuron.reset_complete()
                Logs.log_to_app(Deuron.core.logs, "Reset complete", true)
            end

            if Deuron.stopped
                break
            end
        end

        exit()

        nothing
    catch ex
        bt = catch_backtrace()
        msg = sprint(showerror, ex, bt)

        Logs.log_to_err(Deuron.core.logs, msg)

        # Set flag and exit
        Deuron.set_crashed()
        exit()

        # println("### Exception ###:\n", ex)
        # stacktrace(catch_backtrace())
        # Base.exit(3)
    end

    # Not sure if this will ever trigger. TODO research...
    if trace ≠ nothing
        println(trace)
        println("\n################ Exception ####################")
        println("Goodbye.")
        Base.exit(2)
    end
end

function exit()
    msg = "============================================="
    Logs.log_to_app(Deuron.core.logs, msg, true)
    println("\n", msg)

    # Do exit stuff
    if Deuron.suspended
        Logs.log_to_app(Deuron.core.logs, "Suspending...", true)
        Logs.log_to_app(Deuron.core.logs, "Suspend complete", true)
        Config.set_exit_state(Deuron.core.config, "Paused")
    elseif Deuron.crashed
        Logs.log_to_app(Deuron.core.logs, "Crashed!!! See error.log", true)
        Config.set_exit_state(Deuron.core.config, "Crashed")
    else
        msg = "Exiting..."
        Logs.log_to_app(Deuron.core.logs, msg, true)
        Config.set_exit_state(Deuron.core.config, "Terminated")
    end

    Config.save(Deuron.core.config)

    Deuron.stop()

    msg = "Goodbye."
    Logs.log_to_app(Deuron.core.logs, msg, true)

    Logs.close_logs(Deuron.core.logs)

    println(msg)
    Base.exit(1)
end

end # End module -----------------


text:
"""
    Gamma

* Parametrization
parametrized in \\alpha and \\theta

* Score

* Fisher Information

* `time_varying_params` map.

* Default link
"""
Gamma

function score!(score_til::Matrix{T}, y::T, ::Type{Gamma}, param::Matrix{T}, t::Int) where T
    score_til[t, 1] = log(y) - digamma(param[t, 1]) - log(param[t, 2])
    score_til[t, 2] = y/param[t, 2]^2 - param[t, 1]/param[t, 2]
    return
end

function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{Gamma}, param::Matrix{T}, t::Int) where T
    aux.fisher[1, 1] = trigamma(param[t, 1])
    aux.fisher[2, 2] = param[t, 1] / param[t, 2]^2
    aux.fisher[2, 1] = 1/param[t, 2]
    aux.fisher[1, 2] = 1/param[t, 2]
    return
end

function log_likelihood(::Type{Gamma}, y::Vector{T}, param::Matrix{T}, n::Int) where T
    loglik = 0.0
    for t in 1:n
        loglik += (param[t, 1] - 1)*log(y[t]) - y[t]/param[t, 2] - loggamma(param[t, 1]) - param[t, 1]*log(param[t, 2])
    end
    return -loglik
end

# Links
function link!(param_tilde::Matrix{T}, ::Type{Gamma}, param::Matrix{T}, t::Int) where T
    param_tilde[t, 1] = link(LogLink, param[t, 1], zero(T))
    param_tilde[t, 2] = link(LogLink, param[t, 2], zero(T))
    return
end
function unlink!(param::Matrix{T}, ::Type{Gamma}, param_tilde::Matrix{T}, t::Int) where T
    param[t, 1] = unlink(LogLink, param_tilde[t, 1], zero(T))
    param[t, 2] = unlink(LogLink, param_tilde[t, 2], zero(T))
    return
end
function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{Gamma}, param::Matrix{T}, t::Int) where T
    aux.jac[1] = jacobian_link(LogLink, param[t, 1], zero(T))
    aux.jac[2] = jacobian_link(LogLink, param[t, 2], zero(T))
    return
end

# utils
function update_dist(::Type{Gamma}, param::Matrix{T}, t::Int) where T
    small_threshold!(param, SMALL_NUM, t)
    return Gamma(param[t, 1], param[t, 2])
end

function params_sdm(d::Gamma)
    return Distributions.params(d)
end

function num_params(::Type{Gamma})
    return 2
end


text:
using Base.BinaryPlatforms
import MPIABI

function augment_platform!(platform)
    abi = MPIABI.abi
    BinaryPlatforms.add_tag!(platform.tags, "mpi", string(abi))
    return platform
end


text:
let
    pomdp = BabyPOMDP()

    bu = DiscreteUpdater(pomdp)
    b0 = initialize_belief(bu, initialstate_distribution(pomdp))

    # these values were gotten from FIB.jl
    # alphas = [-29.4557 -36.5093; -19.4557 -16.0629]
    alphas = [ -16.0629 -19.4557; -36.5093 -29.4557]
    policy = AlphaVectorPolicy(pomdp, alphas, ordered_actions(pomdp))

    @test Set(alphapairs(policy)) == Set([[-16.0629, -36.5093]=>false, [-19.4557, -29.4557]=>true])
    @test Set(alphavectors(policy)) == Set([[-16.0629, -36.5093], [-19.4557, -29.4557]])

    # initial belief is 100% confidence in baby not being hungry
    @test isapprox(value(policy, b0), -16.0629)
    @test isapprox(value(policy, [1.0,0.0]), -16.0629)
    @test isapprox(actionvalues(policy, b0), [-16.0629, -19.4557])
    @test length(actionvalues(policy, b0)) == length(actions(pomdp))

    # because baby isn't hungry, policy should not feed (return false)
    @test action(policy, b0) == false

    # SparseCat belief
    sparse_b0 = SparseCat([s for s in ordered_states(pomdp) if pdf(b0, s) != 0.],
                        b0.b[b0.b .> 0.])
    @test isapprox(value(policy, sparse_b0), -16.0629)
    @test isapprox(actionvalues(policy, sparse_b0), [-16.0629, -19.4557])
    @test action(policy, sparse_b0) == false

    # Bool_distribution (if it works for this, it should work for an arbitrary distribution)
    bd = initialstate_distribution(pomdp)::BoolDistribution
    @test action(policy, bd) == false

    # try pushing new vector
    push!(policy, [0.0,0.0], true)

    @test value(policy, b0) == 0.0
    @test action(policy, b0) == true
end


text:
module ReinforcementLearningCore

using ReinforcementLearningBase

const RLCore = ReinforcementLearningCore

@doc """
[ReinforcementLearningCore.jl](https://github.com/JuliaReinforcementLearning/ReinforcementLearningCore.jl) (**RLCore**)
provides some standard and reusable components defined by [**RLBase**](https://github.com/JuliaReinforcementLearning/ReinforcementLearningBase.jl), hoping that they are useful for people to implement and experiment with different kinds of algorithms.
""" RLCore

export RLCore

include("extensions/extensions.jl")
include("utils/utils.jl")
include("core/core.jl")
include("components/components.jl")

end # module


text:
import JLD2, LinearAlgebra, JSON

# erlangDParams = Dict()
JLD2.@load pwd()*"/src/erlangParamsData/erlangDParams.jld2" erlangDParams

JLD2.@load pwd()*"/src/CMEParamsData/CMEParams.jld2" CMEParams

"""
ME constructor method

    ME(
        a::Array{<:Real,2},
        S::Array{<:Real,2},
        s::Array{<:Real,1},
        D::Array{<:Real}=[0],
    )

Inputs:
 - `a` a 1 by p Array of reals
 - `S` a p by p Array of reals
 - `s` a p by 1 Array of reals
 - `D` an optional argument, if empty then me.D is the identity,
    else is a p by p matrix
 Throws an error if the dimensions are inconsistent.
"""
struct ME
    a::Array{<:Real,2}
    S::Array{<:Real,2}
    s::Union{Array{<:Real,1},Array{<:Real,2}}
    D::Array{<:Real,2}
    function ME(
        a::Array{<:Real,2},
        S::Array{<:Real,2},
        s::Union{Array{<:Real,1},Array{<:Real,2}};
        D::Array{<:Real,2}=zeros(1,1),
    )

        if D==zeros(1,1)
            D = Array{Float64}(LinearAlgebra.I(size(S,1)))
        end
        s1 = size(a,1)
        s2 = size(a,2)
        s3 = size(S,1)
        s4 = size(S,2)
        s5 = size(s,1)
        s6 = size(s,2)
        s7 = size(D,1)
        s8 = size(D,2)
        checksquare(D,"D")
        checksquare(S,"S")
        test = (s1!=1) || (s6!=1) || any(([s2;s3;s4;s7;s8].-s5).!=0)
        if test
            error("Dimensions of ME representation not consistent")
        else
            return new(a,S,s,D)
        end
    end
end

pdf(me::SFFM.ME) = x->(me.a*exp(me.S*x)*me.s)[1]
pdf(a::Array{<:Real,2}, me::SFFM.ME) =
    (length(a)==size(me.S,1)) ? (x->(a*exp(me.S*x)*me.s)[1]) : throw(
        DomainError("a and me.S must have compatible size"))

pdf(me::SFFM.ME, x::Real) = pdf(me)(x)
pdf(a::Array{<:Real,2}, me::SFFM.ME, x::Real) = pdf(a,me)(x)

pdf(me::SFFM.ME, x::Array{<:Real}) = pdf(me).(x)
pdf(a::Array{<:Real,2}, me::SFFM.ME, x::Array{<:Real}) = pdf(a,me).(x)

ccdf(me::SFFM.ME) = x->sum(me.a*exp(me.S*x))
ccdf(a::Array{<:Real,2}, me::SFFM.ME) = (length(a)==size(me.S,1)) ? (x->sum(a*exp(me.S*x))) : throw(
    DomainError("a and me.S must have compatible size"))

ccdf(me::SFFM.ME, x::Real) = ccdf(me)(x)
ccdf(a::Array{<:Real,2}, me::SFFM.ME, x::Real) = ccdf(a,me)(x)

ccdf(me::SFFM.ME, x::Array{<:Real}) = ccdf(me).(x)
ccdf(a::Array{<:Real,2}, me::SFFM.ME, x::Array{<:Real}) = ccdf(a,me).(x)

cdf(me::SFFM.ME) = x->1-ccdf(me,x)
cdf(a::Array{<:Real,2}, me::SFFM.ME) = x->1-ccdf(a,me,x)

cdf(me::SFFM.ME, x::Real) = cdf(me)(x)
cdf(a::Array{<:Real,2}, me::SFFM.ME, x::Real) = cdf(a,me)(x)

cdf(me::SFFM.ME, x::Array{<:Real}) = cdf(me).(x)
cdf(a::Array{<:Real,2}, me::SFFM.ME, x::Array{<:Real}) = cdf(a,me).(x)

"""

"""
function MakeME(params; mean::Real = 1)
    N = 2*params["n"]+1
    α = zeros(1,N)
    α[1] = params["c"]
    a = params["a"]
    b = params["b"]
    ω =  params["omega"]
    for k in 1:params["n"]
        kω = k*ω
        α[2*k] = (1/2)*( a[k]*(1+kω) - b[k]*(1-kω) )/(1+kω^2)
        α[2*k+1] = (1/2)*( a[k]*(1-kω) + b[k]*(1+kω) )/(1+kω^2)
    end
    α = α./sum(α)
    Q = zeros(N,N)
    Q[1,1] = -1
    for k in 1:params["n"]
        kω = k*ω
        idx = 2*k:(2*k+1)
        Q[idx,idx] = [-1 -kω; kω -1]
    end
    Q = Q.*sum(-α*Q^-1)./mean
    q = -sum(Q,dims=2)
    return SFFM.ME(α,Q,q;D=params["D"])
end

function MakeErlang(order; mean::Real = 1)
    α = zeros(1,order) # inital distribution
    α[1] = 1
    λ = order/mean
    Q = zeros(order,order)
    Q = Q + LinearAlgebra.diagm(0=>repeat(-[λ],order), 1=>repeat([λ],order-1))
    q = -sum(Q,dims=2)
    D = Array{Float64}(LinearAlgebra.I(order))[end:-1:1,:]
    return SFFM.ME(α,Q,q;D=D)
end

orbit(t,me::SFFM.ME; norm = 1) = begin
    orbits = zeros(length(t),length(me.a))
    for i in 1:length(t)
        num = me.a*exp(me.S*t[i])
        if norm == 1
            denom = sum(num)
        elseif norm == 2
            denom = exp(me.S[1,1]*t[i])#sum(num)
        else
            denom = 1
        end
        orbits[i,:] = num./denom
    end
    return orbits
end

function renewalProperties(me::SFFM.ME)
    density(t) = begin
        Q = me.S
        q = me.s
        α = me.a
        e = ones(size(q))
        (α*exp((Q+q*α)*t)*q)[1]
    end
    mean(t) = begin
        Q = me.S
        q = me.s
        α = me.a
        e = ones(size(q))
        temp1 = α*-Q^-1
        temp2 = temp1*e
        temp3 = temp1./temp2
        temp4 = Q + (q*α)
        ((t./temp2) - α*(I - exp(temp4*t))*(temp4 + e*temp3)^-1*q)[1]
    end
    ExpectedOrbit(t) = begin
        Q = me.S
        q = me.s
        α = me.a
        e = ones(size(q))
        (α*exp((Q + q*α)*t))[1]
    end
    return (density=density,mean=mean,ExpectedOrbit=ExpectedOrbit)
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM