text:
@testset "Specifying Initial state - Continuous Wave" begin
	# ----------------------------------------
	# Manually specifying the initial state for ContinuousWave
  problem = WaveguideProblem(TwoLevelChain(1), ContinuousWave(1.0), 10.0)
  solA = solve(problem, reltol = 1e-10)[2]

  # Project on the Rydberg manifold
  σ_Ryd = transition(NLevelBasis(2), 2, 2)

	ρ₀ = nlevelstate(NLevelBasis(2), 1)
  solB = solve(problem, ρ₀ = ρ₀, reltol = 1e-10)[2]

  @test sum( tracedistance.(solA, solB) ) < 1e-4
end

@testset "Specifying Initial state - Fock" begin
	# ----------------------------------------
	# Variant of "Systems/Perfect Absorption" test
  problem = WaveguideProblem(TwoLevelChain(1), WavePacket(ExpMode(), Fock(1)), (-10.0, 0.0))
  solA = solve(problem, reltol = 1e-10)[2]

  # Project on the Rydberg manifold
  σ_Ryd = transition(NLevelBasis(2), 2, 2)

	ρ₀ = fockstate(FockBasis(1), 1) ⊗ nlevelstate(NLevelBasis(2), 1)
  solB = solve(problem, ρ₀ = ρ₀, reltol = 1e-10)[2]

  @test sum( tracedistance.(solA, solB) ) < 1e-4
end

@testset "Overwriting Nouts" begin
	problem = WaveguideProblem(TwoLevelChain(1), WavePacket(ExpMode(), Fock(1)), ExpMode(), 10.0)
	basisF, basisA = FockBasis(1), NLevelBasis(2)
	ρ₀ = fockstate(basisF, 1) ⊗ nlevelstate(basisA, 1) ⊗ fockstate(basisF, 0)
	@test_logs (:warn, "`Nouts` and `ρ₀` simultaneously set. Ignoring `Nouts`") begin
		solve(problem, ρ₀ = ρ₀, Nouts = [3])
		true
	end
end


text:
function shock_loading(m::HetDSGE, TTT_jump::Matrix{Float64})
    exo = m.exogenous_shocks
    endo = m.endogenous_states

    _RRR = zeros(n_backward_looking_states(m), n_shocks_exogenous(m))
    RRR  = zeros(n_model_states(m), n_shocks_exogenous(m))

    _RRR[endo[:b′_t], exo[:b_sh]] .= 1.
    _RRR[endo[:g′_t], exo[:g_sh]] .= 1.
    _RRR[endo[:z′_t], exo[:z_sh]] .= 1.
    _RRR[endo[:μ′_t], exo[:μ_sh]] .= 1.
    _RRR[endo[:λ_w′_t], exo[:λ_w_sh]] .= 1.
    _RRR[endo[:λ_f′_t], exo[:λ_f_sh]] .= 1.
    _RRR[endo[:rm′_t], exo[:rm_sh]] .= 1.

    # Loading on states and jumps
    RRR[1:n_backward_looking_states(m), :] = _RRR
    RRR[n_backward_looking_states(m)+1:end, :] = TTT_jump*_RRR

    return RRR
end

ext: jl


text:
struct PersistenceDiagramRowIterator{D}
    diagram::D
end

function Base.iterate(it::PersistenceDiagramRowIterator, i=1)
    if i > length(it.diagram)
        return nothing
    else
        int = it.diagram[i]
        dim = get(it.diagram.meta, :dim, missing)
        threshold = get(it.diagram.meta, :threshold, missing)
        return (birth=int.birth, death=int.death, dim=dim, threshold=threshold), i + 1
    end
end

Base.IteratorSize(::PersistenceDiagramRowIterator) = Base.HasLength()
Base.length(it::PersistenceDiagramRowIterator) = length(it.diagram)

Tables.istable(::Type{<:PersistenceDiagram}) = true
Tables.rowaccess(::Type{<:PersistenceDiagram}) = true
function Tables.rows(diagram::PersistenceDiagram)
    return PersistenceDiagramRowIterator(diagram)
end
function Tables.schema(it::PersistenceDiagramRowIterator)
    diagram = it.diagram
    D = hasproperty(diagram, :dim) ? Int : Missing
    T = hasproperty(diagram, :threshold) ? Float64 : Missing
    return Tables.Schema((:birth, :death, :dim, :threshold), (Float64, Float64, D, T))
end

Tables.materializer(::PersistenceDiagram) = PersistenceDiagram


text:
handle_broadcast_name(names, sizes, finalsize) = map(zip(names, sizes)) do (name, sz)
    sz = (sz..., ntuple(() -> 1, length(finalsize) - length(sz))...)
    if sz == finalsize
        return suffix -> "$(name)$suffix"
    else
        return suffix -> begin
            outname = "$(name)$suffix"
            for (s, sf) in zip(sz, finalsize)
                outname = "{($sf / $s){$outname}}"
            end

            return outname
        end
    end
end

handle_broadcast_name(lhs::String, rhs::String, lsize, rsize) =
    (lsize == (1, 1)) && (rsize == (1, 1)) ? (suffix -> "$(lhs)$suffix", suffix -> "$(rhs)$suffix") :
    lsize == (1, 1) ? (suffix -> "{$(rsize[1] * rsize[2]){$(lhs)$suffix}}",
                       suffix -> "$(rhs)$suffix") :
    rsize == (1, 1) ? (suffix -> "$(lhs)$suffix",
                       suffix -> "{$(lsize[1] * lsize[2]){$(rhs)$suffix}}") :
    (suffix -> "$(lhs)$suffix", suffix -> "$(rhs)$suffix")

function write_bcast_instantiation(buffer, prefix, outsize, body)
    write(buffer, """
        genvar $(join(ntuple(i -> "$(prefix)_i_$i", length(outsize)), ", "));

        generate
        """)
    for (i, sz) in enumerate(outsize)
        write(buffer, repeat(" ", i - 1) * "for ($(prefix)_i_$i = 0; $(prefix)_i_$i < $sz; $(prefix)_i_$i = $(prefix)_i_$i + 1) begin : $(prefix)_gen_$i\n")
    end
    write(buffer, repeat(" ", length(outsize)) * "localparam $(prefix)_i = ")
    write(buffer, join(ntuple(i -> "$(prefix)_i_$i * $(join(outsize[(i + 1):end], "*"))", length(outsize) - 1), " + "))
    write(buffer, " + $(prefix)_i_$(length(outsize));\n")
    for line in split(body, "\n"; keepempty = true)
        write(buffer, repeat(" ", length(outsize)) * line * "\n")
    end
    for i in length(outsize):-1:1
        write(buffer, repeat(" ", i - 1) * "end\n")
    end
    write(buffer, "endgenerate\n")
end


text:
# script for sending and decoding requests for image computations

"""
    parse_npy_header(hdr::AbstractString)

Function to parse npy header information to extract dimensions and array dtypes.
Customized implementation of `NPY.parseheader` to handle numpy structured arrays which EE always returns.
See https://numpy.org/devdocs/reference/generated/numpy.lib.format.html#format-version-1-0 for spec.
"""
function parse_npy_header(hdr::AbstractString)
    # rework of the NPY.parseheader function to allow for structured arrays
    # https://numpy.org/doc/stable/user/basics.rec.html#structured-datatypes
    s = NPZ.parsechar(hdr, '{')

    dict = Dict{String,Any}()
    # loop through and get the three key info
    # description, fortran_order, and shape
    for _ = 1:3
        s = strip(s)
        key, s = NPZ.parsestring(s)
        s = strip(s)
        s = NPZ.parsechar(s, ':')
        s = strip(s)
        if key == "descr"
            x, s = split(s, ", 'for")
            # clean up the string info to be able to easily parse
            y = split(
                replace(
                    replace(
                        replace(
                            replace(replace(x, "[" => ""), "]" => ""),
                            "(" => "",
                        ),
                        ")" => "",
                    ),
                    "'" => "",
                ),
                ",",
            )

            bnames = []
            dtypes = []

            # loop through the band info to get datatypes
            for z in y
                if occursin("<", z)
                    push!(dtypes, replace(strip(z), "<" => ""))
                elseif occursin(">", z)
                    push!(dtypes, replace(strip(z), ">" => ""))
                else
                    push!(bnames, strip(z))
                end
            end

            dict[key] = Dict(bnames .=> dtypes)

            # kinda hacky way of reconstructing the header info after splitting...
            s = ", 'for$(s)"

        elseif key == "fortran_order"
            dict[key], s = NPZ.parsebool(s)
        elseif key == "shape"
            dict[key], s = NPZ.parsetuple(s)
        else
            error("parsing header failed: bad dictionary key")
        end
        s = strip(s)
        if s[firstindex(s)] == '}'
            break
        end
        s = NPZ.parsechar(s, ',')
    end
    return dict
end

"""
    decode_npy_response(response::PyObject)

Function to take a response object of a npy bytes and convert to Julia array
"""
function decode_npy_response(response::PyObject)
    foo = response.content

    hdr = parse_npy_header(strip(split(foo, "\n")[1])[11:end])

    banddtypes = map(x -> NPZ.Numpy2Julia[x], values(hdr["descr"]))

    sametyped = all(x -> x == banddtypes[1], banddtypes)

    if ~sametyped
        @error(
            "Resulting bands have different data types, creating arrays with different band dtypes is currently not support ...Try explicitly casting the image to the same data type (ie `ee.Image().toFloat()`)."
        )
    end

    dtype = banddtypes[1]

    # create an empty in-memory buffer to temporarily write data to
    buffer = IOBuffer()
    write(buffer, foo)
    # convert data to correct dtype
    bar = reinterpret(dtype, take!(buffer))
    close(buffer) # close the io buffer for clean memory management

    # get the data shape of the bytes
    height, width = hdr["shape"]
    channels = length(banddtypes)

    # find the number of pixels we expect based on img
    data_size::Int64 = width * height * channels
    # find the offset of what was returned vs what is expected
    offset::Int64 = size(bar, 1) - data_size + 1

    arr::Array{dtype,1} = dtype.(bar)[offset:end]

    if channels > 1
        # reshape array from 1-d to correct 3-d dimensions
        arr3d::Array{dtype,3} = reshape(arr, (channels, width, height))
        # reorder the dimensions so it is W,H,C
        img = permutedims(arr3d, (2, 3, 1))
    else
        arr2d::Array{dtype,2} = reshape(arr, (width, height))
        img = arr2d #rotl90(arr2d)
    end

    return img
end

"""
    computepixels(session::EESession, pixelgrid::PixelGrid, image::EE.AbstractEEObject; format::AbstractString = "NPY")

Function to take an EarthEngine computed image and return an Array with geographic information (i.e. GeoArray).
This signature will return all of the bands within the image.
Currently on the "NPY" format is available.
See https://developers.google.com/earth-engine/reference/rest/v1beta/projects.image/computePixels
"""
function computepixels(
    session::EESession,
    pixelgrid::PixelGrid,
    image::EE.AbstractEEObject;
    format::AbstractString = "NPY",
)
    # get the band information
    bands = computevalue(session, bandNames(image))

    # pass to computepixels with band info
    computepixels(session, pixelgrid, image, bands; format = format)

end

"""
    computepixels(session::EESession, pixelgrid::PixelGrid, image::EE.AbstractEEObject, bands::AbstractVector{String}; format::AbstractString = "NPY")

Function to take an EarthEngine computed image and return an Array with geographic information (i.e. GeoArray).
Currently on the "NPY" format is available.
See https://developers.google.com/earth-engine/reference/rest/v1beta/projects.image/computePixels
"""
function computepixels(
    session::EESession,
    pixelgrid::PixelGrid,
    image::EE.AbstractEEObject,
    bands::AbstractVector{String};
    format::AbstractString = "NPY",
)
    endpoint = "image:computePixels"

    serialized = ee.serializer.encode(image, for_cloud_api = true)

    payload = Dict(
        :expression => serialized,
        :grid => typetodict(pixelgrid),
        :fileFormat => format,
        :bandIds => bands,
    )

    response = _sendrequest(session, endpoint, payload)

    result = decode_npy_response(response)

    affine = extract_affinemap(pixelgrid)

    output = GeoArray(result[:,end:-1:1], affine, pixelgrid.crsCode)

    return output
end


text:
# Autogenerated wrapper script for TOPCOM_jll for armv6l-linux-gnueabihf-cxx03
export chiro2allfinetriangs, chiro2alltriangs, chiro2circuits, chiro2cocircuits, chiro2dual, chiro2finetriang, chiro2finetriangs, chiro2mintriang, chiro2nallfinetriangs, chiro2nalltriangs, chiro2nfinetriangs, chiro2ntriangs, chiro2placingtriang, chiro2triangs, cocircuits2facets, libCHECKREG, libTOPCOM, points2allfinetriangs, points2alltriangs, points2chiro, points2facets, points2finetriang, points2finetriangs, points2flips, points2nallfinetriangs, points2nalltriangs, points2nfinetriangs, points2nflips, points2ntriangs, points2placingtriang, points2triangs, points2volume, santos_22_triang, santos_dim4_triang, santos_triang

using GMP_jll
using cddlib_jll
JLLWrappers.@generate_wrapper_header("TOPCOM")
JLLWrappers.@declare_executable_product(chiro2allfinetriangs)
JLLWrappers.@declare_executable_product(chiro2alltriangs)
JLLWrappers.@declare_executable_product(chiro2circuits)
JLLWrappers.@declare_executable_product(chiro2cocircuits)
JLLWrappers.@declare_executable_product(chiro2dual)
JLLWrappers.@declare_executable_product(chiro2finetriang)
JLLWrappers.@declare_executable_product(chiro2finetriangs)
JLLWrappers.@declare_executable_product(chiro2mintriang)
JLLWrappers.@declare_executable_product(chiro2nallfinetriangs)
JLLWrappers.@declare_executable_product(chiro2nalltriangs)
JLLWrappers.@declare_executable_product(chiro2nfinetriangs)
JLLWrappers.@declare_executable_product(chiro2ntriangs)
JLLWrappers.@declare_executable_product(chiro2placingtriang)
JLLWrappers.@declare_executable_product(chiro2triangs)
JLLWrappers.@declare_executable_product(cocircuits2facets)
JLLWrappers.@declare_library_product(libCHECKREG, "libCHECKREG.so.0")
JLLWrappers.@declare_library_product(libTOPCOM, "libTOPCOM.so.0")
JLLWrappers.@declare_executable_product(points2allfinetriangs)
JLLWrappers.@declare_executable_product(points2alltriangs)
JLLWrappers.@declare_executable_product(points2chiro)
JLLWrappers.@declare_executable_product(points2facets)
JLLWrappers.@declare_executable_product(points2finetriang)
JLLWrappers.@declare_executable_product(points2finetriangs)
JLLWrappers.@declare_executable_product(points2flips)
JLLWrappers.@declare_executable_product(points2nallfinetriangs)
JLLWrappers.@declare_executable_product(points2nalltriangs)
JLLWrappers.@declare_executable_product(points2nfinetriangs)
JLLWrappers.@declare_executable_product(points2nflips)
JLLWrappers.@declare_executable_product(points2ntriangs)
JLLWrappers.@declare_executable_product(points2placingtriang)
JLLWrappers.@declare_executable_product(points2triangs)
JLLWrappers.@declare_executable_product(points2volume)
JLLWrappers.@declare_executable_product(santos_22_triang)
JLLWrappers.@declare_executable_product(santos_dim4_triang)
JLLWrappers.@declare_executable_product(santos_triang)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, cddlib_jll)
    JLLWrappers.@init_executable_product(
        chiro2allfinetriangs,
        "bin/chiro2allfinetriangs",
    )

    JLLWrappers.@init_executable_product(
        chiro2alltriangs,
        "bin/chiro2alltriangs",
    )

    JLLWrappers.@init_executable_product(
        chiro2circuits,
        "bin/chiro2circuits",
    )

    JLLWrappers.@init_executable_product(
        chiro2cocircuits,
        "bin/chiro2cocircuits",
    )

    JLLWrappers.@init_executable_product(
        chiro2dual,
        "bin/chiro2dual",
    )

    JLLWrappers.@init_executable_product(
        chiro2finetriang,
        "bin/chiro2finetriang",
    )

    JLLWrappers.@init_executable_product(
        chiro2finetriangs,
        "bin/chiro2finetriangs",
    )

    JLLWrappers.@init_executable_product(
        chiro2mintriang,
        "bin/chiro2mintriang",
    )

    JLLWrappers.@init_executable_product(
        chiro2nallfinetriangs,
        "bin/chiro2nallfinetriangs",
    )

    JLLWrappers.@init_executable_product(
        chiro2nalltriangs,
        "bin/chiro2nalltriangs",
    )

    JLLWrappers.@init_executable_product(
        chiro2nfinetriangs,
        "bin/chiro2nfinetriangs",
    )

    JLLWrappers.@init_executable_product(
        chiro2ntriangs,
        "bin/chiro2ntriangs",
    )

    JLLWrappers.@init_executable_product(
        chiro2placingtriang,
        "bin/chiro2placingtriang",
    )

    JLLWrappers.@init_executable_product(
        chiro2triangs,
        "bin/chiro2triangs",
    )

    JLLWrappers.@init_executable_product(
        cocircuits2facets,
        "bin/cocircuits2facets",
    )

    JLLWrappers.@init_library_product(
        libCHECKREG,
        "lib/libCHECKREG.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libTOPCOM,
        "lib/libTOPCOM.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        points2allfinetriangs,
        "bin/points2allfinetriangs",
    )

    JLLWrappers.@init_executable_product(
        points2alltriangs,
        "bin/points2alltriangs",
    )

    JLLWrappers.@init_executable_product(
        points2chiro,
        "bin/points2chiro",
    )

    JLLWrappers.@init_executable_product(
        points2facets,
        "bin/points2facets",
    )

    JLLWrappers.@init_executable_product(
        points2finetriang,
        "bin/points2finetriang",
    )

    JLLWrappers.@init_executable_product(
        points2finetriangs,
        "bin/points2finetriangs",
    )

    JLLWrappers.@init_executable_product(
        points2flips,
        "bin/points2flips",
    )

    JLLWrappers.@init_executable_product(
        points2nallfinetriangs,
        "bin/points2nallfinetriangs",
    )

    JLLWrappers.@init_executable_product(
        points2nalltriangs,
        "bin/points2nalltriangs",
    )

    JLLWrappers.@init_executable_product(
        points2nfinetriangs,
        "bin/points2nfinetriangs",
    )

    JLLWrappers.@init_executable_product(
        points2nflips,
        "bin/points2nflips",
    )

    JLLWrappers.@init_executable_product(
        points2ntriangs,
        "bin/points2ntriangs",
    )

    JLLWrappers.@init_executable_product(
        points2placingtriang,
        "bin/points2placingtriang",
    )

    JLLWrappers.@init_executable_product(
        points2triangs,
        "bin/points2triangs",
    )

    JLLWrappers.@init_executable_product(
        points2volume,
        "bin/points2volume",
    )

    JLLWrappers.@init_executable_product(
        santos_22_triang,
        "bin/santos_22_triang",
    )

    JLLWrappers.@init_executable_product(
        santos_dim4_triang,
        "bin/santos_dim4_triang",
    )

    JLLWrappers.@init_executable_product(
        santos_triang,
        "bin/santos_triang",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    find_longest_path(graph_result::BioGraph.GraphResult, optimizer_factory; is_weighted::Bool, source_node::Int64, sink_node::Int64, has_path::String)

Find longest path in graph. Input:

- `GraphResult`
- JuMP optimizer_factory such as `CPLEX.Optimizer`
- `has_path`: optional string that indicates the path must have in longest path.
- `is_weighted`: if `true` find shortest path which is weighted
- `source_node`, `sink_node`: find longest path which has source and sink nodes.
"""
function find_longest_path(graph_result::GraphResult, optimizer_factory; is_weighted::Bool=true, source_node::Int64=0, sink_node::Int64=0, has_path::String="")

    #has_cycle::Bool=false,

    g = copy(graph_result.g)
    e_o_dict, e_i_dict = _edge_label_array_to_dict(graph_result.edge_label)
    l_dict, l_s_dict, l_i_dict = _label_map_array_to_dict(graph_result.node_label)

    start_nodes = []
    end_nodes = []

    add_vertices!(g, 1)

    source_nodes::Array{Int64} = []
    for node in graph_result.source_node
        push!(source_nodes, node.node)
    end

    sink_nodes::Array{Int64} = []
    for node in graph_result.sink_node
        push!(sink_nodes, node.node)
    end

    println("Making graph")
    if sink_node > 0
        add_edge!(g, sink_node, nv(g))
    else
        for i in 1:nv(g)
            if i in sink_nodes
                add_edge!(g, i, nv(g))
            end
        end
    end

    array_e = []
    dict_e = Dict()
    dict_e_contraint = Dict()
    iter = 1
    for e in edges(g)
        push!(array_e, e)
        e_name = string(e.src) * "::" * string(e.dst)
        dict_e[e_name] = iter
        dict_e_contraint[iter] = 0
        iter += 1
    end

    w_dict = _weight_array_to_dict(graph_result.weight)

    weight_e = []
    for e in array_e
        if is_weighted
            if graph_result.weight_inherited
                push!(weight_e, w_dict[e.src])
            else
                if e.dst == nv(g)
                    push!(weight_e, w_dict[e.src])
                else
                    e_name = l_dict[e.src] * "::" * l_dict[e.dst]
                    overlap = e_o_dict[e_name]
                    overlap_int = parse(Int64, chop(overlap))
                    push!(weight_e, w_dict[e.src]-overlap_int)
                end
            end
        else
            push!(weight_e, 1)
        end
    end

    println("Add constraint")
    first_node_constraint = []
    if source_node > 0
        for i in 1:length(array_e)
            if array_e[i].src == source_node
                push!(first_node_constraint, i)
            end
        end
    else
        for i in 1:length(array_e)
            if array_e[i].src in source_nodes
                push!(first_node_constraint, i)
            end
        end
    end

    end_node_constraint = []
    for i in 1:length(array_e)
        if array_e[i].dst == nv(g)
            push!(end_node_constraint, i)
        end
    end

    in_dict = Dict()
    out_dict = Dict()

    for i in 1:nv(g) - 1
        in_dict[i] = []
        out_dict[i] = []
    end

    for j in 1:length(array_e)
        if !(array_e[j].dst in source_nodes) && array_e[j].dst < nv(g)
            push!(in_dict[array_e[j].dst], j)
        end
        if !(array_e[j].src in source_nodes) && array_e[j].src < nv(g)
            push!(out_dict[array_e[j].src], j)
        end
    end

    cycle_constraint = []
    cycle_length_constraint = []

    path_constraint = []

    if has_path != ""
        dict_g = Dict()
        for i in 1:nv(g)
            dict_g[i] = 0
        end
        for p_dummy in graph_result.path
            if has_path == p_dummy.name
                p_arr = p_dummy.path
                for i in 1:(length(p_arr)-1)
                    e_name = string(p_arr[i]) * "::" * string(p_arr[i+1])
                    push!(path_constraint, dict_e[e_name])
                    dict_e_contraint[dict_e[e_name]] += 1
                    dict_g[p_arr[i+1]] += 1
                end
            end
        end
    end

    g_opt = SimpleDiGraph(nv(g))

    obj = nv(g)

    println("Solving")
    i = 1
    if has_path == ""
        while true
            model = Model();
            set_optimizer(model, optimizer_factory);
            set_silent(model)
            @variable(model, x[1:ne(g)], binary = true);
            @objective(model, Max, sum(x[i] * weight_e[i] for i = 1:ne(g)));
            @constraint(model, sum(x[i] for i in first_node_constraint) == 1);
            @constraint(model, sum(x[i] for i in end_node_constraint) == 1);
            for i in 1:nv(g) - 1
                @constraint(model, sum(x[j] for j in in_dict[i]) <= 1);
                @constraint(model, sum(x[j] for j in in_dict[i]) - sum(x[j] for j in out_dict[i]) == 0);
            end
            for i in 1:length(cycle_constraint)
                @constraint(model, sum(x[j] for j in cycle_constraint[i]) <= cycle_length_constraint[i] - 1);
            end
            optimize!(model);
            g_opt = SimpleDiGraph(nv(g))
            for i in 1:ne(g)
                if value(x[i]) >= 0.5
                    add_edge!(g_opt, array_e[i].src, array_e[i].dst)
                end
            end
            obj = objective_value(model)
            println("inter: ", i, ", objective: ", obj, ", num new cycle: ", length(simplecycles(g_opt)))
            i += 1
            for cycle in simplecycles_iter(g_opt, 100)
                cycle_dump = []
                push!(cycle_length_constraint, length(cycle))
                for j in 1:length(array_e)
                    if array_e[j].src in cycle && array_e[j].dst in cycle
                        push!(cycle_dump, j)
                    end
                end
                push!(cycle_constraint, cycle_dump)
            end
            length(simplecycles(g_opt)) > 0 || break
        end
    else
        while true
            model = Model();
            set_optimizer(model, optimizer_factory);
            set_silent(model)
            @variable(model, x[1:ne(g)], integer = true, lower_bound=0);
            @objective(model, Max, sum(x[i] * weight_e[i] for i = 1:ne(g)));
            @constraint(model, sum(x[i] for i in first_node_constraint) == 1);
            @constraint(model, sum(x[i] for i in end_node_constraint) == 1);
            for i in 1:ne(g)
                if dict_e_contraint[i] == 0
                    @constraint(model, x[i]  <= 1);
                else
                    @constraint(model, x[i]  == dict_e_contraint[i]);
                end
            end
            for i in 1:nv(g) - 1
                if dict_g[i] == 0
                    @constraint(model, sum(x[j] for j in in_dict[i]) <= 1);
                else
                    @constraint(model, sum(x[j] for j in in_dict[i]) == dict_g[i]);
                end
                @constraint(model, sum(x[j] for j in in_dict[i]) - sum(x[j] for j in out_dict[i]) == 0);
            end
            for i in 1:length(cycle_constraint)
                @constraint(model, sum(x[j] for j in cycle_constraint[i]) <= cycle_length_constraint[i] - 1);
            end
            optimize!(model);
            g_opt = SimpleDiGraph(nv(g))
            for i in 1:ne(g)
                if value(x[i]) >= 0.8
                    add_edge!(g_opt, array_e[i].src, array_e[i].dst)
                end
            end
            obj = objective_value(model)
            g_opt_com = weakly_connected_components(g_opt)
            lone_nodes = 0
            for com in g_opt_com
                if length(com) > 1
                    sg_source = collect(intersect(Set(source_nodes), Set(com)))
                    if length(sg_source) == 0
                        cycle_dump = []
                        push!(cycle_length_constraint, length(com));
                        for j in 1:length(array_e)
                            if array_e[j].src in com && array_e[j].dst in com
                                push!(cycle_dump, j)
                            end
                        end
                        if cycle_dump != []
                            push!(cycle_constraint, cycle_dump)
                        end
                    end
                else
                    lone_nodes+=1
                end
            end
            println("objective: ", obj, ", num of lone cycle: ", length(g_opt_com)-lone_nodes - 1)
            length(g_opt_com)-lone_nodes > 1 || break
        end
    end

    path = _g_opt_to_path(g_opt, source_nodes)
    label_path = _label_path(graph_result, path)
    return LongestPath(graph_result, g_opt, path, label_path, obj)
end

"""
    get_gfa(g_result::BioGraph.LongestPath; outfile::String)

Get GFA output of `LongestPath`.
"""
function get_gfa(longest_path::LongestPath, outfile::String)
    g_result = longest_path.graph
    lp_result = longest_path.label_path
    e_o_dict, e_i_dict = _edge_label_array_to_dict(g_result.edge_label)
    open(outfile, "w") do io
        write(io, "H\tVN:Z:1.0" * "\n")
        for node in g_result.node_label
            s_dummy = split(node.label, "\t")
            write(io, "S\t" * s_dummy[1] * "\t" * node.sequence * "\t" * node.info * "\n")
        end
        for edge in g_result.edge_label
            write(io, "L\t" * edge.src * "\t" * edge.dst * "\t" * edge.overlap * "\t" * edge.info * "\n")
        end
        PathName = "LongestPath"
        SegmentNames = ""
        for i in 1:length(lp_result)
            lp_dummy = split(lp_result[i], "\t")
            if i < length(lp_result)
                SegmentNames = SegmentNames * lp_dummy[1] * lp_dummy[2] * ","
            else
                SegmentNames = SegmentNames * lp_dummy[1] * lp_dummy[2]
            end
        end
        Overlaps = ""
        for i in 1:length(lp_result) - 1
            e_name = lp_result[i] * "::" * lp_result[i + 1]
            if i < length(lp_result) - 1
                Overlaps = Overlaps * e_o_dict[e_name] * ","
            else
                Overlaps = Overlaps * e_o_dict[e_name]
            end
        end
        write(io, "P\t" * PathName * "\t" * SegmentNames * "\t" * Overlaps)
    end
end

"""
    get_fasta(g_result::BioGraph.LongestPath; outfile::String)

Get FastA output of `LongestPath`.
"""
function get_fasta(longest_path::LongestPath, outfile::String="")
    g_result = longest_path.graph
    lp_result = longest_path.label_path
    l_dict, l_s_dict, l_i_dict = _label_map_array_to_dict(g_result.node_label)
    e_o_dict, e_i_dict = _edge_label_array_to_dict(g_result.edge_label)
    lp_dict = Dict()

    for i in 1:length(l_dict)
        lp_dummy = split(l_dict[i], "\t")
        if lp_dummy[2] == "-"
            lp_dict[l_dict[i]] = string(reverse_complement(LongDNASeq(l_s_dict[i])))
        else
            lp_dict[l_dict[i]] = l_s_dict[i]
        end
    end
    fasta_output = ""
    for i in 1:length(lp_result)
        if i == length(lp_result)
            fasta_output = fasta_output*lp_dict[lp_result[i]]
        else
            e_name = lp_result[i] * "::" * lp_result[i+1]
            overlap = e_o_dict[e_name]
            overlap_int = parse(Int64, chop(overlap))
            fasta_output = fasta_output*chop(lp_dict[lp_result[i]], tail = overlap_int)
        end
    end
    if outfile != ""
        open(outfile, "w") do io
            write(io, ">linear_path\n")
            write(io, fasta_output)
        end
    end
    return fasta_output
end


text:
include("Graph.jl")
include("Exp2.jl")

using Printf

function runMain(ags1, ags2; resultDir = "../result/")
    lg = open(resultDir*"log.txt", "a")
    G = readUnLabeledGraph(ags1)
    println(lg, ags1, " ", G.n, " ", G.m)
    ds = Dict{String, String}()
    ds["1"] = "RoleSim : "
    ds["2"] = "StructSim : "
    ds["3"] = "ForestSimExact : "
    ds["4"] = "ForestSimApprox : "
    df = Dict{String, Function}()
    df["1"] = testRoleSim
    df["2"] = testStructSim
    df["3"] = testForestSim
    df["4"] = testForestSimApprox
    for x in split(ags2, ',')
        println(lg, ds[x], df[x](G, 10))
    end
    close(lg)
end

runMain("Diseasome", "1,2,3,4")
runMain("EmailUniv", "1,2,3,4")
runMain("Hamster", "1,2,3,4")
runMain("GridWorm", "1,2,3,4")
runMain("GrQc", "1,2,3,4")
runMain("Erdos992", "1,2,3,4")
runMain("Reality", "1,2,3,4")
runMain("Dmela", "1,2,3,4")
runMain("HepPh", "2,3,4")
runMain("PagesCompany", "2,3,4")
runMain("AstroPh", "2,3,4")
runMain("Gplus", "2,3,4")
runMain("Brightkite", "2,4")
runMain("BlogCatalog", "2,4")
runMain("Douban", "2,4")
runMain("MathSciNet", "4")
runMain("Flickr", "4")
runMain("IMDB", "4")
runMain("YoutubeSnap", "4")
runMain("Flixster", "4")


text:
#push!(LOAD_PATH,"../src/")

using Documenter, DataDrivenDiffEq
using Flux, SymbolicRegression
using Literate

include("pages.jl")

ENV["GKSwstype"] = "100"

# Evaluate the example directory

src = joinpath(@__DIR__, "src")
lit = joinpath(@__DIR__, "examples")
excludes = []#["symbolic_regression.jl"]
tutorials = []

for (root, _, files) ∈ walkdir(lit), file ∈ files
  file ∈ excludes && continue
  fname, fext = splitext(file)

  fext == ".jl" || continue
  ipath = joinpath(root, file)
  opath = joinpath(splitdir(replace(ipath, lit=>src))[1], "examples")
  script = Literate.script(ipath, opath, execute = false, comments = false)
  code = strip(read(script, String))
  mdpost(str) = replace(str, "@__CODE__" => code)
  Literate.markdown(ipath, opath)
  Literate.markdown(ipath, opath, execute = false, postprocess = mdpost)
  if fname == "0_getting_started"
    pushfirst!(tutorials,  relpath(joinpath(opath, fname*".md"), src))
  else
    push!(tutorials, relpath(joinpath(opath, fname*".md"), src))
  end
end

# Create the docs
makedocs(
    sitename="DataDrivenDiffEq.jl",
    authors="Julius Martensen, Christopher Rackauckas",
    modules=[DataDrivenDiffEq],
    clean=true,doctest=false,
    format = Documenter.HTML(analytics = "UA-90474609-3",
                             assets = ["assets/favicon.ico"],
                             canonical="https://datadriven.sciml.ai/stable/"),
    pages=pages
)

deploydocs(
   repo = "github.com/SciML/DataDrivenDiffEq.jl.git";
   push_preview = true
)


text:
function addStrings(s1::String, s2::String) :: Int64
    chars = ['0','1','2','3','4','5','6','7','8','9']
    dd = Dict([(chars[n+1],n) for n in 0:9])

    result = 0
    len1 = length(s1)
    for ind in 1:len1
        result += dd[s1[ind]] * (10 ^ len1)/(10 ^ ind)
    end

    len2 = length(s2)
    for ind in 1:len2
        result += dd[s2[ind]] * (10 ^ len2)/(10 ^ ind)
    end

    result
end

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
module SemidefiniteOptInterface

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

using Compat
using Compat.LinearAlgebra # for diag

abstract type AbstractSDOptimizer <: MOI.AbstractOptimizer end

include("interface.jl")

const SVF = MOI.SingleVariable
const VVF = MOI.VectorOfVariables
const VF  = Union{SVF, VVF}
const SAF{T} = MOI.ScalarAffineFunction{T}
const ASF{T} = Union{SVF, SAF{T}}

const ZS = Union{MOI.EqualTo, MOI.Zeros}
const NS = Union{MOI.GreaterThan, MOI.Nonnegatives}
const PS = Union{MOI.LessThan, MOI.Nonpositives}
const DS = MOI.PositiveSemidefiniteConeTriangle
const SupportedSets = Union{ZS, NS, PS, DS}

const VI = MOI.VariableIndex
const CI{F, S} = MOI.ConstraintIndex{F, S}

mutable struct SOItoMOIBridge{T, SIT <: AbstractSDOptimizer} <: MOI.AbstractOptimizer
    sdoptimizer::SIT
    setconstant::Dict{Int64, T}
    blkconstant::Dict{Int, T}
    objconstant::T
    objsign::Int
    objshift::T
    nconstrs::Int
    nblocks::Int
    blockdims::Vector{Int}
    free::BitSet
    varmap::Vector{Vector{Tuple{Int, Int, Int, T, T}}} # Variable Index vi -> blk, i, j, coef, shift # x = sum coef * block(X, blk)[i, j] + shift
    zeroblock::Dict{CI, Int}
    constrmap::Dict{CI, UnitRange{Int}} # Constraint Index ci -> cs
    double::Vector{CI} # created when there are two cones for same variable
    function SOItoMOIBridge{T}(sdoptimizer::SIT) where {T, SIT}
        new{T, SIT}(sdoptimizer, Dict{Int64, T}(), Dict{Int, T}(),
            zero(T), 1, zero(T), 0, 0,
            Int[],
            BitSet(),
            Vector{Tuple{Int, Int, Int, T}}[],
            Dict{CI, Int}(),
            Dict{CI, UnitRange{Int}}(),
            CI[])
    end
end
varmap(optimizer::SOItoMOIBridge, vi::VI) = optimizer.varmap[vi.value]
function setvarmap!(optimizer::SOItoMOIBridge{T}, vi::VI, v::Tuple{Int, Int, Int, T, T}) where T
    setvarmap!(optimizer, vi, [v])
end
function setvarmap!(optimizer::SOItoMOIBridge{T}, vi::VI, vs::Vector{Tuple{Int, Int, Int, T, T}}) where T
    optimizer.varmap[vi.value] = vs
end

SDOIOptimizer(sdoptimizer::AbstractSDOptimizer, T=Float64) = SOItoMOIBridge{T}(sdoptimizer)

include("load.jl")

function MOI.get(optimizer::SOItoMOIBridge, attr::MOI.SolverName)
    return MOI.get(optimizer.sdoptimizer, attr)
end

function MOI.is_empty(optimizer::SOItoMOIBridge)
    isempty(optimizer.double) &&
    isempty(optimizer.setconstant) &&
    isempty(optimizer.blkconstant) &&
    iszero(optimizer.objconstant) &&
    optimizer.objsign == 1 &&
    iszero(optimizer.objshift) &&
    iszero(optimizer.nconstrs) &&
    iszero(optimizer.nblocks) &&
    isempty(optimizer.blockdims) &&
    isempty(optimizer.free) &&
    isempty(optimizer.varmap) &&
    isempty(optimizer.zeroblock) &&
    isempty(optimizer.constrmap)
end
function MOI.empty!(optimizer::SOItoMOIBridge{T}) where T
    for s in optimizer.double
        MOI.delete(m, s)
    end
    MOI.empty!(optimizer.sdoptimizer)
    optimizer.double = CI[]
    optimizer.setconstant = Dict{Int64, T}()
    optimizer.blkconstant = Dict{Int, T}()
    optimizer.objconstant = zero(T)
    optimizer.objsign = 1
    optimizer.objshift = zero(T)
    optimizer.nconstrs = 0
    optimizer.nblocks = 0
    optimizer.blockdims = Int[]
    optimizer.free = BitSet()
    optimizer.varmap = Vector{Tuple{Int, Int, Int, T}}[]
    optimizer.zeroblock = Dict{CI, Int}()
    optimizer.constrmap = Dict{CI, UnitRange{Int}}()
end

function setconstant!(optimizer::SOItoMOIBridge, ci::CI, s) end
function setconstant!(optimizer::SOItoMOIBridge, ci::CI, s::MOI.AbstractScalarSet)
    optimizer.setconstant[ci.value] = MOIU.getconstant(s)
end
function set_constant(optimizer::SOItoMOIBridge,
                      ci::CI{<:MOI.AbstractScalarFunction,
                             <:MOI.AbstractScalarSet})
    return optimizer.setconstant[ci.value]
end
function set_constant(optimizer::SOItoMOIBridge{T}, ci::CI) where T
    return zeros(T, length(optimizer.constrmap[ci]))
end
function addblkconstant(optimizer::SOItoMOIBridge, ci::CI{<:Any, <:Union{NS, PS}}, x)
    blk = -ci.value
    return x .+ optimizer.blkconstant[blk]
end
addblkconstant(optimizer::SOItoMOIBridge, ci::CI, x) = x

function MOI.supports(
    optimizer::SOItoMOIBridge{T},
    ::Union{MOI.ObjectiveSense,
            MOI.ObjectiveFunction{<:Union{MOI.SingleVariable,
                                          MOI.ScalarAffineFunction{T}}}}) where T
    return true
end

# Zeros and Nonpositives supports could be removed thanks to variable bridges
# * `VectorOfVariables`-in-`Zeros` would return a `VectorAffineFunction` with
#   zero constant and no variable created.
# * `VectorOfVariables`-in-`Nonpositives` would create variables in
#   `Nonnegatives` and return a `VectorAffineFunction` containing `-` the
#    variables.
function MOI.supports_constraint(
    ::SOItoMOIBridge, ::Type{MOI.VectorOfVariables},
    ::Type{<:Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
                   MOI.PositiveSemidefiniteConeTriangle}})
    return true
end
# This support could be remove thanks to variable bridges.
# The VectorizeVariableBridge would redirect to the above case and then the
# resulting function would be shifted by the constant.
function MOI.supports_constraint(
    ::SOItoMOIBridge{T}, ::Type{MOI.SingleVariable},
    ::Type{<:Union{MOI.EqualTo{T}, MOI.GreaterThan{T}, MOI.LessThan{T}}}) where T
    return true
end
function MOI.supports_constraint(
    ::SOItoMOIBridge{T}, ::Type{MOI.ScalarAffineFunction{T}},
    ::Type{MOI.EqualTo{T}}) where T
    return true
end

function MOI.copy_to(dest::SOItoMOIBridge, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end
MOIU.supports_allocate_load(::SOItoMOIBridge, copy_names::Bool) = !copy_names

MOI.optimize!(m::SOItoMOIBridge) = MOI.optimize!(m.sdoptimizer)

# Objective

function MOI.get(m::SOItoMOIBridge, ::MOI.ObjectiveValue)
    m.objshift + m.objsign * getprimalobjectivevalue(m.sdoptimizer) + m.objconstant
end

# Attributes

const SolverStatus = Union{MOI.TerminationStatus, MOI.PrimalStatus, MOI.DualStatus}
MOI.get(m::SOItoMOIBridge, s::SolverStatus) = MOI.get(m.sdoptimizer, s)

MOI.get(m::SOItoMOIBridge, ::MOI.ResultCount) = 1

function _getblock(M, blk::Integer, s::Type{<:Union{NS, ZS}})
    return diag(block(M, blk))
end
function _getblock(M, blk::Integer, s::Type{<:PS})
    return -diag(block(M, blk))
end
# Vectorized length for matrix dimension d
sympackedlen(d::Integer) = (d*(d+1)) >> 1
function _getblock(M::AbstractMatrix{T}, blk::Integer, s::Type{<:DS}) where T
    B = block(M, blk)
    d = Compat.LinearAlgebra.checksquare(B)
    n = sympackedlen(d)
    v = Vector{T}(undef, n)
    k = 0
    for j in 1:d
        for i in 1:j
            k += 1
            v[k] = B[i, j]
        end
    end
    @assert k == n
    return v
end
function getblock(M, blk::Integer, s::Type{<:MOI.AbstractScalarSet})
    vd = _getblock(M, blk, s)
    @assert length(vd) == 1
    return vd[1]
end
function getblock(M, blk::Integer, s::Type{<:MOI.AbstractVectorSet})
    return _getblock(M, blk, s)
end

getvarprimal(m::SOItoMOIBridge, blk::Integer, S) = getblock(getX(m.sdoptimizer), blk, S)
function getvardual(m::SOItoMOIBridge, blk::Integer, S)
    z = getZ(m.sdoptimizer)
    b = getblock(z, blk, S)
    return getblock(getZ(m.sdoptimizer), blk, S)
end

function MOI.get(m::SOItoMOIBridge{T}, ::MOI.VariablePrimal, vi::VI) where T
    X = getX(m.sdoptimizer)
    x = zero(T)
    for (blk, i, j, coef, shift) in varmap(m, vi)
        x += shift
        if blk != 0
            x += block(X, blk)[i, j] * sign(coef)
        end
    end
    return x
end
function MOI.get(m::SOItoMOIBridge, vp::MOI.VariablePrimal, vi::Vector{VI})
    return MOI.get.(m, vp, vi)
end

function _getattribute(m::SOItoMOIBridge, ci::CI{<:ASF}, f)
    cs = m.constrmap[ci]
    @assert length(cs) == 1
    return f(m, first(cs))
end
function _getattribute(m::SOItoMOIBridge, ci::CI{<:VVF}, f)
    return f.(m, m.constrmap[ci])
end

function MOI.get(m::SOItoMOIBridge, a::MOI.ConstraintPrimal,
                 ci::CI{F, S}) where {F, S}
    if ci.value >= 0
        return set_constant(m, ci)
    else
        # Variable Function-in-S with S different from Zeros and EqualTo and not a double variable constraint
        blk = -ci.value
        return addblkconstant(m, ci, getvarprimal(m, blk, S))
    end
end

function MOI.get(m::SOItoMOIBridge, ::MOI.ConstraintDual, ci::CI{<:VF, S}) where S<:SupportedSets
    if ci.value < 0
        return getvardual(m, -ci.value, S)
    else
        dual = _getattribute(m, ci, getdual)
        if haskey(m.zeroblock, ci) # ZS
            return dual + getvardual(m, m.zeroblock[ci], S)
        else # var constraint on unfree constraint
            return dual
        end
    end
end

function getdual(m::SOItoMOIBridge{T}, c::Integer) where T
    if c == 0
        return zero(T)
    else
        return -gety(m.sdoptimizer)[c]
    end
end
function MOI.get(m::SOItoMOIBridge, ::MOI.ConstraintDual, ci::CI)
    return _getattribute(m, ci, getdual)
end

include("sdpa.jl")
include("mock.jl")

end # module


text:
"""
Snape's Riddle - Author JK Rowling
Danger lies before you, while safety lies behind,
Two of us will help you, whichever you would find,
One among us seven will let you move ahead,
Another will transport the drinker back instead,
Two among our number hold only nettle wine,
Three of us are killers, waiting hidden in line.
Choose, unless you wish to stay here for evermore,
To help you in your choice, we give you these clues four:
First, however slyly the poison tries to hide
You will always find some on nettle wine’s left side;
Second, different are those who stand at either end,
But if you would move onwards, neither is your friend;
Third, as you see clearly, all are different size,
Neither dwarf nor giant holds death in their insides;
Fourth, the second left and the second on the right
Are twins once you taste them, though different at first sight.
"""

# Codes:
# 1 is netted wine
# 2 move ahead
# 3 move back
# 4 kill

# Seven Potions a,b,c,d,e,f,g each have 1 of 4 possible values
# a is left of b which is left of c etc.
Ω,℧ = ABparse(["a, b, c, d, e, f, g  ∈  [1,2,3,4]"])

# Two among our number hold only nettle wine, #code 1
Ω,℧ = ABparse(["{{i}} == 1 {{2}}"], Ω=Ω,℧=℧); Ω[℧]
# Equivalent to ABparse(["a|b|c|d|e|f|g == 1 {2}"], Ω=Ω,℧=℧)

# One among us seven will let you move ahead, #code 2
Ω,℧ = ABparse(["{{i}} == 2 {{1}}"], Ω=Ω,℧=℧)
#Ω,℧ = ABparse(["a|b|c|d|e|f|g == 2 {1}"], Ω=Ω,℧=℧)

#Another will transport the drinker back instead, #code 3
Ω,℧ = ABparse(["{{i}} == 3 {{1}}"], Ω=Ω,℧=℧)

# Three of us are killers #code 4
Ω,℧ = ABparse(["{{i}} == 4 {{3}}"], Ω=Ω,℧=℧)
# This last specification does not reduce the feasible set

Ω[℧] # 420 possible outcomes

# First, however slyly the poison tries to hide
# You will always find some on nettle wine’s left side;
# So if b is poison then a is nettle wine
Ω,℧ = ABparse(["{{j+1}} == 1 ==> {{j}} == 4"], Ω=Ω,℧=℧); Ω[℧]

# Equivalent to:
# ABparse(["b == 1 ==> a == 4","c == 1 ==> b == 4","d == 1 ==> c == 4",
#          "e == 1 ==> d == 4","f == 1 ==> e == 4","g == 1 ==> f == 4"], Ω=Ω,℧=℧)
#
# Wine cannot be in the first bottle since poison is to the left
Ω,℧ = ABparse(["a != 1"], Ω=Ω,℧=℧)

Ω[℧] # 60 possible outcomes

# Second, different are those who stand at either end,
Ω,℧ = ABparse(["a != g"], Ω=Ω,℧=℧)

# But if you would move onwards, neither is your friend;
Ω,℧ = ABparse(["a|g == 2 {0}"], Ω=Ω,℧=℧)

Ω[℧] # 30 possible outcomes

# But if you would move onwards, neither is your friend;
# Third, as you see clearly, all are different size,
# Neither dwarf nor giant holds death in their insides;
# c dwarf, f are giant
Ω,℧ = ABparse(["c|f == 4 {0}"], Ω=Ω,℧=℧)

Ω[℧] # 5 possible outcomes

# Fourth, the second left and the second on the right
# Are twins once you taste them, though different at first sight.
Ω,℧ = ABparse(["b == f"], Ω=Ω, ℧=℧)

Ω[℧] # 1 possible outcome
# You pick the 3rd (move ahead) and the 7th (move back)


text:
ctr_name(::MOI.EqualTo) = "FX"
ctr_name(::MOI.LessThan) = "UB"
ctr_name(::MOI.GreaterThan) = "LB"
ctr_name(::MOI.Integer) = "IN"
ctr_name(::MOI.ZeroOne) = "BN"

function fill_constraint_names(model)
    for (F,S) in MOI.get(model, MOI.ListOfConstraints())
        if F <: MOI.ScalarAffineFunction#, MOI.ScalarQuadraticFunction]
            # check
        elseif F <: MOI.VectorOfVariables
            # error
        elseif F <: MOI.SingleVariable
            for index in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
                name = MOI.get(model, MOI.ConstraintName(), index)
                if name == ""
                    f = MOI.get(model, MOI.ConstraintFunction(), index)
                    name = MOI.get(model, MOI.VariableName(), f.variable)
                    MOI.set(model, MOI.ConstraintName(), index, name*"_"*ctr_name(S))
                end
            end
        else
            #error
        end
    end
end

function remove_intervals(model)
    for (F,S) in MOI.get(model, MOI.ListOfConstraints())
        if S <: MOI.Interval
            for index in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
                cname = MOI.get(model, MOI.ConstraintName(), index)
                func = MOI.get(model, MOI.ConstraintFunction(), index)
                set = MOI.get(model, MOI.ConstraintSet(), index)
                c0 = MOI.transform(model, index, MOI.LessThan(set.upper))
                MOI.set(model, MOI.ConstraintName(), c0, cname) # choose a nice name
                c = MOI.add_constraint(model, func, MOI.GreaterThan(set.lower))
                MOI.set(model, MOI.ConstraintName(), c, cname * "_LB") # choose a nice name
            end
        end
    end
end


text:
using TracedSample, StatsBase, Random, Test

@testset "Addressed calls with global RNG" begin

for fn in (sample, wsample)
   Random.seed!(0)
   traced = fn(:addr, [1, 2, 3], weights([1, 2, 3]))
   Random.seed!(0)
   untraced = fn([1, 2, 3], weights([1, 2, 3]))
   @test traced == untraced

   Random.seed!(0)
   traced = fn(:addr, [1, 2, 3], weights([1, 2, 3]), 3; replace=false)
   Random.seed!(0)
   untraced = fn([1, 2, 3], weights([1, 2, 3]), 3; replace=false)
   @test traced == untraced
end

for fn in (sample!, wsample!)
   Random.seed!(0)
   traced = fn(:addr, [1, 2, 3], weights([1, 2, 3]), zeros(5))
   Random.seed!(0)
   untraced = fn([1, 2, 3], weights([1, 2, 3]), zeros(5))
   @test traced == untraced
end

end

@testset "Addressed calls with custom RNG" begin

rng = MersenneTwister()

for fn in (sample, wsample)
   Random.seed!(rng, 42)
   traced = fn(rng, :addr, [1, 2, 3], weights([1, 2, 3]))
   Random.seed!(rng, 42)
   untraced = fn(rng, [1, 2, 3], weights([1, 2, 3]))
   @test traced == untraced

   Random.seed!(rng, 42)
   traced = fn(rng, :addr, [1, 2, 3], weights([1, 2, 3]), 3; replace=false)
   Random.seed!(rng, 42)
   untraced = fn(rng, [1, 2, 3], weights([1, 2, 3]), 3; replace=false)
   @test traced == untraced
end

for fn in (sample!, wsample!)
   Random.seed!(rng, 42)
   traced = fn(rng, :addr, [1, 2, 3], weights([1, 2, 3]), zeros(5))
   Random.seed!(rng, 42)
   untraced = fn(rng, [1, 2, 3], weights([1, 2, 3]), zeros(5))
   @test traced == untraced
end

end


text:
using MPIMeasurements
using HDF5

scanner = MPIScanner("HeadScanner.toml")
daq = getDAQ(scanner)

params = toDict(daq.params)

params["studyName"]="TestTobi"
params["studyDescription"]="A very cool measurement"
params["scannerOperator"]="Tobi"
params["dfStrength"]=[5e-3]  #[5e-3]
params["acqNumAverages"]=1000
#params["calibFieldToVolt"]=[1]

params["acqFFSequence"] = "HeadScanner102Triangle"
#params["acqFFValues"] = repeat( cat(1,x,reverse(x[2:end-1])),inner=5)
params["acqNumPeriodsPerFrame"]=length(params["acqFFValues"])

robot = getRobot(scanner)
S,S2 = measurementContReadAndSave(daq, robot, params, controlPhase=true, showFT=true)

h5write("systemMatrixNew.h5","/matrix",S)
h5write("phantom.h5","/matrix",S2)


text:
"""
    create_network(config)

Create a network with the Barabási-Albert model.

# Arguments

- `config`: Config file

# Example
```julia-repl
julia> using Random

julia> Random.seed!(0);

julia> create_network(Config(network = cfg(net(agent_count = 100))))
{100, 206} directed simple Int64 graph

```

See also: [`update_network!`](@ref)
"""
function create_network(
    agent_list::AbstractArray,
    config::Config
)
    # this algorithm is modelled after the python networkx implementation:
    # https://github.com/networkx/networkx/blob/master/networkx/generators/random_graphs.py#L655

    g = SimpleDiGraph(config.network.agent_count)
    pref_attach_list = collect(1:config.network.agent_count)

    for source in 1:config.network.agent_count

        shuffle!(pref_attach_list)
        m0 = agent_list[source].desired_input_count

        if m0 >= config.network.agent_count
            m0 = config.network.agent_count - 1
        elseif m0 <= 0
            break
        end

        targets = Array{Int64}(undef, 0)

        for i in pref_attach_list

            if length(targets) >= m0
                break
            end

            if !(i in targets) && i != source
                push!(targets, i)
            end
        end

        for e in zip(targets, fill(source, m0))
            add_edge!(g, e[1], e[2])
        end

        append!(pref_attach_list, targets)
    end

    return g
end

"""
    update_network!(state, config)

Update a network in a simulation.

# Arguments
- `state`: a tuple of the current graph and agent_list
- `config`: Config object as provided by Config

See also: [`create_network`](@ref), [`Agent`](@ref), [`generate_opinion`](@ref), [`generate_inclin_interact`](@ref), , [`generate_check_regularity`](@ref)
"""
function update_network!(
    state::Tuple{AbstractGraph, AbstractArray},
    config::Config
)
    graph, agent_list = state
    pref_attach_list = [
        src(e) for e in edges(graph) if agent_list[src(e)].active
    ]

    for _ in 1:config.network.growth_rate
        push!(
            agent_list,
            Agent(
                nv(graph) + 1,
                generate_opinion(),
                generate_inclin_interact(),
                generate_check_regularity(),
                generate_desired_input_count(
                    config.agent_props.mean_desired_input_count
                )
            )
        )

        add_vertex!(graph)
        targets = Array{Int64}(undef, 0)

        shuffle!(pref_attach_list)
        for i in pref_attach_list
            if (
                length(targets)
                >= ceil(Int, last(agent_list).desired_input_count / 2)
            )
                break
            end

            if !(i in targets) && i != nv(graph)
                push!(targets, i)
            end
        end

        for t in targets
            add_edge!(graph, t, nv(graph))
        end
    end

    return state
end

# suppress output of include()
;


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "UDUNITS"
version = v"2.2.28"

# Collection of sources required to complete build
sources = [
    # .tar.gz archive is missing cmake files: https://github.com/Unidata/UDUNITS-2/issues/108
    ArchiveSource("https://artifacts.unidata.ucar.edu/repository/downloads-udunits/udunits-$(version).zip",
                  "e09d31db68f9a840a0663c7e9909101957733ff0310761b9906f4722e0d92c44"),
    DirectorySource("bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
apk add --upgrade cmake --repository=http://dl-cdn.alpinelinux.org/alpine/v3.14/main # requires CMake 3.19

cd $WORKSPACE/srcdir/udunits-*

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/freebsd.patch

mkdir build
cd build

cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release
make -j${nprocs}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libudunits2", :libudunits2),
    ExecutableProduct("udunits2", :udunits2)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Expat_jll", uuid="2e619515-83b5-522b-bb60-26c02a35a201"); compat="2.2.10")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "Node pop plots" begin

nodePopPlot( mpSim, 12, "A junior", "B junior", "Reserve junior",
    "A senior", "B senior", "Master", "Career" )
nodeEvolutionPlot( mpSim, 12, "A junior", "B junior", "Reserve junior",
    "A senior", "B senior", "Master", "Career" )
#=
nodePopPlot( mpSim, 12, "A junior", "B junior", "Reserve junior",
    "A senior", "B senior", "Master", "Career", showPlot=false,
    filename="plot/node pop test" )
nodeEvolutionPlot( mpSim, 12, "A junior", "B junior", "Reserve junior",
    "A senior", "B senior", "Master", "Career", showPlot=false,
    filename="plot/node evolution test" )
=#
end  # @testset "Node pop plots"


text:
struct HomogeneousMesh{VecArrayType,VecType<:AbstractVec,CellType<:AbstractCell2Node,CNF,FaceType<:AbstractFace2Node,NBF,NF} <: AbstractMesh
    nodes::VecArrayType
    cells::Vector{CellType} # c2n
    faces::Vector{FaceType} # f2n
    bfaces::Vector{FaceType} # bf2n
    #c2c::Vector{UInt}
    c2f::Vector{NTuple{CNF,Int}}
  #  f2c::Vector{Face2Cell}
  #  bf2c::Vector{BFace2Cell}
    cv::Vector{Float64}
    f2cloops::Face2CellLoop{NBF,NF,VecType}
end

vec_type(::Type{HomogeneousMesh{V,C,CNF,FT,NBF,NF}}) where {V,C,CNF,FT,NBF,NF} = V
cell_type(::Type{HomogeneousMesh{V,C,CNF,FT,NBF,NF}}) where {V,C,CNF,FT,NBF,NF} = C
cell_nnodes(::Type{HomogeneousMesh{V,C,CNF,FT,NBF,NF}}) where {V,C,CNF,FT,NBF,NF} = CNF
face_type(::Type{HomogeneousMesh{V,C,CNF,FT,NBF,NF}}) where {V,C,CNF,FT,NBF,NF} = FT
nbfaces(::Type{HomogeneousMesh{V,C,CNF,FT,NBF,NF}}) where {V,C,CNF,FT,NBF,NF} = NBF
nfaces(::Type{HomogeneousMesh{V,C,CNF,FT,NBF,NF}}) where {V,C,CNF,FT,NBF,NF} = NF

for op in (:vec_type,:cell_type,:cell_nnodes,:face_type,:nbfaces,:nfaces)
    @eval begin
        @inline ($op)(a::HomogeneousMesh) = ($op)(typeof(a))
    end
end

function HomogeneousMesh(inputfile::AbstractString)
    nodes, cells, bc_cell, bfaces, bc_face = readneutral(inputfile)
    cv = (x-> 1/(volume(x,nodes))).(cells)
    faces, c2f, bf2c, bfn, bfc, bcv, bccenter = cell_connectivity(cells,bfaces,bc_face,nodes)
    f2c, fn, fc, cvf, ccenter = face_connectivity(faces,c2f,cells,nodes)
    CellType = eltype(cells)
    CNF = nfaces(CellType)
    f2cloops = Face2CellLoop{length(bfaces),length(faces),eltype(nodes),typeof(fn)}(bf2c,bc_face,bfn,bfc,bcv,bccenter,f2c,fn,fc,cvf,ccenter,get_AoL(f2c,fn,ccenter))
    return HomogeneousMesh{typeof(nodes),eltype(nodes),CellType,CNF,eltype(faces),length(bfaces),length(faces)}(nodes,cells,faces,bfaces,c2f,cv,f2cloops)
end

function get_AoL(f2c,fn,cc)
    AoL = zeros(length(fn))
    @inbounds for i in linearindices(f2c)
        n = fn[i]
        c = cc[i]
        AoL[i] = (n⋅n)/(n⋅(c[2] - c[1]))
    end
    return AoL
end

HomogeneousMesh(d::Dict) = HomogeneousMesh(d[:meshfilename])

function cell_connectivity(cells::Vector{<:TriangleCell},bfaces,bc_face,nodes)
    NC = length(cells)
    NBF = length(bfaces)

    faces = Vector{Face2Node{2}}(0)
    c2f = Vector{Tuple{Int,Int,Int}}(NC)
    bf2c = Vector{BFace2Cell}(NBF)
    bfn = Vec2DArray{Float64}(NBF)
    bfc = Vec2DArray{Float64}(NBF)
    bccenter = Vec2DArray{Float64}(NBF)
    bcv = Vector{Float64}(NBF)

    nnbc = 0
    for (i,cell) in enumerate(cells)
        t1,t2,t3 = cell_faces(cell)

        t1notinfaces=true #future work: make this a tnotinfaces = zeros(Bool,NN)
        t2notinfaces=true
        t3notinfaces=true

        c2f_i1 = 0
        c2f_i2 = 0
        c2f_i3 = 0

        nnbc < NBF && for (j,el) in enumerate(bfaces)
            if t1 == el #future work make this a for t in tn, t == el
                c2f_i1 = -j
                bf2c[j] = BFace2Cell((i,))
                bfn[j] = area(t1,nodes)
                bfc[j] = center(t1,nodes)
                bcv[j] = 1/volume(cell,nodes)
                bccenter[j] = center(cell,nodes)
                t1notinfaces = false
                nnbc+=1
            elseif t2 == el
                c2f_i2 = -j
                bf2c[j] = BFace2Cell((i,))
                bfn[j] = area(t2,nodes)
                bfc[j] = center(t2,nodes)
                bcv[j] = 1/volume(cell,nodes)
                bccenter[j] = center(cell,nodes)
                t2notinfaces = false
                nnbc+=1
            elseif t3 == el
                c2f_i3 = -j
                bf2c[j] = BFace2Cell((i,))
                bfn[j] =  area(t3,nodes)
                bfc[j] = center(t3,nodes)
                bcv[j] = 1/volume(cell,nodes)
                bccenter[j] = center(cell,nodes)
                t3notinfaces = false
                nnbc+=1
            end
        end

        for (j,el) in enumerate(faces)
            if (t1[2], t1[1]) == el
                c2f_i1 = j
                t1notinfaces = false
            elseif (t2[2], t2[1]) == el
                c2f_i2 = j
                t2notinfaces = false
            elseif (t3[2], t3[1]) == el
                c2f_i3 = j
                t3notinfaces = false
            end
        end

        if t1notinfaces
            push!(faces,t1)
            c2f_i1 = length(faces)
        end

        if t2notinfaces
            push!(faces,t2)
            c2f_i2 = length(faces)
        end

        if t3notinfaces
            push!(faces,t3)
            c2f_i3 = length(faces)
        end

        @assert (c2f_i1 != 0) & (c2f_i2 != 0) & (c2f_i3 != 0)
        c2f[i] = (c2f_i1, c2f_i2, c2f_i3)
    end
    return faces, c2f, bf2c, bfn, bfc, bcv, bccenter
end

function face_connectivity(faces::Vector{<:Face2Node},c2f,cells::Vector{T},nodes) where {T<:AbstractCell2Node}
    NF = length(faces)
    NC = length(c2f)
    #VecType = nodetype(T)
    f2c = Vector{Face2Cell}(NF)
    fn = Vec2DArray{Float64}(NF)
    fc = Vec2DArray{Float64}(NF)
    cv = Vector{NTuple{2,Float64}}(NF)
    ccenter = Vector{NTuple{2,Vec2D{Float64}}}(NF)

    for i=1:NF
        notfound = true
        j=0
        cell1 = UInt(0)
        cell2 = UInt(0)

        while notfound
            j+=1
            if i in c2f[j]
                notfound = false
                cell1 = UInt(j)
            end
        end

        for l=(j+1):NC
            if i in c2f[l]
                cell2 = UInt(l)
                break
            end
        end

        @assert (cell1 != 0) & (cell2 != 0)
        f2c[i] = Face2Cell((cell1, cell2))
        fn[i] =  area(faces[i],nodes)
        fc[i] =  center(faces[i],nodes)
        cv[i] =  (1/volume(cells[cell1],nodes), 1/volume(cells[cell2],nodes))
        ccenter[i] =  (center(cells[cell1],nodes), center(cells[cell2],nodes))

    end

    return f2c, fn, fc, cv, ccenter
end


text:
# Julia wrapper for header: /usr/include/libavutil/mem.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_malloc,
    av_malloc_array,
    av_realloc,
    av_free,
    av_mallocz,
    av_mallocz_array,
    av_strdup,
    av_freep,
    av_memcpy_backptr

function av_malloc(size::Csize_t)
    ccall((:av_malloc,libavutil),Ptr{Void},(Csize_t,),size)
end

function av_malloc_array(nmemb::Csize_t,size::Csize_t)
    ccall((:av_malloc_array,libavutil),Ptr{Void},(Csize_t,Csize_t),nmemb,size)
end

function av_realloc(ptr,size::Csize_t)
    ccall((:av_realloc,libavutil),Ptr{Void},(Ptr{Void},Csize_t),ptr,size)
end

function av_free(ptr)
    ccall((:av_free,libavutil),Void,(Ptr{Void},),ptr)
end

function av_mallocz(size::Csize_t)
    ccall((:av_mallocz,libavutil),Ptr{Void},(Csize_t,),size)
end

function av_mallocz_array(nmemb::Csize_t,size::Csize_t)
    ccall((:av_mallocz_array,libavutil),Ptr{Void},(Csize_t,Csize_t),nmemb,size)
end

function av_strdup(s)
    ccall((:av_strdup,libavutil),Ptr{UInt8},(Ptr{UInt8},),s)
end

function av_freep(ptr)
    ccall((:av_freep,libavutil),Void,(Ptr{Void},),ptr)
end

function av_memcpy_backptr(dst,back::Integer,cnt::Integer)
    ccall((:av_memcpy_backptr,libavutil),Void,(Ptr{UInt8},Cint,Cint),dst,back,cnt)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


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

makedocs(;
    modules=[MyPackage2],
    sitename="Syslab库使用手册",
    # format = LaTeX(), #ToDo：目前生成PDF失败，待研究
    pages=[
        "Home" => "index.md",
        "Math" => "math.md",
        "Bar" => [
            "Bar1" => "bar1.md"
            "Bar2" => "bar2.md"
        ],
        "All Functions" => "functions.md"
    ],
)


text:
## This script was used to create keplerian_derivatives.jl

using Pkg
Pkg.activate("examples")
Pkg.instantiate()

using GLOM_RV_Example; GLOM_RV = GLOM_RV_Example
using SymEngine
using UnitfulAstro; using Unitful

@vars K, P, M0, h, k, γ, π_sym, t, cosEA, sinEA
@funs EA
e = sqrt(h^2+k^2)
EAf = EA(P, M0, h, k, t)
sinE = sin(EAf)
cosE = cos(EAf)
sinω = h/e
cosω = k/e
e_mod = 1 - e ^ 2
factor = 1 - e * cosE

function kep_EA_simplify(kep_sym::Basic)
    kep_sym = kep_sym(diff(EAf, P) => -2 * π_sym * t / (P * P * factor))
    kep_sym = kep_sym(diff(EAf, M0) => -1 / factor)
    kep_sym = kep_sym(diff(EAf, h) => sinω * sinE / factor)
    kep_sym = kep_sym(diff(EAf, k) => cosω * sinE / factor)
    return kep_sym
end
syms = [K, P, M0, h, k, γ]
n_kep = length(syms)

kep = K / factor * (e_mod * cosE * cosω - sqrt(e_mod) * sinE * sinω) + γ
dkeps = [kep_EA_simplify(diff(kep, sym)) for sym in syms]
ddkeps = [[kep_EA_simplify(diff(dkep, sym)) for dkep in dkeps] for sym in syms]

function kep_simplify!(kep::Basic, to_be_replaced::Basic, to_replace::Basic)
	kep = kep(to_be_replaced => to_replace)
	for i in 1:n_kep
		dkeps[i] = dkeps[i](to_be_replaced => to_replace)
		for j in 1:n_kep
			ddkeps[i][j] = ddkeps[i][j](to_be_replaced => to_replace)
		end
	end
	return kep
end

@vars esq, e, jsq, j, qmod

kep = kep_simplify!(kep, cosE, cosEA)
kep = kep_simplify!(kep, sinE, sinEA)
kep = kep_simplify!(kep, h^2+k^2, esq)
kep = kep_simplify!(kep, sqrt(esq), e)
kep = kep_simplify!(kep, 1 - esq, jsq)
kep = kep_simplify!(kep, sqrt(jsq), j)
kep = kep_simplify!(kep, 1 - e * cosEA, qmod)

begin
	file_loc = "src/keplerian_derivatives.jl"
	io = open(file_loc, "w")

	# begin to write the function including assertions that the amount of hyperparameters are correct
	write(io, """using UnitfulAstro; using Unitful

	\"\"\"
	kep_deriv() function created by examples/create_keplerian_derivatives.jl.
	Derivative of a Keplerian radial velocity signal using h and k instead of
	e and ω

	Parameters:
	K (Unitful.Velocity): velocity semi-amplitude
	P (Unitful.Time): period
	M0 (real): initial mean anomaly
	h (real): eccentricity * sin(argument of periastron)
	k (real): eccentricity * cos(argument of periastron)
	γ (Unitful.Velocity): velocity offset
	t (Unitful.Time): time
	dorder (vector of integers): A vector of how many partial derivatives you want
		to take with respect to each variable in the input order

	Returns:
	Unitful Quantity: The derivative specified with dorder
	\"\"\"
	function kep_deriv(
		K::Unitful.Velocity,
		P::Unitful.Time,
		M0::Real,
		h::Real,
		k::Real,
		γ::Unitful.Velocity,
		t::Unitful.Time,
		dorder::Vector{<:Integer})

		validate_kepler_dorder(dorder)

		esq = h * h + k * k
		e = sqrt(esq)
		jsq = 1 - esq
		j = sqrt(jsq)
		EAval = ecc_anomaly(t, P, M0, e)
		cosEA = cos(EAval)
		sinEA = sin(EAval)
		qmod = 1 - e * cosEA

	""")
	max_derivs = 3

	tot_inds =[[i] for i in 0:n_kep]
	append!(tot_inds, collect(Iterators.flatten([[[i,j] for i in 1:j] for j in 1:n_kep])))

	for inds in tot_inds
		dorder = zeros(Int, n_kep)
		if inds[1] != 0
			dorder[inds[1]] = 1
			func = dkeps[inds[1]]
		else
			func = kep
		end
		if length(inds) == 2
			dorder[inds[2]] += 1
			func = ddkeps[inds[1]][inds[2]]
		end

		func_str = SymEngine.toString(func)
		func_str = replace(func_str, "π_sym"=>"π")
		write(io, "    if dorder ==" * string(dorder) * "\n")
		write(io, "        func = " * func_str * "\n    end\n\n")
	end

	write(io, "    return float(func)\n\nend\n")

	write(io, """kep_deriv(ks::kep_signal, t::Unitful.Time, dorder::Vector{<:Integer}) =\n    kep_deriv(ks.K, ks.P, ks.M0, ks.h, ks.k, ks.γ, t, dorder)""")

	close(io)
end


text:
# Autogenerated wrapper script for racket_jll for x86_64-linux-gnu
export racket, raco

JLLWrappers.@generate_wrapper_header("racket")
JLLWrappers.@declare_executable_product(racket)
JLLWrappers.@declare_executable_product(raco)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        racket,
        "bin/racket",
    )

    JLLWrappers.@init_executable_product(
        raco,
        "bin/raco",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@recipe function plot(p::AbstractParticles)
    seriestype --> :histogram
    @series p.particles
end

# @recipe f(::Type{<:AbstractParticles}, p::AbstractParticles) = p.particles # Does not seem to be needed

function handle_args(p)
    length(p.args) < 2 && throw(ArgumentError("This function is called with at least two arguments (x, y, ..."))
    x,y = p.args[1:2]
    y isa AbstractArray{<:AbstractParticles} || throw(ArgumentError("The second argument must be a vector of some kind of Particles"))
    x,y
end

function quantiles(y,q::Number)
    m = mean.(y)
    q > 0.5 && (q = 1-q)
    lower = -(quantile.(y,q)-m)
    upper = quantile.(y,1-q)-m
    lower,upper
end

function quantiles(y,q)
    m = mean.(y)
    lower = -(quantile.(y,q[1])-m)
    upper = quantile.(y,q[2])-m
    lower,upper
end

@userplot Errorbarplot
@recipe function plt(p::Errorbarplot)
    x,y = handle_args(p)
    q = length(p.args) >= 3 ? p.args[3] : 0.025
    m = mean.(y)
    lower,upper = quantiles(y, q)
    label --> "Mean with $q quantile"
    yerror := (lower,upper)
    x,m
end

@userplot MCplot
@recipe function plt(p::MCplot)
    x,y = handle_args(p)
    label --> ""
    alpha --> 1/log(nparticles(y))
    m = Matrix(y)'
    x,m
end

@userplot Ribbonplot
@recipe function plt(p::Ribbonplot)
    x,y = handle_args(p)
    q = length(p.args) >= 3 ? p.args[3] : 0.025
    label --> "Mean with $q quantile"
    m = mean.(y)
    lower,upper = quantiles(y, q)
    ribbon := [lower,upper]
    x,m
end

"""
errorbarplot(x,y,[q=0.025])

Plots a vector of particles with error bars at quantile `q`.
If `q::Tuple`, then you can specify both lower and upper quantile, e.g., `(0.01, 0.99)`.
"""
errorbarplot

"""
mcplot(x,y,[q=0.025])

Plots all trajectories represented by a vector of particles
"""
mcplot

"""
ribbonplot(x,y,[q=0.025])

Plots a vector of particles with a ribbon covering quantiles `q, 1-q`.
If `q::Tuple`, then you can specify both lower and upper quantile, e.g., `(0.01, 0.99)`.
"""
ribbonplot

@recipe function plt(y::MvParticles, q=0.025)
    label --> "Mean with ($q, $(1-q)) quantiles"
    lower,upper = quantiles(y, q)
    ribbon := [lower,upper]
    mean.(y)
end

@recipe function plt(func::Function, x::MvParticles, q=0.025)
    y = func.(x)
    label --> "Mean with ($q, $(1-q)) quantiles"
    xerror := quantiles(x, q)
    yerror := quantiles(y, q)
    mean.(x), mean.(y)
end

@recipe function plt(x::MvParticles, y::MvParticles, q=0.025; points=false)
    my = mean.(y)
    mx = mean.(x)
    if points
        @series begin
            seriestype --> :scatter
            primary := false
            alpha --> 0.1
            Matrix(x), Matrix(y)
        end
    else
        yerror := quantiles(y, q)
        xerror := quantiles(x, q)
        label --> "Mean with $q quantile"
    end
    mx, my
end

@recipe function plt(x::MvParticles, y::AbstractArray, q=0.025)
    mx = mean.(x)
    lower,upper = quantiles(x, q)
    xerror := (lower,upper)
    mx, y
end

@recipe function plt(x::AbstractArray, y::MvParticles, q=0.025)
    ribbon := [quantiles(y, q)...]
    label --> "Mean with ($q, $(1-q)) quantiles"
    x, mean.(y)
end


text:
# import TrajectoryOptimization: dynamics
# import TrajectoryOptimization: Model, LQRCost, Problem, Objective, rollout!, iLQRSolverOptions,
#     AbstractSolver, jacobian!, _backwardpass!, _backwardpass_sqrt!, AugmentedLagrangianSolverOptions, ALTROSolverOptions,
#     goal_constraint, update_constraints!, update_active_set!, jacobian!, update_problem,
#     line_trajectory, total_time, generate_jacobian, _check_dynamics, AnalyticalModel, _test_jacobian,
#     f_augmented!, add_slack_controls, add_min_time_controls
#
# ## Augment dynamics
# model = Dynamics.pendulum
# model_d = discretize_model(model,:rk4)
#
# # Add slack controls
# model_inf = add_slack_controls(model_d)
#
# @test model_inf.n == model_d.n
# @test model_inf.m == model_d.m + model_d.n
#
# x = rand(model_d.n)
# u = rand(model_d.m)
# u_inf = rand(model_inf.m); u_inf[1:model_d.m] .= u
# x_d = zeros(model_d.n)
# x_inf = zeros(model_inf.n)
# Z_d = zeros(model_d.n,model_d.n+model_d.m+1)
# Z_inf = zeros(model_inf.n,model_inf.n+model_inf.m+1)
#
# evaluate!(x_d,model_d,x,u,1.0)
# evaluate!(x_inf,model_inf,x,u_inf,1.0)
#
# evaluate!(Z_d,rand(model_d.n),model_d,x,u,1.0)
# model_inf.∇f(Z_inf,x,u_inf,1.0)
#
# @test x_inf == x_d + u_inf[model_d.m+1:end]
# @test Z_d == Z_inf[1:model_d.n,[(1:model_d.n+model_d.m)...,2model_d.n+model_d.m+1]]
# @test Z_inf[1:model_d.n,(1:model_d.n) .+ (model_d.n+model_d.m)] == Diagonal(1.0I,model_d.n)
#
# # Add minimum time controls
# model_mintime = add_min_time_controls(model_d)
#
# @test model_mintime.n == model_d.n + 1
# @test model_mintime.m == model_d.m + 1
#
# dt = 0.75
# x = rand(model_d.n)
# u = rand(model_d.m)
# u_mintime = sqrt(dt)*ones(model_mintime.m); u_mintime[1:model_d.m] .= u
# x_d = zeros(model_d.n)
# x_mintime = rand(model_mintime.n); x_mintime[1:model_d.n] .= x
# x₊_mintime = zeros(model_mintime.n)
# Z_d = zeros(model_d.n,model_d.n+model_d.m+1)
# Z_mintime = zeros(model_mintime.n,model_mintime.n+model_mintime.m)
#
# evaluate!(x_d,model_d,x,u,dt)
# evaluate!(x₊_mintime,model_mintime,x_mintime,u_mintime,Inf) # note we can use Inf here since dt input is not used
#
# evaluate!(Z_d,rand(model_d.n),model_d,x,u,dt)
# model_mintime.∇f(Z_mintime,x_mintime,u_mintime,Inf)
#
# @test isapprox(x₊_mintime[1:model_d.n],x_d)
#
# @test isapprox(x₊_mintime[1:model_d.n],x_d)
# @test x₊_mintime[end] == u_mintime[end]
# @test isapprox(Z_d[1:model_d.n,1:model_d.n+model_d.m],Z_mintime[1:model_d.n,[(1:model_d.n)...,((1:model_d.m) .+ (model_d.n+1))...]])
# @test isapprox(Z_mintime[1:model_d.n,model_d.n+model_d.m+2],Z_d[:,end]*2.0*u_mintime[end])
# @test Z_mintime[model_d.n+1,model_d.n+model_d.m+2] == 1.0
#
# # Add infeasible and minimum time controls
# model_altro = add_min_time_controls(model_inf)
#
# dt = 0.75
# x = rand(model_d.n)
# u = rand(model_d.m)
# u_altro = sqrt(dt)*ones(model_altro.m); u_altro[1:model_d.m] .= u
# u_altro
# x_d = zeros(model_d.n)
# x_altro = rand(model_altro.n); x_altro[1:model_d.n] .= x
# x₊_altro = zeros(model_altro.n)
# Z_d = zeros(model_d.n,model_d.n+model_d.m+1)
# Z_altro = zeros(model_altro.n,model_altro.n+model_altro.m)
#
# @test model_altro.n == model_d.n + 1
# @test model_altro.m == model_d.m + model_d.n + 1
#
# evaluate!(x_d,model_d,x,u,dt)
# evaluate!(x₊_altro,model_altro,x_altro,u_altro,Inf) # note we can use Inf here since dt input is not used
#
# evaluate!(Z_d,rand(model_d.n),model_d,x,u,dt)
# model_altro.∇f(Z_altro,x_altro,u_altro,Inf)
#
# @test isapprox(x₊_altro[1:model_d.n],x_d + u_altro[(1:model_d.n) .+ model_d.n])
# @test x₊_altro[end] == u_altro[end]
# @test isapprox(Z_d[1:model_d.n,1:model_d.n+model_d.m],Z_altro[1:model_d.n,[(1:model_d.n)...,((1:model_d.m) .+ (model_d.n+1))...]])
# @test isapprox(Z_altro[1:model_d.n,model_altro.n+model_altro.m],Z_d[:,end]*2.0*u_altro[end])
# @test Z_altro[model_d.n+1,model_altro.n+model_altro.m] == 1.0

import TrajectoryOptimization: dynamics
import TrajectoryOptimization: Model, LQRCost, Problem, Objective, rollout!, iLQRSolverOptions,
    AbstractSolver, jacobian!, _backwardpass!, _backwardpass_sqrt!, AugmentedLagrangianSolverOptions, ALTROSolverOptions,
    goal_constraint, update_constraints!, update_active_set!, jacobian!, update_problem,
    line_trajectory, total_time, generate_jacobian, _check_dynamics, AnalyticalModel, _test_jacobian,
    f_augmented!, add_slack_controls, add_min_time_controls, Nominal

## Augment dynamics
model = Dynamics.pendulum
model_d = discretize_model(model,:rk4)

# Add slack controls
model_inf = add_slack_controls(model_d)

@test model_inf.n == model_d.n
@test model_inf.m == model_d.m + model_d.n

x = rand(model_d.n)
u = rand(model_d.m)
u_inf = rand(model_inf.m); u_inf[1:model_d.m] .= u
x_d = zeros(model_d.n)
x_inf = zeros(model_inf.n)
Z_d = zeros(model_d.n,model_d.n+model_d.m+1)
Z_inf = zeros(model_inf.n,model_inf.n+model_inf.m+1)

evaluate!(x_d,model_d,x,u,1.0)
evaluate!(x_inf,model_inf,x,u_inf,1.0)

jacobian!(Z_d,model_d,x,u,1.0)
model_inf.∇f(Z_inf,x,u_inf,1.0)

@test x_inf == x_d + u_inf[model_d.m+1:end]
@test Z_d == Z_inf[1:model_d.n,[(1:model_d.n+model_d.m)...,2model_d.n+model_d.m+1]]
@test Z_inf[1:model_d.n,(1:model_d.n) .+ (model_d.n+model_d.m)] == Diagonal(1.0I,model_d.n)

# Add minimum time controls
model_mintime = add_min_time_controls(model_d)

@test model_mintime.n == model_d.n + 1
@test model_mintime.m == model_d.m + 1

dt = 0.75
x = rand(model_d.n)
u = rand(model_d.m)
u_mintime = sqrt(dt)*ones(model_mintime.m); u_mintime[1:model_d.m] .= u
x_d = zeros(model_d.n)
x_mintime = rand(model_mintime.n); x_mintime[1:model_d.n] .= x
x₊_mintime = zeros(model_mintime.n)
Z_d = zeros(model_d.n,model_d.n+model_d.m+1)
Z_mintime = zeros(model_mintime.n,model_mintime.n+model_mintime.m)

evaluate!(x_d,model_d,x,u,dt)
evaluate!(x₊_mintime,model_mintime,x_mintime,u_mintime,Inf) # note we can use Inf here since dt input is not used

jacobian!(Z_d,model_d,x,u,dt)
model_mintime.∇f(Z_mintime,x_mintime,u_mintime,Inf)

@test isapprox(x₊_mintime[1:model_d.n],x_d)

@test isapprox(x₊_mintime[1:model_d.n],x_d)
@test x₊_mintime[end] == u_mintime[end]
@test isapprox(Z_d[1:model_d.n,1:model_d.n+model_d.m],Z_mintime[1:model_d.n,[(1:model_d.n)...,((1:model_d.m) .+ (model_d.n+1))...]])
@test isapprox(Z_mintime[1:model_d.n,model_d.n+model_d.m+2],Z_d[:,end]*2.0*u_mintime[end])
@test Z_mintime[model_d.n+1,model_d.n+model_d.m+2] == 1.0

# Add infeasible and minimum time controls
model_altro = add_min_time_controls(model_inf)

dt = 0.75
x = rand(model_d.n)
u = rand(model_d.m)
u_altro = sqrt(dt)*ones(model_altro.m); u_altro[1:model_d.m] .= u
x_d = zeros(model_d.n)
x_altro = rand(model_altro.n); x_altro[1:model_d.n] .= x
x₊_altro = zeros(model_altro.n)
Z_d = zeros(model_d.n,model_d.n+model_d.m+1)
Z_altro = zeros(model_altro.n,model_altro.n+model_altro.m)

@test model_altro.n == model_d.n + 1
@test model_altro.m == model_d.m + model_d.n + 1

evaluate!(x_d,model_d,x,u,dt)
evaluate!(x₊_altro,model_altro,x_altro,u_altro,Inf) # note we can use Inf here since dt input is not used

jacobian!(Z_d,model_d,x,u,dt)
model_altro.∇f(Z_altro,x_altro,u_altro,Inf)

@test isapprox(x₊_altro[1:model_d.n],x_d + u_altro[(1:model_d.n) .+ model_d.n])
@test x₊_altro[end] == u_altro[end]
@test isapprox(Z_d[1:model_d.n,1:model_d.n+model_d.m],Z_altro[1:model_d.n,[(1:model_d.n)...,((1:model_d.m) .+ (model_d.n+1))...]])
@test isapprox(Z_altro[1:model_d.n,model_altro.n+model_altro.m],Z_d[:,end]*2.0*u_altro[end])
@test Z_altro[model_d.n+1,model_altro.n+model_altro.m] == 1.0

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
## Caculate Trace of Hat matrix in LMM (only for PLSGeneral type in MixedModels)

function Λ0(m::LinearMixedModel)
	vv = SparseMatrixCSC{Float64,Int}[]
	for i in 1:length(m.λ)
		nl = length(m.facs[i].pool)
		ll = m.λ[i]
		#p = size(ll,1)
		p = size(ll)[1]
		if p == 1
			isa(ll,MixedModels.PDScalF) || error("1×1 λ section should be a PDScalF type")
			push!(vv,ll.s .* speye(nl))
		else
			isa(ll, MixedModels.PDDiagF) ? sfll = sparse(full(ll.d)): sfll = sparse(full(ll.ch.UL))
			push!(vv,blkdiag([sfll for _ in 1:nl]...))
		end
	end
	Triangular(blkdiag(vv...),:L,false)
end

function traceHat(m::LinearMixedModel)
	X=m.X.m
	XtX = Symmetric(X'X,:L)
	Zt=MixedModels.zt(m)
	Λt=Λ0(m).data

	Ztc = CholmodSparse(Λt*Zt,0)
	L = cholfact(Ztc,1.,true)

	perm=L.Perm .+ one(eltype(L.Perm))
	RX=cholfact(MixedModels.symcontents(XtX),:L)
	RZX = Λt*Zt*X

	copy!(RZX,solve(L, RZX[perm,:], CHOLMOD_L))

	LAPACK.potrf!('L',BLAS.syrk!('L','T',-1.,RZX, 1.,
		copy!(RX.UL, MixedModels.symcontents(XtX))))

	Cl= solve(L,solve(L,Ztc,CHOLMOD_P),CHOLMOD_L)
	#Cl= solve(L,Ztc,CHOLMOD_L)
	Cr=Ac_ldiv_B(tril(RX.UL), X' - full(Ac_mul_B(Cl,RZX))' )

	rho=sumabs2(Cl.nzval)+sumabs2(vec(Cr))
	rho
end

function vinv(m::LinearMixedModel)

	Zt=MixedModels.zt(m)
	#Λ=SelectMixedModels.Λ0(m)

    L1 = cholfact(CholmodSparse(SelectMixedModels.Λ0(m).data,0),true)
    T=solve(L1,CholmodSparse(speye(size(L1)[1]),0), CHOLMOD_A)
	L2 = cholfact(CholmodSparse(sparse(T) + Zt *Zt'),true)
	S = solve(L2,CholmodSparse(speye(size(L2)[1]),0), CHOLMOD_A)
	return  speye(size(m)[1]) - Zt' * sparse(S) * Zt
end

max_stars_repo_name: ilovealldata/SelectMixedModels

max_issues_repo_name: ilovealldata/SelectMixedModels

max_forks_repo_name: ilovealldata/SelectMixedModels


text:
@testset "various types" begin
    for var in (vars..., tuples...)
        @inferred randobs(var)
        @inferred randobs(var, 4)
        @test typeof(randobs(var)) == typeof(getobs(var, 1))
        @test typeof(randobs(var, 4)) == typeof(getobs(var, 1:4))
        @test nobs(randobs(var, 4)) == nobs(getobs(var, 1:4))
    end
end

@testset "check that sampled obs exist" begin
    for tX in (X, DataSubset(X))
        X_rnd = @inferred(randobs(tX, 30))
        for i = 1:30
            @testset "random obs $i" begin
                found = false
                for j = 1:150
                    if all(X_rnd[:,i] .== X[:,j])
                        found = true
                    end
                end
                @test found
            end
        end
    end
end

# test if obs in tuple match each other
@testset "Tuple" begin
    for i = 1:30
        @test 0 < randobs(CustomType()) <= 150
        x1, y1 = randobs((X1,Y1))
        @test all(x1 .== y1)
        x1, y1 = randobs((X1,Y1), 5)
        @test all(x1' .== y1)
    end
end


text:
@testset "pentland" begin
    img = generate_surface(SynthSphere(), 0.5, [0,0,1], radius = 50, background = false)
    Z = retrieve_surface(Pentland(), img)
    @test Z[76,76] ≈ 3.4232738028272176
    @test isapprox(Z[1,1], 1.160438047437746e-14; atol = 2e-15)
    @test isapprox(Z[151,151], 6.374017915783071e-15; atol = 2e-15)
    @test Z[90,90] ≈ 3.4196712241892793
    @test Z[101,36] ≈ 2.342218759464893
    @test isapprox(Z[22, 136], 9.89746106229774e-15; atol = 2e-15)
    @test Z[42,42] ≈ 2.9211132774577244
end


text:
const flipmap = [0x02, 0x01, 0x00, 0x03]
const onebitmap = [0x00 => 0x01, 0x00 => 0x02, 0x00 => 0x03,
    0xff => 0xff, 0x01 => 0x02, 0x01 => 0x03, 0xff => 0xff, 0xff => 0xff, 0x02 => 0x03]

"""
    get_genotypes!(buf, p, v; ldbuf)

Computes unphased biallelic dosage of ALT allele.

- `buf`: stores genotype values.
- `p`: a `Pgen` object.
- `v`: a `Variant` object.
- `ldbuf`: most recent non-LD-compressed genotypes.

Returns:
- `buf`: resulting genotype values
- `variant_record`: current record for variant
- `offset`: end of dosage record on the current variant record track.
"""
function get_genotypes!(buf::Vector{UInt8}, p::Pgen, v::Variant;
    ldbuf::Union{Nothing, Vector{UInt8}}=nothing)
    compression_type = v.record_type & 0x07
    @assert v.record_type & 0x08 == 0 "Multiallelic case unsupported"
    if (compression_type == 0x02 || compression_type == 0x03)
        if ldbuf === nothing
            # load most recent non-LD-compressed dosage
            get_genotypes!(buf, p, p.header.most_recent_non_ld[v.index])
        else
            buf .= ldbuf
        end
    end
    if p.variant_record_cache !== nothing
        seek(p.io, v.offset)
        p.variant_record_cache .= read(p.io, v.length)
        variant_record = p.variant_record_cache
    else
        variant_record = @view p.data[Int(v.offset + 1) : Int(v.offset + v.length)]
    end
    if compression_type == 0x00
        offset = _get_genotypes_no_compression!(buf, p, variant_record)
    elseif compression_type == 0x01
        offset = _get_genotypes_1bit!(buf, p, variant_record)
    elseif compression_type == 0x02 || compression_type == 0x03
        offset = _get_genotypes_difflist!(buf, p, variant_record)
        if compression_type == 0x03
            for i in 1:p.header.n_samples
                buf[i] = flipmap[buf[i] + 0x01]
            end
        end
    elseif compression_type == 0x04 || compression_type == 0x06 || compression_type == 0x07
        if compression_type == 0x04
            fill!(buf, 0x00)
        elseif compression_type == 0x06
            fill!(buf, 0x02)
        else
            fill!(buf, 0x03)
        end
        offset = _get_genotypes_difflist!(buf, p, variant_record)
    else
        @error "invalid compression type"
    end
    buf, variant_record, offset
end

"""
    get_genotypes(p, v)

Computes unphased biallelic dosage of ALT allele.

- `p`: a `Pgen` object.
- `v`: a `Variant` object.

Returns:
- `buf`: resulting genotype values
- `variant_record`: current record for variant
- `offset`: end of dosage record on the current variant record track.
"""
function get_genotypes(p::Pgen, v::Variant)
    buf = Vector{UInt8}(undef, p.header.n_samples)
    get_genotypes!(buf, p, v)
end

# for compression mode "0x00", no compression
function _get_genotypes_no_compression!(buf::Vector{UInt8}, p::Pgen, variant_record::AbstractVector{UInt8})
    n_samples = p.header.n_samples
    n_bytes = (n_samples + 3) >> 2
    genotypes_raw_cache = @view variant_record[1:n_bytes]
    @inbounds for i in 1:n_samples
        ip3 = i + 3
        buf[i] = (genotypes_raw_cache[ip3 >> 2] >> ((ip3 & 0x03) << 1)) & 0x03
    end
    n_bytes
end

function _get_genotypes_1bit!(buf::Vector{UInt8}, p::Pgen, variant_record::AbstractVector{UInt8})
    n_samples = p.header.n_samples
    n_bytes = (n_samples + 7) >> 3
    falseval, trueval = onebitmap[variant_record[1]]
    bv = BitsVector(pointer(variant_record, 2), 1, n_samples)
    @inbounds for i in 1:n_samples
        buf[i] = bv[i] == 0x01 ? trueval : falseval
    end
    dl, offset = parse_difflist(variant_record, UInt(n_bytes + 1), p.header.bytes_per_sample_id, true)
    _get_difflist_genotypes!(buf, p, dl)
    offset
end

function _get_genotypes_difflist!(buf::Vector{UInt8}, p::Pgen, variant_record::AbstractVector{UInt8})
    dl, offset = parse_difflist(variant_record, zero(UInt), p.header.bytes_per_sample_id, true)
    _get_difflist_genotypes!(buf, p, dl)
    offset
end

function _get_difflist_genotypes!(buf::Vector{UInt8}, p::Pgen, dl::DiffList)
    ngroups = (dl.len + 63) ÷ 64
    for gid in 1:ngroups
        parse_difflist_sampleids!(p.difflist_cache, p.difflist_cache_incr, dl, gid)
        for (idx, sampleid) in enumerate(p.difflist_cache)
            totalidx = 64 * (gid - 1) + idx
            if totalidx > dl.len
                break
            end
            buf[sampleid] = dl.genotypes[totalidx]
        end
    end
end


text:
"""
    AbstractClock{T<:Real}

Supertype of clocks with time of type `T`.
"""
abstract type AbstractClock{T<:Real} end

"""
    currenttime(c::AbstractClock)

Return current time of clock `c`.
"""
function currenttime end

"""
    limit(c::AbstractClock)

Return the limit of clock `c`. For a ContinuousClock `c`,
the max time might larger than `limit(c)`.
"""
function limit end

"""
    start(c::AbstractClock)

Return the start time of clock `c`.
"""
function start end

"""
    init!(c::AbstractClock)

Update current time to start time.
"""
function init! end

"""
    DiscreteClock([start], timelist)
    DiscreteClock(stop)

Clock for discrete-time process, time of which muse be increased with given step and
can't
be updated manually. The `timelist` must be a non-empty and monotonically increasing
`AbstractVector`.  If the `start` is not specified, the first item of `timelist` will be
deleted and set as `start`. During iteration, the current time will be updated automatically
and returned as iteration item. When the iteration finished without `break`, [`init!`](@ref)
will be applied.  `DiscreteClock(stop)` will create a clock with `start=0` and
`timelist=Base.OneTo(stop)`

# Examples

```jldoctest
julia> c = DiscreteClock(0:3);

julia> currenttime(c)
0

julia> [(t, currenttime(c)) for t in c]
3-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (2, 2)
 (3, 3)

julia> currenttime(c)
0

julia> c = DiscreteClock(3); # similar to DiscreteClock(0:3)

julia> (currenttime(c), collect(c))
(0, [1, 2, 3])
```
"""
struct DiscreteClock{T<:Real,I<:AbstractVector{T}} <: AbstractClock{T}
    current::Array{T,0}
    start::T
    timelist::I
    function DiscreteClock(start::T, timelist::I) where {T,I<:AbstractVector{T}}
        return new{T,I}(fill(start), start, timelist)
    end
end
function DiscreteClock(timelist::AbstractVector)
    isempty(timelist) && throw(ArgumentError("timelist must contain at least one element"))
    return DiscreteClock(timelist[1], timelist[2:end])
end
function DiscreteClock(stop::Integer)
    stop > 0 || throw(ArgumentError("stop time must be > 0"))
    return DiscreteClock(zero(stop), Base.OneTo(stop))
end

# iterator interfaces
Base.IteratorSize(::Type{<:DiscreteClock}) = Base.HasLength()
Base.length(c::DiscreteClock) = length(c.timelist)
Base.eltype(::Type{<:DiscreteClock{T}}) where {T} = T
Base.iterate(c::DiscreteClock) = _itr_update!(c, iterate(c.timelist))
function Base.iterate(c::DiscreteClock, state)
    return _itr_update!(c, iterate(c.timelist, state))
end
_itr_update!(c::DiscreteClock, ::Nothing) = (init!(c); nothing)
function _itr_update!(c::DiscreteClock{T,I}, ret::Tuple{T,Any}) where {T,I}
    c.current[] = ret[1]
    return ret
end

# clock interfaces
currenttime(c::DiscreteClock) = c.current[]
limit(c::DiscreteClock) = last(c.timelist)
init!(c::DiscreteClock) = c.current[] = c.start
start(c::DiscreteClock) = c.start

"""
    ContinuousClock{T, I<:Union{Nothing, DiscreteClock}} <: AbstractClock{T}
    ContinuousClock(stop, [start=zero(stop)]; [max_epoch=nothing])

A clock for continuous-time process. Unlike the [`DiscreteClock`](@ref DiscreteClock),
during iteration, the current time will not be update automatically, but update by
[`increase!`](@ref increase!) manually. Besides the epoch of current iteration instead of
current time will be returned. If the `max_epoch` is specified, the iteration will break
when epoch reach to the `max_epoch`, even `currenttime(c) < limit(c)`, and break in this way the
[`init!(c)`](@ref init!) will not be applied.

# Examples

```jldoctest
julia> c = ContinuousClock(3.0; max_epoch=2);

julia> for epoch in c
           increase!(c, 1)
           println(currenttime(c), '\t', epoch)
       end
1.0	1
2.0	2

julia> for epoch in c
           increase!(c, 1)
           println(currenttime(c), '\t', epoch)
       end
3.0	1

julia> for epoch in c
           increase!(c, 1)
           println(currenttime(c), '\t', epoch)
       end
1.0	1
2.0	2
```
"""
struct ContinuousClock{T<:Real,I<:Union{Nothing,Integer}} <: AbstractClock{T}
    current::Array{T,0}
    start::T
    stop::T
    epoch::I
    function ContinuousClock(
        start::T,
        stop::T,
        epoch::I,
    ) where {T<:Real,I<:Union{Nothing,Integer}}
        start > stop && throw(ArgumentError("stop must be larger than start"))
        return new{T,I}(fill(start), start, stop, epoch)
    end
end
function ContinuousClock(
    stop::Real,
    start::Real=zero(stop);
    max_epoch::Union{Nothing,Integer}=nothing,
)
    return ContinuousClock(promote(start, stop)..., max_epoch)
end

# iterator interfaces
Base.IteratorSize(::Type{<:ContinuousClock}) = Base.SizeUnknown()
Base.eltype(::ContinuousClock{T}) where {T} = T
Base.iterate(c::ContinuousClock) =
    currenttime(c) < limit(c) ? _itr(c.epoch) : (init!(c); nothing)
Base.iterate(c::ContinuousClock, state) =
    currenttime(c) < limit(c) ? _itr(c.epoch, state) : (init!(c); nothing)

_itr(::Nothing, i::Int=1) = (i, i + 1)
_itr(lim::Integer, i::Integer=one(lim)) = ifelse(i > lim, nothing, (i, i + 1))

# clock interfaces
currenttime(c::ContinuousClock) = c.current[]
limit(c::ContinuousClock) = c.stop
start(c::ContinuousClock) = c.start
init!(c::ContinuousClock) = c.current[] = start(c)

"""
    increase!(c::ContinuousClock, t::Real)

Update current time of clock `c` to `currenttime(c) + t`.
"""
increase!(c::ContinuousClock{T}, t::Real) where {T<:Real} = increase!(c, convert(T, t))::T
increase!(c::ContinuousClock{T}, t::T) where {T<:Real} = c.current[] += t

# vim:tw=92:ts=4:sw=4:et


text:
using CSV
using Plots
using Random
using StatsFuns
using DataFrames
using Distributions
using LinearAlgebra
using TripletEmbeddings
using GreenVideoExperiments

# Y = Matrix(transpose(load("TaskA", "ground_truth", "30Hz")[begin:30:end,:data]))
Y = CSV.read(expanduser("~/Documents/Research/PRL2019/data/TaskA.csv"), DataFrame, header=false)
Y = Matrix(transpose(Y.Column1))[:,begin:30:end]

plt = plot(legend=:bottomleft)
plot!(plt, transpose(Y), label="Embedding", color=:black)

# for μ in [0.7, 0.8, 0.9]
for μ in [0.7]
    triplets = Triplets(Y; f = x -> logistic(10x))
    Ŷ = Embedding(size(Y))

    misclassifications = fit!(tSTE(α=30), triplets, Ŷ)
    Ŷ, tr = procrustes(Ŷ, Y)

    @show norm(Y - Ŷ).^2

    plot!(plt, transpose(Ŷ), label="Embedding μ = $(μ)")
end

display(plt)

################
### Debiased ###
################
let
    triplets = LabeledTriplets(Y; f = x -> logistic(10x))

    Ŷ = Embedding(size(Y))
    loss = ρSTE(ρ = Dict([-1, +1] .=> [0.3, 0.3]))
    misclassifications = fit!(loss, triplets, Ŷ; max_iterations=1000)

    Ŷ, tr = procrustes(Ŷ, Y)

    @show norm(Y - Ŷ).^2

    plot!(plt, transpose(Ŷ), label="Debiased")
end


text:
# all tests

using Test
import MathOptInterface
const MOI = MathOptInterface

import GLPK
glpk = MOI.OptimizerWithAttributes(
    GLPK.Optimizer,
    MOI.Silent() => true,
    "tol_int" => 1e-10,
    "tol_bnd" => 1e-10,
    "mip_gap" => 1e-10,
)

println("starting PajaritoExtras tests")
@testset "PajaritoExtras tests" begin
    include("JuMP_tests.jl")
    TestJuMP.runtests(glpk)
end