text:
include("setup_2d.jl")

println("dwt by lifting (N=",N,"x",N,"), ", L, " levels")
f(x0, L, wl) = for i = 1:tn; dwt(x0, wl, L); end
f(x0, L, wl);
@time f(x0, L, wl);

#@profile for i=1:tn dwt!(x0, L, wl, true); end
#let
#    @profile dwt(x0, wl, L);
#    Profile.print()
#end


text:
## Krylov.jl

struct KrylovJL{F,I,A,K} <: AbstractKrylovSubspaceMethod
    KrylovAlg::F
    gmres_restart::I
    window::I
    args::A
    kwargs::K
end

function KrylovJL(args...; KrylovAlg = Krylov.gmres!,
                  gmres_restart=0, window=0,
                  kwargs...)

    return KrylovJL(KrylovAlg, gmres_restart, window,
                    args, kwargs)
end

KrylovJL_CG(args...;kwargs...) =
    KrylovJL(args...; KrylovAlg=Krylov.cg!, kwargs...)
KrylovJL_GMRES(args...;kwargs...) =
    KrylovJL(args...; KrylovAlg=Krylov.gmres!, kwargs...)
KrylovJL_BICGSTAB(args...;kwargs...) =
    KrylovJL(args...; KrylovAlg=Krylov.bicgstab!, kwargs...)
KrylovJL_MINRES(args...;kwargs...) =
    KrylovJL(args...; KrylovAlg=Krylov.minres!, kwargs...)

function get_KrylovJL_solver(KrylovAlg)
    KS =
    if     (KrylovAlg === Krylov.lsmr!      ) Krylov.LsmrSolver
    elseif (KrylovAlg === Krylov.cgs!       ) Krylov.CgsSolver
    elseif (KrylovAlg === Krylov.usymlq!    ) Krylov.UsymlqSolver
    elseif (KrylovAlg === Krylov.lnlq!      ) Krylov.LnlqSolver
    elseif (KrylovAlg === Krylov.bicgstab!  ) Krylov.BicgstabSolver
    elseif (KrylovAlg === Krylov.crls!      ) Krylov.CrlsSolver
    elseif (KrylovAlg === Krylov.lsqr!      ) Krylov.LsqrSolver
    elseif (KrylovAlg === Krylov.minres!    ) Krylov.MinresSolver
    elseif (KrylovAlg === Krylov.cgne!      ) Krylov.CgneSolver
    elseif (KrylovAlg === Krylov.dqgmres!   ) Krylov.DqgmresSolver
    elseif (KrylovAlg === Krylov.symmlq!    ) Krylov.SymmlqSolver
    elseif (KrylovAlg === Krylov.trimr!     ) Krylov.TrimrSolver
    elseif (KrylovAlg === Krylov.usymqr!    ) Krylov.UsymqrSolver
    elseif (KrylovAlg === Krylov.bilqr!     ) Krylov.BilqrSolver
    elseif (KrylovAlg === Krylov.cr!        ) Krylov.CrSolver
    elseif (KrylovAlg === Krylov.craigmr!   ) Krylov.CraigmrSolver
    elseif (KrylovAlg === Krylov.tricg!     ) Krylov.TricgSolver
    elseif (KrylovAlg === Krylov.craig!     ) Krylov.CraigSolver
    elseif (KrylovAlg === Krylov.diom!      ) Krylov.DiomSolver
    elseif (KrylovAlg === Krylov.lslq!      ) Krylov.LslqSolver
    elseif (KrylovAlg === Krylov.trilqr!    ) Krylov.TrilqrSolver
    elseif (KrylovAlg === Krylov.crmr!      ) Krylov.CrmrSolver
    elseif (KrylovAlg === Krylov.cg!        ) Krylov.CgSolver
    elseif (KrylovAlg === Krylov.cg_lanczos!) Krylov.CgLanczosShiftSolver
    elseif (KrylovAlg === Krylov.cgls!      ) Krylov.CglsSolver
    elseif (KrylovAlg === Krylov.cg_lanczos!) Krylov.CgLanczosSolver
    elseif (KrylovAlg === Krylov.bilq!      ) Krylov.BilqSolver
    elseif (KrylovAlg === Krylov.minres_qlp!) Krylov.MinresQlpSolver
    elseif (KrylovAlg === Krylov.qmr!       ) Krylov.QmrSolver
    elseif (KrylovAlg === Krylov.gmres!     ) Krylov.GmresSolver
    elseif (KrylovAlg === Krylov.fom!       ) Krylov.FomSolver
    end

    return KS
end

function init_cacheval(alg::KrylovJL, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)

    KS = get_KrylovJL_solver(alg.KrylovAlg)

    memory = (alg.gmres_restart == 0) ? min(20, size(A,1)) : alg.gmres_restart

    solver = if(
        alg.KrylovAlg === Krylov.dqgmres! ||
        alg.KrylovAlg === Krylov.diom!    ||
        alg.KrylovAlg === Krylov.gmres!   ||
        alg.KrylovAlg === Krylov.fom!
       )
        KS(A, b, memory)
    elseif(
           alg.KrylovAlg === Krylov.minres! ||
           alg.KrylovAlg === Krylov.symmlq! ||
           alg.KrylovAlg === Krylov.lslq!   ||
           alg.KrylovAlg === Krylov.lsqr!   ||
           alg.KrylovAlg === Krylov.lsmr!
          )
        (alg.window != 0) ? KS(A,b; window=alg.window) : KS(A, b)
    else
        KS(A, b)
    end

    solver.x = u

    return solver
end

function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
    if cache.isfresh
        solver = init_cacheval(alg, cache.A, cache.b, cache.u, cache.Pl, cache.Pr, cache.maxiters, cache.abstol, cache.reltol, cache.verbose)
        cache = set_cacheval(cache, solver)
    end

    M = cache.Pl
    N = cache.Pr

    M = (M === Identity()) ? I : InvPreconditioner(M)
    N = (N === Identity()) ? I : InvPreconditioner(N)

    atol    = float(cache.abstol)
    rtol    = float(cache.reltol)
    itmax   = cache.maxiters
    verbose = cache.verbose ? 1 : 0

    args   = (cache.cacheval, cache.A, cache.b)
    kwargs = (atol=atol, rtol=rtol, itmax=itmax, verbose=verbose,
              alg.kwargs...)

    if cache.cacheval isa Krylov.CgSolver
        N !== I  &&
            @warn "$(alg.KrylovAlg) doesn't support right preconditioning."
        Krylov.solve!(args...; M=M,
                      kwargs...)
    elseif cache.cacheval isa Krylov.GmresSolver
        Krylov.solve!(args...; M=M, N=N,
                      kwargs...)
    elseif cache.cacheval isa Krylov.BicgstabSolver
        Krylov.solve!(args...; M=M, N=N,
                      kwargs...)
    elseif cache.cacheval isa Krylov.MinresSolver
        N !== I  &&
            @warn "$(alg.KrylovAlg) doesn't support right preconditioning."
        Krylov.solve!(args...; M=M,
                      kwargs...)
    else
        Krylov.solve!(args...; kwargs...)
    end

    return SciMLBase.build_linear_solution(alg,cache.u,nothing,cache)
end

## IterativeSolvers.jl

struct IterativeSolversJL{F,I,A,K} <: AbstractKrylovSubspaceMethod
    generate_iterator::F
    gmres_restart::I
    args::A
    kwargs::K
end

function IterativeSolversJL(args...;
                            generate_iterator = IterativeSolvers.gmres_iterable!,
                            gmres_restart=0, kwargs...)

    return IterativeSolversJL(generate_iterator, gmres_restart,
                              args, kwargs)
end

IterativeSolversJL_CG(args...; kwargs...) =
    IterativeSolversJL(args...;
                       generate_iterator=IterativeSolvers.cg_iterator!,
                       kwargs...)
IterativeSolversJL_GMRES(args...;kwargs...) =
    IterativeSolversJL(args...;
                       generate_iterator=IterativeSolvers.gmres_iterable!,
                       kwargs...)
IterativeSolversJL_BICGSTAB(args...;kwargs...) =
    IterativeSolversJL(args...;
                       generate_iterator=IterativeSolvers.bicgstabl_iterator!,
                       kwargs...)
IterativeSolversJL_MINRES(args...;kwargs...) =
    IterativeSolversJL(args...;
                       generate_iterator=IterativeSolvers.minres_iterable!,
                       kwargs...)

function init_cacheval(alg::IterativeSolversJL, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    restart = (alg.gmres_restart == 0) ? min(20, size(A,1)) : alg.gmres_restart

    kwargs = (abstol=abstol, reltol=reltol, maxiter=maxiters,
              alg.kwargs...)

    iterable = if alg.generate_iterator === IterativeSolvers.cg_iterator!
        Pr !== Identity() &&
          @warn "$(alg.generate_iterator) doesn't support right preconditioning"
        alg.generate_iterator(u, A, b, Pl;
                              kwargs...)
    elseif alg.generate_iterator === IterativeSolvers.gmres_iterable!
        alg.generate_iterator(u, A, b; Pl=Pl, Pr=Pr, restart=restart,
                              kwargs...)
    elseif alg.generate_iterator === IterativeSolvers.bicgstabl_iterator!
        Pr !== Identity() &&
          @warn "$(alg.generate_iterator) doesn't support right preconditioning"
        alg.generate_iterator(u, A, b, alg.args...; Pl=Pl,
                              abstol=abstol, reltol=reltol,
                              max_mv_products=maxiters*2,
                              alg.kwargs...)
    else # minres, qmr
        alg.generate_iterator(u, A, b, alg.args...;
                              abstol=abstol, reltol=reltol, maxiter=maxiters,
                              alg.kwargs...)
    end
    return iterable
end

function SciMLBase.solve(cache::LinearCache, alg::IterativeSolversJL; kwargs...)
    if cache.isfresh || !(typeof(alg) <: IterativeSolvers.GMRESIterable)
        solver = init_cacheval(alg, cache.A, cache.b, cache.u, cache.Pl, cache.Pr, cache.maxiters, cache.abstol, cache.reltol, cache.verbose)
        cache = set_cacheval(cache, solver)
    end
    purge_history!(cache.cacheval, cache.u, cache.b)

    cache.verbose && println("Using IterativeSolvers.$(alg.generate_iterator)")
    i = 0
    for iter in enumerate(cache.cacheval)
        i += 1
        cache.verbose && println("Iter: $(iter[1]), residual: $(iter[2])")
        # TODO inject callbacks KSP into solve cb!(cache.cacheval)
    end
    cache.verbose && println()

    return SciMLBase.build_linear_solution(alg,cache.u,nothing,cache; iters = i)
end

purge_history!(iter, x, b) = nothing
function purge_history!(iter::IterativeSolvers.GMRESIterable, x, b)
  iter.k = 1
  iter.x  = x
  fill!(x,false)
  iter.b  = b

  iter.residual.current = IterativeSolvers.init!(iter.arnoldi, iter.x, iter.b, iter.Pl, iter.Ax, initially_zero = true)
  IterativeSolvers.init_residual!(iter.residual, iter.residual.current)
  iter.β = iter.residual.current
  nothing
end

ext: jl


text:
module RelationalData
export Relation, TABLE_DUM, TABLE_DEE, restrict, extend, rename, project, naturaljoin, ⨝

"""
    Heading

A `Heading` defines names and types for a relation.
"""
struct Heading
  names
  t

  Heading() = new((), Tuple{})
  Heading(empty::Tuple{}) = new((), Tuple{})

  function Heading(h::NamedTuple{names, T} where {names, T<:Tuple{Vararg{DataType}}})
    sorted = (sort([keys(h)...])...,)
    head = NamedTuple{sorted}(h)
    new(sorted, Tuple{values(head)...})
  end

  function Heading(h::NamedTuple)
    sorted = (sort([keys(h)...])...,)
    head = NamedTuple{sorted}(h)
    new(sorted, Tuple{typeof(head).types...})
  end

  Heading(h::Heading, types::DataType...) = new(h.names, Tuple{types...})
end

"""
    shapeto(nt::NamedTuple, h::Heading)

Returns the named tuple converted to conform to the heading with attributes sorted in the canonical order
and types converted as appropriate (only safe conversions). Throws if conversion is impossible.
"""
function shapeto(nt::NamedTuple, h::Heading)
  length(nt) == length(h.names) || throw(DomainError(nt, "Does not match $h"))
  # For some reason we can't get the order right for the types here
  ont = NamedTuple{h.names}(nt)
  convert(NamedTuple{h.names, h.t}, ont)
end

function shapeto(empty::Tuple{}, h::Heading)
  length(h.names) == 0 || throw(DomainError(empty, "Does not match $h"))
  NamedTuple{h.names}(empty)
end

"""
    Relation

A Relation is a set of named tuples conforming to a Heading
"""
struct Relation{names, T}
  body::Set{NamedTuple}

  function Relation(ht::NamedTuple{names, T} where {names, T<:Tuple{Vararg{DataType}}})
    heading = Heading(ht)
    new{heading.names, heading.t}(Set{NamedTuple{heading.names, heading.t}}())
  end
  function Relation(ht::NamedTuple{names, T} where {names, T<:Tuple{Vararg{DataType}}},
     s::AbstractSet{V} where {V<:NamedTuple})
     heading = Heading(ht)
     new{heading.names, heading.t}(Set(shapeto.(s, Ref(heading))))
  end
  function Relation(s::AbstractSet{NamedTuple{names, T}}) where {names, T}
    heading = Heading((; zip(names, fieldtypes(T))...))
    new{heading.names, heading.t}(Set(shapeto.(s, Ref(heading))))
  end
  function Relation(nts::NamedTuple{names, T}...) where {names, T}
    heading = Heading((; zip(names, fieldtypes(T))...))
    new{heading.names, heading.t}(Set(shapeto.(nts, Ref(heading))))
  end
  function Relation(ht::NamedTuple{names, T},
   values::V...) where {V<:Tuple} where {names, T<:Tuple{Vararg{DataType}}}
   namedValues = NamedTuple{names}.(values)
   heading = Heading(ht)
   # Seems we have to reorder first, then apply types
   reorderedValues = NamedTuple{heading.names}.(namedValues)
   new{heading.names, heading.t}(Set(NamedTuple{heading.names, heading.t}.(reorderedValues)))
  end
  function Relation(ht::NamedTuple{names, T}, itr) where {names, T<:Tuple{Vararg{DataType}}}
    namedValues = NamedTuple{names}.(itr)
    heading = Heading(ht)
    # Seems we have to reorder first, then apply types
    reorderedValues = NamedTuple{heading.names}.(namedValues)
    new{heading.names, heading.t}(Set(NamedTuple{heading.names, heading.t}.(reorderedValues)))
  end
  Relation() = new{(), Tuple{}}(Set{NamedTuple{(), Tuple{}}}())
  Relation(empty::Tuple...) = new{(), Tuple{}}(Set([NamedTuple{(), Tuple{}}(empty[1])]))
  # trusted constructor
  Relation(names::Tuple{Vararg{Symbol}}, t::Type{T}, body::AbstractSet{V}) where {T<:Tuple, V<:NamedTuple} = new{names, t}(body)
end

function _promote(h::Heading, nt::NamedTuple)
  length(nt) == length(h.names) || throw(DomainError(nt, "Does not match $h"))
  # For some reason we can't get the order right for the types here
  ont = NamedTuple{h.names}(nt)
  htypes = promote_type.(fieldtypes(h.t), fieldtypes(typeof(ont).parameters[2]))
  Heading(h, htypes...)
end

function Relation(itr)
  next = iterate(itr)
  if next === nothing
    return Relation()
  end
  heading = Heading(next[1])
  next = iterate(itr, next[2])
  while next !== nothing
    (nt, state) = next
    heading = _promote(heading, nt)
    next = iterate(itr, state)
  end
  Relation(heading.names, heading.t, Set(NamedTuple{heading.names, heading.t}.(NamedTuple{heading.names}.(itr))))
end

const TABLE_DUM = Relation()
const TABLE_DEE = Relation(())

Base.iterate(r::Relation, i...) = iterate(r.body, i...)
Base.IteratorSize(r::Relation) = Base.HasLength()
Base.length(r::Relation) = length(r.body)
Base.IteratorEltype(r::Relation) = Base.HasEltype()
Base.eltype(r::Relation{names, T}) where {names, T} = NamedTuple{names, T}

Base.:(==)(q::Relation{names, T}, r::Relation{names, T}) where {names, T} = q.body == r.body
Base.:(==)(q::Relation{qnames, Q}, r::Relation{rnames, R}) where {qnames, Q, rnames, R} = false

Base.filter(f, r::Relation) = Relation(filter(f, r.body))
restrict(r::Relation, f) = filter(f, r)

"""
    extend(r::Relation, s::Symbol, f)

Extends every tuple of a relation with the given symbol assigned the value obtained by applying function f to the tuple.
"""
extend(r::Relation, s::Symbol, f) = Relation([merge(t, (; s => f(t))) for t in r])

function rename(r::Relation{names, T}, p::Pair{Symbol,Symbol}...) where {names, T}
  replacements = Dict(p...)
  renamed = map(n -> get(replacements, n, n), names)
  Relation(Set([NamedTuple{renamed, T}(values(nt)) for nt in r]))
end

project(r::Relation, names::Symbol...) = Relation(Set(NamedTuple{(names...,)}.(r)))

function ⨝(r1::Relation{names1, T1}, r2::Relation{names2, T2}) where {names1, T1, names2, T2}
  common = (intersect(names1, names2)...,)
  Relation(Set([merge(nt1, nt2) for nt1 in r1 for nt2 in r2 if NamedTuple{common}(nt1) == NamedTuple{common}(nt2)]))
end

function naturaljoin(r1::Relation, relations::Relation...)
  length(relations) == 1 && return r1 ⨝ relations[1]
  r1 ⨝ naturaljoin(relations...)
end

function Base.union(r1::Relation{names, T}, relations::Relation{names, T}...) where {names, T}
  Relation(names, T, union(r1.body, [r.body for r in relations]...))
end

Base.map(f, r::Relation) = Relation([f(t) for t in r])

end # module

max_stars_repo_name: tobega/RelationalData

max_issues_repo_name: tobega/RelationalData

max_forks_repo_name: tobega/RelationalData


text:
using Test, F1Method

using LinearAlgebra, DiffEqBase, ForwardDiff

# Set up:
# - overload `SteadyStateProblem` constructor
# - overload `solve` function
# - define solver algorithm (basic Newton here)
# - define type for that algorithm (here `MyAlg`)
include("simple_setup.jl")

@testset "quasi-Rosenbrock derivative" begin
    include("rosenbrock.jl")
end


text:
const CerebellarVemusLobuleI = AnatomicalStructure("cerebellar vermus lobule I")
const CerebellarVemusLobuleII = AnatomicalStructure("cerebellar vermus lobule II")
const CerebellarVemusLobuleIII = AnatomicalStructure("cerebellar vermus lobule III")
const CerebellarVemusLobuleIV = AnatomicalStructure("cerebellar vermus lobule IV")
const CerebellarVemusLobuleV = AnatomicalStructure("cerebellar vermus lobule V")
const CerebellarVemusLobuleVI = AnatomicalStructure("cerebellar vermus lobule VI")
const CerebellarVemusLobuleVII = AnatomicalStructure("cerebellar vermus lobule VII")
const CerebellarVemusLobuleVIII = AnatomicalStructure("cerebellar vermus lobule VIII")
const CerebellarVemusLobuleIX = AnatomicalStructure("cerebellar vermus lobule IX")
const CerebellarVemusLobuleX = AnatomicalStructure("cerebellar vermus lobule X")

const FlocculonodularLobe = AnatomicalStructure("flocculonodular lobe")

# TODO sub structures need more about cerebellar lob components
# children = Flocculonodular lobe
const Vestibulocerebellum = BilateralStructure("vestibulocerebellum")
const Archicerebellum = Vestibulocerebellum

# TODO sub structures for Cerebrocerebellum
const Cerebrocerebellum = BilateralStructure("cerebrocerebellum" )

"""
    Vermis

* Catani and Thiebaut de Schotten, 2008
* Dell'Acqua et al., 2013
* Meola et al., 2016a; Meola et al., 2016c; Wakana et al., 2004
"""
const Vermis = AnatomicalStructure("vermis")

const CerebellumCortex = BilateralStructure("cerebellum cortex")
#const CerebellumWhiteMatter = AnatomicalStructure("cerebellum white matter")

#= TODO this needs a different separation than white matter vs cortex
const Cerebellum = AnatomicalStructure(
    name = "Cerebellum",
    children = (
        CerebellumCortex,
        CerebellumWhiteMatter
    )
)
=#

# TODO double check paravermis
const Paravermis = BilateralStructure("paravermis")

# children Vermis and paravermis
const Spinocerebellum = BilateralStructure("spinocerebellum")
const Paleocerebellum = Spinocerebellum


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ResourceSkusResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{ResourceSku} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function ResourceSkusResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(ResourceSkusResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(ResourceSkusResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type ResourceSkusResult

const _property_map_ResourceSkusResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_ResourceSkusResult = Dict{Symbol,String}(Symbol("value")=>"Vector{ResourceSku}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ ResourceSkusResult }) = collect(keys(_property_map_ResourceSkusResult))
Swagger.property_type(::Type{ ResourceSkusResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ResourceSkusResult[name]))}
Swagger.field_name(::Type{ ResourceSkusResult }, property_name::Symbol) =  _property_map_ResourceSkusResult[property_name]

function check_required(o::ResourceSkusResult)
    (getproperty(o, Symbol("value")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ResourceSkusResult }, name::Symbol, val)
end


text:
module SodiumSeal

@static if VERSION >= v"1.3"
    using libsodium_jll: libsodium
else
    let deps = joinpath(@__DIR__, "..", "deps", "deps.jl")
        isfile(deps) || error(
            "SodiumSeal is not installed properly. ",
            """Run Pkg.build("SodiumSeal"), restart Julia and try again.""",
        )
        include(deps)
    end
end

export
    KeyPair,
    seal,
    unseal

using Base64: base64encode, base64decode

const SEALBYTES = Ref{Csize_t}()
const PUBLICKEYBYTES = Ref{Csize_t}()
const SECRETKEYBYTES = Ref{Csize_t}()

### Errors ###

struct SodiumError
    fun::Symbol
    code::Int
end

Base.showerror(io::IO, e::SodiumError) = print(io, "$(e.fun) returned error code $(e.code)")

macro check(allowed, ex)
    fun = ex.args[2].args[1]
    return quote
        code = $(esc(ex))
        code in $allowed || throw(SodiumError($fun, code))
        code
    end
end

### Key pairs ###

"""
    KeyPair([public[, secret]])

Construct a new `KeyPair` with existing keys, or generate a new one.
If the keys are `AbstractString`s, they are assumed to be Base64-encoded.
If you are only interested in encrypting, you need not supply the secret key.
"""
struct KeyPair
    public::Vector{Cuchar}
    secret::Vector{Cuchar}
end

KeyPair(public) = KeyPair(public, [])
KeyPair(public::AbstractString) = KeyPair(public, "")
KeyPair(public::AbstractString, secret::AbstractString) =
    KeyPair(base64decode(public), base64decode(secret))
function KeyPair()
    public = Vector{Cuchar}(undef, PUBLICKEYBYTES[])
    secret = Vector{Cuchar}(undef, SECRETKEYBYTES[])

    @check 0 ccall(
        (:crypto_box_keypair, libsodium),
        Cint,
        (Ptr{Cuchar}, Ptr{Cuchar}),
        public, secret,
    )

    return KeyPair(public, secret)
end

Base.show(io::IO, k::KeyPair) = print(io, "KeyPair(...)")

### Encryption ###

"""
    seal(plaintext, keypair::KeyPair) -> Union{String, Vector{UInt8}}

Encrypt some data.

If `plaintext` is an `AbstractString`, it is assumed to be Base64-encoded,
and the output is also a Base64-encoded `String`.
Otherwise, it is a `Vector{UInt8}`.
"""
seal(plaintext::AbstractString, keypair::KeyPair) =
    base64encode(seal(base64decode(plaintext), keypair))
function seal(plaintext, keypair::KeyPair)
    len = length(plaintext)
    dest = Vector{Cuchar}(undef, SEALBYTES[] + len)

    @check 0 ccall(
        (:crypto_box_seal, libsodium),
        Cint,
        (Ptr{Cuchar}, Ptr{Cuchar}, Culonglong, Ptr{Cuchar}),
        dest, plaintext, len, keypair.public,
    )

    return dest
end

### Decryption ###

"""
    unseal(ciphertext, keypair::KeyPair) -> Union{String, Vector{UInt8}}

Decrypt some data.

If `ciphertext` is an `AbstractString`, it is assumed to be Base64-encoded,
and the output is also a Base64-encoded `String`.
Otherwise, it is a `Vector{UInt8}`.
"""
unseal(ciphertext::AbstractString, keypair::KeyPair) =
    base64encode(unseal(base64decode(ciphertext), keypair))
function unseal(ciphertext, keypair)
    len = length(ciphertext)
    dest = Vector{Cuchar}(undef, len - SEALBYTES[])

    @check 0 ccall(
        (:crypto_box_seal_open, libsodium),
        Cint,
        (Ptr{Cuchar}, Ptr{Cuchar}, Culonglong, Ptr{Cuchar}, Ptr{Cuchar}),
        dest, ciphertext, len, keypair.public, keypair.secret,
    )

    return dest
end

### Init ###

function __init__()
    code = @check (0, 1) ccall((:sodium_init, libsodium), Cint, ())
    if code == 0
        SEALBYTES[] = ccall((:crypto_box_sealbytes, libsodium), Csize_t, ())
        PUBLICKEYBYTES[] = ccall((:crypto_box_publickeybytes, libsodium), Csize_t, ())
        SECRETKEYBYTES[] = ccall((:crypto_box_secretkeybytes, libsodium), Csize_t, ())
    end
end

end


text:
# function SCIPcreate()
#     s = SCIP_t(Array(Ptr{SCIP}, 1))
#     SCIPcreate(s)
#     finalizer(s, s->SCIPfree(s))
#     return s
# end

# Convenience constructors/destructors
function SCIPcreate()
    scip = SCIP_t(Array(Ptr{SCIP}, 1))
    SCIPcreate(scip)
    finalizer(scip, scip->SCIPfree(scip))
    return scip
end

function SCIPcreateVarBasic(scip::SCIP_t, name::String, lb::SCIP_Real, ub::SCIP_Real, obj::SCIP_Real, vartype::SCIP_VARTYPE)
    var = SCIP_VAR_t(Array(Ptr{SCIP_VAR}, 1))
    SCIPcreateVarBasic(scip, var, name, lb, ub, obj, vartype)
    #finalizer(var, var->SCIPreleaseVar(scip, var))
    return var
end

function SCIPcreateConsBasicLinear(scip::SCIP_t, name::String, vars::Vector{SCIP_VAR_t}, vals::Vector{SCIP_Real}, lhs::SCIP_Real, rhs::SCIP_Real)
    nvars = length(vars)
    @assert length(vals) == nvars
    cons = SCIP_CONS_t(Array(Ptr{SCIP_CONS}, 1))
    v = SCIP_VAR_t([pointer(var) for var in vars])
    SCIPcreateConsBasicLinear(scip, cons, name, nvars, v, vals, lhs, rhs)
    #finalizer(cons, cons->SCIPreleaseCons(scip, cons))
    return cons
end


text:
using InfiniteOpt, Ipopt, Plots
using KNITRO

# Parmeters
# Time span
t0 = 0
tf = 1
# Initial values
x1 = 1

# Model
model = InfiniteModel(KNITRO.Optimizer)

# infinite_parameter
@infinite_parameter(model, t ∈ [t0, tf], num_supports = 101)

# variable
@variable(model, x ≥ 0, Infinite(t))

@variable(model, u, Infinite(t), start = 0)

# objective
@objective(model, Min, x(1)^2 + ∫(u^2, t))

# constraint
@constraint(model, x(0) == x1)
@constraint(model, x_constr, ∂(x, t) == x + u)

# Optimization

print(model)

optimize!(model)

# Visulization
ts = value(t)
u_value = value(u)
x_value = value(x)

display(plot(ts, u_value))
display(plot(ts, x_value))

using DiffEqFlux, DifferentialEquations
const solveeq = DifferentialEquations.solve

tspan = (t0, tf)
ann = FastChain(FastDense(1, 64, tanh), FastDense(64, 1))
θ = initial_params(ann)
function sir_nn(du, u, p, t)
    du[1] = u[1] + ann([t], p)[1]
end
u0 = [x1]
ts = Float32.(collect(0.0:0.01:tspan[2]))
prob = ODEProblem(sir_nn, u0, tspan, θ)
sol_init = solveeq(prob, Vern9(), abstol = 1e-10, reltol = 1e-10)
plot(sol_init)
function predict_adjoint(θ)
    Array(solveeq(prob, Vern9(), p = θ, saveat = ts, sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP(true))))
end
function loss_adjoint(θ)
    x = predict_adjoint(θ)
    x[1, end]^2 + 0.01 * sum(abs2, [first(ann([t], θ)) for t in ts])
end
l = loss_adjoint(θ)

const losses = []
callback(θ, l, pred) = begin
    push!(losses, l)
    if length(losses) % 50 == 0
        println(losses[end])
        p = plot(solveeq(remake(prob, p = θ), Tsit5(), saveat = ts), lw = 3)
        #plot!(p, ts, [first(ann([t], θ)) for t in ts], label = "u(t)", lw = 3)
        display(p)
    end
    false
end
# cb = function (θ, l)
#     println(l)
#     p = plot(solveeq(remake(prob, p = θ), Tsit5(), saveat = ts), lw = 3)
#     #plot!(p, ts, [first(ann([t], θ)) for t in ts], label = "u(t)", lw = 3)
#     display(p)
#     return false
# end
# Display the ODE with the current parameter values.
cb(θ, l)
loss1 = loss_adjoint(θ)
res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.05), cb = cb, maxiters = 100)
res2 = DiffEqFlux.sciml_train(loss_adjoint, res1.u,
    BFGS(initial_stepnorm = 0.01), cb = cb, maxiters = 100,
    allow_f_increases = false)

using DataFrames
using CSV
p_value = res2.u
results_saving_opt = DataFrame()
results_saving_opt.t = ts
results_saving_opt.I = x_value
results_saving_opt.u = u_value
results_saving_opt.I_nn = Array(solveeq(remake(prob, p = p_value), Tsit5(), saveat = ts))[1, :]
results_saving_opt.u_nn = [first(ann([t], p_value)) for t in ts]
CSV.write("results_saving_Example3_2.csv", results_saving_opt)

plot(results_saving_opt.t, results_saving_opt.I, lw = 2, label = "I(t) by direct collection")
plot!(results_saving_opt.t, results_saving_opt.I_nn, lw = 2, label = "I(t) by deep learning")
xlabel!("t(day)")
ylabel!("I(t)")
savefig("Example3_2_I.png")

plot(results_saving_opt.t[2:end-1], results_saving_opt.u[2:end-1], lw = 2, label = "u(t) by direct collection")
plot!(results_saving_opt.t[2:end-1], results_saving_opt.u_nn[2:end-1], lw = 2, label = "u(t) by deep learning")
xlabel!("t(day)")
ylabel!("u(t)")
savefig("Example3_2_control.png")

#
opt_loss = objective_value(model)
# opt_loss = 1.75
opt_loss_nn = loss_adjoint(p_value)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-DEV
#     language: julia
#     name: julia-1.7
# ---

# %% [markdown]
# # 正弦函数の13次多項式による近似
#
# 区間 $[-\pi/4, \pi/4]$ で $\sin x$ を効率よく近似する13次多項式 $F_{13}(x)$ を計算する. 以下で $F_{10}(x)$ は10次多項式.
#
# * `Remez_abs`: $\sin(x) - F_{13}(x)$ のsupノルムを最小化
# * `Remez_rel`: $\dfrac{\sin(x) - F_{13}(x)}{\sin(x)}$ のsupノルムを最小化)
# * `Remez_rel1`: $\dfrac{\sin(x) - x - x^3 F_{10}(x)}{\sin(x)}$ のsupノルムを最小化)
# * `Remez_rel2`: $\dfrac{\sin(x) - x - x^3 F_{10}(x)}{x}$ のsupノルムを最小化)
#
# 最後の `Remez_rel2` の場合の $x + x^3 F_{10}(x)$ が [Julia の trig.jl](https://github.com/JuliaLang/julia/blob/master/base/special/trig.jl#L54) で使われている13次多項式に非常に近い.
#
# さらに, 低次の係数からFloat64での値に順次固定して求めて行くと, さらに近くなる.

# %% [markdown] toc=true
#
Table of Contents

#

1  Remez.jlの使用例

2  低次の係数からFloat64での値に順次固定して求めてみる

2.1  weight 1/x (minimaximize the pseudo relative error)

2.2  weight 1/sin(x) (minimaximize the exact relative error)

3  以下は汚い試行錯誤の残骸

# %% [markdown]
# https://github.com/simonbyrne/Remez.jl

# %%
using Remez
using Quadmath
using Polynomials
using Plots

relerr(x, y) = iszero(x - y) ? x - y : (x - y)/y

# %%
?ratfn_minimax

# %% [markdown]
# ## Remez.jlの使用例

# %%
@time N0, D0, E0, X0 = ratfn_minimax(sin, (-big(π)/4, big(π)/4), 13, 0)
sleep(0.1)
c0 = @. Float128(round(N0; digits=50))
c0[begin+1:2:end] |> display
sin_Remez_abs = Polynomial(c0)

# %%
function chebyshev_approx(d, F::Polynomial{CT}) where CT
    P = F
    for k in degree(F):-1:d+1
        degree(P) < k && continue
        a = coeffs(P)[end]
        T = ChebyshevT([zeros(CT, k); a/CT(2)^(k-1)])
        P = P - T
    end
    P
end

function maclaurin_sinpio4(d, CT=BigFloat)
    s = [isodd(k) ? (-1)^((k-1)÷2)*(big(π)/4)^k/factorial(CT(k)) : zero(CT) for k in 0:d]
    Polynomial(s)
end

function maclaurin_sin(d, CT=BigFloat)
    s = [isodd(k) ? (-1)^((k-1)÷2)/factorial(CT(k)) : zero(CT) for k in 0:d]
    Polynomial(s)
end

function ordinaryscale(F::Polynomial{CT}) where CT
    c = coeffs(F) ./ [(CT(π)/4)^k for k in 0:degree(F)]
    Polynomial(c)
end

S1 = chebyshev_approx(13, maclaurin_sinpio4(61)) |> ordinaryscale
S2 = chebyshev_approx(13, maclaurin_sinpio4(63)) |> ordinaryscale
S3 = chebyshev_approx(13, maclaurin_sinpio4(65)) |> ordinaryscale
@show [
    S1 - S2
    S2 - S3
]
sleep(0.1)

s = Float128.(coeffs(S2))
s[begin+1:2:end] |> display
sin_genkuroki_abs = Polynomial(s)

# %%
coeffs(sin_genkuroki_abs - sin_Remez_abs)[begin+1:2:end]

# %%
coeffs(sin_Remez_abs)[begin+1:2:end] .|> Float64

# %%
@time N1, D1, E1, X1 = ratfn_minimax(sin, (-big(π)/4, big(π)/4), 13, 0, (x,y)->1/(eps(BigFloat)+abs(y)))
sleep(0.1)
c1 = @. Float128(round(N1; digits=50))
c1[begin+1:2:end] |> display
sin_Remez_rel = Polynomial(c1)

# %%
coeffs(sin_Remez_rel)[begin+1:2:end] .|> Float64

# %%
sinmxox3(x) = iszero(x) ? -1/big(6) : (sin(x) - x)/x^3 # sin(x) minus X over x^3
w_sinmxox3(x, y) = iszero(x) ? x : abs(x^3/sin(x))

@time N2, D2, E2, X2 = ratfn_minimax(sinmxox3, (-big(π)/4, big(π)/4), 10, 0, w_sinmxox3)
sleep(0.1)
c2 = @. Float128(round(N2; digits=50))
c2 = Float128[0; 1; 0; c2]
c2[begin+1:2:end] |> display
sin_Remez_rel1 = Polynomial(c2)

# %%
coeffs(sin_Remez_rel1)[begin+1:2:end] .|> Float64

# %% [markdown]
# * https://github.com/JuliaLang/julia/blob/master/base/special/trig.jl#L54
# * https://github.com/JuliaMath/openlibm/blob/master/src/k_sin.c#L51
# * https://github.com/freebsd/freebsd-src/blob/main/lib/msun/src/k_sin.c#L50

# %%
DS0 =  1.0
DS1 = -1.66666666666666324348e-01
DS2 =  8.33333333332248946124e-03
DS3 = -1.98412698298579493134e-04
DS4 =  2.75573137070700676789e-06
DS5 = -2.50507602534068634195e-08
DS6 =  1.58969099521155010221e-10

sin_julia = Polynomial(Float128[0.0, DS0, 0.0, DS1, 0.0, DS2, 0.0, DS3, 0.0, DS4, 0.0, DS5, 0.0, DS6])

# %%
coeffs(sin_Remez_rel1 - sin_julia)[begin+1:2:end]

# %%
[0.0, DS0, 0.0, DS1, 0.0, DS2, 0.0, DS3, 0.0, DS4, 0.0, DS5, 0.0, DS6][begin+1:2:end]

# %%
Float64.(coeffs(sin_Remez_rel1)[begin+1:2:end])

# %%
(Float64[0.0, DS0, 0.0, DS1, 0.0, DS2, 0.0, DS3, 0.0, DS4, 0.0, DS5, 0.0, DS6] -
    Float64.(coeffs(sin_Remez_rel1)))[begin+1:2:end]

# %%
sinmxox3(x) = iszero(x) ? -1/big(6) : (sin(x) - x)/x^3 # sin(x) minus X over x^3

@time N3, D3, E3, X3 = ratfn_minimax(sinmxox3, (-big(π)/4, big(π)/4), 10, 0, (x,y)->x^2)
sleep(0.1)
c3 = @. Float128(round(N3; digits=50))
c3 = Float128[0; 1; 0; c3]
c3[begin+1:2:end] |> display
sin_Remez_rel2 = Polynomial(c3)

# %%
coeffs(sin_Remez_rel2 - sin_julia)[begin+1:2:end]

# %%
Float64.(coeffs(sin_Remez_rel1)[begin+1:2:end])

# %%
[0.0, DS0, 0.0, DS1, 0.0, DS2, 0.0, DS3, 0.0, DS4, 0.0, DS5, 0.0, DS6][begin+1:2:end]

# %%
([0.0, DS0, 0.0, DS1, 0.0, DS2, 0.0, DS3, 0.0, DS4, 0.0, DS5, 0.0, DS6] -
    Float64.(coeffs(sin_Remez_rel2)))[begin+1:2:end]

# %%
x = range(eps(Float128), π/4; length=1000)
y_Remez_abs  = sin_Remez_abs.(x)
y_Remez_rel  = sin_Remez_rel.(x)
y_Remez_rel1 = sin_Remez_rel1.(x)
y_Remez_rel2 = sin_Remez_rel2.(x)
y_julia      = sin_julia.(x)
y            = sin.(x)

error_Remez_abs  = y_Remez_abs  - y .|> Float64
error_Remez_rel  = y_Remez_rel  - y .|> Float64
error_Remez_rel1 = y_Remez_rel1 - y .|> Float64
error_Remez_rel2 = y_Remez_rel2 - y .|> Float64
error_julia      = y_julia      - y .|> Float64

relerror_Remez_abs  = relerr.(y_Remez_abs , y) .|> Float64
relerror_Remez_rel  = relerr.(y_Remez_rel , y) .|> Float64
relerror_Remez_rel1 = relerr.(y_Remez_rel1, y) .|> Float64
relerror_Remez_rel2 = relerr.(y_Remez_rel2, y) .|> Float64
relerror_julia      = relerr.(y_julia     , y) .|> Float64

@show extrema(error_Remez_abs)
@show extrema(error_Remez_rel)
@show extrema(error_Remez_rel1)
@show extrema(error_Remez_rel2)
@show extrema(error_julia)
println()
@show extrema(relerror_Remez_abs)
@show extrema(relerror_Remez_rel)
@show extrema(relerror_Remez_rel1)
@show extrema(relerror_Remez_rel2)
@show extrema(relerror_julia)
;

# %%
plot(; legend=:bottom)
plot!(; title="10¹⁶ × absolute errors", titlefontsize=11)
plot!(x, 10^16*error_julia; label="Julia", ls=:solid)
plot!(x, 10^16*error_Remez_rel;  label="Remez_rel", ls=:dash)
plot!(x, 10^16*error_Remez_rel1; label="Remez_rel1", ls=:dot, lw=1.8)
plot!(x, 10^16*error_Remez_rel1; label="Remez_rel2", ls=:dash)
plot!(x, 10^16*error_Remez_abs;  label="Remez_abs", ls=:dashdot)

# %%
plot(; legend=:bottom)
plot!(; title="10¹⁶ × relative errors", titlefontsize=11)
plot!(x, 10^16*relerror_julia; label="Julia", ls=:solid)
plot!(x, 10^16*relerror_Remez_rel;  label="Remez_rel", ls=:dash)
plot!(x, 10^16*relerror_Remez_rel1; label="Remez_rel1", ls=:dot, lw=1.8)
plot!(x, 10^16*relerror_Remez_rel1; label="Remez_rel2", ls=:dash)
plot!(x, 10^16*relerror_Remez_abs;  label="Remez_abs", ls=:dashdot)

# %%
plot(; legend=:bottom)
plot!(; title="10¹⁶ × relative errors", titlefontsize=11)
plot!(x, 10^16*relerror_julia; label="Julia", ls=:solid)
plot!(x, 10^16*relerror_Remez_rel;  label="Remez_rel",  ls=:dash)
plot!(x, 10^16*relerror_Remez_rel1; label="Remez_rel1", ls=:dot, lw=1.8)
plot!(x, 10^16*relerror_Remez_rel1; label="Remez_rel2", ls=:dash)
plot!(x, 10^16*relerror_Remez_abs;  label="Remez_abs",  ls=:dashdot)
plot!(; ylim=(-0.08, 0.06))

# %%
plot(; legend=:bottomleft)
plot!(; title="10¹⁶ × relative errors", titlefontsize=11)
plot!(x, 10^16*relerror_julia; label="Julia", ls=:solid)
plot!(x, 10^16*relerror_Remez_rel2; label="Remez_rel2", ls=:dash)

# %%
sin_Remez_rel2_F64 = convert(Polynomial{Float64}, sin_Remez_rel2)
sin_julia_F64 = Polynomial([0.0, DS0, 0.0, DS1, 0.0, DS2, 0.0, DS3, 0.0, DS4, 0.0, DS5, 0.0, DS6])

x = range(eps(Float128), π/4; length=1000)
y_Remez_rel2_F64 = sin_Remez_rel2_F64.(x)
y_julia_F64 = sin_julia_F64.(x)

relerror_Remez_rel2_F64 = relerr.(y_Remez_rel2_F64, y) .|> Float64
relerror_julia_F64 = relerr.(y_julia_F64, y) .|> Float64

plot(; legend=:bottomleft)
plot!(; title="10¹⁶ × relative errors", titlefontsize=11)
plot!(x, 10^16*relerror_julia_F64; label="Julia", ls=:solid)
plot!(x, 10^16*relerror_Remez_rel2_F64; label="Remez_rel2_F64", ls=:dash)

# %% [markdown]
# ## 低次の係数からFloat64での値に順次固定して求めてみる

# %% [markdown]
# ### weight 1/x (minimaximize the pseudo relative error)

# %%
a1 = big"1.0"

# %%
r(x) = √abs(x)
F(x) = iszero(x) ? -1/typeof(x)(6) : (sin(x) - x*a1)/x^3
f(x) = F(r(x))

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 5, 0, (x,y)->abs(x))
sleep(0.1)

c3 = Float64.(N)
@show a3 = c3[1]
c3

# %%
F(x) = iszero(x) ? 1/big(120) : (sin(x) - x*(a1 + x^2*a3))/x^5
f(x) = F(r(x))

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 4, 0, (x,y)->abs(x)^2)
sleep(0.1)

c5 = Float64.(N)
@show a5 = c5[1]
c5

# %%
F(x) = iszero(x) ? -1/factorial(big(7)) : (sin(x) - x*(a1 + x^2*(a3 + x^2*a5)))/x^7
f(x) = F(r(x))

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 3, 0, (x,y)->abs(x)^3)
sleep(0.1)

c7 = Float64.(N)
@show a7 = c7[1]
c7

# %%
F(x) = iszero(x) ? 1/factorial(big(9)) : (sin(x) - x*(a1 + x^2*(a3 + x^2*(a5 + x^2*a7))))/x^9
f(x) = F(r(x))

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 2, 0, (x,y)->abs(x)^4)
sleep(0.1)

c9 = Float64.(N)
@show a9 = c9[1]
c9

# %%
F(x) = iszero(x) ? 1/factorial(big(11)) : (sin(x) - x*(a1 + x^2*(a3 + x^2*(a5 + x^2*(a7 + x^2*a9)))))/x^11
f(x) = F(r(x))

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 1, 0, (x,y)->abs(x)^5)
sleep(0.1)

c11 = Float64.(N)
@show a11 = c11[1]
c11

# %%
F(x) = iszero(x) ? 1/factorial(big(13)) :
    (sin(x) - x*(a1 + x^2*(a3 + x^2*(a5 + x^2*(a7 + x^2*(a9 + x^2*a11))))))/x^13
f(x) = F(r(x))

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 0, 0, (x,y)->abs(x)^6)
sleep(0.1)

c13 = Float64.(N)
@show a13 = c13[1]
c13

# %%
a = Float64[0.0, a1, 0.0, a3, 0.0, a5, 0.0, a7, 0.0, a9, 0.0, a11, 0.0, a13]
a[begin+1:2:end]

# %%
a_julia = [0.0, DS0, 0.0, DS1, 0.0, DS2, 0.0, DS3, 0.0, DS4, 0.0, DS5, 0.0, DS6]
a_julia[begin+1:2:end]

# %%
(a - a_julia)[begin+1:2:end]

# %%
sin_a = Polynomial(a)
sin_j = Polynomial(a_julia)

x = range(eps(Float128), Float128(π)/4; length=1000)
y_a  = sin_a.(x)
y_j  = sin_j.(x)
y = sin.(x)

relerror_a = relerr.(y_a, y) .|> Float64
relerror_j = relerr.(y_j, y) .|> Float64
relerror_aj = relerr.(y_a, y_j) .|> Float64

@show extrema(relerror_a)
@show extrema(relerror_j)
@show extrema(relerror_aj)

plot(; legend=:bottomleft)
plot!(; title="10¹⁶ × relative errors", titlefontsize=11)
plot!(x, 10^16*relerror_j; label="Julia", ls=:solid, lw=1.2)
plot!(x, 10^16*relerror_a; label="Mine",  ls=:dash, lw=1.2)

# %%
sin_a = Polynomial(a)
sin_j = Polynomial(a_julia)

x = range(eps(), π/4; length=1000)
y_a  = sin_a.(x)
y_j  = sin_j.(x)
y = sin.(Float128.(x))

relerror_a = relerr.(y_a, y) .|> Float64
relerror_j = relerr.(y_j, y) .|> Float64
relerror_aj = relerr.(y_a, y_j) .|> Float64

@show extrema(relerror_a)
@show extrema(relerror_j)
@show extrema(relerror_aj)

P = plot(; legend=:bottomleft)
plot!(; title="10¹⁶ × relative errors (Julia)", titlefontsize=11)
plot!(x, 10^16*relerror_j; label="", lw=0.5)

Q = plot(; legend=:bottomleft)
plot!(; title="10¹⁶ × relative errors (Mine)", titlefontsize=11)
plot!(x, 10^16*relerror_a; label="", lw=0.5)

plot(P, Q; size=(1000, 300))

# %% [markdown]
# ### weight 1/sin(x) (minimaximize the exact relative error)

# %%
a1 = big"1.0"

# %%
r(x) = √abs(x)
F(x) = iszero(x) ? -1/typeof(x)(6) : (sin(x) - x*a1)/x^3
f(x) = F(r(x))
W(x, y) = iszero(x) ? x : x^3/sin(x)
w(x, y) = W(r(x), y)

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 5, 0, w)
sleep(0.1)

@show N
c3 = Float64.(N)
@show a3 = c3[1]
c3

# %%
F(x) = iszero(x) ? -1/typeof(x)(6) : (sin(x) - x*(a1 + x^2*a3))/x^5
f(x) = F(r(x))
W(x, y) = iszero(x) ? x : x^5/sin(x)
w(x, y) = W(r(x), y)

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 4, 0, w)
sleep(0.1)

c5 = Float64.(N)
@show a5 = c5[1]
c5

# %%
F(x) = iszero(x) ? -1/factorial(big(7)) : (sin(x) - x*(a1 + x^2*(a3 + x^2*a5)))/x^7
f(x) = F(r(x))
W(x, y) = iszero(x) ? x : x^7/sin(x)
w(x, y) = W(r(x), y)

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 3, 0, w)
sleep(0.1)

c7 = Float64.(N)
@show a7 = c7[1]
c7

# %%
F(x) = iszero(x) ? 1/factorial(big(9)) : (sin(x) - x*(a1 + x^2*(a3 + x^2*(a5 + x^2*a7))))/x^9
f(x) = F(r(x))
W(x, y) = iszero(x) ? x : x^9/sin(x)
w(x, y) = W(r(x), y)

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 2, 0, w)
sleep(0.1)

c9 = Float64.(N)
@show a9 = c9[1]
c9

# %%
F(x) = iszero(x) ? 1/factorial(big(9)) : (sin(x) - x*(a1 + x^2*(a3 + x^2*(a5 + x^2*(a7 + x^2*a9)))))/x^11
f(x) = F(r(x))
W(x, y) = iszero(x) ? x : x^11/sin(x)
w(x, y) = W(r(x), y)

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 1, 0, w)
sleep(0.1)

c11 = Float64.(N)
@show a11 = c11[1]
c11

# %%
F(x) = iszero(x) ? 1/factorial(big(9)) :
    (sin(x) - x*(a1 + x^2*(a3 + x^2*(a5 + x^2*(a7 + x^2*(a9 + x^2*a11))))))/x^13
f(x) = F(r(x))
W(x, y) = iszero(x) ? x : x^13/sin(x)
w(x, y) = W(r(x), y)

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 0, 0, w)
sleep(0.1)

c13 = Float64.(N)
@show a13 = c13[1]
c13

# %%
b = Float64[0.0, a1, 0.0, a3, 0.0, a5, 0.0, a7, 0.0, a9, 0.0, a11, 0.0, a13]
b[begin+1:2:end]

# %%
sin_b = Polynomial(b)

x = range(eps(Float128), Float128(π)/4; length=1000)
y_a  = sin_a.(x)
y_b  = sin_b.(x)
y_j  = sin_j.(x)
y = sin.(x)

relerror_a = relerr.(y_a, y) .|> Float64
relerror_b = relerr.(y_b, y) .|> Float64
relerror_j = relerr.(y_j, y) .|> Float64
relerror_aj = relerr.(y_a, y_j) .|> Float64
relerror_bj = relerr.(y_b, y_j) .|> Float64

@show extrema(relerror_a)
@show extrema(relerror_b)
@show extrema(relerror_j)
@show extrema(relerror_aj)
@show extrema(relerror_bj)

plot(; legend=:bottomleft)
plot!(; title="10¹⁶ × relative errors", titlefontsize=11)
plot!(x, 10^16*relerror_j; label="Julia", ls=:solid, lw=1.2)
plot!(x, 10^16*relerror_a; label="My a",  ls=:dash, lw=1.2)
plot!(x, 10^16*relerror_b; label="My b",  ls=:dashdot, lw=1.2)

# %%
sin_a = Polynomial(a)
sin_b = Polynomial(b)
sin_j = Polynomial(a_julia)

x = range(eps(), π/4; length=1000)
y_a  = sin_a.(x)
y_b  = sin_b.(x)
y_j  = sin_j.(x)
y = sin.(Float128.(x))

relerror_a = relerr.(y_a, y) .|> Float64
relerror_b = relerr.(y_b, y) .|> Float64
relerror_j = relerr.(y_j, y) .|> Float64
relerror_aj = relerr.(y_a, y_j) .|> Float64
relerror_bj = relerr.(y_b, y_j) .|> Float64

@show extrema(relerror_a)
@show extrema(relerror_b)
@show extrema(relerror_j)
@show extrema(relerror_aj)
@show extrema(relerror_bj)

P = plot(; legend=:bottomleft)
plot!(; title="10¹⁶ × relative errors (Julia)", titlefontsize=11)
plot!(x, 10^16*relerror_j; label="", lw=0.5)

Q = plot(; legend=:bottomleft)
plot!(; title="10¹⁶ × relative errors (My b)", titlefontsize=11)
plot!(x, 10^16*relerror_b; label="", lw=0.5)

plot(P, Q; size=(1000, 300))

# %% [markdown]
# ## 以下は汚い試行錯誤の残骸

# %%
r(x) = √abs(x)
f(x) = iszero(x) ? -1/big(6) : (sin(r(x)) - r(x))/r(x)^3
w_f(x, y) = iszero(x) ? x : abs(r(x)^3/sin(r(x)))

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 5, 0, w_f)
sleep(0.1)

c = eltype(N3)[0; 1; 0; [isodd(k) ? N[(k+1)÷2] : 0 for k in 1:11]]
c = Float128.(c)
c[begin+1:2:end] |> display
sin_Remez_rel1′ = Polynomial(c)

# %%
sin_Remez_rel1′ - sin_Remez_rel1

# %%
r(x) = √abs(x)
f(x) = iszero(x) ? -1/big(6) : (sin(r(x)) - r(x))/r(x)^3
w_f(x, y) = iszero(x) ? x : abs(r(x)^3/sin(r(x)))

@time N, D, E, X = ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 5, 0, (x,y)->abs(x))
sleep(0.1)

c = eltype(N)[0; 1; 0; [isodd(k) ? N[(k+1)÷2] : 0 for k in 1:11]]
c = Float128.(c)
c[begin+1:2:end] |> display
sin_Remez_rel2′ = Polynomial(c)

# %%
sin_Remez_rel2′ - sin_Remez_rel2

# %%
using Remez: ratfn_leastsquares, ratfn_eval, find_extrema, winnow_extrema, ratfn_equal_deviation

function my_ratfn_minimax(f, interval, n, d,
                       w = (x,y)->BigFloat(1);
        epsbits = precision(BigFloat),
        threshold = 2^(-epsbits/3)
    )
    # We start off by finding a least-squares approximation. This
    # doesn't need to be perfect, but if we can get it reasonably good
    # then it'll save iterations in the refining stage.
    #
    # Least-squares approximations tend to look nicer in a minimax
    # sense if you evaluate the function at a big pile of Chebyshev
    # nodes rather than uniformly spaced points. These values will
    # also make a good grid to use for the initial search for error
    # extrema, so we'll keep them around for that reason too.

    # Construct the grid.
    lo = BigFloat(minimum(interval))
    hi = BigFloat(maximum(interval))

    local grid
    let
        mid = (hi+lo)/2
        halfwid = (hi-lo)/2
        nnodes = 16 * (n+d+1)
        grid = [ mid - halfwid * cospi(big(i)/big(nnodes)) for i=0:nnodes ]
    end

    # Find the initial least-squares approximation.
    (nc, dc) = ratfn_leastsquares(f, grid, n, d, w)

    # Threshold of convergence. We stop when the relative difference
    # between the min and max (winnowed) error extrema is less than
    # this.
    #
    # This is set to the cube root of machine epsilon on a more or
    # less empirical basis, because the rational-function case will
    # not converge reliably if you set it to only the square root.
    # (Repeatable by using the --test mode.) On the assumption that
    # input and output error in each iteration can be expected to be
    # related by a simple power law (because it'll just be down to how
    # many leading terms of a Taylor series are zero), the cube root
    # was the next thing to try.
    #epsbits = precision(BigFloat)
    #threshold = 2^(-epsbits/3)

    # Main loop.
    maxiters = 10^3
    iter = 0
    max_err = zero(BigFloat)
    extrema = Vector{Tuple{BigFloat, BigFloat}}()
    while iter < maxiters
        iter += 1
        iszero(mod(iter, 100)) && print(iter, " ")
        # Find all the error extrema we can.
        function compute_error(x)
            real_y = f(x)
            approx_y = ratfn_eval(nc, dc, x)
            return (approx_y - real_y) * w(x, real_y)
        end
        extrema = find_extrema(compute_error, grid)

        # Winnow the extrema down to the right number, and ensure they
        # have alternating sign.
        extrema = winnow_extrema(extrema, n+d+2)

        # See if we've finished.
        min_err = minimum([abs(y) for (x,y) = extrema])
        max_err = maximum([abs(y) for (x,y) = extrema])
        variation = (max_err - min_err) / max_err
        if variation < threshold
            return nc, dc, max_err, extrema
        end

        # If not, refine our function by equalising the error at the
        # extrema points, and go round again.
        (nc, dc) = ratfn_equal_deviation(f, map(x->x[1], extrema),
                                         n, d, max_err, w)
    end

    println("maxiters = $maxiters")
    return nc, dc, max_err, extrema
end

# %%
r(x) = √abs(x)
f(x) = iszero(x) ? -1/big(6) : (sin(r(x)) - r(x))/r(x)^3
w_f(x, y) = iszero(x) ? x : abs(r(x)^3/sin(r(x)))

setprecision(80)

@time N, D, E, X = my_ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 5, 0, w_f;
    threshold = 2^(-precision(BigFloat)/7))
sleep(0.1)

c = eltype(N)[0; 1; 0; [isodd(k) ? N[(k+1)÷2] : 0 for k in 1:11]]
c[begin+1:2:end] |> display

sin_Remez = Polynomial(c)
sin_Remez |> display

setprecision(256);

# %%
r(x) = √abs(x)
f(x) = iszero(x) ? -1/big(6) : (sin(r(x)) - r(x))/r(x)^3
w_f(x, y) = iszero(x) ? x : abs(r(x)^3/sin(r(x)))

setprecision(64)

@time N, D, E, X = my_ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 5, 0, w_f;
    threshold = 2^(-precision(BigFloat)/7))
sleep(0.1)

c = eltype(N)[0; 1; 0; [isodd(k) ? N[(k+1)÷2] : 0 for k in 1:11]]
c[begin+1:2:end] |> display

sin_Remez = Polynomial(c)
sin_Remez |> display

setprecision(256);

# %%
r(x) = √abs(x)
f(x) = iszero(x) ? -1/big(6) : (sin(r(x)) - r(x))/r(x)^3
w_f(x, y) = iszero(x) ? x : abs(r(x)^3/sin(r(x)))

setprecision(60)

@time N, D, E, X = my_ratfn_minimax(f, (big"0.0", (big(π)/4)^2), 5, 0, w_f;
    threshold = 2^(-precision(BigFloat)/12))
sleep(0.1)

c = eltype(N)[0; 1; 0; [isodd(k) ? N[(k+1)÷2] : 0 for k in 1:11]]
c[begin+1:2:end] |> display
c[begin+1:2:end] .|> Float64 |> display

sin_Remez = Polynomial(c)
sin_Remez |> display

setprecision(256);

# %%
[0.0, DS0, 0.0, DS1, 0.0, DS2, 0.0, DS3, 0.0, DS4, 0.0, DS5, 0.0, DS6][begin+1:2:end] |> display

# %%
g1(x) = iszero(x) ? x : (sin(x) - x)/x
@time N, D, E, X = ratfn_minimax(g1, (-big(π)/4, big(π)/4), 12, 0)
sleep(0.1)
c = @. Float128(round(N; digits=50))
c = Float128[0; 1+c[1]; c[2:end]]
c[begin+1:2:end] |> display
sin_Remez_g1 = Polynomial(c)
sin_Remez_g1 |> display
y_g1 = sin_Remez_g1.(x)
relerror_Remez_g1 = relerr.(y_g1, y) .|> Float64

g2(x) = iszero(x) ? x : (sin(x) - x)/x^2
@time N, D, E, X = ratfn_minimax(g2, (-big(π)/4, big(π)/4), 11, 0, (x,y)->abs(x))
sleep(0.1)
c = @. Float128(round(N; digits=50))
c = Float128[0; 1; c]
c[begin+1:2:end] |> display
sin_Remez_g2 = Polynomial(c)
sin_Remez_g2 |> display
y_g2 = sin_Remez_g2.(x)
relerror_Remez_g2 = relerr.(y_g2, y) .|> Float64

g3(x) = iszero(x) ? -1/typeof(x)(6) : (sin(x) - x)/x^3
@time N, D, E, X = ratfn_minimax(g3, (-big(π)/4, big(π)/4), 10, 0, (x,y)->x^2)
sleep(0.1)
c = @. Float128(round(N; digits=50))
c = Float128[0; 1; 0; c]
c[begin+1:2:end] |> display
sin_Remez_g3 = Polynomial(c)
sin_Remez_g3 |> display
y_g3 = sin_Remez_g3.(x)
relerror_Remez_g3 = relerr.(y_g3, y) .|> Float64

plot(; legend=:bottomleft)
plot!(; title="10¹⁶ × relative errors (gₖ(x)=(sin(x)-x)/xᵏ)", titlefontsize=11)
plot!(x, 10^16*relerror_julia; label="Julia", ls=:solid)
plot!(x, 10^16*relerror_Remez_g1;  label="Remez_g1",  ls=:dash)
plot!(x, 10^16*relerror_Remez_g2;  label="Remez_g2",  ls=:dot, lw=1.8)
plot!(x, 10^16*relerror_Remez_g3;  label="Remez_g3",  ls=:dashdot)

# %%
cg2 = Float64.(coeffs(sin_Remez_g2))[begin+1:2:end]

# %%
cg3 = Float64.(coeffs(sin_Remez_g3))[begin+1:2:end]

# %%
cg2 == cg3

# %%
cjulia = [0.0, DS0, 0.0, DS1, 0.0, DS2, 0.0, DS3, 0.0, DS4, 0.0, DS5, 0.0, DS6][begin+1:2:end]

# %%
cg3 - cjulia

# %%

max_stars_repo_path: 0004/Remez.jl example - degree-13 polynomial approximations of sine.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using DataFrames

## ACTmath
ACTmath = DataFrame(
    [0 0 0
    1 1 1
    2 1 3
    3 3 13
    4 9 42
    5 18 59
    6 59 95
    7 67 131
    8 91 158
    9 144 161
    10 149 194
    11 192 164
    12 192 166
    13 192 197
    14 201 177
    15 204 158
    16 217 169
    17 181 132
    18 184 158
    19 170 151
    20 201 134
    21 147 137
    22 163 122
    23 147 110
    24 140 116
    25 147 132
    26 126 104
    27 113 104
    28 100 114
    29 106 97
    30 107 107
    31 91 88
    32 83 80
    33 73 79
    34 72 70
    35 75 61
    36 50 48
    37 37 47
    38 38 29
    39 23 32
    40 15 12],
    ["scale","xcount", "ycount"]
)

## KBneatX
KBneatX = DataFrame(
    [8 3
    21 6
    31 10
    7 2
    18 5
    36 12
    17 5
    8 2
    8 5
    20 6
    12 4
    11 5
    27 9
    7 2
    12 8
    7 2
    16 5
    27 9
    16 5
    7 3
    16 5
    22 7
    13 6
    9 4
    15 4
    29 7
    14 5
    26 7
    14 6
    20 5
    10 3
    17 5
    18 6
    19 4
    14 3
    12 5
    11 4
    20 7
    17 7
    20 5
    13 5
    10 5
    12 4
    15 5
    16 6
    8 3
    11 3
    32 10
    11 6
    19 5
    17 6
    25 8
    26 8
    10 2
    4 1
    20 8
    11 4
    20 4
    14 4
    15 5
    20 6
    25 9
    20 7
    8 4
    22 8
    10 2
    15 4
    25 7
    12 6
    8 2
    19 6
    14 5
    11 5
    30 10
    15 4
    25 7
    21 6
    7 0
    31 10
    10 4
    8 6
    8 4
    10 3
    16 6
    10 2
    17 4
    16 4
    17 4
    18 8
    14 6
    20 5
    15 3
    7 1
    17 6
    9 2
    12 5
    10 5
    12 3
    14 4
    20 6
    14 4
    20 5
    21 8
    9 4
    19 7
    12 3
    13 4
    14 5
    18 6
    28 8
    14 2
    17 3
    17 5
    18 6
    20 5
    22 6
    11 4
    15 7
    2 0
    11 4
    11 3
    15 5
    7 1
    10 4
    16 5
    10 3
    16 6
    8 3
    30 12
    22 6
    28 9
    13 5
    16 6
    27 8
    12 4
    10 3
    16 7
    32 10
    6 3
    18 6
    18 4
    17 5
    16 7
    12 5
    7 2
    27 11
    24 8
    27 9
    20 5
    33 10
    32 11
    21 7
    17 4
    20 4
    13 3
    15 7
    23 6
    6 1
    10 1
    6 2
    11 5
    14 4
    9 2
    13 6
    18 4
    25 8
    23 9
    14 4
    12 4
    14 4
    8 4
    11 4
    30 10
    21 5
    18 5
    11 2
    11 4
    9 2
    19 5
    14 2
    16 3
    20 7
    12 6
    23 8
    14 4
    8 2
    21 7
    18 6
    25 9
    10 5
    16 4
    9 1
    20 6
    25 9
    18 5
    10 5
    5 2
    11 3
    8 4
    17 8
    18 7
    17 5
    13 5
    15 3
    15 3
    22 6
    14 4
    7 0
    12 2
    11 3
    22 8
    9 4
    24 8
    14 5
    24 7
    31 10
    14 4
    8 3
    10 4
    10 2
    34 12
    16 5
    20 6
    12 6
    6 3
    29 9
    13 3
    18 7
    14 6
    14 6
    12 4
    10 2
    13 3
    13 6
    13 2
    13 4
    18 6
    16 4
    23 9
    22 7
    13 6
    11 5
    21 5
    19 6
    3 2
    26 9
    18 6
    20 6
    27 9
    12 6
    15 3
    22 7
    10 4
    6 3
    20 7
    6 2
    11 2
    17 7
    25 8
    6 3
    20 7
    7 3
    25 6
    13 5
    28 9
    12 2
    10 3
    12 4
    26 9
    16 5
    14 3
    20 7
    8 2
    8 1
    15 3
    17 6
    10 3
    14 6
    27 9
    7 3
    6 1
    14 5
    14 5
    15 7
    28 7
    9 2
    10 6
    15 5
    11 4
    20 9
    6 2
    17 4
    25 5
    10 3
    17 2
    13 4
    15 5
    9 3
    11 4
    10 2
    13 5
    29 11
    15 3
    20 8
    11 4
    21 8
    11 3
    30 11
    29 7
    9 2
    9 4
    27 8
    15 5
    20 7
    7 3
    17 3
    26 9
    20 7
    20 6
    15 3
    18 3
    9 0
    21 7
    15 4
    15 4
    12 5
    22 8
    12 5
    11 3
    19 6
    19 4
    25 8
    13 3
    28 10
    19 5
    20 5
    17 6
    8 3
    20 5
    20 5
    24 9
    14 5
    18 7
    15 4
    9 1
    22 6
    17 5
    21 6
    22 8
    19 5
    13 3
    12 5
    8 2
    10 1
    12 5
    14 2
    8 2
    10 4
    14 9
    8 1
    19 8
    15 8
    9 1
    15 5
    23 6
    10 4
    11 6
    14 4
    5 1
    13 2
    14 5
    11 4
    7 1
    16 6
    10 3
    28 9
    24 9
    17 7
    16 6
    14 3
    14 4
    22 7
    10 3
    13 3
    25 8
    28 8
    28 9
    21 6
    18 7
    27 9
    22 7
    19 8
    11 2
    26 9
    18 4
    13 3
    14 6
    28 9
    16 5
    11 2
    14 3
    12 4
    7 2
    18 3
    13 1
    13 3
    13 1
    20 7
    8 3
    10 4
    13 4
    11 5
    11 5
    11 4
    24 8
    14 3
    15 4
    8 3
    8 2
    20 5
    12 3
    15 6
    24 6
    9 2
    13 5
    28 9
    11 3
    11 5
    12 4
    12 2
    14 4
    10 4
    23 7
    22 6
    21 6
    11 4
    12 3
    9 1
    19 5
    28 9
    14 4
    12 2
    8 3
    20 8
    13 4
    23 7
    13 4
    7 2
    18 6
    14 6
    14 5
    24 8
    11 3
    12 3
    15 5
    7 2
    16 6
    26 8
    11 2
    15 6
    8 3
    21 5
    23 7
    18 5
    27 8
    17 7
    21 8
    11 2
    13 4
    17 4
    20 5
    29 11
    28 8
    32 11
    20 7
    16 5
    27 11
    16 6
    28 9
    21 8
    20 5
    9 1
    12 4
    27 6
    22 8
    17 4
    25 10
    12 2
    11 5
    25 7
    18 6
    16 3
    16 7
    27 6
    27 7
    12 2
    15 6
    20 7
    5 1
    17 4
    26 10
    19 6
    17 6
    12 6
    11 4
    20 8
    32 10
    12 5
    13 4
    24 8
    21 7
    17 6
    15 3
    21 6
    16 4
    13 5
    8 1
    12 5
    12 4
    27 10
    9 4
    19 5
    9 2
    13 3
    12 4
    19 5
    31 10
    10 3
    10 3
    30 8
    20 9
    17 5
    11 4
    6 3
    28 10
    21 8
    11 3
    10 4
    23 7
    18 6
    10 2
    13 4
    20 8
    25 9
    25 7
    16 5
    10 2
    32 11
    7 2
    34 11
    15 4
    14 3
    27 8
    28 8
    8 5
    7 4
    16 5
    22 6
    22 6
    9 4
    14 3
    24 7
    20 8
    18 5
    23 6
    8 3
    9 6
    8 1
    18 6
    9 3
    17 6
    9 4
    10 1
    31 11
    16 5
    29 8
    27 10
    12 5
    28 8
    23 8
    30 11
    6 2
    13 3
    12 4
    24 7
    18 7
    12 2
    15 6
    5 2
    17 6
    12 4
    11 7
    15 4
    10 4
    12 6
    24 9
    34 11
    8 3
    30 10
    15 5
    20 7
    12 4
    16 5
    5 2
    18 5
    11 2
    9 3
    21 7
    10 3
    21 6
    4 2
    7 2
    9 2
    11 2
    12 4
    9 4
    7 2
    15 4
    11 5
    8 3
    13 4
    24 7
    27 10
    12 5
    10 5
    4 1
    10 4
    13 4
    12 4
    9 3
    17 6
    24 7
    20 7
    25 5
    12 5
    14 4
    17 4
    22 7
    13 4
    21 4
    20 7
    10 3
    23 6
    10 2
    13 4
    10 4
    26 8
    9 0
    9 5
    11 4
    26 8
    16 6
    29 10
    21 5
    17 3
    15 6
    7 4
    20 7
    12 2
    8 2
    21 6
    25 8
    12 3
    10 4
    25 9
    16 6
    14 3
    18 7
    19 7
    25 7
    15 5
    11 5
    12 3
    11 6
    9 2
    14 7
    25 7
    17 5
    15 7
    12 5
    17 6
    19 8
    23 7
    12 2
    14 5
    13 4
    13 6
    10 5
    12 5
    26 9
    7 4
    26 7
    12 3
    20 7
    32 11
    14 4
    13 6
    9 1
    8 3
    10 4
    14 4
    7 3
    14 3
    20 6
    7 1
    18 8
    8 4
    17 7
    11 4
    12 3
    11 5
    12 4
    13 4
    8 4
    16 4
    20 4
    22 9
    12 3
    9 1
    22 6
    30 10
    11 2
    8 3
    15 4
    11 4
    9 4
    9 3
    8 3
    16 7
    14 3
    22 7
    22 7
    6 1
    17 6
    22 7
    17 5
    18 6
    8 5
    15 5
    15 5
    16 6
    17 5
    32 11
    20 8
    12 7
    8 3
    24 9
    29 9
    8 6
    7 2
    27 9
    20 7
    22 6
    18 5
    11 3
    20 6
    6 2
    11 3
    10 4
    11 4
    15 5
    23 7
    16 3
    22 6
    23 7
    26 8
    17 6
    18 6
    18 8
    19 6
    24 9
    15 6
    16 6
    30 10
    11 4
    15 6
    12 5
    5 3
    25 7
    11 4
    10 4
    10 4
    24 8
    12 4
    8 2
    13 4
    17 5
    24 8
    14 6
    18 4
    21 6
    12 2
    8 2
    13 4
    4 1
    11 4
    12 5
    24 8
    33 12
    11 3
    16 5
    19 7
    13 3
    12 6
    13 6
    13 5
    20 8
    16 4
    25 9
    14 3
    21 8
    29 11
    29 10
    23 9
    14 5
    19 5
    22 7
    14 3
    6 1
    8 4
    17 5
    17 4
    20 5
    10 2
    12 4
    10 2
    9 6
    7 1
    18 8
    18 6
    9 3
    10 3
    15 5
    7 0
    13 6
    16 6
    23 8
    7 2
    6 2
    16 5
    16 6
    12 6
    12 3
    14 3
    16 5
    16 3
    9 2
    25 8
    10 3
    15 3
    15 4
    10 3
    13 5
    31 9
    19 6
    19 8
    25 8
    11 3
    25 8
    6 3
    20 8
    29 9
    16 6
    24 7
    5 0
    14 3
    22 7
    3 0
    17 5
    9 3
    13 5
    17 5
    9 3
    19 4
    23 7
    16 4
    26 8
    20 7
    14 4
    12 3
    29 9
    27 9
    7 2
    9 4
    22 8
    9 2
    14 6
    13 3
    22 7
    14 5
    19 6
    13 3
    17 6
    9 3
    25 10
    15 4
    14 6
    10 5
    30 10
    8 3
    19 4
    31 11
    9 3
    8 3
    9 2
    5 1
    29 10
    18 6
    11 3
    28 9
    5 2
    7 2
    22 7
    13 5
    17 10
    28 9
    13 5
    12 3
    26 8
    15 4
    10 4
    13 4
    14 5
    17 5
    15 3
    17 5
    11 2
    29 10
    11 2
    25 9
    11 3
    12 6
    21 10
    22 7
    23 9
    7 1
    14 6
    28 10
    20 5
    15 6
    9 2
    17 6
    16 3
    21 7
    9 3
    16 7
    16 6
    21 7
    5 3
    21 7
    27 9
    15 2
    16 5
    12 4
    7 3
    4 3
    10 4
    12 4
    18 6
    11 4
    12 5
    16 4
    16 6
    23 8
    25 8
    21 7
    16 6
    29 8
    6 3
    11 3
    19 4
    20 6
    12 4
    23 10
    26 8
    4 2
    22 6
    6 2
    10 5
    19 5
    12 5
    23 5
    11 2
    11 4
    14 6
    11 5
    11 5
    11 4
    11 4
    14 3
    14 5
    16 4
    18 6
    12 3
    23 7
    12 5
    8 1
    9 3
    18 6
    7 4
    12 4
    9 2
    23 6
    5 2
    10 3
    16 6
    7 3
    22 7
    11 2
    16 4
    10 1
    9 3
    16 5
    17 6
    21 8
    6 3
    13 4
    9 1
    10 4
    21 7
    16 4
    4 1
    10 5
    25 7
    12 3
    13 3
    13 5
    4 3
    17 7
    25 7
    13 4
    7 0
    14 6
    12 4
    22 5
    17 5
    24 8
    12 2
    9 0
    13 5
    17 4
    15 5
    17 4
    7 3
    19 7
    30 10
    17 6
    16 4
    15 4
    22 8
    20 6
    24 6
    10 3
    19 6
    24 7
    8 3
    21 9
    14 3
    9 2
    20 7
    9 5
    14 7
    17 5
    16 5
    7 0
    21 7
    25 8
    10 3
    14 4
    17 4
    7 2
    12 6
    6 3
    31 11
    16 5
    19 5
    13 3
    20 6
    18 6
    18 4
    12 3
    22 6
    13 3
    14 3
    24 10
    9 2
    17 7
    16 7
    26 9
    18 5
    14 4
    15 5
    15 6
    11 4
    17 5
    23 7
    6 2
    13 4
    14 5
    20 4
    9 1
    22 9
    13 3
    25 8
    28 7
    34 12
    9 1
    16 3
    24 8
    14 4
    14 3
    10 3
    13 6
    21 7
    22 6
    20 4
    6 2
    9 4
    12 4
    16 7
    11 5
    17 5
    4 2
    16 6
    18 6
    13 4
    18 3
    6 3
    9 3
    8 4
    8 4
    8 2
    17 5
    30 10
    16 6
    14 4
    14 7
    25 8
    16 4
    7 1
    23 9
    18 6
    14 4
    21 6
    12 3
    13 4
    10 2
    27 9
    11 3
    15 6
    10 3
    21 7
    11 3
    12 4
    22 7
    11 4
    26 8
    18 4
    23 7
    15 3
    16 6
    12 2
    26 9
    9 3
    6 3
    18 5
    11 6
    10 6
    15 6
    6 3
    3 1
    26 8
    21 8
    21 6
    18 4
    14 3
    12 4
    8 2
    17 5
    12 3
    11 3
    14 4
    8 4
    20 7
    32 11
    11 2
    6 2
    11 3
    9 4
    19 6
    21 5
    17 6
    22 7
    13 3
    27 9
    12 4
    8 4
    15 4
    9 3
    7 1
    3 1
    19 6
    15 8
    20 8
    8 3
    20 7
    7 2
    12 5
    14 4
    15 4
    17 7
    12 3
    14 4
    7 2
    19 6
    12 5
    11 5
    22 7
    17 5
    13 5
    15 7
    11 4
    17 8
    23 6
    16 2
    17 6
    27 9
    10 4
    29 9
    24 9
    7 2
    21 7
    6 3
    12 6
    14 6
    18 5
    13 5
    21 9
    20 6
    11 3
    28 9
    8 3
    15 6
    33 10
    10 4
    24 7
    26 7
    18 4
    12 5
    26 8
    19 9
    16 6
    14 5
    10 3
    8 1
    8 3
    7 4
    9 2
    28 8
    10 4
    30 10
    7 4
    19 8
    8 3
    11 6
    8 4
    16 8
    16 6
    16 3
    22 8
    23 8
    8 2
    11 3
    11 5
    29 9
    16 6
    12 3
    12 6
    23 7
    11 1
    9 4
    9 3
    14 3
    6 2
    11 3
    22 6
    10 4
    6 1
    13 2
    15 5
    25 8
    12 4
    9 3
    18 7
    28 9
    12 3
    25 7
    17 7
    14 5
    14 6
    19 6
    33 11
    12 4
    15 4
    21 7
    18 7
    19 6
    21 7
    8 3
    26 8
    14 6
    19 7
    25 7
    11 2
    11 5
    13 4
    13 6
    18 5
    20 5
    23 9
    17 5
    10 4
    7 1
    27 9
    7 0
    12 4
    15 7
    11 3
    8 3
    16 7
    22 8
    26 9
    14 3
    17 6
    12 4
    11 5
    16 4
    7 2
    11 3
    18 6
    25 9
    24 8
    23 10
    11 5
    13 4
    23 8
    13 5
    8 3
    23 7
    14 4
    19 6
    10 3
    12 2
    22 5
    5 0
    26 8
    14 3
    10 3
    34 12
    14 3
    8 3
    23 5
    32 10
    19 7
    26 8
    13 5
    10 3
    6 4
    15 6
    10 4
    9 4
    12 3
    13 6
    18 7
    21 7
    14 4
    18 8
    14 6
    16 6
    29 7
    15 5
    12 6
    10 4
    20 6
    21 7
    18 6
    16 4
    11 7
    20 5
    29 8
    11 4
    14 5
    11 4
    16 5
    27 11
    12 5
    12 3
    25 9
    16 6
    11 6
    7 3
    19 7
    17 5
    9 3
    18 6
    17 4
    22 7
    19 6
    15 5
    8 2
    13 4
    10 3
    18 7
    26 9
    11 2
    9 3
    22 10
    17 8
    14 5
    35 11
    11 2
    26 7
    12 5
    16 5
    16 5
    29 11
    12 5
    17 4
    15 7
    12 3
    33 11
    14 5
    9 2
    9 1
    23 7
    13 5
    17 6
    11 6
    8 3
    8 3
    27 8
    27 9
    10 3
    14 5
    9 3
    22 5
    18 8
    23 6
    12 5
    7 3
    10 2
    32 11
    24 9
    11 2
    12 2
    15 6
    11 4
    8 4
    30 10
    11 4
    16 6
    8 2
    14 5
    10 2
    13 2
    13 3
    10 4
    18 5
    12 3
    12 2
    8 2
    10 1
    8 2
    9 3
    14 4
    14 5
    27 8
    10 3
    12 4
    19 8
    11 4
    17 5
    16 7
    24 7
    29 8
    10 4
    18 6
    13 3
    9 3
    6 2
    26 7
    15 6
    11 6
    13 5
    6 3
    9 1
    23 7
    14 2
    18 3
    34 11
    12 5
    12 4
    9 3
    17 7
    12 6
    18 7
    10 3
    21 7
    21 6
    20 7
    17 4
    18 4
    16 5
    12 4
    32 12
    14 4
    17 5
    17 6
    19 5
    14 7
    15 7
    13 5
    32 11
    11 3
    12 3
    22 9
    6 1
    7 3
    26 8
    22 8
    11 3
    30 10
    25 8
    11 3
    8 2
    17 6
    21 4
    15 4
    12 5
    18 5
    17 4
    14 4
    6 1
    6 3
    21 4
    9 3
    13 5
    14 5
    23 8
    22 6
    16 2
    12 4
    21 5
    18 5
    17 7
    14 3
    15 4
    16 5
    18 5
    15 5
    16 7
    12 3
    19 5
    9 3
    9 2
    12 6
    24 8
    8 3
    33 12
    9 3
    10 3
    16 4
    18 5
    12 2
    14 5
    19 7
    5 1
    13 4
    6 2
    10 2
    24 7
    15 5
    18 4
    14 3
    9 3
    15 5
    30 9
    11 5
    8 4
    3 0
    9 2
    25 8
    9 4
    19 6
    17 6
    11 6
    17 6
    9 4
    11 3
    14 7
    9 3
    14 6
    13 5
    15 5
    12 3
    18 6
    14 6
    22 7
    27 6
    13 3
    13 4
    20 6
    17 7
    9 4
    16 6
    21 5
    32 10
    25 6
    27 8
    23 7
    13 5
    24 7
    24 8
    8 1
    21 7
    12 3
    26 8
    8 4
    14 4
    21 6
    33 11
    22 6
    10 2
    16 5
    20 6
    29 9
    23 7
    24 7
    16 6
    12 6
    30 10
    9 4
    11 1
    14 3
    19 7
    14 5
    19 6
    20 7
    10 4
    12 5
    22 8
    14 5
    20 7
    8 2
    8 4],
    ["total", "anchor"]
)

KBneatY = DataFrame(
    [29 9
    16 5
    18 6
    24 7
    19 5
    17 5
    15 6
    15 3
    15 5
    19 8
    16 4
    8 1
    22 7
    21 5
    3 0
    12 4
    10 3
    25 5
    20 6
    25 9
    18 6
    14 5
    32 10
    26 7
    11 6
    12 3
    12 4
    13 4
    12 3
    16 6
    32 11
    6 1
    16 5
    27 9
    12 4
    26 6
    18 6
    10 4
    35 12
    16 5
    5 3
    16 6
    29 8
    26 8
    33 10
    17 7
    14 3
    18 5
    31 10
    13 4
    32 10
    11 3
    14 4
    30 10
    31 10
    11 2
    11 3
    19 6
    18 3
    10 3
    23 7
    20 5
    33 10
    11 3
    16 3
    22 7
    15 5
    19 6
    9 3
    19 4
    18 7
    9 4
    22 3
    15 6
    18 5
    15 7
    18 7
    19 4
    16 6
    9 2
    22 4
    18 6
    28 9
    9 4
    20 8
    10 5
    26 10
    33 11
    9 3
    26 6
    9 3
    23 6
    14 4
    17 5
    5 2
    8 1
    23 10
    19 7
    17 5
    18 6
    21 7
    18 5
    24 8
    12 5
    7 1
    13 3
    25 9
    16 6
    23 8
    20 8
    35 12
    23 6
    34 11
    22 7
    21 5
    11 4
    11 4
    22 7
    15 4
    15 4
    32 12
    15 4
    17 5
    22 7
    18 4
    28 9
    21 6
    28 9
    25 8
    6 3
    24 8
    24 8
    13 3
    20 5
    16 7
    15 3
    16 7
    12 3
    23 9
    10 2
    21 7
    20 8
    15 5
    13 4
    14 4
    30 10
    19 8
    12 4
    20 6
    20 4
    31 10
    18 6
    23 6
    9 3
    25 7
    11 3
    6 2
    31 9
    22 4
    27 7
    16 5
    15 4
    25 7
    19 5
    31 9
    21 7
    24 8
    19 7
    12 5
    18 5
    13 6
    19 5
    6 1
    35 11
    15 5
    27 7
    23 5
    17 6
    15 5
    19 7
    22 7
    11 3
    17 7
    36 12
    13 6
    24 7
    16 5
    12 4
    14 4
    27 8
    23 7
    25 9
    29 9
    23 6
    11 3
    29 9
    11 4
    22 5
    19 7
    5 1
    18 5
    25 7
    27 9
    22 7
    17 4
    26 7
    20 7
    20 7
    23 5
    19 5
    21 10
    20 7
    17 6
    16 6
    19 5
    15 4
    10 1
    9 2
    24 7
    27 8
    13 6
    12 4
    10 3
    15 4
    20 5
    12 2
    19 5
    24 8
    9 2
    13 4
    31 10
    27 8
    29 10
    9 2
    15 2
    31 10
    23 10
    23 9
    13 4
    21 7
    11 2
    14 3
    15 4
    10 3
    35 12
    14 7
    14 5
    23 7
    19 3
    10 3
    8 3
    26 8
    16 7
    17 2
    27 9
    12 5
    17 6
    20 5
    27 9
    24 7
    21 7
    11 4
    20 7
    15 5
    18 4
    13 2
    24 7
    16 4
    11 5
    29 11
    20 7
    15 5
    13 4
    6 2
    20 6
    18 7
    12 5
    18 7
    11 2
    16 3
    14 2
    26 8
    20 5
    16 4
    9 3
    18 5
    20 6
    28 8
    6 2
    22 7
    10 1
    32 10
    7 0
    27 9
    26 10
    31 9
    30 10
    12 5
    26 8
    9 2
    21 7
    8 3
    11 4
    24 7
    17 7
    17 5
    28 8
    12 7
    27 8
    24 8
    12 3
    19 4
    11 4
    28 10
    12 5
    31 11
    10 3
    34 11
    26 8
    20 4
    26 8
    16 4
    13 3
    16 7
    18 3
    17 8
    17 7
    17 5
    27 9
    13 6
    23 9
    10 4
    11 3
    9 5
    19 7
    19 6
    21 4
    17 5
    20 7
    21 5
    20 7
    33 9
    14 5
    30 9
    25 7
    17 6
    20 7
    26 6
    23 8
    15 6
    10 4
    14 4
    18 4
    16 4
    12 1
    30 8
    16 4
    33 11
    12 6
    22 7
    23 9
    26 9
    28 10
    23 8
    27 10
    28 9
    32 9
    31 11
    18 5
    30 10
    17 3
    13 4
    17 5
    14 4
    11 4
    13 4
    6 2
    12 5
    29 9
    16 5
    20 7
    21 6
    23 7
    16 7
    22 6
    21 6
    21 7
    13 6
    21 6
    14 5
    23 8
    13 4
    8 1
    20 6
    20 6
    8 2
    17 7
    16 4
    20 7
    15 5
    15 8
    15 5
    18 6
    15 5
    18 6
    26 10
    19 6
    16 3
    12 2
    19 8
    23 6
    16 5
    28 9
    4 2
    13 4
    13 3
    31 9
    24 8
    13 1
    26 7
    21 8
    11 5
    28 9
    32 12
    20 7
    13 5
    19 8
    17 4
    17 6
    23 7
    25 7
    22 6
    17 4
    17 5
    28 9
    17 6
    15 4
    16 5
    29 10
    22 6
    24 8
    15 4
    17 5
    32 9
    22 9
    18 5
    16 6
    18 6
    14 4
    23 8
    19 7
    12 4
    33 11
    14 2
    17 5
    16 7
    28 7
    21 7
    20 4
    20 8
    10 4
    21 7
    26 7
    14 5
    12 3
    21 6
    26 8
    18 4
    19 6
    19 4
    14 5
    17 7
    26 9
    18 7
    12 6
    22 6
    25 8
    18 3
    24 10
    11 2
    17 3
    23 6
    24 7
    19 6
    14 7
    15 4
    23 3
    22 6
    16 3
    30 9
    19 4
    22 6
    21 5
    30 10
    14 7
    26 8
    23 8
    24 7
    30 9
    23 6
    18 6
    16 8
    13 6
    23 9
    26 10
    6 1
    12 3
    23 8
    24 9
    21 6
    17 3
    23 7
    23 7
    25 9
    16 4
    15 5
    17 5
    14 4
    16 5
    11 5
    11 3
    18 6
    22 6
    26 8
    16 5
    14 5
    11 5
    23 9
    14 4
    20 7
    14 5
    30 9
    26 7
    12 4
    23 10
    23 6
    29 9
    26 8
    33 10
    14 3
    9 2
    7 2
    11 3
    13 4
    13 5
    21 7
    30 9
    25 7
    14 6
    11 4
    14 5
    24 6
    24 9
    16 7
    14 3
    21 6
    10 3
    19 8
    25 9
    18 4
    20 7
    18 6
    14 3
    13 5
    30 10
    16 5
    20 7
    29 9
    10 3
    33 10
    18 5
    22 8
    25 7
    10 3
    30 11
    12 3
    33 10
    5 1
    14 3
    13 5
    22 8
    17 8
    16 5
    12 2
    15 5
    15 6
    19 7
    21 5
    14 6
    28 7
    14 4
    23 7
    8 2
    19 6
    33 10
    27 10
    27 7
    7 1
    12 3
    19 6
    26 9
    28 9
    17 6
    23 7
    31 9
    12 4
    17 6
    6 3
    24 8
    15 4
    30 10
    26 6
    19 4
    12 3
    10 3
    33 11
    20 8
    22 7
    21 6
    15 5
    24 7
    14 6
    7 4
    20 6
    11 4
    26 6
    17 8
    30 9
    17 4
    22 4
    31 9
    9 3
    19 5
    18 6
    15 5
    17 8
    25 8
    13 5
    24 9
    7 2
    19 8
    12 5
    25 8
    16 5
    19 6
    10 6
    8 1
    19 5
    27 10
    25 8
    27 8
    11 4
    6 1
    25 8
    6 2
    6 1
    17 6
    22 7
    20 5
    30 10
    29 10
    8 2
    20 5
    12 5
    13 6
    10 5
    14 5
    20 8
    19 6
    18 5
    18 6
    16 2
    13 6
    9 2
    14 5
    15 5
    27 8
    18 3
    27 8
    21 4
    9 4
    22 9
    20 6
    6 1
    15 6
    20 8
    9 5
    9 3
    13 3
    23 8
    23 6
    11 3
    15 7
    21 6
    18 4
    28 9
    5 1
    10 3
    14 3
    16 3
    22 7
    15 4
    14 4
    26 8
    31 10
    7 3
    29 7
    16 6
    12 4
    20 6
    27 10
    16 4
    10 2
    11 3
    21 7
    25 7
    18 4
    8 0
    14 5
    15 4
    30 9
    29 9
    17 7
    13 4
    30 8
    23 6
    23 8
    13 6
    24 9
    11 3
    18 6
    25 7
    29 10
    18 4
    25 8
    27 9
    6 3
    19 4
    26 9
    7 2
    23 8
    21 8
    22 7
    23 7
    12 4
    14 6
    15 4
    21 6
    13 3
    17 7
    13 3
    14 4
    14 5
    28 8
    25 8
    14 4
    9 2
    22 6
    17 5
    23 7
    20 7
    18 4
    25 10
    19 7
    12 3
    10 6
    22 7
    23 4
    24 5
    19 8
    14 5
    25 10
    20 7
    25 8
    34 12
    14 4
    22 6
    17 7
    27 10
    26 8
    30 12
    15 5
    9 4
    30 8
    18 6
    12 3
    32 9
    27 8
    33 10
    24 9
    13 5
    7 1
    23 9
    14 3
    17 6
    22 8
    20 7
    23 5
    11 3
    20 5
    15 6
    21 7
    33 11
    31 11
    27 8
    17 3
    16 6
    14 4
    19 6
    13 3
    10 4
    14 4
    13 5
    17 5
    18 6
    22 6
    27 7
    19 5
    27 8
    7 2
    27 7
    14 4
    29 10
    25 9
    15 5
    19 4
    29 10
    12 3
    31 10
    13 5
    21 7
    10 0
    8 2
    8 5
    11 4
    30 11
    27 10
    7 3
    16 5
    12 3
    12 3
    14 7
    33 12
    20 8
    25 6
    16 6
    26 7
    9 2
    19 4
    14 3
    15 5
    7 4
    20 5
    21 7
    10 5
    22 7
    22 7
    18 4
    8 1
    15 5
    28 9
    24 9
    12 5
    12 5
    30 9
    13 5
    16 6
    34 10
    20 6
    8 5
    24 9
    11 3
    15 4
    24 10
    21 7
    18 6
    21 7
    14 3
    15 5
    8 3
    25 8
    11 5
    14 6
    25 8
    22 5
    12 3
    18 7
    25 8
    6 2
    19 6
    15 5
    7 5
    17 7
    10 2
    20 9
    11 3
    17 5
    11 6
    19 6
    17 3
    22 6
    17 5
    15 5
    25 8
    23 7
    17 5
    22 6
    9 2
    29 10
    23 5
    10 4
    9 1
    17 3
    19 6
    33 11
    20 7
    16 2
    10 1
    12 4
    14 6
    27 8
    36 12
    29 8
    27 9
    23 7
    23 7
    9 3
    11 2
    21 7
    30 9
    30 10
    11 6
    12 4
    21 7
    24 6
    20 5
    25 7
    25 8
    11 4
    24 7
    23 7
    12 4
    25 8
    15 5
    22 6
    4 0
    11 6
    7 2
    8 2
    12 2
    18 7
    23 7
    17 5
    6 3
    26 9
    18 4
    28 7
    23 7
    21 6
    19 5
    14 7
    25 6
    7 1
    30 11
    11 2
    29 9
    15 5
    12 4
    10 5
    29 8
    24 7
    32 10
    16 5
    23 8
    15 6
    9 2
    27 8
    17 7
    17 5
    12 4
    17 4
    15 4
    27 9
    19 7
    14 4
    18 5
    16 4
    16 4
    12 3
    22 7
    19 6
    21 6
    8 2
    8 1
    21 5
    14 4
    7 2
    8 2
    19 5
    8 4
    12 4
    14 4
    11 4
    10 5
    10 4
    30 9
    8 3
    25 7
    19 6
    28 7
    16 5
    21 5
    12 5
    14 3
    28 8
    17 7
    17 6
    26 7
    16 6
    14 3
    16 5
    10 4
    17 6
    30 9
    20 6
    13 3
    8 0
    15 6
    24 5
    22 7
    20 7
    13 4
    15 7
    18 6
    14 4
    10 2
    28 9
    33 11
    23 7
    21 5
    29 10
    27 8
    16 6
    20 6
    27 10
    18 7
    29 10
    11 3
    28 7
    20 6
    17 4
    13 4
    6 2
    4 0
    17 7
    14 5
    15 3
    10 4
    10 4
    14 4
    15 6
    10 3
    24 6
    19 4
    10 4
    19 6
    13 4
    11 3
    31 7
    25 8
    21 6
    30 8
    19 6
    32 12
    25 9
    35 11
    31 9
    18 4
    17 5
    18 5
    22 6
    16 5
    24 7
    22 7
    8 3
    15 3
    24 7
    11 2
    30 10
    17 6
    24 7
    7 3
    19 7
    28 9
    18 4
    5 0
    33 11
    6 1
    16 6
    7 4
    29 9
    24 9
    20 7
    24 8
    33 9
    33 12
    15 6
    12 6
    11 5
    15 7
    16 3
    12 3
    29 9
    13 4
    28 8
    21 6
    24 8
    16 6
    17 7
    10 4
    19 3
    30 9
    20 7
    29 9
    17 6
    15 5
    24 7
    32 9
    14 5
    16 6
    12 3
    13 2
    33 9
    18 4
    21 6
    24 7
    17 6
    28 8
    13 3
    24 6
    24 6
    16 5
    23 7
    23 8
    18 6
    14 3
    29 11
    27 8
    12 4
    18 5
    9 2
    14 4
    13 5
    12 5
    19 7
    20 5
    15 5
    31 10
    19 6
    24 7
    19 7
    25 9
    11 4
    19 5
    13 3
    15 5
    16 5
    13 4
    22 8
    10 1
    15 5
    19 4
    6 3
    25 7
    7 4
    27 8
    28 9
    11 5
    18 6
    16 4
    13 2
    8 3
    21 7
    30 9
    33 10
    22 7
    17 7
    22 6
    12 4
    15 8
    12 2
    22 6
    13 2
    11 5
    17 6
    15 8
    11 4
    12 3
    13 5
    32 8
    13 5
    23 7
    27 8
    25 8
    25 5
    32 11
    16 4
    28 8
    6 2
    16 2
    26 8
    19 6
    28 10
    33 11
    24 10
    7 2
    10 1
    29 10
    25 7
    27 8
    14 6
    14 6
    26 7
    22 6
    33 9
    21 8
    23 7
    18 7
    14 7
    16 3
    15 3
    17 6
    17 5
    29 8
    26 7
    17 7
    22 6
    30 9
    5 0
    27 8
    27 9
    17 5
    14 4
    21 7
    29 10
    15 4
    19 8
    12 5
    19 5
    29 8
    20 5
    17 5
    15 4
    18 5
    11 6
    27 8
    19 7
    33 11
    14 4
    10 4
    17 6
    12 4
    14 6
    13 4
    24 8
    15 4
    29 11
    10 2
    20 7
    17 6
    9 4
    12 4
    14 5
    17 5
    11 4
    14 4
    23 9
    27 8
    16 5
    17 4
    21 8
    16 6
    12 7
    29 9
    27 7
    17 3
    15 4
    21 7
    23 6
    24 7
    8 4
    24 8
    31 9
    10 3
    29 10
    10 4
    17 5
    8 3
    19 5
    19 6
    31 10
    19 5
    20 5
    27 8
    24 6
    23 8
    25 8
    11 3
    23 7
    22 5
    25 8
    6 1
    30 9
    23 7
    10 4
    10 4
    14 5
    15 5
    20 7
    29 9
    9 2
    12 3
    16 5
    29 10
    21 7
    17 4
    20 7
    9 2
    18 8
    17 6
    17 5
    13 5
    22 7
    16 7
    21 7
    20 5
    10 2
    27 6
    21 8
    10 3
    20 6
    15 6
    9 3
    18 5
    23 6
    25 9
    17 6
    15 5
    26 9
    23 8
    25 7
    22 7
    24 7
    21 7
    16 5
    26 7
    25 7
    19 8
    15 3
    27 9
    25 8
    13 5
    8 4
    12 4
    15 6
    22 7
    19 6
    27 8
    12 3
    20 7
    26 9
    27 10
    15 6
    25 8
    16 7
    12 2
    18 6
    23 8
    34 12
    24 8
    9 0
    15 5
    25 7
    22 5
    9 2
    17 4
    24 9
    31 10
    13 3
    20 4
    15 6
    17 2
    24 8
    9 3
    22 7
    12 3
    17 6
    19 5
    17 5
    20 6
    9 4
    27 9
    26 8
    10 4
    6 2
    28 8
    13 6
    14 5
    32 11
    11 5
    20 7
    29 11
    9 4
    19 7
    12 5
    21 6
    12 6
    18 6
    7 2
    24 9
    5 2
    15 4
    28 10
    26 9
    16 4
    9 5
    13 3
    18 3
    21 5
    10 2
    9 3
    17 4
    15 7
    14 5
    7 1
    20 6
    6 2
    18 3
    17 5
    11 5
    25 8
    20 6
    8 2
    10 4
    33 10
    10 3
    29 9
    17 7
    14 5
    19 4
    23 9
    23 10
    18 5
    23 7
    20 6
    18 5
    17 7
    28 11
    8 2
    33 10
    27 8
    18 6
    30 10
    25 7
    16 7
    17 6
    21 6
    19 7
    6 0
    20 8
    11 3
    10 1
    12 4
    22 6
    30 11
    8 1
    15 4
    13 5
    15 6
    15 2
    18 4
    10 3
    8 2
    21 7
    22 7
    30 10
    15 4
    19 5
    15 2
    23 7
    24 8
    26 9
    5 2
    11 3
    11 2
    16 6
    18 8
    12 4
    15 2
    12 3
    14 4
    14 6
    16 5
    12 3
    26 7
    22 7
    14 6
    10 3
    19 7
    17 5
    24 7
    19 7
    15 6
    24 7
    7 2
    15 3
    15 7
    30 9
    17 7
    9 2
    25 7
    18 5
    23 8
    27 9
    24 9
    21 8
    30 9
    21 8
    10 3
    27 9
    18 6
    17 6
    32 10
    18 5
    26 10
    5 4
    18 6
    26 8
    8 5
    12 3
    14 1
    20 6
    26 7
    26 9
    11 2
    13 6
    7 2
    13 4
    11 3
    22 5
    19 5
    21 6
    14 4
    10 4
    18 7
    16 6
    10 3
    9 2
    9 4
    8 3
    10 4
    13 5
    25 6
    34 11
    19 7
    8 4
    17 5
    25 8
    15 4
    15 5
    17 4
    14 4
    13 3
    17 6
    30 11
    12 3
    16 5
    18 7
    32 11
    13 6
    27 8
    22 7
    25 7
    30 10
    19 3
    22 7
    16 4
    26 8
    19 5
    25 8
    20 6
    22 7
    19 5
    22 9
    27 8
    11 3
    25 8
    31 10
    10 3
    16 3
    9 1
    21 5
    17 7
    28 10
    24 6
    24 8
    21 9
    17 5
    21 7
    32 10
    23 6
    16 6
    32 10
    17 5
    25 8
    26 8
    27 8
    32 11
    9 3
    33 11
    28 10
    24 7
    18 6
    27 8
    11 1
    23 7
    6 1
    22 5
    24 8
    8 3
    23 9
    31 10
    20 6
    11 3
    15 4
    14 6
    28 9
    14 6],
    ["total", "anchor"]
)

KBp48 = [
    .00,
    .01,
    .03,
    .08,
    .22,
    .53,
    1.42,
    2.88,
    4.70,
    7.42,
    10.80,
    14.74,
    19.17,
    23.61,
    28.15,
    32.83,
    37.69,
    42.28,
    46.50,
    50.59,
    54.87,
    58.89,
    62.47,
    66.05,
    69.37,
    72.68,
    75.84,
    78.60,
    81.06,
    83.44,
    85.90,
    88.18,
    90.19,
    92.00,
    93.67,
    95.37,
    96.81,
    97.82,
    98.68,
    99.39,
    99.83
]


text:
module TestDefineSingletons

using DefineSingletons
using Documenter: doctest
using Test

abstract type SuperType end
struct NonparametricSingleton end
struct ParametricSingleton{T} end

@def_singleton singleton_gentype_gensym
@def_singleton singleton_gentype_withname::GenTypeName
@def_singleton singleton_gentype_withsuper1::GenSubTypeName1 <: SuperType
@def_singleton singleton_gentype_withsuper2::(GenSubTypeName2 <: SuperType)
@def_singleton singleton_gentype_withsuper3 isa SuperType
@def_singleton singleton_predefined_nonparametric = NonparametricSingleton()
@def_singleton singleton_predefined_parametric = ParametricSingleton{1}()

issingleton(x) = Base.issingletontype(typeof(x))

@testset "DefineSingletons.jl" begin
    @test issingleton(singleton_gentype_gensym)
    @test issingleton(singleton_gentype_withname)
    @test issingleton(singleton_gentype_withsuper1)
    @test issingleton(singleton_gentype_withsuper2)
    @test issingleton(singleton_gentype_withsuper3)
    @test issingleton(singleton_predefined_nonparametric)
    @test issingleton(singleton_predefined_parametric)

    @test singleton_gentype_withname isa GenTypeName
    @test singleton_gentype_withsuper1 isa GenSubTypeName1
    @test singleton_gentype_withsuper2 isa GenSubTypeName2
    @test singleton_gentype_withsuper3 isa SuperType
    @test GenSubTypeName1 <: SuperType
    @test GenSubTypeName2 <: SuperType

    kw = (; context = :limit => true)
    @test repr(singleton_gentype_gensym; kw...) === "singleton_gentype_gensym"
    @test repr(singleton_gentype_withname; kw...) === "singleton_gentype_withname"
    @test repr(singleton_gentype_withsuper1; kw...) === "singleton_gentype_withsuper1"
    @test repr(singleton_gentype_withsuper2; kw...) === "singleton_gentype_withsuper2"
    @test repr(singleton_gentype_withsuper3; kw...) === "singleton_gentype_withsuper3"
    @test repr(singleton_predefined_nonparametric; kw...) ===
          "singleton_predefined_nonparametric"
    @test repr(singleton_predefined_parametric; kw...) === "singleton_predefined_parametric"

    @test endswith(repr(singleton_gentype_gensym), ".singleton_gentype_gensym")
    @test endswith(repr(singleton_gentype_withname), ".singleton_gentype_withname")
    @test endswith(repr(singleton_gentype_withsuper1), ".singleton_gentype_withsuper1")
    @test endswith(repr(singleton_gentype_withsuper2), ".singleton_gentype_withsuper2")
    @test endswith(repr(singleton_gentype_withsuper3), ".singleton_gentype_withsuper3")
    @test endswith(
        repr(singleton_predefined_nonparametric),
        ".singleton_predefined_nonparametric",
    )
    @test endswith(
        repr(singleton_predefined_parametric),
        ".singleton_predefined_parametric",
    )
end

@testset "doctest" begin
    doctest(DefineSingletons; manual = false)
end

end  # module


text:
# This file includes chaotic maps

export Logistic, Cubic, Tent, Henon, Bernoulli,
    Lorenz, Chua, Rossler,
    trajectory!, normalize

abstract type AbstractOscillator end
abstract type AbstractDiscreteOscillator <: AbstractOscillator  end
abstract type AbstractContinuousOscillator <: AbstractOscillator  end

# ----------------------------- Discrete Time Maps --------------------------------- #

struct Cubic <: AbstractDiscreteOscillator
    a::Float64
    b::Float64
    x::Vector{Float64}
    t::Float64
end
Cubic() = Cubic(4, 3, rand(1), 0.)
(cmap::Cubic)(dx, x, u, t) = ( dx[1] = cmap.a * x[1]^3 - cmap.b * x[1] )

struct Tent <: AbstractDiscreteOscillator
    a::Float64
    b::Float64
    c::Float64
    x::Vector{Float64}
    t::Float64
end
Tent() = Tent(0.6, 0.6, 0.4, rand(1), 0.)
(cmap::Tent)(dx, x, u, t) = ( dx[1] = 0 ≤ x[1] ≤ cmap.a ? x[1] / cmap.b : (1 - x[1]) / cmap.c )

struct Henon <: AbstractDiscreteOscillator
    a::Float64
    b::Float64
    x::Vector{Float64}
    t::Float64
end
Henon() = Henon(1.4, 0.3, rand(2), 0.)

function (cmap::Henon)(dx, x, u, t)
    dx[1] = 1 + x[2] - cmap.a * x[1]^2
    dx[2] = cmap.b * x[1]
end

struct Bernoulli <: AbstractDiscreteOscillator
    a::Float64
    x::Vector{Float64}
    t::Float64
end
Bernoulli() = Bernoulli(1.2, rand(1), 0.)
(cmap::Bernoulli)(dx, x, u, t) = ( dx[1] = x[1] < 0 ? cmap.a * x[1] + 1 : cmap.a * x[1] - 1 )

# ----------------------------- Continuous Time Maps --------------------------------- #

struct Lorenz <: AbstractContinuousOscillator
    σ::Float64
    β::Float64
    ρ::Float64
    γ::Float64
    x::Vector{Float64}
    t::Float64
end
Lorenz() = Lorenz(10, 8 / 3, 28, 1, rand(3), 0)

function (cmap::Lorenz)(dx, x, u, t)
    dx[1] = cmap.σ * (x[2] - x[1])
    dx[2] = x[1] * (cmap.ρ - x[3]) - x[2]
    dx[3] = x[1] * x[2] - cmap.β * x[3]
end

struct Diode
    a::Float64
    b::Float64
    bp::Float64
end
Diode() = Diode(-1.143, -0.714, 1.)
function (diode::Diode)(x)
    if x < -cmap.bp
        cmap.b * x + (cmap.b - cmap.a) * cmap.bp
    elseif -cmap.bp ≤ x ≤ cmap.bp
        cmap.a * x
    elseif x > cmap.bp
        cmap.b * x +  (cmap.a - cmap.b) * bp1
    end
end

struct Chua <: AbstractContinuousOscillator
    α::Float64
    β::Float64
    h::Diode
    x::Vector{Float64}
    t::Float64
end
Chua() = Chua(15, 28, Diode(), rand(3), 0.)

function (cmap::Chua)(dx, x, u, t)
    dx[1] = cmap.α * (x[2] - x[1] - cmap.h(x[1]))
    dx[2] = x[1] - x[2] + x[3]
    dx[3] = -β * x[2]
end

struct Rossler <: AbstractContinuousOscillator
    a::Float64
    b::Float64
    c::Float64
    x::Vector{Float64}
    t::Float64
end
Rossler() = Rossler(0.38, 0.3, 4.82, rand(3), 0.)

function (cmap::Rossler)(dx, x, u, t)
    dx[1] = -x[2] - x[3]
    dx[2] = x[1] + cmap.a * x[2]
    dx[3] = cmap.b + x[3] * (x[1] - cmap.c)
end

# --------------------------------------- Methods  -------------------------------------- #

"""
   $SIGNATURES

Normalizes (zero mean and unity variance)  `x`.
"""
normalize(x) = (x .- mean(x)) / std(x)

"""
    $SIGNATURES

Returns the dimension of the state space of `cmap`.
"""
statedim(cmap::AbstractOscillator) = length(cmap.x)

"""
    $SIGNATURES

Returns a trajectory! of `camp` for a time span of `trange`. `idx` is the indices of trajectory! to be returned.
"""
function trajectory! end

function trajectory!(cmap::AbstractDiscreteOscillator, trange, idx=1)
    sol = solve(DiscreteProblem(cmap, cmap.x, (cmap.t, cmap.t + trange)))
    cmap.t = sol.t[end]
    cmap.x = sol.u[end]
    map(i -> normalize(getindex.(sol.u, i)), idx)
end

function trajectory!(cmap::AbstractContinuousOscillator, trange, tsample, idx=1)
    sol = solve(ODEProblem(cmap, cmap.x, (cmap.t, cmap.t + trange)), saveat=tsample)
    cmap.t = sol.t[end]
    cmap.x = sol.u[end]
    map(i -> normalize(getindex.(sol.u, i)), idx)
end


text:
# Image segmentation experiments, but using sparse reduction + maxflow
using MAT
include("../src/hypergraphs.jl")
include("../src/SparseCard.jl")
include("image_graphs_utils.jl")

## Load optimal solution to the problem, for comparison

for dataset = ["020_smallplant","001_oct"]

for numsuperpix = [200 500]

if numsuperpix == 200
    lambda = 0.002
else
    lambda = 0.008
end

## Load edges from non-superpixel (clique) regions
G = matread("GraphPotentials/Graph_noclique_$(dataset)_$(numsuperpix).mat")
An = G["A"]
tvec = G["tvec"]
svec = G["svec"]
n = size(An,1)

In,Jn,Wn = findnz(An)

## Load hyperedges from superpixel (clique) regions

Edges, EdgesW = load_superpixel_cliques(dataset,lambda,numsuperpix)
K = maximum(length(e) for e in Edges)

## Sparse Reduction Runs

#epsilons = [1.0;0.5;0.1;0.05;0.01;0.005;0.001];0.0005;0.0001]
small = 4
epsilons = 10 .^ range( 0, -small, length = 20)
runs = length(epsilons)
objvals = zeros(runs)
approxcuts = zeros(runs)
runtimes = zeros(runs)
numedges = zeros(runs)
preptime = zeros(runs)
println("")
for i = 1:runs
    epsi = epsilons[i]
    epsilon = epsi*ones(K)
    tic = time()
    Ac,Ic,Jc,Wc = SymmetricCard_reduction(Edges,EdgesW,n,epsilon)
    preptime[i] = time()-tic
    N = size(Ac,1)
    I = [In;Ic]
    J = [Jn;Jc]
    W = [Wn;Wc]
    A = sparse(I,J,W,N,N)

    ## Run max flow
    sparse_edges = n + length(I)

    Svec = vec([svec;zeros(N-n)])
    Tvec = vec([tvec;zeros(N-n)])
    tic = time()
    F = maxflow(A,Svec,Tvec,0.0)
    flowtime = time()-tic
    S = source_nodes_min(F,1e-10)[2:end].-1
    ct = F.cutvalue
    cval = check_cutval(A,Svec,Tvec,S)
    @assert(abs(ct-cval)<1e-8)
    Strue = intersect(S,1:n)
    original_obj = eval_fS_image(Edges,EdgesW,svec,tvec,An,Strue)

    # timeratio = flowtime/opttime
    # approx = original_obj/opt
    # sparsity = sparse_edges/total_edges
    objvals[i] = original_obj
    approxcuts[i] = ct     # cut value in the approximate reduced graph
    numedges[i] = sparse_edges
    runtimes[i] = flowtime

    println("$epsi \t $original_obj \t $sparse_edges \t $flowtime")

end

## Save the output
matwrite("Output/SparseCard_$(dataset)_$(numsuperpix)_$small.mat",Dict("objvals"=>objvals,
"numedges" => numedges, "runtimes"=>runtimes,"preptime"=>preptime,
"epsilons"=>epsilons,"approxcuts"=>approxcuts))

end

end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
function inverseSide(side::HorizonSide)
    return HorizonSide(mod(Int(side)+2, 4))
end

function main(r::Robot)
    stepsToWest = moveSide(r, West)
    stepsToNord = moveSide(r, Nord)

    while !ismarker(r)
        moveSide(r, Sud, true, true)
        if isborder(r, Ost) != true
            moveStep(r, Ost, true)
        end
        moveSide(r, Nord, true, true)
        if isborder(r, Ost) != true
            moveStep(r, Ost, true)
        end
    end
    moveSide(r, West)
    moveSide(r, Nord)

    followByCount(r, Ost, stepsToWest)
    followByCount(r, Sud, stepsToNord)
end

function moveSide(r::Robot, side::HorizonSide, needMark=false, needCheck=false)
    steps = 0
    while isborder(r, side) != true
        moveStep(r, side, needMark)
        if needCheck
            if ismarker(r)
                return
            end
        end
        steps += 1
    end
    return steps
end

function moveStep(r::Robot, side::HorizonSide, needMark::Bool)
    if needMark
        putmarker!(r)
    end
    move!(r, side)
end

function followByCount(r::Robot, side::HorizonSide, count::Int)
    for i = 0:count-1
        moveStep(r, side, false)
    end
end

max_stars_repo_name: TheDarkMoon/mirea-progs

max_issues_repo_name: TheDarkMoon/mirea-progs

max_forks_repo_name: TheDarkMoon/mirea-progs


text:
using HomotopyContinuation, LinearAlgebra, Test, Parameters, Random
import Arblib
using StaticArrays
using ProjectiveVectors: PVector, affine_chart
using TreeViews: TreeViews
const HC = HomotopyContinuation
import SemialgebraicSets

include("test_systems.jl")

function test_treeviews(x)
    @test TreeViews.hastreeview(x)
    @test_nowarn TreeViews.treelabel(devnull, x, MIME"application/prs.juno.inline"())
    for i = 1:TreeViews.numberofnodes(x)
        @test_nowarn TreeViews.nodelabel(devnull, x, i, MIME"application/prs.juno.inline"())
        @test_nowarn TreeViews.treenode(x, i)
    end
end

function test_show_juno(x)
    @test show(stdout, MIME("application/prs.juno.inline"), x) === x
end

@testset "HomotopyContinuation" begin
    include("utils_test.jl")
    include("double_double_test.jl")
    include("binomial_system_test.jl")
    include("norm_test.jl")
    include("voronoi_tree_test.jl")
    include("unique_points_test.jl")
    include("model_kit_test.jl")
    include("linear_algebra_test.jl")
    include("homotopies_test.jl")
    include("systems_test.jl")
    include("tracker_test.jl")
    include("valuation_test.jl")
    include("endgame_tracker_test.jl")
    include("linear_test.jl")
    include("polyhedral_test.jl")
    include("solve_test.jl")
    include("endgame_test.jl")
    include("monodromy_test.jl")
    include("witness_set_test.jl")
    include("certification_test.jl")
    include("semialgebraic_sets_test.jl")

    if "extensive" in ARGS
        include("extensive/extensive_test.jl")
    end
end


text:
import Base64
import CompatHelper
import Dates
import GitHub
import HTTP
import JSON
import Pkg
import Printf
import Test
import TimeZones

const registries_1 = Pkg.Types.RegistrySpec[Pkg.RegistrySpec(name = "General",
                                                             uuid = "23338594-aafe-5451-b93e-139f81909106",
                                                             url = "https://github.com/JuliaRegistries/General.git")]

const registries_2 = Pkg.Types.RegistrySpec[Pkg.RegistrySpec(name = "General",
                                                             uuid = "23338594-aafe-5451-b93e-139f81909106",
                                                             url = "https://github.com/JuliaRegistries/General.git"),
                                            Pkg.RegistrySpec(name = "BioJuliaRegistry",
                                                             uuid = "ccbd2cc2-2954-11e9-1ccf-f3e7900901ca",
                                                             url = "https://github.com/BioJulia/BioJuliaRegistry.git")]

include("testutils.jl")

Test.@testset "CompatHelper.jl" begin
    Test.@testset "CompatHelper.jl unit tests" begin
        @info("Running the CompatHelper.jl unit tests")
        include("unit-tests.jl")
    end

    COMPATHELPER_RUN_INTEGRATION_TESTS = get(ENV, "COMPATHELPER_RUN_INTEGRATION_TESTS", "")::String
    if COMPATHELPER_RUN_INTEGRATION_TESTS == "true"
        Test.@testset "CompatHelper.jl integration tests" begin
            @info("Running the CompatHelper.jl integration tests")
            include("integration-tests.jl")
        end
    end
end


text:
module Atomistic

using AtomsBase
using InteratomicPotentials
using DFTK
using NBodySimulator
using PeriodicTable
using Plots
using PyCall
using Random
using StaticArrays
using Unitful
using UnitfulAtomic
using UnitfulRecipes

import Base: @kwdef, RefValue
import InteratomicPotentials: ArbitraryPotential
import DFTK: Mixing
import NBodySimulator: Body, NullThermostat, SimulationResult, Thermostat
import Plots: Plot

include("exceptions.jl")

export MolecularDynamicsResult, get_system, get_time_range, temperature, reference_temperature, kinetic_energy, potential_energy, total_energy, rdf
include("interfaces/md_result.jl")
export MolecularDynamicsSimulator, simulate
include("interfaces/md_simulator.jl")

export DynamicAtom, DynamicSystem
include("integrations/ab_integration.jl")
export LJPotential, NBSimulator, NBSResult
include("integrations/nbs_integration.jl")
export DFTKPotential, analyze_convergence
include("integrations/dftk_integration.jl")

export plot_temperature, plot_temperature!, plot_energy, plot_energy!, plot_rdf
include("analysis/plotting.jl")
export write_nbs_animation, write_ase_trajectory
include("analysis/visualization.jl")

end

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia