text:
"""
    rankUpdate!(C, A)
    rankUpdate!(C, A, α)
    rankUpdate!(C, A, α, β)

A rank-k update, C := β*C + α*A'A, of a Hermitian (Symmetric) matrix.

`α` and `β` both default to 1.0.  When `α` is -1.0 this is a downdate operation.
The name `rankUpdate!` is borrowed from [https://github.com/andreasnoack/LinearAlgebra.jl]
The order of the arguments
"""
function rankUpdate! end

function rankUpdate!(C::HermOrSym{T,S}, a::StridedVector{T},
        α=true) where {T<:BlasReal,S<:StridedMatrix}
    BLAS.syr!(C.uplo, T(α), a, C.data)
    C  ## to ensure that the return value is HermOrSym
end

function rankUpdate!(C::HermOrSym{T,S}, A::StridedMatrix{T},
        α=true, β=true) where {T<:BlasReal,S<:StridedMatrix}
    BLAS.syrk!(C.uplo, 'N', T(α), A, T(β), C.data)
    C
end

function rankUpdate!(C::HermOrSym{T,S}, A::SparseMatrixCSC{T}, α=true, β=true) where {T,S}
    m, n = size(A)
    @argcheck(m == size(C, 2), DimensionMismatch)
    @argcheck(C.uplo == 'L', ArgumentError)
    Cd = C.data
    isone(β) || rmul!(LowerTriangular(Cd), β)
    rv = rowvals(A)
    nz = nonzeros(A)
    @inbounds for jj in 1:n
        rangejj = nzrange(A, jj)
        lenrngjj = length(rangejj)
        for (k, j) in enumerate(rangejj)
            anzj = α * nz[j]
            rvj = rv[j]
            for i in k:lenrngjj
                kk = rangejj[i]
                Cd[rv[kk], rvj] += nz[kk] * anzj
            end
        end
    end
    C
end

function rankUpdate!(C::Diagonal{T}, A::SparseMatrixCSC{T}, α=true, β=true) where {T <: Number}
    m, n = size(A)
    dd = C.diag
    @argcheck(length(dd) == m, DimensionMismatch)
    isone(β) || rmul!(dd, β)
    nz = nonzeros(A)
    rv = rowvals(A)
    @inbounds for j in 1:n
        nzr = nzrange(A, j)
        if !isempty(nzr)
            length(nzr) == 1 || throw(ArgumentError("A*A' has off-diagonal elements"))
            k = nzr[1]
            dd[rv[k]] += α * abs2(nz[k])
        end
    end
    C
end

function rankUpdate!(C::HermOrSym{T,UniformBlockDiagonal{T}}, A::BlockedSparse{T},
        α=true) where {T}
    Arb = A.rowblocks
    Cdf = C.data.facevec
    (m = length(Arb)) == length(Cdf) ||
        throw(DimensionMismatch("length(A.rowblocks) = $m ≠ $(length(Cdf)) = length(C.data.facevec)"))
    for (b, d) in zip(Arb, Cdf)
        for v in b
            BLAS.syr!('L', α, v, d)
        end
    end
    C
end

rankUpdate!(C::HermOrSym{T,Matrix{T}}, A::BlockedSparse{T}, α=true) where {T} = rankUpdate!(C, A.cscmat, α)


text:
# input.jl

# This script defines the problem being solved in the form of a data structure.

# Modify the data structure manually to change the problem
# In the future we'd want scripts that take data in other forms and output to new input.jl files

###### Example Input ######

println("Input Importing Begun")

function getMolecule()
	println("Processing a Random Molecule")
	return Molecule(rand())
end

println("Input Importing Finished")

ext: jl

max_stars_repo_name: voutasaurus/metalbolicalgorithm

max_issues_repo_name: voutasaurus/metalbolicAlgorithm

max_forks_repo_name: voutasaurus/metalbolicAlgorithm


text:
#!/usr/bin/env julia

using RDatasets, Gadfly

import StatsBase
haskde = :kde in names(StatsBase)

tests = [
    ("points",                                6inch, 3inch),
    ("function_plots",                        6inch, 3inch),
    ("function_layers",                       6inch, 3inch),
    ("colored_points",                        6inch, 3inch),
    ("vstack",                                6inch, 6inch),
    ("hstack",                                6inch, 3inch),
    ("colorful_hist",                         6inch, 3inch),
    ("discrete_histogram",                    6inch, 3inch),
    ("discrete_histogram_horizontal",         6inch, 3inch),
    ("stacked_discrete_histogram",            6inch, 3inch),
    ("stacked_discrete_histogram_horizontal", 6inch, 3inch),
    ("stacked_continuous_histogram",          6inch, 3inch),
    ("dodged_discrete_histogram",             6inch, 3inch),
    ("array_aesthetics",                      6inch, 3inch),
    ("subplot_grid",                          6inch, 3inch),
    ("subplot_grid_histogram",                6inch, 3inch),
    ("labels",                                6inch, 6inch),
    ("percent",                               6inch, 6inch),
    ("timeseries_day",                        6inch, 3inch),
    ("timeseries_month",                      6inch, 3inch),
    ("timeseries_year_1",                     6inch, 3inch),
    ("timeseries_year_2",                     6inch, 3inch),
    ("timeseries_year_3",                     6inch, 3inch),
    ("timeseries_colorful",                   6inch, 3inch),
    ("custom_themes",                         6inch, 3inch),
    ("issue98",                               6inch, 3inch),
    ("issue82",                               6inch, 3inch),
    ("histogram2d",                           6inch, 3inch),
    ("rectbin",                               6inch, 3inch),
    ("density",                               6inch, 3inch),
    ("colorful_density",                      6inch, 3inch),
    ("explicit_colorkey_title",               6inch, 3inch),
    ("explicit_subplot_titles",               6inch, 3inch),
    ("subplot_grid_smooth",                   6inch, 8inch),
    ("histogram_errorbar",                    6inch, 3inch),
    ("issue106",                              6inch, 3inch),
    ("continuous_color_scale_range",          6inch, 3inch),
    ("continuous_scale_range",                6inch, 3inch),
    ("log10_scale_range",                     6inch, 3inch),
    ("histogram_explicit_bins",               6inch, 3inch),
    ("histogram2d_explicit_bins",             6inch, 3inch),
    ("explicit_number_format",                6inch, 3inch),
    ("explicit_xy_ticks",                     6inch, 3inch),
    ("boxplot",                               6inch, 3inch),
    ("subplot_categorical_bar",               6inch, 3inch),
    ("errorbar",                              6inch, 3inch),
    ("issue120",                              6inch, 3inch),
    ("histogram2d_discrete",                  6inch, 3inch),
    ("layer_themes",                          6inch, 3inch),
    ("discrete_color_manual",                 6inch, 3inch),
    ("ordered_line",                          6inch, 3inch),
    ("nan_skipping",                          6inch, 3inch),
    ("hexbin",                                6inch, 3inch),
    ("spy",                                   6inch, 3inch),
    ("issue177",                              6inch, 3inch),
    ("ribbon",                                6inch, 3inch)
]

backends = {
    "svg" => (name, width, height) -> SVG("$(name).svg", width, height),
    "d3"  => (name, width, height) -> D3("$(name).js",   width, height),
    "png" => (name, width, height) -> PNG("$(name).png", width, height),
    "ps"  => (name, width, height) -> PS("$(name).ps",   width, height),
    "pdf" => (name, width, height) -> PDF("$(name).pdf", width, height)
}

function run_tests(output_filename)
    testdir = Pkg.dir("Gadfly", "tests")
    whitelist = Set()
    if !isempty(ARGS)
        union!(whitelist, ARGS)
    else
        union!(whitelist, [name for (name, width, height) in tests])
    end

    for (name, width, height) in tests
        # 0.2 is currently in a weird state where kde is mid-migration and has
        # just ceased to exist. Until I figure that out, I'm just going to
        # disable density plots on 0.2. Sorry.
        if !haskde && (name == "density" || name == "colorful_density")
            continue
        end

        if !in(name, whitelist)
            continue
        end

        for (backend_name, backend) in backends
            println(STDERR, "Rendering $(name) on $(backend_name) backend.")
            try
                p = evalfile(joinpath(testdir, "$(name).jl"))
                draw(backend(name, width, height), p)
            catch
                println(STDERR, "FAILED!")
                rethrow()
            end
        end
    end

    d3src = Pkg.dir("Compose", "data", "d3.v3.min.js")
    gadflysrc = Pkg.dir("Gadfly", "src", "gadfly.js")

    output = open(output_filename, "w")
    print(output,
        """
        <!DOCTYPE html>

        <meta charset="utf-8" />

            <title>Gadfly Test Plots</title>

        """)

    for (name, width, height) in tests
        if !in(name, whitelist)
            continue
        end

        println(output, "
", name, "
")
        print(output, """

""")
        print(output, """""")
        print(output, """""")
        print(output, """""")
        print(output, """\n""")
    end

    print(output,
        """

        """)

    close(output)
end

run_tests("test.html")


text:
__precompile__()
module MathOptInterfaceGurobi

export GurobiOptimizer

using Gurobi
const GRB = Gurobi
using MathOptInterface
const MOI = MathOptInterface
using LinQuadOptInterface
const LQOI = LinQuadOptInterface

const SUPPORTED_OBJECTIVES = [
    LQOI.Linear,
    LQOI.Quad
]
const SUPPORTED_CONSTRAINTS = [
    (LQOI.Linear, LQOI.EQ),
    (LQOI.Linear, LQOI.LE),
    (LQOI.Linear, LQOI.GE),
    # (Linear, IV),
    (LQOI.Quad, LQOI.EQ),
    (LQOI.Quad, LQOI.LE),
    (LQOI.Quad, LQOI.GE),
    (LQOI.SinVar, LQOI.EQ),
    (LQOI.SinVar, LQOI.LE),
    (LQOI.SinVar, LQOI.GE),
    (LQOI.SinVar, LQOI.IV),
    (LQOI.SinVar, MOI.ZeroOne),
    (LQOI.SinVar, MOI.Integer),
    (LQOI.VecVar, LQOI.SOS1),
    (LQOI.VecVar, LQOI.SOS2),
    (LQOI.VecVar, MOI.Nonnegatives),
    (LQOI.VecVar, MOI.Nonpositives),
    (LQOI.VecVar, MOI.Zeros),
    (LQOI.VecLin, MOI.Nonnegatives),
    (LQOI.VecLin, MOI.Nonpositives),
    (LQOI.VecLin, MOI.Zeros)
]

mutable struct GurobiOptimizer <: LQOI.LinQuadOptimizer
    LQOI.@LinQuadOptimizerBase
    env
    params::Dict{String,Any}
    GurobiOptimizer(::Void) = new()
end

LQOI.LinQuadModel(::Type{GurobiOptimizer},env) = GRB.Model(env::GRB.Env,"defaultname")

function GurobiOptimizer(;kwargs...)

    env = GRB.Env()
    m = GurobiOptimizer(nothing)
    m.env = env
    m.params = Dict{String,Any}()
    MOI.empty!(m)
    for (name,value) in kwargs
        m.params[string(name)] = value
        GRB.setparam!(m.inner, string(name), value)
    end
    return m
end

function MOI.empty!(m::GurobiOptimizer)
    MOI.empty!(m,m.env)
    for (name,value) in m.params
        GRB.setparam!(m.inner, name, value)
    end
end

LQOI.lqs_supported_constraints(s::GurobiOptimizer) = SUPPORTED_CONSTRAINTS
LQOI.lqs_supported_objectives(s::GurobiOptimizer) = SUPPORTED_OBJECTIVES
#=
    inner wrapper
=#

#=
    Main
=#

# LinQuadSolver # Abstract type
# done above

# LQOI.lqs_setparam!(env, name, val)
# TODO fix this one
LQOI.lqs_setparam!(m::GurobiOptimizer, name, val) = GRB.setparam!(m.inner, string(name), val)

# LQOI.lqs_setlogfile!(env, path)
# TODO fix this one
LQOI.lqs_setlogfile!(m::GurobiOptimizer, path) = GRB.setlogfile(m.env, path::String)

# LQOI.lqs_getprobtype(m)
# TODO - consider removing, apparently useless

#=
    Constraints
=#

cintvec(v::Vector) = convert(Vector{Int32}, v)

_getsense(m::GurobiOptimizer, ::MOI.EqualTo{Float64}) = Cchar('=')
_getsense(m::GurobiOptimizer, ::MOI.LessThan{Float64}) = Cchar('<')
_getsense(m::GurobiOptimizer, ::MOI.GreaterThan{Float64}) = Cchar('>')
_getsense(m::GurobiOptimizer, ::MOI.Zeros)        = Cchar('=')
_getsense(m::GurobiOptimizer, ::MOI.Nonpositives) = Cchar('<')
_getsense(m::GurobiOptimizer, ::MOI.Nonnegatives) = Cchar('>')
_getboundsense(m::GurobiOptimizer, ::MOI.Nonpositives) = Cchar('>')
_getboundsense(m::GurobiOptimizer, ::MOI.Nonnegatives) = Cchar('<')

# LQOI.lqs_chgbds!(m, colvec, valvec, sensevec)
# TODO - improve single type
function LQOI.lqs_chgbds!(instance::GurobiOptimizer, colvec, valvec, sensevec)
    lb_len = count(x->x==Cchar('L'), sensevec)
    LB_val = Array{Float64}(0)
    sizehint!(LB_val, lb_len)
    LB_col = Array{Cint}(0)
    sizehint!(LB_col, lb_len)

    ub_len = count(x->x==Cchar('U'), sensevec)
    UB_val = Array{Float64}(0)
    sizehint!(UB_val, ub_len)
    UB_col = Array{Cint}(0)
    sizehint!(UB_col, ub_len)

    for i in eachindex(valvec)
        if sensevec[i] == Cchar('L')
            push!(LB_col, colvec[i])
            push!(LB_val, valvec[i])
        elseif sensevec[i] == Cchar('U')
            push!(UB_col, colvec[i])
            push!(UB_val, valvec[i])
        end
    end

    if lb_len > 0
        GRB.set_dblattrlist!(instance.inner, "LB", LB_col, LB_val)
    end

    if ub_len > 0
        GRB.set_dblattrlist!(instance.inner, "UB", UB_col, UB_val)
    end

    GRB.update_model!(instance.inner)
    nothing
end

# LQOI.lqs_getlb(m, col)
LQOI.lqs_getlb(instance::GurobiOptimizer, col) = (GRB.update_model!(instance.inner);GRB.get_dblattrlist( instance.inner, "LB", GRB.ivec(col))[1])
# LQOI.lqs_getub(m, col)
LQOI.lqs_getub(instance::GurobiOptimizer, col) = GRB.get_dblattrlist( instance.inner, "UB", GRB.ivec(col))[1]

# LQOI.lqs_getnumrows(m)
LQOI.lqs_getnumrows(instance::GurobiOptimizer) = GRB.num_constrs(instance.inner)

# LQOI.lqs_addrows!(m, rowvec, colvec, coefvec, sensevec, rhsvec)
LQOI.lqs_addrows!(instance::GurobiOptimizer, rowvec, colvec, coefvec, sensevec, rhsvec) = (GRB.add_constrs!(instance.inner, rowvec, colvec, coefvec, sensevec, rhsvec);GRB.update_model!(instance.inner))

# LQOI.lqs_getrhs(m, rowvec)
LQOI.lqs_getrhs(instance::GurobiOptimizer, row) = GRB.get_dblattrlist( instance.inner, "RHS", GRB.ivec(row))[1]

# colvec, coef = LQOI.lqs_getrows(m, rowvec)
# TODO improve
function LQOI.lqs_getrows(instance::GurobiOptimizer, idx)
    A = GRB.get_constrs(instance.inner, idx, 1)'
    return A.rowval-1, A.nzval
end

# LQOI.lqs_getcoef(m, row, col) #??
# TODO improve
function LQOI.lqs_getcoef(instance::GurobiOptimizer, row, col) #??
    return getcoeff(model::Model, row::Integer, col::Integer)
    # A = GRB.get_rows(m, row, row)'
    # cols = A.rowval
    # vals = A.nzval

    # pos = findfirst(cols, col)
    # if pos > 0
    #     return vals[pos]
    # else
    #     return 0.0
    # end
end

# LQOI.lqs_chgcoef!(m, row, col, coef)
# TODO SPLIT THIS ONE
function LQOI.lqs_chgcoef!(instance::GurobiOptimizer, row, col, coef)
    if row == 0
        GRB.set_dblattrlist!(instance.inner, "Obj", Cint[col], Float64[coef])
    elseif col == 0
        GRB.set_dblattrlist!(instance.inner, "RHS", Cint[row], Float64[coef])
    else
        GRB.chg_coeffs!(instance.inner, row, col, coef)
        #TODO fix this function in gurobi
    end
end

# LQOI.lqs_delrows!(m, row, row)
LQOI.lqs_delrows!(instance::GurobiOptimizer, rowbeg, rowend) = GRB.del_constrs!(instance.inner, cintvec(collect(rowbeg:rowend)))

# LQOI.lqs_chgctype!(m, colvec, typevec)
# TODO fix types
LQOI.lqs_chgctype!(instance::GurobiOptimizer, colvec, typevec) = GRB.set_charattrlist!(instance.inner, "VType", GRB.ivec(colvec), GRB.cvec(typevec))

# LQOI.lqs_chgsense!(m, rowvec, sensevec)
# TODO fix types
LQOI.lqs_chgsense!(instance::GurobiOptimizer, rowvec, sensevec) = GRB.set_charattrlist!(instance.inner, "Sense", GRB.ivec(rowvec), GRB.cvec(sensevec))

const VAR_TYPE_MAP = Dict{Symbol,Cchar}(
    :CONTINUOUS => Cchar('C'),
    :INTEGER => Cchar('I'),
    :BINARY => Cchar('B')
)
LQOI.lqs_vartype_map(m::GurobiOptimizer) = VAR_TYPE_MAP

# LQOI.lqs_addsos(m, colvec, valvec, typ)
LQOI.lqs_addsos!(instance::GurobiOptimizer, colvec, valvec, typ) = (GRB.add_sos!(instance.inner, typ, colvec, valvec);GRB.update_model!(instance.inner))
# LQOI.lqs_delsos(m, idx, idx)
LQOI.lqs_delsos!(instance::GurobiOptimizer, idx1, idx2) = (GRB.del_sos!(instance.inner, cintvec(collect(idx1:idx2)));GRB.update_model!(instance.inner))

const SOS_TYPE_MAP = Dict{Symbol,Symbol}(
    :SOS1 => :SOS1,#Cchar('1'),
    :SOS2 => :SOS2#Cchar('2')
)
LQOI.lqs_sertype_map(m::GurobiOptimizer) = SOS_TYPE_MAP

# LQOI.lqs_getsos(m, idx)
# TODO improve getting processes
function LQOI.lqs_getsos(instance::GurobiOptimizer, idx)
    A, types = GRB.get_sos_matrix(instance.inner)
    line = A[idx,:] #sparse vec
    cols = line.nzind
    vals = line.nzval
    typ = types[idx] == Cint(1) ? :SOS1 : :SOS2
    return cols, vals, typ
end

# LQOI.lqs_getnumqconstrs(m)
LQOI.lqs_getnumqconstrs(instance::GurobiOptimizer) = GRB.num_qconstrs(instance.inner)

# LQOI.lqs_addqconstr(m, cols,coefs,rhs,sense, I,J,V)
LQOI.lqs_addqconstr!(instance::GurobiOptimizer, cols,coefs,rhs,sense, I,J,V) = GRB.add_qconstr!(instance.inner, cols, coefs, I, J, V, sense, rhs)

# LQOI.lqs_chgrngval
LQOI.lqs_chgrngval!(instance::GurobiOptimizer, rows, vals) = GRB.chg_rhsrange!(instance.inner, cintvec(rows), -vals)

const CTR_TYPE_MAP = Dict{Symbol,Cchar}(
    :RANGE => Cchar('R'),
    :LOWER => Cchar('L'),
    :UPPER => Cchar('U'),
    :EQUALITY => Cchar('E')
)
LQOI.lqs_ctrtype_map(m::GurobiOptimizer) = CTR_TYPE_MAP

#=
    Objective
=#

# LQOI.lqs_copyquad(m, intvec,intvec, floatvec) #?
function LQOI.lqs_copyquad!(instance::GurobiOptimizer, I, J, V)
    GRB.delq!(instance.inner)
    for i in eachindex(V)
        if I[i] == J[i]
            V[i] /= 2
        end
    end
    GRB.add_qpterms!(instance.inner, I, J, V)
    return nothing
end

# LQOI.lqs_chgobj(m, colvec,coefvec)
function LQOI.lqs_chgobj!(instance::GurobiOptimizer, colvec, coefvec)
    nvars = GRB.num_vars(instance.inner)
    obj = zeros(Float64, nvars)

    for i in eachindex(colvec)
        obj[colvec[i]] = coefvec[i]
    end

    GRB.set_dblattrarray!(instance.inner, "Obj", 1, num_vars(instance.inner), obj)
    GRB.update_model!(instance.inner)
    nothing
end

# LQOI.lqs_chgobjsen(m, symbol)
# TODO improve min max names
function LQOI.lqs_chgobjsen!(instance::GurobiOptimizer, symbol)
    if symbol in [:minimize, :Min]
        GRB.set_sense!(instance.inner, :minimize)
    else
        GRB.set_sense!(instance.inner, :maximize)
    end
    GRB.update_model!(instance.inner)
end

# LQOI.lqs_getobj(instance.inner)
LQOI.lqs_getobj(instance::GurobiOptimizer) = GRB.get_dblattrarray( instance.inner, "Obj", 1, num_vars(instance.inner)   )

# lqs_getobjsen(m)
function LQOI.lqs_getobjsen(instance::GurobiOptimizer)
    s = GRB.model_sense(instance.inner)
    if s in [:maximize, :Max]
        return MOI.MaxSense
    else
        return MOI.MinSense
    end
end

#=
    Variables
=#

# LQOI.lqs_getnumcols(m)
LQOI.lqs_getnumcols(instance::GurobiOptimizer) = (GRB.update_model!(instance.inner); GRB.num_vars(instance.inner))

# LQOI.lqs_newcols!(m, int)
LQOI.lqs_newcols!(instance::GurobiOptimizer, int) = (GRB.add_cvars!(instance.inner, zeros(int));GRB.update_model!(instance.inner))

# LQOI.lqs_delcols!(m, col, col)
LQOI.lqs_delcols!(instance::GurobiOptimizer, col, col2) = (GRB.del_vars!(instance.inner, col);GRB.update_model!(instance.inner))

# LQOI.lqs_addmipstarts(m, colvec, valvec)
function LQOI.lqs_addmipstarts!(instance::GurobiOptimizer, colvec, valvec)
    x = zeros(GRB.num_vars(instance.inner))
    for i in eachindex(colvec)
        x[colvec[i]] = valvec[i]
    end
    GRB.loadbasis(instance.inner, x)
end
#=
    Solve
=#

# LQOI.lqs_mipopt!(m)
LQOI.lqs_mipopt!(instance::GurobiOptimizer) = LQOI.lqs_lpopt!(instance)

# LQOI.lqs_qpopt!(m)
LQOI.lqs_qpopt!(instance::GurobiOptimizer) = LQOI.lqs_lpopt!(instance)

# LQOI.lqs_lpopt!(m)
LQOI.lqs_lpopt!(instance::GurobiOptimizer) = (GRB.update_model!(instance.inner);GRB.optimize(instance.inner))

# LQOI.lqs_terminationstatus(m)
function LQOI.lqs_terminationstatus(instance::GurobiOptimizer)

    stat = get_status(instance.inner)

    if stat == :loaded
        return MOI.OtherError
    elseif stat == :optimal
        return MOI.Success
    elseif stat == :infeasible
        if hasdualray(instance)
            return MOI.Success
        else
            return MOI.InfeasibleNoResult
        end
    elseif stat == :inf_or_unbd
        return MOI.InfeasibleOrUnbounded
    elseif stat == :unbounded
        if hasprimalray(instance)
            return MOI.Success
        else
            return MOI.UnboundedNoResult
        end
    elseif stat == :cutoff
        return MOI.ObjectiveLimit
    elseif stat == :iteration_limit
        return MOI.IterationLimit
    elseif stat == :node_limit
        return MOI.NodeLimit
    elseif stat == :time_limit
        return MOI.TimeLimit
    elseif stat == :solution_limit
        return MOI.SolutionLimit
    elseif stat == :interrupted
        return MOI.Interrupted
    elseif stat == :numeric
        return MOI.NumericalError
    elseif stat == :suboptimal
        return MOI.OtherLimit
    elseif stat == :inprogress
        return MOI.OtherError
    elseif stat == :user_obj_limit
        return MOI.ObjectiveLimit
    end
    return MOI.OtherError
end

function LQOI.lqs_primalstatus(instance::GurobiOptimizer)

    stat = get_status(instance.inner)

    if stat == :optimal
        return MOI.FeasiblePoint
    elseif stat == :solution_limit
        return MOI.FeasiblePoint
    elseif stat in [:inf_or_unbd, :unbounded] && hasprimalray(instance)
        return MOI.InfeasibilityCertificate
    elseif stat == :suboptimal
        return MOI.FeasiblePoint
    else
        return MOI.UnknownResultStatus
    end
end
function LQOI.lqs_dualstatus(instance::GurobiOptimizer)
    stat = get_status(instance.inner)

    if GRB.is_mip(instance.inner) || GRB.is_qcp(instance.inner)
        return MOI.UnknownResultStatus
    else
        if stat == :optimal
            return MOI.FeasiblePoint
        elseif stat == :solution_limit
            return MOI.FeasiblePoint
        elseif stat in [:inf_or_unbd, :infeasible] && hasdualray(instance)
            return MOI.InfeasibilityCertificate
        elseif stat == :suboptimal
            return MOI.FeasiblePoint
        else
            return MOI.UnknownResultStatus
        end
    end
end

# LQOI.lqs_getx!(m, place)
LQOI.lqs_getx!(instance::GurobiOptimizer, place) = GRB.get_dblattrarray!(place, instance.inner, "X", 1)

# LQOI.lqs_getax!(m, place)
function LQOI.lqs_getax!(instance::GurobiOptimizer, place)
    GRB.get_dblattrarray!(place, instance.inner, "Slack", 1)
    rhs = GRB.get_dblattrarray(instance.inner, "RHS", 1, num_constrs(instance.inner))
    for i in eachindex(place)
        place[i] = -place[i]+rhs[i]
    end
    nothing
end
# LQOI.lqs_getdj!(m, place)
LQOI.lqs_getdj!(instance::GurobiOptimizer, place) = GRB.get_dblattrarray!(place, instance.inner, "RC", 1)

# LQOI.lqs_getpi!(m, place)
LQOI.lqs_getpi!(instance::GurobiOptimizer, place) = GRB.get_dblattrarray!(place, instance.inner, "Pi", 1)

# LQOI.lqs_getobjval(m)
LQOI.lqs_getobjval(instance::GurobiOptimizer) = GRB.get_objval(instance.inner)

# LQOI.lqs_getbestobjval(m)
LQOI.lqs_getbestobjval(instance::GurobiOptimizer) = GRB.get_objval(instance.inner)

# LQOI.lqs_getmiprelgap(m)
function LQOI.lqs_getmiprelgap(instance::GurobiOptimizer)
    L = GRB.get_objval(instance.inner)
    U = GRB.get_objbound(instance.inner)
    return abs(U-L)/U
end

# LQOI.lqs_getitcnt(m)
LQOI.lqs_getitcnt(instance::GurobiOptimizer)  = GRB.get_iter_count(instance.inner)

# LQOI.lqs_getbaritcnt(m)
LQOI.lqs_getbaritcnt(instance::GurobiOptimizer) = GRB.get_barrier_iter_count(instance.inner)

# LQOI.lqs_getnodecnt(m)
LQOI.lqs_getnodecnt(instance::GurobiOptimizer) = GRB.get_node_count(instance.inner)

# LQOI.lqs_dualfarkas(m, place)
LQOI.lqs_dualfarkas!(instance::GurobiOptimizer, place) = GRB.get_dblattrarray!(place, instance.inner, "FarkasDual", 1)

function hasdualray(instance::GurobiOptimizer)
    try
        GRB.get_dblattrarray(instance.inner, "FarkasDual", 1, GRB.num_constrs(instance.inner))
        return true
    catch
        return false
    end
end

# LQOI.lqs_getray(m, place)
LQOI.lqs_getray!(instance::GurobiOptimizer, place) = GRB.get_dblattrarray!(place, instance.inner, "UnbdRay", 1)

function hasprimalray(instance::GurobiOptimizer)
    try
        GRB.get_dblattrarray(instance.inner, "UnbdRay", 1, GRB.num_vars(instance.inner))
        return true
    catch
        return false
    end
end

MOI.free!(m::GurobiOptimizer) = GRB.free_model(m.inner)

"""
    writeproblem(m: :MOI.AbstractOptimizer, filename::String)
Writes the current problem data to the given file.
Supported file types are solver-dependent.
"""
writeproblem(m::GurobiOptimizer, filename::String, flags::String="") = GRB.write_model(m.inner, filename)

# blocked
MOI.addconstraint!(m::GurobiOptimizer, func::LQOI.Linear, set::LQOI.IV) = error("not supported")
MOI.addconstraints!(m::GurobiOptimizer, func::Vector{LQOI.Linear}, set::Vector{LQOI.IV}) = error("not supported")

MOI.canget(m::GurobiOptimizer, any, c::LQOI.LCI{LQOI.IV}) = false
MOI.canmodifyconstraint(m::GurobiOptimizer, c::LQOI.LCI{LQOI.IV}, chg) = false
MOI.candelete(m::GurobiOptimizer, c::LQOI.LCI{LQOI.IV}) = false

end # module


text:
@testset "Primatives" begin

    @test True :: TrueOrFalse == True
    @test False :: TrueOrFalse == False

end


text:
using Feather
using DataFrames
using ProgressMeter
# using CSV

function loadSave(filename)
    return Feather.read(filename) #save (domain) or save2 (article)
end

function loadBias(filename)
    # CSV.read is necessary to update to new dataframes
    # bias = CSV.read(filename, header=false)
    bias = readtable(filename, header=false)
    code = r"^(([^:/?#]+):)?(//([^/?#]*))?([^?#]*)(\?([^#]*))?(#(.*))?"

    biasnames = String[]
    p = Progress(size(bias,1), 0.1)
    for i in 1:size(bias,1)
        if ~isna(bias[i, 5])
            url = match(code, AbstractString(bias[i, 5]))
            bias[i, 5] = replace(url[4], "www.", "")
            push!(biasnames, bias[i, 5])
        else
            push!(biasnames, "NA")
        end
    end
    return bias, biasnames
end

"""
    buildWeb(df, by_domain=true)

Build a data structure for storing the links in an indexed form.

The first key is the type of link one of: ["a", "link", "script", "img", "mut"]
where a,link,script, and img are html tag types, and mut is "mutual a links".

df is a dataframe with columns sdom, ddom, type as the first three columns.

if by_domain is false, then we look for the 4th and 5th columns which have the filtered
article identifiers.
"""
function buildWeb(df, by_domain=true)
    cats = Dict{String, Dict}()
    for key in ["a", "link", "script", "img", "mut"]
        cats[key] = Dict{String, Dict}()
    end

    p = Progress(size(df, 1), 0.1)
    for i in 1:size(df, 1)
       # discriminating between save and save2 feather files
        if by_domain
          # by domain
          sdom = get(df[i, 1]) # sdom
          ddom = get(df[i, 2]) # ddom
        else
          # by article
          sdom = get(df[i, 4]) # rsrc
          ddom = get(df[i, 5]) # rdest
        end
        if isascii(sdom) && isascii(ddom)
            cat = cats[get(df[i, 3])] # linktype
            if ~haskey(cat, sdom)
                cat[sdom] = Dict{String, Int}()
            end
            if ~haskey(cat[sdom], ddom)
                cat[sdom][ddom] = 1
            else
                cat[sdom][ddom] += 1
            end
        end
        next!(p)
    end

    p = Progress(length(cats["a"]), 0.1)
    for src in keys(cats["a"])
        for dst in keys(cats["a"][src])
            if haskey(cats["a"], dst) && haskey(cats["a"][dst], src)
                if ~haskey(cats["mut"], src)
                     cats["mut"][src] = Dict{String, Int}()
                end
                cats["mut"][src][dst] = cats["a"][src][dst]
            end
        end
        next!(p)
    end
    return cats
end

max_stars_repo_name: jpfairbanks/NewsGraphing

max_issues_repo_name: jpfairbanks/NewsGraphing

max_forks_repo_name: jpfairbanks/NewsGraphing


text:
function convergence_test(tol=1e-3)
	global pa
	FBD.fbd!(pa);
	@test Misfits.error_after_scaling(pa.pfibd.optm.cal.g,pa.pfibd.optm.obs.g)[1] < tol
	@test Misfits.error_after_scaling(pa.pfibd.optm.cal.s[1:pa.plsbd.om.nts-1],pa.pfibd.optm.obs.s[1:pa.plsbd.om.nts-1])[1] < tol
	@test Misfits.error_after_scaling(pa.plsbd.optm.cal.g,pa.plsbd.optm.obs.g)[1] < tol
	@test Misfits.error_after_scaling(pa.plsbd.optm.cal.s,pa.plsbd.optm.obs.s)[1] < tol

end

FBD.__init__(optg="optim", opts="optim")
pa=FBD.random_problem();
convergence_test(1e-3)

FBD.__init__()
pa=FBD.random_problem();
convergence_test(1e-3)


text:
immutable PoissonProcess
    lambda::Float64
end

type PoissonState
    t::Vector{Float64}
    n::Int64
    k::Int64
    y::Vector{Float64}
    dt::Vector{Float64}

    PoissonState(t, k) = new(t, length(t), k, zeros(k), diff(t))
end

function next(p::PoissonProcess, s::PoissonState, i)
    if i == 1
        return s.y
    end

    s.y .= s.y .+ rand(Poisson(p.lambda * s.dt[i-1]), s.k)
    s.y
end

function cumsim(p::PoissonProcess, t, k)
    s = PoissonState(t, k)
    result = zeros(s.n, s.k)
    for i in 1:s.n
        y = next(p, s, i)
        result[i,:] .= y[:]
    end
    result
end

function sim(p::PoissonProcess, t, k)
    s = PoissonState(t, k)
    for i in 1:(s.n-1)
        next(p, s, i)
    end
    next(p, s, s.n)
end

distribution(p::PoissonProcess, t) = Poisson(p.lambda * t)


text:
module Crash

using StochasticPrograms
using MathProgBase

abstract type CrashMethod end

struct None <: CrashMethod end

function (crash::None)(sp::JuMP.Model,solver::MathProgBase.AbstractMathProgSolver)
    return rand(sp.numCols)
end

struct EVP <: CrashMethod end

function (crash::EVP)(sp::JuMP.Model,solver::MathProgBase.AbstractMathProgSolver)
    evp = StochasticPrograms.EVP(sp,solver)
    status = solve(evp)
    status != :Optimal && error("Could not solve EVP model during crash procedure. Aborting.")
    return evp.colVal[1:sp.numCols]
end

struct Scenario{S <: AbstractScenarioData} <: CrashMethod
    scenario::S

    function (::Type{Scenario})(scenario::S) where S <: AbstractScenarioData
        return new{S}(scenario)
    end
end

function (crash::Scenario)(sp::JuMP.Model,solver::MathProgBase.AbstractMathProgSolver)
    ws = WS(sp,crash.scenario,solver)
    status = solve(ws)
    status != :Optimal && error("Could not solve wait-and-see model during crash procedure. Aborting.")
    return ws.colVal[1:sp.numCols]
end

struct Custom{T <: Real} <: CrashMethod
    x₀::Vector{T}

    function (::Type{Custom})(x₀::Vector{T}) where T <: Real
        return new{T}(x₀)
    end
end

function (crash::Custom)(sp::JuMP.Model,solver::MathProgBase.AbstractMathProgSolver)
    return x₀
end

end


text:
#
# This script is intended to see whether we have all the interface we need
# whether we need additional structure, and see how things would look like
# from a developer and user perspective.
#

using Distributed

using JuLIP, ACE, ACEfit, ACEatoms

using ACEfit: eval_obs, vec_obs

# load some example observation types
include("obsexamples.jl")
@everywhere include("obsexamples.jl")

OE = ObsExamples.ObsPotentialEnergy
OF = ObsExamples.ObsForces

## [1] create a dataset

function create_dataset(N)
   sw = StillingerWeber()
   data = Dat[]
   for _=1:N
      at = bulk(:Si, cubic=true) * 3
      rattle!(at, 0.1)
      dat = Dat(at, "bulk",
                [ eval_obs(OE, sw, at),
                  eval_obs(OF, sw, at) ] )
      push!(data, dat)
   end
   return data
end

data = create_dataset(10)

## [2] create a model, in this case a linear model represented by a basis
Bsel = SimpleSparseBasis(3, 8)
B1p = ACE.Utils.RnYlm_1pbasis(r0 = rnn(:Si), rcut = cutoff(StillingerWeber()))
basis = ACE.SymmetricBasis(ACE.Invariant(), B1p, Bsel)
model = ACEatoms.ACESitePotential(Dict(:Si => ACE.LinearACEModel(basis)))

# we now need to explain how to evaluate the observations on a basis; so we
# need to implement `energy` and `forces`. But this is already done for us
# in ACEatoms ...

## [3] assemble a loss function

loss(model, dat::Dat) = sum(
      sum(abs2, vec_obs(o) - vec_obs( eval_obs(typeof(o), model, dat.config) ))
      for o in ACEfit.observations(dat) )

function Loss(model, data)
   L = 0.0
   ACEfit.siterate( (i, dat) -> (L += loss(model, dat)), data )
   return L
end

println("Loss before fit:  ", Loss(model, data))

# note that the gradient of loss would always be obtained via AD, we never
# do this by hand.

## [4] assemble a lsq system
# this is a little more challenging and still needs a few fixes...
# at the moment, I'm "hacking" it so it will do something sensible but it
# is not a good solution yet.

import ACEbase

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# temporary hack to get get_basis and set_params! into ACEfit
ACEfit.get_basis(m::typeof(model)) = ACEatoms.basis(m)
ACEfit.eval(:(function set_params! end))
ACEfit.set_params!(m::typeof(model), args...) = ACE.set_params!(m, args...)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

ACEfit.llsq!(model, data, :serial; solver = ACEfit.LSQR())
#ACEfit.llsq!(model, data, :dist; solver = ACEfit.LSQR())

# and we can confirm that this actually makes the loss small :)
println("Loss after fit:  ", Loss(model, data))


text:
#unum-sqrs.jl

#square and square root

#computes the square of a unum x.  For ubounds we must do a check.
@universal function sqr(x::Ubound)
  signcode::UInt16 = 0
  is_neg_def(x.lower) && (signcode += 1)
  is_neg_def(x.upper) && (signcode += 2)

  #parse through the signcode possibilities
  if signcode == 0
    lower_result = resolve_lower(x.lower * x.lower)
    upper_result = resolve_upper(x.upper * x.upper)

    is_ulp(lower_result) && is_ulp(upper_result) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  elseif signcode == 1
    lower_result = resolve_upper(x.lower * x.lower)
    upper_result = resolve_upper(x.upper * x.upper)

    B(zero(U), max(lower_result, upper_result))
  #signcode 2 is impossible
  elseif signcode == 3
    lower_result = resolve_lower(x.upper * x.upper)
    upper_result = resolve_upper(x.lower * x.lower)

    is_ulp(lower_result) && is_ulp(upper_result) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  end
end

export sqr

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using ColorBitstring, Test

@testset "ColorBitstring" begin
    io = IOBuffer()
    printbits(io, Float16(1))
    @test String(take!(io)) == "\e[31m0\e[0m\e[32m01111\e[0m\e[34m0000000000\e[0m"
    printbits(io, Float32(1))
    @test String(take!(io)) == "\e[31m0\e[0m\e[32m01111111\e[0m\e[34m00000000000000000000000\e[0m"
    printbits(io, 1.0)
    @test String(take!(io)) == "\e[31m0\e[0m\e[32m01111111111\e[0m\e[34m0000000000000000000000000000000000000000000000000000\e[0m"
    printbits(io, 1)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

type Mortar <: BoundaryProblem
    dimension :: Int
    rotate_normals :: Bool
    adjust :: Bool
    dual_basis :: Bool
    use_forwarddiff :: Bool
    distval :: Float64
    store_fields :: Vector{Symbol}
end

function Mortar()
    default_fields = []
    return Mortar(-1, false, false, false, false, Inf, default_fields)
end

function get_unknown_field_name(problem::Problem{Mortar})
    return "reaction force"
end

function get_formulation_type(problem::Problem{Mortar})
    if problem.properties.use_forwarddiff
        return :forwarddiff
    else
        return :incremental
    end
end

function assemble!(problem::Problem{Mortar}, time::Float64)
    if problem.properties.dimension == -1
        problem.properties.dimension = dim = size(first(problem.elements), 1)
        info("assuming dimension of mesh tie surface is $dim")
        info("if this is wrong set is manually using problem.properties.dimension")
    end
    dimension = Val{problem.properties.dimension}
    use_forwarddiff = Val{problem.properties.use_forwarddiff}
    assemble!(problem, time, dimension, use_forwarddiff)
end


text:
# 1. 生成网络函数
function generateNet(net0)
    (Ntotal,k) = size(net0);
    net = zeros(Int64,Ntotal,k);
    for hang=1:Ntotal
        for lie=1:k
            hl = convert(Int64, net0[hang,lie]);
            net[hang,lie] = hl;
        end
    end
    return net
end

# 2. 网络初始化函数
function initNet(net)
    (Ntotal,k) = size(net0);
    S = net;
    strategy = zeros(Int64,Ntotal,1);
    a = rand(1:Ntotal,1);
    strategy[a] = 1;
    return S,strategy
end

# 3. 更新策略函数
function updateStr(S,strategy,death_num,b)

    w = 0.01;
    (Ntotal,k) = size(S);

    # 死亡节点的邻居序列
    death_neigh = S[death_num,:];
    death_neigh_strategy = zeros(Int64,1,4);
    # 死亡节点邻居的策略序列
    death_neigh_strategy = strategy[death_neigh];
    sum_dnei_str = sum(death_neigh_strategy);

    # 更新策略

    # 如果死亡邻居全为D或全为C，则学习该策略
    if sum_dnei_str == 0
        new_strategy = 0;
    elseif sum_dnei_str == k
        new_strategy = 1;
    else
        # 说明死亡节点邻居中既有C也有D，要分别计算邻居的收益
        # 计算死亡节点邻居的收益，就需要找出死亡节点邻居的邻居
        death_neigh_Payoff=zeros(1,k);
        for nn=1:k
            neighNow = death_neigh[nn];
            death_nei_nei = S[neighNow,:];
            # 计算出邻居的邻居的策略序列
            death_nei_nei_strategyC =  sum(strategy[death_nei_nei]);
            # 如果邻居节点i为C策略
            if strategy[neighNow] == 1
                death_neigh_Payoff[nn] = -k * 1 + b * death_nei_nei_strategyC;
            else
                death_neigh_Payoff[nn] = b * death_nei_nei_strategyC;
            end
        end
        death_neigh_fitness = 1 - w + w * death_neigh_Payoff;
        fitc = sum(death_neigh_fitness .* death_neigh_strategy);
        fitTotal = sum(death_neigh_fitness);
        Pc=fitc / fitTotal;
        probC = rand();
        if probC<=Pc
            new_strategy = 1;
        else
            new_strategy = 0;
        end

    end

    return new_strategy

end

function Round(net,b)
    S,strategy = initNet(net); #println(strategy); ok2
    (Ntotal,k) = size(net);
    nc = 1;
    while ((nc>0) && (nc<Ntotal))
        death_num = rand(1:Ntotal,1);
        old_strategy = strategy[death_num][1];  # 原策略
        # 更新策略
        new_strategy = updateStr(S,strategy,death_num,b);
        strategy[death_num] = new_strategy;
        nc = nc + new_strategy - old_strategy;
    end
    return nc
end

function main(net0,Ncwin,b)
    net = generateNet(net0); #println(net); ok1
    (Ntotal,k) = size(net0);
    for i=1:10000
        @time nc = Round(net,b);
        if (nc == Ntotal)
            Ncwin += 1;
        end
    end
end

net0 = readdlm("N500K4.txt" );
w = 0.01;
c = 1;
Ncwin = 0;
b = 3.9;
@time main(net0,Ncwin,b);

max_stars_repo_name: dev-juliacn/sailing

max_issues_repo_name: dev-juliacn/sailing

max_forks_repo_name: dev-juliacn/sailing


text:
module LittleMC

# Because there is no package to do just what I want, this is
# a simple Metropolis-Hastings sampler to just get things done.
# inspired by github.com/dfm/emcee/mh.py

export MC, sample, start, runstats

using Distributions
using PDMats
using HDF5

# The main object that contains all of the necessary MCMC information
type MC{T <: AbstractMvNormal} #Could try making this immutable later
    f::Function # Function to evaluate, generally the log-posterior
    nsamples::Int
    p0::Vector{Float64} # Starting parameters
    proposal::T # Jump proposal covariance matrix
    nparams::Int
    samples::Matrix{Float64}
    lnprobs::Vector{Float64}
    naccepted::Int
    csv::IOStream
end

# Initialization function
function MC(f::Function, nsamples::Int, p0::Vector{Float64}, propcov::Matrix{Float64}, csv::IOStream)
    # Do some custom initialization
    nparams = length(p0)

    # Make sure we didn't screw up with matrices
    @assert(size(propcov)[1] == size(propcov)[2], "propcov must be square matrix")
    @assert(nparams == size(propcov)[1] , "propcov mismatched with nparams")

    samples = Array(Float64, nparams, nsamples)
    lnprobs = Array(Float64, nsamples)

    # Turn propcov matrix into a distribution
    proposal = MvNormal(propcov)

    MC(f, nsamples, p0, proposal, nparams, samples, lnprobs, 0, csv)
end

function sample(mc::MC, p0::Vector{Float64}, lnprob0::Float64)

    # Get parameter proposal
    q = p0 .+ rand(mc.proposal)

    # Evaluate the proposal
    lnprob_proposed = mc.f(q)

    diff = lnprob_proposed - lnprob0

    # M-H acceptance ratio
    if diff < 0
        diff = exp(diff) - rand()
    end

    if diff >= 0
        # proposal accepted
        mc.naccepted += 1
        return (q, lnprob_proposed)
    else
        # proposal rejected
        return (p0, lnprob0)
    end

    return (p, lnprob)
end

# Start running the MCMC
function start(mc::MC)
    # Evaluate at the current parameters
    p0 = mc.p0
    lnprob0 = mc.f(p0)

    # always accept this as a starting point
    mc.samples[:, 1] = p0
    mc.lnprobs[1] = lnprob0

    writecsv(mc.csv, p0')

    for i=2:mc.nsamples
        (p0, lnprob0) = sample(mc, p0, lnprob0)
        mc.samples[:, i] = p0
        mc.lnprobs[i] = lnprob0
        writecsv(mc.csv, p0')
        flush(mc.csv) # Make sure this gets written.
    end
end

# Calculate the Metropolis acceptance fraction
function acceptance(mc::MC)
    return mc.naccepted/mc.nsamples
end

# Calculate the covariance of the samples
# function covariance(mc::MC)
# end

# Wrapper function for all the analysis tools
function runstats(mc::MC)
    @printf("Acceptance fraction %.2f\n", acceptance(mc))
    # more to come
end

# Optionally write the samples to an HDF5 file
function write(mc::MC, fname::AbstractString)
    fid = h5open(fname, "w")
    fid["samples"] = mc.samples

    # Add the acceptance fraction as an attribute on the "samples" dataset
    attrs(fid["samples"])["acceptance"] = acceptance(mc)

    close(fid)
end

end #module

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leaflettilelayer

"""
    ''_leaflettilelayer(;kwargs...)

A LeafletTileLayer component.

Keyword arguments:
- `id` (String; optional)
- `attribution` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `opacity` (Real; optional)
- `url` (String; optional)
- `zIndex` (Real; optional)
"""
function ''_leaflettilelayer(; kwargs...)
        available_props = Symbol[:id, :attribution, :loading_state, :opacity, :url, :zIndex]
        wild_props = Symbol[]
        return Component("''_leaflettilelayer", "LeafletTileLayer", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
"""
    gradient!{T <: AbstractFloat}(
        output::AbstractArray,
        f::Function,
        x::AbstractArray{T},
        buffer::AbstractArray{T},
        ::ForwardMode,
    )

# Description

Evaluate the gradient of `f` at `x` using forward-mode finite differences.
Store the results into `output`. Work with a user-provided `buffer` to ensure
that we can work without copies, but also without mutating `x`. In Rust jargon,
we take ownership of both `output` and `buffer`, but do not require ownership
of `x`. We just need read-only access to `x`.

# Arguments

* `output::AbstractArray`: An array that will be mutated to contain the
    gradient of `f` at `x`.
* `f::Function`: The function to be differentiated.
* `x::AbstractArray{T}`: The value of `x` at which to evaluate the gradient of
    `f`.
* `buffer::AbstractArray{T}`: A buffer that is equivalent to `similar(x)`. Used
    as a temporary copy of `x` that can be mutated safely during the execution
    of the function.
* `::ForwardMode`: An instance of the `ForwardMode` type.

# Returns

* `nothing::Void`: This function is called for its side effects.

# Examples

```julia
import FiniteDiff: gradient!, ForwardMode
x = [0.0, 0.0]
output, buffer = Array(Float64, 2), Array(Float64, 2)
gradient!(
    output,
    x -> sin(x[1]) + 2 * sin(x[2]),
    x,
    buffer,
    ForwardMode(),
)
```
"""
function gradient!{T <: AbstractFloat}(
    output::AbstractArray,
    f::Function,
    x::AbstractArray{T},
    buffer::AbstractArray{T},
    ::ForwardMode,
)::Void
    # Validate that all inputs have the same number of dimensions.
    n = length(x)
    if length(output) != n || length(buffer) != n
        throw(DomainError())
    end

    # Cache the value of f(x) to avoid multiple calls to f, which is pure.
    f_x = f(x)

    # Copy x into the user-supplied buffer.
    copy!(buffer, x)

    # Iterate over the dimensions of the input.
    for i in 1:n
        # Make a copy of the true value of x[i].
        x_i = x[i]

        # Temporarily change buffer[i] to a new value.
        ϵ = step_size(ForwardMode(), x_i)
        buffer[i] = x_i + ϵ

        # Evaluate f at the new value of x stored in the buffer.
        f_xp = f(buffer)

        # Restore the true value of x[i] in the buffer.
        buffer[i] = x_i

        # Store the i-th component of the gradient into the output array.
        # TODO: Remove this
        # @printf("%d\t%s\t%s\t%s\t%s\n", i, f_xp, f_x, f_xp - f_x, ϵ)
        output[i] = (f_xp - f_x) / ϵ
    end

    # Return nothing. The result is already stored in output.
    return
end

"""
    gradient!{T <: AbstractFloat}(
        output::AbstractArray,
        f::Function,
        x::AbstractArray{T},
        buffer::AbstractArray{T},
        ::BackwardMode,
    )

# Description

Evaluate the gradient of `f` at `x` using backward-mode finite differences.
Store the results into `output`. Work with a user-provided `buffer` to ensure
that we can work without copies, but also without mutating `x`. In Rust jargon,
we take ownership of both `output` and `buffer`, but do not require ownership
of `x`. We just need read-only access to `x`.

# Arguments

* `output::AbstractArray`: An array that will be mutated to contain the
    gradient of `f` at `x`.
* `f::Function`: The function to be differentiated.
* `x::AbstractArray{T}`: The value of `x` at which to evaluate the gradient of
    `f`.
* `buffer::AbstractArray{T}`: A buffer that is equivalent to `similar(x)`. Used
    as a temporary copy of `x` that can be mutated safely during the execution
    of the function.
* `::BackwardMode`: An instance of the `BackwardMode` type.

# Returns

* `nothing::Void`: This function is called for its side effects.

# Examples

```julia
import FiniteDiff: gradient!, BackwardMode
x = [0.0, 0.0]
output, buffer = Array(Float64, 2), Array(Float64, 2)
gradient!(
    output,
    x -> sin(x[1]) + 2 * sin(x[2]),
    x,
    buffer,
    BackwardMode(),
)
```
"""
function gradient!{T <: AbstractFloat}(
    output::AbstractArray,
    f::Function,
    x::AbstractArray{T},
    buffer::AbstractArray{T},
    ::BackwardMode,
)::Void
    # Validate that all inputs have the same number of dimensions.
    n = length(x)
    if length(output) != n || length(buffer) != n
        throw(DomainError())
    end

    # Cache the value of f(x) to avoid multiple calls to f, which is pure.
    f_x = f(x)

    # Copy x into the user-supplied buffer.
    copy!(buffer, x)

    # Iterate over dimensions of the input.
    for i in 1:n
        # Make a copy of the true value of x[i].
        x_i = buffer[i]

        # Temporarily change buffer[i] to a new value.
        ϵ = step_size(BackwardMode(), x_i)
        buffer[i] = x_i - ϵ

        # Evaluate f at the new value of x stored in the buffer.
        f_xm = f(buffer)

        # Restore the true value of x[i] in the buffer.
        buffer[i] = x_i

        # Store the i-th component of the gradient into the output array.
        output[i] = (f_x - f_xm) / ϵ
    end

    # Return nothing. The result is already stored in output.
    return
end

"""
    gradient!{T <: AbstractFloat}(
        output::AbstractArray,
        f::Function,
        x::AbstractArray{T},
        buffer::AbstractArray{T},
        ::CentralMode,
    )

# Description

Evaluate the gradient of `f` at `x` using central-mode finite differences.
Store the results into `output`. Work with a user-provided `buffer` to ensure
that we can work without copies, but also without mutating `x`. In Rust jargon,
we take ownership of both `output` and `buffer`, but do not require ownership
of `x`. We just need read-only access to `x`.

# Arguments

* `output::AbstractArray`: An array that will be mutated to contain the
    gradient of `f` at `x`.
* `f::Function`: The function to be differentiated.
* `x::AbstractArray{T}`: The value of `x` at which to evaluate the gradient of
    `f`.
* `buffer::AbstractArray{T}`: A buffer that is equivalent to `similar(x)`. Used
    as a temporary copy of `x` that can be mutated safely during the execution
    of the function.
* `::CentralMode`: An instance of the `CentralMode` type.

# Returns

* `nothing::Void`: This function is called for its side effects.

# Examples

```julia
import FiniteDiff: gradient!, CentralMode
x = [0.0, 0.0]
output, buffer = Array(Float64, 2), Array(Float64, 2)
gradient!(
    output,
    x -> sin(x[1]) + 2 * sin(x[2]),
    x,
    buffer,
    CentralMode(),
)
```
"""
function gradient!{T <: AbstractFloat}(
    output::AbstractArray,
    f::Function,
    x::AbstractArray{T},
    buffer::AbstractArray{T},
    ::CentralMode,
)::Void
    # Validate that all inputs have the same number of dimensions.
    n = length(x)
    if length(output) != n || length(buffer) != n
        throw(DomainError())
    end

    # Cache the value of f(x) to avoid multiple calls to f, which is pure.
    f_x = f(x)

    # Copy x into the user-supplied buffer.
    copy!(buffer, x)

    # Iterate over dimensions of the input.
    for i in 1:n
        # Make a copy of the true value of x[i].
        x_i = buffer[i]

        # Temporarily change buffer[i] to a new value.
        ϵ = step_size(CentralMode(), x_i)
        buffer[i] = x_i + ϵ

        # Evaluate f at the new value of x stored in the buffer.
        f_xp = f(buffer)

        # Temporarily change buffer[i] to a new value.
        buffer[i] = x_i - ϵ

        # Evaluate f at the new value of x stored in the buffer.
        f_xm = f(buffer)

        # Restore the true value of x[i] in the buffer.
        buffer[i] = x_i

        # Store the i-th component of the gradient into the output array.
        # TODO: Remove this
        # @printf("%d\t%s\t%s\t%s\t%s\n", i, f_xp, f_xm, f_xp - f_xm, ϵ)
        output[i] = (f_xp - f_xm) / (ϵ + ϵ)
    end

    # Return nothing. The result is already stored in output.
    return
end

"""
    gradient!{T <: AbstractFloat}(
        output::AbstractArray,
        f::Function,
        x::AbstractArray{T},
        buffer::AbstractArray{Complex{T}},
        ::ComplexMode,
    )

# Description

Evaluate the gradient of `f` at `x` using complex-mode finite differences.
Store the results into `output`. Work with a user-provided `buffer` to ensure
that we can work without copies, but also without mutating `x`. In Rust jargon,
we take ownership of both `output` and `buffer`, but do not require ownership
of `x`. We just need read-only access to `x`.

**NOTE**: The buffer must have complex elements to allow the computation of
complex-mode finite differences.

**NOTE**: This mode of finite differences will work correctly only when:

* `f` supports complex inputs.
* `f` is an analytic function in the complex analysis sense of the word.

# Arguments

* `output::AbstractArray`: An array that will be mutated to contain the
    gradient of `f` at `x`.
* `f::Function`: The function to be differentiated.
* `x::AbstractArray{T}`: The value of `x` at which to evaluate the gradient of
    `f`.
* `buffer::AbstractArray{Complex{T}}`: A buffer that is equivalent to
    `similar(x, Complex{eltype(x)})`. Used as a temporary copy of `x` that can
    be mutated safely during the execution of the function.
* `::ComplexMode`: An instance of the `ComplexMode` type.

# Returns

* `nothing::Void`: This function is called for its side effects.

# Examples

```julia
import FiniteDiff: gradient!, ComplexMode
x = [0.0, 0.0]
output, buffer = Array(Float64, 2), Array(Complex{Float64}, 2)
gradient!(
    output,
    x -> sin(x[1]) + 2 * sin(x[2]),
    x,
    buffer,
    ComplexMode(),
)
```
"""
function gradient!{T <: AbstractFloat}(
    output::AbstractArray,
    f::Function,
    x::AbstractArray{T},
    buffer::AbstractArray{Complex{T}},
    ::ComplexMode,
)::Void
    # Validate that all inputs have the same number of dimensions.
    n = length(x)
    if length(output) != n || length(buffer) != n
        throw(DomainError())
    end

    # Copy x into the user-supplied buffer.
    copy!(buffer, x)

    # Iterate over dimensions of the input.
    for i in 1:n
        # Make a copy of the true value of x[i].
        x_i = x[i]

        # Temporarily change buffer[i] to a new value.
        ϵ = step_size(ComplexMode(), x_i)
        buffer[i] = x_i + ϵ * im

        # Evaluate f at the new value of x stored in the buffer.
        f_xim = f(buffer)

        # Restore the true value of x[i] in the buffer.
        buffer[i] = x_i

        # Store the i-th component of the gradient into the output array.
        output[i] = imag(f_xim) / ϵ
    end

    # Return nothing. The result is already stored in output.
    return
end

"""
    gradient!{T <: AbstractFloat}(
        output::AbstractArray,
        f::Function,
        x::AbstractArray{T},
        buffer::AbstractArray{T},
    )

# Description

Evaluate the gradient of `f` at `x` using finite differences. Defaults to
central mode finite differences. See the documentation for that mode for
more details.

# Arguments

* `output::AbstractArray`: An array that will be mutated to contain the
    gradient of `f` at `x`.
* `f::Function`: The function to be differentiated.
* `x::AbstractArray{T}`: The value of `x` at which to evaluate the gradient of
    `f`.
* `buffer::AbstractArray{T}`: A buffer that is equivalent to `similar(x)`. Used
    as a temporary copy of `x` that can be mutated safely during the execution
    of the function.

# Returns

* `nothing::Void`: This function is called for its side effects.

# Examples

```julia
import FiniteDiff: gradient!
x = [0.0, 0.0]
output, buffer = Array(Float64, 2), Array(Float64, 2)
gradient!(
    output,
    x -> sin(x[1]) + 2 * sin(x[2]),
    x,
    buffer,
)
```
"""
function gradient!{T <: AbstractFloat}(
    output::AbstractArray,
    f::Function,
    x::AbstractArray{T},
    buffer::AbstractArray{T},
)::Void
    gradient!(output, f, x, buffer, CentralMode())
    return
end

"""
    gradient{T <: AbstractFloat}(
        f::Function,
        x::AbstractArray{T},
        mode::Mode = CentralMode(),
    )

# Description

Evaluate the gradient of `f` at `x` using finite differences. Defaults to
central mode finite differences. See the documentation for that mode for
more details.

**NOTE**: Allocates new memory to store the output as well as memory to use as
a buffer.

# Arguments

* `f::Function`: The function to be differentiated.
* `x::AbstractArray{T}`: The value of `x` at which to evaluate the gradient of
    `f`.
* `m::Mode`: An instance of the `Mode` type. This will determine the mode
    of finite differences that will be used.

# Returns

* `output::AbstractArray`: The gradient.

# Examples

```julia
import FiniteDiff: gradient!, ForwardMode
x = [0.0, 0.0]
output = gradient!(x -> sin(x[1]) + 2 * sin(x[2]), x, ForwardMode())
```
"""
function gradient{T <: AbstractFloat}(
    f::Function,
    x::AbstractArray{T},
    mode::Mode = CentralMode(),
)::AbstractArray{T}
    if mode != ComplexMode()
        output, buffer = similar(x), similar(x)
    else
        output, buffer = similar(x), similar(x, Complex{T})
    end
    gradient!(output, f, x, buffer, mode)
    return output
end

"""
    gradient(f::Function, mode::Mode = CentralMode(); mutates::Bool = false)

# Description

Construct a new function that will evaluate the gradient of `f` at any value
of `x`. The user can specify the mode as a positional argument and can also
indicate whether to return a mutating or non-mutating function using a keyword
argument.

# Arguments

* `f::Function`: The function to be differentiated.
* `m::Mode`: An instance of the `Mode` type. This will determine the mode
    of finite differences that will be used.

# Keyword Arguments

* `mutates::Bool = false`: Determine whether the resulting function will mutate
    its inputs or will be a pure function.

# Returns

* `g′::Function`: A function to evaluate the gradient of `f` at a new point
    `x`.

# Examples

```julia
import FiniteDiff: gradient!, ForwardMode
x = [0.0, 0.0]
g′ = gradient(x -> sin(x[1]) + 2 * sin(x[2]), ForwardMode(), mutates = false)
output = g′(x)
```
"""
function gradient(
    f::Function,
    mode::Mode = CentralMode();
    mutates::Bool = false,
)::Function
    if mutates
        return (output, x, buffer) -> gradient!(output, f, x, buffer, mode)
    else
        return (x, ) -> gradient(f, x, mode)
    end
end


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent
using Random

##==----------------------------------------------------------------------------------------------------------------------

isunitary(g::AbstractGate) = Qaintessent.matrix(g) * Qaintessent.matrix(Base.adjoint(g)) ≈ I

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "gate helper functions" begin
    @testset "kron" begin
        @test kron(I,Z) ≈  kron(Matrix(I, (2,2)), matrix(Z))
    end

    @testset "gate addition rotation gates" begin
        for g in [RxGate, RyGate, RzGate, PhaseShiftGate, EntanglementXXGate, EntanglementYYGate, EntanglementZZGate]
            θ = 2π*rand()
            ϕ = 2π*rand()

            g1 = g(θ)
            g2 = g(ϕ)

            @test g1 + g2 ≈ g(θ+ϕ)
        end
    end

    @testset "gate addition general rotation gate" begin

        θ = 2π*rand()
        m = rand(3)
        m = m./norm(m)

        ϕ = 2π*rand()
        n = rand(3)
        n = n./norm(n)

        g1 = RotationGate(θ, m)
        g2 = RotationGate(ϕ, n)

        @test g1 + g2 ≈ RotationGate(m*θ + ϕ*n)
    end

    @testset "gate phase multiplication rotation gates" begin
        for g in [RxGate, RyGate, RzGate, PhaseShiftGate, EntanglementXXGate, EntanglementYYGate, EntanglementZZGate]
            θ = 2π*rand()
            α = rand()

            g1 = g(θ)

            @test α*g1 ≈ g(θ*α)
        end
    end

    @testset "gate multiplication general rotation gate" begin

        θ = 2π*rand()
        m = rand(3)
        m = m./norm(m)

        α = rand()

        g1 = RotationGate(θ, m)

        @test α*g1 ≈ RotationGate(m*θ*α)
    end
end

@testset ExtendedTestSet "quantum gates" begin

    θ = 2π*rand()
    ϕ = 2π*rand()
    n = randn(3); n /= norm(n)

    @testset "basic quantum gates" begin
        for g in [X, Y, Z, HadamardGate(), SGate(), SdagGate(), TGate(), TdagGate(), RxGate(θ), RyGate(θ), RzGate(θ), RotationGate(θ, n), PhaseShiftGate(ϕ), SwapGate(), EntanglementXXGate(θ), EntanglementYYGate(θ), EntanglementZZGate(θ), controlled_not()]
            @test isunitary(g)
            gdag = adjoint(g)

            @test Qaintessent.matrix(gdag) == adjoint(Qaintessent.matrix(g))
            @test LinearAlgebra.ishermitian(g) == (Qaintessent.matrix(gdag) == Qaintessent.matrix(g))
            @test Qaintessent.sparse_matrix(g) == Qaintessent.matrix(g)
        end
    end

    @testset "general rotation gate" begin
        @test Qaintessent.matrix(RotationGate(θ, [1, 0, 0])) ≈ Qaintessent.matrix(RxGate(θ))
        @test Qaintessent.matrix(RotationGate(θ, [0, 1, 0])) ≈ Qaintessent.matrix(RyGate(θ))
        @test Qaintessent.matrix(RotationGate(θ, [0, 0, 1])) ≈ Qaintessent.matrix(RzGate(θ))

        @test RotationGate(θ * n) ≈ RotationGate(θ, n)
        @test Qaintessent.matrix(RotationGate(0.3*θ, n)) * Qaintessent.matrix(RotationGate(0.7*θ, n)) ≈ Qaintessent.matrix(RotationGate(θ, n))
        @test Qaintessent.sparse_matrix(RotationGate(θ, n)) ≈ Qaintessent.matrix(RotationGate(θ, n))
    end

    @testset "general rotation gate exceptions" begin
        @test_throws ErrorException("Rotation axis vector must have length 3.") RotationGate(θ, [1, 0, 0, 0])
        @test_throws ErrorException("Norm of rotation axis vector must be 1.")  RotationGate(θ, [1, 2, 0])
        @test_throws ErrorException("Rotation axis vector must have length 3.") RotationGate([0, θ, 0, 1])
    end

    @testset "general unitary gate" begin

        # test MatrixGate
        N = 3
        d = 2
        A = randn(ComplexF64, d^N, d^N)
        U, _ = qr(A)
        U = Array(U)
        gateU = MatrixGate(U)
        gateUdag = adjoint(gateU)
        @test Qaintessent.matrix(gateU) ≈ U
        @test isunitary(gateU)
        @test LinearAlgebra.ishermitian(gateU) == (Qaintessent.matrix(gateUdag) == Qaintessent.matrix(gateU))
        @test LinearAlgebra.ishermitian(gateU) == (Qaintessent.sparse_matrix(gateUdag) == Qaintessent.sparse_matrix(gateU))

        U = Qaintessent.sparse_matrix(X)
        gateU = MatrixGate(U)
        gateUdag = adjoint(gateU)
        @test Qaintessent.sparse_matrix(gateU) ≈ U
        @test isunitary(gateU)
        @test true == (Qaintessent.sparse_matrix(gateUdag) == Qaintessent.sparse_matrix(gateU))
    end

    # not unitary
    @testset "general unitary gate exceptions" begin

        N = 3
        d = 2
        @test_throws ErrorException("Quantum gate must be unitary") MatrixGate(randn(ComplexF64, d^N, d^N))
    end

    @testset "general controlled gate" begin
        # test MatrixGate
        N = 4
        M = 2
        d = 2
        A = randn(ComplexF64, d^M, d^M)
        U, _ = qr(A)
        U = Array(U)
        gateU = MatrixGate(U)
        cgateU = ControlledGate(gateU, M)

        cgateUdag = adjoint(cgateU)
        @test LinearAlgebra.ishermitian(cgateU) == (Qaintessent.sparse_matrix(cgateU) == Qaintessent.sparse_matrix(cgateUdag))
        @test LinearAlgebra.ishermitian(cgateU) == (Qaintessent.matrix(cgateU) == Qaintessent.matrix(cgateUdag))
    end
end

@testset ExtendedTestSet "abstract gate tests" begin
    @test_throws Exception adjoint(AbstractGate())
    @test_throws Exception num_wires(AbstractGate())
    @test_throws Exception matrix(AbstractGate())
    @test_throws Exception sparse_matrix(AbstractGate())
end


text:
using OrdinaryDiffEq, DiffEqDevTools, ParameterizedFunctions, Plots, ODE, ODEInterfaceDiffEq, LSODA, Sundials
gr() #gr(fmt=:png)
using LinearAlgebra

f = @ode_def Orego begin
  dy1 = p1*(y2+y1*(1-p2*y1-y2))
  dy2 = (y3-(1+y1)*y2)/p1
  dy3 = p3*(y1-y3)
end p1 p2 p3

p = [77.27,8.375e-6,0.161]
prob = ODEProblem(f,[1.0,2.0,3.0],(0.0,30.0),p)
sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
abstols = 1.0 ./ 10.0 .^ (4:11)
reltols = 1.0 ./ 10.0 .^ (1:8);

plot_prob = ODEProblem(f,[1.0,2.0,3.0],(0.0,400.0),p)
sol = solve(plot_prob,CVODE_BDF())
plot(sol,yscale=:log10)

#sol = solve(prob,ode23s()); println("Total ODE.jl steps: $(length(sol))")
#using GeometricIntegratorsDiffEq
#try
#    sol = solve(prob,GIRadIIA3(),dt=1/10)
#catch e
#    println(e)
#end

sol = solve(prob,ARKODE(),abstol=1e-5,reltol=1e-1);

sol = solve(prob,ARKODE(nonlinear_convergence_coefficient = 1e-3),abstol=1e-5,reltol=1e-1);

sol = solve(prob,ARKODE(order=3),abstol=1e-5,reltol=1e-1);

sol = solve(prob,ARKODE(order=3,nonlinear_convergence_coefficient = 1e-5),abstol=1e-5,reltol=1e-1);

sol = solve(prob,ARKODE(order=5),abstol=1e-5,reltol=1e-1);

setups = [
          #Dict(:alg=>ROCK2())    #Unstable
          #Dict(:alg=>ROCK4())    #needs more iterations
          #Dict(:alg=>ESERK5()),
          ]

sol = solve(prob,EXPRB53s3(),dt=2.0^(-8));
sol = solve(prob,EPIRK4s3B(),dt=2.0^(-8));
sol = solve(prob,EPIRK5P2(),dt=2.0^(-8));

sol = solve(prob,PDIRK44(),dt=2.0^(-8));

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);
setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>FBDF()),
          Dict(:alg=>QNDF()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>rodas()),
          Dict(:alg=>radau()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>ROS34PW1a()),
          Dict(:alg=>lsoda()),
          ]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2,numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2,numruns=10)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Kvaerno3()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>KenCarp3()),
          Dict(:alg=>lsoda()),
    # Dict(:alg=>SDIRK2()), # Removed because it's bad
          Dict(:alg=>radau())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose = false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2,numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2,numruns=10)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>KenCarp5()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>KenCarp3()),
          Dict(:alg=>ARKODE(order=5)),
          Dict(:alg=>ARKODE(nonlinear_convergence_coefficient = 1e-6)),
          Dict(:alg=>ARKODE(nonlinear_convergence_coefficient = 1e-5,order=3))
          ]
names = ["Rosenbrock23" "KenCarp5" "KenCarp4" "KenCarp3" "ARKODE5" "ARKODE4" "ARKODE3"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>ImplicitEulerExtrapolation()),
          #Dict(:alg=>ImplicitDeuflhardExtrapolation()), # Diverges
          #Dict(:alg=>ImplicitHairerWannerExtrapolation()), # Diverges
          Dict(:alg=>ABDF2()),
          Dict(:alg=>FBDF()),
          Dict(:alg=>QNDF()),
          Dict(:alg=>Exprb43()),
          Dict(:alg=>Exprb32()),
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:13)
reltols = 1.0 ./ 10.0 .^ (4:10)

setups = [
          Dict(:alg=>FBDF()),
          Dict(:alg=>QNDF()),
          Dict(:alg=>Rodas4P()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>ddebdf()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>rodas()),
          Dict(:alg=>radau()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>lsoda()),
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2,numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2,numruns=10)
plot(wp)

setups = [Dict(:alg=>GRK4A()),
          Dict(:alg=>Rodas5()),
          Dict(:alg=>Kvaerno5()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>KenCarp5()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>lsoda()),
          Dict(:alg=>radau())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2,numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2,numruns=10)
plot(wp)

#setups = [Dict(:alg=>Hairer4()),
          #Dict(:alg=>Hairer42()),
          #Dict(:alg=>Rodas3()),
          #Dict(:alg=>Kvaerno4()),
          #Dict(:alg=>Cash4())
#]
#wp = WorkPrecisionSet(prob,abstols,reltols,setups;
#                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
#plot(wp)

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
using Base.Cartesian
using BenchmarkTools

@generated function unpack_n(v::Array{T, N}) where {M, U, T <: NTuple{M, U}, N}
    quote
        @nexprs $M i -> out_i = similar(v,U)
        @inbounds @nloops $N j v begin
            @nexprs $M i -> ((@nref $N out_i j) = (@nref $N v j)[i])
        end
        return @ntuple $M out
    end
end

@generated function unpack_t(v::Array{T,N}) where {T <: Tuple, N}
    TT = T.types
    M = length(TT)
    quote
        @nexprs $M i -> out_i = similar(v,$TT[i])
        @inbounds @nloops $N j v begin
            @nexprs $M i -> ((@nref $N out_i j) = (@nref $N v j)[i])
        end
        return @ntuple $M out
    end
end

f(a, b) = a+b, a*b, a-b

for sz in [10 50 100 200]
    packed = f.(rand(sz,sz,1), rand(1,1,sz))
    r = (@belapsed unpack_t($packed))/(@belapsed unpack_n($packed))
    println("$sz^3 : $r")
end


text:
####################################################################################################
function displayIteration(i, funceval, residual, itlinear = 0)
	(i == 0) && println("\n Newton Iterations \n   Iterations      Func-count      f(x)      Linear-Iterations\n")
	if length(itlinear)==1
		@printf("%9d %16d %14.4e %9d\n", i, funceval, residual, itlinear);
	else
		if itlinear isa Tuple{Int64,Int64}
			@printf("%9d %16d %14.4e (%6d, %6d)\n", i, funceval, residual, itlinear[1], itlinear[2]);
		else
			# used for nested linear solves
			@printf("%9d %16d %14.4e  ", i, funceval, residual);println(itlinear);
		end
	end
end
####################################################################################################
function computeEigenvalues(contparams::ContinuationPar, n_unstable::Int64, J)
	# the next line is to ensure we compute enough eigenvalues to probe stability
	nev_ = max(n_unstable + 5, contparams.nev)
	eiginfo = contparams.newtonOptions.eigsolver(J, nev_)
	isstable, n_unstable, n_imag = is_stable(contparams, eiginfo[1])
	return eiginfo, isstable, n_unstable, n_imag
end

computeEigenvalues(iter::PALCIterable, state::PALCStateVariables) = computeEigenvalues(iter.contParams, state.n_unstable[2], iter.J(getu(state), getp(state)))

function computeEigenvalues(iter::PALCIterable, state::PALCStateVariables)
	n_unstable = state.n_unstable[2]
	nev_ = max(n_unstable + 5, iter.contParams.nev)
	J = iter.J(getu(state), getp(state))
	eiginfo = iter.contParams.newtonOptions.eigsolver(J, nev_)
	isstable, n_unstable, n_imag = is_stable(iter.contParams, eiginfo[1])
	return eiginfo, isstable, n_unstable, n_imag
end

function computeEigenvalues!(iter::PALCIterable, state::PALCStateVariables)
	# we compute the eigenelements
	n_unstable = state.n_unstable[2]
	nev_ = max(n_unstable + 5, iter.contParams.nev)
	J = iter.J(getu(state), getp(state))
	eigvals, eigvecs, flag, it_number = iter.contParams.newtonOptions.eigsolver(J, nev_)
	# we update the state
	isstable, n_unstable, n_imag = is_stable(iter.contParams, eigvals)
	updatestability!(state, n_unstable, n_imag)
	if isnothing(state.eigvals) == false
		state.eigvals = eigvals
	end
	if iter.contParams.saveEigenvectors
		state.eigvecs = eigvecs
	end
	return it_number
end

####################################################################################################
"""
Compute a Jacobian by Finite Differences
"""
function finiteDifferences(F, x::AbstractVector; δ = 1e-9)
	N = length(x)
	Nf = length(F(x))
	J = zeros(eltype(x), Nf, N)
	finiteDifferences!(F, J, x; δ = δ)
	return J
end

function finiteDifferences!(F, J, x::AbstractVector; δ = 1e-9)
	f = F(x)
	x1 = copy(x)
	for i ∈ eachindex(x)
		x1[i] += δ
		J[:, i] .= (F(x1) .- F(x)) / δ
		x1[i] -= δ
	end
	return J
end

####################################################################################################
"""
Save solution / data in JLD2 file
- `filename` is for example "example.jld2"
- `sol` is the solution
- `p` is the parameter
- `i` is the index of the solution to be saved
"""
function saveToFile(filename, sol, p, i::Int64, br::ContResult, contParam::ContinuationPar)
	# create a group in the JLD format
	jldopen(filename*".jld2", "a+") do file
		mygroup = JLD2.Group(file, "sol-$i")
		mygroup["sol"] = sol
		mygroup["param"] = p
	end

	jldopen(filename*"-branch.jld2", "w") do file
		file["branch"] = br
		file["contParam"] = contParam
	end
end
####################################################################################################
using BlockArrays, SparseArrays
"""
Function waiting to be accepted to BlockArrays.jl
"""
function blockToSparse(J::AbstractBlockArray)
	nl, nc = size(J.blocks)
	# form the first line of blocks
	res = J[Block(1,1)]
	for j=2:nc
		res = hcat(res, J[Block(1,j)])
	end
	# continue with the other lines
	for i=2:nl
		line = J[Block(i,1)]
		for j=2:nc
			line = hcat(line, J[Block(i,j)])
		end
		res = vcat(res,line)
	end
	return res
end
####################################################################################################
# this trick is extracted from KrylovKit. It allows for the Jacobian to be specified as a matrix (sparse / dense) or as a function.
apply(A::AbstractMatrix, x::AbstractVector) = A * x
apply(f, x) = f(x)

# apply!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
# apply!(y, f, x) = copyto!(y, f(x))

# empty eigenvectors to save memory
_empty(a::AbstractVector{T}, ::Type{U}=T) where {T,U} = Vector{U}()
_empty(a::AbstractMatrix{T}, ::Type{U}=T) where {T,U} = similar(a, (0,0))


text:
"""
    PyCode(code::String, filename::String, mode::Symbol)

A Python code object, representing the compiled contents of `code`.

The `filename` is used for exception printing. The mode must be `:exec` or `:eval`.

See also [`@py_cmd`](@ref) and [`@pyv_cmd`](@ref).
"""
mutable struct PyCode
    ref::PyRef
    code::String
    filename::String
    mode::Symbol
    PyCode(code::String, filename::String, mode::Symbol) = begin
        mode in (:exec, :eval) || error("invalid mode $(repr(mode))")
        new(PyRef(), code, filename, mode)
    end
end
export PyCode

ispyreftype(::Type{PyCode}) = true
pyptr(co::PyCode) = begin
    ptr = co.ref.ptr
    if isnull(ptr)
        ptr =
            co.ref.ptr = C.Py_CompileString(
                co.code,
                co.filename,
                co.mode == :exec ? C.Py_file_input :
                co.mode == :eval ? C.Py_eval_input : error("invalid mode $(repr(co.mode))"),
            )
    end
    ptr
end
Base.unsafe_convert(::Type{CPyPtr}, x::PyCode) = checknull(pyptr(x))


text:
function mdl_source!(voltage_r,
                     voltage_i,
                        current_r,
                        current_i,
                     device::StaticSource,
                     sys::DynamicSystem)

        #Load device parameters
        V_R = device.V_R
        V_I = device.V_I
        X_th = device.X_th

        #I = ( (V_R + V_I*1im) - (V_tR + V_tI*1im) )/(X_th*1im)
        #I_r = (V_I - V_tI)/X_th #in system pu flowing out
        #I_i = -(V_R - V_tR)/X_th #in system pu flowing out

        #update current
        current_r[1] += (V_I - voltage_i[1])/X_th #in system pu flowing out
        current_i[1] += -(V_R - voltage_r[1])/X_th #in system pu flowing out

        return
end