text:
function feature_string()
    llvmlib_path = VERSION ≥ v"1.6.0-DEV.1429" ? Base.libllvm_path() : only(filter(lib->occursin(r"LLVM\b", basename(lib)), Libdl.dllist()))
    libllvm = Libdl.dlopen(llvmlib_path)
    gethostcpufeatures = Libdl.dlsym(libllvm, :LLVMGetHostCPUFeatures)
    features_cstring = ccall(gethostcpufeatures, Cstring, ())
    features = filter(ext -> (ext ≠ "" && (m = match(r"\d", ext); isnothing(m) ? true : m.offset != 2)), split(unsafe_string(features_cstring), ','))
    features, features_cstring
end

const FEATURE_SET = Set{String}()
_has_feature(str::String) = str ∈ FEATURE_SET

archstr() = Sys.ARCH === :i686 ? "x86_64_" : string(Sys.ARCH) * '_'

# feature_name(ext) = archstr() * replace(ext[2:end], r"\." => "_")
feature_name(ext) = archstr() * ext[2:end]
process_feature(ext) = (feature_name(ext), first(ext) == '+')

has_feature(_) = False()
@noinline function set_featue(feature::String, has::Bool)
    featqn = QuoteNode(Symbol(feature))
    if has
        @eval has_feature(::Val{$featqn}) = True()
    else
        @eval has_feature(::Val{$featqn}) = False()
    end
end

function set_features!()
    features, features_cstring = feature_string()
    for ext ∈ features
      feature, has = process_feature(ext)
      has && push!(FEATURE_SET, feature)
      set_featue(feature, has)
    end
    Libc.free(features_cstring)
end
set_features!()

function reset_features!()
    features, features_cstring = feature_string()
    for ext ∈ features
        feature, has = process_feature(ext)
        if _has_feature(feature) ≠ has
            @debug "Defining $(has ? "presence" : "absense") of feature $feature."
            set_featue(feature, has)
        end
    end
    Libc.free(features_cstring)
end

register_size(::Type{T}) where {T} = register_size()
register_size(::Type{T}) where {T<:Union{Signed,Unsigned}} = simd_integer_register_size()

function define_cpu_name()
    cpu = QuoteNode(Symbol(Sys.CPU_NAME::String))
    @eval cpu_name() = Val{$cpu}()
end
define_cpu_name()


text:
#==
Instumenting Debugger Pass:

The purpose of this pass is to modify the IR code to instert debug statements.
One is inserted before each other statement in the IR.
Rough pseudocode for a Debug statment:
```
if should_break  # i.e. this breakpoint is active
    variables = filter(isdefined, slots)
    call break_action(variables)  # launch the debugging REPL etc.
end
```
The reality is a bit more complicated, as we actually workouit what is defined at
compile time.
==#

"""
    extended_insert_statements!(code, codelocs, stmtcount, newstmts)

Like `Cassette.insert_statements` but the `newstmts` function takes
3 arguments:
 - `statement`: the IR statement it will be replacing
 - `dest_i`: the first index that this will be inserted into, after renumbering to insert prior replacements.
     - in the `newst,ts` function `dest_i` should be used to calculate SSAValues and goto addresses
 - `src_i`: the index of the original IR statement that this will be replacing
     - You may wish to use this in the `newstmts` function if the code actually depends on which
       statement number of original IR is being replaced.
"""
function extended_insert_statements!(code, codelocs, stmtcount, newstmts)
    ssachangemap = fill(0, length(code))
    labelchangemap = fill(0, length(code))
    worklist = Tuple{Int,Int}[]
    @inbounds for i in 1:length(code)
        stmt = code[i]
        nstmts = stmtcount(stmt, i)
        if nstmts !== nothing
            addedstmts = nstmts - 1
            push!(worklist, (i, addedstmts))
            ssachangemap[i] = addedstmts
            if i < length(code)
                labelchangemap[i + 1] = addedstmts
            end
        end
    end
    Core.Compiler.renumber_ir_elements!(code, ssachangemap, labelchangemap)
    @inbounds for (src_i, addedstmts) in worklist
        dest_i = src_i + ssachangemap[src_i] - addedstmts # correct the index for accumulated offsets
        stmts = newstmts(code[dest_i], dest_i, src_i)
        @assert(length(stmts) == (addedstmts + 1), "$(length(stmts)) == $(addedstmts + 1)")
        code[dest_i] = stmts[end]
        for j in 1:(length(stmts) - 1) # insert in reverse to maintain the provided ordering
            insert!(code, dest_i, stmts[end - j])
            insert!(codelocs, dest_i, codelocs[dest_i])
        end
    end
end

"""
    solve_isdefined_map(ci, nargs)
Returns a vector which maps from statement index
to a list of slot_ids that are definately defined at that index
"""
solve_isdefined_map(reflection) = solve_isdefined_map(reflection.code_info, reflection.method.nargs)
function solve_isdefined_map(ci, nargs)
    @assert nargs >= 1  # must be at least one for #self
    cfg = Core.Compiler.compute_basic_blocks(ci.code)
    domtree = Core.Compiler.construct_domtree(cfg)

    isdefined_on = Vector{Vector{Int}}(undef, length(ci.code))
    function proc_block!(cur_block_ii, defined_slots)
        block = cfg.blocks[cur_block_ii]
        domnode = Core.Compiler.getindex(domtree.nodes, cur_block_ii)
        for stmt_ii in block.stmts.start : block.stmts.stop
            # We want to know what has been defined before we run this statement
            # i.e. what can occur on RHS of assigment or in a call
            isdefined_on[stmt_ii] = defined_slots

            # Now update `defined_slots` for future statements
            stmt = ci.code[stmt_ii]
            if isexpr(stmt, :(=)) && stmt.args[1] isa Core.SlotNumber
                id = stmt.args[1].id
                if id ∉ defined_slots
                    defined_slots = copy(defined_slots)
                    defined_slots = push!(defined_slots, id)
                end
            end
        end
        Core.Compiler.foreach(domnode.children) do block_ii
            proc_block!(block_ii, defined_slots)
        end
    end

    # Initial defined slots start from 2 rather than 1 too skip #self#
    proc_block!(1, collect(2:nargs))
    return isdefined_on
end

"""
    call_expr([mod], func, args...)

This function returns the IR expression for calling the given function, with the
given args. It is maked with `nooverdub` which will stop Cassette recursing into it.
"""
call_expr(mod::Module, func::Symbol, args...) = call_expr(GlobalRef(mod, func), args...)
call_expr(f, args...) = Expr(:call, Expr(:nooverdub, f), args...)

"""
    enter_debug_statements(slotnames, slot_created_ons, method::Method, ind::Int, orig_ind::Int)

This returns the IR code for a debug statement (as decribed at the top of this file).
This basically means creating code that checks if we `should_break` at this statement,
and if so works out what variable are defined, then passing those to the `break_action` call which will
show the debugging prompt.

 - slotnames: the names of the slots from the CodeInfo
 - slot_created_on: a vector saying where the variables were declared (as returned by `created_on`
 - method: the method being instruments
 - ind: the actual index in the code IR this is being  inserted at. This is where the SSAValues start from
 - orig_ind: the index in the original code IR for where this is being inserted. (before other debug statements were inserted above)
"""
@inline function enter_debug_statements(
    slotnames, isdefined_map, method::Method,
    stmt, ind::Int, orig_ind::Int
    )

    stmt_count = 6
    defined_slotids = isdefined_map[orig_ind]
    statements = Vector{Any}(undef, stmt_count)
    statements[1] = call_expr(MagneticReadHead, :should_break, method, orig_ind)
    statements[2] = Expr(:gotoifnot, Core.SSAValue(ind), ind + stmt_count - 1)  # go to last statement (i.e. the original stmt)
    statements[3] = Tuple(slotnames[defined_slotids])
    statements[4] = call_expr(Core, :tuple, Core.SlotNumber.(defined_slotids)...)
    names_ssa = Core.SSAValue(ind+2)
    values_ssa = Core.SSAValue(ind+3)
    statements[5] = call_expr(
        MagneticReadHead, :break_action,
        method,
        orig_ind,
        names_ssa, values_ssa
    )
    statements[6] = stmt  # last put in the original statement -- this is where we jump to
    #Core.println(statements)
    return statements
end

"""
    instrument!(::Type{<:HandEvalCtx}, reflection::Cassette.Reflection)
This is the transform for the debugger cassette pass.
it is the main method of this file, and calls all the ones defined earlier.
"""
@inline function instrument!(::Type{<:HandEvalCtx}, reflection::Cassette.Reflection)
    ir = reflection.code_info

    isdefined_map = solve_isdefined_map(reflection)
    extended_insert_statements!(
        ir.code, ir.codelocs,
        (stmt, ii) -> begin
            # We instrument every new line, before the line.
            # So the first statement of the line is the one we replace
            stmt isa Expr || return nothing
            ii>1 && @inbounds(ir.codelocs[ii]==ir.codelocs[ii-1]) && return nothing
            return 6
        end,
        (stmt, i, orig_i) -> enter_debug_statements(
            ir.slotnames, isdefined_map, reflection.method,
            stmt, i, orig_i
        );
    )
    return ir
end

const handeval_pass = Cassette.@pass instrument!

ext: jl


text:
using FractionalDiffEq, Plots, SpecialFunctions
analytical(x)=exp(1+x)*erfc(sqrt(1+x))
tspan = LinRange(-1, 1, 100)
prob = FIEProblem([1, 1], [0.5, 0], 1, tspan)
sol = solve(prob, 20, SpectralUltraspherical())
solanalytical = analytical.(tspan)
plot(sol, title="Second kind Abel integral equation", label="Numerical")
plot!(tspan, solanalytical, ls=:dash, label="Analytical")


text:
## Reads in data for stackoverflow hypergraph from
# https://www.cs.cornell.edu/~arb/data/stackoverflow-answers/
# and stores it in .mat format with main clusters

include("../src/hypergraph-clustering-utils.jl")
## Read label names
pathname = joinpath(dirname(dirname(@__FILE__)),"data")
dataname = "stackoverflow-answers"
open("$pathname/$dataname/label-names-$dataname.txt") do f
    for line in eachline(f)
        push!(names, line)
    end
end
## Read labels, save in matrix format
using SparseArrays

Labels = Vector{Vector{Int64}}()
pathname = joinpath(dirname(dirname(@__FILE__)),"data")
open("$pathname/$dataname/node-labels-$dataname.txt") do f
    for line in eachline(f)
        labels = parse.(Int64,split(line,','))
        push!(Labels, labels)
    end
end
N = length(names)
LM = elist2incidence(Labels,N)

## Read in the hypergraph
Edges = Vector{Vector{Int64}}()
pathname = joinpath(dirname(dirname(@__FILE__)),"data")
maxN = 0
open("$pathname/$dataname/hyperedges-$dataname.txt") do f
    global maxN
    for line in eachline(f)
        edge = [parse(Int64, v) for v in split(line, ',')]
        sort!(edge)
        if 2 <= length(edge)
            push!(Edges,edge)
            n = maximum(edge)
            if n > maxN
                maxN = n
            end
        end
    end
end

Hmat = elist2incidence(Edges,maxN)
order = round.(Int64,vec(sum(Hmat,dims=2)))
d = vec(sum(Hmat,dims=1))
volA = sum(d)

## Get main clusters, following conditions outlined by Veldt et al. [KDD 2020]
Tconds = Vector{Float64}()
Tsizes = Vector{Int64}()
labels = Vector{Int64}()
num = 1
for lab = 1:length(names)
    global num
    T = findnz(LM[:,lab])[1]
    nT = length(T)
    if nT < 10000 && nT > 2000
        condT, volT, cutT = tl_cond(Hmat,T,d,1.0,volA,order)
        if condT < .2
            println("$num\t  $nT \t $condT \t"*names[lab])
            num+= 1
            push!(Tconds,condT)
            push!(labels,lab)
            push!(Tsizes,nT)
        end
    end
end

LabelMat = LM[:,labels]
LabelNames = names[labels]

## Save in .mat format
matwrite("Stackoverflow_data.mat", Dict("H"=>Hmat,"LabelMat"=>LabelMat,
                                "LabelNames"=>LabelNames,"ClusSizes"=>Tsizes,
                                "ClusConds"=>Tconds))

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
# This file was generated, do not modify it. # hide
scitype(["cat", "mouse", "dog"])

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@testset "sparse arrays" begin
    S = sparse([1],[2],[3])
    G = Bennu.GeneralSparseMatrixCSC(S)
    H = adapt(Array, G)

    for A = (G, H)
        @test size(S) == size(A)
        @test SparseArrays.getcolptr(S) == SparseArrays.getcolptr(A)
        @test rowvals(S) == rowvals(A)
        @test nonzeros(S) == nonzeros(A)
    end
end


text:
module Scarf

export Instance, DP_sS
using Distributions, SpecialFunctions

mutable struct Instance{T <: Real}
    holding_cost::T
    backorder_cost::T
    setup_cost::T
    production_cost::T
    lead_time::Int
    gamma::T
    backlog::Bool
    demand_forecasts::Vector{Normal{T}}
    lt_demand_forecasts::Vector{Normal{T}}
    H::Int
    s::Array{T,1}
    S::Array{T,1}
end

function Instance(h,b,K,c,CV,LT,demands, gamma = 1. ; backlog = one(CV))
    T = typeof(CV)
    dists = Normal{T}[]
    @assert b > c && b > h
    for d in demands
        push!(dists, Normal(d, CV*d))
    end
    lt_dists = Normal{T}[]
    for t in 1:(length(dists)-LT)
        mean_sum = zero(T)
        mean_var = zero(T)
        for i in 0:LT
            mean_sum += mean(dists[t+i])
            mean_var += var(dists[t+i])
        end
        push!(lt_dists, Normal(mean_sum, sqrt(mean_var)))
    end
    backlog = min(one(backlog), max(zero(backlog), backlog))
    S = fill(-Inf, length(dists))
    s = fill(-Inf, length(dists))
    Instance{T}(h,b,K,c,LT, gamma, backlog, dists,lt_dists,length(dists), S, s)
end

mutable struct Pwla{T, S <: StepRangeLen}
    breakpoints::Vector{T}
    range::S
end

function Pwla(s)
    stepsize = Float64(s)
    Pwla(Float64[], stepsize:stepsize:stepsize)
end

Base.push!(p::Pwla, x) = push!(p.breakpoints, x)

function (pwla::Pwla)(y)
    if length(pwla.breakpoints) > 0
        return (@view pwla.breakpoints[end:-1:1])[Int(round((y - pwla.range[1])/step(pwla.range))+1)]
    else
        return zero(y)
    end
end

function production_cost(instance::Instance, q)
    if q > 0
        return instance.production_cost*q + instance.setup_cost
    else
        return zero(q)
    end
end

function L(instance::Instance, y, t::Int)
    h = instance.holding_cost
    b = instance.backorder_cost
    μ = mean(instance.lt_demand_forecasts[t])
    v = var(instance.lt_demand_forecasts[t])
    e = MathConstants.e
    π = MathConstants.pi
    fh = zero(y)
    fb = zero(y)
    if y >= 0
        fh = (((e^(-((y - μ)^2)/2v) -e^(-(μ^2)/2v))*sqrt(v))/(sqrt(2π))) + 1//2*(y-μ)*(erf((y-μ)/(sqrt(2)*sqrt(v))) + erf(μ/(sqrt(2)*sqrt(v))))
        fh *= h
        fb = ((e^(-((y - μ)^2)/2v))*v + sqrt(π/2) * (μ - y) * sqrt(v) * erfc((y-μ)/(sqrt(2)*sqrt(v))))/(sqrt(2π)*sqrt(v))
        fb *= b
    else
        fb = (e^(-(μ^2)/2v)*v + sqrt(π/2) * (μ - y) * sqrt(v) * (1 + erf(μ/(sqrt(2)*sqrt(v)))))/(sqrt(2π)*sqrt(v))
        fb *= b
    end
    return fb + fh
end

function G(instance::Instance, y, t::Int)
    instance.production_cost*y + L(instance, y, t)
end

function C(instance::Instance, x, t::Int, pwla::Pwla)
    S = instance.S[t]
    s = instance.s[t]
    q = x <= s ? S - x : 0.0
    return production_cost(instance, q) + L(instance, x + q, t) + instance.gamma*pwla(x + q)
end

function expected_future_cost(instance::Instance, y, t::Int, pwla::Pwla)
    if t >= instance.H-instance.lead_time
        return zero(y)
    else
        df = instance.demand_forecasts[t]
        ub = instance.backlog ? quantile(df, 0.99999) : y
        ξ = step(pwla.range):step(pwla.range):ub
        x = y .- ξ
        p = cdf.(df, ξ .+ step(ξ)/2) .- cdf.(df, ξ .- step(ξ)/2)
        c(x) = C(instance, x, t+1, pwla)
        return sum(c.(x) .* p) + (1 - cdf(df, y))*c(zero(y))*(1-instance.backlog)
    end
end

function DP_sS(instance::Instance{T}, stepsize::T = one(T)) where T <: Real
    H = instance.H
    λ = instance.lead_time
    maxdemand = max(stepsize, maximum(mean.(instance.demand_forecasts)))
    critical_ratio = 1# (instance.backorder_cost-instance.holding_cost)/instance.backorder_cost
    EOQ = sqrt(2*maxdemand*instance.setup_cost/(critical_ratio*instance.holding_cost))
    ub = 2*(EOQ+maxdemand*λ)
    upperbound = ub + stepsize - ub%stepsize
    C_tplus1 = Pwla(stepsize)
    for t in H-λ:-1:1
        C_t = Pwla(stepsize)
        descending = true
        y = upperbound
        EFC = expected_future_cost(instance, y, t, C_tplus1)
        g = G(instance, y, t) + EFC
        push!(C_t, EFC)
        instance.S[t] = y
        ming = g
        while g <= ming + instance.setup_cost || descending
            y -= stepsize
            EFC = expected_future_cost(instance, y, t, C_tplus1)
            push!(C_t, EFC)
            g = G(instance, y, t) + EFC
            if g < ming
                ming = g
                instance.S[t] = y
            elseif descending
                descending = false
            end
        end
        instance.s[t] = y
        C_t.range = y:stepsize:upperbound
        if C_t(upperbound) < C_t(y)
            #@warn "Upperbound is too low at iteration $t: $(C_t(upperbound)) vs$(C_t(y))"
        end
        C_tplus1 = C_t
    end
end

end #module


text:
#=
    roles

    Tiny[Huge] is a collective value, gathering all finite Float64 values <= [>=] itself.
    Tiny[Huge] is largest inadmissble small [smallest inadmissable large] discrete magnitude.

    constraints

     isfinite(Tiny) & !issubnormal(Tiny) & isfinite(Huge) & !!issubnormal(Tiny)
     ldexp(0.5,-510) < Tiny < ldexp(0.5,-508)
     ldexp(0.5,+509) < Huge < ldexp(0.5,+512)
     0.0 == 1.0-Tiny*Huge == Huge-1.0/Tiny == Tiny-1.0/Huge

     assignments

     Tiny = parse(Float64, "0x1p-510") # ldexp(0.5, -509) == 2.9833_3629_2480_0830e-154
     Huge = parse(Float64, "0x1p+510") # ldexp(0.5, +511) == 3.3519_5198_2485_6493e+153
=#

@inline clrabs{T<:Float65}(x::T) = clearstate(abs(x))

@inline isTiny{T<:AbstractFloat}(x::T) = ((abs(x) <= Tiny(T)) & isnotzero(x))
@inline isHuge{T<:AbstractFloat}(x::T) = ((abs(x) >= Huge(T)) & isnotinf(x))
# projection takes domain(Float64) onto  halfdomain(Float64)
#@inline isProjectedTiny{T<:AbstractFloat}(x::T) = (abs(x) <= TinyProjected(T)) & isnotzero(x)
#@inline isProjectedHuge{T<:AbstractFloat}(x::T) = (abs(x) >= HugeProjected(T)) & isnotinf(x)
# stateful takes projection into one of two coextant, equivalued halfdomain(Float64) bistablities
#@inline isStatefulTiny{T<:AbstractFloat}(x::T)  = (abs(x) <= TinyProjected(T)) & isnotzero(x)
#@inline isStatefulHuge{T<:AbstractFloat}(x::T)  = (abs(x) >= HugeProjected(T)) & isnotinf(x)

@inline isTiny{T<:Float65}(x::T) = (abs(x) <= TinyProjected(Float64)) & isnotzero(x)
@inline isHuge{T<:Float65}(x::T) = (abs(x) >= HugeProjected(Float64)) & isnotinf(x)
# projection takes domain(Float64) onto  halfdomain(Float64)
@inline isProjectedTiny{T<:Float65}(x::T) = isTiny(clrabs(x))
@inline isProjectedHuge{T<:Float65}(x::T) = isHuge(clrabs(x))
# stateful takes projection into one of two coextant, equivalued halfdomain(Float64) bistablities
@inline isStatefulTiny{T<:Float65}(x::T)  = getstate(x) && isTiny(clrabs(x))
@inline isStatefulHuge{T<:Float65}(x::T)  = getstate(x) && isHuge(clrabs(x))

# Bias is used in projective and reflective mappings (projection and reflection behave as dual ops)
Bias(::Type{Float64}) = 511                     # exponent_bias(Float64) >> 1
Bias(::Type{UInt64})  = 0x00000000000001ff      # use with bitops

# Tiny and Huge, assigned
#
Tiny(::Type{Float64}) = 2.983336292480083e-154  # collective value representing all smallest normal magnitudes
Huge(::Type{Float64}) = 3.3519519824856493e153  # collective value representing all largest  normal magnitudes
Tiny(::Type{UInt64})  = 0x2010000000000000      # use with bitops (Tiny as largest  inadmissble small magnitude)
Huge(::Type{UInt64})  = 0x5fd0000000000000      # use with bitops (Huge as smallest inadmissble large magnitude)

# Tiny and Huge, projected
#
TinyProjected(::Type{Float64}) = 4.450147717014403e-308 # ldexp(0.5,-1020) == frexp(pushout(Tiny(Float64))
HugeProjected(::Type{Float64}) = 0.5000000000000000     # ldexp(0.5,+0)    == frexp(pushout(Huge(Float64))
TinyStateful( ::Type{Float64}) = 8.0000000000000000     # ldexp(0.5,+4)    == set_ebit(TinyProjected(Float64))
HugeStateful( ::Type{Float64}) = 8.98846567431158e+307  # ldexp(0.5,+1024) == set_ebit(HugeProjected(Float64))
                                                        #            use with bitops
TinyProjected(::Type{UInt64})  = 0x0020000000000000     # reinterpret(UInt64, TinyProjected(Float64))
HugeProjected(::Type{UInt64})  = 0x3fe0000000000000     # reinterpret(UInt64, HugeProjected(Float64))
TinyStateful( ::Type{UInt64})  = 0x4020000000000000     # reinterpret(UInt64, TinyStateful(Float64))
HugeStateful( ::Type{UInt64})  = 0x7fe0000000000000     # reinterpret(UInt64, HugeStateful(Float64))


text:
abstract type AbstractAmplifier end

struct SimpleAmplifier <: AbstractAmplifier
    name::String
    in_min::typeof(0.0V)
    in_max::typeof(0.0V)
    gain::Float64
    offset::typeof(0.0V)
    current_max::typeof(0.0mA)
end

name(amp::SimpleAmplifier) = amp.name
in_min(amp::SimpleAmplifier) = amp.in_min
out_min(amp::SimpleAmplifier) = amp.in_min*gain + offset
in_max(amp::SimpleAmplifier) = amp.in_max
out_max(amp::SimpleAmplifier) = amp.in_max*gain + offset
in_span(amp::SimpleAmplifier) = in_max(amp) - in_min(amp)
out_span(amp::SimpleAmplifier) = out_max(amp) - out_min(amp)
current_max(amp::SimpleAmplifier) = amp.current_max
gain(amp::SimpleAmplifier) = amp.gain
offset(amp::SimpleAmplifier) = amp.offset

function validate_input(amp::SimpleAmplifier, input::HasVoltageUnits)
    if input < in_min(amp) || input > in_max(amp)
        error("Voltage $command is outside of the acceptable command range of the amplifier")
    end
end

function amplify(amp::SimpleAmplifier, input::HasVoltageUnits)
    validate_input(amp, input)
    return gain(amp) * input + offset(amp)
end

AMPLIFIERS = Dict{String, Any}()

function retrieve_amplifier(amp_name::String)
    if !haskey(AMPLIFIERS, amp_name)
        error("Unknown amplifier")
    end
    return AMPLIFIERS[amp_name]
end

function register_amplifier(amp::T) where {T<:AbstractAmplifier}
    if haskey(AMPLIFIERS, name(amp))
        error("An amplifier named $(name(amp)) is already registered.")
    end
    AMPLIFIERS[name(amp)] = amp
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* File containg models of columns: distillation, stripping, absorbers
#* rectifier, ....
#*
#* The default nomenclature is:
#*		Type_Column_reboilertype_condensertyper
#*
#* where:
#*	Type = refluxed or reboiled or section
#*	Column = Stripping, Absorption, Rectifier, Distillation
#*	Reboiler type (if exists) = kettle or thermosyphon
#*	Condenser type (if exists) = with subccoling or without subcooling
#*
#*-----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: column.mso 511 2008-05-12 17:25:33Z paula $
#*---------------------------------------------------------------------
type Section_Column
	Section_Column()=begin
		new(
			SectionColumnBasic(),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed stream",
				:PosX=>0,
				:PosY=>0.55
			)),
			[
				:(FeedTray.F= TRAYS([1:_base_1.NumberOfTrays]).Inlet.F),
				:(FeedTray.T = TRAYS([1:_base_1.NumberOfTrays]).Inlet.T),
				:(FeedTray.P = TRAYS([1:_base_1.NumberOfTrays]).Inlet.P),
				:(FeedTray.z = TRAYS([1:_base_1.NumberOfTrays]).Inlet.z),
				:(FeedTray.v = TRAYS([1:_base_1.NumberOfTrays]).Inlet.v),
				:(FeedTray.h = TRAYS([1:_base_1.NumberOfTrays]).Inlet.h),
				:(0*"mol/h"= TRAYS([1:_base_1.NumberOfTrays]).Inlet.F),
				:(300*"K" = TRAYS([1:_base_1.NumberOfTrays]).Inlet.T),
				:(1*"atm" = TRAYS([1:_base_1.NumberOfTrays]).Inlet.P),
				:(0.1 = TRAYS([1:_base_1.NumberOfTrays]).Inlet.z),
				:(0 = TRAYS([1:_base_1.NumberOfTrays]).Inlet.v),
				:(0*"J/mol" = TRAYS([1:_base_1.NumberOfTrays]).Inlet.h),
			],
			[
				"FeedTrayTop Inlet Flow","FeedTrayTop Inlet Temperature","FeedTrayTop Inlet Pressure","FeedTrayTop Inlet Composition","FeedTrayTop Inlet Vapour Fraction","FeedTrayTop Inlet Enthalpy","Inlet Tray - Flow Sealed","Inlet Tray -  Temperature","Inlet Tray -  Pressure","Inlet Tray -  Composition","Inlet Tray -  Vapour Fraction","Inlet Tray -  Enthalpy",
			],
			[:FeedTray,]
		)
	end
	_base_1::SectionColumnBasic
	FeedTray::stream
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Section_Column
function set(in::Section_Column)
	_base_1.NumberOfFeeds = 1

end
function setEquationFlow(in::Section_Column)
	#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	# Equating Feed Tray Variables to Trays Variables
	#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	if i == FeedTrayLocation(1)
		addEquation(1)
		addEquation(2)
		addEquation(3)
		addEquation(4)
		addEquation(5)
		addEquation(6)
	else
		addEquation(7)
		addEquation(8)
		addEquation(9)
		addEquation(10)
		addEquation(11)
		addEquation(12)
	end
end
function atributes(in::Section_Column,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/SectionColumn"
	fields[:Brief]="Model of a column section."
	fields[:Info]="== Model of a column section containing ==
* NumberOfTrays TRAYS.
* One Feed Inlet.

== Specify ==
* the feed stream of each tray (Inlet);
* the Murphree eficiency for each tray Emv;
* the InletLiquid stream of the top tray;
* the InletVapour stream of the bottom tray.

== Initial Conditions ==
* the TRAYS temperature (OutletLiquid.T);
* the TRAYS liquid level (Level) OR the TRAYS liquid flow (OutletLiquid.F);
* (NoComps - 1) OutletLiquid (OR OutletVapour) compositions for each tray.
"
	drive!(fields,_)
	return fields
end
Section_Column(_::Dict{Symbol,Any})=begin
	newModel=Section_Column()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# worker mode
const _GW_RUNME_WORKERMODE_FlAG = "-w"
is_worker_mode(gwt::GWTask) = get!(gwt, _GW_RUNME_WORKERMODE_FlAG, false)

function _parse_args!(gwt::GWTask, argv)
    # worker mode
    set!(gwt, _GW_RUNME_WORKERMODE_FlAG, any(argv .== _GW_RUNME_WORKERMODE_FlAG))
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
module LogMath

include("math.jl")
include("LogNum.jl")

end # module


text:
"""
`PDBML <: Format`

Protein Data Bank Markup Language (PDBML), a representation of PDB data in XML format.
"""
immutable PDBML <: Format end

function _get_text(elem, name)
    sub = find_element(elem, name)
    if sub !== nothing
        return(content(sub))
    end
    throw(ErrorException("There is not $name for $elem"))
end

function _get_ins_code(elem)
    sub = find_element(elem, "pdbx_PDB_ins_code")
    if sub !== nothing
        return(content(sub))
    else
        return("")
    end
end

function _get_atom_iterator(document::LightXML.XMLDocument)
    pdbroot = root(document)
    child_elements(get_elements_by_tagname(pdbroot, "atom_siteCategory")[1])
end

"""
`parse(pdbml, ::Type{PDBML}; chain="all", model="all", group="all", atomname="all", onlyheavy=false, label=true, occupancyfilter=false)`

Reads a `LightXML.XMLDocument` representing a pdb file.
Returns a list of `PDBResidue`s (view `MIToS.PDB.PDBResidues`).
Setting `chain`, `model`, `group`, `atomname` and `onlyheavy` values
can be used to select of a subset of all residues. If not set, all residues are returned.
If the keyword argument `label` (default: `true`) is `false`,
the **auth_** attributes will be use instead of the **label_** attributes for `chain`, `atom` and residue `name` fields.
The **auth_** attributes are alternatives provided by an author in order to match the identification/values
used in the publication that describes the structure.
If the keyword argument `occupancyfilter` (default: `false`) is `true`, only the atoms with the best occupancy are returned.
"""
function parse(pdbml::LightXML.XMLDocument, ::Type{PDBML}; chain::ASCIIString = "all",
               model::ASCIIString = "all", group::ASCIIString = "all", atomname::ASCIIString="all",
               onlyheavy::Bool=false, label::Bool=true, occupancyfilter::Bool=false)

    residue_dict = OrderedDict{PDBResidueIdentifier, Vector{PDBAtom}}()

    prefix = label ? "label" : "auth"
    chain_attribute = string(prefix, "_asym_id")
    atom_attribute = string(prefix, "_atom_id")
    comp_attribute = string(prefix, "_comp_id")

    atoms = _get_atom_iterator(pdbml)
    for atom in atoms

        atom_group = _get_text(atom, "group_PDB")
        atom_model = _get_text(atom, "pdbx_PDB_model_num")
        atom_chain = _get_text(atom, chain_attribute)
        atom_name = _get_text(atom, atom_attribute)
        element = _get_text(atom, "type_symbol")

        if  (group=="all" || group==atom_group) && (chain=="all" || chain==atom_chain) &&
                (model=="all" || model==atom_model) && (atomname=="all" || atomname==atom_name) && (!onlyheavy || element!="H")

            PDBe_number = _get_text(atom, "label_seq_id")

            #  Residue_No  _atom_site.auth_seq_id
            #  Ins_Code    _atom_site.pdbx_PDB_ins_code
            PDB_number = string(_get_text(atom, "auth_seq_id"), _get_ins_code(atom))
            name = _get_text(atom, comp_attribute)
            x = float(_get_text(atom, "Cartn_x"))
            y = float(_get_text(atom, "Cartn_y"))
            z = float(_get_text(atom, "Cartn_z"))
            occupancy = float(_get_text(atom, "occupancy"))
            B = _get_text(atom, "B_iso_or_equiv")

            residue_id = PDBResidueIdentifier(PDBe_number, PDB_number, name, atom_group, atom_model, atom_chain)
            atom_data  = PDBAtom(Coordinates(x,y,z), atom_name, element, occupancy, B)

            value = get!(residue_dict, residue_id, PDBAtom[])
            push!(value, atom_data)
        end
    end
    _generate_residues(residue_dict, occupancyfilter)
end

# Download PDB
# ============

function _inputnameforgzip(outfile)
    if endswith(outfile, ".gz")
        return(outfile)
    end
    string(outfile, ".gz")
end

"""
Download a gzipped PDB file from PDB database.
Requires a four character `pdbcode`.
By default the `format` is xml and uses the `baseurl` http://www.rcsb.org/pdb/files/.
`outfile` is the path/name of the output file.
"""
function downloadpdb(pdbcode::AbstractString; format::ASCIIString="xml", outfile::AbstractString="default", baseurl::ASCIIString="http://www.rcsb.org/pdb/files/")
    if length(pdbcode)== 4
        filename = string(uppercase(pdbcode), ".", lowercase(format),".gz")
        outfile = outfile == "default" ? filename : _inputnameforgzip(outfile)
        sepchar = endswith(baseurl,"/") ? "" : "/";
        download(string(baseurl,sepchar,filename) , outfile)
    else
        throw(string(pdbcode, " is not a correct PDB code"))
    end
end

# RESTful PDB interface
# =====================

"""
Access general information about a PDB entry (e.g., Header information) using the RESTful interface of the PDB database (describePDB).
Returns a Dict for the four character `pdbcode`.
"""
function getpdbdescription(pdbcode::ASCIIString)
    if length(pdbcode)== 4
        query = string("http://www.rcsb.org/pdb/rest/describePDB?structureId=", lowercase(pdbcode))
        xmlroot = root(parse_file(download(query)))
        description = get_elements_by_tagname(xmlroot, "PDB")[1]
        return( attributes_dict(description) )
    else
        throw(string(pdbcode, " is not a correct PDB code"))
    end
end


text:
using Compat.Test, Compat.SparseArrays, Nothings, Compat

@test 1==1


text:
using ResourceManagement
using DataFrames
using Test
using CSV

# dflabor = ReadLaborTracker("src\\TEAM_PLANNED_FWD24.csv");
dflabor = ReadLaborTracker("C:\\Users\\junqueg\\OneDrive - dematic.com\\MyDocuments_BackUp\\15. Programming\\Projects\\ResourceManagement.jl\\src\\TEAM_PLANNED_FWD24_NOV.csv");
# dfAvail = ReadAvailHours("src\\UTILREPORT_FWD_NOV.csv");
dfAvail = ReadAvailHours("C:\\Users\\junqueg\\OneDrive - dematic.com\\MyDocuments_BackUp\\15. Programming\\Projects\\ResourceManagement.jl\\src\\UTILREPORT_FWD_NOV.csv");
# dfRev = ReadLaborTracker("src\\TEAM_ACTUALPLAN_REV18.csv");
dfRev = ReadLaborTracker("C:\\Users\\junqueg\\Documents\\My Documents\\15. Programming\\Projects\\ResourceManagement.jl\\src\\TEAM_ACTUALPLAN_REV18.csv");

dfcost = ReadCostTracker("src\\costracker.csv");
dfcostnew =  ReadCostTracker("src\\costracker_0606.csv");
## costs
p149529 = Project("149529", "Walmart");
fetchAndWriteProjectFinances!(dfcost, "149529", "430300", p149529)

p150547 = Project("150547", "Richemont");
fetchAndWriteProjectFinances!(dfcost, "150547", "430300", p150547)

p158070 = Project("158070", "Walmart");
fetchAndWriteProjectFinances!(dfcost, "158070", "430300", p158070)

p152242 = Project("152242", "Amazon");
fetchAndWriteProjectFinances!(dfcost, "152242", "430300", p152242)

p152385 = Project("152385", "Glazers");
fetchAndWriteProjectFinances!(dfcost, "152385", "430300", p152385)

p154662 = Project("154662", "Burlington")
fetchAndWriteProjectFinances!(dfcostnew, "154662", "430300", p154662)

# dfcost2 = CSV.read("src\\costracker.csv", DataFrame)

## Build Program
program = Program();
program = p149529 + p150547
###

A = DisciplineLabor()
B = DisciplineLabor();
C = DisciplineLabor();
A.Budget.Hours = 160.0;
B.Budget.Hours = 200.0;
C.Budget.Hours = 240.0;
T=TeamLabor()
T = A + B + C

Tony = DisciplineLabor("430300", "Higa Anthony");
fetchAndWritePlannedHours!(dflabor, "HIGA ANTHONY", Tony, :fwd);
fetchAndWritePlannedHours!(dfRev, "HIGA ANTHONY", 18, Tony, :rev);
fetchAndWriteRevActualHours!(dfRev, "HIGA ANTHONY", 18, Tony);

Julie = DisciplineLabor("430300", "BRIONES JULITA");
fetchAndWritePlannedHours!(dflabor, "BRIONES JULITA", Julie, :fwd);
fetchAndWritePlannedHours!(dfRev, "BRIONES JULITA", Julie, :rev);
fetchAndWriteRevActualHours!(dfRev, "BRIONES JULITA",Julie);

Brad = DisciplineLabor("430300", "BORDEN BRADLEY")
fetchAndWritePlannedHours!(dflabor, "BORDEN BRADLEY" , Brad, :fwd);
fetchAndWritePlannedHours!(dfRev, "BORDEN BRADLEY", 18, Brad, :rev);
fetchAndWriteRevActualHours!(dfRev, "BORDEN BRADLEY", 18, Brad);

writeAvailableFwdHours!(Tony, dfAvail, 24);
writeAvailableFwdHours!(Julie, dfAvail, 24);
writeAvailableFwdHours!(Brad, dfAvail, 24);

Team1 = Tony + Julie
Team2 = Tony + Julie + Brad

Vp = getFwdPlannedHours(Tony, "");
Vp1 = getFwdPlannedHours(Tony, "153804");
Vp2 = getFwdPlannedHours(Tony, "154558");
Va = Tony.FwdHoursAvailable;
Vpr = getRevPlannedHours(Tony, "");
Vp1r = getRevPlannedHours(Tony, "153804");
Vp2r = getRevPlannedHours(Tony, "154558");
Tactual = getRevActualHours(Tony, "");

TU1 = getFwdUtilization(Tony, "");
TU2 = getFwdUtilization(Tony, "153804");
TU3 = getFwdUtilization(Tony, "154558");

J1 = getFwdPlannedHours(Julie, "");
J1R = getRevPlannedHours(Julie, "");
Jactual = getRevActualHours(Julie, "");

B1 = getFwdPlannedHours(Brad, "");
B1R = getRevPlannedHours(Brad, "");
Bactual = getRevActualHours(Brad, "");

Tcap = getCapacity(Tony);
Jcap = getCapacity(Julie);
Bcap = getCapacity(Brad);

T1cap = getCapacity(Team1);
T2cap = getCapacity(Team2);

TeamDump(Team1)
TeamDump(Team2)

#### Macros
struct MacrosType{T} end
MacrosType(s::String) = MacrosType{Symbol(s)}();

# macro definition
macro prettyPrint(ex)
    return _prettyPrint(ex);
end

function _prettyPrint(ex)
    args = ex.args;
    if ((ex.head == :->) && (length(args))==2)
        cmd = args[2].args[2];
        cmd = MacrosType("$cmd");
        p = eval(args[1])

        ProcessMacros(p, cmd);

    else println("error")
    end

end

function  ProcessMacros(p::DisciplineLabor , cmd::MacrosType{:All})

    println("""
                    Name: $(p.Name)

                    """)
end

function ProcessMacros(::MacrosType{:Sum})
    fwdsum=sum(eval(args[1]).FwdHoursForecast);
    fwdcostsum=sum(eval(args[1]).FwdCostsForecast);
    revsum=sum(eval(args[1]).RevHoursForecast);
    revcostsum=sum(eval(args[1]).RevCostsForecast);

    println("""
                    Name: $(eval(args[1]).Name)
                    Dept: $(eval(args[1]).Dept)
                    Project: $(eval(args[1]).Project)
                    =============================
                    Rate: $(eval(args[1]).Rate)
                    =============================
                    FwdHoursForecast: $fwdsum
                    FwdCostsForecast: $fwdcostsum
                    =============================
                    RevHoursForecast: $revsum
                    RevCostsForecast: $revcostsum
                    =============================
    """)
end

function ProcessMacros(::MacrosType{:Max})
    fwdmax=maximum(eval(args[1]).FwdHoursForecast);
    fwdcostmax=maximum(eval(args[1]).FwdCostsForecast);
    revmax=maximum(eval(args[1]).RevHoursForecast);
    revcostmax=maximum(eval(args[1]).RevCostsForecast);

    println("""
                    Name: $(eval(args[1]).Name)
                    Dept: $(eval(args[1]).Dept)
                    Project: $(eval(args[1]).Project)
                    =============================
                    Rate: $(eval(args[1]).Rate)
                    =============================
                    FwdHoursForecast: $fwdmax
                    FwdCostsForecast: $fwdcostmax
                    =============================
                    RevHoursForecast: $revmax
                    RevCostsForecast: $revcostmax
                    =============================
    """
    )

end

@prettyPrint Tony -> All

@testset "ResourceManagement.jl" begin
    # Write your tests here.

    @test 1 == 1;

    @test Statistics.mean(2.9, 3.1) == 3.0;

    # @test A.Budget.Hours == 160.0;
    # @test B.Budget.Hours == 200.0;
    # @test C.Budget.Hours == 240.0;

    # @test T.Budget.Hours == 600.0;

    @test size(dflabor) == (137, 79);

    @test sum(getFwdPlannedHours(Tony, "")) == 1838;
    @test sum(getFwdPlannedHours(Tony, "153804")) == 240.0
    # @test sum(getRevPlannedHours(Tony, "153804")) == 240.0
    # @test sum(getRevPlannedHours(Tony, "150547")) == 18.0
    @test sum(getFwdPlannedHours(Tony, "154558")) == 200.0
    # @test sum(getRevPlannedHours(Tony, "154558")) == 300.0
    # @test sum(getRevPlannedHours(Tony, "")) == 2233;
    # @test sum(getRevActualHours(Tony, "")) == 2937;

    # @test  sum(Va .*(getFwdUtilization(Tony, "")))== 1446;
    @test sum(Tony.FwdHoursAvailable)[1] == 4016
    @test sum(getFwdAvailableMonthHours(Tony))[1] == 4016
    @test getName(Tony) == "Higa Anthony";
    # @test getProjects(Tony) == ["152242", "153804" , "154558", "154662", "IND-43", "NZ430300"];
    @test sum(getCapacity(Tony)) == 2178;
    @test sum(getCapacity(Tony)) == (sum(getFwdAvailableMonthHours(Tony)) - sum(getFwdPlannedHours(Tony, "")))

    @test sum(getFwdPlannedHours(Julie, "")) == 3078;
    @test sum(getFwdPlannedHours(Julie, "150547")) == 277;
    @test sum(getFwdPlannedHours(Julie, "156257")) == 137;
    @test sum(getFwdAvailableMonthHours(Julie))== 4016
    # @test getProjects(Julie) == ["150547", "152242", "154662", "156257", "NZ430300"]
    # @test sum(getRevPlannedHours(Julie, "")) == 2802;
    # @test sum(getRevActualHours(Julie, "")) == 3216;

    @test sum(getFwdPlannedHours(Brad, "")) == 3821;
    @test sum(getFwdPlannedHours(Brad, "154644")) == 344;
    # @test getProjects(Brad) == ["154644", "154662", "158070"]
    # @test sum(getRevPlannedHours(Brad, "")) == 2727;
    # @test sum(getRevActualHours(Brad, "")) == 3369;

    @test getFwdPlannedHours(Team1, "") == getFwdPlannedHours(Tony, "") + getFwdPlannedHours(Julie, "");
    # @test getRevPlannedHours(Team1, "") == getRevPlannedHours(Tony, "") + getRevPlannedHours(Julie, "");
    @test getFwdPlannedHours(Team2, "") == getFwdPlannedHours(Tony, "") + getFwdPlannedHours(Julie, "") + getFwdPlannedHours(Brad, "");
    # @test getRevPlannedHours(Team2, "") == getRevPlannedHours(Tony, "") + getRevPlannedHours(Julie, "") + getRevPlannedHours(Brad, "");

    @test getFwdAvailableMonthHours(Team1) == getFwdAvailableMonthHours(Tony) + getFwdAvailableMonthHours(Julie);

    # @test getRevActualHours(Team1, "") == getRevActualHours(Tony, "") + getRevActualHours(Julie, "");
    # @test getRevActualHours(Team2, "") == getRevActualHours(Tony, "") + getRevActualHours(Julie, "") + getRevActualHours(Brad, "");

end


text:
using DataStructures

mutable struct Node
    ip::Int
    neighbours::Array{Tuple{Node,Int}}
    risk::Union{Float64, Nothing}
    predecessor::Union{Node, Nothing}
    probability::Float64
end

function initialize_single_source!(graph, start)
	for n in graph
		n.risk = typemax(Float64)
	end
	start.risk = 0
end

function relax!(from_node, to_node, cost)
	risk = from_node.risk + (cost / to_node.probability)
	if to_node.risk > risk
		to_node.risk = risk
		to_node.predecessor = from_node
	end
end

function find_min!(ls::Array{Node, 1})::Node
	min = ls[1]
	index = 1
	for i in 1 : length(ls)
		if ls[i].risk < min.risk
			min = ls[i]
			index = i
		end
	end
	deleteat!(ls, i)
	return min
end

function find(ls::Array{Node, 1}, node)::Int
	for i in 1 : length(ls)
		if ls[i] == node
			return i
		end
	end
	return 0
end

function dijkstra!(graph, start)
	initialize_single_source!(graph, start)
	que::Array{Node, 1} = []
	for i in graph
		push!(que, i)
	end
	while length(que) > 0
		u = find_min!(que)
		for j in u.neighbours
			relax!(u, j[1], j[2])
		end
	end
end

function bellman_ford!(graph, start)
	initialize_single_source!(graph, start)

	for k in 1 : length(graph) - 1
		for i in graph
			for j in i.neighbours
				relax!(i, j[1], j[2])
			end
		end
	end
	for i in graph
		for j in i.neighbours
			if j[1].risk > i.risk + j[2]
				return false
			end
		end
	end
	return true
end


text:
include("fake_imgs.jl")
include("slicetiming.jl")


text:
# LSP indexing solution
# https://github.com/julia-vscode/julia-vscode/issues/800#issuecomment-650085983
if isdefined(@__MODULE__, :LanguageServer)
    include("../src/SpinDoctor.jl")
    using .SpinDoctor
end

using SpinDoctor
using LinearAlgebra
using QuadGK
using Test
using Aqua

include("testutils.jl")

# Run through typical workflow
@testset "Workflow" begin
    include("workflow.jl")
end

@testset "Aqua" begin
    Aqua.test_all(
        SpinDoctor;
        ambiguities = false,
        project_toml_formatting = false, # https://github.com/JuliaTesting/Aqua.jl/issues/72
    )
end

#= TODO: Restructure tests
# Datatypes
@testset "Setups" begin
    include("datatypes/sequences.jl")
    include("datatypes/experiment.jl")
    include("datatypes/femesh.jl")
    include("datatypes/gradient.jl")
    include("datatypes/model.jl")
    include("datatypes/create_model.jl")
    include("datatypes/setup.jl")
    include("datatypes/get_coefficients.jl")
end

# Geometry
@testset "Geometry" begin
    include("geometry/call_tetgen.jl")
    include("geometry/convexhull.jl")
    include("geometry/create_cells.jl")
    include("geometry/create_directions.jl")
    include("geometry/create_fibonacci_sphere.jl")
    include("geometry/create_geometry.jl")
    include("geometry/create_surfaces_cylinder.jl")
    include("geometry/create_surfaces_neuron.jl")
    include("geometry/create_surfaces_sphere.jl")
    include("geometry/deform_domain.jl")
    include("geometry/get_volumes.jl")
    include("geometry/get_mesh_surface.jl")
    include("geometry/get_mesh_surfacenormals.jl")
    include("geometry/gmesh2fem.jl")
    include("geometry/read_cells.jl")
    include("geometry/read_surfaces.jl")
    include("geometry/read_tetgen.jl")
    include("geometry/save_cells.jl")
    include("geometry/save_surfaces.jl")
    include("geometry/save_tetgen.jl")
    include("geometry/split_mesh.jl")
end

# Matrix assembly
@testset "Matrix assembly" begin
    include("matrix_assembly/assemble_mass_matrix.jl")
    include("matrix_assembly/assemble_stiffness_matrix.jl")
    include("matrix_assembly/assemble_flux_matrices.jl")
    include("matrix_assembly/assemble_flux_matrix.jl")
    include("matrix_assembly/couple_flux_matrix.jl")
end

# Matrix formalism
@testset "Matrix formalism" begin
    include("matrix_formalism/eig2length.jl")
    include("matrix_formalism/compute_laplace_eig.jl")
    include("matrix_formalism/solve_mf.jl")
end

# Analytical
@testset "Analytical" begin
    include("analytical/alpha_func.jl")
    include("analytical/compute_bc.jl")
    include("analytical/compute_beta.jl")
    include("analytical/compute_int_I.jl")
    include("analytical/compute_int_J.jl")
    include("analytical/compute_int_K.jl")
    include("analytical/compute_JY.jl")
    include("analytical/compute_nu.jl")
    include("analytical/compute_v.jl")
    include("analytical/find_alpha.jl")
    include("analytical/solve_analytical.jl")
end

# BTPDE
@testset "BTPDE" begin
    include("btpde/solve_btpde.jl")
    include("btpde/solve_btpde_midpoint.jl")
    include("btpde/solve_karger.jl")
end

# ADC
@testset "ADC" begin
    include("adc/solve_hadc.jl")
end

# Postprocessing
@testset "Postprocess" begin
    include("postprocess/fit_adc.jl")
    include("postprocess/fit_tensors.jl")
    include("postprocess/savefield.jl")
end

# Plot
@testset "Plot" begin
    include("plot/plot_mesh.jl")
end
=#


text:
include("../utils.jl")

using Manifolds: induced_basis
using ManifoldsBase: TFVector

@testset "Sphere" begin
    M = Sphere(2)
    @testset "Sphere Basics" begin
        @test repr(M) == "Sphere(2, ℝ)"
        @test typeof(get_embedding(M)) === Euclidean{Tuple{3},ℝ}
        @test representation_size(M) == (3,)
        @test injectivity_radius(M) == π
        @test injectivity_radius(M, ExponentialRetraction()) == π
        @test injectivity_radius(M, ProjectionRetraction()) == π / 2
        @test base_manifold(M) === M
        @test !is_point(M, [1.0, 0.0, 0.0, 0.0])
        @test !is_vector(M, [1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0])
        @test_throws DomainError is_point(M, [2.0, 0.0, 0.0], true)
        @test !is_point(M, [2.0, 0.0, 0.0])
        @test !is_vector(M, [1.0, 0.0, 0.0], [1.0, 0.0, 0.0])
        @test_throws DomainError is_vector(M, [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], true)
        @test injectivity_radius(M, [1.0, 0.0, 0.0], ProjectionRetraction()) == π / 2
    end
    types = [Vector{Float64}]
    TEST_FLOAT32 && push!(types, Vector{Float32})
    TEST_STATIC_SIZED && push!(types, MVector{3,Float64})

    basis_types = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))
    test_atlases = (Manifolds.StereographicAtlas(), Manifolds.RetractionAtlas())
    for T in types
        @testset "Type $T" begin
            pts = [
                convert(T, [1.0, 0.0, 0.0]),
                convert(T, [0.0, 1.0, 0.0]),
                convert(T, [0.0, 0.0, 1.0]),
            ]
            test_manifold(
                M,
                pts,
                test_reverse_diff=isa(T, Vector),
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                test_mutating_rand=isa(T, Vector),
                point_distributions=[Manifolds.uniform_distribution(M, pts[1])],
                tvector_distributions=[
                    Manifolds.normal_tvector_distribution(M, pts[1], 1.0),
                ],
                basis_types_vecs=(DiagonalizingOrthonormalBasis([0.0, 1.0, 2.0]),),
                basis_types_to_from=basis_types,
                test_vee_hat=false,
                retraction_methods=[ProjectionRetraction(), ExponentialRetraction()],
                inverse_retraction_methods=[ProjectionInverseRetraction()],
                is_tangent_atol_multiplier=1,
                test_atlases=test_atlases,
                test_inplace=true,
            )
            @test isapprox(-pts[1], exp(M, pts[1], log(M, pts[1], -pts[1])))
        end
    end

    @testset "Distribution tests" begin
        usd_mvector = Manifolds.uniform_distribution(M, @MVector [1.0, 0.0, 0.0])
        @test isa(rand(usd_mvector), MVector)

        gtsd_mvector =
            Manifolds.normal_tvector_distribution(M, (@MVector [1.0, 0.0, 0.0]), 1.0)
        @test isa(rand(gtsd_mvector), MVector)
    end

    @testset "Embedding test" begin
        p = [1.0, 0.0, 0.0]
        X = [0.0, 1.0, 1.0]
        @test embed(M, p) == p
        q = similar(p)
        embed!(M, q, p)
        @test q == p
        @test embed(M, p, X) == X
        Y = similar(X)
        embed!(M, Y, p, X)
        @test Y == X
    end

    @testset "DefaultOrthonormalBasis edge cases" begin
        B = DefaultOrthonormalBasis(ℝ)
        n = manifold_dimension(M)
        p = [1; zeros(n)]
        for i in 1:n  # p'x ± 1
            vcoord = [j == i for j in 1:n]
            v = [0; vcoord]
            @test get_coordinates(M, p, v, B) ≈ vcoord
            @test get_vector(M, p, vcoord, B) ≈ v
            @test get_coordinates(M, -p, v, B) ≈ vcoord
            @test get_vector(M, -p, vcoord, B) ≈ v
        end
        p = [0; 1; zeros(n - 1)] # p'x = 0
        for i in 1:n
            vcoord = [j == i for j in 1:n]
            v = get_vector(M, p, vcoord, B)
            @test is_vector(M, p, v)
            @test get_coordinates(M, p, v, B) ≈ vcoord
        end
    end

    @testset "log edge case" begin
        n = manifold_dimension(M)
        x = normalize(randn(n + 1))
        v = log(M, x, -x)
        @test norm(v) ≈ π
        @test isapprox(dot(x, v), 0; atol=1e-12)
        vexp = normalize(project(M, x, [1, zeros(n)...]))
        @test v ≈ π * vexp

        x = [1, zeros(n)...]
        v = log(M, x, -x)
        @test norm(v) ≈ π
        @test isapprox(dot(x, v), 0; atol=1e-12)
        vexp = normalize(project(M, x, [0, 1, zeros(n - 1)...]))
        @test v ≈ π * vexp
    end

    @testset "Complex Sphere" begin
        M = Sphere(2, ℂ)
        @test repr(M) == "Sphere(2, ℂ)"
        @test typeof(get_embedding(M)) === Euclidean{Tuple{3},ℂ}
        @test representation_size(M) == (3,)
        p = [1.0, 1.0im, 1.0]
        q = project(M, p)
        @test is_point(M, q)
        Y = [2.0, 1.0im, 20.0]
        X = project(M, q, Y)
        @test is_vector(M, q, X, true; atol=10^(-14))
    end

    @testset "Quaternion Sphere" begin
        M = Sphere(2, ℍ)
        @test repr(M) == "Sphere(2, ℍ)"
        @test typeof(get_embedding(M)) === Euclidean{Tuple{3},ℍ}
        @test representation_size(M) == (3,)
        p = [Quaternion(1.0), Quaternion(1.0im), Quaternion(0.0, 0.0, -1.0, 0.0)]
        q = project(M, p)
        @test is_point(M, q)
        Y = [Quaternion(2.0), Quaternion(1.0im), Quaternion(0.0, 0.0, 20.0, 0.0)]
        X = project(M, q, Y)
        @test is_vector(M, q, X, true; atol=10^(-14))
    end

    @testset "Array Sphere" begin
        M = ArraySphere(2, 2; field=ℝ)
        @test repr(M) == "ArraySphere(2, 2; field = ℝ)"
        @test typeof(get_embedding(M)) === Euclidean{Tuple{2,2},ℝ}
        @test representation_size(M) == (2, 2)
        p = ones(2, 2)
        q = project(M, p)
        @test is_point(M, q)
        Y = [1.0 0.0; 0.0 1.1]
        X = project(M, q, Y)
        @test is_vector(M, q, X)
        M = ArraySphere(2, 2; field=ℂ)

        @test repr(M) == "ArraySphere(2, 2; field = ℂ)"
        @test typeof(get_embedding(M)) === Euclidean{Tuple{2,2},ℂ}
        @test representation_size(M) == (2, 2)
    end

    @testset "StereographicAtlas" begin
        M = Sphere(2)
        A = Manifolds.StereographicAtlas()
        p = randn(3)
        p ./= norm(p)
        for k in [1, -1]
            p *= k
            i = Manifolds.get_chart_index(M, A, p)
            @test i === (p[1] < 0 ? :south : :north)
            a = get_parameters(M, A, i, p)
            q = get_point(M, A, i, a)
            @test isapprox(M, p, q)

            p2 = randn(3)
            p2 ./= norm(p2)
            p3 = randn(3)
            p3 ./= norm(p3)
            X2 = log(M, p, p2)
            X3 = log(M, p, p3)
            B = induced_basis(M, A, i, TangentSpace)

            X2B = get_coordinates(M, p, X2, B)
            X3B = get_coordinates(M, p, X3, B)

            @test inner(M, p, X2, X3) ≈ dot(X2B, local_metric(M, p, B) * X3B)

            X2back = get_vector(M, p, X2B, B)
            X3back = get_vector(M, p, X3B, B)
            @test isapprox(M, p, X2, X2back)
            @test isapprox(M, p, X3, X3back)

            @testset "transition_map" begin
                other_chart = i === :south ? :north : :south
                a_other = Manifolds.transition_map(M, A, i, other_chart, a)
                @test isapprox(M, p, get_point(M, A, other_chart, a_other))

                a_other2 = allocate(a_other)
                Manifolds.transition_map!(M, a_other2, A, i, A, other_chart, a)
                @test isapprox(M, p, get_point(M, A, other_chart, a_other2))

                a_other2 = allocate(a_other)
                Manifolds.transition_map!(M, a_other2, A, i, other_chart, a)
                @test isapprox(M, p, get_point(M, A, other_chart, a_other2))
            end
        end
    end

    @testset "Metric conversion is the identity" begin
        p = [1.0, 0.0, 0.0]
        X = [0.0, 1.0, 1.0]
        Y = change_representer(M, EuclideanMetric(), p, X)
        @test Y == X
        Z = change_metric(M, EuclideanMetric(), p, X)
        @test Z == X
        @test local_metric(M, p, DefaultOrthonormalBasis()) == Diagonal([1.0, 1.0])
    end

    @testset "Distributions" begin
        sphere = Sphere(2)
        haar_measure = uniform_distribution(sphere)
        pts = rand(haar_measure, 5)
        @test all(p -> is_point(sphere, p), pts)
    end
end


text:
#Подсчитать число прямоугольных перегородок (отрезков нет)

function ugol(r,side1,side2) #функция ведет робота в угол side1 side2
    while !isborder(r,side1) || !isborder(r,side2)
        if !isborder(r,side1)
            move!(r,side1)
        else
            move!(r,side2)
        end
    end
end

function back(r, side)#функция возвращает робота в начало строки в side
    c = 1
    while c != 0
        if !isborder(r, side)
            move!(r, side)
        else
            c = obhod(r,side)
            if c == 0
                break
            end
        end
    end
end

function reverse(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+2, 4))
        end
    end
end

function next_side(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+3,4))
        end
    end
end

function obhod(r,side)#функция обходит перегородки чтобы идти по прямой
    k = 0
    while (isborder(r,side) && !isborder(r,next_side(side)))
        move!(r,next_side(side))
        k += 1
    end
    if isborder(r,next_side(side))&&isborder(r,side)
        for i in 1:k
            move!(r,reverse(next_side(side)))
        end
        return 0
    end
    if !isborder(r,side)
        move!(r,side)
    end
    if k != 0
        while isborder(r,reverse(next_side(side)))
            move!(r,side)
        end
        for i in 1:k
            move!(r,reverse(next_side(side)))
        end
    end
    return 1
end

function num_pr(r)
    ugol(r, West, Sud)
    n = 0
    flag = 1
    while !isborder(r, Nord)
        c = 1
        while c != 0
            if !isborder(r, Ost)
                move!(r, Ost)
                if(flag == 0 && isborder(r, Nord))
                    flag = 1
                    n += 1
                elseif(flag == 1 && !isborder(r, Nord))
                    flag = 0
                end
            else
                c = obhod(r,Ost)
                if c == 0
                    break
                end
            end
        end
        move!(r, Nord)
        back(r,West)
    end
print(n)
end

max_stars_repo_name: Aminulik/mirea-progs

max_issues_repo_name: Aminulik/mirea-progs

max_forks_repo_name: Aminulik/mirea-progs


text:
@testset "optimize!" begin
    fm = _build_dmip_formulation()
    result = @inferred Cerberus.optimize!(fm, CONFIG)
    @test result.primal_bound ≈ 0.1 / 2.1
    @test result.dual_bound ≈ 0.1 / 2.1
    @test length(result.best_solution) == 3
    @test result.best_solution ≈ [1.0, 2 / 2.1, 0.0]
    @test result.termination_status == Cerberus.OPTIMAL
    @test result.total_node_count == 3
    @test result.total_simplex_iters == 0
end

@testset "process_node!" begin
    # A feasible model
    let
        fm = _build_dmip_formulation()
        state = _CurrentState()
        node = Cerberus.Node()
        result = @inferred Cerberus.process_node!(state, fm, node, CONFIG)
        @test result.status == Cerberus.OPTIMAL_LP
        @test result.cost ≈ 0.5 - 2.5 / 2.1
        @test length(result.x) == 3
        @test result.x ≈ [0.5, 2.5 / 2.1, 0.0]
        @test result.simplex_iters == 0
        @test result.depth == 0
        result.int_infeas == 0

        true_basis = @inferred Cerberus.get_basis(state)
        _test_is_equal_to_dmip_basis(true_basis)
    end

    # A feasible model with no incrementalism
    let
        fm = _build_dmip_formulation()
        no_inc_config = Cerberus.AlgorithmConfig(
            warm_start_strategy = Cerberus.NO_WARM_STARTS,
            model_reuse_strategy = Cerberus.NO_MODEL_REUSE,
        )
        state = _CurrentState()
        node = Cerberus.Node()
        result =
            @inferred Cerberus.process_node!(state, fm, node, no_inc_config)

        @test result.status == Cerberus.OPTIMAL_LP
        @test result.cost ≈ 0.5 - 2.5 / 2.1
        @test length(result.x) == 3
        @test result.x ≈ [0.5, 2.5 / 2.1, 0.0]
        @test result.simplex_iters == 0
        @test result.depth == 0
        result.int_infeas == 0
    end

    # An infeasible model
    let
        fm = _build_dmip_formulation()
        state = _CurrentState()
        # A bit hacky, but force infeasibility by branching both up and down.
        node = Cerberus.Node(
            [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))],
            [Cerberus.BoundUpdate(_CVI(1), _GT(1.0))],
            Cerberus.AffineConstraint{_LT}[],
            Cerberus.AffineConstraint{_GT}[],
            2,
        )
        result = @inferred Cerberus.process_node!(state, fm, node, CONFIG)
        @test result.status == Cerberus.INFEASIBLE_LP
        @test result.cost == Inf
        @test !isdefined(result, :x)
        @test result.simplex_iters == 0
        @test result.depth == 2
        @test result.int_infeas == 0
    end

    @testset "prune by parent bound" begin
        fm = _build_dmip_formulation()
        node = Cerberus.Node()
        node.dual_bound = 1.0
        state = _CurrentState(primal_bound = 0.0)
        result = @inferred Cerberus.process_node!(state, fm, node, CONFIG)
        @test result.status == Cerberus.PRUNED_BY_PARENT_BOUND
        @test isnan(result.cost)
        @test !isdefined(result, :x)
        @test result.simplex_iters == 0
        @test result.depth == 0
        result.int_infeas == 0
    end
end

@testset "_num_int_infeasible" begin
    config = Cerberus.AlgorithmConfig()
    let form = _build_dmip_formulation()
        state = Cerberus.CurrentState()
        let x_int = [1.0, 3.2, 0.0]
            state.current_solution = x_int
            @test Cerberus._num_int_infeasible(state, form, CONFIG) == 0
        end
        let x_int = [1.0 - 0.9CONFIG.int_tol, 3.2, 0.0 + 0.9CONFIG.int_tol]
            state.current_solution = x_int
            @test Cerberus._num_int_infeasible(state, form, CONFIG) == 0
        end
        let x_int = [1.0 - 2CONFIG.int_tol, 3.2, 0.0 + 1.1CONFIG.int_tol]
            state.current_solution = x_int
            @test Cerberus._num_int_infeasible(state, form, CONFIG) == 2
        end
    end

    let form = _build_gi_dmip_formulation()
        state = Cerberus.CurrentState()
        x_int = [0.6, 0.0, 2.0]
        state.current_solution = x_int
        @test Cerberus._num_int_infeasible(state, form, CONFIG) == 0
        x_int[2] = 0.9
        @test Cerberus._num_int_infeasible(state, form, CONFIG) == 1
        x_int[2] = 1.0
        x_int[3] = 3.0
        @test Cerberus._num_int_infeasible(state, form, CONFIG) == 0
        x_int[2] = 1.0
        x_int[3] = 2.9
        @test Cerberus._num_int_infeasible(state, form, CONFIG) == 1
    end
end

@testset "_store_basis_if_desired!" begin
    fm = _build_dmip_formulation()
    state = Cerberus.CurrentState()
    node = Cerberus.Node()
    Cerberus.process_node!(state, fm, node, CONFIG)
    fc = Cerberus.Node(
        Cerberus.BoundUpdate{_LT}[],
        [Cerberus.BoundUpdate(_CVI(1), _GT(1.0))],
        Cerberus.AffineConstraint{_LT}[],
        Cerberus.AffineConstraint{_GT}[],
        2,
    )
    oc = Cerberus.Node(
        [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))],
        Cerberus.BoundUpdate{_GT}[],
        Cerberus.AffineConstraint{_LT}[],
        Cerberus.AffineConstraint{_GT}[],
        2,
    )
    cost = 12.3
    result = Cerberus.NodeResult(
        Cerberus.OPTIMAL_LP,
        cost,
        [1.2, 2.3, 3.4],
        1492,
        12,
        13,
    )
    let no_incrementalism_config = Cerberus.AlgorithmConfig(
            warm_start_strategy = Cerberus.NO_WARM_STARTS,
            model_reuse_strategy = Cerberus.NO_MODEL_REUSE,
        )
        @inferred Cerberus._store_basis_if_desired!(
            state,
            (oc, fc),
            no_incrementalism_config,
        )
        @test isempty(state.warm_starts)
    end
    let warm_start_only_config = Cerberus.AlgorithmConfig(
            warm_start_strategy = Cerberus.WARM_START_WHENEVER_POSSIBLE,
            model_reuse_strategy = Cerberus.NO_MODEL_REUSE,
        )
        @inferred Cerberus._store_basis_if_desired!(
            state,
            (oc, fc),
            warm_start_only_config,
        )
        @test length(state.warm_starts) == 2
        @test haskey(state.warm_starts, fc)
        fc_basis = state.warm_starts[fc]
        _test_is_equal_to_dmip_basis(fc_basis)
        @test haskey(state.warm_starts, oc)
        oc_basis = state.warm_starts[oc]
        _test_is_equal_to_dmip_basis(oc_basis)
        # TODO: Test more than one child
    end
    empty!(state.warm_starts)
    let warm_start_off_dives_config = Cerberus.AlgorithmConfig(
            warm_start_strategy = Cerberus.WARM_START_WHEN_BACKTRACKING,
            model_reuse_strategy = Cerberus.REUSE_MODEL_ON_DIVES,
        )
        @inferred Cerberus._store_basis_if_desired!(
            state,
            (oc, fc),
            warm_start_off_dives_config,
        )
        @test length(state.warm_starts) == 1
        @test haskey(state.warm_starts, oc)
        oc_basis = state.warm_starts[oc]
        _test_is_equal_to_dmip_basis(oc_basis)
    end
end

@testset "update_state!" begin
    fm = _build_dmip_formulation()
    starting_pb = 12.3
    simplex_iters_per = 18
    depth = 7
    state = _CurrentState(primal_bound = starting_pb)
    @test Cerberus._is_root_node(Cerberus.pop_node!(state.tree))
    node = Cerberus.Node()

    @test isempty(state.warm_starts)

    # 1. Prune by infeasibility
    let nr = Cerberus.NodeResult(node)
        nr.status = Cerberus.INFEASIBLE_LP
        nr.cost = Inf
        nr.simplex_iters = simplex_iters_per
        nr.int_infeas = 0
        @inferred Cerberus.update_state!(state, fm, node, nr, CONFIG)
        @test isempty(state.tree)
        @test state.total_node_count == 1
        @test state.primal_bound == starting_pb
        @test state.dual_bound == -Inf
        @test length(state.best_solution) == 0
        @test all(isnan, values(state.best_solution))
        @test state.total_simplex_iters == simplex_iters_per
        @test isempty(state.warm_starts)
    end
    state.on_a_dive = true
    state.rebuild_model = true

    # 2. Prune by bound
    let nr = Cerberus.NodeResult(node)
        nr.status = Cerberus.OPTIMAL_LP
        nr.cost = 13.5
        frac_soln = [0.2, 3.4, 0.6]
        nr.x = frac_soln
        state.current_solution = nr.x
        nr.simplex_iters = simplex_iters_per
        nr.depth = depth
        nr.int_infeas = Cerberus._num_int_infeasible(state, fm, CONFIG) == 2
        @inferred Cerberus.update_state!(state, fm, node, nr, CONFIG)
        @test isempty(state.tree)
        @test state.total_node_count == 2
        @test state.primal_bound == starting_pb
        @test state.dual_bound == -Inf
        @test length(state.best_solution) == 0
        @test all(isnan, values(state.best_solution))
        @test state.total_simplex_iters == 2 * simplex_iters_per
        @test isempty(state.warm_starts)
    end
    state.on_a_dive = true
    state.rebuild_model = true

    # 3. Prune by integrality
    new_pb = 11.1
    int_soln = [1.0, 3.4, 0.0]
    let nr = Cerberus.NodeResult(node)
        nr.status = Cerberus.OPTIMAL_LP
        nr.cost = new_pb
        nr.x = copy(int_soln)
        state.current_solution = nr.x
        nr.simplex_iters = simplex_iters_per
        nr.depth = depth
        nr.int_infeas = Cerberus._num_int_infeasible(state, fm, CONFIG)
        @test nr.int_infeas == 0
        @inferred Cerberus.update_state!(state, fm, node, nr, CONFIG)
        @test isempty(state.tree)
        @test state.total_node_count == 3
        @test state.primal_bound == new_pb
        @test state.dual_bound == -Inf
        @test state.best_solution == int_soln
        @test state.total_simplex_iters == 3 * simplex_iters_per
        @test isempty(state.warm_starts)
    end
    state.on_a_dive = true
    state.rebuild_model = true

    # 4. Branch
    db = 10.1
    let nr = Cerberus.process_node!(state, fm, node, CONFIG)
        frac_soln_2 = [0.0, 2.9, 0.6]
        nr.cost = db
        nr.x = frac_soln_2
        state.current_solution = nr.x
        nr.simplex_iters = simplex_iters_per
        nr.depth = depth
        nr.int_infeas = Cerberus._num_int_infeasible(state, fm, CONFIG)
        @test nr.int_infeas == 1
        @inferred Cerberus.update_state!(state, fm, node, nr, CONFIG)
        @test Cerberus.num_open_nodes(state.tree) == 2
        @test state.total_node_count == 4
        @test state.primal_bound == new_pb
        @test state.dual_bound == -Inf
        @test state.best_solution == int_soln
        @test state.total_simplex_iters == 4 * simplex_iters_per
        @inferred Cerberus.update_dual_bound!(state)
        @test state.dual_bound == db
        let fc = Cerberus.pop_node!(state.tree)
            @test fc.lt_bounds == [Cerberus.BoundUpdate(_CVI(3), _LT(0.0))]
            @test isempty(fc.gt_bounds)
            @test fc.dual_bound == db
            @test length(state.warm_starts) == 1
            @test !haskey(state.warm_starts, fc)
        end
        let oc = Cerberus.pop_node!(state.tree)
            @test isempty(oc.lt_bounds)
            @test oc.gt_bounds == [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))]
            @test oc.dual_bound == db
            @test haskey(state.warm_starts, oc)
            oc_basis = state.warm_starts[oc]
            _test_is_equal_to_dmip_basis(oc_basis)
        end
    end

    # 5. Prune by parent bound
    let nr = Cerberus.NodeResult(node)
        nr.status = Cerberus.PRUNED_BY_PARENT_BOUND
        nr.depth = depth
        @inferred Cerberus.update_state!(state, fm, node, nr, CONFIG)
        @test isempty(state.tree)
        @test state.total_node_count == 5
        @test state.primal_bound == new_pb
        @test state.dual_bound == db
        @test state.best_solution == int_soln
        @test state.total_simplex_iters == 4 * simplex_iters_per
        @test length(state.warm_starts) == 1
    end
end


text:
module TestMerge

include("preamble.jl")
using BangBang.NoBang: SingletonDict
using InitialValues: Init, asmonoid

function signature(args)
    #! format: off
    strs = (
        if a isa Function
            string(a)
        elseif a isa NamedTuple
            "::NamedTuple"
        else
            "::$(typeof(a))"
        end for a in args
    )
    #! format: on
    return join(strs, ", ")
end

@testset "widen" begin
    @testset "merge!!($(signature(args)))" for (args, desired) in Any[
        ((Dict(:a => 1), Dict(:b => 0.5)), Dict(:a => 1.0, :b => 0.5)),
        ((Dict(:a => 1), Dict("b" => 2)), Dict(:a => 1, "b" => 2)),
        ((Dict(:a => 1), Dict("b" => 0.5)), Dict(:a => 1.0, "b" => 0.5)),
        (((a = 1,), (b = 2,)), (a = 1, b = 2)),
        (((a = 1,), Dict(:b => 2)), (a = 1, b = 2)),
        (((a = 1,), Dict(:b => 2), (b = 3,)), (a = 1, b = 3)),
        ((Empty(Dict), Dict(:a => 1)), Dict(:a => 1)),
        ((Empty(Dict), SingletonDict(:a => 1)), Dict(:a => 1)),
        ((SingletonDict(:a => 1), Dict(:b => 2)), Dict(:a => 1, :b => 2)),
        ((SingletonDict(:a => 1), Dict(:a => 2)), Dict(:a => 2)),
    ]
        d = deepcopy(args[1])
        @test merge!!(args...) ==ₜ desired
        @test args[1] == d
    end
    @testset "merge!!($(signature(args)))" for (args, desired) in Any[
        ((+, Dict(:a => 1), Dict(:a => 0.5)), Dict(:a => 1.5)),
        ((+, (a = 1,), Dict(:a => 0.5)), (a = 1.5,)),
    ]
        d = deepcopy(args[2])
        @test merge!!(args...) ==ₜ desired
        @test args[2] == d
    end
end

@testset "mutate" begin
    @testset "merge!!($(signature(args)))" for (args, desired) in Any[
        ((Dict(:a => 1), Dict(:b => 2)), Dict(:a => 1, :b => 2)),
        ((Dict(:a => 1), (b = 2,)), Dict(:a => 1, :b => 2)),
    ]
        @test merge!!(args...) === args[1] == desired
    end
    @testset "merge!!($(signature(args)))" for (args, desired) in Any[
        ((+, Dict(:a => 1), Dict(:a => 2)), Dict(:a => 3)),
        ((+, Dict(:a => 1), (a = 2,)), Dict(:a => 3)),
        ((+, Dict(:a => 1), (a = 2,), Dict(:a => 3)), Dict(:a => 6)),
    ]
        @test merge!!(args...) === args[2] == desired
    end
end

@testset "Init" begin
    for m!! in [merge!!, mergewith!!(+)]
        d = Dict(:a => 1)
        @test m!!(d, Init(m!!)) === d
        @test m!!(Init(m!!), d) !== d
        @test m!!(Init(m!!), d) ==ₜ d
        @test m!!(Init(m!!), (x for x in d if true)) ==ₜ d
        @test m!!(Init(m!!), Init(m!!)) === Init(m!!)
        @test asmonoid(m!!) === m!!
    end
end

end  # module


text:
function target(directory::AbstractString,
                target_descr::Matrix)
    open("$directory/target", "w") do file
        for i in 1:size(target_descr,1)
            if target_descr[i,1] != ""
                write(file, string(target_descr[i,1]))
                if target_descr[i,2] != ""
                    write(file, " = ", string(target_descr[i,2]))
                    if target_descr[i,3] != ""
                        write(file, " ! ", string(target_descr[i,3]))
                    end
                end
            else
                write(file, "--------")
            end
            write(file, '\n')
        end
    end
end

separator(file) = write(file, "--------\n")
property(file, p, v, c) = write(file, "$p = $v ! $c\n")

function target(directory::AbstractString,
                heading::AbstractString,
                Z::Real, nelc::Integer,
                partial_waves::Vector{Term},
                targets::BSR.StateDict,
                perturbers::BSR.StateDict)
    open("$directory/target", "w") do file
        ntarg = sum([length(targets[t]) for t in keys(targets)])
        partial_waves = sort(collect(keys(perturbers)))
        # nlsp = length(partial_waves)
        kpert = sum([length(perturbers[p]) for p in keys(perturbers)])
        write(file, heading, "\n")
        separator(file)
        property(file, "coupling", "LS", "coupling scheme")
        property(file, "nz", Z, "nuclear charge")
        property(file, "nelc", nelc, "number of bound electrons")
        separator(file)
        property(file, "ntarg", ntarg, "number of targets")
        separator(file)
        for k in keys(targets)
            for t in targets[k]
                write(file, t, ".c\n")
            end
        end
        separator(file)
        property(file, "nlsp", nlsp, "number of partial waves")
        separator(file)
        for pw in partial_waves
            write(file, "$(string(pw)) $(pw.L) $(multiplicity(pw)) $(pw.parity)\n")
        end
        separator(file)
        property(file, "kpert", kpert, "number of perturbers")
        separator(file)
        for ti = eachindex(partial_waves)
            for p in perturbers[partial_waves[ti]]
                write(file, "$ti $p\n")
            end
        end
    end
    readall("$directory/target")
end

function load_target(directory::AbstractString)
    open(joinpath(directory, "target")) do file
        target = Dict()
        read_parm = T -> begin
            line = split(readline(file), "=")
            val = strip(split(line[2], "!")[1])
            target[symbol(strip(line[1]))] = T<:AbstractString ? val : parse(T, val)
        end
        expect_line = () -> begin
            line = readline(file)
            line[1]==line[2]=='-' || error("Expected divider")
        end
        readline(file)
        expect_line()
        coupling = read_parm(ASCIIString)
        read_parm(Int)
        read_parm(Int)
        expect_line()
        ntarg = read_parm(Int)
        expect_line()
        target[:targets] = []
        for i = 1:ntarg
            line = split(readline(file))
            t = Dict(:name => line[1],
                     :parity => parse(Int, line[5]),
                     :E => parse(Float64, line[6]),
                     :ncfg => parse(Int, line[7]))
            if coupling == "LS"
                t[:L] = parse(Int, line[3])
                t[:S] = parse(Int, line[4])
            else
                t[:J] = parse(Int, line[3])//2
            end
            push!(target[:targets], t)
        end
        expect_line()
        read_parm(Int)
        read_parm(Int)
        read_parm(Int)
        expect_line()
        nlsp = read_parm(Int)
        expect_line()
        target[:partial_waves] = map(1:nlsp) do i
            line = split(readline(file))
            a,b = parse(Int, line[2]),parse(Int, line[3])
            pw = Dict(:name => line[1],
                      :parity => parse(Int, line[4]))
            if coupling == "LS"
                pw[:L] = a
                pw[:S] = b
            else
                pw[:J] = a//2
            end
            pw
        end
        expect_line()
        kpert = read_parm(Int)
        expect_line()
        target[:perturbers] = Dict()
        for i = 1:kpert
            line = split(readline(file))
            pk = parse(Int, line[1])
            target[:perturbers][pk] = vcat(get(target[:perturbers], pk, [])..., line[2])
        end
        expect_line()
        readline(file)
        expect_line()
        for i = 1:nlsp
            line = split(readline(file))
            nch = parse(Int, line[5])
            target[:partial_waves][i][:channels] = map(1:nch) do j
                line = split(readline(file))
                Dict(:name => line[1],
                     :l => parse(Int, line[2]),
                     :target => parse(Int, line[3]),
                     :last_cfg => parse(Int, line[5]),
                     :K => (parse(Int, line[6])-1)//2)
            end
            expect_line()
        end
        target
    end
end

@bsr_exe bsr_prep
@bsr_exe bsr_conf

export knot_set, load_knot_set, target, load_target, bsr_prep, bsr_conf


text:
#=
All units are cgs.
=#

module StarStruc

using Plots
using LaTeXStrings
using Logging

@Logging.configure(level=INFO)

# provides c0, G, m_H, k_B, sigma_sb, a_rad
include("constants.jl")

type Star
    M::AbstractFloat
    m::AbstractFloat
    mf::AbstractFloat
    X::AbstractFloat
    Y::AbstractFloat
    grid::Array{AbstractFloat}
    r::Array{AbstractFloat}
    l::Array{AbstractFloat}
    P::Array{AbstractFloat}
    T::Array{AbstractFloat}
end

Star(M) = Star(M, 0.01 * Msun, 0.8 * M, 0.7, 0.28,
               zeros(1), zeros(1), zeros(1), zeros(1), zeros(1))
Star(M, X, Y) = Star(M, 0.01 * Msun, 0.8 * M, X, Y,
                     zeros(1), zeros(1), zeros(1), zeros(1), zeros(1))

include("opacity.jl")
include("shootf.jl")
include("nr.jl")
include("tests.jl")

function starplot(star::Star)
    p1 = plot(star.grid / Msun, star.r / Rsun, label="r",
              xlabel=L"Mass [$M_\odot$]", ylabel=L"Radius [$R_\odot$]")
    p2 = plot(star.grid / Msun, star.l / Lsun, label="l",
              xlabel=L"Mass [$M_\odot$]", ylabel=L"Luminosity [$L_\odot$]")
    p3 = plot(star.grid / Msun, star.P, label="P", yscale=:log10,
              xlabel=L"Mass [$M_\odot$]", ylabel=L"Pressure [dyne cm$^{-2}$]")
    p4 = plot(star.grid / Msun, star.T, label="T", yscale=:log10,
              xlabel=L"Mass [$M_\odot$]", ylabel="Temperature [K]")
    plot(p1, p2, p3, p4, layout=(2, 2))
end

export shootf
export shootf!
export profiles
export Msun
export Rsun
export Star
export starplot

end # module


text:
# Autogenerated wrapper script for lttng_tools_jll for x86_64-linux-gnu
export liblttng_ctl, lttng, lttng_consumerd, lttng_crash, lttng_relayd, lttng_sessiond

using URCU_jll
using Popt_jll
using XML2_jll
using lttng_ust_jll
JLLWrappers.@generate_wrapper_header("lttng_tools")
JLLWrappers.@declare_library_product(liblttng_ctl, "liblttng-ctl.so.0")
JLLWrappers.@declare_executable_product(lttng)
JLLWrappers.@declare_executable_product(lttng_consumerd)
JLLWrappers.@declare_executable_product(lttng_crash)
JLLWrappers.@declare_executable_product(lttng_relayd)
JLLWrappers.@declare_executable_product(lttng_sessiond)
function __init__()
    JLLWrappers.@generate_init_header(URCU_jll, Popt_jll, XML2_jll, lttng_ust_jll)
    JLLWrappers.@init_library_product(
        liblttng_ctl,
        "lib/liblttng-ctl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

    JLLWrappers.@init_executable_product(
        lttng_consumerd,
        "lib/lttng/libexec/lttng-consumerd",
    )

    JLLWrappers.@init_executable_product(
        lttng_crash,
        "bin/lttng-crash",
    )

    JLLWrappers.@init_executable_product(
        lttng_relayd,
        "bin/lttng-relayd",
    )

    JLLWrappers.@init_executable_product(
        lttng_sessiond,
        "bin/lttng-sessiond",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    links_from_webserver_index(url)

Extracts all the content links from a webservers directory index page.
These follow a pretty standard form.
This one is tested so far on Apache/2.2.15
"""
function get_dataurls_from_webserver_index(datapage_url)
    datapage = getpage(datapage_url)

    data_links = eachmatch(sel"td a", datapage.root)
    @assert( text_only(data_links[1]) == "Parent Directory", repr(data_links[1]))

    data_hrefs = getindex.(attrs.(data_links[2:end]), "href")
    data_urls = joinpath.(datapage_url, data_hrefs)
end

"""
    citation_text(doi)

Uses the DOI formatted citation service to generate citation text for a given DOI.
This works for DOI's issued by: CrossRef, DataCite, and mEDRA.

See https://citation.crosscite.org/docs.html#sec-4-1
"""
function citation_text(doi)
    # GOLDPLATE: this could support so much more for different styles, but we don't need it
    url = startswith(doi, "http") ? doi : joinpath("https://doi.org/", doi)
    resp = HTTP.get(url, ["Accept"=>"text/x-bibliography"]; forwardheaders=true)
    resp.body |> String |> strip
end


text:
using Test
using Plots
using LinearAlgebra

using GCMMesh
using GCMMesh.Mesh

using GCMMesh.Graphs

FT = Float64

x1min, x1max = FT(0), FT(1)
x2min, x2max = FT(0), FT(1)

limits = (x1min, x2min, x1max, x2max)

function is_in_fiedler_vector_space(
    eigvec,
    graph::Graph{I,FT},
    evec_tol = FT(1.0E-6),
) where {I,FT}
    nv = graph.nverts
    vl, vu = FT(1), FT(1) # not used
    il, iu = 1, nv
    st, en = 2, 2
    abstol = FT(1.0E-14)    # abs tolerance for syevr! LAPACK solver
    eig_reltol = FT(1.0E-4) # relative tolerance for finding fiedler eigenvalues with algebraic
    # multiplicity > 1
    lmat = Graphs.laplacian_matrix_full(graph)
    eigvals, eigvecs =
        LinearAlgebra.LAPACK.syevr!('V', 'I', 'L', lmat, vl, vu, il, iu, abstol)

    fdlr_eig = eigvals[2]
    for i = 3:nv
        if abs(fdlr_eig - eigvals[i]) < eig_reltol
            en = i
        else
            break
        end
    end
    diff = zeros(FT, nv)
    for i = st:en
        diff .+= (eigvec' * eigvecs[:, i]) .* eigvecs[:, i]
    end
    return norm(diff .- eigvec) < evec_tol
end

@testset "graph for a simple 4 x 4 rectangular grid" begin
    @testset "face connectivity graph for a 4×4 element quad mesh" begin
        mesh = equispaced_rectangular_mesh(limits..., 4, 4, (false, false))
        fgraph = build_face_graph(mesh)
        off = fgraph.edge_offset

        @test fgraph.nverts == 4 * 4
        @test fgraph.edge_data[off[1]:off[2]-1] == [2, 5]
        @test fgraph.edge_data[off[2]:off[3]-1] == [1, 3, 6]
        @test fgraph.edge_data[off[3]:off[4]-1] == [2, 4, 7]
        @test fgraph.edge_data[off[4]:off[5]-1] == [3, 8]
        @test fgraph.edge_data[off[5]:off[6]-1] == [1, 6, 9]
        @test fgraph.edge_data[off[6]:off[7]-1] == [2, 5, 7, 10]
        @test fgraph.edge_data[off[7]:off[8]-1] == [3, 6, 8, 11]
        @test fgraph.edge_data[off[8]:off[9]-1] == [4, 7, 12]
        @test fgraph.edge_data[off[9]:off[10]-1] == [5, 10, 13]
        @test fgraph.edge_data[off[10]:off[11]-1] == [6, 9, 11, 14]
        @test fgraph.edge_data[off[11]:off[12]-1] == [7, 10, 12, 15]
        @test fgraph.edge_data[off[12]:off[13]-1] == [8, 11, 16]
        @test fgraph.edge_data[off[13]:off[14]-1] == [9, 14]
        @test fgraph.edge_data[off[14]:off[15]-1] == [10, 13, 15]
        @test fgraph.edge_data[off[15]:off[16]-1] == [11, 14, 16]
        @test fgraph.edge_data[off[16]:off[17]-1] == [12, 15]
    end

    @testset "face connectivity graph for a 4×4 element quad mesh, periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 4, 4, (true, true))
        fgraph = build_face_graph(mesh)
        off = fgraph.edge_offset

        @test fgraph.nverts == 4 * 4
        @test fgraph.edge_data[off[1]:off[2]-1] == [2, 4, 5, 13]
        @test fgraph.edge_data[off[2]:off[3]-1] == [1, 3, 6, 14]
        @test fgraph.edge_data[off[3]:off[4]-1] == [2, 4, 7, 15]
        @test fgraph.edge_data[off[4]:off[5]-1] == [1, 3, 8, 16]
        @test fgraph.edge_data[off[5]:off[6]-1] == [1, 6, 8, 9]
        @test fgraph.edge_data[off[6]:off[7]-1] == [2, 5, 7, 10]
        @test fgraph.edge_data[off[7]:off[8]-1] == [3, 6, 8, 11]
        @test fgraph.edge_data[off[8]:off[9]-1] == [4, 5, 7, 12]
        @test fgraph.edge_data[off[9]:off[10]-1] == [5, 10, 12, 13]
        @test fgraph.edge_data[off[10]:off[11]-1] == [6, 9, 11, 14]
        @test fgraph.edge_data[off[11]:off[12]-1] == [7, 10, 12, 15]
        @test fgraph.edge_data[off[12]:off[13]-1] == [8, 9, 11, 16]
        @test fgraph.edge_data[off[13]:off[14]-1] == [1, 9, 14, 16]
        @test fgraph.edge_data[off[14]:off[15]-1] == [2, 10, 13, 15]
        @test fgraph.edge_data[off[15]:off[16]-1] == [3, 11, 14, 16]
        @test fgraph.edge_data[off[16]:off[17]-1] == [4, 12, 13, 15]
    end

    @testset "vertex connectivity graph for a 4×4 element quad mesh" begin
        mesh = equispaced_rectangular_mesh(limits..., 4, 4, (false, false))
        vgraph = build_vertex_graph(mesh)
        off = vgraph.edge_offset

        @test vgraph.nverts == 4 * 4
        @test vgraph.edge_data[off[1]:off[2]-1] == [2, 5, 6]
        @test vgraph.edge_data[off[2]:off[3]-1] == [1, 3, 5, 6, 7]
        @test vgraph.edge_data[off[3]:off[4]-1] == [2, 4, 6, 7, 8]
        @test vgraph.edge_data[off[4]:off[5]-1] == [3, 7, 8]
        @test vgraph.edge_data[off[5]:off[6]-1] == [1, 2, 6, 9, 10]
        @test vgraph.edge_data[off[6]:off[7]-1] == [1, 2, 3, 5, 7, 9, 10, 11]
        @test vgraph.edge_data[off[7]:off[8]-1] == [2, 3, 4, 6, 8, 10, 11, 12]
        @test vgraph.edge_data[off[8]:off[9]-1] == [3, 4, 7, 11, 12]
        @test vgraph.edge_data[off[9]:off[10]-1] == [5, 6, 10, 13, 14]
        @test vgraph.edge_data[off[10]:off[11]-1] == [5, 6, 7, 9, 11, 13, 14, 15]
        @test vgraph.edge_data[off[11]:off[12]-1] == [6, 7, 8, 10, 12, 14, 15, 16]
        @test vgraph.edge_data[off[12]:off[13]-1] == [7, 8, 11, 15, 16]
        @test vgraph.edge_data[off[13]:off[14]-1] == [9, 10, 14]
        @test vgraph.edge_data[off[14]:off[15]-1] == [9, 10, 11, 13, 15]
        @test vgraph.edge_data[off[15]:off[16]-1] == [10, 11, 12, 14, 16]
        @test vgraph.edge_data[off[16]:off[17]-1] == [11, 12, 15]
    end

    @testset "vertex connectivity graph for a 4×4 element quad mesh" begin
        mesh = equispaced_rectangular_mesh(limits..., 4, 4, (true, true))
        vgraph = build_vertex_graph(mesh)
        off = vgraph.edge_offset

        @test vgraph.nverts == 4 * 4
        @test vgraph.edge_data[off[1]:off[2]-1] == [2, 4, 5, 6, 8, 13, 14, 16]
        @test vgraph.edge_data[off[2]:off[3]-1] == [1, 3, 5, 6, 7, 13, 14, 15]
        @test vgraph.edge_data[off[3]:off[4]-1] == [2, 4, 6, 7, 8, 14, 15, 16]
        @test vgraph.edge_data[off[4]:off[5]-1] == [1, 3, 5, 7, 8, 13, 15, 16]
        @test vgraph.edge_data[off[5]:off[6]-1] == [1, 2, 4, 6, 8, 9, 10, 12]
        @test vgraph.edge_data[off[6]:off[7]-1] == [1, 2, 3, 5, 7, 9, 10, 11]
        @test vgraph.edge_data[off[7]:off[8]-1] == [2, 3, 4, 6, 8, 10, 11, 12]
        @test vgraph.edge_data[off[8]:off[9]-1] == [1, 3, 4, 5, 7, 9, 11, 12]
        @test vgraph.edge_data[off[9]:off[10]-1] == [5, 6, 8, 10, 12, 13, 14, 16]
        @test vgraph.edge_data[off[10]:off[11]-1] == [5, 6, 7, 9, 11, 13, 14, 15]
        @test vgraph.edge_data[off[11]:off[12]-1] == [6, 7, 8, 10, 12, 14, 15, 16]
        @test vgraph.edge_data[off[12]:off[13]-1] == [5, 7, 8, 9, 11, 13, 15, 16]
        @test vgraph.edge_data[off[13]:off[14]-1] == [1, 2, 4, 9, 10, 12, 14, 16]
        @test vgraph.edge_data[off[14]:off[15]-1] == [1, 2, 3, 9, 10, 11, 13, 15]
        @test vgraph.edge_data[off[15]:off[16]-1] == [2, 3, 4, 10, 11, 12, 14, 16]
        @test vgraph.edge_data[off[16]:off[17]-1] == [1, 3, 4, 9, 11, 12, 13, 15]
    end
end

@testset "graph for a cube panel mesh" begin
    @testset "face connectivity graph for a 6-element cube panel mesh" begin
        mesh = cube_panel_mesh(1, FT)
        fgraph = build_face_graph(mesh)
        off = fgraph.edge_offset

        @test fgraph.nverts == 6
        @test fgraph.edge_data[off[1]:off[2]-1] == [2, 3, 4, 5]
        @test fgraph.edge_data[off[2]:off[3]-1] == [1, 3, 5, 6]
        @test fgraph.edge_data[off[3]:off[4]-1] == [1, 2, 4, 6]
        @test fgraph.edge_data[off[4]:off[5]-1] == [1, 3, 5, 6]
        @test fgraph.edge_data[off[5]:off[6]-1] == [1, 2, 4, 6]
        @test fgraph.edge_data[off[6]:off[7]-1] == [2, 3, 4, 5]
    end

    @testset "vertex connectivity graph for a 6-element cube panel mesh" begin
        mesh = cube_panel_mesh(1, FT)
        vgraph = build_vertex_graph(mesh)
        off = vgraph.edge_offset

        @test vgraph.nverts == 6
        @test vgraph.edge_data[off[1]:off[2]-1] == [2, 3, 4, 5]
        @test vgraph.edge_data[off[2]:off[3]-1] == [1, 3, 5, 6]
        @test vgraph.edge_data[off[3]:off[4]-1] == [1, 2, 4, 6]
        @test vgraph.edge_data[off[4]:off[5]-1] == [1, 3, 5, 6]
        @test vgraph.edge_data[off[5]:off[6]-1] == [1, 2, 4, 6]
        @test vgraph.edge_data[off[6]:off[7]-1] == [2, 3, 4, 5]
    end

    @testset "face connectivity graph for 24 element cube panel mesh with 4 elements per panel" begin
        mesh = cube_panel_mesh(2, FT)
        fgraph = build_face_graph(mesh)
        off = fgraph.edge_offset

        @test fgraph.nverts == 24
        @test fgraph.edge_data[off[1]:off[2]-1] == [2, 3, 13, 17]
    end

    @testset "face connectivity graph for 24 element cube panel mesh with 4 elements per panel" begin
        mesh = cube_panel_mesh(2, FT)
        vgraph = build_vertex_graph(mesh)
        off = vgraph.edge_offset

        @test vgraph.nverts == 24
        @test vgraph.edge_data[off[1]:off[2]-1] == [2, 3, 4, 13, 15, 17, 18]
    end
end

@testset "graph laplacian matrix x vector" begin
    @testset "for a 4×4 element quad mesh" begin
        mesh = equispaced_rectangular_mesh(limits..., 4, 4, (false, false))
        fgraph = build_face_graph(mesh)
        nverts = fgraph.nverts
        vin = ones(FT, nverts)
        vout = Vector{FT}(undef, nverts)
        laplacian_x_v!(vout, vin, fgraph)
        @test vout == zeros(FT, nverts)
    end

    @testset "for a 24 element cube panel mesh" begin
        mesh = cube_panel_mesh(2, FT)
        fgraph = build_face_graph(mesh)
        nverts = fgraph.nverts
        vin = ones(FT, nverts)
        vout = Vector{FT}(undef, nverts)
        laplacian_x_v!(vout, vin, fgraph)
        @test vout == zeros(FT, nverts)
    end
end

@testset "test split graph" begin
    @testset "for a 4×4 element quad mesh" begin
        mesh = equispaced_rectangular_mesh(limits..., 4, 4, (false, false))
        graph = build_face_graph(mesh)
        nv = graph.nverts
        nv1 = Int(cld(nv, 2))
        g1, g2 = split_graph(graph, nv1)
        @test g1.nverts == nv1         # check if the graph size is right
        @test g2.nverts == nv - nv1
        @test isempty(intersect(g1.vertices, g2.vertices)) # no common vertices
        @test sort(graph.vertices) == sort(union(g1.vertices, g2.vertices))
    end

    @testset "for a 24 element cube panel mesh" begin
        mesh = cube_panel_mesh(2, FT)
        graph = build_vertex_graph(mesh)
        nv = graph.nverts
        nv1 = Int(cld(nv, 2))
        g1, g2 = split_graph(graph, nv1)
        @test g1.nverts == nv1         # check if the graph size is right
        @test g2.nverts == nv - nv1
        @test isempty(intersect(g1.vertices, g2.vertices)) # no common vertices
        @test sort(graph.vertices) == sort(union(g1.vertices, g2.vertices))
    end

end

@testset "test Fiedler vector calculation using iterative methods" begin
    @testset "for a 4×4 element quad mesh" begin
        mesh = equispaced_rectangular_mesh(limits..., 4, 4, (false, false))
        graph = build_face_graph(mesh)
        fdlr_direct = fiedler_vector(graph, Graphs.Direct())
        fdlr_power = fiedler_vector(graph, Graphs.PowerIt())
        fdlr_lanczos = fiedler_vector(graph, Graphs.Lanczos())
        @test is_in_fiedler_vector_space(fdlr_direct, graph)
        @test is_in_fiedler_vector_space(fdlr_power, graph)
        @test is_in_fiedler_vector_space(fdlr_lanczos, graph)
    end
end

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
module GradientTests

using DiffTests, ForwardDiff, ReverseDiff, Test

include(joinpath(dirname(@__FILE__), "../utils.jl"))

function test_unary_gradient(f, x)
    test = ForwardDiff.gradient!(DiffResults.GradientResult(x), f, x)

    # without GradientConfig

    test_approx(ReverseDiff.gradient(f, x), DiffResults.gradient(test))

    out = similar(x)
    ReverseDiff.gradient!(out, f, x)
    test_approx(out, DiffResults.gradient(test))

    result = DiffResults.GradientResult(x)
    ReverseDiff.gradient!(result, f, x)
    test_approx(DiffResults.value(result), DiffResults.value(test))
    test_approx(DiffResults.gradient(result), DiffResults.gradient(test))

    # with GradientConfig

    cfg = ReverseDiff.GradientConfig(x)

    test_approx(ReverseDiff.gradient(f, x, cfg), DiffResults.gradient(test))

    out = similar(x)
    ReverseDiff.gradient!(out, f, x, cfg)
    test_approx(out, DiffResults.gradient(test))

    result = DiffResults.GradientResult(x)
    ReverseDiff.gradient!(result, f, x, cfg)
    test_approx(DiffResults.value(result), DiffResults.value(test))
    test_approx(DiffResults.gradient(result), DiffResults.gradient(test))

    # with GradientTape

    seedx = rand(eltype(x), size(x))
    tp = ReverseDiff.GradientTape(f, seedx)

    test_approx(ReverseDiff.gradient!(tp, x), DiffResults.gradient(test))

    out = similar(x)
    ReverseDiff.gradient!(out, tp, x)
    test_approx(out, DiffResults.gradient(test))

    result = DiffResults.GradientResult(x)
    ReverseDiff.gradient!(result, tp, x)
    test_approx(DiffResults.value(result), DiffResults.value(test))
    test_approx(DiffResults.gradient(result), DiffResults.gradient(test))

    # with compiled GradientTape

    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
        ctp = ReverseDiff.compile(tp)

        test_approx(ReverseDiff.gradient!(ctp, x), DiffResults.gradient(test))

        out = similar(x)
        ReverseDiff.gradient!(out, ctp, x)
        test_approx(out, DiffResults.gradient(test))

        result = DiffResults.GradientResult(x)
        ReverseDiff.gradient!(result, ctp, x)
        test_approx(DiffResults.value(result), DiffResults.value(test))
        test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
    end
end

function test_ternary_gradient(f, a, b, c)
    test_val = f(a, b, c)
    test_a = ForwardDiff.gradient(x -> f(x, b, c), a)
    test_b = ForwardDiff.gradient(x -> f(a, x, c), b)
    test_c = ForwardDiff.gradient(x -> f(a, b, x), c)

    # without GradientConfig

    ∇a, ∇b, ∇c = ReverseDiff.gradient(f, (a, b, c))
    test_approx(∇a, test_a)
    test_approx(∇b, test_b)
    test_approx(∇c, test_c)

    ∇a, ∇b, ∇c = map(similar, (a, b, c))
    ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c))
    test_approx(∇a, test_a)
    test_approx(∇b, test_b)
    test_approx(∇c, test_c)

    ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
    ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c))
    test_approx(DiffResults.value(∇a), test_val)
    test_approx(DiffResults.value(∇b), test_val)
    test_approx(DiffResults.value(∇c), test_val)
    test_approx(DiffResults.gradient(∇a), test_a)
    test_approx(DiffResults.gradient(∇b), test_b)
    test_approx(DiffResults.gradient(∇c), test_c)

    # with GradientConfig

    cfg = ReverseDiff.GradientConfig((a, b, c))

    ∇a, ∇b, ∇c = ReverseDiff.gradient(f, (a, b, c), cfg)
    test_approx(∇a, test_a)
    test_approx(∇b, test_b)
    test_approx(∇c, test_c)

    ∇a, ∇b, ∇c = map(similar, (a, b, c))
    ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c), cfg)
    test_approx(∇a, test_a)
    test_approx(∇b, test_b)
    test_approx(∇c, test_c)

    ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
    ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c), cfg)
    test_approx(DiffResults.value(∇a), test_val)
    test_approx(DiffResults.value(∇b), test_val)
    test_approx(DiffResults.value(∇c), test_val)
    test_approx(DiffResults.gradient(∇a), test_a)
    test_approx(DiffResults.gradient(∇b), test_b)
    test_approx(DiffResults.gradient(∇c), test_c)

    # with GradientTape

    tp = ReverseDiff.GradientTape(f, (rand(eltype(a), size(a)), rand(eltype(b), size(b)), rand(eltype(c), size(c))))

    ∇a, ∇b, ∇c = ReverseDiff.gradient!(tp, (a, b, c))
    test_approx(∇a, test_a)
    test_approx(∇b, test_b)
    test_approx(∇c, test_c)

    ∇a, ∇b, ∇c = map(similar, (a, b, c))
    ReverseDiff.gradient!((∇a, ∇b, ∇c), tp, (a, b, c))
    test_approx(∇a, test_a)
    test_approx(∇b, test_b)
    test_approx(∇c, test_c)

    ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
    ReverseDiff.gradient!((∇a, ∇b, ∇c), tp, (a, b, c))
    test_approx(DiffResults.value(∇a), test_val)
    test_approx(DiffResults.value(∇b), test_val)
    test_approx(DiffResults.value(∇c), test_val)
    test_approx(DiffResults.gradient(∇a), test_a)
    test_approx(DiffResults.gradient(∇b), test_b)
    test_approx(DiffResults.gradient(∇c), test_c)

    # with compiled GradientTape

    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
        ctp = ReverseDiff.compile(tp)

        ∇a, ∇b, ∇c = ReverseDiff.gradient!(ctp, (a, b, c))
        test_approx(∇a, test_a)
        test_approx(∇b, test_b)
        test_approx(∇c, test_c)

        ∇a, ∇b, ∇c = map(similar, (a, b, c))
        ReverseDiff.gradient!((∇a, ∇b, ∇c), ctp, (a, b, c))
        test_approx(∇a, test_a)
        test_approx(∇b, test_b)
        test_approx(∇c, test_c)

        ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
        ReverseDiff.gradient!((∇a, ∇b, ∇c), ctp, (a, b, c))
        test_approx(DiffResults.value(∇a), test_val)
        test_approx(DiffResults.value(∇b), test_val)
        test_approx(DiffResults.value(∇c), test_val)
        test_approx(DiffResults.gradient(∇a), test_a)
        test_approx(DiffResults.gradient(∇b), test_b)
        test_approx(DiffResults.gradient(∇c), test_c)
    end
end

# issue https://github.com/JuliaDiff/ReverseDiff.jl/issues/140
nested_array_mul_140(x) = sum(sum(x[1] * [[x[2], x[3]]]))
test_println("Issue #140", nested_array_mul_140)
test_unary_gradient(nested_array_mul_140, [1.0, 2.0, 1.0, -2.4, 4.0])

for f in DiffTests.MATRIX_TO_NUMBER_FUNCS
    test_println("MATRIX_TO_NUMBER_FUNCS", f)
    test_unary_gradient(f, rand(5, 5))
end

for f in DiffTests.VECTOR_TO_NUMBER_FUNCS
    test_println("VECTOR_TO_NUMBER_FUNCS", f)
    test_unary_gradient(f, rand(5))
end

for f in DiffTests.TERNARY_MATRIX_TO_NUMBER_FUNCS
    test_println("TERNARY_MATRIX_TO_NUMBER_FUNCS", f)
    test_ternary_gradient(f, rand(5, 5), rand(5, 5), rand(5, 5))
end

end # module