text:
#
# `depth_first_constructive` is a very simple recursive implementation
# but it is more difficult to cap the costs so scales very badly
#

function depth_first_constructive(T::Vector{<: ITensor}) where {LabelSetT}
  indsT = [inds(Tₙ) for Tₙ in T]
  return depth_first_constructive(DimT, indsT)
end

function depth_first_constructive(::Type{DimT}, T::Vector{IndexSetT}) where {IndexSetT <: IndexSet, DimT}
  labels, dims = label_dims(DimT, T)
  nlabels = length(dims)
  if nlabels ≤ 16
    return depth_first_constructive(UInt16, labels, dims)
  elseif nlabels ≤ 32
    return depth_first_constructive(UInt32, labels, dims)
  elseif nlabels ≤ 64
    return depth_first_constructive(UInt64, labels, dims)
  elseif nlabels ≤ 128
    return depth_first_constructive(UInt128, labels, dims)
  else
    return depth_first_constructive(BitSet, labels, dims)
  end
end

function depth_first_constructive(::Type{LabelSetT}, labels::Vector, dims::Vector) where {LabelSetT}
  return depth_first_constructive(map(label -> bitset(LabelSetT, label), labels), dims)
end

function depth_first_constructive(::Type{LabelSetT}, ::Type{DimT}, T::Vector{<: ITensor}) where {LabelSetT, DimT}
  indsT = [inds(Tₙ) for Tₙ in T]
  return depth_first_constructive(LabelSetT, DimT, indsT)
end

function depth_first_constructive(::Type{LabelSetT}, ::Type{DimT}, T::Vector{IndexSetT}) where {IndexSetT <: IndexSet, LabelSetT, DimT}
  labels, dims = label_dims(DimT, T)
  return depth_first_constructive(map(label -> bitset(LabelSetT, label), labels), dims)
end

function depth_first_constructive(T::Vector, ind_dims::Vector)
  optimal_cost = Ref(typemax(eltype(ind_dims)))
  optimal_sequence = Vector{Pair{Int, Int}}(undef, length(T)-1)
  _depth_first_constructive!(optimal_sequence, optimal_cost, Pair{Int, Int}[], T, ind_dims, collect(1:length(T)), 0)
  return pair_sequence_to_tree(optimal_sequence, length(T))
end

function _depth_first_constructive!(optimal_sequence, optimal_cost, sequence, T, ind_dims, remaining, cost)
  if length(remaining) == 1
    # Only should get here if the contraction was the best
    # Otherwise it would have hit the `continue` below
    @assert cost ≤ optimal_cost[]
    optimal_cost[] = cost
    optimal_sequence .= sequence
  end
  for aᵢ in 1:length(remaining)-1, bᵢ in aᵢ+1:length(remaining)
    a = remaining[aᵢ]
    b = remaining[bᵢ]
    current_cost, Tᵈ = contraction_cost(T[a], T[b], ind_dims)
    new_cost = cost + current_cost
    if new_cost ≥ optimal_cost[]
      continue
    end
    new_sequence = push!(copy(sequence), a => b)
    new_T = push!(copy(T), Tᵈ)
    new_remaining = deleteat!(copy(remaining), (aᵢ, bᵢ))
    push!(new_remaining, length(new_T))
    _depth_first_constructive!(optimal_sequence, optimal_cost, new_sequence, new_T, ind_dims, new_remaining, new_cost)
  end
end


text:
#@testset "Fluxi.Basics_01" begin

    using MetabolicFluxAnalysis
    lca = MetabolicFluxAnalysis.parseLC(" 3.2 A + 1   B + C + 2   D  ")

    s1 = MetabolicFluxAnalysis.Species("s1",4)
    s2 = MetabolicFluxAnalysis.Species("s2",0)
    s3 = MetabolicFluxAnalysis.Species("a2",3)
    ss = [s1;s2;s3]
    #print( sort(ss) )

    sas1 = MetabolicFluxAnalysis.SAS("sa",1)
    sas2 = MetabolicFluxAnalysis.SAS("sa",2)
    sas3 = MetabolicFluxAnalysis.SAS("sb",1)
    sas4 = MetabolicFluxAnalysis.SAS("sa",2)

    isequal(sas2,sas4)

    hash(sas2)==hash(sas4)
    ss1 = Set( [sas1;sas2;sas3;sas4] )
    @test length(ss1)==3

    sas1a = MetabolicFluxAnalysis.SAS("sa",1)
    sas2a = MetabolicFluxAnalysis.SAS("sa",2)
    sas1b = MetabolicFluxAnalysis.SAS("sa",1)
    sas2b = MetabolicFluxAnalysis.SAS("sa",2)

    ss1 = Set( [sas1a;sas2a] )
    ss2 = Set( [sas1b;sas2b] )
    @test ss1==ss2

    ras1 = MetabolicFluxAnalysis.RAS("f1",MetabolicFluxAnalysis.SAS("sa",1),1)
    ras2 = MetabolicFluxAnalysis.RAS("f1",MetabolicFluxAnalysis.SAS("sb",1),1)
    ras3 = MetabolicFluxAnalysis.RAS("f1",MetabolicFluxAnalysis.SAS("sa",1),1)
    ps1 = Set([ras1;ras2;ras3])
    @test length(ps1)==2

    pas1 = MetabolicFluxAnalysis.PAS("f1",MetabolicFluxAnalysis.SAS("sa",1),1)
    pas2 = MetabolicFluxAnalysis.PAS("f1",MetabolicFluxAnalysis.SAS("sb",1),1)
    pas3 = MetabolicFluxAnalysis.PAS("f1",MetabolicFluxAnalysis.SAS("sa",1),1)
    ps1 = Set([pas1;pas2;pas3])
    @test length(ps1)==2

    m1 = MetabolicFluxAnalysis.ASTEMUMassRatio( MetabolicFluxAnalysis.SpeciesEMU( [MetabolicFluxAnalysis.SAS("A",1);MetabolicFluxAnalysis.SAS("A",3)] ) , 1.0 )
    m2 = MetabolicFluxAnalysis.ASTEMUMassRatio( MetabolicFluxAnalysis.SpeciesEMU( [MetabolicFluxAnalysis.SAS("A",1);MetabolicFluxAnalysis.SAS("A",3)] ) , 1.0 )
    @test hash(m1) == hash(m2)
    @test isequal( m1.semu , m2.semu )
    @test isequal(m1,m2)

#end

ext: jl


text:
function give_all_states(N::Int, L::Int, Lp::Int)

    # states = Array{String}( undef, (Lp + 1)^N )
    states = String[]

    # 以 Lp + 1 进制表示，生成最大值
    max_s = string(Lp, base = Lp + 1)
    for i in 2:N
        max_s = max_s * string(Lp, base = Lp + 1)
    end

    # 从0到max_s，找出所有符合粒子总数条件的状态
    for i in 0:parse(Int, max_s, base = Lp + 1)
        # 进制转换为 Lp + 1
        the_state = string(i, base = Lp + 1, pad = N)
        # 粒子数符合条件则记录
        if total_num_particle(the_state, N, Lp) == L
            push!(states, the_state)
        end
    end
    # parse(Int, "444", base =5)
    # string(4444, base = 5)
    states
    # max_s
end

@memoize function total_num_particle(the_state::String, N::Int, Lp::Int)
    num = 0
    for i in 1:N
        # num = num + parse(Int, the_state[i], base = Lp+1)
        num = num + num_of_position(the_state[i], Lp)
    end
    num
end

@memoize function num_of_position(the_state_of_i::Char, Lp::Int)
    parse(Int, the_state_of_i, base = Lp + 1)
end

function fill_H_μ_U!(H, the_model::BH_Model)

    for j in 1:length(the_model.states)
        the_sum = 0.0
        for i in 1:the_model.N
            the_num = num_of_position(the_model.states[j][i], the_model.Lp)
            the_sum = the_sum - (0.5*the_model.U + the_model.μ)*the_num + 0.5*the_model.U*(the_num^2)
        end
        H[j, j] = H[j, j] + the_sum
    end

    true
end

function fill_H_J!(H, the_model::BH_Model)
    for j in 1:length(the_model.states)
        for i in 1:the_model.N

            index, quan = give_index_cr_an(the_model.states, j, i, i + 1, the_model.N, the_model.Lp)
            if index != 0
                H[index, j] = H[index, j] - the_model.J*quan
                # println("$i-cr, $(i + 1)-an $(the_model.states[j]) -> $(the_model.states[index]) | ($j, $index), $quan")  # debug
            end
            index, quan = give_index_cr_an(the_model.states, j, i + 1, i, the_model.N, the_model.Lp)
            if index != 0
                H[index, j] = H[index, j] - the_model.J*quan
                # println("$(i + 1)-cr, $i-an $(the_model.states[j]) -> $(the_model.states[index]) |($j, $index), $quan")  # debug
            end
        end

    end
    true

end

# i产生, j湮灭
@memoize function give_index_cr_an(states::Vector{String}, label::Int, i::Int, j::Int, N::Int, Lp::Int)
    the_state_s = states[label]
    i, j = adjust_index_for_pbc(i, j, N)
    i_num = num_of_position(the_state_s[end - i + 1], Lp)
    j_num = num_of_position(the_state_s[end - j + 1], Lp)
    if j_num == 0 || i_num == Lp
        return 0, 0.0
    end
    the_state_n = parse(Int, the_state_s, base = Lp + 1)
    the_state_n = the_state_n + (Lp + 1)^(i - 1) - (Lp + 1)^(j - 1)
    the_state_s = string(the_state_n, base = Lp + 1, pad = N)
    # findall(x -> x == the_state_s, states), i_num, j_num, states[label], the_state_s # debug

    findall(x -> x == the_state_s, states)[1], sqrt(i_num + 1)*sqrt(j_num) # release

end

function fill_H_V!(H, the_model::BH_Model)

    for j in 1:length(the_model.states)
        the_sum = 0.0
        for i in 1:the_model.N
            x, y = adjust_index_for_pbc(i, i + 1, the_model.N)
            the_sum = the_sum + the_model.V*num_of_position(the_model.states[j][x], the_model.Lp)*num_of_position(the_model.states[j][y], the_model.Lp)
        end
        H[j, j] = H[j, j] + the_sum
    end

    # for i in 1:length(the_model.states)
    #     i, j = adjust_index_for_pbc(i, i + 1, the_model.N)
    #     H[i, i] = H[i, i] + the_model.V*total_num_particle(the_model.states[i], the_model.N, the_model.Lp)*total_num_particle(the_model.states[j], the_model.N, the_model.Lp)
    # end
    true

end

function fill_H_ε!(H, the_model::BH_Model)
    for j in 1:length(the_model.states)

        for i in 1:the_model.N
            index_1 , v1, index_2, v2 = give_index_cr_and_an_2times(the_model.states, j, i, the_model.Lp, the_model.N)
            if index_1 != 0
                H[index_1, j] = H[index_1, j] - 0.5*the_model.ε*v1
            end
            if index_2 != 0
                H[index_2, j] = H[index_2, j] - 0.5*the_model.ε*v2
            end
        end

    end
    true

end

@memoize function give_index_cr_and_an_2times(states::Vector{String}, label::Int, i::Int, Lp::Int, N::Int)
    the_state_s = states[label]
    i_num = num_of_position(the_state_s[end - i + 1], Lp)
    index_1 , v1, index_2, v2 = 0, 0.0, 0, 0.0
    # 产生两次
    if i_num <= Lp - 2
        the_state_1_n = parse(Int, the_state_s, base = Lp + 1)
        the_state_1_n = the_state_1_n + 2*(Lp + 1)^(i - 1)
        the_state_1_s = string(the_state_1_n, base = Lp + 1, pad = N)
        index_1, v1 = findall(x -> x == the_state_1_s, states)[1], sqrt(i_num + 1)*sqrt(i_num + 2)
        # println("cr 2 times, $the_state_s -> $the_state_1_s, ($index_1, $label) = $v1") # debug
    end
    # 湮灭两次
    if i_num >= 2
        the_state_2_n = parse(Int, the_state_s, base = Lp + 1)
        the_state_2_n = the_state_2_n - 2*(Lp + 1)^(i - 1)
        the_state_2_s = string(the_state_2_n, base = Lp + 1, pad = N)
        index_2, v2 = findall(x -> x == the_state_2_s, states)[1], sqrt(i_num)*sqrt(i_num - 1)
        # println("an 2 times, $the_state_s -> $the_state_2_s, ($index_2, $label) = $v2") # debug
    end
    index_1 , v1, index_2, v2
end

# 根据周期性边界条件调整临近作用坐标
@memoize function adjust_index_for_pbc(i::Int, j::Int, max::Int)
    if i > max
        # print("adjust_index_for_pbc, ($i, $j) -> ") #debug
        i = 1
    end
    if j > max
        # print("adjust_index_for_pbc, ($i, $j) -> ") #debug
        j = 1
    end
    # println("($i, $j)") #debug
    i, j
end

# 求特征向量的粒子数
function give_the_num_of_particles_of_eigvector(vec::Vector{Float64}, the_model::BH_Model)
    probability = vec.^2
    nhil = length(the_model.states)
    num_of_p = Array{Float64}(undef, nhil)
    for i in 1:nhil
        num_of_p[i] = total_num_particle(the_model.states[i], the_model.N, the_model.Lp)
    end
    sum(probability .* num_of_p)
end

function c_ff(the_model::BH_Model, vec::Vector{Float64})

    # 周期性边界条件
    pbc = Array{Int}(undef, 2*the_model.N)
    pbc[1:the_model.N] .= [i for i in 1:the_model.N]
    pbc[1 + the_model.N:end] .= [i for i in 1:the_model.N]

    # 希尔伯特空间维数
    nhil = length(the_model.states)

    # 大矩阵
    mat = Array{Float64}(undef, nhil, nhil)

    # 关联函数计算结果存放
    result = Array{Float64}(undef, the_model.N + 1)

    # 开始计算关联函数
    for r0 in 0:the_model.N-1
        # 置零
        mat .= 0.0
        for j in 1:nhil
            the_sum = 0.0
            for i in 1:the_model.N
                # 周期性边界条件
                x, y = i, pbc[i + r0]
                the_sum = the_sum + num_of_position(the_model.states[j][x], the_model.Lp)*num_of_position(the_model.states[j][y], the_model.Lp)
            end
            mat[j, j] = mat[j, j] + the_sum
        end
        result[r0 + 1] = vec' * mat * vec
    end

    result[end] = result[1]

    result ./ the_model.N
end


text:
# Practical Machine learning
# Reinforcement learning - Q learning example
# Chapter 12

using NNGraph, DeepQLearning
reload("DeepQLearning")

dqn = DeepQLearning.DQN(10,100,5)

s0 = randNNMat(10,1)
a = DeepQLearning.forward(dqn, s0)
DeepQLearning.act(dqn,s0)
DeepQLearning.learn(dqn, 0.)

s1 = randNNMat(10,1)
a = DeepQLearning.forward(dqn, s1)
DeepQLearning.act(dqn,s1)
DeepQLearning.learn(dqn, 0.)


text:
"""
Produce a temporary name that will not conflict with its Julian or Vulkan version.
"""
function tmp_argname(name, type)
    "_" * name
end

function wrap!(w_api, sdef::SDefinition)
    new_sdef = structure(sdef)
    has_bag(sdef.name) && setindex!(w_api.bags, create_bag(sdef), bagtype(sdef.name))
    wrap_structure!(w_api, new_sdef)
    wrap_constructor!(w_api, new_sdef, sdef)
end

function wrap_structure!(w_api, new_sdef)
    if new_sdef.name ∉ keys(w_api.structs)
        w_api.structs[new_sdef.name] = new_sdef # add to wrapped structs
    else
        @warn "Wrapped type $(new_sdef.name) but was already processed"
    end
end

optional_create_info_types = Dict(
    "VkImage" => "VkImageCreateInfo"
)

function structure(sdef)
    sname = sdef.name
    new_fields = OrderedDict()
    if is_handle(sname)
        new_fields["handle"] = sname
        abstract_type = "Handle"
    elseif dfmatch(vulkan_structs, :name, sname).returnedonly
        abstract_type = "ReturnedOnly"
        for (name, type) ∈ sdef.fields
            drop_field(name, type, sname) && continue
            new_name, new_type = field_transform(name, type, sname)
            new_fields[new_name] = new_type
        end
    else
        abstract_type = "VulkanStruct"
        new_fields["vks"] = sname
        has_bag(sname) && setindex!(new_fields, bagtype(sname), "bag")
    end
    new_name = name_transform(sdef)
    SDefinition(new_name, is_handle(sdef.name), new_fields, abstract_type, nothing)
end

# parameters to keep in a constructor rather than in a struct
# const parameters_dict = Dict(
#     "pNext" => KeywordArgument("next", "C_NULL"),
#     "flags" => KeywordArgument("flags", "0"),
#     "pAllocator" => KeywordArgument("pAllocator", "C_NULL"),
# )

drop_field(name, type, sname) = name == "sType" || is_count_variable(name, sname)
keeps_original_layout(sdef) = !has_bag(sdef.name) && "sType" ∉ keys(sdef.fields)

instance_name(sdef) = convert(SnakeCaseLower, sdef.name)

function create_bag(sdef)
    sname = sdef.name
    fields = OrderedDict()
    for (name, type) ∈ sdef.fields
        new_name, new_type = field_transform(name, type, sname)
        tmp_name = tmp_argname(name, type)
        if is_ptr(type)
            eltype = inner_type(type)
            has_bag(eltype) && setindex!(fields, "", "bag_" * new_name)
            name == "pNext" && setindex!(fields, "", "bag_next")
            eltype == "Cstring" && setindex!(fields, "", new_name * "_ptrarray")
        end
        if is_ptr(type) || type == "Cstring"
            field_name = (type ≠ "Cstring" && is_vulkan_struct(inner_type(type))) || name == "pNext" ? new_name * "_ref" : new_name
            fields[field_name] = ""
        end
    end
    SDefinition(bagtype(sname), false, fields, "Bag", nothing)
end

"Convert a Ptr{T} to a fetch_known_type(T)."
function convert_nested_type(type)
    eltype = inner_type(type)
    is_ptr(type) && is_vulkan_type(eltype) && return fieldtype_transform("", eltype)
    is_ptr(type) && return "Ptr{$(fieldtype_transform("", eltype))}"
    is_ntuple(type) && return replace(type, eltype => fieldtype_transform("", eltype))
    error("Nested type $type is neither a pointer nor an array (pointer to pointer) type")
end

function type_dependencies(sdef::SDefinition)
    type_deps = filter(x -> x ≠ sdef.name, map((x, y) -> isempty(x) ? [y] : x, type_dependencies.(values(sdef.fields)), values(sdef.fields)))
    isempty(type_deps) ? String[] : unique(vcat(type_deps...))
end

function field_transform(name, type, sname)
    new_type = fieldtype_transform(name, type, sname)
    fieldname_transform(name, type) => new_type isa Converted ? new_type.final_type : new_type
end

convention_exceptions = Dict(
    "textureCompressionASTC_LDR" => "texture_compression_astc_ldr",
    "textureCompressionASTC_HDR" => "texture_compression_astc_hdr",
    "formatA4R4G4B4" => "format_a4r4g4b4",
    "formatA4B4G4R4" => "format_a4b4g4r4",
)

function fieldname_transform(name, type)
    name ∈ keys(convention_exceptions) && return convention_exceptions[name]
    startswith(name, "pfn") && return convert(SnakeCaseLower, CamelCaseUpper(name[4:end])).value
    if !isnothing(match(r"^\s*p+[A-Z]", name))
        cc = CamelCaseUpper(lstrip(name, 'p'))
    else
        cc = detect_convention(name, instance=true)
    end
    convert(VulkanGen.julia_convention[:variable], cc).value
end

function fieldtype_transform(name, type, sname = nothing)
    (isnothing(type) || type == "Any") && return "Any"
    startswith(type, "PFN_") && return "Function"
    occursin("version", lowercase(name)) && type == "UInt32" && return Converted(type, "VersionNumber")
    is_ntuple(type) && inner_type(type) == "UInt8" && return Converted(type, "String")
    @return_value_if_key_exists type_conversions type
    !isnothing(sname) && !isempty(name) && is_array_variable(name, sname) && return "Array{$(fieldtype_transform("", inner_type(type), sname))}"
    is_vulkan_type(type) && return remove_vk_prefix(type)
    is_ptr(type) && return convert_nested_type(type)
    type
end

name_transform(str, ::Type{SDefinition}) = str[3:end]


text:
#################### Posterior Statistics ####################

function autocor(c::AbstractChains; lags::Vector=[1, 5, 10, 50],
                 relative::Bool=true)
  if relative
    lags *= step(c)
  elseif any(lags .% step(c) .!= 0)
    throw(ArgumentError("lags do not correspond to thinning interval"))
  end
  labels = map(x -> "Lag " * string(x), lags)
  vals = mapslices(x -> autocor(x, lags)', c.value, dims=[1, 2])
  ChainSummary(vals, c.names, labels, header(c))
end

function cor(c::AbstractChains)
  ChainSummary(cor(combine(c)), c.names, c.names, header(c))
end

function changerate(c::AbstractChains)
  n, p, m = size(c.value)
  r = zeros(Float64, p, 1, 1)
  r_mv = 0.0
  delta = Array{Bool}(undef, p)
  for k in 1:m
    prev = c.value[1, :, k]
    for i in 2:n
      for j in 1:p
        x = c.value[i, j, k]
        dx = x != prev[j]
        r[j] += dx
        delta[j] = dx
        prev[j] = x
      end
      r_mv += any(delta)
    end
  end
  vals = round.([r; r_mv] / (m * (n - 1)), digits=3)
  ChainSummary(vals, [c.names; "Multivariate"], ["Change Rate"], header(c))
end

describe(c::AbstractChains; args...) = describe(stdout, c; args...)

function describe(io::IO, c::AbstractChains;
                  q::Vector=[0.025, 0.25, 0.5, 0.75, 0.975], etype=:bm, args...)
  ps_stats = summarystats(c; etype=etype, args...)
  ps_quantiles = quantile(c, q=q)
  println(io, ps_stats.header)
  print(io, "Empirical Posterior Estimates:\n")
  show(io, ps_stats)
  print(io, "Quantiles:\n")
  show(io, ps_quantiles)
end

function hpd(x::Vector{T}; alpha::Real=0.05) where {T<:Real}
  n = length(x)
  m = max(1, ceil(Int, alpha * n))

  y = sort(x)
  a = y[1:m]
  b = y[(n - m + 1):n]
  _, i = findmin(b - a)

  [a[i], b[i]]
end

function hpd(c::AbstractChains; alpha::Real=0.05)
  pct = first(showoff([100.0 * (1.0 - alpha)]))
  labels = ["$(pct)% Lower", "$(pct)% Upper"]
  vals = permutedims(
    mapslices(x -> hpd(vec(x), alpha=alpha), c.value, dims=[1, 3]),
    [2, 1, 3]
  )
  ChainSummary(vals, c.names, labels, header(c))
end

function quantile(c::AbstractChains; q::Vector=[0.025, 0.25, 0.5, 0.75, 0.975])
  labels = map(x -> string(100 * x) * "%", q)
  vals = permutedims(
    mapslices(x -> quantile(vec(x), q), c.value, dims=[1, 3]),
    [2, 1, 3]
  )
  ChainSummary(vals, c.names, labels, header(c))
end

function summarystats(c::AbstractChains; etype=:bm, args...)
  f = x -> [mean(x), std(x), sem(x), mcse(vec(x), etype; args...)]
  labels = ["Mean", "SD", "Naive SE", "MCSE", "ESS"]
  vals = permutedims(
    mapslices(x -> f(x), c.value, dims=[1, 3]),
    [2, 1, 3]
  )
  stats = [vals  min.((vals[:, 2] ./ vals[:, 4]).^2, size(c.value, 1))]
  ChainSummary(stats, c.names, labels, header(c))
end


text:
function header(f::IOStream)
    sdf_magic = simple_str(read(f, 4))
    endianness = read(f, Int32)
    file_version = read(f, Int32)
    file_revision = read(f, Int32)
    code_name = simple_str(read(f, ID_LENGTH))
    first_block_location = read(f, Int64)
    summary_location = read(f, Int64)
    summary_size = read(f, Int32)
    nblocks = read(f, Int32)
    block_header_length = read(f, Int32)
    step = read(f, Int32)
    time = read(f, Float64)
    jobid1 = read(f, Int32)
    jobid2 = read(f, Int32)
    string_length = read(f, Int32)
    code_io_version = read(f, Int32)

    Header(
        sdf_magic,
        endianness,
        file_version,
        file_revision,
        code_name,
        first_block_location,
        summary_location,
        summary_size,
        nblocks,
        block_header_length,
        step,
        time,
        jobid1,
        jobid2,
        string_length,
        code_io_version,
    )
end

header(filename::AbstractString) = open(header, filename)

function Base.read(f, block::BlockHeader{T,D,BLOCKTYPE_CONSTANT}) where {T,D}
    val = read(f, T)
    ConstantBlockHeader(block, val)
end

function Base.read(f, block::BlockHeader{T,N,BLOCKTYPE_PLAIN_MESH}) where {T,N}
    dims = Array{Int32, 1}(undef, N)

    mults, labels, units, geometry, minval, maxval = read_mesh_common!(f,N)
    read!(f, dims)

    PlainMeshBlockHeader(
        block,
        mults,
        labels,
        units,
        geometry,
        minval,
        maxval,
        dims,
    )
end

function Base.read(f, block::BlockHeader{T,N,BLOCKTYPE_POINT_MESH}) where {T,N}
    mults, labels, units, geometry, minval, maxval = read_mesh_common!(f,N)
    np = read(f, Int64)

    PointMeshBlockHeader(
        block,
        mults,
        labels,
        units,
        geometry,
        minval,
        maxval,
        np,
    )
end

function Base.read(f, block::BlockHeader{T,D,BLOCKTYPE_PLAIN_VARIABLE}) where {T,D}
    dims = Array{Int32, 1}(undef, D)

    mult, units, mesh_id = read_variable_common!(f)
    read!(f, dims)
    stagger = Stagger(read(f, Int32))

    PlainVariableBlockHeader(block, mult, units, mesh_id, Tuple(dims), stagger)
end

function Base.read(f, block::BlockHeader{T,D,BLOCKTYPE_POINT_VARIABLE}) where {T,D}
    mult, units, mesh_id = read_variable_common!(f)
    np = read(f, Int64)

    PointVariableBlockHeader(block, mult, units, mesh_id, np)
end

function Base.read(f, block::BlockHeader{T,D,BLOCKTYPE_RUN_INFO}) where {T,D}
    code_version = read(f, Int32)
    code_revision = read(f, Int32)
    commit_id = simple_str(read(f, ID_LENGTH))
    sha1sum = simple_str(read(f, ID_LENGTH))
    compile_machine = simple_str(read(f, ID_LENGTH))
    compile_flags = simple_str(read(f, ID_LENGTH))
    defines = read(f, Int64)
    compile_date = read(f, Int32)
    run_date = read(f, Int32)
    io_date = read(f, Int32)
    # TODO save these fields

    RunInfoBlockHeader(block)
end

function Base.read(f, block::BlockHeader{T,D,BLOCKTYPE_CPU_SPLIT}) where {T,D}
    CPUSplitBlockHeader(block)
end

function read_variable_common!(f)
    mult = read(f, Float64)
    units = simple_str(read(f, ID_LENGTH))
    mesh_id = simple_str(read(f, ID_LENGTH))

    mult, units, mesh_id
end

function read_mesh_common!(f, n)
    minval = Array{Float64, 1}(undef, n)
    maxval = Array{Float64, 1}(undef, n)

    mults = ntuple(n) do i
        read(f, Float64)
    end

    labels = ntuple(n) do i
        simple_str(read(f, ID_LENGTH))
    end
    units = ntuple(n) do i
        simple_str(read(f, ID_LENGTH))
    end
    geometry = Geometry(read(f, Int32))
    read!(f, minval)
    read!(f, maxval)

    mults, labels, units, geometry, minval, maxval
end

simple_str(s) = replace(String(rstrip(String(s))), "\0" => "")


text:
using Clang.Generators
using BlueZ_jll

include_dir = joinpath(BlueZ_jll.artifact_dir ,"include")

options = load_options(joinpath(@__DIR__, "wrap.toml"))

args = get_default_args()
push!(args, "-I$include_dir")

bluetooth = joinpath(include_dir, "bluetooth")

headers = [joinpath(bluetooth, header) for header in ["bluetooth.h"]]

ctx = create_context(headers, args, options)
build!(ctx)


text:
# This file was generated, do not modify it. # hide
y, X = unpack(iris, ==(:Species), !=(:PetalLength))
first(X, 1) |> pretty

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
### Types
abstract type ApproximationScheme end

struct RosenfeldFMT{T} <: ApproximationScheme
    A::T
    B::T
    G::T
end

struct PercusYevick{T} <: ApproximationScheme
    α::T
    β::T
    δ::T
end

struct VerletWeis{T} <: ApproximationScheme
    coreliquid::HardSpheres{T}
    subscheme::PercusYevick{T}
    α::T
end

struct SharmaSharma{S, T} <: ApproximationScheme
    coreliquid::HardSpheres{T}
    subscheme::S
end

struct NonLinearSharmaSharma{S, T} <: ApproximationScheme
    coreliquid::HardSpheres{T}
    subscheme::S
end

struct MSA{S, T} <: ApproximationScheme
    coreliquid::HardSpheres{T}
    subscheme::S
    α₀::T
    α₁::T
    α₂::T
    α₃::T
    β₀::T
    β₁::T
    β₂::T
    β₃::T
    a₁::T
    a₃::T
    b₁::T
    b₃::T
end

### Constructors
function RosenfeldFMT(liquid::HardDisks)
    η  = volume_fraction(liquid)
    η² = η * η
    η³ = η² * η
    η⁴ = η³ * η
    η₋⁻¹ = 1 / (1 - η)
    η₋⁻² = η₋⁻¹ * η₋⁻¹

    Z  = (1 + 0.128η² + 0.027η³ + 0.06η⁴) * η₋⁻²
    Z′ = ((0.256η + 0.081η² + 0.24η³) * η₋⁻²) + (2Z * η₋⁻¹)
    χ  = Z + η * Z′

    G = sqrt(Z′ / 2)
    A = (1 + (2η - 1) * χ + 2η * G) / η
    B = (-1 + (1 - η) * χ - 3η * G) / η

    return RosenfeldFMT{typeof(A)}(A, B, G)
end

function PercusYevick(liquid::HardSpheres)
    η = volume_fraction(liquid)
    η₀⁴ = @fastmath (1 - η)^4
    η₁² = (1 + 2η)^2
    η₂² = (2 +  η)^2

    α = -η₁² / η₀⁴
    β = 3η * η₂² / 2η₀⁴
    δ = -η * η₁² / 2η₀⁴

    return PercusYevick{typeof(α)}(α, β, δ)
end

function VerletWeis(liquid::HardSpheres)
    η  = volume_fraction(liquid)
    κ  = 1 - η / 16
    η′ = η * κ
    α  = ∛(κ)

    coreliquid = HardSpheres(η′)
    subscheme = PercusYevick(coreliquid)

    return VerletWeis{typeof(α)}(coreliquid, subscheme, α)
end

function SharmaSharma{S}(liquid::AttractiveHardSpheres{TT}) where
    {S<:ApproximationScheme, TT}

    η = volume_fraction(liquid)
    coreliquid = HardSpheres(η)
    subscheme = S(coreliquid)
    SS = typeof(subscheme)

    return SharmaSharma{SS, TT}(coreliquid, subscheme)
end

function NonLinearSharmaSharma{S}(liquid::AttractiveHardSpheres{TT}) where
    {S <: ApproximationScheme, TT}

    η = volume_fraction(liquid)
    coreliquid = HardSpheres(η)
    subscheme = S(coreliquid)
    SS = typeof(subscheme)

    return NonLinearSharmaSharma{SS, TT}(coreliquid, subscheme)
end

function MSA{S}(liquid::DipolarHardSpheres{TT}; tol = √(eps(TT))) where
    {S<:ApproximationScheme, TT}

    η = volume_fraction(liquid)
    T = temperature(liquid)
    coreliquid = HardSpheres(η)
    subscheme = S(coreliquid)
    SS = typeof(subscheme)
    κ = dhs_msa_parameter(liquid, tol)

    T⁻¹ = 1 / T
    ξ   = κ * η

    ξ₁² = (1 - 2ξ)^2
    ξ₁⁴ = ξ₁² * ξ₁²
    ξ₂⁴ = (1 +  ξ)^4
    ξ₃² = (1 + 4ξ)^2

    r₁ = ξ₃² / ξ₁⁴
    r₂ = ξ₁² / ξ₂⁴
    r₃ = 8 * (1 + ξ)^2 / ξ₁⁴
    r₄ = (2 - ξ)^2 / ξ₂⁴

    α₀′ = 2κ * (-r₁ + r₂)
    α₁′ = 3ξ * κ * ( r₃ + r₄)
    α₃′ = -ξ * κ * (2r₁ + r₂)
    β₁′ = 3ξ * κ / 8 * (2r₃ - r₄)
    β₃′ = -ξ * κ / 4 * (4r₁ - r₂)

    α₀ = α₀′ / 9  + α₁′ / 12  + α₃′ / 18
    α₁ = α₀′ / 90 + α₁′ / 108 + α₃′ / 144
    α₂ = α₀′ +  α₁′ +  α₃′
    α₃ = α₀′ + 2α₁′ + 4α₃′

    β₀ = -T⁻¹ / 9
    β₁ =  β₁′ / 270 + β₃′ / 360 + β₀ / 10
    β₂ =  β₁′ +  β₃′ - T⁻¹
    β₃ = 4β₁′ + 6β₃′ + β₂

    a₁ = 2 * (α₁′ + 8β₁′) / 3
    a₃ = 8 * (α₃′ + 4β₃′)
    b₁ = 2 * (α₁′ - 4β₁′) / 3
    b₃ = 8 * (α₃′ - 2β₃′)

    return MSA{SS, TT}(coreliquid, subscheme,
                       α₀, α₁, α₂, α₃, β₀, β₁, β₂, β₃, a₁, a₃, b₁, b₃)
end


text:
include("setup.jl")

tests = ["radon", "fbp", "geometry"]

if !isempty(ARGS)
    tests = ARGS
end

@testset "Tests for the MartaCT package" begin
    for t ∈ tests
        include("test_$t.jl")
    end
end


text:
export ingredients
export WidthOverDocs

"""
   ingredients(path) allows Pluto notebooks to include a module from another file.
It's like include, but for inside a pluto notebook.
Source: https://github.com/fonsp/Pluto.jl/issues/115
"""
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

""" Provides checkbox to toggle full width versus narrow with column for LiveDocs """
function WidthOverDocs(enabled::Bool=false)  # From PlutoThemes.jl
	checked = wide ? "checked" : ""
	init = wide ? "toggle_width(document.getElementById('width-over-livedocs'))" : ""
	return HTML("""
<!-- https://github.com/fonsp/Pluto.jl/issues/400#issuecomment-695040745 -->
<input
	type="checkbox"
	id="width-over-livedocs"
	name="width-over-livedocs"
    onclick="window.plutoOptIns.toggle_width(this)"
	$(checked)>
<label for="width-over-livedocs">
	Full Width Mode
</label>
<style>
	body.width-over-docs #helpbox-wrapper {
    	display: none !important;
	}
	body.width-over-docs main {
		max-width: calc(100% - 4rem);
		margin-right: 2rem;
	}
</style>
<script>
	const toggle_width = function(t) {
		t.checked
		? document.body.classList.add("width-over-docs")
		: document.body.classList.remove("width-over-docs") }
	window.plutoOptIns = window.plutoOptIns || {}
	window.plutoOptIns.toggle_width = toggle_width
	$(init)
</script>
""")
end


text:
"""
    c = @sc_ircode f(args...)

Julia SSA IR explore.

```julia
c                  # view IR in the REPL
display(c)         # (ditto)

c.cfg              # control-flow graph (CFG) visualizer
c.cfg_only         # CFG without IR in node label
c.dom              # dominator tree visualizer
c.dom_only         # dominator tree without IR in node label
display(c.cfg)     # display CFG

c.llvm             # create LLVM IR explore
c.native           # create native code explore
c.att              # (ditto)
c.intel            # create native code explore in intel syntax
edit(c.native)
abspath(c.native)
```

... and so on; type `c.` + TAB to see the full list.

Since visualizers such as `c.cfg` and `c.cfg_only` work via the standard
`show` mechanism, they interoperable well with other packages like FileIO.jl
and DisplayAs.jl.

```julia
using FileIO
save("PATH/TO/IMAGE.png", c.cfg_only)
save("PATH/TO/IMAGE.svg", c.cfg_only)
save("PATH/TO/IMAGE.pdf", c.cfg_only)

using DisplayAs
c.cfg_only |> DisplayAs.SVG
```
"""
:(@sc_ircode)

struct IRCodeView <: AbstractCode
    ir::Core.Compiler.IRCode
    f::Any
    atype::Any
    rtype::Union{Type,Nothing}
    args::Any
    kwargs::Any
end

macro sc_ircode(args...)
    gen_call_with_extracted_types_and_kwargs(__module__, sc_ircode, args)
end

function sc_ircode(f, argument_type; kwargs...)
    @nospecialize
    # TODO: handle multiple returns?
    (ir, rtype), = code_ircode(f, argument_type; kwargs...)
    args = (f, argument_type)
    return IRCodeView(ir, f, argument_type, rtype, args, kwargs)
end

function sc_ircode(mi::Core.Compiler.MethodInstance; kwargs...)
    @nospecialize
    mth = mi.def
    if mth isa Method
        ftype = Base.tuple_type_head(mth.sig)
        if Base.issingletontype(ftype)
            f = ftype.instance
        else
            f = ftype  # ?
        end
        atype = Base.tuple_type_tail(mth.sig)
    else
        f = "f?"
        atype = "Tuple{?}"
    end

    args = (mi,)
    ir, rtype = code_ircode(args...; kwargs...)
    return IRCodeView(ir, f, atype, rtype, args, kwargs)
end

function sc_ircode(ir::Core.Compiler.IRCode; kwargs...)
    f = "f?"
    atype = "Tuple{?}"
    rtype = nothing
    args = (ir,)
    return IRCodeView(ir, f, atype, rtype, args, kwargs)
end

function sc_ircode(ci::Core.Compiler.CodeInfo; kwargs...)
    ir = CompilerUtils.ircode_from_codeinfo(ci)
    f = ci.parent.def.name
    atype = Base.tuple_type_tail(ci.parent.def.sig)
    rtype = ci.rettype
    args = (ci,)
    return IRCodeView(ir, f, atype, rtype, args, kwargs)
end

function Base.summary(io::IO, llvm::IRCodeView)
    @unpack f, atype = Fields(llvm)
    print(io, "IRCodeView of ", f, " with ", atype)
    return
end

function show_ir_config(ir; debuginfo = :default)
    if debuginfo === :source
        return Base.IRShow.default_config(ir; verbose_linetable = true)
    else
        return Base.IRShow.default_config(ir; verbose_linetable = false)
    end
end

function _show_ir(io::IO, ir::Core.Compiler.IRCode; options...)
    if isempty(options)
        show(io, MIME"text/plain"(), ir)
        return
    end
    Base.IRShow.show_ir(io, ir, show_ir_config(ir; options...))
end

function Base.show(io::IO, ::MIME"text/plain", ircv::IRCodeView; options...)
    @unpack ir, rtype = Fields(ircv)
    summary(io, ircv)
    println(io)
    _show_ir(io, ir; options...)
    print(io, "⇒ ", rtype)
    println(io)
    return
end

Base.propertynames(::IRCodeView) = (
    :ir,
    :rtype,
    # explores:
    :llvm,
    :native,
    :intel,
    :att,
    # visualizers:
    :cfg,
    :cfg_only,
    :dom,
    :dom_only,
)

function Base.getproperty(ircv::IRCodeView, name::Symbol)
    @unpack args = Fields(ircv)
    if name === :llvm
        return sc_llvm(args...)
    elseif name === :native || name === :att
        return sc_native(args...)
    elseif name === :intel
        return sc_intel(args...)
    elseif name === :cfg
        return IRCodeCFGDot(ircv, true)
    elseif name === :cfg_only
        return IRCodeCFGDot(ircv, false)
    elseif name === :dom
        return IRCodeDomTree(ircv, true)
    elseif name === :dom_only
        return IRCodeDomTree(ircv, false)
    end
    return getfield(ircv, name)
end

abstract type AbstractLazyDot <: AbstractCode end

function dot_to_iobuffer(dot)
    io = IOBuffer()
    print_dot(io, dot)
    seekstart(io)
    return io
end

function run_dot(output::IO, input::IO, options)
    cmd = getcmd(:dot)
    cmd = `$cmd -Gfontname=monospace -Nfontname=monospace -Efontname=monospace $options`
    @debug "Run: $cmd"
    run(pipeline(cmd, stdout = output, stderr = stderr, stdin = input))
    return
end

# https://www.iana.org/assignments/media-types/text/vnd.graphviz
Base.show(io::IO, ::MIME"text/vnd.graphviz", dot::AbstractLazyDot) = print_dot(io, dot)

Base.show(io::IO, ::MIME"image/png", dot::AbstractLazyDot) =
    run_dot(io, dot_to_iobuffer(dot), `-Tpng`)
Base.show(io::IO, ::MIME"image/svg+xml", dot::AbstractLazyDot) =
    run_dot(io, dot_to_iobuffer(dot), `-Tsvg`)
Base.show(io::IO, ::MIME"application/pdf", dot::AbstractLazyDot) =
    run_dot(io, dot_to_iobuffer(dot), `-Tpdf`)

struct IRCodeCFGDot <: AbstractLazyDot
    ircv::IRCodeView
    include_code::Bool
end

function escape_dot_label(io::IO, str)
    for c in str
        if c in "\\{}<>|\"\n"
            # https://graphviz.org/doc/info/attrs.html#k:escString
            print(io, '\\', c)
        else
            print(io, c)
        end
    end
end

function Base.summary(io::IO, dot::IRCodeCFGDot)
    @unpack ircv = Fields(dot)
    @unpack f, atype = Fields(ircv)
    print(io, "CFG of $f on $atype")
end

function find_syncregions(ir, bb)
    ids = Int[]
    for i in bb.stmts
        inst = ir.stmts[i][:inst]
        if inst isa Expr && inst.head === :syncregion
            push!(ids, i)
        end
    end
    return ids
end

print_stmt(io, x) = print(io, x)

function print_stmt(io, ex::Expr)
    if Meta.isexpr(ex, :enter, 1)
        print(io, "enter #", ex.args[1])
    elseif ex.head in (:leave, :pop_exception)
        print(io, ex.head)
        for a in ex.args
            print(io, ' ')
            print(io, a)
        end
    else
        print(io, ex)
    end
end

function print_stmt(io::IO, x::Core.GotoNode)
    @unpack label = x
    print(io, "goto #", label)
end

function print_stmt(io::IO, x::Core.GotoIfNot)
    @unpack cond, dest = x
    print(io, "goto #", dest, " if not ", cond)
end

function print_bb_stmts_for_dot_label(io, ir, bb)
    for s in bb.stmts
        ln = sprint() do io
            print(io, "%", s, " = ")
            print_stmt(io, ir.stmts.inst[s])
        end
        escape_dot_label(io, ln)
        print(io, "\\l")
    end
end

print_dot(dot) = print_dot(stdout, dot)
function print_dot(io::IO, dot::IRCodeCFGDot)
    @unpack ircv, include_code = Fields(dot)
    @unpack ir = Fields(ircv)

    function bblabel(i)
        inst = ir.stmts.inst[ir.cfg.blocks[i].stmts[end]]
        if inst isa Core.ReturnNode
            if isdefined(inst, :val)
                return "#$(i)⏎"
            else
                return "#$(i)⚠"
            end
        end
        ids = find_syncregions(ir, ir.cfg.blocks[i])
        if !isempty(ids)
            return "#$i SR(" * join(("%$i" for i in ids), ", ") * ")"
        end
        return string("#", i)
    end

    graphname = summary(dot)
    print(io, "digraph \"")
    escape_dot_label(io, graphname)
    println(io, "\" {")
    indented(args...) = print(io, "    ", args...)
    indented("label=\"")
    escape_dot_label(io, graphname)
    println(io, "\";")
    for (i, bb) in enumerate(ir.cfg.blocks)
        indented(i, " [shape=record")

        # Print code
        if include_code
            print(io, ", label=\"{$(bblabel(i)):\\l")
        else
            print(io, ", label=\"{$(bblabel(i))}\", tooltip=\"")
        end
        print_bb_stmts_for_dot_label(io, ir, bb)
        if include_code
            print(io, "}\"")
        else
            print(io, '"')
        end
        println(io, "];")

        # Print edges
        term = ir.stmts[bb.stmts[end]][:inst]
        if term isa DetachNode && length(bb.succs) == 2 && term.label in bb.succs
            det, = (i for i in bb.succs if i != term.label)
            cont = term.label
            indented(i, " -> ", cont, " [label = \" C($(term.syncregion))\"];\n")
            indented(i, " -> ", det, " [label = \" D($(term.syncregion))\"];\n")
        elseif term isa ReattachNode && length(bb.succs) == 1 && bb.succs[1] == term.label
            indented(i, " -> ", term.label, " [label = \" R($(term.syncregion))\"];\n")
        elseif term isa SyncNode && length(bb.succs) == 1
            indented(i, " -> ", bb.succs[1], " [label = \" S($(term.syncregion))\"];\n")
        else
            if term isa Expr && term.head === :enter
                attr = "label = \" E\""
            elseif term isa Expr && term.head === :leave
                attr = "label = \" L\""
            else
                attr = ""
            end
            for s in bb.succs
                attr2 = i == s ? "dir = back " : ""
                indented(i, " -> ", s, " [", attr2, attr, "]", ";\n")
            end
        end
    end
    println(io, '}')
end

struct IRCodeDomTree <: AbstractLazyDot
    ircv::IRCodeView
    include_code::Bool
    domtree::Core.Compiler.DomTree
end

function IRCodeDomTree(ircv::IRCodeView, include_code::Bool)
    @unpack ir = Fields(ircv)
    domtree = Core.Compiler.construct_domtree(ir.cfg.blocks)
    return IRCodeDomTree(ircv, include_code, domtree)
end

function Base.summary(io::IO, d::IRCodeDomTree)
    @unpack f, atype = Fields(Fields(d).ircv)
    print(io, "Dominator tree for $f on $atype")
end

# https://github.com/JuliaDebug/Cthulhu.jl/issues/26
AbstractTrees.treekind(::IRCodeDomTree) = AbstractTrees.IndexedTree()
AbstractTrees.childindices(d::IRCodeDomTree, i::Int) = d[i].children
AbstractTrees.childindices(::IRCodeDomTree, ::IRCodeDomTree) = (1,)
AbstractTrees.parentlinks(::IRCodeDomTree) = AbstractTrees.StoredParents()
AbstractTrees.printnode(io::IO, i::Int, ::IRCodeDomTree) = print(io, i)
Base.getindex(d::IRCodeDomTree, i) = Fields(d).domtree.nodes[i]

function Base.show(io::IO, ::MIME"text/plain", d::IRCodeDomTree)
    summary(io, d)
    println(io)
    AbstractTrees.print_tree(io, 1; roottree = d)
end

function print_dot(io::IO, dot::IRCodeDomTree)
    @unpack ircv, domtree, include_code = Fields(dot)
    @unpack ir = Fields(ircv)

    graphname = summary(dot)
    print(io, "digraph \"")
    escape_dot_label(io, graphname)
    println(io, "\" {")
    indented(args...) = print(io, "    ", args...)
    indented("label=\"")
    escape_dot_label(io, graphname)
    println(io, "\";")

    @assert length(domtree.nodes) == length(ir.cfg.blocks)
    for (i, (node, bb)) in enumerate(zip(domtree.nodes, ir.cfg.blocks))
        indented(i, " [shape=record")

        # Print code
        if include_code
            print(io, ", label=\"{$i:\\l")
        else
            print(io, ", label=\"{$i}\", tooltip=\"")
        end
        print_bb_stmts_for_dot_label(io, ir, bb)
        if include_code
            print(io, "}\"")
        else
            print(io, '"')
        end
        println(io, "];")

        # Print edges
        for s in node.children
            indented(i, " -> ", s, ";\n")
        end
    end
    println(io, '}')
end


text:
#=
An optimization oracle is a linear optimization oracle for the set S, i.e., it solves
    z^\ast(c) :=    min c^Tw
                    s.t. w \in S
The solution of the above problem returned by the oracle is referred to as w^\ast(\cdot) in
the paper

A specific implementation of an optimization oracle should take a single input argument c,
and should return the pair (z^\ast(c), w^\ast(c))

This file contains the code for the assignment problem optimization oracle, which
relies on the Hungarian solver package.
=#

using Hungarian

"""
    assignment_oracle(c)

Optimization oracle for the assignment problem. Uses Hungarian algorithm. Assumes
that the input c is a vector of length d^2 and must be reshaped into a matrix
"""
function assignment_oracle(c::Vector{Float64})
    d_sq = length(c)
    d = isqrt(d_sq)
    C = reshape(c, (d,d))
    assignment, z = hungarian(C)

    W = zeros(d, d)
    for i = 1:d
        W[i, assignment[i]] = 1
    end
    w = vec(W)
    return (z, w)
end

"""
    assignment_hamming_reduce(c)

Utility function that converts c to a cost matrix with 0 on the optimal edges and 1 on all
other possible edges. Therefore, the SPO loss becomes equal to the hamming loss between the
predicted assignment and the optimal one. Returns the new cost vector c associated with
the Hamming loss.
"""
function assignment_hamming_reduce(c::Vector{Float64})
    d_sq = length(c)
    d = isqrt(d_sq)
    C = reshape(c, (d,d))
    assignment, z = hungarian(C)

    new_C = ones(d, d)
    for i = 1:d
        new_C[i, assignment[i]] = 0
    end
    new_c = vec(new_C)
    return new_c
end

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
abstract type AbstractDDSM <: Name end
struct DDSMMMAP <: AbstractDDSM end
struct DDSM <: AbstractDDSM end

const DDSM_PATH = Ref{String}(joinpath(DATADIR, "DDSM", "data.h5"))
const DDSM_FILE = Ref{HDF5.File}(HDF5.File(-1, "nnn"))

function set_DDSM_PATH!(path)
    DDSM_PATH[] = path
end

function open_ddsm()
    if !isfile(DDSM_PATH[])
        error("File $(DDSM_PATH[]) does not exists. Use DatasetProvider.set_DDSM_PATH! function to set the proper path.")
    end
    fid = DDSM_FILE[]
    if !isopen(fid)
        DDSM_FILE[] = fid = h5open(DDSM_PATH[], "r")
    end
    return fid
end

close_ddsm() = close(DDSM_FILE[])

# download options
task(::Type{<:AbstractDDSM}) = MultiClass
format(::Type{<:AbstractDDSM}) = GrayImages
source(::Type{<:AbstractDDSM}) = "https://www.kaggle.com/skooch/ddsm-mammography?select=test10_data"
make_datadep(::Type{<:AbstractDDSM}) = nothing

nclasses(::Type{<:AbstractDDSM}) = 4
classes(::Type{<:AbstractDDSM}) = collect(0:4)
nattributes(::Type{<:AbstractDDSM}) = (299, 299, 3)
ninstances(::Type{<:AbstractDDSM}) = (55890, 0, 15364)

# dataset description
name(::Type{<:AbstractDDSM}) = "DDSM Mammography"

function author(::Type{<:AbstractDDSM})
    return ["Eric A. Scuccimarra"]
end

# loading
function loadraw(N::Type{<:AbstractDDSM}, type)
    hassubset(N, type)
    split = type == :test ? "test" : "train"

    fid = open_ddsm()
    if N <: DDSMMMAP
        data = HDF5.readmmap(fid[split]["data"]) |> HDF5Array
    else
        @warn "due to the size of the DDSM dataset, it may take several minutes to load data into memory"
        data = HDF5.read(fid[split]["data"])
    end
    targets = fid[split]["labels"][:]
    close(fid)
    return data, targets
end


text:
using Test, LinearAlgebra
using DiffEqOperators
using ModelingToolkit

@testset "utility functions" begin
    @test DiffEqOperators.unit_indices(2) == (CartesianIndex(1,0), CartesianIndex(0,1))
    @test DiffEqOperators.add_dims(zeros(2,2), ndims(zeros(2,2)) + 2) == [6. 6.; 0. 0.; 0. 0.]
    @test DiffEqOperators.perpindex(collect(1:5), 3) == [1, 2, 4, 5]
    @test DiffEqOperators.perpsize(zeros(2,2,3,2), 3) == (2, 2, 2)
end

@testset "finite-difference weights from fornberg(1988) & fornberg(2020)" begin
    order = 2; z = 0.0; x = [-1, 0, 1.0];
    @test DiffEqOperators.calculate_weights(order, z, x) == [1,-2,1]  # central difference of second-derivative with unit-step

    order = 1; z = 0.0; x = [-1., 1.0];
    @test DiffEqOperators.calculate_weights(order, z, x) == [-0.5,0.5] # central difference of first-derivative with unit step

    order = 1; z = 0.0; x = [0, 1];
    @test DiffEqOperators.calculate_weights(order, z, x) == [-1, 1] # forward difference

    order = 1; z = 1.0; x = [0, 1];
    @test DiffEqOperators.calculate_weights(order, z, x) == [-1, 1] # backward difference

    # forward-diff of third derivative with order of accuracy == 3
    order = 3; z = 0.0; x = [0,1,2,3,4,5]
    @test DiffEqOperators.calculate_weights(order, z, x) == [-17/4,	71/4	,−59/2,	49/2,	−41/4,	7/4]

    order = 3; z = 0.0; x = collect(-3:3)
    d, e = DiffEqOperators.calculate_weights(order, z, x;dfdx = true)
    @test d ≈ [-167/18000, -963/2000, -171/16,0,171/16,963/2000,167/18000]
    @test e ≈ [-1/600,-27/200,-27/8,-49/3,-27/8,-27/200,-1/600]

    order = 3; z = 0.0; x = collect(-4:4)
    d, e = DiffEqOperators.calculate_weights(order, z, x;dfdx = true)
    @test d ≈ [-2493/5488000, -12944/385875, -87/125 ,-1392/125,0,1392/125,87/125,12944/385875,2493/5488000]
    @test e ≈ [-3/39200,-32/3675,-6/25,-96/25,-205/12, -96/25, -6/25,-32/3675,-3/39200]
end

@testset "count differentials 1D" begin
    @parameters t x
    @variables u(..)
    Dt = Differential(t)

    Dx = Differential(x)
    eq  = Dt(u(t,x)) ~ -Dx(u(t,x))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))

    Dxx = Differential(x)^2
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))

    Dxxxx = Differential(x)^4
    eq  = Dt(u(t,x)) ~ -Dxxxx(u(t,x))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 4
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
end

@testset "count differentials 2D" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dt = Differential(t)

    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
    @test first(DiffEqOperators.differential_order(eq.rhs, y.val)) == 2
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
    @test isempty(DiffEqOperators.differential_order(eq.lhs, y.val))
end

@testset "count with mixed terms" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dx = Differential(x)
    Dy = Differential(y)
    Dt = Differential(t)

    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + Dx(Dy(u(t,x,y)))
    @test DiffEqOperators.differential_order(eq.rhs, x.val) == Set([2, 1])
    @test DiffEqOperators.differential_order(eq.rhs, y.val) == Set([2, 1])
end

@testset "Kuramoto–Sivashinsky equation" begin
    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dx2 = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    α = 1
    β = 4
    γ = 1
    eq = Dt(u(x,t)) + u(x,t)*Dx(u(x,t)) + α*Dx2(u(x,t)) + β*Dx3(u(x,t)) + γ*Dx4(u(x,t)) ~ 0
    @test DiffEqOperators.differential_order(eq.lhs, x.val) == Set([4, 3, 2, 1])
end


text:
using GeometricIntegrators.Equations

ϑ₁(t, q) = q[4] + A₁(t,q)
ϑ₂(t, q) = q[5] + A₂(t,q)
ϑ₃(t, q) = q[6] + A₃(t,q)

function ϑ(t, q, p)
    p[1] = ϑ₁(t,q)
    p[2] = ϑ₂(t,q)
    p[3] = ϑ₃(t,q)
    p[4] = zero(eltype(q))
    p[5] = zero(eltype(q))
    p[6] = zero(eltype(q))
    nothing
end

# ϕ₀(x) = E₀*sin(2π*x)
# ϕ(t,q) = ϕ₀(q[3])

# E₁(t,q) = zero(eltype(q))
# E₂(t,q) = zero(eltype(q))
# E₃(t,q) = - 2π*E₀*cos(2π*q[3])

hamiltonian(t,q) = 0.5 * (q[4]^2 + q[5]^2 + q[6]^2) + ϕ(t, q)

angular_momentum(t,q) = q[1] * ϑ₂(t,q) - q[2] * ϑ₁(t,q)

function charged_particle_3d_pᵢ(qᵢ, tᵢ=0)
    pᵢ = zero(qᵢ)

    if ndims(qᵢ) == 1
        ϑ(tᵢ, qᵢ, pᵢ)
    else
        for i in 1:size(qᵢ,2)
            @views ϑ(tᵢ, qᵢ[:,i], pᵢ[:,i])
        end
    end
    pᵢ
end

charged_particle_3d_iode_ϑ(t, q, v, p) = ϑ(t, q, p)

function charged_particle_3d_iode_f(t, q, v, f)
    f[1] = dA₁dx₁(t,q) * v[1] + dA₂dx₁(t,q) * v[2] + dA₃dx₁(t,q) * v[3] + E₁(t,q)
    f[2] = dA₁dx₂(t,q) * v[1] + dA₂dx₂(t,q) * v[2] + dA₃dx₂(t,q) * v[3] + E₂(t,q)
    f[3] = dA₁dx₃(t,q) * v[1] + dA₂dx₃(t,q) * v[2] + dA₃dx₃(t,q) * v[3] + E₂(t,q)
    f[4] = v[1] - q[4]
    f[5] = v[2] - q[5]
    f[6] = v[3] - q[6]
    nothing
end

function charged_particle_3d_iode_g(t, q, λ, g)
    g[1] = dA₁dx₁(t,q) * λ[1] + dA₂dx₁(t,q) * λ[2] + dA₃dx₁(t,q) * λ[3]
    g[2] = dA₁dx₂(t,q) * λ[1] + dA₂dx₂(t,q) * λ[2] + dA₃dx₂(t,q) * λ[3]
    g[3] = dA₁dx₃(t,q) * λ[1] + dA₂dx₃(t,q) * λ[2] + dA₃dx₃(t,q) * λ[3]
    g[4] = λ[1]
    g[5] = λ[2]
    g[6] = λ[3]
    nothing
end

function charged_particle_3d_iode_v(t, q, v)
    v[1] = q[4]
    v[2] = q[5]
    v[3] = q[6]
    v[4] = E₁(t,q) + q[5] * B₃(t,q) - q[6] * B₂(t,q)
    v[5] = E₂(t,q) + q[6] * B₁(t,q) - q[4] * B₂(t,q)
    v[6] = E₃(t,q) + q[4] * B₂(t,q) - q[5] * B₁(t,q)
    nothing
end

function charged_particle_3d_iode(q₀=q₀)
    IODE(charged_particle_3d_iode_ϑ, charged_particle_3d_iode_f,
         charged_particle_3d_iode_g, q₀, charged_particle_3d_pᵢ(q₀);
         h=hamiltonian, v=charged_particle_3d_iode_v)
end


text:
module DensityRatioEstimation

using Requires, Parameters, Statistics, LinearAlgebra

include("moment_matching.jl")
export MMDAnalytical, MMDNumerical

export estimate_ratio

function __init__()
    # The path name `glue` comes from here: https://github.com/JuliaLang/Pkg.jl/issues/1285#issuecomment-525891481
    @require JuMP="4076af6c-e467-56ae-b986-b466b2749572" @require Ipopt="b6b21f68-93f8-5de0-b562-5493be1d77c9" include("glue/opt.jl")
end

end


text:
# get the first layer parameters for influence analysis
open("visualize.dot", "w") do io
      println(io, mx.to_graphviz(lstm))
  end
  run(pipeline(`dot -Tsvg visualize.dot`, stdout="visualize.svg"))

beta = copy(model.arg_params[:fullyconnected0_weight])
z = maximum(abs(beta),2);
cax3 = matshow(z', interpolation="nearest")
colorbar(cax3)
xlabels = [string(i) for i=1:40]
xticks(0:39, xlabels)


text:
using VDPTag
using Plots
using POMDPToolbox
using Reel
using ProgressMeter

frames = Frames(MIME("image/png"), fps=2)

mdp = VDPTagMDP()

hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(1))
# hist = simulate(hr, mdp, RandomPolicy(mdp, rng=MersenneTwister(2)))
hist = simulate(hr, mdp, ToNextML(mdp))

gr()
@showprogress "Creating gif..." for s in state_hist(hist)
    push!(frames, plot(mdp, s))
end

filename = string(tempname(), "_vdprun.gif")
write(filename, frames)
println(filename)
run(`setsid gifview $filename`)


text:
using Base.Test
using Clustering

D = [0 1 2 3
     1 0 1 2
     2 1 0 1
     3 2 1 0]

@assert size(D) == (4, 4)

a = [1, 1, 2, 2]
c = [2, 2]

@test_approx_eq silhouettes(a, c, D) [1.5/2.5, 0.5/1.5, 0.5/1.5, 1.5/2.5]

a = [1, 2, 1, 2]
c = [2, 2]

@test_approx_eq silhouettes(a, c, D) [0.0, -0.5, -0.5, 0.0]


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function genbo2(foldername, filename, u, xij, qi, qj, ti, tj, ai, aj, rho, mu, eta, kappa, natomtype)

if length(u)>0
    if (u[1]==0)
        m = 0;
        n = 0;
    else
        n,m = size(u);
        # if (n==0) || (n == 1)
        #     m = 0;
        # elseif (n == natomtype+1)
        #     if u[end,1] == 0
        #         m = 0;
        #     end
        # else
        #     error("Number of functions for the two-body bond order potential must be equal to 1 plus the number of atom types");
        # end
    end
else
    n = 0;
    m = 0;
end

gen = 1;
ifile = 0;
stropu = "";
strcpu = "";
strgpu = "";
potnum = 0;
fn = "PaircDensity";

for i = 1:m
    ui = u[:,i];
    if (ui[n] != 0)
        potnum = potnum + 1;
        stropu1, strcpu1, strgpu1 = genpair(foldername, filename*string(potnum), ui[1:(n-1)], xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
        stropu1 = modifystr2(stropu1, natomtype);
        strcpu1 = modifystr2(strcpu1, natomtype);
        strgpu1 = modifystr2(strgpu1, natomtype);
        stropu2, strcpu2, strgpu2 = gendensity(foldername, fn*string(potnum), ui[n], rho, mu, eta, kappa, gen, ifile);
        stropu = stropu * stropu1 * "\n" * stropu2 * "\n";
        strcpu = strcpu * strcpu1 * "\n" * strcpu2 * "\n";
        strgpu = strgpu * strgpu1 * "\n" * strgpu2 * "\n";
    end
end

if potnum==0
    ifile = 0;
    gen = 0;
    stropu, strcpu, strgpu = genpair(foldername, filename, u, xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
    tmpopu, tmpcpu, tmpgpu = gendensity(foldername, fn, u, mu, mu, eta, kappa, gen, ifile);
    stropu = stropu  * "\n" * tmpopu * "\n";
    strcpu = strcpu  * "\n" * tmpcpu * "\n";
    strgpu = strgpu  * "\n" * tmpgpu * "\n";
    cppfiles(foldername, filename, stropu, strcpu, strgpu, 1);
else
    sp0, sp1, sp2, sp3 = getpotstr(2);
    sp0 = replace(sp0, "int ng)" => "int ng, int potnum)");
    sp1 = replace(sp1, "int);" => "int, int);");
    sp2 = replace(sp2, "int);" => "int, int);");
    tmpopu, tmpcpu, tmpgpu = genpotnum(filename, potnum, sp0, sp1, sp2, sp3);
    stropu = stropu * tmpopu * "\n";
    strcpu = strcpu * tmpcpu * "\n";
    strgpu = strgpu * tmpgpu * "\n";

    # here
    tmpgpu = replace(tmpgpu, "(T *u," => "Gradient(T *u, T *du, T *u_xij,");
    tmpgpu = replace(tmpgpu, "(u," => "Gradient(u, du, u_xij,");
    tmpgpu = replace(tmpgpu, "(double *" => "Gradient(double *, double *, double*");
    tmpgpu = replace(tmpgpu, "(float *" => "Gradient(float *, float *, float*");
    strgpu = strgpu * "\n" * tmpgpu;

    sp0, sp1, sp2, sp3 = getpotstr(0);
    sp0 = replace(sp0, "int ng)" => "int ng, int potnum)");
    sp1 = replace(sp1, "int);" => "int, int);");
    sp2 = replace(sp2, "int);" => "int, int);");
    tmpopu, tmpcpu, tmpgpu = genpotnum(fn, potnum, sp0, sp1, sp2, sp3);
    stropu = stropu * tmpopu;
    strcpu = strcpu * tmpcpu;
    strgpu = strgpu * tmpgpu;

    # here
    tmpgpu = replace(tmpgpu, "(T *u," => "Gradient(T *u, T *du, T *u_rho,");
    tmpgpu = replace(tmpgpu, "(u," => "Gradient(u, du, u_rho,");
    tmpgpu = replace(tmpgpu, "(double *" => "Gradient(double *, double *, double*");
    tmpgpu = replace(tmpgpu, "(float *" => "Gradient(float *, float *, float*");
    strgpu = strgpu * "\n" * tmpgpu;

    cppfiles(foldername, filename, stropu, strcpu, strgpu, 0);
end

end

function modifystr2(oldstr, natomtype)

newstr = oldstr;
for i = 1:natomtype
    s = "u[" * string(i-1) * " * i*" * string(natomtype) * "]";
    if contains(newstr, s)
        news = "if (tj[i] == " * string(i) * ") \n\t\t\t" * "u[i]";
        newstr = replace(newstr, s => news);
    end
end

return newstr

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using DelimitedFiles

"""
    get_cvi_data(data_file::String)

Get the cvi data specified by the data_file path.
"""
function get_cvi_data(data_file::String)
    # Parse the data
    data = readdlm(data_file, ',')
    data = permutedims(data)
    train_x = data[1:2, :]
    train_y = convert(Array{Int64}, data[3, :])

    return train_x, train_y
end # get_cvi_data(data_file::String)

"""
    sort_cvi_data(data::Array{N, 2}, labels::Array{M, 1}) where {N<:Real, M<:Int}

Sorts the CVI data by the label index, assuring that clusters are provided incrementally.
"""
function sort_cvi_data(data::Array{N, 2}, labels::Array{M, 1}) where {N<:Real, M<:Int}
    index = sortperm(labels)
    data = data[:, index]
    labels = labels[index]

    return data, labels
end # sort_cvi_data(data::Array{N, 2}, labels::Array{M, 1}) where {N<:Real, M<:Int}

"""
    relabel_cvi_data(labels::Array{M, 1}) where {M<:Int}

Relabels the vector to present new labels in incremental order.
"""
function relabel_cvi_data(labels::Array{M, 1}) where {M<:Int}
    # Get the unique labels and their order of appearance
    unique_labels = unique(labels)
    n_unique_labels = length(unique_labels)
    n_labels = length(labels)

    # Create a new ordered list of unique labels
    new_unique_labels = [x for x in 1:n_unique_labels]

    # Map the old unique labels to the new ones
    label_mapping = Dict(zip(unique_labels, new_unique_labels))

    # Create a new labels vector with ordered labels
    new_labels = zeros(Int, n_labels)
    for ix = 1:n_labels
        new_labels[ix] = label_mapping[labels[ix]]
    end

    return new_labels
end # relabel_cvi_data(labels::Array{M, 1}) where {M<:Int}


text:
module PLMC

using StatsBase
using Distributions
using Clustering
using ClusterComplex
using Random
using ComputationalHomology
using LogExpFunctions: logsumexp
using Requires

import Clustering: assignments, counts, nclusters
import ClusterComplex: models

export PLMClusteringResult, modelclass,
       refinedmdl, ibdiff, mdldiff,
       Agglomeration,
       agglomerate, agglomerateph, testmerges, plmc,
       #re-export,
       MvNormal, MixtureModel

include("types.jl")
include("mdl.jl")
include("ib.jl")
include("spectral.jl")
include("plots.jl")
include("utils.jl")

# Agglomerative Clustering

function show_vector(io, v::AbstractVector)
    print(io, "[")
    for (i, elt) in enumerate(v)
        i > 1 && print(io, ", ")
        print(io, elt)
    end
    print(io, "]")
end

function testmerges(scorefunc, mtree::Vector{Vector{Vector{Int}}},
                    mcr::ModelClusteringResult, X::AbstractMatrix)
    c = length(mtree)
    MDL = fill(Inf, c)
    LL = fill(Inf, c)
    Rₘₐₓ = fill(Inf, c)
    for i in 1:c
        MDL[i], LL[i], Rₘₐₓ[i] = scorefunc(mtree[i], mcr, X)
        @debug begin
            io = IOBuffer()
            show_vector(io, mtree[i])
            strcls = String(take!(io))
            "$i -> $strcls\nMDL = $(MDL[i])\nLL = $(LL[i])\nRₘₐₓ = $(Rₘₐₓ[i])"
        end
    end
    return MDL, LL, Rₘₐₓ
end
"""
    agglomerate(flt::Filtration) -> Agglomeration

Construct agglomerative pairwise merge tree from the filtration `flt`.
"""
function agglomerate(flt::Filtration; individual=true, kwargs...)
    mtree   = Vector{Vector{Int}}[]
    mergers = Vector{Vector{Int}}[]
    mvals   = valtype(flt)[]

    # use individual clusters as PLM clusters
    if individual
        insert!(mtree, 1, [[i] for i in 1:size(complex(flt),0)])
        push!(mvals, zero(eltype(mvals)))
    end
    cplx = complex(flt)

    ST = Vector{Int}[]
    for (v, splxs) in flt
        merged = false
        for (d,sid) in splxs
            simplex = cplx[sid, d]
            # join simplex points
            # println("$v => $simplex")
            pts = values(simplex)
            if length(pts) == 1
                # println("[] <= $(pts[1])")
                push!(ST, pts)
            else
                j = pts[1]
                jj = findfirst(a-> j ∈ a, ST)
                for k in pts[2:end]
                    kk = findfirst(a-> k ∈ a, ST)
                    # println("$j <= $k : $jj <= $kk")
                    jj == kk && continue
                    # println("$jj <= $kk : $(ST[jj]) <= $(ST[kk])")
                    push!(mergers, [deepcopy(ST[jj]), deepcopy(ST[kk])])
                    append!(ST[jj], ST[kk])
                    deleteat!(ST, kk)
                    merged = true
                end
            end
        end
        if merged
            push!(mtree, deepcopy(ST))
            push!(mvals, v)
            # println("$(mtree[end]) => $v")
        end
    end
    return Agglomeration(mtree, mergers, mvals)
end

"""
    agglomerate(flt::Filtration, mcr::ModelClusteringResult, X::AbstractMatrix) -> (PLMClusteringResult, Agglomeration)

Return agglomerative piecewise clustering from the filtration `flt` by finding a corrspondent minimum MDL value.
"""
function agglomerate(flt::Filtration, mcr::ModelClusteringResult,
                     X::AbstractMatrix; individual=true)
    agg = agglomerate(flt, individual=individual)
    return plmc(agg, mcr, X, filtration=flt), agg
end

function agglomerateph(flt::Filtration; individual=true)
    mtree   = Vector{Vector{Int}}[]
    mergers = Vector{Vector{Int}}[]
    mvals   = valtype(flt)[]

    # use individual clusters as PLM clusters
    if individual
        insert!(mtree, 1, [[i] for i in 1:size(complex(flt),0)])
        push!(mvals, zero(eltype(mvals)))
    end

    cplx = complex(flt)
    ph = persistenthomology(PersistentCocycleReduction{Float64}, flt)
    dgm = diagram(ph)
    dims = sort!(collect(keys(dgm)))
    fvals = values(ph)

    ST = Vector{Int}[]
    # i = 16
    for i in 1:length(fvals)-1
        cv, nv = fvals[i:i+1]
        itr = collect(Iterators.flatten(filter(x-> cv < death(x) <=nv, dgm[d]) for d in dims))
        gs = ComputationalHomology.generator.(itr)
        # g = gs[2]
        for g in gs
            ComputationalHomology.dim(g) > 0 && continue

            splxs = unique(Iterators.flatten(values(cplx[sid, 0]) for sid in keys(g)))
            mrg = deepcopy(mtree[end])
            mrgidx = map(c->any(s-> s ∈ c, splxs), mrg)
            push!(mergers, mrg[mrgidx])
            merged = vcat(mrg[mrgidx]...)
            deleteat!(mrg, mrgidx)
            append!(mrg, [merged])

            push!(mtree, mrg)
            push!(mvals, cv)
            @debug "Merging" at=mvals[end] state=mtree[end]' generator=g
        end
    end
    return Agglomeration(mtree, mergers, mvals)
end

"""
    agglomerate(flt::Filtration, mcr::ModelClusteringResult, X::AbstractMatrix) -> (PLMClusteringResult, Agglomeration)

Return agglomerative piecewise clustering from the filtration `flt` by finding a corrspondent minimum MDL value.
"""
function agglomerateph(flt::Filtration, mcr::ModelClusteringResult,
                       X::AbstractMatrix; individual=true)
    agg = agglomerateph(flt, individual=individual)
    return plmc(agg, mcr, X, filtration=flt), agg
end

function mergecost(nodes::Vector{Vector{Int}}, measure::Function,
                   mcr::ModelClusteringResult, X::AbstractMatrix, β::Int;
                   kwargs...)
    logps = hcat((logpdf(p, X) for p in models(mcr))...)
    assign = assignments(mcr)
    return mergecost(nodes, measure, logps, assign, β; kwargs...)
end

function mergecost(nodes::Vector{Vector{Int}}, measure::Function,
                   logps::AbstractMatrix{T}, assign::Vector{Int}, β::Int;
                   kwargs...) where {T <: AbstractFloat}
    n = length(nodes)
    # println("nodes: ", n)
    ds = [let nds = nodes[[i,j]]; nds=>measure(nds, logps, assign; kwargs...); end for i in 1:n, j in 1:n if i < j]
    dsi = sortperm(ds, by=last)
    mi = min(β, length(dsi))
    return ds[dsi[1:mi]]
end

"""
    agglomerate(measure::Function, mcr::ModelClusteringResult, X::AbstractMatrix; β=1) -> Agglomeration

Construct agglomerative pairwise clustering through the beam search using `measure` function. `β` is a length of the beam.
"""
function agglomerate(measure::Function, mcr::ModelClusteringResult,
                     X::AbstractMatrix{T}; kwargs...) where {T <: AbstractFloat}
    k = nclusters(mcr)
    args = Dict(kwargs...)
    β = get(args, :β, 1)
    βsize = get(args, :βsize, β)

    assign = assignments(mcr)
    priorC = map(c->c/length(assign), counts(mcr))
    logps = hcat((logpdf(p, X) for p in models(mcr))...)
    logpx = sum(exp.(logps' .+ log.(priorC)), dims=1)' |> vec

    # initialize beam storage
    beam = Agglomeration[]
    push!(beam, Agglomeration([[i] for i in 1:k]))
    bcount = β
    cls = last(first(beam))
    terminate = false
    while length(cls) > 1 && !terminate
        # fill beam structure with possible options
        while bcount > 0
            tmp = Agglomeration[]
            while length(beam) > 0
                agg = pop!(beam)
                mcosts = mergecost(last(agg), measure, logps, assign, β; marginal=logpx, kwargs...)
                for mrg in mcosts
                    @debug "Merging" at=mrg aggregation=agg
                    if isinf(mrg.second)
                        push!(tmp, agg)
                        terminate = true
                        break
                    else
                        newagg = deepcopy(agg)
                        push!(newagg, mrg)
                        push!(tmp, newagg)
                    end
                end
            end
            unique!(a->last(a), tmp)
            beam = tmp
            bcount -= 1
        end
        bcount = 1

        # prune beam structure
        idxs = sortperm(map(count, beam))
        mi = min(β^βsize, length(idxs))
        beam = beam[idxs[1:mi]]
        # map(a->a.clusters[end] => count(a), beam)
        cls = last(first(beam))
    end
    return first(beam)
end

"""
    plmc(agg::Agglomeration, mcr::ModelClusteringResult, X::AbstractMatrix) -> PLMClusteringResult

Construct a piecewise clustering from the agglomerative merge tree `agg` by finding a corrspondent minimum MDL value.
"""
function plmc(agg::Agglomeration, mcr::ModelClusteringResult,
              X::AbstractMatrix, filtration=nothing; kwargs...)
    args = Dict(kwargs...)
    scorefunc = get(args, :score, refinedmdl)
    findscore = get(args, :find, findglobalmin)
    mdlidx = get(args, :aggidx, 0)

    scores = nothing
    if mdlidx == 0
        scores, _ = testmerges(scorefunc, agg.clusters, mcr, X)
        mdlidx = findscore(scores)
    end

    ϵ = Inf
    scplx = SimplicialComplex()
    if filtration !== nothing
        ϵ = agg.costs[mdlidx]
        scplx = complex(filtration, ϵ)
    end
    return PLMClusteringResult(mcr, agg.clusters[mdlidx], scplx, ϵ), agg, scores
end

"""
    plmc(AlgorithmType, flt::Filtration, mcr::ModelClusteringResult, X::AbstractMatrix) -> PLMClusteringResult

Construct a piecewise clustering of the dataset `X` using `AlgorithmType` algorithm given the cover represented by a filtration `flt` and a mode-based clustering `mcr`.
"""
function plmc(::Type{Topological}, X::AbstractMatrix, flt::Filtration,
              mcr::ModelClusteringResult; kwargs...)
    agg = agglomerate(flt; kwargs...)
    plmc(agg, mcr, X, flt; kwargs...)
end

function plmc(::Type{PHomology}, X::AbstractMatrix, flt::Filtration,
              mcr::ModelClusteringResult; kwargs...)
    agg = agglomerateph(flt; kwargs...)
    plmc(agg, mcr, X, flt; kwargs...)
end

function plmc(::Type{MDL}, X::AbstractMatrix, flt::Filtration,
              mcr::ModelClusteringResult; kwargs...)
    agg = agglomerate(mdldiff, mcr, X; kwargs...)
    plmc(agg, mcr, X; kwargs...)
end

function plmc(::Type{InformationBottleneck}, X::AbstractMatrix, flt::Filtration,
              mcr::ModelClusteringResult; kwargs...)
    agg = agglomerate(ibdiff, mcr, X; kwargs...)
    plmc(agg, mcr, X; kwargs...)
end

function __init__()
    @require Optim="429524aa-4258-5aef-a3af-852621145aeb" @eval using Optim: optimize
end

end

max_stars_repo_name: wildart/PLMC

max_issues_repo_name: wildart/PLMC

max_forks_repo_name: wildart/PLMC


text:
# Lebesgue measure

export Lebesgue

struct Lebesgue{X} <: PrimitiveMeasure end

function Pretty.tile(::Lebesgue{X}) where {X}
    result = Pretty.literal("Lebesgue(")
    result *= Pretty.tile(X)
    result *= Pretty.literal(")")
end

Lebesgue(X) = Lebesgue{X}()

sampletype(::Lebesgue{ℝ}) = Float64
sampletype(::Lebesgue{ℝ₊}) = Float64
sampletype(::Lebesgue{𝕀}) = Float64

testvalue(::Lebesgue{ℝ}) = 0.0
testvalue(::Lebesgue{𝕀}) = 0.5
testvalue(::Lebesgue{ℝ₊}) = 1.0
testvalue(::Lebesgue{<:Real}) = 0.0

logdensity(::Lebesgue, x) = zero(x)

Base.:∘(::typeof(basemeasure), ::Type{Lebesgue}) = Lebesgue


text:
# THIS FILE TAKEN FROM GraphLayout AND EDITED SO IT WORKS

"""
Compute graph layout using stress majorization

Inputs:

    δ: Matrix of pairwise distances
    p: Dimension of embedding (default: 2)
    w: Matrix of weights. If not specified, defaults to
           w[i,j] = δ[i,j]^-2 if δ[i,j] is nonzero, or 0 otherwise
    X0: Initial guess for the layout. Coordinates are given in rows.
        If not specified, default to random matrix of Gaussians

Additional optional keyword arguments control the convergence of the algorithm
and the additional output as requested:

    maxiter:   Maximum number of iterations. Default: 400size(X0, 1)^2
    abstols:   Absolute tolerance for convergence of stress.
               The iterations terminate if the difference between two
               successive stresses is less than abstol.
               Default: √(eps(eltype(X0))
    reltols:   Relative tolerance for convergence of stress.
               The iterations terminate if the difference between two
               successive stresses relative to the current stress is less than
               reltol. Default: √(eps(eltype(X0))
    abstolx:   Absolute tolerance for convergence of layout.
               The iterations terminate if the Frobenius norm of two successive
               layouts is less than abstolx. Default: √(eps(eltype(X0))
    verbose:   If true, prints convergence information at each iteration.
               Default: false
    returnall: If true, returns all iterates and their associated stresses.
               If false (default), returns the last iterate

Output:

    The final layout X, with coordinates given in rows, unless returnall=true.

Reference:

    The main equation to solve is (8) of:

    @incollection{
        author = {Emden R Gansner and Yehuda Koren and Stephen North},
        title = {Graph Drawing by Stress Majorization}
        year={2005},
        isbn={978-3-540-24528-5},
        booktitle={Graph Drawing},
        seriesvolume={3383},
        series={Lecture Notes in Computer Science},
        editor={Pach, J\'anos},
        doi={10.1007/978-3-540-31843-9_25},
        publisher={Springer Berlin Heidelberg},
        pages={239--250},
    }
"""
function my_layout_stressmajorize_adj(δ, p::Int=2, w=nothing, X0=randn(size(δ, 1), p);
        maxiter = 400size(X0, 1)^2, abstols=√(eps(eltype(X0))),
        reltols=√(eps(eltype(X0))), abstolx=√(eps(eltype(X0))),
        verbose = false, returnall = false)

    @assert size(X0, 2)==p

    if w==nothing
        w = δ.^-2
        w[.!isfinite.(w)] .= 0
    end

    @assert size(X0, 1)==size(δ, 1)==size(δ, 2)==size(w, 1)==size(w, 2)
    Lw = weightedlaplacian(w)
    pinvLw = pinv(Lw)
    newstress = stress(X0, δ, w)
    Xs = Matrix[X0]
    stresses = [newstress]
    itcounter = 0
    for iter = 1:maxiter
        itcounter = iter
        #TODO the faster way is to drop the first row and col from the iteration
        X = pinvLw * (LZ(X0, δ, w)*X0)
        @assert all(isfinite.(X))
        newstress, oldstress = stress(X, δ, w), newstress
        verbose && info("""Iteration $iter
        Change in coordinates: $(vecnorm(X - X0))
        Stress: $newstress (change: $(newstress-oldstress))
        """)
        push!(Xs, X)
        push!(stresses, newstress)
        abs(newstress - oldstress) < reltols * newstress && break
        abs(newstress - oldstress) < abstols && break
        norm(X - X0) < abstolx && break
        X0 = X
    end
    itcounter == maxiter && warn("Maximum number of iterations reached without convergence")
    returnall ? (Xs, stresses) : Xs[end]
end

"""
Stress function to majorize

Input:
    X: A particular layout (coordinates in rows)
    d: Matrix of pairwise distances
    w: Weights for each pairwise distance

See (1) of Reference
"""
function stress(X, d=fill(1.0, size(X, 1), size(X, 1)), w=nothing)
    s = 0.0
    n = size(X, 1)
    if w==nothing
        w = d.^-2
        w[.!isfinite.(w)] = 0
    end
    @assert n==size(d, 1)==size(d, 2)==size(w, 1)==size(w, 2)
    for j=1:n, i=1:j-1
        s += w[i, j] * (norm(X[i,:] - X[j,:]) - d[i,j])^2
    end
    @assert isfinite.(s)
    s
end

function _checksquare(M::Matrix)
    r,c = size(M)
    return r
end
"""
Compute weighted Laplacian given ideal weights w

Lʷ defined in (4) of the Reference
"""
function weightedlaplacian(w)
    n = _checksquare(w)
    T = eltype(w)
    Lw = zeros(T, n, n)
    for i=1:n
        D = zero(T)
        for j=1:n
            i==j && continue
            Lw[i, j] = -w[i, j]
            D += w[i, j]
        end
        Lw[i, i]=D
    end
    Lw
end

"""
Computes L^Z defined in (5) of the Reference

Input: Z: current layout (coordinates)
       d: Ideal distances (default: all 1)
       w: weights (default: d.^-2)
"""
function LZ(Z, d, w)
    n = size(Z, 1)
    L = zeros(n, n)
    for i=1:n
        D = 0.0
        for j=1:n
            i==j && continue
            nrmz = norm(Z[i,:] - Z[j,:])
            nrmz==0 && continue
            δ = w[i, j] * d[i, j]
            L[i, j] = -δ/nrmz
            D -= -δ/nrmz
        end
        L[i, i] = D
    end
    @assert all(isfinite.(L))
    L
end


text:
module DiffEqParamEstim
using DiffEqBase, LsqFit, PenaltyFunctions,
      RecursiveArrayTools, ForwardDiff, Calculus, Distributions,
      LinearAlgebra, DiffEqSensitivity, Dierckx,
      SciMLBase

import PreallocationTools
STANDARD_PROB_GENERATOR(prob,p) = remake(prob;u0=eltype(p).(prob.u0),p=p)
STANDARD_PROB_GENERATOR(prob::EnsembleProblem,p) = EnsembleProblem(
                             remake(prob.prob;u0=eltype(p).(prob.prob.u0),p=p),
                             output_func = prob.output_func,
                             prob_func = prob.prob_func,
                             reduction = prob.reduction,
                             u_init = prob.u_init)
STANDARD_MS_PROB_GENERATOR = function(prob, p, k)
    t0, tf = prob.tspan
    P, N = length(prob.p), length(prob.u0)
    K = Int( (length(p) - P) / N)
    τ = range(t0, tf, length=K+1)
    remake( prob; u0=p[(1+(k-1)*N):(k*N)], p=p[end-P+1:end],
            tspan=(τ[k], τ[k+1]))
end

include("cost_functions.jl")
include("lm_fit.jl")
include("build_loss_objective.jl")
include("build_lsoptim_objective.jl")
include("kernels.jl")
include("two_stage_method.jl")
include("multiple_shooting_objective.jl")

end # module


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct ReactivePowerPI <: ReactivePowerControl
        Kp_q::Float64
        Ki_q::Float64
        ωf::Float64
        V_ref::Float64
        Q_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of a Proportional-Integral Reactive Power controller for a specified power reference

# Arguments
- `Kp_q::Float64`: Proportional Gain, validation range: `(0, nothing)`
- `Ki_q::Float64`: Integral Gain, validation range: `(0, nothing)`
- `ωf::Float64`: filter frequency cutoff, validation range: `(0, nothing)`
- `V_ref::Float64`: Voltage Set-point, validation range: `(0, nothing)`
- `Q_ref::Float64`: Reactive Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the ReactivePowerPI model are:
	σq_oc: Integrator state of the PI Controller,
	q_oc: Measured reactive power of the inverter model
- `n_states::Int`: ReactivePowerPI has two states
"""
mutable struct ReactivePowerPI <: ReactivePowerControl
    "Proportional Gain"
    Kp_q::Float64
    "Integral Gain"
    Ki_q::Float64
    "filter frequency cutoff"
    ωf::Float64
    "Voltage Set-point"
    V_ref::Float64
    "Reactive Power Set-point"
    Q_ref::Float64
    ext::Dict{String, Any}
    "The states of the ReactivePowerPI model are:
	σq_oc: Integrator state of the PI Controller,
	q_oc: Measured reactive power of the inverter model"
    states::Vector{Symbol}
    "ReactivePowerPI has two states"
    n_states::Int
end

function ReactivePowerPI(Kp_q, Ki_q, ωf, V_ref=1.0, Q_ref=1.0, ext=Dict{String, Any}(), )
    ReactivePowerPI(Kp_q, Ki_q, ωf, V_ref, Q_ref, ext, [:σq_oc, :q_oc], 2, )
end

function ReactivePowerPI(; Kp_q, Ki_q, ωf, V_ref=1.0, Q_ref=1.0, ext=Dict{String, Any}(), states=[:σq_oc, :q_oc], n_states=2, )
    ReactivePowerPI(Kp_q, Ki_q, ωf, V_ref, Q_ref, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function ReactivePowerPI(::Nothing)
    ReactivePowerPI(;
        Kp_q=0,
        Ki_q=0,
        ωf=0,
        V_ref=0,
        Q_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`ReactivePowerPI`](@ref) `Kp_q`."""
get_Kp_q(value::ReactivePowerPI) = value.Kp_q
"""Get [`ReactivePowerPI`](@ref) `Ki_q`."""
get_Ki_q(value::ReactivePowerPI) = value.Ki_q
"""Get [`ReactivePowerPI`](@ref) `ωf`."""
get_ωf(value::ReactivePowerPI) = value.ωf
"""Get [`ReactivePowerPI`](@ref) `V_ref`."""
get_V_ref(value::ReactivePowerPI) = value.V_ref
"""Get [`ReactivePowerPI`](@ref) `Q_ref`."""
get_Q_ref(value::ReactivePowerPI) = value.Q_ref
"""Get [`ReactivePowerPI`](@ref) `ext`."""
get_ext(value::ReactivePowerPI) = value.ext
"""Get [`ReactivePowerPI`](@ref) `states`."""
get_states(value::ReactivePowerPI) = value.states
"""Get [`ReactivePowerPI`](@ref) `n_states`."""
get_n_states(value::ReactivePowerPI) = value.n_states

"""Set [`ReactivePowerPI`](@ref) `Kp_q`."""
set_Kp_q!(value::ReactivePowerPI, val) = value.Kp_q = val
"""Set [`ReactivePowerPI`](@ref) `Ki_q`."""
set_Ki_q!(value::ReactivePowerPI, val) = value.Ki_q = val
"""Set [`ReactivePowerPI`](@ref) `ωf`."""
set_ωf!(value::ReactivePowerPI, val) = value.ωf = val
"""Set [`ReactivePowerPI`](@ref) `V_ref`."""
set_V_ref!(value::ReactivePowerPI, val) = value.V_ref = val
"""Set [`ReactivePowerPI`](@ref) `Q_ref`."""
set_Q_ref!(value::ReactivePowerPI, val) = value.Q_ref = val
"""Set [`ReactivePowerPI`](@ref) `ext`."""
set_ext!(value::ReactivePowerPI, val) = value.ext = val


text:
"""
    similar_rxns(model::MetNet, iders = eachindex(model.rxns); verbose = true)

Findall the reaction that transforms the same elements
if we define: r1 = 2A -> 2B, r2 = B <-> A, r3 = A <- B, r4 = A + B ->
the method will returns [[r1, r2, r3]]
This method only use the information in the binary S.
Bounds are not relevant either.

"""
function similar_rxns(model::MetNet,
        iders = eachindex(model.rxns); verbose = true)

    idxs = [rxnindex(model, ider) for ider in iders]
    # collecting react and prods hashs
    # (prods_hash, react_hash) => reaction idxs
    # (react_hash, prods_hash) => reaction idxs
    hash_table = Dict{Tuple{UInt64,UInt64}, Vector{Int}}()
    for rxni in idxs
        rhash = hash(Set(rxn_reacts(model, rxni)))
        phash = hash(Set(rxn_prods(model, rxni)))
        rxnis = get!(hash_table, (rhash, phash), Int[])
        push!(rxnis, rxni)
    end

    # Form pairs
    similars = Vector{Vector{Int}}()
    for (fkey, frxnis) in hash_table

        bkey = fkey |> reverse
        brxnis = get(hash_table, bkey, Int[])

        sims_ = union(frxnis, brxnis)
        length(sims_) > 1 && push!(similars, sims_)

        # deleting
        delete!(hash_table, fkey)
        delete!(hash_table, bkey)

    end

    if verbose
        for rxnis in similars
            println("Similars ", length(rxnis))
            for rxni in rxnis
                println(model.rxns[rxni], " [", rxni, "]: ", rxn_str(model, rxni))
            end
            println()
        end
    end
    return similars
end


text:
# Example of Jacobian usage for HINT layer (Kruse et al, 2020)
# Author: Gabrio Rizzuti, grizzuti3@gatech.edu
# Date: October 2020

using InvertibleNetworks, LinearAlgebra, Test

# Initialize HINT layer
nx = 32
ny = 32
n_channel = 32
n_hidden = 64
batchsize = 2
logdet = true
N = CouplingLayerHINT(n_channel, n_hidden; logdet=logdet, permute="full")
θ = get_params(N)

## f: X × θ ↦ Y, ∂f: ΔX × Δθ ↦ ΔY

# Compute Jacobian
X = randn(Float32, nx, ny, n_channel, batchsize)
J = Jacobian(N, X; io_mode="X×θ↦Y")

# Evaluate Jacobian
ΔX = randn(Float32, size(X))
Δθ = Array{Parameter, 1}(undef, length(θ)); for i = 1:length(θ) Δθ[i] = Parameter(randn(Float32, size(θ[i]))); end;
ΔY = J*(ΔX, Δθ)

# Evaluate adjoint Jacobian
JT = adjoint(J)
ΔY_ = randn(Float32, size(X))
ΔX_, Δθ_ =  JT*ΔY_

# Adjoint Test
a = dot(ΔX, ΔX_)+dot(Δθ, Δθ_)
b = dot(ΔY, ΔY_)
@test isapprox(a, b; rtol=1f-3)

## f(x; ⋅) : θ ↦ Y, ∂f(x): Δθ ↦ ΔY

# Compute Jacobian
X = randn(Float32, nx, ny, n_channel, batchsize)
J = Jacobian(N, X; io_mode="θ↦Y")

# Evaluate Jacobian
ΔX = randn(Float32, size(X))
Δθ = Array{Parameter, 1}(undef, length(θ)); for i = 1:length(θ) Δθ[i] = Parameter(randn(Float32, size(θ[i]))); end;
ΔY = J*Δθ

# Evaluate adjoint Jacobian
JT = adjoint(J)
ΔY_ = randn(Float32, size(X))
Δθ_ =  JT*ΔY_

# Adjoint Test
a = dot(Δθ, Δθ_)
b = dot(ΔY, ΔY_)
@test isapprox(a, b; rtol=1f-3)


text:
export Segment, SegmentParameters

struct SigmaRho
  sigma :: Float64
  rho   :: Float64
  function SigmaRho(sigma,rho)
    if ~isnan(rho) && ~isnan(sigma)
      throw(ArgumentError("Cannot specify both rho and sigma"))
    end
    new(sigma,rho)
  end
end
SigmaRho(;sigma = NaN, rho=NaN) = SigmaRho(sigma,rho)
function Base.show(io::IO, x::SigmaRho)
  if isnan(x.sigma) && isnan(x.rho)
    return nothing
  end
  print(io,"+ ")
  if ~isnan(x.sigma)
    @printf(io," sigma=%.9e",x.sigma)
  end
  if ~isnan(x.rho)
    @printf(io," rho=%.9e",x.rho)
  end
  println(io)
  return nothing
end

mutable struct WxWyWz
  xyz :: Array{Float64,1}
  isdefault :: Bool
  function WxWyWz(wx, wy, wz)
    if norm([wx,wy,wz])<1e-9
      throw(ArgumentError("norm([wx,wy,wz])<1e-9"))
    end
    new([wx,wy,wz],isnan(wx) && isnan(wy) && isnan(wz))
  end
end
WxWyWz(;wx=NaN, wy=NaN, wz=NaN) =
  WxWyWz(wx,wy,wz)

function Base.show(io::IO, x::WxWyWz)
  if ~x.isdefault
    print(io,"+ ")
    if ~isnan(x.xyz[1])
      @printf(io," wx=%.9e",x.xyz[1])
    end
    if ~isnan(x.xyz[2])
      @printf(io," wy=%.9e",x.xyz[2])
    end
    if ~isnan(x.xyz[3])
      @printf(io," wz=%.9e",x.xyz[3])
    end
    println(io)
  end
  return nothing
end

struct WH
  w :: Float64
  h :: Float64
  nhinc :: Int
  nwinc :: Int
  rh :: Float64
  rw :: Float64
  function WH(w, h, nhinc, nwinc, rh, rw)
    if nhinc<0
      throw(ArgumentError("nhinc must be a positive integer"))
    end
    if nwinc<0
      throw(ArgumentError("nwinc must be a positive integer"))
    end
    if rh<0
      throw(ArgumentError("rh must be positive"))
    end
    if rw<0
      throw(ArgumentError("rw must be positive"))
    end
    return new(w, h, nhinc, nwinc, rh, rw)
  end
end
WH(;w=NaN, h=NaN, nhinc=0, nwinc=0, rh=NaN, rw=NaN) =
  WH(w, h, nhinc, nwinc, rh, rw)
function Base.show(io::IO, x::WH)
  if isnan(x.w) && isnan(x.h) && x.nhinc==0 && x.nwinc==0 && isnan(x.rh) && isnan(x.rw)
    return nothing
  end
  print(io,"+ ")
  if ~isnan(x.w)
    @printf(io," w=%.9e",x.w)
  end
  if ~isnan(x.h)
    @printf(io," h=%.9e",x.h)
  end
  if x.nhinc>0
    print(io," nhinc=",x.nhinc)
  end
  if x.nwinc>0
    print(io," nwinc=",x.nwinc)
  end
  if ~isnan(x.rh)
    @printf(io," rh=%.9e",x.rh)
  end
  if ~isnan(x.rw)
    @printf(io," rw=%.9e",x.rw)
  end
  println(io)
  return nothing
end

"""
    SegmentParameters(<keyword arguments>)
    SegmentParameters(parameters::SegmentParameters, <keyword arguments>)

Object to hold parameters for `Segment` and `Default`

**Keyword Arguments**

- `w`               -- segment width
- `h`               -- segment height
- `sigma`           -- conductivity
- `rho`             -- resistivity
- `wx`, `wy`, `wz`  -- segment orientation.  vector pointing along width of segment's cross section.
- `nhinc`, `nwinc`  -- integer number of filaments in height and width
- `rh`, `rw`        -- ratio in the height and width

The second form replaces values in `parameters` with any keyword parameters present.
"""
struct SegmentParameters
  wh        :: WH
  sigmarho  :: SigmaRho
  wxwywz    :: WxWyWz
end
SegmentParameters(;w=NaN, h=NaN, sigma=NaN, rho=NaN, wx=NaN, wy=NaN, wz=NaN,
                  nhinc=0, nwinc=0, rh=NaN, rw=NaN) =
  SegmentParameters(WH(w,h,nhinc,nwinc,rh,rw), SigmaRho(sigma,rho), WxWyWz(wx,wy,wz))
function SegmentParameters(sp::SegmentParameters; w=NaN, h=NaN, sigma=NaN, rho=NaN, wx=NaN,
                           wy=NaN, wz=NaN, nhinc=0, nwinc=0, rh=NaN, rw=NaN)
  new_w = isnan(w) ? sp.wh.w : w
  new_h = isnan(h) ? sp.wh.h : h
  new_nhinc = nhinc==0 ? sp.wh.nhinc : nhinc
  new_nwinc = nwinc==0 ? sp.wh.nwinc : nwinc
  new_rh = isnan(rh) ? sp.wh.rh : rh
  new_rw = isnan(rw) ? sp.wh.rw : rw
  new_sigma = isnan(sigma) ? sp.sigmarho.sigma : sigma
  new_rho = isnan(rho) ? sp.sigmarho.rho : rho
  new_wx = isnan(wx) ? sp.wxwywz.xyz[1] : wx
  new_wy = isnan(wy) ? sp.wxwywz.xyz[2] : wy
  new_wz = isnan(wz) ? sp.wxwywz.xyz[3] : wz
  new_wh = WH(new_w,new_h,new_nhinc,new_nwinc,new_rh,new_rw)
  new_sigmarho = SigmaRho(new_sigma, new_rho)
  new_wxwywz = WxWyWz(new_wx,new_wy,new_wz)
  SegmentParameters(new_wh, new_sigmarho, new_wxwywz)
end

w(sp::SegmentParameters) = sp.wh.w
h(sp::SegmentParameters) = sp.wh.h
nhinc(sp::SegmentParameters) = sp.wh.nhinc
nwinc(sp::SegmentParameters) = sp.wh.nwinc
rh(sp::SegmentParameters) = sp.wh.rh
rw(sp::SegmentParameters) = sp.wh.rw
sigma(sp::SegmentParameters) = sp.sigmarho.sigma
rho(sp::SegmentParameters) = sp.sigmarho.rho
wx(sp::SegmentParameters) = sp.wxwywz.xyz[1]
wy(sp::SegmentParameters) = sp.wxwywz.xyz[2]
wz(sp::SegmentParameters) = sp.wxwywz.xyz[3]
wxyz(sp::SegmentParameters) = sp.wxwywz.xyz

"""
    Segment([name], n1::Node, n2::Node, [parameters::SegmentParameters])
    Segment([name], n1::Node, n2::Node, <keyword arguments>)

`Segment` objects `show` a FastHenry segment command.

Keyword arguments are the same as for `SegmentParameters`.

note: When rotating segments, the vector [wx, wy, wz] will be rotated.  Default
values will be used if not specified.  wx, wy, wz will show after `transform[!]` is
called on a `Segment`.
"""
struct Segment <: Element
  name      :: Symbol
  node1     :: Node
  node2     :: Node
  wh        :: WH
  sigmarho  :: SigmaRho
  wxwywz    :: WxWyWz
  function Segment(n, n1::Node, n2::Node, wh::WH, sigmarho::SigmaRho, wxwywz::WxWyWz)
    new_wxwywz = deepcopy(wxwywz)
#    initialixe_wxwywz!(n1,n2,new_wxwywz)
    new(Symbol(n), n1, n2, wh, sigmarho, new_wxwywz)
  end
end
Segment(n1::Node, n2::Node; w=NaN, h=NaN, sigma=NaN, rho=NaN,
        wx=NaN, wy=NaN, wz=NaN, nhinc=0, nwinc=0, rh=NaN, rw=NaN) =
            Segment(Symbol(""), n1, n2, WH(w,h,nhinc,nwinc,rh,rw), SigmaRho(sigma,rho), WxWyWz(wx,wy,wz))
Segment(name, n1::Node, n2::Node; w=NaN, h=NaN, sigma=NaN, rho=NaN,
        wx=NaN, wy=NaN, wz=NaN, nhinc=0, nwinc=0, rh=NaN, rw=NaN) =
            Segment(name, n1, n2, WH(w,h,nhinc,nwinc,rh,rw), SigmaRho(sigma,rho), WxWyWz(wx,wy,wz))
Segment(n1::Node, n2::Node, sp::SegmentParameters) =
            Segment(Symbol(""), n1, n2, sp.wh, sp.sigmarho, sp.wxwywz)
Segment(name, n1::Node, n2::Node, sp::SegmentParameters) =
            Segment(name, n1, n2, sp.wh, sp.sigmarho, sp.wxwywz)


text:
using Revise
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ImplicitFreeSurface
using Statistics
using Printf
using LinearAlgebra, SparseArrays
using Oceananigans.Solvers: constructors, unpack_constructors

function geostrophic_adjustment_simulation(free_surface, topology; arch = Oceananigans.CPU())

    Lh = 100kilometers
    Lz = 400meters

    grid = RectilinearGrid(arch,
        size = (64, 3, 1),
        x = (0, Lh), y = (0, Lh), z = (-Lz, 0),
        topology = topology)

    coriolis = FPlane(f = 1e-4)

    model = HydrostaticFreeSurfaceModel(grid = grid,
        coriolis = coriolis,
        free_surface = free_surface)

    gaussian(x, L) = exp(-x^2 / 2L^2)

    U = 0.1 # geostrophic velocity
    L = grid.Lx / 40 # gaussian width
    x₀ = grid.Lx / 4 # gaussian center

    vᵍ(x, y, z) = -U * (x - x₀) / L * gaussian(x - x₀, L)

    g = model.free_surface.gravitational_acceleration
    η = model.free_surface.η

    η₀ = coriolis.f * U * L / g # geostrohpic free surface amplitude

    ηᵍ(x) = η₀ * gaussian(x - x₀, L)

    ηⁱ(x, y) = 2 * ηᵍ(x)

    set!(model, v = vᵍ)
    set!(η, ηⁱ)

    gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed
    wave_propagation_time_scale = model.grid.Δxᶜᵃᵃ / gravity_wave_speed
    simulation = Simulation(model, Δt = 2wave_propagation_time_scale, stop_iteration = 300)

    return simulation
end

function run_and_analyze(simulation)
    η = simulation.model.free_surface.η
    u, v, w = simulation.model.velocities
    Δt = simulation.Δt

    ηx = Field(∂x(η))
    compute!(ηx)

    f = simulation.model.free_surface
    @views u₀ = interior(u)[:, 1, 1]
    @views v₀ = interior(v)[:, 1, 1]
    @views η₀ = interior(η)[:, 1, 1]
    @views ηx₀ = interior(ηx)[:, 1, 1]

    if f isa SplitExplicitFreeSurface
        solver_method = "SplitExplicitFreeSurface"
    else
        solver_method = string(simulation.model.free_surface.solver_method)
    end

    simulation.output_writers[:fields] = JLD2OutputWriter(simulation.model, (η, ηx, u, v, w),
        schedule = TimeInterval(Δt),
        prefix = "solution_$(solver_method)",
        force = true)

    progress_message(sim) = @info @sprintf("[%.2f%%], iteration: %d, time: %.3f, max|w|: %.2e",
        100 * sim.model.clock.time / sim.stop_time, sim.model.clock.iteration,
        sim.model.clock.time, maximum(abs, sim.model.velocities.u))

    simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(10))

    run!(simulation)

    compute!(ηx)

    @views u₁ = interior(u)[:, 1, 1]
    @views v₁ = interior(v)[:, 1, 1]
    @views η₁ = interior(η)[:, 1, 1]
    @views ηx₁ = interior(ηx)[:, 1, 1]

    @show mean(Array(η₀))
    @show mean(Array(η₁))

    Δη = Array(η₁) - Array(η₀)

    return (; η₀, η₁, Δη, ηx₀, ηx₁, u₀, u₁, v₀, v₁)
end

# fft_based_free_surface = ImplicitFreeSurface()
pcg_free_surface = ImplicitFreeSurface(solver_method = :PreconditionedConjugateGradient);
matrix_free_surface = ImplicitFreeSurface(solver_method = :HeptadiagonalIterativeSolver);
splitexplicit_free_surface = SplitExplicitFreeSurface()

topology_types = [(Bounded, Periodic, Bounded), (Periodic, Periodic, Bounded)]
topology_types = [topology_types[1]]

archs = [Oceananigans.CPU(), Oceananigans.GPU()]
archs = [archs[1]]

free_surfaces = [pcg_free_surface, matrix_free_surface, splitexplicit_free_surface];
simulations = [geostrophic_adjustment_simulation(free_surface, topology_type, arch = arch) for free_surface in free_surfaces, topology_type in topology_types, arch in archs];
data = [run_and_analyze(sim) for sim in simulations];
# run_and_analyze(simulations[3])

using GLMakie
using JLD2

file1 = jldopen("solution_PreconditionedConjugateGradient.jld2")
file2 = jldopen("solution_HeptadiagonalIterativeSolver.jld2")
file3 = jldopen("solution_SplitExplicitFreeSurface.jld2")

grid = file1["serialized/grid"]

x = grid.xᶜᵃᵃ[1:grid.Nx]
xf = grid.xᶠᵃᵃ[1:grid.Nx+1]
y = grid.yᵃᶜᵃ[1:grid.Ny]

iterations = parse.(Int, keys(file1["timeseries/t"]))
iterations = iterations[1:200]

iter = Node(0) # Node or Observable depending on Makie version
mid = Int(floor(grid.Ny / 2))
η0 = file1["timeseries/1/0"][:, mid, 1]
η1 = @lift(Array(file1["timeseries/1/"*string($iter)])[:, mid, 1])
η2 = @lift(Array(file2["timeseries/1/"*string($iter)])[:, mid, 1])
η3 = @lift(Array(file3["timeseries/1/"*string($iter)])[:, mid, 1])
u1 = @lift(Array(file1["timeseries/3/"*string($iter)])[:, mid, 1])
u2 = @lift(Array(file2["timeseries/3/"*string($iter)])[:, mid, 1])
u3 = @lift(Array(file3["timeseries/3/"*string($iter)])[:, mid, 1])
fig = Figure(resolution = (1400, 1000))
options = (; ylabelsize = 22,
    xlabelsize = 22, xgridstyle = :dash, ygridstyle = :dash, xtickalign = 1,
    xticksize = 10, ytickalign = 1, yticksize = 10, xlabel = "y [m]", xlims = extrema(x))
ax1 = Axis(fig[1, 1]; options..., ylabel = "η [m]", ylims = (-5e-4, 5e-4))
ax2 = Axis(fig[1, 2]; options..., ylabel = "u [m/s]", ylims = (-5e-5, 5e-5))

ηlines0 = lines!(ax1, x, η0, color = :black)
ηlines1 = lines!(ax1, x, η1, color = :red)
ηlines2 = lines!(ax1, x, η2, color = :blue)
ηlines3 = lines!(ax1, x, η3, color = :orange)
axislegend(ax1,
    [ηlines0, ηlines1, ηlines2, ηlines3],
    ["Initial Condition", "PCG", "Matrix", "Split-Explicit"])
ylims!(ax1, (-5e-4, 5e-3))
u0 = Array(file3["timeseries/3/"*string(0)])[:, mid, 1]
xf = length(u0) == length(xf) ? xf : x
ulines1 = lines!(ax2, xf, u1, color = :red)
ulines2 = lines!(ax2, xf, u2, color = :blue)
ulines3 = lines!(ax2, xf, u3, color = :orange)
axislegend(ax2,
    [ulines1, ulines2, ulines3],
    ["PCG", "Matrix", "Split-Explicit"])
ylims!(ax2, (-2e-4, 2e-4))
GLMakie.record(fig, "free_surface_bounded.mp4", iterations, framerate = 12) do i
    @info "Plotting iteration $i of $(iterations[end])..."
    iter[] = i
end


text:
a = 5

if a == 5
a = 10
end

println(a)


text:
module OutlierDetectionInterface
    using Reexport:@reexport
    using CategoricalArrays:CategoricalValue

    @reexport using MLJModelInterface
    const MMI = MLJModelInterface

    const CLASS_NORMAL = "normal"
    const CLASS_OUTLIER = "outlier"
    const DEFAULT_THRESHOLD = 0.9

    # those do not get reexported from MLJModelInterface
    const _process_model_def = MMI._process_model_def
    const _model_constructor = MMI._model_constructor
    const _model_cleaner = MMI._model_cleaner

    # detector types
    export Detector,
           SupervisedDetector,
           UnsupervisedDetector

    # interface types
    export DetectorModel,
           Scores,
           Labels,
           Data,
           Fit

    # constants
    export CLASS_NORMAL,
           CLASS_OUTLIER,
           DEFAULT_THRESHOLD

    # helpers
    export ncolons

    # macros
    export @detector,
           @default_frontend,
           @default_metadata

    # macro helpers
    export _process_model_def,
           _model_constructor,
           _model_cleaner

    include("base.jl")
    include("interface.jl")
    include("macros.jl")
end


text:
#!/usr/bin/julia

using Plots
using PlyIO

function plot_ply_2d(filename)
    ply = load_ply(filename)
    graph = Plots.scatter(ply["vertex"]["x"], ply["vertex"]["y"], title=filename)
    graph
end

function plot_ply_3d(filename)
    ply = load_ply(filename)
    graph = Plots.scatter(ply["vertex"]["x"], ply["vertex"]["y"], ply["vertex"]["z"], title=filename)
    graph
end

function convex_hull_2d()
    filename = "output.ply"
    plot_ply_2d(filename)
end

function raw_data()
    filename = "data/lettuce.ply"
    plot_ply_3d(filename)
end

function data_soil_removed()
    filename = "soil_removed.ply"
    plot_ply_3d(filename)
end

function data_3d()
    filename = "output.ply"
    plot_ply_3d(filename)
end

function data_soil_removed_o3d()
    filename = "soil_removed_o3d.ply"
    plot_ply_3d(filename)
end

# graph = raw_data()
graph = data_soil_removed()
# graph = data_soil_removed_o3d()
# graph = convex_hull_2d()
# graph = data_3d()
display(graph)


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# add call to list that is sorted by call priorities and then by arrival times within priorities
function addCallToQueueSortPriorityThenTime!(queuedCallList::Vector{Call}, call::Call)
	# find where to insert new call into list
	# maintain sorting by priority (lower value means a higher priority),
	# and then sorting by arrival time within priority
	# calls nearer to end of list are higher priority and arrived sooner
	i = length(queuedCallList) + 1
	while i > 1 && (call.priority > queuedCallList[i-1].priority ||
		(call.priority == queuedCallList[i-1].priority && call.arrivalTime > queuedCallList[i-1].arrivalTime))
		i -= 1
	end
	insert!(queuedCallList, i, call)

	if checkMode
		# check that calls are ordered first by priority, then by time
		for i = length(queuedCallList): -1 : 2
			@assert(queuedCallList[i].priority < queuedCallList[i-1].priority ||
				(queuedCallList[i].priority == queuedCallList[i-1].priority &&
				queuedCallList[i].arrivalTime <= queuedCallList[i-1].arrivalTime))
		end
	end
end