text:
"""
Line Outage Distribution Factors (LODFs) are a sensitivity measure of how a change in
a line’s flow affects the flows on other lines in the system.
"""
struct LODF{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Float64}
    data::Array{Float64, 2}
    axes::Ax
    lookup::L
end

function _buildlodf(branches, nodes, dist_slack::Array{Float64} = [0.1])
    linecount = length(branches)
    ptdf, a = PowerSystems._buildptdf(branches, nodes, dist_slack)
    H = gemm('N', 'N', ptdf, a)
    ptdf_denominator = H
    for iline in 1:linecount
        if (1.0 - ptdf_denominator[iline, iline]) < 1.0E-06
            ptdf_denominator[iline, iline] = 0.0
        end
    end
    (Dem, dipiv, dinfo) = getrf!(
        Matrix{Float64}(LinearAlgebra.I, linecount, linecount) -
        Array(LinearAlgebra.Diagonal(ptdf_denominator)),
    )
    lodf = gemm('N', 'N', H, getri!(Dem, dipiv))
    lodf =
        lodf - Array(LinearAlgebra.Diagonal(lodf)) -
        Matrix{Float64}(LinearAlgebra.I, linecount, linecount)
    return lodf
end

"""
Builds the LODF matrix from a group of branches and nodes. The return is a LOLDF array indexed with the branch name.

# Keyword arguments
- `dist_slack::Vector{Float64}`: Vector of weights to be used as distributed slack bus.
    The distributed slack vector has to be the same length as the number of buses
"""
function LODF(branches, nodes, dist_slack::Vector{Float64} = [0.1])

    #Get axis names
    line_ax = [branch.name for branch in branches]
    lodf = _buildlodf(branches, nodes, dist_slack)

    axes = (line_ax, line_ax)
    look_up = (_make_ax_ref(line_ax), _make_ax_ref(line_ax))

    return LODF(lodf, axes, look_up)
end

"""
Builds the LODF matrix from a system. The return is a LOLDF array indexed with the branch name.

# Keyword arguments
- `dist_slack::Vector{Float64}`: Vector of weights to be used as distributed slack bus.
    The distributed slack vector has to be the same length as the number of buses
"""
function LODF(sys::System, dist_slack::Vector{Float64} = [0.1])
    branches = sort!(
        collect(get_components(ACBranch, sys)),
        by = x -> (get_number(get_arc(x).from), get_number(get_arc(x).to)),
    )
    nodes = sort!(collect(get_components(Bus, sys)), by = x -> get_number(x))
    return LODF(branches, nodes, dist_slack)
end


text:
push!(LOAD_PATH, "./src")
using ParametricOptInterface
using MathOptInterface
using GLPK
import Random
#using SparseArrays
using TimerOutputs

const MOI = MathOptInterface
const POI = ParametricOptInterface
SOLVER = GLPK

if SOLVER == GLPK
    MAX_ITER_PARAM = "it_lim"
elseif SOLVER == Gurobi
    MAX_ITER_PARAM = "IterationLimit"
elseif SOLVER == Xpress
    MAX_ITER_PARAM = "LPITERLIMIT"
end

struct PMedianData
    num_facilities::Int
    num_customers::Int
    num_locations::Int
    customer_locations::Vector{Float64}
end

# This is the LP relaxation.
function generate_poi_problem(model, data::PMedianData, add_parameters::Bool)
    NL = data.num_locations
    NC = data.num_customers

    ###
    ### 0 <= facility_variables <= 1
    ###

    facility_variables = MOI.add_variables(model, NL)

    for v in facility_variables
        MOI.add_constraint(model, v, MOI.Interval(0.0, 1.0))
    end

    ###
    ### assignment_variables >= 0
    ###

    assignment_variables = reshape(MOI.add_variables(model, NC * NL), NC, NL)
    for v in assignment_variables
        MOI.add_constraint(model, v, MOI.GreaterThan(0.0))
        # "Less than 1.0" constraint is redundant.
    end

    ###
    ### Objective function
    ###

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            [
                MOI.ScalarAffineTerm(
                    data.customer_locations[i] - j,
                    assignment_variables[i, j]
                )
                for i in 1:NC for j in 1:NL
            ],
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    ###
    ### assignment_variables[i, j] <= facility_variables[j]
    ###

    for i in 1:NC, j in 1:NL
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                [
                    MOI.ScalarAffineTerm(1.0, assignment_variables[i, j]),
                    MOI.ScalarAffineTerm(-1.0, facility_variables[j])
                ],
                0.0,
            ),
            MOI.LessThan(0.0),
        )
    end

    ###
    ### sum_j assignment_variables[i, j] = 1
    ###

    for i in 1:NC
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                [
                    MOI.ScalarAffineTerm(1.0, assignment_variables[i, j])
                    for j in 1:NL
                ],
                0.0,
            ),
            MOI.EqualTo(1.0),
        )
    end

    ###
    ### sum_j facility_variables[j] == num_facilities
    ###

    if add_parameters
        d, cd = MOI.add_constrained_variable(model, POI.Parameter(data.num_facilities))
    end

    if add_parameters
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.(1.0, vcat(facility_variables,d)),
                0.0,
            ),
            MOI.EqualTo{Float64}(0),
        )
    else
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.(1.0, facility_variables),
                0.0,
            ),
            MOI.EqualTo{Float64}(data.num_facilities),
        )
    end

    return assignment_variables, facility_variables
end

function solve_moi(data::PMedianData, optimizer; vector_version, params, add_parameters=false)
    model = optimizer()
    for (param, value) in params
        MOI.set(model, param, value)
    end
    @timeit "generate" x, y = if vector_version
        generate_poi_problem_vector(model, data, add_parameters)
    else
        generate_poi_problem(model, data, add_parameters)
    end
    @timeit "solve" MOI.optimize!(model)
    return MOI.get(model, MOI.ObjectiveValue())
end

function POI_OPTIMIZER()
    return POI.Optimizer(SOLVER.Optimizer())
end

function MOI_OPTIMIZER()
    return SOLVER.Optimizer()
end

function solve_moi_loop(data::PMedianData; vector_version, max_iters=Inf, time_limit_sec=Inf, loops)
    params = []
    if isfinite(time_limit_sec)
        push!(params, (MOI.TimeLimitSec(), time_limit_sec))
    end
    if isfinite(max_iters)
        push!(params, (MOI.RawOptimizerAttribute(MAX_ITER_PARAM), max_iters))
    end
    push!(params, (MOI.Silent(),true))
    s_type = vector_version ? "vector" : "scalar"

    @timeit(
        "$(SOLVER) MOI $(s_type)",
        for _ in 1:loops
            solve_moi(
                data, MOI_OPTIMIZER; vector_version=vector_version, params=params
            )
        end
    )
end

function solve_poi_no_params_loop(data::PMedianData; vector_version, max_iters=Inf, time_limit_sec=Inf, loops)
    params = []
    if isfinite(time_limit_sec)
        push!(params, (MOI.TimeLimitSec(), time_limit_sec))
    end
    if isfinite(max_iters)
        push!(params, (MOI.RawOptimizerAttribute(MAX_ITER_PARAM), max_iters))
    end
    push!(params, (MOI.Silent(),true))
    s_type = vector_version ? "vector" : "scalar"
    @timeit(
        "$(SOLVER) POI NO PARAMS $(s_type)",
        for _ in 1:loops
            solve_moi(
                data,POI_OPTIMIZER; vector_version=vector_version, params=params
            )
        end
    )
end

function solve_poi_loop(data::PMedianData; vector_version, max_iters=Inf, time_limit_sec=Inf, loops=1)
    params = []
    if isfinite(time_limit_sec)
        push!(params, (MOI.TimeLimitSec(), time_limit_sec))
    end
    if isfinite(max_iters)
        push!(params, (MOI.RawOptimizerAttribute(MAX_ITER_PARAM), max_iters))
    end
    push!(params, (MOI.Silent(),true))
    s_type = vector_version ? "vector" : "scalar"

    @timeit(
        "$(SOLVER) POI $(s_type)",
        for _ in 1:loops
            solve_moi(
                data,POI_OPTIMIZER; vector_version=vector_version, params=params, add_parameters = true
            )
        end
    )
end

function run_benchmark(;
    num_facilities, num_customers, num_locations, time_limit_sec, max_iters, loops
)
    Random.seed!(10)
    reset_timer!()
    data = PMedianData(num_facilities, num_customers, num_locations, rand(num_customers) .* num_locations)
    GC.gc()
    solve_moi_loop(data, vector_version=false, max_iters=max_iters, time_limit_sec=time_limit_sec, loops=loops)
    GC.gc()
    solve_poi_no_params_loop(data, vector_version=false, max_iters=max_iters, time_limit_sec=time_limit_sec, loops=loops)
    GC.gc()
    solve_poi_loop(data, vector_version=false, max_iters=max_iters, time_limit_sec=time_limit_sec, loops=loops)
    GC.gc()
    print_timer()
    println()
end

run_benchmark(
            num_facilities = 100,
            num_customers = 100,
            num_locations = 100,
            time_limit_sec = 0.0001,
            max_iters = 1,
            loops = 1
        )

run_benchmark(
    num_facilities = 100,
    num_customers = 100,
    num_locations = 100,
    time_limit_sec = 0.0001,
    max_iters = 1,
    loops = 100
)

ext: jl


text:
using Plots

function main()
    x = 0:0.1:6
    y = sin.(x)

    plot(x, y)
end

main()


text:
export Slicing, Slicer, SlicedEinsum

struct Slicer
    log2_sizes::Vector{Float64}   # the size dict after slicing
    legs::Dict{Int,Float64}   # sliced leg and its original size
    max_size::Int       # maximum number of sliced legs
end
Slicer(log2_sizes::AbstractVector{Float64}, max_size::Int) = Slicer(collect(log2_sizes), Dict{Int,Float64}(), max_size)
Base.length(s::Slicer) = length(s.legs)
function Base.replace!(slicer::Slicer, pair::Pair)
    worst, best = pair
    @assert haskey(slicer.legs, worst)
    @assert !haskey(slicer.legs, best)
    slicer.log2_sizes[worst] = slicer.legs[worst]       # restore worst size
    slicer.legs[best] = slicer.log2_sizes[best]  # add best to legs
    slicer.log2_sizes[best] = 0.0
    delete!(slicer.legs, worst)                  # remove worst from legs
    return slicer
end

function Base.push!(slicer, best)
    @assert length(slicer) < slicer.max_size
    @assert !haskey(slicer.legs, best)
    slicer.legs[best] = slicer.log2_sizes[best]  # add best to legs
    slicer.log2_sizes[best] = 0.0
    return slicer
end

struct Slicing{LT}
    legs::Vector{LT}   # sliced leg and its original size
end
Base.:(==)(se::Slicing, se2::Slicing) = se.legs == se2.legs

Slicing(s::Slicer, inverse_map) = Slicing([inverse_map[l] for (l, s) in s.legs])
Base.length(s::Slicing) = length(s.legs)

struct SlicedEinsum{LT, Ein} <: AbstractEinsum
    slicing::Slicing{LT}
    eins::Ein
end
Base.:(==)(se::SlicedEinsum, se2::SlicedEinsum) = se.slicing == se2.slicing && se.eins == se2.eins

# Iterate over tensor network slices, its iterator interface returns `slicemap` as a Dict
# slice and fill tensors with
# * take_slice(x, label_of_x, slicemap)
# * fill_slice!(x, label_of_x, x_slice, slicemap)
struct SliceIterator{LT,IT<:CartesianIndices}
    ixsv::Vector{Vector{LT}}
    iyv::Vector{LT}
    sliced_labels::Vector{LT}
    indices::IT
    size_dict_sliced::Dict{LT,Int}
end

function SliceIterator(se::SlicedEinsum, size_dict::Dict{LT}) where LT
    iyv = OMEinsum.getiyv(se.eins.eins)
    ixsv = OMEinsum.getixsv(se.eins)
    n = length(se.slicing)

    size_dict_sliced = copy(size_dict)
    sliced_sizes = Vector{Int}(undef, n)
    sliced_labels = Vector{LT}(undef, n)
    for i = 1:n
        l = se.slicing.legs[i]
        sliced_sizes[i] = size_dict[l]
        sliced_labels[i] = l
        size_dict_sliced[l] = 1
    end
    indices = CartesianIndices((sliced_sizes...,))
    SliceIterator(ixsv, iyv, sliced_labels, indices, size_dict_sliced)
end
Base.length(si::SliceIterator) = length(si.indices)
Base.eltype(::Type{SliceIterator{LT,IT}}) where {LT,IT} = Dict{LT,Int}

# returns `slicemap` as a Dict
function Base.iterate(si::SliceIterator)
    ci, cistate = iterate(si.indices)
    slicemap = Dict(zip(si.sliced_labels, ones(Int,length(si.sliced_labels))))
    slicemap, (1,(ci,cistate),slicemap)
end
function Base.iterate(si::SliceIterator, state)
    i, (ci,cistate), slicemap = state
    if i >= length(si.indices)
        return nothing  # NOTE: ci is same as cistate
    else
        ci, cistate = iterate(si.indices, cistate)
        for (l, v) in zip(si.sliced_labels, ci.I)
            slicemap[l] = v
        end
        return slicemap, (i+1, (ci,cistate), slicemap)
    end
end
function Base.getindex(si::SliceIterator, indices...)
    ci = si.indices[indices...]
    slicemap = Dict(zip(si.sliced_labels, ci.I))
    return slicemap
end

function take_slice(x, ix, slicemap::Dict)
    slices = map(l->haskey(slicemap, l) ? slicemap[l] : Colon(), ix)
    if all(x->x isa Integer, slices)
        return copy(view(x,slices...))
    else
        return x[slices...]
    end
end
function fill_slice!(x, ix, chunk, slicemap::Dict)
    if ndims(x) == 0
        x .+= chunk  # to avoid CUDA `getindex!`.
    else
        slices = map(l->haskey(slicemap, l) ? slicemap[l] : Colon(), ix)
        view(x, slices...) .+= chunk
    end
end

function (se::SlicedEinsum{LT,ET})(@nospecialize(xs::AbstractArray...); size_info = nothing, kwargs...) where {LT, ET}
    length(se.slicing) == 0 && return se.eins(xs...; size_info=size_info, kwargs...)
    size_dict = size_info===nothing ? Dict{OMEinsum.labeltype(se),Int}() : copy(size_info)
    OMEinsum.get_size_dict!(se, xs, size_dict)

    it = SliceIterator(se, size_dict)
    res = OMEinsum.get_output_array(xs, getindex.(Ref(size_dict), it.iyv))
    eins_sliced = drop_slicedim(se.eins, se.slicing)
    for (k, slicemap) in enumerate(it)
        @debug "computing slice $k/$(length(it))"
        xsi = ntuple(i->take_slice(xs[i], it.ixsv[i], slicemap), length(xs))
        resi = einsum(eins_sliced, xsi, it.size_dict_sliced; kwargs...)
        fill_slice!(res, it.iyv, resi, slicemap)
    end
    return res
end

function drop_slicedim(ne::NestedEinsum, slicing::Slicing)
    isleaf(ne) && return ne
    ixs = map(ix->filter(∉(slicing.legs), ix), getixsv(ne.eins))
    iy = filter(∉(slicing.legs), getiyv(ne.eins))
    NestedEinsum(map(arg->drop_slicedim(arg, slicing), ne.args), similar_eincode(ne.eins, ixs, iy))
end
similar_eincode(::DynamicEinCode, ixs, iy) = DynamicEinCode(ixs, iy)
similar_eincode(::StaticEinCode, ixs, iy) = StaticEinCode{(Tuple.(ixs)...,), (iy...,)}()

# forward some interfaces
function OMEinsum.timespacereadwrite_complexity(code::SlicedEinsum, size_dict)
    size_dict_sliced = copy(size_dict)
    for l in code.slicing.legs
        size_dict_sliced[l] = 1
    end
    tc, sc, rw = timespacereadwrite_complexity(code.eins, size_dict_sliced)
    sliceoverhead = sum(log2.(getindex.(Ref(size_dict), code.slicing.legs)))
    tc + sliceoverhead, sc, rw+sliceoverhead
end
function OMEinsum.flop(code::SlicedEinsum, size_dict)
    size_dict_sliced = copy(size_dict)
    for l in code.slicing.legs
        size_dict_sliced[l] = 1
    end
    fl = flop(code.eins, size_dict_sliced)
    fl * prod(getindex.(Ref(size_dict), code.slicing.legs))
end

OMEinsum.flatten(se::SlicedEinsum) = OMEinsum.flatten(se.eins)
OMEinsum.labeltype(::SlicedEinsum{LT}) where LT = LT
OMEinsum.get_size_dict!(se::SlicedEinsum, xs, size_info::Dict) = OMEinsum.get_size_dict!(se.eins, xs, size_info)
OMEinsum.getixsv(se::SlicedEinsum) = OMEinsum.getixsv(se.eins)
OMEinsum.getiyv(se::SlicedEinsum) = OMEinsum.getiyv(se.eins)
OMEinsum.label_elimination_order(code::SlicedEinsum) = OMEinsum.label_elimination_order(code.eins)


text:
glorot_uniform(rng::Random.AbstractRNG, dims...) = begin; (rand(rng, Float32, dims...) .- 0.5f0) .* sqrt(24.0f0/sum(dims)) end;
glorot_normal(rng::Random.AbstractRNG, dims...) = randn(rng, Float32, dims...) .* sqrt(2.0f0/sum(dims))

module ExpUtils

import ..JuliaRL
using ..Resampling
using Reproduce

export algorithm_settings!, build_algorithm_dict

function algorithm_settings!(s::ArgParseSettings)
    @add_arg_table s begin
        "--inc"
        action = :store_true
        help = "set to use incremental td"
        # Algorithms
        "--is"
        action = :store_true
        help = "set to include base off-policy importance sampling"
        "--normis"
        action = :store_true
        help = "set to include base off-policy importance sampling"
        "--incnormis"
        action = :store_true
        help = "set to include base off-policy importance sampling"
        "--wsnormis"
        action = :store_true
        help = "set to include base off-policy importance sampling"
        "--wsavgnormis"
        action = :store_true
        help = "set to include base off-policy importance sampling"
        "--ir"
        action = :store_true
        help = "Set to include ir"
        "--bcir"
        action = :store_true
        help = "Set to include incremental buffer corrected ir"
        "--wisbatch"
        action = :store_true
        help = "Set to include wisbatch in experiments"
        "--wisbuffer"
        action = :store_true
        help = "Set to include wisbuffer in experiments"
        "--wisoptimal"
        action = :store_true
        help = "Set to run wis optimal experiments"
        "--wisrupam"
        action = :store_true
        help = "Set to run wis rupam experiments"
        "--res_wisrupam"
        action = :store_true
        help = "Set to run wis rupam experiments"
        "--init_u"
        nargs='+'
        arg_type = Float64
        "--normalize_eta"
        action= :store_true
        help = "Set to normalize eta for WIS-TD(0)"
        "--sarsa"
        action = :store_true
        help = "Set to run sarsa"
        "--expectedsarsa"
        action = :store_true
        help = "Set to run expected sarsa"
        "--vtrace"
        action = :store_true
        help = "Set to run clipped is"
        "--clip_value"
        nargs='+'
        arg_type = Float64
        "--clip_value_perc"
        nargs='+'
        arg_type = Float64

    end
    return s
end

function build_algorithm_dict(parsed; max_is=1.0, size_features = 0)
    algo_dict = Dict{String, Resampling.LearningUpdate}()
    sample_dict = Dict{String, String}()
    value_type_dict = Dict{String, String}()

    tdupdate = parsed["inc"] ? IncTD() : BatchTD()

    if parsed["is"]
        algo_dict["IS"] = tdupdate
        sample_dict["IS"] = "ER"
        value_type_dict["IS"] = "State"
    end
    if parsed["normis"]
        algo_dict["NormIS"] = tdupdate
        sample_dict["NormIS"] = "ER"
        value_type_dict["NormIS"] = "State"
    end
    if parsed["incnormis"]
        algo_dict["IncNormIS"] = IncNormIS(tdupdate)
        sample_dict["IncNormIS"] = "ER"
        value_type_dict["IncNormIS"] = "State"
    end
    if parsed["wsnormis"]
        algo_dict["WSNormIS"] = WSNormIS(0.9f0, tdupdate)
        sample_dict["WSNormIS"] = "ER"
        value_type_dict["WSNormIS"] = "State"
    end
    if parsed["wsavgnormis"]
        algo_dict["WSAvgNormIS"] = WSNormIS()
        sample_dict["WSAvgNormIS"] = "ER"
        value_type_dict["WSAvgNormIS"] = "State"
    end
    if parsed["ir"]
        algo_dict["IR"] = tdupdate
        sample_dict["IR"] = "IR"
        value_type_dict["IR"] = "State"
    end
    if parsed["bcir"]
        algo_dict["BCIR"] = tdupdate
        sample_dict["BCIR"] = "IR"
        value_type_dict["BCIR"] = "State"
    end

    if parsed["wisbatch"]
        algo_dict["WISBatch"] = WISBatchTD(tdupdate)
        sample_dict["WISBatch"] = "ER"
        value_type_dict["WISBatch"] = "State"
    end
    if parsed["wisbuffer"]
        algo_dict["WISBuffer"] = tdupdate
        sample_dict["WISBuffer"] = "ER"
        value_type_dict["WISBuffer"] = "State"
    end

    if parsed["wisoptimal"]
        algo_dict["WISOptimal"] = WISBatchTD(tdupdate)
        sample_dict["WISOptimal"] = "Optimal"
        value_type_dict["WISOptimal"] = "State"
    end

    if parsed["vtrace"]
        for clip_value in parsed["clip_value"]
            algo_dict["VTrace_$(clip_value)"] = VTrace(clip_value, tdupdate)
            sample_dict["VTrace_$(clip_value)"] = "ER"
            value_type_dict["VTrace_$(clip_value)"] = "State"
        end
        for clip_value_perc in parsed["clip_value_perc"]
            algo_dict["VTrace_$(clip_value_perc*max_is)"] = VTrace(clip_value_perc*max_is, tdupdate)
            sample_dict["VTrace_$(clip_value_perc*max_is)"] = "ER"
            value_type_dict["VTrace_$(clip_value_perc*max_is)"] = "State"
        end
    end

    if parsed["sarsa"]
        algo_dict["Sarsa"] = BatchSarsa()
        sample_dict["Sarsa"] = "ER"
        value_type_dict["Sarsa"] = "StateAction"
    end

    if parsed["expectedsarsa"]
        algo_dict["ExpectedSarsa"] = BatchExpectedSarsa()
        sample_dict["ExpectedSarsa"] = "ER"
        value_type_dict["ExpectedSarsa"] = "StateAction"
    end

    if parsed["wisrupam"]
        if length(algo_dict) != 0
            throw("WISRupam must be done independently")
        end
        for u_0 in parsed["init_u"]
            algo_dict["WISRupam_$(u_0)"] = Resampling.WISBatchTD_Rupam(u_0, parsed["normalize_eta"])
            sample_dict["WISRupam_$(u_0)"] = "ER"
            value_type_dict["WISRupam_$(u_0)"] = "State"
        end
    end

    if parsed["res_wisrupam"]
        # if length(algo_dict) != 0
        #     throw("WISRupam must be done independently")
        # end
        for u_0 in parsed["init_u"]
            algo_dict["RES_WISRupam_$(u_0)"] = Resampling.WISBatchTD_Rupam(u_0, parsed["normalize_eta"])
            sample_dict["RES_WISRupam_$(u_0)"] = "IR"
            value_type_dict["RES_WISRupam_$(u_0)"] = "State"
        end
    end

    return algo_dict, sample_dict, value_type_dict
end

export MarkovChainUtil
include("exp_utils/MarkovChain.jl")

export FourRoomsUtil
include("exp_utils/FourRooms.jl")

export FourRoomsContUtil
include("exp_utils/FourRoomsCont.jl")

export FluxUtil
include("exp_utils/Flux.jl")

end


text:
using Documenter, SyntheticGrids

makedocs(
    modules = [SyntheticGrids],
    format = :html,
    pages = [
        "Home" => "index.md",
        "Types" => "Types.md",
        "Model" => "Model.md",
        "Functions" => [
		"Functions.md",
		"Private.md",
	],
    ],
    repo = "https://github.com/invenia/SyntheticGrids.jl/blob/{commit}{path}#L{line}",
    sitename = "SyntheticGrids.jl",
    authors = "Invenia Technical Computing",
    assets = ["assets/invenia.css"],
)

deploydocs(
    repo = "github.com/invenia/SyntheticGrids.jl",
    deps   = Deps.pip("python-markdown-math"),
    julia  = "0.6",
    target = "build",
    make = nothing,
)


text:
######################################################## model
F(z::BorderedArray,θ::AbstractVector,N::Int,M::Int) = F(z.u,(θ=θ,p=z.p),N::Int,M::Int)
function F(u::AbstractVector{T},parameters::NamedTuple,N::Int,M::Int) where T<:Number

	@unpack θ,p = parameters
	f = first(u)*first(p)*first(θ)
	F = zeros(typeof(f),length(u))

	F[1] = sin(p)^2 - ( θ[1]*sin(p)^2 + 1 )*u[1]
	for i ∈ 2:N
		F[i] = u[i-1] - u[i]
	end

	for i ∈ 2:M
		F[1+mod(i-1,N)] -= u[1+mod(i-1,N)]*θ[i]
	end

	return F
end

function scaling_backward(N::Int,M::Int)
	f(u,p) = F(u,p,N,M)

	X = StateSpace(N,-π:0.01:π,[0.0])
	θ = SizedVector{M}(ones(M))

	println("N = $N M = $M")
	@time ∇loss(f,θ,X)
	return @elapsed ∇loss(f,θ,X)
end

function scaling_forward(N::Int,M::Int)
	f(u,p) = F(u,p,N,M)

	X = StateSpace(N,-π:0.01:π,[0.0])
	θ = SizedVector{M}(ones(M))

	println("N = $N M = $M")
	@time loss(f,θ,X)
	return @elapsed loss(f,θ,X)
end

function scaling_continuation(N::Int,M::Int)
	f(u,p) = F(u,p,N,M)

	X = StateSpace(N,-π:0.01:π,[0.0])
	parameters = (θ=SizedVector{M}(ones(M)),p=minimum(X.parameter))
	hyperparameters = getParameters(X)

	println("N = $N M = $M")
	@time deflationContinuation(f,X.roots,parameters,hyperparameters)
	return @elapsed deflationContinuation(f,X.roots,parameters,hyperparameters)
end


text:
module SurrealNumbers

import Base: convert, promote, promote_rule, hash, show, read, delete!, empty!,
    size, sign, round, floor, ceil, mod, trunc,
    isinteger, isinf, isnan, isfinite, isodd, iseven, iszero,
    <, <=, zero, one, ==, <, -, +, *, /,
    summarysize

export Surreal, SurrealFinite, SurrealShort, SurrealDyadic, SurrealDAGstats
export SurrealZero, SurrealOne, SurrealMinusOne, SurrealTwo, SurrealMinusTwo, SurrealThree, ϕ

# eventually these will be removed, but they are convenient for testing
export ExistingSurreals, ExistingConversions, ExistingCanonicals, ExistingProducts, ExistingSums, ExistingNegations, Count
export clearcache

# export SurrealAlt2
# export ExistingSurrealsAlt2, ExistingProductsAlt2, ExistingSumsAlt2

export ∅, ϕ, <=, zero, one, ==, ≅, ≇, <, -, +, *, ≀, ↑, ↓, dag_add,
    dali, convert, overbarString, promote, hash,
    pf, pff, spf, show, surreal2node, surreal2dot, surreal2dag, surreal2tex, read, expand,
    generation, canonicalise, iscanonical, parents, isancestor, ≺, ≻, ⪯, ⪰,
    unique2!,
    sign, round, floor, ceil, trunc, mod,
    isinteger, isinf, isnan, isfinite, isodd, iseven, isdivisible,
    dag_stats, nodes, edges, paths, tree_nodes, breadth, width, surrealDAG,
    uniqueness, uniqueness_max, uniqueness_failure

abstract type Surreal <: Real end

@static if VERSION < v"0.7.0"
    const stdout = STDOUT
end

# include("Dyadic.jl")
include("SurrealFinite.jl")
# include("SurrealTrans.jl")
# include("SurrealAlt2.jl")

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

let nextidx = 0
    global nextproc
    function nextproc()
        p = -1
        if p == -1
            p = workers()[(nextidx % nworkers()) + 1]
            nextidx += 1
        end
        p
    end
end

spawnat(p, thunk) = sync_add(remotecall(thunk, p))

spawn_somewhere(thunk) = spawnat(nextproc(),thunk)

"""
    @spawn

Create a closure around an expression and run it on an automatically-chosen process,
returning a [`Future`](@ref) to the result.

# Examples
```julia-repl
julia> addprocs(3);

julia> f = @spawn myid()
Future(2, 1, 5, Nullable{Any}())

julia> fetch(f)
2

julia> f = @spawn myid()
Future(3, 1, 7, Nullable{Any}())

julia> fetch(f)
3
```
"""
macro spawn(expr)
    thunk = esc(:(()->($expr)))
    :(spawn_somewhere($thunk))
end

"""
    @spawnat

Create a closure around an expression and run the closure
asynchronously on process `p`. Returns a [`Future`](@ref) to the result.
Accepts two arguments, `p` and an expression.

# Examples
```julia-repl
julia> addprocs(1);

julia> f = @spawnat 2 myid()
Future(2, 1, 3, Nullable{Any}())

julia> fetch(f)
2
```
"""
macro spawnat(p, expr)
    thunk = esc(:(()->($expr)))
    :(spawnat($(esc(p)), $thunk))
end

"""
    @fetch

Equivalent to `fetch(@spawn expr)`.
See [`fetch`](@ref) and [`@spawn`](@ref).

# Examples
```julia-repl
julia> addprocs(3);

julia> @fetch myid()
2

julia> @fetch myid()
3

julia> @fetch myid()
4

julia> @fetch myid()
2
```
"""
macro fetch(expr)
    thunk = esc(:(()->($expr)))
    :(remotecall_fetch($thunk, nextproc()))
end

"""
    @fetchfrom

Equivalent to `fetch(@spawnat p expr)`.
See [`fetch`](@ref) and [`@spawnat`](@ref).

# Examples
```julia-repl
julia> addprocs(3);

julia> @fetchfrom 2 myid()
2

julia> @fetchfrom 4 myid()
4
```
"""
macro fetchfrom(p, expr)
    thunk = esc(:(()->($expr)))
    :(remotecall_fetch($thunk, $(esc(p))))
end

# extract a list of modules to import from an expression
extract_imports!(imports, x) = imports
function extract_imports!(imports, ex::Expr)
    if Meta.isexpr(ex, (:import, :using))
        return push!(imports, ex.args[1])
    elseif Meta.isexpr(ex, :let)
        return extract_imports!(imports, ex.args[1])
    elseif Meta.isexpr(ex, (:toplevel, :block))
        for i in eachindex(ex.args)
            extract_imports!(imports, ex.args[i])
        end
    end
    return imports
end
extract_imports(x) = extract_imports!(Symbol[], x)

"""
    @everywhere [procs()] expr

Execute an expression under `Main` on all `procs`.
Errors on any of the processes are collected into a
`CompositeException` and thrown. For example:

    @everywhere bar = 1

will define `Main.bar` on all processes.

Unlike [`@spawn`](@ref) and [`@spawnat`](@ref),
`@everywhere` does not capture any local variables.
Instead, local variables can be broadcast using interpolation:

    foo = 1
    @everywhere bar = \$foo

The optional argument `procs` allows specifying a subset of all
processes to have execute the expression.

Equivalent to calling `remotecall_eval(Main, procs, expr)`.
"""
macro everywhere(ex)
    return :(@everywhere $procs() $(esc(ex))) # interpolation needs to work around hygiene bugs (#22307)
end

macro everywhere(procs, ex)
    imps = [Expr(:import, m) for m in extract_imports(ex)]
    return quote
        $(isempty(imps) ? nothing : Expr(:toplevel, imps...)) # run imports locally first
        let ex = $(Expr(:quote, ex)), procs = $(esc(procs))
            remotecall_eval(Main, procs, ex)
        end
    end
end

"""
    remotecall_eval(m::Module, procs, expression)

Execute an expression under module `m` on the processes
specified in `procs`.
Errors on any of the processes are collected into a
`CompositeException` and thrown.

See also `@everywhere`.
"""
function remotecall_eval(m::Module, procs, ex)
    @sync begin
        run_locally = 0
        for pid in procs
            if pid == myid()
                run_locally += 1
            else
                @async remotecall_wait(Core.eval, pid, m, ex)
            end
        end
        yield() # ensure that the remotecall_fetch have had a chance to start

        # execute locally last as we do not want local execution to block serialization
        # of the request to remote nodes.
        for _ in 1:run_locally
            @async Core.eval(m, ex)
        end
    end
    nothing
end

# optimized version of remotecall_eval for a single pid
# and which also fetches the return value
function remotecall_eval(m::Module, pid::Int, ex)
    return remotecall_fetch(Core.eval, pid, m, ex)
end

# Statically split range [1,N] into equal sized chunks for np processors
function splitrange(N::Int, np::Int)
    each = div(N,np)
    extras = rem(N,np)
    nchunks = each > 0 ? np : extras
    chunks = Vector{UnitRange{Int}}(nchunks)
    lo = 1
    for i in 1:nchunks
        hi = lo + each - 1
        if extras > 0
            hi += 1
            extras -= 1
        end
        chunks[i] = lo:hi
        lo = hi+1
    end
    return chunks
end

function preduce(reducer, f, R)
    N = length(R)
    chunks = splitrange(N, nworkers())
    all_w = workers()[1:length(chunks)]

    w_exec = Task[]
    for (idx,pid) in enumerate(all_w)
        t = Task(()->remotecall_fetch(f, pid, reducer, R, first(chunks[idx]), last(chunks[idx])))
        schedule(t)
        push!(w_exec, t)
    end
    reduce(reducer, [wait(t) for t in w_exec])
end

function pfor(f, R)
    [@spawn f(R, first(c), last(c)) for c in splitrange(length(R), nworkers())]
end

function make_preduce_body(var, body)
    quote
        function (reducer, R, lo::Int, hi::Int)
            $(esc(var)) = R[lo]
            ac = $(esc(body))
            if lo != hi
                for $(esc(var)) in R[(lo+1):hi]
                    ac = reducer(ac, $(esc(body)))
                end
            end
            ac
        end
    end
end

function make_pfor_body(var, body)
    quote
        function (R, lo::Int, hi::Int)
            for $(esc(var)) in R[lo:hi]
                $(esc(body))
            end
        end
    end
end

"""
    @parallel

A parallel for loop of the form :

    @parallel [reducer] for var = range
        body
    end

The specified range is partitioned and locally executed across all workers. In case an
optional reducer function is specified, `@parallel` performs local reductions on each worker
with a final reduction on the calling process.

Note that without a reducer function, `@parallel` executes asynchronously, i.e. it spawns
independent tasks on all available workers and returns immediately without waiting for
completion. To wait for completion, prefix the call with [`@sync`](@ref), like :

    @sync @parallel for var = range
        body
    end
"""
macro parallel(args...)
    na = length(args)
    if na==1
        loop = args[1]
    elseif na==2
        reducer = args[1]
        loop = args[2]
    else
        throw(ArgumentError("wrong number of arguments to @parallel"))
    end
    if !isa(loop,Expr) || loop.head !== :for
        error("malformed @parallel loop")
    end
    var = loop.args[1].args[1]
    r = loop.args[1].args[2]
    body = loop.args[2]
    if na==1
        thecall = :(pfor($(make_pfor_body(var, body)), $(esc(r))))
    else
        thecall = :(preduce($(esc(reducer)), $(make_preduce_body(var, body)), $(esc(r))))
    end
    thecall
end


text:
using GridPotentialFlow
using Test
using Literate

const GROUP = get(ENV, "GROUP", "All")

ENV["GKSwstype"] = "nul"

notebookdir = "../examples"
docdir = "../docs/src/manual"
litdir = "./literate"
testdir = "./"
solverdir = "./solver"

for (root, dirs, files) in walkdir(litdir)
  for file in files
    if endswith(file,".jl")
      (GROUP == "All" || GROUP == "Notebooks") && Literate.notebook(joinpath(root, file),notebookdir)
      (GROUP == "All" || GROUP == "Documentation") && Literate.markdown(joinpath(root, file),docdir)
      (GROUP == "All" || GROUP == "Scripts") && (Literate.script(joinpath(root, file),testdir); include(file))
    end
  end
end


text:
import ..Traits
using ..Traits: Trait
export Err,
       LogEnergy,
       Mem,
       Intervene,
       Rng,
       Scope,
       Cond,
       Propose,
       RandMutate,
       IgnoreCondition,
       Seen,
       Solve

# # Primitive Traits
struct Err end
struct LogEnergy end
struct Mem end
struct Intervene end
struct Rng end
struct Scope end
struct Cond end
struct Propose end
struct RandMutate end
struct IgnoreCondition end
struct Seen end
struct Solve end

function symtotrait(x::Symbol)
  if x == :err
    Err
  elseif x == :logenergy
    LogEnergy
  elseif x == :mem
    Mem
  elseif x == :intervene
    Intervene
  elseif x == :rng
    Rng
  elseif x == :scope
    Scope
  elseif x == :condition
    Cond
  elseif x == :ignorecondition
    IgnoreCondition
  elseif x == :propose
    Propose
  elseif x == :randmutate
    RandMutate
  elseif x == :seen
    Seen
  elseif x == :solve
    Solve
  else
    error("Unknown trait: $x")
  end
end

@generated function Traits.traits(k::Tags{K, V}) where {K, V}
  traits_ = map(symtotrait, K)
  Trait{Union{traits_...}}()
end

@generated function Traits.traits(k::Type{Tags{K, V}}) where {K, V}
  traits_ = map(symtotrait, K)
  Trait{Union{traits_...}}()
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Oceananigans.Solvers
using Oceananigans.Operators
using Oceananigans.Architectures
using Oceananigans.Grids: with_halo
using Oceananigans.Fields: Field, ZReducedField

import Oceananigans.Solvers: solve!
import Oceananigans.Architectures: architecture

"""
    struct PCGImplicitFreeSurfaceSolver{V, S, R}

The preconditioned conjugate gradient iterative implicit free-surface solver.

$(TYPEDFIELDS)
"""
struct PCGImplicitFreeSurfaceSolver{V, S, R}
    "The vertically-integrated lateral areas"
    vertically_integrated_lateral_areas :: V
    "The preconditioned conjugate gradient solver"
    preconditioned_conjugate_gradient_solver :: S
    "The right hand side of the free surface evolution equation"
    right_hand_side :: R
end

architecture(solver::PCGImplicitFreeSurfaceSolver) =
    architecture(solver.preconditioned_conjugate_gradient_solver)

"""
    PCGImplicitFreeSurfaceSolver(grid, settings)

Return a solver based on a preconditioned conjugate gradient method for the elliptic equation

```math
[∇ ⋅ H ∇ - Az / (g Δt²)] ηⁿ⁺¹ = (∇ʰ ⋅ Q★ - Az ηⁿ / Δt) / (g Δt)
```

representing an implicit time discretization of the linear free surface evolution equation
for a fluid with variable depth `H`, horizontal areas `Az`, barotropic volume flux `Q★`, time
step `Δt`, gravitational acceleration `g`, and free surface at time-step `n` `ηⁿ`.
"""
function PCGImplicitFreeSurfaceSolver(grid::AbstractGrid, gravitational_acceleration::Number, settings)
    # Initialize vertically integrated lateral face areas
    ∫ᶻ_Axᶠᶜᶜ = Field{Face, Center, Nothing}(with_halo((3, 3, 1), grid))
    ∫ᶻ_Ayᶜᶠᶜ = Field{Center, Face, Nothing}(with_halo((3, 3, 1), grid))

    vertically_integrated_lateral_areas = (xᶠᶜᶜ = ∫ᶻ_Axᶠᶜᶜ, yᶜᶠᶜ = ∫ᶻ_Ayᶜᶠᶜ)

    compute_vertically_integrated_lateral_areas!(vertically_integrated_lateral_areas)

    right_hand_side = Field{Center, Center, Nothing}(grid)

    # Set maximum iterations to Nx * Ny if not set
    settings = Dict{Symbol, Any}(settings)
    maximum_iterations = get(settings, :maximum_iterations, grid.Nx * grid.Ny)
    settings[:maximum_iterations] = maximum_iterations

    # Set preconditioner to default preconditioner if not specified
    preconditioner = get(settings, :preconditioner_method, implicit_free_surface_precondition!)
    settings[:preconditioner_method] = preconditioner

    solver = PreconditionedConjugateGradientSolver(implicit_free_surface_linear_operation!;
                                                   template_field = right_hand_side,
                                                   settings...)

    return PCGImplicitFreeSurfaceSolver(vertically_integrated_lateral_areas, solver, right_hand_side)
end

build_implicit_step_solver(::Val{:PreconditionedConjugateGradient}, grid, gravitational_acceleration, settings) =
    PCGImplicitFreeSurfaceSolver(grid, gravitational_acceleration, settings)

#####
##### Solve...
#####

function solve!(η, implicit_free_surface_solver::PCGImplicitFreeSurfaceSolver, rhs, g, Δt)
    #=
    # Somehow take an explicit step first to produce an improved initial guess
    # for η for the iterative solver.
    event = explicit_ab2_step_free_surface!(free_surface, model, Δt, χ)
    wait(device(model.architecture), event)
    =#

    ∫ᶻA = implicit_free_surface_solver.vertically_integrated_lateral_areas
    solver = implicit_free_surface_solver.preconditioned_conjugate_gradient_solver

    # solve!(x, solver, b, args...) solves A*x = b for x.
    solve!(η, solver, rhs, ∫ᶻA.xᶠᶜᶜ, ∫ᶻA.yᶜᶠᶜ, g, Δt)

    return nothing
end

function compute_implicit_free_surface_right_hand_side!(rhs,
                                                        implicit_solver::PCGImplicitFreeSurfaceSolver,
                                                        g, Δt, ∫ᶻQ, η)

    solver = implicit_solver.preconditioned_conjugate_gradient_solver
    arch = architecture(solver)
    grid = solver.grid

    event = launch!(arch, grid, :xy,
                    implicit_free_surface_right_hand_side!,
                    rhs, grid, g, Δt, ∫ᶻQ, η,
		            dependencies = device_event(arch))

    return event
end

""" Compute the divergence of fluxes Qu and Qv. """
@inline flux_div_xyᶜᶜᶜ(i, j, k, grid, Qu, Qv) = δxᶜᵃᵃ(i, j, k, grid, Qu) + δyᵃᶜᵃ(i, j, k, grid, Qv)

@kernel function implicit_free_surface_right_hand_side!(rhs, grid, g, Δt, ∫ᶻQ, η)
    i, j = @index(Global, NTuple)
    Az = Azᶜᶜᶜ(i, j, 1, grid)
    δ_Q = flux_div_xyᶜᶜᶜ(i, j, 1, grid, ∫ᶻQ.u, ∫ᶻQ.v)
    @inbounds rhs[i, j, 1] = (δ_Q - Az * η[i, j, 1] / Δt) / (g * Δt)
end

"""
Returns `L(ηⁿ)`, where `ηⁿ` is the free surface displacement at time step `n`
and `L` is the linear operator that arises
in an implicit time step for the free surface displacement `η`.

(See the docs section on implicit time stepping.)
"""
function implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)
    grid = L_ηⁿ⁺¹.grid
    arch = architecture(L_ηⁿ⁺¹)

    fill_halo_regions!(ηⁿ⁺¹, arch)

    event = launch!(arch, grid, :xy, _implicit_free_surface_linear_operation!,
                    L_ηⁿ⁺¹, grid,  ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    fill_halo_regions!(L_ηⁿ⁺¹, arch)

    return nothing
end

# Kernels that act on vertically integrated / surface quantities
@inline ∫ᶻ_Ax_∂x_ηᶠᶜᶜ(i, j, k, grid, ∫ᶻ_Axᶠᶜᶜ, η) = @inbounds ∫ᶻ_Axᶠᶜᶜ[i, j, k] * ∂xᶠᶜᶜ(i, j, k, grid, η)
@inline ∫ᶻ_Ay_∂y_ηᶜᶠᶜ(i, j, k, grid, ∫ᶻ_Ayᶜᶠᶜ, η) = @inbounds ∫ᶻ_Ayᶜᶠᶜ[i, j, k] * ∂yᶜᶠᶜ(i, j, k, grid, η)

"""
Compute the horizontal divergence of vertically-uniform quantity using
vertically-integrated face areas `∫ᶻ_Axᶠᶜᶜ` and `∫ᶻ_Ayᶜᶠᶜ`.
"""
@inline Az_∇h²ᶜᶜᶜ(i, j, k, grid, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, η::ZReducedField) =
    (δxᶜᵃᵃ(i, j, k, grid, ∫ᶻ_Ax_∂x_ηᶠᶜᶜ, ∫ᶻ_Axᶠᶜᶜ, η) +
     δyᵃᶜᵃ(i, j, k, grid, ∫ᶻ_Ay_∂y_ηᶜᶠᶜ, ∫ᶻ_Ayᶜᶠᶜ, η))

"""
    _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)

Return the left side of the "implicit η equation"

```math
(∇ʰ⋅H∇ʰ - 1/gΔt²) ηⁿ⁺¹ = 1 / (gΔt) ∇ʰ ⋅ Q★ - 1 / (gΔt²) ηⁿ
----------------------
        ≡ L_ηⁿ⁺¹
```

which is derived from the discretely summed barotropic mass conservation equation,
and arranged in a symmetric form by multiplying by horizontal areas Az:

```
δⁱÂʷ∂ˣηⁿ⁺¹ + δʲÂˢ∂ʸηⁿ⁺¹ - 1/gΔt² Az ηⁿ⁺¹ = 1 / (gΔt) (δⁱÂʷu̅ˢᵗᵃʳ + δʲÂˢv̅ˢᵗᵃʳ) - 1 / gΔt² Az ηⁿ
```

where  ̂ indicates a vertical integral, and
       ̅ indicates a vertical average
"""
@kernel function _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)
    i, j = @index(Global, NTuple)
    Az = Azᶜᶜᶜ(i, j, 1, grid)
    @inbounds L_ηⁿ⁺¹[i, j, 1] = Az_∇h²ᶜᶜᶜ(i, j, 1, grid, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, ηⁿ⁺¹) - Az * ηⁿ⁺¹[i, j, 1] / (g * Δt^2)
end

"""
    _implicit_free_surface_precondition!(P_r, grid, r, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)

Return the diagonally dominant inverse preconditioner applied to the residuals consistently with
 `M = D⁻¹(I - (A - D)D⁻¹) ≈ A⁻¹` where `I` is the Identity matrix, D is the matrix
containing the diagonal of A, and A is the linear operator applied to η

```math
P_r = M * r
```
which expanded in components is
```math
P_rᵢⱼ = rᵢⱼ / Acᵢⱼ - 1 / Acᵢⱼ ( Ax⁻ / Acᵢ₋₁ rᵢ₋₁ⱼ + Ax⁺ / Acᵢ₊₁ rᵢ₊₁ⱼ + Ay⁻ / Acⱼ₋₁ rᵢⱼ₋₁+ Ay⁺ / Acⱼ₊₁ rᵢⱼ₊₁ )
```

where `Ac`, `Ax⁻`, `Ax⁺`, `Ay⁻` and `Ay⁺` are the coefficients of
`ηᵢⱼ`, `ηᵢ₋₁ⱼ`, `ηᵢ₊₁ⱼ`, `ηᵢⱼ₋₁` and `ηᵢⱼ₊₁` in `_implicit_free_surface_linear_operation!`
"""
function implicit_free_surface_precondition!(P_r, r, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)
    grid = ∫ᶻ_Axᶠᶜᶜ.grid
    arch = architecture(P_r)

    fill_halo_regions!(r, arch)

    event = launch!(arch, grid, :xy, _implicit_free_surface_precondition!,
                    P_r, grid, r, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    fill_halo_regions!(P_r, arch)

    return nothing
end

# Kernels that calculate coefficients for the preconditioner
@inline Ax⁻(i, j, grid, ax) = @inbounds   ax[i, j, 1] / Δxᶠᶜᶜ(i, j, 1, grid)
@inline Ay⁻(i, j, grid, ay) = @inbounds   ay[i, j, 1] / Δyᶜᶠᶜ(i, j, 1, grid)
@inline Ax⁺(i, j, grid, ax) = @inbounds ax[i+1, j, 1] / Δxᶠᶜᶜ(i+1, j, 1, grid)
@inline Ay⁺(i, j, grid, ay) = @inbounds ay[i, j+1, 1] / Δyᶜᶠᶜ(i, j+1, 1, grid)

@inline Ac(i, j, grid, g, Δt, ax, ay) = - ( Ax⁻(i, j, grid, ax)
                                          + Ax⁺(i, j, grid, ax)
                                          + Ay⁻(i, j, grid, ay)
                                          + Ay⁺(i, j, grid, ay)
                                          + Azᶜᶜᶜ(i, j, 1, grid) / (g * Δt^2) )

@inline heuristic_inverse_times_residuals(i, j, r, grid, g, Δt, ax, ay) = @inbounds 1 / Ac(i, j, grid, g, Δt, ax, ay) * ( r[i, j, 1] -
                            Ax⁻(i, j, grid, ax) / Ac(i-1, j, grid, g, Δt, ax, ay) * r[i-1, j, 1] -
                            Ax⁺(i, j, grid, ax) / Ac(i+1, j, grid, g, Δt, ax, ay) * r[i+1, j, 1] -
                            Ay⁻(i, j, grid, ay) / Ac(i, j-1, grid, g, Δt, ax, ay) * r[i, j-1, 1] -
                            Ay⁺(i, j, grid, ay) / Ac(i, j+1, grid, g, Δt, ax, ay) * r[i, j+1, 1] )

@kernel function _implicit_free_surface_precondition!(P_r, grid, r, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)
    i, j = @index(Global, NTuple)
    @inbounds P_r[i, j, 1] = heuristic_inverse_times_residuals(i, j, r, grid, g, Δt, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ)
end


text:
# Loosely based on MiniZinc's fzn_lex_lesseq_float for MIP:
# https://github.com/MiniZinc/libminizinc/blob/master/share/minizinc/linear/fzn_lex_lesseq_float.mzn
# Major difference: that reformulation is only for two vectors,
# this implementation is for an arbitrary number of vectors (there is no
# dedicated set for just two vectors).

"""
Bridges `CP.LexicographicallyGreaterThan` to indicators.
"""
struct LexicographicallyGreaterThan2IndicatorBridge{T} <: MOIBC.AbstractBridge
    vars_eq::Matrix{MOI.VariableIndex}
    vars_gt::Matrix{MOI.VariableIndex}
    vars_eq_bin::Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    vars_gt_bin::Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    cons_one_gt::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
    cons_move::Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_indic_eq::Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}}}
    cons_indic_gt::Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.GreaterThan{T}}}}
end

function MOIBC.bridge_constraint(
    ::Type{LexicographicallyGreaterThan2IndicatorBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.LexicographicallyGreaterThan,
) where {T}
    return MOIBC.bridge_constraint(
        LexicographicallyGreaterThan2IndicatorBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{LexicographicallyGreaterThan2IndicatorBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.LexicographicallyGreaterThan,
) where {T}
    @assert s.row_dim >= 1
    @assert s.column_dim >= 2

    vars_eq = Matrix{MOI.VariableIndex}(undef, s.column_dim - 1, s.row_dim)
    vars_gt = Matrix{MOI.VariableIndex}(undef, s.column_dim - 1, s.row_dim)
    vars_eq_bin = Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}(undef, s.column_dim - 1, s.row_dim)
    vars_gt_bin = Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}(undef, s.column_dim - 1, s.row_dim)
    cons_one_gt = Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}(undef, s.column_dim - 1)
    cons_move = Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}(undef, s.column_dim - 1, s.row_dim - 1)
    cons_indic_eq = Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}}}(undef, s.column_dim - 1, s.row_dim)
    cons_indic_gt = Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.GreaterThan{T}}}}(undef, s.column_dim - 1, s.row_dim)

    f_scalars = MOIU.scalarize(f)
    f_matrix = reshape(f_scalars, s.column_dim, s.row_dim)
    # First index: column; second index: row. Rationale: columns are sorted
    # lexicographically, not rows.

    # Add the constraint between the columns i and i+1.
    for i in 1:(s.column_dim - 1)
        vars_eq[i, :], vars_eq_bin[i, :] = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:s.row_dim])
        vars_gt[i, :], vars_gt_bin[i, :] = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:s.row_dim])

        cons_one_gt[i] = MOI.add_constraint(
            model,
            sum(one(T) .* vars_gt[i, :]),
            MOI.LessThan(one(T)),
        )

        cons_move[i, :] = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
            MOI.add_constraint(
                model,
                one(T) * vars_eq[i, j - 1] - one(T) * vars_eq[i, j] - one(T) * vars_gt[i, j],
                MOI.EqualTo(zero(T)),
            )
            for j in 2:s.row_dim
        ]

        cons_indic_eq[i, :] = MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}}[
            MOI.add_constraint(
                model,
                MOIU.vectorize(
                    [
                        one(T) * vars_eq[i, j],
                        f_matrix[i, j] - f_matrix[i + 1, j],
                    ]
                ),
                MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(zero(T))),
            )
            for j in 1:s.row_dim
        ]

        cons_indic_gt[i, :] = MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.GreaterThan{T}}}[
            MOI.add_constraint(
                model,
                MOIU.vectorize(
                    [
                        one(T) * vars_gt[i, j],
                        f_matrix[i, j] - f_matrix[i + 1, j],
                    ]
                ),
                MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(zero(T))),
            )
            for j in 1:s.row_dim
        ]
    end

    return LexicographicallyGreaterThan2IndicatorBridge(vars_eq, vars_gt, vars_eq_bin, vars_gt_bin, cons_one_gt, cons_move, cons_indic_eq, cons_indic_gt)
end

function MOI.supports_constraint(
    ::Type{LexicographicallyGreaterThan2IndicatorBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.LexicographicallyGreaterThan},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{LexicographicallyGreaterThan2IndicatorBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{LexicographicallyGreaterThan2IndicatorBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
        (MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}),
        (MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.GreaterThan{T}}),
    ]
end

function MOI.get(b::LexicographicallyGreaterThan2IndicatorBridge{T}, ::MOI.NumberOfVariables) where {T}
    return length(b.vars_eq) + length(b.vars_gt)
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return length(b.vars_eq_bin) + length(b.vars_gt_bin)
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T}
    return length(b.cons_one_gt)
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return length(b.cons_move)
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}},
    },
) where {T}
    return length(b.cons_indic_eq)
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.GreaterThan{T}},
    },
) where {T}
    return length(b.cons_indic_gt)
end

function MOI.get(b::LexicographicallyGreaterThan2IndicatorBridge{T}, ::MOI.ListOfVariableIndices) where {T}
    return vcat(vec(b.vars_eq), vec(b.vars_gt))
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return vcat(vec(b.vars_eq_bin), vec(b.vars_gt_bin))
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T}
    return copy(b.cons_one_gt)
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return vec(b.cons_move)
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}},
    },
) where {T}
    return vec(b.cons_indic_eq)
end

function MOI.get(
    b::LexicographicallyGreaterThan2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.GreaterThan{T}},
    },
) where {T}
    return vec(b.cons_indic_gt)
end


text:
# Utils for foveation and spacial querying of the database

export
  whosNear2D,
  whosNear3D,
  foveateQueryToPoint

  """
      whosNear2D(cg::CloudGraph, session::AbstractString; x,y,yaw,dist,angle )

  Find vertices near the point specified and return dictionary of symbol to Neo4j ID pairs.
  """
  function whosNear2D(cg::CloudGraph, session::AbstractString;
        x=nothing,
        y=nothing,
        yaw=nothing,
        dist::Float64=3.0,
        angle::Float64=pi/12.0  )
    #
    if x==nothing && y==nothing && yaw==nothing
      error("please give some info on where you want to run the spacial query, xyYaw.")
    end
    # Build the query
    query = "match (n:$(session)) where exists(n.MAP_est) "
    query = x==nothing ? query : query*"and abs(n.MAP_est[0] - $(x)) < $(dist) "
    query = y==nothing ? query : query*"and abs(n.MAP_est[1] - $(y)) < $(dist) "
    query = yaw==nothing ? query : query*"and abs(n.MAP_est[2] - $(yaw)) < $(angle) "
    query = query*"return n.label, id(n)"

    cph, = executeQuery(cg, query)

    symneo = Dict{Symbol, Int}()
    for data in cph.results[1]["data"]
      symneo[Symbol(data["row"][1])] = data["row"][2]
    end
    # neoid = cph.results[1]["data"][1]["row"][1]

    symneo
  end

  """
      whosNear3D(cg::CloudGraph, session::AbstractString; x,y,z,roll, pitch,yaw,dist,angle )

  Find vertices near the point specified and return dictionary of symbol to Neo4j ID pairs.
  """
  function whosNear3D(cg::CloudGraph, session::AbstractString;
        x=nothing,
        y=nothing,
        z=nothing,
        roll=nothing,
        pitch=nothing,
        yaw=nothing,
        dist::Float64=3.0,
        angle::Float64=pi/12.0  )
    #
    if x==nothing && y==nothing && z==nothing && roll==nothing && pitch==nothing && yaw==nothing
      error("please give some info on where you want to run the spacial query, xyzRPY.")
    end
    # Build the query
    query = "match (n:$(session)) where exists(n.MAP_est) "
    query = x==nothing ? query : query*"and abs(n.MAP_est[0] - $(x)) < $(dist) "
    query = y==nothing ? query : query*"and abs(n.MAP_est[1] - $(y)) < $(dist) "
    query = z==nothing ? query : query*"and abs(n.MAP_est[2] - $(z)) < $(dist) "
    query = roll==nothing ? query : query*"and abs(n.MAP_est[3] - $(roll)) < $(angle) "
    query = pitch==nothing ? query : query*"and abs(n.MAP_est[4] - $(pitch)) < $(angle) "
    query = yaw==nothing ? query : query*"and abs(n.MAP_est[5] - $(yaw)) < $(angle) "
    query = query*"return n.label, id(n)"

    cph, = executeQuery(cg, query)

    symneo = Dict{Symbol, Int}()
    for data in cph.results[1]["data"]
      symneo[Symbol(data["row"][1])] = data["row"][2]
    end

    symneo
  end

function foveateQueryToPoint{T <: AbstractString}(cg::CloudGraph,
        sessions::Vector{T};
        point::Vector{Float64}=[0.0;0.0],
        minrange::Number=1.5, maxrange::Number=5, fovrad::Number=0.3  )
  #
  neoids, syms = Vector{Int}(), Vector{Symbol}()
  loadtx = transaction(cg.neo4j.connection)
  nms = length(sessions)
  query = "MATCH (n) WHERE ("
  cnt = 0
  for session in sessions
    query *= "(n:$(session)) "
    cnt+=1
    if cnt < nms
      query *= " or "
    else
      query *= ") and "
    end
  end
  query *=
  "// Region filtering
  	($(point[1])-n.MAP_est[0])*($(point[1])-n.MAP_est[0]) +
  	($(point[2])-n.MAP_est[1])*($(point[2])-n.MAP_est[1]) > $(minrange^2)
  AND
  	($(point[1])-n.MAP_est[0])*($(point[1])-n.MAP_est[0]) +
  	($(point[2])-n.MAP_est[1])*($(point[2])-n.MAP_est[1]) < $(maxrange^2)
  AND
  // Frustum cutoff filtering
  	(abs(atan2( ($(point[2])-n.MAP_est[1]),
  	($(point[1])-n.MAP_est[0])) - n.MAP_est[2] ) < $(fovrad)
  OR
  	abs(atan2( ($(point[2])-n.MAP_est[1]),
  	($(point[1])-n.MAP_est[0])) - n.MAP_est[2] + 6.28) < $(fovrad)
  OR
  	abs(atan2( ($(point[2])-n.MAP_est[1]),
  	($(point[1])-n.MAP_est[0])) - n.MAP_est[2] - 6.28) < $(fovrad))
  RETURN id(n), n.label"
  cph = loadtx(query, submit=true)
  for res in cph.results[1]["data"]
    push!(neoids, res["row"][1]  )
    push!(syms, Symbol(res["row"][2])  )
  end
  return neoids, syms
end

function foveateQueryToPoint(cg::CloudGraph,
        session::AbstractString;
        point::Vector{Float64}=[0.0;0.0],
        minrange::Number=1.5, maxrange::Number=5, fovrad::Number=0.3  )
  #
  foveateQueryToPoint(cg::CloudGraph,
          [session],
          point=point,
          minrange=minrange, maxrange=maxrange, fovrad=fovrad  )
end


text:
module Day05

using AdventOfCode2019

function day05(input::String = readInput(joinpath(@__DIR__, "..", "data", "day05.txt")))
    data = parse.(Int, split(input, ","))

    c1 = Channel{Int}(2)
    put!(c1, 1)
    out1 = AdventOfCode2019.IntCode.run_program!(copy(data), c1, nothing, nothing)
    close(c1)

    c2 = Channel{Int}(2)
    put!(c2, 5)
    out2 = AdventOfCode2019.IntCode.run_program!(copy(data), c2, nothing, nothing)
    close(c2)

    return [out1[end], out2[end]]
end

end # module


text:
# Saving and loading happens in a nested fashion:
# save(filename, mc) calls
#   save_mc(filename, mc) calls
#       save_model(filename, mc.model) calls
#           save_lattice(filename, mc.lattice) (TODO, not used yet)
#       save_measurements(filename, mc) calls
#           save_measurement(filename, measurement)
#
# loading follows the same structure
# > Each level (beyond the outermost save) also has an argument `entryname::String`.
# > Each save_something() should write a "VERSION" and a "type" to the file. The
#   former simplifies updates down the line, the latter allows for dispatch.
# > Each load_something should be dispatched on a type, e.g.
#   `load_model(data, ::Type(HubbardModelAttractive))`

# Loading:
# load(filename[, groups])
#   - loads from filename, maybe specific group
#   - jumps into _load at fitting part
# _load(data, ::Type{...})
#   - loads from Dict or JLD2.JLDFile or JLD2.Group (data) to some type (for dispatch)

# TODOs:
# - rename save_x to _save(filename, ::X, entryname) (mirror _load)
# - make _load less volatile to changes, i.e. change Type -> Val{Symbol}

"""
    save(filename, mc; overwrite=false, rename=true)

Saves the given MonteCarlo simulation `mc` to a JLD-file `filename`.

If `rename = true` the filename will be adjusted if it already exists. If
`overwrite = true` it will be overwritten. In this case a temporary backup
will be created. If neither are true an error will be thrown.
"""
function save(
        filename, mc::MonteCarloFlavor;
        overwrite = false, rename = true, compress = true,
        backend = endswith(filename, "jld2") ? JLD2 : JLD, kwargs...
    )
    # endswith(filename, ".jld") || (filename *= ".jld")

    # handle ranming and overwriting
    isfile(filename) && !overwrite && !rename && throw(ErrorException(
        "Cannot save because \"$filename\" already exists. Consider setting " *
        "`rename = true` to adjust the filename or `overwrite = true`" *
        " to overwrite the file."
    ))
    if isfile(filename) && !overwrite && rename
        filename = _generate_unique_filename(filename)
    end

    temp_filename = ""
    if isfile(filename) && overwrite
        parts = splitpath(filename)
        parts[end] = "." * parts[end]
        temp_filename = _generate_unique_filename(joinpath(parts...))
        mv(filename, temp_filename)
    end

    mode = isfile(filename) ? "r+" : "w"
    file = FileWrapper(
        backend.jldopen(filename, mode, compress=compress; kwargs...), filename
    )

    write(file, "VERSION", 1)
    save_mc(file, mc, "MC")
    save_rng(file)
    close(file.file)

    if overwrite && !isempty(temp_filename) && isfile(temp_filename)
        rm(temp_filename)
    end

    return filename
end

# Something like
# existing_file.jld -> existing_file_aJ3c.jld
function _generate_unique_filename(filename)
    isfile(filename) || return filename
    x = rand("0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz")
    s = "_$(Char(x))"
    parts = split(filename, '.')
    filename = join(parts[1:end-1], '.') * s
    while isfile(filename * '.' * parts[end])
        x = rand("0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz")
        s = string(Char(x))
        filename = filename * s
    end
    filename * '.' * parts[end]
end

to_tag(f::FileWrapper) = to_tag(f.file)
function to_tag(data::Union{JLDFile, JLD2.Group, Dict{String, Any}})
    haskey(data, "tag") && return Val(Symbol(data["tag"]))
    haskey(data, "type") && return to_tag(data["type"])
    error("Failed to get tag from $data")
end
to_tag(::Type{<: AbstractLattice}) = Val(:Generic)
to_tag(::Type{<: Model}) = Val(:Generic)

# get all files in directory recursively
function to_files(path_or_filename)
    if isfile(path_or_filename)
        return [path_or_filename]
    elseif isdir(path_or_filename)
        return vcat(to_files.(readdir(path_or_filename, join=true))...)
    else
        error("$path_or_filename is neither a valid directory nor file path.")
    end
end

function load(
        paths_or_filenames::Vector{String};
        prefix = "", postfix = r"jld|jld2", simplify = false, silent = false,
        parallel = true, on_error = e -> @error(exception = e)
    )
    # Normalize input to filepaths (recursively)
    files = String[]
    for path_or_file in paths_or_filenames
        _files = to_files(path_or_file)
        filter!(_files) do filepath
            _, filename = splitdir(filepath)
            startswith(filename, prefix) && endswith(filename, postfix)
        end
        append!(files, _files)
    end

    println(
        "Loading $(length(files)) Simulations",
        parallel && nprocs() > 1 ? " on $(nworkers()) workers" : ""
    )
    flush(stdout)

    # Might be worth shuffling files for more equal load times?
    mcs = ProgressMeter.@showprogress pmap(files, distributed = parallel, on_error = on_error) do f
        mc = load(f)
        simplify && simplify_measurements!(mc)
        mc
    end

    return mcs
end

"""
    load(filename[, groups...])
    load(path_or_collection; prefix = "", postfix = "jld|jld2", simplify = false)

Loads one or many MonteCarlo simulations from the given file path, directory or
`Vector` thereof.

If the given argument is directory or `Vector` of directories and file paths all
the directories will be expanded recursively. You can use `prefix` and `postfix`
to filter valid filenames. (These are only applied to the filename, not the
full path.)  If `simplify = true` a conversion of measurements to `ValueWrapper`
will be attempted. (This is useful to reduce the memory requirements.)
"""
function load(filename::String, groups::String...; kwargs...)
    if isfile(filename)
        data = if endswith(filename, "jld2")
            FileWrapper(JLD2.jldopen(filename, "r"), filename)
        else
            FileWrapper(JLD.load(filename), filename)
        end
        output = try
            if haskey(data, "MC") && !("MC" in groups)
                _load(data, "MC", groups...) else _load(data, groups...)
            end
        finally
            endswith(filename, "jld2") && close(data.file)
        end

        return output
    elseif isdir(filename)
        return load([filename]; kwargs...)
    else
        error("$filename is not a valid file or directory.")
    end
end
_load(data, g1::String, g2::String, gs::String...) = _load(data[g1], g2, gs...)
function _load(data, g::String)
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load $(data.path) version $(data["VERSION"])"))
    end

    haskey(data[g], "RNG") && load_rng!(data)

    _load(data[g], to_tag(data[g]))
end
_load(data) = _load(data, to_tag(data))

"""
    resume!(filename[; kwargs...])

Resumes a Monte Carlo simulation from a savefile generated by `run!(mc)`. Takes
the same keyword arguments as `run!`. Returns the simulation and the state
returned by `run!`.

See also: [`run!`](@ref)
"""
function resume!(filename; kwargs...)
    data = if endswith(filename, "jld2")
        FileWrapper(JLD2.jldopen(filename, "r"), filename)
    else
        FileWrapper(JLD.load(filename), filename)
    end

    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load $filename version $(data["VERSION"])"))
    end

    mc = _load(data["MC"], to_tag(data["MC"]))
    resume_init!(mc)
    load_rng!(data)
    endswith(filename, "jld2") && close(data)

    state = run!(mc; kwargs...)
    mc, state
end

function save_mc(
        filename::String, mc::MonteCarloFlavor, entryname::String="MC";
        backend = endswith(filename, "jld2") ? JLD2 : JLD,
        kwargs...
    )
    mode = isfile(filename) ? "r+" : "w"
    file = backend.jldopen(filename, mode; kwargs...)
    save_mc(file, mc, entryname)
    close(file)
    nothing
end
to_tag(::Type{<: Union{UnknownType, JLD2.UnknownType}}) = Val{:UNKNOWN}
function _load(data, ::Val{:UNKNOWN})
    @info "Failed to load (Unknowntype)"
    @info "Available fields: $(keys(data))"
    @info "You may be missing external packages."
    nothing
end

#     save_model(filename, model, entryname)
#
# Save (minimal) information necessary to reconstruct the given `model` in a
# jld-file `filename` under group `entryname`.
#
# By default the full model object is saved. When saving a simulation, the
# entryname defaults to `MC/Model`.
function save_model(
        filename::String, model, entryname::String;
        backend = endswith(filename, "jld2") ? JLD2 : JLD, kwargs...
    )
    mode = isfile(filename) ? "r+" : "w"
    file = backend.jldopen(filename, mode; kwargs...)
    save_model(file, model, entryname)
    close(file)
    nothing
end
function save_model(file::JLDFile, model, entryname::String)
    write(file, entryname * "/VERSION", 0)
    write(file, entryname * "/tag", "Generic")
    write(file, entryname * "/data", model)
    nothing
end

"""
    _load(data, ::Type{...})

Loads `data` where `data` is either a `JLD2.JLDFile`, `JLD2.Group` or a `Dict`.

The default `_load` will check that `data["VERSION"] == 0` and simply return
`data["data"]`. You may implement `_load(data, ::Type{<: MyType})` to add
specialized loading behavior.
"""
function _load(data, ::Val{:Generic}) where T
    data["VERSION"] == 0 || throw(ErrorException(
        "Version $(data["VERSION"]) incompatabile with default _load for $T."
    ))
    data["data"]
end

#     save_lattice(filename, lattice, entryname)
#
# Save (minimal) information necessary to reconstruct the given `lattice` in a
# jld-file `filename` under group `entryname`.
#
# By default the full lattice object is saved. When saving a simulation, the
# entryname defaults to `MC/Model/Lattice`.
function save_lattice(
        filename::String, lattice::AbstractLattice, entryname::String;
        backend = endswith(filename, "jld2") ? JLD2 : JLD, kwargs...
    )
    mode = isfile(filename) ? "r+" : "w"
    file = backend.jldopen(filename, mode; kwargs...)
    save_lattice(file, lattice, entryname)
    close(file)
    nothing
end
function save_lattice(file::JLDFile, lattice::AbstractLattice, entryname::String)
    write(file, entryname * "/VERSION", 0)
    write(file, entryname * "/tag", "Generic")
    write(file, entryname * "/data", lattice)
    nothing
end

const _GLOBAL_RNG = VERSION < v"1.3.0" ? Random.GLOBAL_RNG : Random.default_rng()

"""
    save_rng(filename [; rng = _GLOBAL_RNG, entryname = "RNG"])

Saves the current state of Julia's random generator (`Random.GLOBAL_RNG`) to the
given `filename`.
"""
function save_rng(
        filename::String; rng = _GLOBAL_RNG, entryname::String="RNG",
        backend = endswith(filename, "jld2") ? JLD2 : JLD, kwargs...
    )
    mode = isfile(filename) ? "r+" : "w"
    file = backend.jldopen(filename, mode; kwargs...)
    save_rng(file, rng=rng, entryname=entryname)
    close(file)
end
function save_rng(file::JLDFile; rng = _GLOBAL_RNG, entryname::String="RNG")
    try
        write(file, entryname, rng)
    catch e
        error("Error while saving RNG state: ", e)
    end
end

"""
    load_rng!(data[; rng = _GLOBAL_RNG, entryname = "RNG"])

Loads an RNG from a given `data` dictinary generated by `JLD.load()` to `rng`.
"""
function load_rng!(data; rng = _GLOBAL_RNG, entryname::String="RNG")
    try
        copy!(rng, data[entryname])
    catch e
        error("Error while restoring RNG state: ", e)
    end
end


text:
# importing
import Base.repr
import Base.in
import Base.(==)
import Base.(<)
import Base.(>)
import Base.isless
import Base.getindex
import Base.push!

# just a little type alias
typealias Handler Function

# The children should be sorted with static endpoints first
abstract Endpoint

const ROOT_TAG = "__pivot__"

# For debugging. Makes endpoints a little nicer looking in the terminal.
function repr(e::Endpoint, offset=0)
    val = (typeof(e) == StaticEndpoint) ?  "S" : "D"
    repeat("\t", offset) * "$(e.tag)-$val$offset-$(length(e.handlermap))-($((isempty(e.children))?"":"\n")" * join(map((ep) -> repr(ep, offset+1), e.children), "\n") * ")"
end

"""
StaticEndpoint
"""
type StaticEndpoint <: Endpoint
    tag
    children::Vector{Endpoint}
    handlermap::Dict{Verb, Handler}
    StaticEndpoint() = new(ROOT_TAG, [], Dict())
    StaticEndpoint(name) = new(name, [], Dict())
end

"""
DynamicEndpoint
matches any string and captures what it compares during a comparison check.
"""
type DynamicEndpoint <: Endpoint
    tag
    captured
    children::Vector{Endpoint}
    handlermap::Dict{Verb, Handler}
    DynamicEndpoint(name) = new(name, nothing, [], Dict())
end

(==)(tag::AbstractString, ep::StaticEndpoint) = tag == ep.tag
(==)(ep::StaticEndpoint, tag::AbstractString) = tag == ep.tag
(==)(tag::AbstractString, ep::DynamicEndpoint) = (ep.captured = tag; true)
(==)(ep::DynamicEndpoint, tag::AbstractString) = (ep.captured = tag; true)

"""
static endpoint should be less that a dynamic endpoint.
"""
isless(::StaticEndpoint, ::DynamicEndpoint) = true
isless(::DynamicEndpoint, ::StaticEndpoint) = false
isless(::StaticEndpoint, ::StaticEndpoint) = false
isless(::DynamicEndpoint, ::DynamicEndpoint) = false

"""
check if the tag is named
"""
Base.in(tag::AbstractString, ep::Endpoint) = in(tag, ep.children)

"""
getindex
returns the child with the specified tag name.
"""
function getindex(ep::Endpoint, tag::AbstractString)
    for cep in ep.children
        cep == tag && return cep
    end

    error("No endpoint named $tag.")
end

function getindex(ep::Endpoint,
                  tags::AbstractVector)
    isempty(tags) && return ep
    leaf = ep[shift!(tags)]
    while !isempty(tags)
        leaf = leaf[shift!(tags)]
    end
    leaf
end

function getindex(ep::Endpoint,
                  tag::AbstractString,
                  captured::Dict{Any, Any})
    ret = ep[tag]
    if typeof(ret) == DynamicEndpoint
        captured[ret.tag[2:end]] = ret.captured
    end
    (ret, captured)
end

"""
getindex returns the child endpoint of the parent endpoint that
matches the string
"""
function getindex(ep::Endpoint,
                  tags::AbstractVector,
                  captured::Dict{Any, Any})
    isempty(tags) && return (ep, captured)
    leaf, captured = ep[shift!(tags), captured]
    while !isempty(tags)
        leaf, captured = leaf[shift!(tags), captured]
    end
    leaf, captured
end

function push!(ep::Endpoint, o::Endpoint)
    if in(o.tag, ep)
        return ep[o.tag]
    end
    push!(ep.children, o)
    sort!(ep.children)
    o
end

"""
# push!
converts the tag into an endpoint, then pushes it the the chileren
of ep.
"""
function push!(ep::Endpoint, tag::AbstractString;
               dynamic_prefix=':')
    push!(ep, StaticEndpoint(tag))
    ep[tag]
end

"""
# push!

Takes a taglist::Vector{AbstractString}, shifts them one-by-one and pushes
the token onto the current endpoint, and returns
"""
function push!(ep::Endpoint, taglist::Vector)
    while !isempty(taglist)
        tag = shift!(taglist)
        ep = push!(ep, tag)
    end
    ep
end

"""
# buildtree

builds a routing tree from a taglist
"""
function buildtree(taglist::Vector; dynamic_identifer=':', root=StaticEndpoint())
    push!(root, taglist)
    return root
end

"""
# finalize!

this function converts the static endpoint nodes into dynamic endpoint
nodes.

This must be called after all of the routing is defined
"""
function finalize!(root::Endpoint; dynamic_prefix=':')
    if startswith(root.tag, dynamic_prefix) && typeof(root) == StaticEndpoint
        children = root.children
        hm = root.handlermap
        root = DynamicEndpoint(root.tag)
        root.children = children
        root.handlermap = hm
    end
    root.children = map(finalize!, root.children)
    sort!(root.children)
    return root
end


text:
# src/arithmetic.jl
# Arithmetic of Rationals

import Base: sign,
             signbit,
             abs,
             abs2,
             copysign,
             flipsign,
             BitSigned,
             +,
             -,
             *,
             /,
             //,
             ^,
             inv,
             fma,
             power_by_squaring,
             div,
             fld,
             cld,
             checked_mul

"""```
sign(x::AbstractRational)
    === oftype(x, sign(numerator(x)))
```
Return the sign of the rational `x`.
"""
sign(x::AbstractRational) = oftype(x, sign(numerator(x)))

"""```
signbit(x::AbstractRational)
    === signbit(numerator(x))
```
Return the sign bit of the rational `x`.
"""
signbit(x::AbstractRational) = signbit(numerator(x))

"""```
abs(x::AbstractRatio)
    === signbit(x) ? -x : x
```
Return the absolute value of the ratio `x`.
"""
abs(x::AbstractRatio) = signbit(x) ? -x : x

"""```
abs2(x::AbstractRatio)
    === x * x
```
Return the square of the absolute value of the ratio `x`.
"""
abs2(x::AbstractRatio) = x * x

"""```
copysign(x::AbstractRatio, y)
    === signbit(y) ? -abs(x) : abs(x)
```
Return the absolute value of the ratio `x` with the sign of `y`.
"""
copysign(x::AbstractRatio, y) = signbit(y) ? -abs(x) : abs(x)

"""```
flipsign(x::AbstractRatio, y)
    === signbit(y) ? -x : x
```
Return the ratio `x` with the sign flipped according to the sign of `y`.
"""
flipsign(x::AbstractRatio, y) = signbit(y) ? -x : x

"""```
+x::AbstractRatio
```
Return the ratio `x`.
"""
+(x::AbstractRatio) = build_ratio(typeof(x), +numerator(x), denominator(x))

"""```
-x::AbstractRatio
```
Return the negative of ratio `x`.
"""
-(x::AbstractRatio) = build_ratio(typeof(x), -numerator(x), denominator(x))

"""```
-x::AbstractRational{T <: BitSigned}
```
Return the negative of rational `x` with underlying bit signed integer type.
"""
function -(x::AbstractRational{T}) where {T<:BitSigned}
    if numerator(x) == typemin(T)
        throw(OverflowError("rational numerator is typemin(T)"))
    end
    return build_ratio(typeof(x), -numerator(x), denominator(x))
end

"""```
-x::AbstractRational{T <: Unsigned}
```
Return the negative of rational `x` with underlying unsigned integer type.
"""
function -(x::AbstractRational{T}) where {T<:Unsigned}
    if numerator(x) != zero(T)
        throw(OverflowError("cannot negate unsigned number"))
    end
    return x
end

"""```
x::AbstractRational * y::AbstractRational
```
Return the product of rational `x` with rational `y`.
"""
function *(x::AbstractRational, y::AbstractRational)
    xn, yd = divgcd(numerator(x), denominator(y))
    xd, yn = divgcd(denominator(x), numerator(y))
    return build_ratio(typeof(x), checked_mul(xn, yn), checked_mul(xd, yd))
end

"""```
x::AbstractRational / y
```
Divide the rational `x` by `y`.
"""
function /(x::AbstractRational, y)
    return build_ratio(typeof(x), numerator(x) * denominator(y), denominator(x) * numerator(y))
end

"""```
x / y::AbstractRational
```
Divide `x` by the rational `y`.
"""
function /(x, y::AbstractRational)
    return build_ratio(typeof(y), numerator(x) * denominator(y), denominator(x) * numerator(y))
end

"""```
x::AbstractRational // y
```
Rationally divide the rational `x` by `y`.
"""
function //(x::AbstractRational, y)
    xn, yn = divgcd(numerator(x), y)
    return xn // checked_mul(x.den, yn)
end

"""```
x // y::AbstractRational
```
Rationally divide `x` by the rational `y`.
"""
function //(x, y::AbstractRational)
    xn, yn = divgcd(x, numerator(y))
    return checked_mul(xn, denominator(y)) // yn
end

"""```
x::AbstractRational // y::AbstractRational
```
Rationally divide the rational `x` by the rational `y`.
"""
function //(x::AbstractRational, y::AbstractRational)
    xn, yn = divgcd(numerator(x), numerator(y))
    xd, yd = divgcd(denominator(x), denominator(y))
    return checked_mul(xn, yd) // checked_mul(xd, yn)
end

"""```
x::AbstractRatio ^ n::Integer
```
Exponentiate the ratio `x` by the integer `n`.
"""
function ^(x::AbstractRatio, n::Integer)
    return n >= zero(n) ? power_by_squaring(x, n) : power_by_squaring(inv(x), -n)
end

"""```
x ^ y::AbstractRatio
```
Exponentiate `x` by the ratio `y`.
"""
^(x, y::AbstractRatio) = x^(numerator(y) / denominator(y))

"""```
x::AbstractFloat ^ y::AbstractRatio
```
Exponentiate the floating point value `x` by the ratio `y`.
"""
^(x::T, y::AbstractRatio) where {T<:AbstractFloat} = x^convert(T, y)

"""```
z::Complex ^ y::AbstractRatio
```
Exponentiate the complex number `z` by the ratio `y`.
"""
^(z::Complex{T}, y::AbstractRatio) where {T} = z^convert(typeof(one(T)^y), y)

"""```
z::Complex{<:AbstractRatio} ^ n::Bool
```
Exponentiate the complex ratio `z` by the boolean `n`.
"""
^(z::Complex{<:AbstractRatio}, n::Bool) = n ? z : one(z)

"""```
z::Complex{<:AbstractRatio} ^ n::Integer
```
Exponentiate the complex ratio `z` by the integer `n`.
"""
function ^(z::Complex{<:AbstractRatio}, n::Integer)
    return n >= zero(n) ? power_by_squaring(z, n) : power_by_squaring(inv(z), -n)
end

"""```
inv(x::AbstractRatio)
    === build_ratio(typeof(x), denominator(x), numerator(x))
```
Invert the ratio `x`.
"""
inv(x::AbstractRatio) = build_ratio(typeof(x), denominator(x), numerator(x))

"""```
fma(x::AbstractRatio, y::AbstractRatio, z::AbstractRatio)
    === (x * y) + z
```
Return the fused Mulitply-Add of the ratios `x`, `y`, and `z`.
"""
fma(x::AbstractRatio, y::AbstractRatio, z::AbstractRatio) = (x * y) + z

"""```
div(x::AbstractRational, y::Integer, r::RoundingMode)
```
Integer divide the rational `x` by the integer `y` and round accoring to rounding mode `r`.
"""
function div(x::AbstractRational, y::Integer, r::RoundingMode)
    xn, yn = divgcd(numerator(x), y)
    return div(xn, checked_mul(denominator(x), yn), r)
end

"""```
div(x::Integer, y::AbstractRational, r::RoundingMode)
```
Integer divide the integer `x` by the rational `y` and round accoring to rounding mode `r`.
"""
function div(x::Integer, y::AbstractRational, r::RoundingMode)
    xn, yn = divgcd(x, numerator(y))
    return div(checked_mul(xn, denominator(y)), yn, r)
end

"""```
div(x::AbstractRational, y::AbstractRational, r::RoundingMode)
```
Integer divide the rational `x` by the rational `y` and round accoring to rounding mode `r`.
"""
function div(x::AbstractRational, y::AbstractRational, r::RoundingMode)
    xn, yn = divgcd(numerator(x), numerator(y))
    xd, yd = divgcd(denominator(x), denominator(y))
    return div(checked_mul(xn, yd), checked_mul(xd, yn), r)
end

"""```
div(x::AbstractRatio, y)
    === div(x, y, RoundToZero)
```
Integer divide the rational `x` by `y`.
"""
div(x::AbstractRatio, y) = div(x, y, RoundToZero)

"""```
fld(x::AbstractRatio, y)
    === div(x, y, RoundDown)
```
Integer divide the rational `x` by `y` and round down.
"""
fld(x::AbstractRatio, y) = div(x, y, RoundDown)

"""```
cld(x::AbstractRatio, y)
    === div(x, y, RoundUp)
```
Integer divide the rational `x` by `y` and round up.
"""
cld(x::AbstractRatio, y) = div(x, y, RoundUp)

"""```
div(x, y::AbstractRatio)
    === div(x, y, RoundToZero)
```
Integer divide `x` by the rational `y`.
"""
div(x, y::AbstractRatio) = div(x, y, RoundToZero)

"""```
fld(x, y::AbstractRatio)
    === div(x, y, RoundDown)
```
Integer divide `x` by the rational `y` and round down.
"""
fld(x, y::AbstractRatio) = div(x, y, RoundDown)

"""```
cld(x, y::AbstractRatio)
    === div(x, y, RoundUp)
```
Integer divide `x` by the rational `y` and round up.
"""
cld(x, y::AbstractRatio) = div(x, y, RoundUp)

"""```
div(x::AbstractRatio, y::AbstractRatio)
    === div(x, y, RoundToZero)
```
Integer divide the rational `x` by the rational `y`.
"""
div(x::AbstractRatio, y::AbstractRatio) = div(x, y, RoundToZero)

"""```
fld(x::AbstractRatio, y::AbstractRatio)
    === div(x, y, RoundDown)
```
Integer divide the rational `x` by the rational `y` and round down.
"""
fld(x::AbstractRatio, y::AbstractRatio) = div(x, y, RoundDown)

"""```
cld(x::AbstractRatio, y::AbstractRatio)
    === div(x, y, RoundUp)
```
Integer divide the rational `x` by the rational `y` and round up.
"""
cld(x::AbstractRatio, y::AbstractRatio) = div(x, y, RoundUp)


text:
module CURAND

using ..CuArrays
using ..CuArrays: libcurand

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

using CUDAnative

using CEnum

# core library
include("libcurand_common.jl")
include("error.jl")
include("libcurand.jl")

# low-level wrappers
include("wrappers.jl")

# high-level integrations
include("random.jl")

const handles_lock = ReentrantLock()
const created_generators = Dict{Tuple{UInt,Int},RNG}()
const active_generators = Vector{Union{Nothing,RNG}}()

function generator()
    tid = Threads.threadid()
    if @inbounds active_generators[tid] === nothing
        ctx = context()
        key = (objectid(ctx), tid)
        lock(handles_lock) do
            active_generators[tid] = get!(created_generators, key) do
                rng = RNG()
                Random.seed!(rng)
                rng
            end
        end
    end
    @inbounds active_generators[tid]
end

function __init__()
    resize!(active_generators, Threads.nthreads())
    fill!(active_generators, nothing)

    CUDAnative.atcontextswitch() do tid, ctx
        # we don't eagerly initialize handles, but do so lazily when requested
        active_generators[tid] = nothing
    end
end

end

const seed! = CURAND.seed!
const rand = CURAND.rand
const randn = CURAND.randn
const rand_logn = CURAND.rand_logn
const rand_poisson = CURAND.rand_poisson

@deprecate curand CuArrays.rand
@deprecate curandn CuArrays.randn


text:
@testset "matrix functions" begin
    @testset "LinearAlgebra.exp!(A::Matrix) frule" begin
        n = 10
        # each normalization hits a specific branch
        @testset "A::Matrix{$T}, opnorm(A,1)=$nrm" for T in (Float64, ComplexF64),
            nrm in (0.01, 0.1, 0.5, 1.5, 3.0, 6.0, 12.0)

            A = randn(T, n, n)
            A *= nrm / opnorm(A, 1)
            tols = nrm == 0.1 ? (atol=1e-8, rtol=1e-8) : NamedTuple()
            test_frule(LinearAlgebra.exp!, A; tols...)
        end
        @testset "imbalanced A" begin
            A = Float64[0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0]
            test_frule(LinearAlgebra.exp!, A)
        end
        @testset "hermitian A, T=$T" for T in (Float64, ComplexF64)
            A = Matrix(Hermitian(randn(T, n, n)))
            test_frule(LinearAlgebra.exp!, A)
            test_frule(LinearAlgebra.exp!, A ⊢ Matrix(Hermitian(randn(T, n, n))))
        end
    end

    @testset "exp(A::Matrix) rrule" begin
        n = 10
        # each normalization hits a specific branch
        @testset "A::Matrix{$T}, opnorm(A,1)=$nrm" for T in (Float64, ComplexF64),
            nrm in (0.01, 0.1, 0.5, 1.5, 3.0, 6.0, 12.0)

            A = randn(T, n, n)
            A *= nrm / opnorm(A, 1)
            # rrule is not inferrable, but pullback should be
            tols = nrm == 0.1 ? (atol=1e-8, rtol=1e-8) : NamedTuple()
            test_rrule(exp, A; check_inferred=false, tols...)
            Y, back = rrule(exp, A)
            @maybe_inferred back(rand_tangent(Y))
        end
        @testset "imbalanced A" begin
            A = Float64[0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0]
            test_rrule(exp, A; check_inferred=false)
        end
        @testset "hermitian A, T=$T" for T in (Float64, ComplexF64)
            A = Matrix(Hermitian(randn(T, n, n)))
            test_rrule(exp, A; check_inferred=false)
            test_rrule(
                exp,
                A;
                check_inferred=false,
                output_tangent=Matrix(Hermitian(randn(T, n, n))),
            )
        end
    end
end


text:
using ProgressMeter, JLD

include("Q1-Ising2D.jl")

cd(dirname(@__FILE__))

datapath = "../../Data/Q1/"

βList = hcat(0.0:0.005:1.0)
Llist = [5, 10, 15]
runnum = 50
Data = zeros(4, length(βList), runnum, length(Llist))

PRG = Progress(length(βList) * runnum)
for n ∈ 1:runnum
    for i ∈ 1:length(βList)
        for l ∈ 1:length(Llist)
            L = Llist[l]
            β = βList[i]
            Parameters = Dict(
                :dim => L,
                :MCLSize => 1000,
                :NLSize => 1000,
                :SkipNum => 200,
                :ProgBar => false)
            Data[:, i, n, l] = [Ising2D.IsingModel(β; Parameters...)...]
        end
        next!(PRG)
    end
end
save(datapath * "PTData.jld", "Data", Data, "βList", βList, "Llist", Llist)


text:
function printDocstrDict(mod::Module, binding::Symbol)
    print(Docs.meta(mod)[Docs.Binding(mod, binding)].docs)
end

max_stars_repo_name: AnirudhPraveen/ExamplePackage

max_issues_repo_name: AnirudhPraveen/ExamplePackage

max_forks_repo_name: AnirudhPraveen/ExamplePackage


text:
module BenchmarkingEconomicEfficiency

    """
    BenchmarkingEconomicEfficiency
    A Julia package for economic efficiency measurement using Data Envelopment Analysis (DEA).
    [BenchmarkingEconomicEfficiency repository](https://github.com/javierbarbero/BenchmarkingEconomicEfficiency.jl).
    """

    using SparseArrays
    using LinearAlgebra

    using Printf: @sprintf
    using Statistics: std
    using StatsBase: CoefTable

    import StatsBase: nobs, mean

    using Reexport
    @reexport using DataEnvelopmentAnalysis

    export
    # Types
    CostDDFDEAModel, RevenueDDFDEAModel,
    CostAdditiveDEAModel, RevenueAdditiveDEAModel, ProfitAdditiveDEAModel,
    CostRussellDEAModel, RevenueRussellDEAModel, ProfitRussellDEAModel,
    ProfitERGDEAModel,
    ProfitModifiedDDFDEAModel,
    CostHolderDEAModel, RevenueHolderDEAModel, ProfitHolderDEAModel,
    CostReverseDDFDEAModel, RevenueReverseDDFDEAModel, ProfitReverseDDFDEAModel,
    CostGDADEAModel, RevenueGDADEAModel, ProfitGDADEAModel,

    # Economic models
    deacostddf, dearevenueddf,
    deacostadd, dearevenueadd, deaprofitadd,
    deacostrussell, dearevenuerussell, deaprofitrussell,
    deaprofiterg,
    deaprofitmddf,
    deacostholder, dearevenueholder,deaprofitholder,
    deacostrddf, dearevenuerddf, deaprofitrddf,
    deacostgda, dearevenuegda, deaprofitgda

    # Include code of functions
    include("deacostddf.jl")
    include("dearevenueddf.jl")
    include("deacostadd.jl")
    include("dearevenueadd.jl")
    include("deaprofitadd.jl")
    include("deacostrussell.jl")
    include("dearevenuerussell.jl")
    include("deaprofitrussell.jl")
    include("deaprofiterg.jl")
    include("deaprofitmddf.jl")
    include("deacostholder.jl")
    include("dearevenueholder.jl")
    include("deaprofitholder.jl")
    include("deacostrddf.jl")
    include("dearevenuerddf.jl")
    include("deaprofitrddf.jl")
    include("deacostgda.jl")
    include("dearevenuegda.jl")
    include("deaprofitgda.jl")

end


text:
export ToricHomotopy, ToricHomotopyCache

"""
    ToricHomotopy(support, lifting, coefficients)

Construct the homotopy ``H(x, s) = Σ_α c_α e^{(α̂⋅γ̂) s} x^α``.
"""
mutable struct ToricHomotopy{S<:SPSystem, T} <: AbstractHomotopy
    system::S
    nterms::Vector{Int}
    support::Matrix{Int32}
    lifting::Vector{Int32}
    coeffs::Vector{T}
    s_weights::Vector{Float64}
    perms::Vector{Vector{Int}}
end

function ToricHomotopy(supports::Vector{<:Matrix}, coefficients::Vector{Vector{T}}) where T
    dummy_lifting = Vector{Int32}[collect(Int32(1):Int32(size(A, 2))) for A in supports]
    ToricHomotopy(supports, dummy_lifting, coefficients)
end
function ToricHomotopy(supports::Vector{<:Matrix}, liftings::Vector{Vector{Int32}}, coefficients::Vector{Vector{T}}) where T
    @assert all(length.(coefficients) .== size.(supports,2))
    @assert length(liftings) == length(supports)

    system = SPSystem(supports, coefficients)
    # StaticPolynomials possibly changes the order of the terms.
    # Therefore the coefficient vectors are maybe no more correct.
    # We correct this by applying the permutation applied to the columns
    # of the support

    nterms = length.(liftings)
    M = sum(nterms)
    n = length(liftings)
    support = Matrix{Int32}(undef, n, M)
    lifting = Vector{Int32}(undef, M)
    coeffs = Vector{T}(undef, M)
    perms = Vector{Vector{Int}}()
    k = 1
    for (i, f) in enumerate(system.system.polys)
        p = SP.permutation(f)
        for j in p
            @. support[:,k] = supports[i][:,j]
            lifting[k] = liftings[i][j]
            coeffs[k] = coefficients[i][j]
            k += 1
        end
        push!(perms, p)
    end
    s_weights = zeros(sum(nterms))

    ToricHomotopy(system, nterms, support, lifting, coeffs, s_weights, perms)
end

function update_cell!(H::ToricHomotopy, cell::MixedSubdivisions.MixedCell)
    n = length(H.nterms)
    k = 1
    s_max = 1.0
    s_min = Inf
    @inbounds for (i, m) in enumerate(H.nterms)
        βᵢ = cell.β[i]
        for _ in 1:m
            s_k = Float64(H.lifting[k]) - βᵢ
            for l in 1:n
                s_k += H.support[l,k] * cell.normal[l]
            end
            H.s_weights[k] = s_k

            s_max = max(s_max, s_k)
            # avoid zeroish values
            if s_k > 1e-8
                s_min = min(s_min, s_k)
            end

            k += 1
        end
    end
    # We normalize such that 5 is the highest power
    s_max /= 5
    for k in eachindex(H.s_weights)
        H.s_weights[k] /= s_max
    end

    s_min /= s_max

    s_min
end

function update_lifting!(H::ToricHomotopy, lifting::Vector{<:Vector{<:Integer}})
    k = 1
    for (i, p) in enumerate(H.perms), j in p
        H.lifting[k] = lifting[i][j]
        k += 1
    end
    H
end

"""
    ToricHomotopyCache

An simple cache for `ToricHomotopyCache`.
"""
mutable struct ToricHomotopyCache{C<:AbstractSystemCache, T} <: AbstractHomotopyCache
    system::C
    coeffs::Vector{T}
    coeffs_dt::Vector{T}
    views_coeffs::Vector{SubArray{T,1,Vector{T},Tuple{UnitRange{Int}}, true}}
    views_coeffs_dt::Vector{SubArray{T,1,Vector{T},Tuple{UnitRange{Int}}, true}}
    s::Float64
    ds::Float64
    active_coeffs::ActiveCoeffs
end

function cache(H::ToricHomotopy{S,T}, x, s) where {S,T}
    system = cache(H.system, x, s)
    U = promote_type(T, typeof(s))
    coeffs = convert.(U, H.coeffs)
    coeffs_dt = copy(H.coeffs)
    views_coeffs = [view(coeffs, 1:H.nterms[1])]
    views_coeffs_dt = [view(coeffs_dt, 1:H.nterms[1])]
    k = H.nterms[1] + 1
    for i in 2:length(H.nterms)
        push!(views_coeffs, view(coeffs, k:k+H.nterms[i]-1))
        push!(views_coeffs_dt, view(coeffs_dt, k:k+H.nterms[i]-1))
        k += H.nterms[i]
    end
    s = ds = NaN
    active_coeffs = COEFFS_UNKNOWN
    ToricHomotopyCache(system, coeffs, coeffs_dt, views_coeffs, views_coeffs_dt, s, ds, active_coeffs)
end

Base.size(H::ToricHomotopy) = size(H.system)

function update_coeffs!(cache::ToricHomotopyCache, H::ToricHomotopy, s)
    if s == cache.s
        if cache.active_coeffs != COEFFS_EVAL
            set_coefficients!(H.system, cache.views_coeffs)
            cache.active_coeffs == COEFFS_EVAL
        end
        return nothing
    end

    cache.s = s
    @inbounds for k in eachindex(cache.coeffs)
        cache.coeffs[k] = exp(H.s_weights[k] * s) * H.coeffs[k]
    end
    set_coefficients!(H.system, cache.views_coeffs)
    cache.active_coeffs == COEFFS_EVAL
    nothing
end

function update_coeffs_dt!(cache::ToricHomotopyCache, H::ToricHomotopy, s)
    if s == cache.ds
        if cache.active_coeffs != COEFFS_DT
            set_coefficients!(H.system, cache.views_coeffs_dt)
            cache.active_coeffs == COEFFS_DT
        end
        return nothing
    end

    if s == cache.s
        cache.ds = s
        @inbounds for k in eachindex(cache.coeffs)
            cache.coeffs_dt[k] = H.s_weights[k] * cache.coeffs[k]
        end
    else
        cache.ds = s
        @inbounds for k in eachindex(cache.coeffs)
            cache.coeffs_dt[k] = H.s_weights[k] * exp(H.s_weights[k] * s) * H.coeffs[k]
        end
    end
    set_coefficients!(H.system, cache.views_coeffs_dt)
    cache.active_coeffs == COEFFS_DT
    nothing
end

function evaluate!(u, H::ToricHomotopy, x, s, c::ToricHomotopyCache)
    update_coeffs!(c, H, real(s))
    @inbounds evaluate!(u, H.system, x, c.system)
    u
end

function dt!(u, H::ToricHomotopy, x, s, c::ToricHomotopyCache)
    update_coeffs_dt!(c, H, real(s))
    @inbounds evaluate!(u, H.system, x, c.system)
    u
end

function jacobian!(U, H::ToricHomotopy, x, s, c::ToricHomotopyCache)
    update_coeffs!(c, H, real(s))
    @inbounds jacobian!(U, H.system, x, c.system)
end

function evaluate_and_jacobian!(u, U, H::ToricHomotopy, x, s, c::ToricHomotopyCache)
    update_coeffs!(c, H, real(s))
    @inbounds evaluate_and_jacobian!(u, U, H.system, x, c.system)
    nothing
end

function evaluate(H::ToricHomotopy, x, s, c::ToricHomotopyCache)
    update_coeffs!(c, H, real(s))
    evaluate(H.system, x, c.system)
end

function dt(H::ToricHomotopy, x, s, c::ToricHomotopyCache)
    update_coeffs_dt!(c, H, real(s))
    evaluate(H.system, x, c.system)
end

function jacobian(H::ToricHomotopy, x, s, c::ToricHomotopyCache)
    update_coeffs!(c, H, real(s))
    jacobian(H.system, x, c.system)
end

(H::ToricHomotopy)(x, t, c=cache(H, x, t)) = evaluate(H, x, t, c)


text:
# Generate some synthetic data
#
using DelimitedFiles,Random
using NLS_Fit

Random.seed!(1234)

n = 50
model = Gaussian_Peak()
θ = Float64[2,n/2,3]

X=Float64[1:n;]
Y=eval_y(model,X,θ) + 0.1*(rand(n) .- 0.5)

writedlm("simple_gaussian.txt",hcat(X,Y))


text:
"""
LDA Slater exchange (DOI: 10.1017/S0305004100016108 and 10.1007/BF01340281)
"""
function lda_x!(ρ; α=1, E=nothing, Vρ=nothing)
    tmp = @. -α * cbrt(3/π * ρ)

    if E !== nothing
        E .= 3/4 .* tmp
    end
    if Vρ !== nothing
        Vρ .= tmp
    end

    (E=E, Vρ=Vρ)
end


text:
function array_rep_all(db)
    zeros(Int, length(db))
end

function indices_set_rarity(db, set, rarity)
    findall(x -> x["set"] == set && x["rarity"] == rarity, db)
end

function indices_rarity(db, rarity::Int)::Vector{Int}
    findall(x -> x["rarity"] == rarity, db)
end

function indices_set(db, set)
    findall(x -> x["set"] == set, db)
end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
f(i)=i=="dotnet ef"&&print(" "^13,"/\\__
"," "^7,raw"---==/    \\
 ___  ___   |.    \|\
| __|| __|  |  )   \\\
| _| | _|   \_/ |  //|\
|___||_|       /   \\/\\
")

max_stars_repo_name: pgrepds/Codegolf

max_issues_repo_name: pgrepds/Codegolf

max_forks_repo_name: pgrepds/Codegolf


text:
struct NonBitsUnion
    a::Union{String,Real}
end

@testset "Structs with non-bits unions" begin
    @test begin
        b = JlrsReflect.reflect([NonBitsUnion])
        sb = JlrsReflect.StringWrappers(b)

        sb[JlrsReflect.basetype(NonBitsUnion)] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.NonBitsUnion")]
        pub struct NonBitsUnion<'frame, 'data> {
            pub a: ::jlrs::wrappers::ptr::ValueRef<'frame, 'data>,
        }"""
    end
end


text:
using Distributions
using DataFrames

include("utils/mctools.jl")

include("models/mcs_def.jl")


text:
using StochasticDiffEq, DiffEqDevTools, LinearAlgebra, Random
Random.seed!(100)
dts = (1/2) .^ (10:-1:2) #14->7 good plot

const σ_const = 0.87
const μ = 1.01

u0 = rand(2)
A = [-2.0 1.0;1.0 -2.0]
B = Diagonal([σ_const for i in 1:2])

function f_commute(du,u,p,t)
  mul!(du,A,u)
  du .+= 1.01u
end

function f_commute_analytic(u0,p,t,W)
 tmp = (A+1.01I-2*(B^2))*t + B*sum(W)
 exp(tmp)*u0
end

function σ(du,u,p,t)
  du[1,1] = σ_const*u[1]
  du[2,1] = σ_const*u[2]
  du[1,2] = σ_const*u[1]
  du[2,2] = σ_const*u[2]
  du[1,3] = σ_const*u[1]
  du[2,3] = σ_const*u[2]
  du[1,4] = σ_const*u[1]
  du[2,4] = σ_const*u[2]
end

ff_commute = SDEFunction(f_commute,σ,analytic=f_commute_analytic)

prob = SDEProblem(ff_commute,σ,u0,(0.0,1.0),noise_rate_prototype=rand(2,4))

sol = solve(prob,RKMilCommute(),dt=1/2^(8))
sol = solve(prob,RKMil_General(ii_approx=IICommutative()),dt=1/2^(8))
sol = solve(prob,EM(),dt=1/2^(10))

dts = (1/2) .^ (10:-1:3) #14->7 good plot
sim2 = test_convergence(dts,prob,EM(),trajectories=Int(1e2))
sim2 = test_convergence(dts,prob,RKMilCommute(),trajectories=Int(2e2))
sim2 = test_convergence(dts,prob,RKMil_General(ii_approx=IICommutative()),trajectories=Int(2e2))

abs(sim2.𝒪est[:final] - 1) < 0.2


text:
# Run all our tests
include("compat.jl")
include("archive_utils.jl")
include("dependencies.jl")
include("platforms.jl")
include("prefix.jl")
include("products.jl")
include("sources.jl")
include("rootfs.jl")
include("runners.jl")