text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

__precompile__(true)

"""
Profiling support, main entry point is the [`@profile`](@ref) macro.
"""
module Profile

import Base.StackTraces: lookup, UNKNOWN, show_spec_linfo
using Base: iszero

export @profile

"""
    @profile

`@profile <expression>` runs your expression while taking periodic backtraces. These are
appended to an internal buffer of backtraces.
"""
macro profile(ex)
    quote
        try
            status = start_timer()
            if status < 0
                error(error_codes[status])
            end
            $(esc(ex))
        finally
            stop_timer()
        end
    end
end

####
#### User-level functions
####

"""
    init(; n::Integer, delay::Float64)

Configure the `delay` between backtraces (measured in seconds), and the number `n` of
instruction pointers that may be stored. Each instruction pointer corresponds to a single
line of code; backtraces generally consist of a long list of instruction pointers. Default
settings can be obtained by calling this function with no arguments, and each can be set
independently using keywords or in the order `(n, delay)`.
"""
function init(; n::Union{Void,Integer} = nothing, delay::Union{Void,Float64} = nothing)
    n_cur = ccall(:jl_profile_maxlen_data, Csize_t, ())
    delay_cur = ccall(:jl_profile_delay_nsec, UInt64, ())/10^9
    if n === nothing && delay === nothing
        return Int(n_cur), delay_cur
    end
    nnew = (n === nothing) ? n_cur : n
    delaynew = (delay === nothing) ? delay_cur : delay
    init(nnew, delaynew)
end

function init(n::Integer, delay::Float64)
    status = ccall(:jl_profile_init, Cint, (Csize_t, UInt64), n, round(UInt64,10^9*delay))
    if status == -1
        error("could not allocate space for ", n, " instruction pointers")
    end
end

# init with default values
# Use a max size of 1M profile samples, and fire timer every 1ms
if Sys.iswindows()
    __init__() = init(1_000_000, 0.01)
else
    __init__() = init(1_000_000, 0.001)
end

"""
    clear()

Clear any existing backtraces from the internal buffer.
"""
clear() = ccall(:jl_profile_clear_data, Void, ())

const LineInfoDict = Dict{UInt64, Vector{StackFrame}}
const LineInfoFlatDict = Dict{UInt64, StackFrame}

struct ProfileFormat
    maxdepth::Int
    mincount::Int
    noisefloor::Float64
    sortedby::Symbol
    combine::Bool
    C::Bool
    function ProfileFormat(;
        C = false,
        combine = true,
        maxdepth::Int = typemax(Int),
        mincount::Int = 0,
        noisefloor = 0,
        sortedby::Symbol = :filefuncline)
        return new(maxdepth, mincount, noisefloor, sortedby, combine, C)
    end
end

"""
    print([io::IO = STDOUT,] [data::Vector]; kwargs...)

Prints profiling results to `io` (by default, `STDOUT`). If you do not
supply a `data` vector, the internal buffer of accumulated backtraces
will be used.

The keyword arguments can be any combination of:

 - `format` -- Determines whether backtraces are printed with (default, `:tree`) or without (`:flat`)
   indentation indicating tree structure.

 - `C` -- If `true`, backtraces from C and Fortran code are shown (normally they are excluded).

 - `combine` -- If `true` (default), instruction pointers are merged that correspond to the same line of code.

 - `maxdepth` -- Limits the depth higher than `maxdepth` in the `:tree` format.

 - `sortedby` -- Controls the order in `:flat` format. `:filefuncline` (default) sorts by the source
    line, whereas `:count` sorts in order of number of collected samples.

 - `noisefloor` -- Limits frames that exceed the heuristic noise floor of the sample (only applies to format `:tree`).
    A suggested value to try for this is 2.0 (the default is 0). This parameter hides samples for which `n <= noisefloor * √N`,
    where `n` is the number of samples on this line, and `N` is the number of samples for the callee.

 - `mincount` -- Limits the printout to only those lines with at least `mincount` occurrences.
"""
function print(io::IO, data::Vector{<:Unsigned} = fetch(), lidict::LineInfoDict = getdict(data);
        format = :tree,
        C = false,
        combine = true,
        maxdepth::Int = typemax(Int),
        mincount::Int = 0,
        noisefloor = 0,
        sortedby::Symbol = :filefuncline)
    print(io, data, lidict, ProfileFormat(C = C,
            combine = combine,
            maxdepth = maxdepth,
            mincount = mincount,
            noisefloor = noisefloor,
            sortedby = sortedby),
        format)
end

function print(io::IO, data::Vector{<:Unsigned}, lidict::LineInfoDict, fmt::ProfileFormat, format::Symbol)
    cols::Int = Base.displaysize(io)[2]
    if format == :tree
        tree(io, data, lidict, cols, fmt)
    elseif format == :flat
        flat(io, data, lidict, cols, fmt)
    else
        throw(ArgumentError("output format $(repr(format)) not recognized"))
    end
end

"""
    print([io::IO = STDOUT,] data::Vector, lidict::LineInfoDict; kwargs...)

Prints profiling results to `io`. This variant is used to examine results exported by a
previous call to [`retrieve`](@ref). Supply the vector `data` of backtraces and
a dictionary `lidict` of line information.

See `Profile.print([io], data)` for an explanation of the valid keyword arguments.
"""
print(data::Vector{<:Unsigned} = fetch(), lidict::LineInfoDict = getdict(data); kwargs...) =
    print(STDOUT, data, lidict; kwargs...)

"""
    retrieve() -> data, lidict

"Exports" profiling results in a portable format, returning the set of all backtraces
(`data`) and a dictionary that maps the (session-specific) instruction pointers in `data` to
`LineInfo` values that store the file name, function name, and line number. This function
allows you to save profiling results for future analysis.
"""
function retrieve()
    data = fetch()
    return (copy(data), getdict(data))
end

function getdict(data::Vector{UInt})
    uip = unique(data)
    return LineInfoDict(UInt64(ip)=>lookup(ip) for ip in uip)
end

"""
    flatten(btdata, lidict) -> (newdata::Vector{UInt64}, newdict::LineInfoFlatDict)

Produces "flattened" backtrace data. Individual instruction pointers
sometimes correspond to a multi-frame backtrace due to inlining; in
such cases, this function inserts fake instruction pointers for the
inlined calls, and returns a dictionary that is a 1-to-1 mapping
between instruction pointers and a single StackFrame.
"""
function flatten(data::Vector, lidict::LineInfoDict)
    # Makes fake instruction pointers, counting down from typemax(UInt)
    newip = typemax(UInt64) - 1
    taken = Set(keys(lidict))  # make sure we don't pick one that's already used
    newdict = Dict{UInt64,StackFrame}()
    newmap  = Dict{UInt64,Vector{UInt64}}()
    for (ip, trace) in lidict
        if length(trace) == 1
            newdict[ip] = trace[1]
        else
            newm = UInt64[]
            for sf in trace
                while newip ∈ taken && newip > 0
                    newip -= 1
                end
                newip == 0 && error("all possible instruction pointers used")
                push!(newm, newip)
                newdict[newip] = sf
                newip -= 1
            end
            newmap[ip] = newm
        end
    end
    newdata = UInt64[]
    for ip in data
        local ip::UInt64
        if haskey(newmap, ip)
            append!(newdata, newmap[ip])
        else
            push!(newdata, ip)
        end
    end
    return (newdata, newdict)
end

"""
    callers(funcname, [data, lidict], [filename=<filename>], [linerange=<start:stop>]) -> Vector{Tuple{count, lineinfo}}

Given a previous profiling run, determine who called a particular function. Supplying the
filename (and optionally, range of line numbers over which the function is defined) allows
you to disambiguate an overloaded method. The returned value is a vector containing a count
of the number of calls and line information about the caller. One can optionally supply
backtrace `data` obtained from [`retrieve`](@ref); otherwise, the current internal
profile buffer is used.
"""
function callers end

function callers(funcname::String, bt::Vector, lidict::LineInfoDict; filename = nothing, linerange = nothing)
    if filename === nothing && linerange === nothing
        return callersf(li -> li.func == funcname, bt, lidict)
    end
    filename === nothing && throw(ArgumentError("if supplying linerange, you must also supply the filename"))
    if linerange === nothing
        return callersf(li -> li.func == funcname && li.file == filename, bt, lidict)
    else
        return callersf(li -> li.func == funcname && li.file == filename &&
            in(li.line, linerange), bt, lidict)
    end
end

callers(funcname::String; kwargs...) = callers(funcname, retrieve()...; kwargs...)
callers(func::Function, bt::Vector, lidict::LineInfoDict; kwargs...) =
    callers(string(func), bt, lidict; kwargs...)
callers(func::Function; kwargs...) = callers(string(func), retrieve()...; kwargs...)

##
## For --track-allocation
##
# Reset the malloc log. Used to avoid counting memory allocated during
# compilation.

"""
    clear_malloc_data()

Clears any stored memory allocation data when running julia with `--track-allocation`.
Execute the command(s) you want to test (to force JIT-compilation), then call
[`clear_malloc_data`](@ref). Then execute your command(s) again, quit
Julia, and examine the resulting `*.mem` files.
"""
clear_malloc_data() = ccall(:jl_clear_malloc_data, Void, ())

# C wrappers
start_timer() = ccall(:jl_profile_start_timer, Cint, ())

stop_timer() = ccall(:jl_profile_stop_timer, Void, ())

is_running() = ccall(:jl_profile_is_running, Cint, ())!=0

get_data_pointer() = convert(Ptr{UInt}, ccall(:jl_profile_get_data, Ptr{UInt8}, ()))

len_data() = convert(Int, ccall(:jl_profile_len_data, Csize_t, ()))

maxlen_data() = convert(Int, ccall(:jl_profile_maxlen_data, Csize_t, ()))

error_codes = Dict(
    -1=>"cannot specify signal action for profiling",
    -2=>"cannot create the timer for profiling",
    -3=>"cannot start the timer for profiling",
    -4=>"cannot unblock SIGUSR1")

"""
    fetch() -> data

Returns a reference to the internal buffer of backtraces. Note that subsequent operations,
like [`clear`](@ref), can affect `data` unless you first make a copy. Note that the
values in `data` have meaning only on this machine in the current session, because it
depends on the exact memory addresses used in JIT-compiling. This function is primarily for
internal use; [`retrieve`](@ref) may be a better choice for most users.
"""
function fetch()
    len = len_data()
    maxlen = maxlen_data()
    if (len == maxlen)
        warn("""The profile data buffer is full; profiling probably terminated
                before your program finished. To profile for longer runs, call Profile.init
                with a larger buffer and/or larger delay.""")
    end
    return unsafe_wrap(Array, get_data_pointer(), (len,))
end

# Number of backtrace "steps" that are triggered by taking the backtrace, e.g., inside profile_bt
# TODO: may be platform-specific?
const btskip = 0

## Print as a flat list
# Counts the number of times each line appears, at any nesting level
function count_flat(data::Vector{T}) where T<:Unsigned
    linecount = Dict{T,Int}()
    toskip = btskip
    for ip in data
        if toskip > 0
            toskip -= 1
            continue
        end
        if ip == 0
            toskip = btskip
            continue
        end
        linecount[ip] = get(linecount, ip, 0)+1
    end
    iplist = Vector{T}()
    n = Vector{Int}()
    for (k,v) in linecount
        push!(iplist, k)
        push!(n, v)
    end
    return (iplist, n)
end

function parse_flat(iplist, n, lidict::LineInfoFlatDict, C::Bool)
    # Convert instruction pointers to names & line numbers
    lilist = [lidict[ip] for ip in iplist]
    # Keep only the interpretable ones
    # The ones with no line number might appear multiple times in a single
    # backtrace, giving the wrong impression about the total number of backtraces.
    # Delete them too.
    keep = .!Bool[x == UNKNOWN || x.line == 0 || (x.from_c && !C) for x in lilist]
    n = n[keep]
    lilist = lilist[keep]
    return (lilist, n)
end

function flat(io::IO, data::Vector, lidict::LineInfoFlatDict, cols::Int, fmt::ProfileFormat)
    if !fmt.C
        data = purgeC(data, lidict)
    end
    iplist, n = count_flat(data)
    if isempty(n)
        warning_empty()
        return
    end
    lilist, n = parse_flat(iplist, n, lidict, fmt.C)
    print_flat(io, lilist, n, cols, fmt)
    nothing
end

function flat(io::IO, data::Vector, lidict::LineInfoDict, cols::Int, fmt::ProfileFormat)
    newdata, newdict = flatten(data, lidict)
    flat(io, newdata, newdict, cols, fmt)
    nothing
end

function print_flat(io::IO, lilist::Vector{StackFrame}, n::Vector{Int},
        cols::Int, fmt::ProfileFormat)
    p = liperm(lilist)
    lilist = lilist[p]
    n = n[p]
    if fmt.combine
        j = 1
        for i = 2:length(lilist)
            if lilist[i] == lilist[j]
                n[j] += n[i]
                n[i] = 0
            else
                j = i
            end
        end
        keep = n .> 0
        n = n[keep]
        lilist = lilist[keep]
    end
    if fmt.sortedby == :count
        p = sortperm(n)
        n = n[p]
        lilist = lilist[p]
    end
    wcounts = max(6, ndigits(maximum(n)))
    maxline = 0
    maxfile = 6
    maxfunc = 10
    for li in lilist
        maxline = max(maxline, li.line)
        maxfile = max(maxfile, length(string(li.file)))
        maxfunc = max(maxfunc, length(string(li.func)))
    end
    wline = max(5, ndigits(maxline))
    ntext = cols - wcounts - wline - 3
    maxfunc += 25
    if maxfile + maxfunc <= ntext
        wfile = maxfile
        wfunc = maxfunc
    else
        wfile = floor(Integer, 2*ntext/5)
        wfunc = floor(Integer, 3*ntext/5)
    end
    println(io, lpad("Count", wcounts, " "), " ", rpad("File", wfile, " "), " ",
        lpad("Line", wline, " "), " ", rpad("Function", wfunc, " "))
    for i = 1:length(n)
        n[i] < fmt.mincount && continue
        li = lilist[i]
        Base.print(io, lpad(string(n[i]), wcounts, " "), " ")
        Base.print(io, rpad(rtruncto(string(li.file), wfile), wfile, " "), " ")
        Base.print(io, lpad(string(li.line), wline, " "), " ")
        fname = string(li.func)
        if !li.from_c && !isnull(li.linfo)
            fname = sprint(show_spec_linfo, li)
        end
        Base.print(io, rpad(ltruncto(fname, wfunc), wfunc, " "))
        println(io)
    end
    nothing
end

## A tree representation
# Identify and counts repetitions of all unique backtraces
function tree_aggregate(data::Vector{UInt64})
    iz = find(iszero, data)  # find the breaks between backtraces
    treecount = Dict{Vector{UInt64},Int}()
    istart = 1 + btskip
    for iend in iz
        tmp = data[iend - 1 : -1 : istart]
        treecount[tmp] = get(treecount, tmp, 0) + 1
        istart = iend + 1 + btskip
    end
    bt = Vector{Vector{UInt64}}()
    counts = Vector{Int}()
    for (k, v) in treecount
        if !isempty(k)
            push!(bt, k)
            push!(counts, v)
        end
    end
    return (bt, counts)
end

tree_format_linewidth(x::StackFrame) = ndigits(x.line) + 6

function tree_format(lilist::Vector{StackFrame}, counts::Vector{Int}, level::Int, cols::Int)
    nindent = min(cols>>1, level)
    ndigcounts = ndigits(maximum(counts))
    ndigline = maximum([tree_format_linewidth(x) for x in lilist])
    ntext = cols - nindent - ndigcounts - ndigline - 5
    widthfile = floor(Integer, 0.4ntext)
    widthfunc = floor(Integer, 0.6ntext)
    strs = Vector{String}(uninitialized, length(lilist))
    showextra = false
    if level > nindent
        nextra = level - nindent
        nindent -= ndigits(nextra) + 2
        showextra = true
    end
    for i = 1:length(lilist)
        li = lilist[i]
        if li != UNKNOWN
            base = " "^nindent
            if showextra
                base = string(base, "+", nextra, " ")
            end
            if li.line == li.pointer
                strs[i] = string(base,
                    rpad(string(counts[i]), ndigcounts, " "),
                    " ",
                    "unknown function (pointer: 0x",
                    hex(li.pointer,2*sizeof(Ptr{Void})),
                    ")")
            else
                fname = string(li.func)
                if !li.from_c && !isnull(li.linfo)
                    fname = sprint(show_spec_linfo, li)
                end
                strs[i] = string(base,
                    rpad(string(counts[i]), ndigcounts, " "),
                    " ",
                    rtruncto(string(li.file), widthfile),
                    ":",
                    li.line == -1 ? "?" : string(li.line),
                    "; ",
                    ltruncto(fname, widthfunc))
            end
        else
            strs[i] = ""
        end
    end
    return strs
end

# Print a "branch" starting at a particular level. This gets called recursively.
function tree(io::IO, bt::Vector{Vector{UInt64}}, counts::Vector{Int},
        lidict::LineInfoFlatDict, level::Int, cols::Int, fmt::ProfileFormat, noisefloor::Int)
    if level > fmt.maxdepth
        return
    end
    # Organize backtraces into groups that are identical up to this level
    if fmt.combine
        # Combine based on the line information
        d = Dict{StackFrame,Vector{Int}}()
        for i = 1:length(bt)
            ip = bt[i][level + 1]
            key = lidict[ip]
            indx = Base.ht_keyindex(d, key)
            if haskey(d, key)
                push!(d[key], i)
            else
                d[key] = [i]
            end
        end
        # Generate counts
        dlen = length(d)
        lilist = Vector{StackFrame}(uninitialized, dlen)
        group = Vector{Vector{Int}}(uninitialized, dlen)
        n = Vector{Int}(uninitialized, dlen)
        i = 1
        for (key, v) in d
            lilist[i] = key
            group[i] = v
            n[i] = sum(counts[v])
            i += 1
        end
    else
        # Combine based on the instruction pointer
        d = Dict{UInt64,Vector{Int}}()
        for i = 1:length(bt)
            key = bt[i][level+1]
            if haskey(d, key)
                push!(d[key], i)
            else
                d[key] = [i]
            end
        end
        # Generate counts, and do the code lookup
        dlen = length(d)
        lilist = Vector{StackFrame}(uninitialized, dlen)
        group = Vector{Vector{Int}}(uninitialized, dlen)
        n = Vector{Int}(uninitialized, dlen)
        i = 1
        for (key, v) in d
            lilist[i] = lidict[key]
            group[i] = v
            n[i] = sum(counts[v])
            i += 1
        end
    end
    # Order the line information
    if length(lilist) > 1
        p = liperm(lilist)
        lilist = lilist[p]
        group = group[p]
        n = n[p]
    end
    # Generate the string for each line
    strs = tree_format(lilist, n, level, cols)
    # Recurse to the next level
    len = Int[length(x) for x in bt]
    for i = 1:length(lilist)
        n[i] < fmt.mincount && continue
        n[i] < noisefloor && continue
        if !isempty(strs[i])
            println(io, strs[i])
        end
        idx = group[i]
        keep = len[idx] .> level+1
        if any(keep)
            idx = idx[keep]
            tree(io, bt[idx], counts[idx], lidict, level + 1, cols, fmt, fmt.noisefloor > 0 ? floor(Int, fmt.noisefloor * sqrt(n[i])) : 0)
        end
    end
    nothing
end

function tree(io::IO, data::Vector{UInt64}, lidict::LineInfoFlatDict, cols::Int, fmt::ProfileFormat)
    if !fmt.C
        data = purgeC(data, lidict)
    end
    bt, counts = tree_aggregate(data)
    if isempty(counts)
        warning_empty()
        return
    end
    level = 0
    len = Int[length(x) for x in bt]
    keep = len .> 0
    tree(io, bt[keep], counts[keep], lidict, level, cols, fmt, 0)
    nothing
end

function tree(io::IO, data::Vector, lidict::LineInfoDict, cols::Int, fmt::ProfileFormat)
    newdata, newdict = flatten(data, lidict)
    tree(io, newdata, newdict, cols, fmt)
    nothing
end

function callersf(matchfunc::Function, bt::Vector, lidict::LineInfoDict)
    counts = Dict{StackFrame, Int}()
    lastmatched = false
    for id in bt
        if id == 0
            lastmatched = false
            continue
        end
        li = lidict[id]
        if lastmatched
            if haskey(counts, li)
                counts[li] += 1
            else
                counts[li] = 1
            end
        end
        lastmatched = matchfunc(li)
    end
    k = collect(keys(counts))
    v = collect(values(counts))
    p = sortperm(v, rev=true)
    return [(v[i], k[i]) for i in p]
end

# Utilities
function rtruncto(str::String, w::Int)
    if length(str) <= w
        return str
    else
        return string("...", str[chr2ind(str, length(str)-w+4):end])
    end
end
function ltruncto(str::String, w::Int)
    if length(str) <= w
        return str
    else
        return string(str[1:chr2ind(str,w-4)], "...")
    end
end

truncto(str::Symbol, w::Int) = truncto(string(str), w)

# Order alphabetically (file, function) and then by line number
function liperm(lilist::Vector{StackFrame})
    comb = Vector{String}(uninitialized, length(lilist))
    for i = 1:length(lilist)
        li = lilist[i]
        if li != UNKNOWN
            comb[i] = @sprintf("%s:%s:%06d", li.file, li.func, li.line)
        else
            comb[i] = "zzz"
        end
    end
    return sortperm(comb)
end

warning_empty() = warn("""
            There were no samples collected. Run your program longer (perhaps by
            running it multiple times), or adjust the delay between samples with
            Profile.init().""")

function purgeC(data::Vector{UInt64}, lidict::LineInfoFlatDict)
    keep = Bool[d == 0 || lidict[d].from_c == false for d in data]
    return data[keep]
end

end # module

max_stars_repo_name: ihnorton/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using NumericalRange
using Test

@testset "NumericalRange.jl" begin
    # Write your tests here.
end

ext: jl

max_stars_repo_name: ThomasChaffey/NumericalRange


text:
using Distributions
using Base.Test
using SkewDist

if isa(Pkg.installed("Calculus"), VersionNumber)
    using Calculus
    function mean_mgf_consistent(dist::ContinuousUnivariateDistribution)
        d = derivative(x->mgf(dist, x))
        m0 = d(0)
        @test_approx_eq derivative(x->mgf(dist, x), 0.0) mean(dist)
    end
    function std_mgf_consistent(dist::ContinuousUnivariateDistribution)
        d1 = derivative(x->mgf(dist,x))
        d2 = derivative(d1)
        @test_approx_eq_eps sqrt(d2(0) - d1(0)^2)  std(dist) 1e-6
    end
else
    mean_mgf_consistent(dist::ContinuousUnivariateDistribution) = true
    std_mgf_consistent(dist::ContinuousUnivariateDistribution) = true
end

function pdf_integrates_to_one(dist::ContinuousUnivariateDistribution)
    @test_approx_eq quadgk(t->pdf(dist, t), minimum(dist), maximum(dist))[1] 1.0
end

function mean_pdf_consistent(dist::ContinuousUnivariateDistribution)
    @test_approx_eq quadgk(t->t*pdf(dist, t), minimum(dist), maximum(dist))[1] mean(dist)
end

function quantile_cdf_consistent(dist::ContinuousUnivariateDistribution)
    for β in 0.05:0.05:0.95
        @test_approx_eq β cdf(dist, quantile(dist, β))
    end
end

# function mean_rand_consistent(dist::ContinuousUnivariateDistribution, n::Int, β::Real)
#     samp = rand(dist, n)
#     m0 = mean(samp)

# end

function test_dist(dist::ContinuousUnivariateDistribution)
    println("Testing distribution $(dist)...")
    pdf_integrates_to_one(dist)
    mean_pdf_consistent(dist)
    quantile_cdf_consistent(dist)
    try
        mean_mgf_consistent(dist)
        std_mgf_consistent(dist)
    catch err
        if !isa(err, MethodError)
            rethrow(err)
        end
    end
end

skewnorm1 = SkewNormal(1.0, 1.0, -1.0)
test_dist(skewnorm1)

skewtdist1 = SkewTDist(1.0, 2.0, 2.0, 3.0)
test_dist(skewtdist1)


text:
function part1(xs)
    sum(diff(xs) .> 0)
end

function valid_conv(xs, w)
    n = length(xs)
    m = length(w)
    rw = reverse(w)

    [sum(xs[i:i+m-1] .* rw) for i=1:(n-m+1)]
end

function part2(xs)
    w = repeat(1:1, 3)
    part1(valid_conv(xs, w))
end

xs = open(ARGS[1]) do fp
    map(x -> parse(Int, x), split(read(fp, String)))
end

println("Part 1: ", part1(xs))
println("Part 2: ", part2(xs))


text:
using ModelingToolkit, OrdinaryDiffEq, Test, LinearAlgebra

@parameters t σ ρ β
@variables x(t) y(t) z(t) k(t)
D = Differential(t)

eqs = [D(D(x)) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z]

@named sys′ = ODESystem(eqs)
sys = ode_order_lowering(sys′)

eqs2 = [0 ~ x * y - k,
    D(D(x)) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z]
@named sys2 = ODESystem(eqs2, t, [x, y, z, k], parameters(sys′))
sys2 = ode_order_lowering(sys2)
# test equation/varible ordering
ModelingToolkit.calculate_massmatrix(sys2) == Diagonal([1, 1, 1, 1, 0])

u0 = [D(x) => 2.0,
    x => 1.0,
    y => 0.0,
    z => 0.0]

p = [σ => 28.0,
    ρ => 10.0,
    β => 8 / 3]

tspan = (0.0, 100.0)
prob = ODEProblem(sys, u0, tspan, p, jac = true)
probexpr = ODEProblemExpr(sys, u0, tspan, p, jac = true)
sol = solve(prob, Tsit5())
solexpr = solve(eval(prob), Tsit5())
@test all(x -> x == 0, Array(sol - solexpr))
#using Plots; plot(sol,vars=(:x,:y))

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(x) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z]

lorenz1 = ODESystem(eqs, name = :lorenz1)
lorenz2 = ODESystem(eqs, name = :lorenz2)

@variables α(t)
@parameters γ
connections = [0 ~ lorenz1.x + lorenz2.y + α * γ]
@named connected = ODESystem(connections, t, [α], [γ], systems = [lorenz1, lorenz2])

u0 = [lorenz1.x => 1.0,
    lorenz1.y => 0.0,
    lorenz1.z => 0.0,
    lorenz2.x => 0.0,
    lorenz2.y => 1.0,
    lorenz2.z => 0.0,
    α => 2.0]

p = [lorenz1.σ => 10.0,
    lorenz1.ρ => 28.0,
    lorenz1.β => 8 / 3,
    lorenz2.σ => 10.0,
    lorenz2.ρ => 28.0,
    lorenz2.β => 8 / 3,
    γ => 2.0]

tspan = (0.0, 100.0)
prob = ODEProblem(connected, u0, tspan, p)
sol = solve(prob, Rodas5())
@test maximum(sol[2, :] + sol[6, :] + 2sol[1, :]) < 1e-12
#using Plots; plot(sol,vars=(:α,Symbol(lorenz1.x),Symbol(lorenz2.y)))


text:
"""
    _screen_args(argtuple, kwargtuple, modsname)

Internal helper to check the presence and shape of modification and argument
arrays in [`screen`](@ref) and pals.
"""
function _screen_args(argtuple, kwargtuple, modsname)

    mods = get(kwargtuple, modsname, nothing)
    args = get(kwargtuple, :args, nothing)

    if isnothing(mods)
        if isnothing(args)
            throw(
                DomainError(args, "at least one of `$modsname` and `args` must be defined"),
            )
        end
        return NamedTuple{(modsname,)}(Ref([[] for _ in args]))
    elseif isnothing(args)
        return (args = [() for _ in mods],)
    else
        if size(mods) != size(args)
            throw(
                DomainError(
                    "$(size(mods)) != $(size(args))",
                    "sizes of `$modsname` and `args` differ",
                ),
            )
        end
        return ()
    end
end

"""
    screen(
        model::MetabolicModel;
        variants::Array{V,N}, # defaults to an array of identities
        analysis,
        args::Array{A,N}, # defaults to an array of empty argument lists
        workers = [myid()],
    )::Array where {V<:AbstractVector,A,N}

Take an array of model-modifying function vectors in `variants`, and execute
the function `analysis` on all variants of the `model` specified by `variants`.
The computation is distributed over worker IDs in `workers`. If `args` are
supplied (as an array of the same size as the `variants`), they are forwarded
as arguments to the corresponding analysis function calls.

The array of variants must contain vectors of single-parameter functions, these
are applied to model in order. The functions must *not* modify the model, but
rather return a modified copy. The copy should be made as shallow as possible,
to increase memory efficiency of the process. Variant generators that modify
the argument model in-place will cause unpredictable results. Refer to the
definition of [`screen_variant`](@ref) for details.

The function `analysis` will receive a single argument (the modified model),
together with arguments from `args` expanded by `...`. Supply an array of
tuples or vectors to pass in multiple arguments to each function. If the
argument values should be left intact (not expanded to multiple arguments),
they must be wrapped in single-item tuple or vector.

The modification and analysis functions are transferred to `workers` as-is; all
packages required to run them (e.g. the optimization solvers) must be loaded
there. Typically, you want to use the macro `@everywhere using
MyFavoriteSolver` from `Distributed` package for loading the solvers.

# Return value

The results of running `analysis` are collected in to the resulting array, in a
way that preserves the shape of the `variants`, similarly as with `pmap`.

The results of `analysis` function must be serializable, preferably made only
from pure Julia structures, because they may be transferred over the network
between the computation nodes. For that reason, functions that return whole
JuMP models that contain pointers to allocated C structures (such as
[`flux_balance_analysis`](@ref) used with `GLPK` or `Gurobi` otimizers) will
generally not in this context.

Note: this function is a thin argument-handling wrapper around
[`_screen_impl`](@ref).

# Example
```
function reverse_reaction(i::Int)
    (model::CoreModel) -> begin
        mod = copy(model)
        mod.S[:,i] .*= -1   # this is unrealistic but sufficient for demonstration
        mod
    end
end

m = load_model(CoreModel, "e_coli_core.xml")

screen(m,
    variants = [
        [reverse_reaction(5)],
        [reverse_reaction(3), reverse_reaction(6)]
    ],
    analysis = mod -> mod.S[:,3])  # observe the changes in S

screen(m,
    variants = [
        [reverse_reaction(5)],
        [reverse_reaction(3), reverse_reaction(6)]
    ],
    analysis = mod -> flux_balance_analysis_vec(mod, GLPK.Optimizer))
```
"""
screen(args...; kwargs...) =
    _screen_impl(args...; kwargs..., _screen_args(args, kwargs, :variants)...)

"""
    _screen_impl(
        model::MetabolicModel;
        variants::Array{V,N},
        analysis,
        args::Array{A,N},
        workers = [myid()],
    )::Array where {V<:AbstractVector,A,N}

The actual implementation of [`screen`](@ref).
"""
function _screen_impl(
    model::MetabolicModel;
    variants::Array{V,N},
    analysis,
    args::Array{A,N},
    workers = [myid()],
)::Array where {V<:AbstractVector,A,N}

    map(fetch, save_at.(workers, :cobrexa_screen_variants_model, Ref(model)))
    map(fetch, save_at.(workers, :cobrexa_screen_variants_analysis_fn, Ref(analysis)))
    map(fetch, get_from.(workers, Ref(:(precache!(cobrexa_screen_variants_model)))))

    res = pmap(
        (vars, args)::Tuple -> screen_variant(
            (@remote cobrexa_screen_variants_model),
            vars,
            (@remote cobrexa_screen_variants_analysis_fn),
            args,
        ),
        CachingPool(workers),
        zip(variants, args),
    )

    map(fetch, remove_from.(workers, :cobrexa_screen_variants_model))
    map(fetch, remove_from.(workers, :cobrexa_screen_variants_analysis_fn))

    return res
end

"""
    screen_variant(model::MetabolicModel, variant::Vector, analysis, args = ())

Helper function for [`screen`](@ref) that applies all single-argument
functions in `variant` to the `model` (in order from "first" to
"last"), and executes `analysis` on the result.

Can be used to test model variants locally.
"""
function screen_variant(model::MetabolicModel, variant::Vector, analysis, args = ())
    for fn in variant
        model = fn(model)
    end
    analysis(model, args...)
end

"""
    screen_variants(model, variants, analysis; workers=[myid()])

A shortcut for [`screen`](@ref) that only works with model variants.
"""
screen_variants(model, variants, analysis; workers = [myid()]) =
    screen(model; variants = variants, analysis = analysis, workers = workers)

"""
    screen_optimize_objective(_, optmodel)::Maybe{Float64}

A variant of [`optimize_objective`](@ref) directly usable in
[`screen_optmodel_modifications`](@ref).
"""
screen_optimize_objective(_, optmodel)::Maybe{Float64} = optimize_objective(optmodel)

"""
    _screen_optmodel_prepare(model, optimizer, common_modifications)

Internal helper for [`screen_optmodel_modifications`](@ref) that creates the
model and applies the modifications.
"""
function _screen_optmodel_prepare(model, optimizer, common_modifications)
    precache!(model)
    optmodel = make_optimization_model(model, optimizer)
    for mod in common_modifications
        mod(model, optmodel)
    end
    return (model, optmodel)
end

"""
    _screen_optmodel_item((mods, args))

Internal helper for [`screen_optmodel_modifications`](@ref) that computes one
item of the screening task.
"""
function _screen_optmodel_item((mods, args))
    (model, optmodel) = @remote cobrexa_screen_optmodel_modifications_data
    for mod in mods
        mod(model, optmodel)
    end
    (@remote cobrexa_screen_optmodel_modifications_fn)(model, optmodel, args...)
end

"""
    screen_optmodel_modifications(
        model::MetabolicModel,
        optimizer;
        common_modifications::VF = [],
        modifications::Array{V,N}, # defaults to an array with no modifications
        args::Array{A,N}, # defaults to an array of empty argument lists
        analysis::Function,
        workers = [myid()],
    )::Array where {V<:AbstractVector,VF<:AbstractVector,A,N}

Screen multiple modifications of the same optimization model.

This function is potentially more efficient than [`screen`](@ref) because it
avoids making variants of the model structure and remaking of the optimization
model. On the other hand, modification functions need to keep the optimization
model in a recoverable state (one that leaves the model usable for the next
modification), which limits the possible spectrum of modifications applied.

Internally, `model` is distributed to `workers` and transformed into the
optimization model using [`make_optimization_model`](@ref).
`common_modifications` are applied to the models at that point. Next, vectors
of functions in `modifications` are consecutively applied, and the result of
`analysis` function called on model are collected to an array of the same
extent as `modifications`. Calls of `analysis` are optionally supplied with
extra arguments from `args` expanded with `...`, just like in [`screen`](@ref).

Both the modification functions (in vectors) and the analysis function here
have 2 base parameters (as opposed to 1 with [`screen`](@ref)): first is the
`model` (carried through as-is), second is the prepared JuMP optimization model
that may be modified and acted upon. As an example, you can use modification
[`change_constraint`](@ref) and analysis [`screen_optimize_objective`](@ref).

Note: This function is a thin argument-handling wrapper around
[`_screen_optmodel_modifications_impl`](@ref).
"""
screen_optmodel_modifications(args...; kwargs...) = _screen_optmodel_modifications_impl(
    args...;
    kwargs...,
    _screen_args(args, kwargs, :modifications)...,
)

"""
    _screen_optmodel_modifications_impl(
        model::MetabolicModel,
        optimizer;
        common_modifications::VF = [],
        modifications::Array{V,N},
        args::Array{A,N},
        analysis::Function,
        workers = [myid()],
    )::Array where {V<:AbstractVector,VF<:AbstractVector,A,N}

The actual implementation of [`screen_optmodel_modifications`](@ref).
"""
function _screen_optmodel_modifications_impl(
    model::MetabolicModel,
    optimizer;
    common_modifications::VF = [],
    modifications::Array{V,N},
    args::Array{A,N},
    analysis::Function,
    workers = [myid()],
)::Array where {V<:AbstractVector,VF<:AbstractVector,A,N}

    map(
        fetch,
        save_at.(
            workers,
            :cobrexa_screen_optmodel_modifications_data,
            Ref(
                :($COBREXA._screen_optmodel_prepare(
                    $model,
                    $optimizer,
                    $common_modifications,
                )),
            ),
        ),
    )
    map(fetch, save_at.(workers, :cobrexa_screen_optmodel_modifications_fn, Ref(analysis)))

    res = pmap(_screen_optmodel_item, CachingPool(workers), zip(modifications, args))

    map(fetch, remove_from.(workers, :cobrexa_screen_optmodel_modifications_data))
    map(fetch, remove_from.(workers, :cobrexa_screen_optmodel_modifications_fn))

    return res
end


text:
using Logging, LinearAlgebra, DataFrames, StatsModels, Random

"""
$(TYPEDSIGNATURES)
"""
function d_criterion(model_matrix::Array{Float64, 2})
    det((model_matrix' * model_matrix) / size(model_matrix, 1)) ^ (1 / size(model_matrix, 2))
end

"""
$(TYPEDSIGNATURES)
"""
function linear_model(data::DataFrame)
    Term(:_) ~ sum(term.(names(data)))
end

"""
$(TYPEDSIGNATURES)

Computes `deltas` for  the exchanges between all pairs  made between experiments
in `model_matrix`  and experiments in `candidate_set`.   The `dispersion_matrix`
array   must   contain  the   dispersion   matrix   for  `model_matrix`.    Both
`model_matrix`  and   `candidate_set`  must  be  model   matrices.   The  arrays
`design_variance`,  `candidate_variance`,  `predictions`  and `deltas`  will  be
overwritten.
"""
function variance_prediction!(dispersion_matrix::Array{Float64, 2},
                              model_matrix::Array{Float64, 2},
                              candidate_set::Array{Float64, 2},
                              design_variance::Array{Float64, 1},
                              candidate_variance::Array{Float64, 1},
                              predictions::Array{Float64, 2},
                              deltas::Array{Float64, 2})
    @inbounds Threads.@threads for i = 1:size(model_matrix, 1)
        design_variance[i] = model_matrix[i, :]' * dispersion_matrix * model_matrix[i, :]
    end

    @inbounds Threads.@threads for i = 1:size(candidate_set, 1)
        candidate_variance[i] = candidate_set[i, :]' * dispersion_matrix * candidate_set[i, :]
    end

    @inbounds Threads.@threads for i = 1:size(predictions, 1)
        @inbounds for j = 1:size(predictions, 2)
            predictions[i, j] = model_matrix[i, :]' * dispersion_matrix * candidate_set[j, :]
        end
    end

    @inbounds Threads.@threads for i = 1:size(deltas, 1)
        @inbounds for j = 1:size(deltas, 2)
            deltas[i, j] = candidate_variance[j] - design_variance[i] -
                           ((design_variance[i] * candidate_variance[j]) -
                            (predictions[i, j] ^ 2))
        end
    end

    return
end

"""
$(TYPEDSIGNATURES)

Generates a random design of size `size` with `factors` numerical factors.
"""
function generate_random_design(factors::Int, size::Int)
    data = Dict("_" => ones(size))

    for i in 1:factors
        data["X$i"] = rand(size)
    end

    DataFrame(data)
end

"""
$(TYPEDSIGNATURES)

Generates  a design  from all  combinations  of `factors`  factors and  `levels`
levels.
"""
function expanded_design(factors::Int, levels::Int)
    expanded_grid = expand_grid(factors, levels)
    data = Dict("_" => ones(levels ^ factors))

    for i in 1:factors
        data["X$i"] = expanded_grid[:, i]
    end

    DataFrame(data)
end

"""
$(TYPEDSIGNATURES)

Generates all combinations of `factors` factors with `levels` levels.
"""
function expand_grid(factors::Int, levels::Int)
    chosen_levels = (((collect(1:levels) .- 1) / (levels - 1)) .* (2)) .- 1
    complete_set  = Array{Float64, 2}(undef, levels ^ factors, factors)
    iterator      = Base.product(repeat([chosen_levels], inner = factors)...)
    set_index     = 1

    for i in iterator
        complete_set[set_index, :] = collect(i)
        set_index += 1
    end

    complete_set
end

"""
$(TYPEDSIGNATURES)
"""
function optimize_design(;factors::Int,
                         levels::Int,
                         experiments::Int,
                         design::DataFrame,
                         candidate_set::DataFrame,
                         iterations::Int,
                         refresh_candidate_set::Bool = true)
    candidate_set_size = size(candidate_set, 1)

    formula = linear_model(design)

    new_schema = schema(formula, design)
    schema_formula = apply_schema(formula, new_schema)

    dummy_response, model_matrix = modelcols(schema_formula, design)
    dummy_reponse, candidate_model_matrix = modelcols(schema_formula, candidate_set)

    design_variance = similar(model_matrix[:, 1])
    candidate_variance = similar(candidate_model_matrix[:, 1])

    predictions = similar(candidate_model_matrix[:, 1],
                          (size(model_matrix, 1),
                           size(candidate_model_matrix, 1)))

    deltas = similar(candidate_model_matrix[:, 1],
                     (size(model_matrix, 1),
                      size(candidate_model_matrix, 1)))

    dispersion_matrix = inv(model_matrix' * model_matrix)

    variance_prediction!(dispersion_matrix,
                         model_matrix,
                         candidate_model_matrix,
                         design_variance,
                         candidate_variance,
                         predictions,
                         deltas)

    current_d_criterion = d_criterion(model_matrix)
    @debug "Starting D-Criterion" current_d_criterion

    for i in 1:iterations
        best_swap = findmax(deltas)
        @debug "Best Swap" best_swap

        if best_swap[1] <= 0.0
            @debug "Max. delta was less than, or equal to, zero"

            if refresh_candidate_set
                @debug "Regenerating candidate set..." i

                candidate_set = generate_random_design(factors, candidate_set_size)

                formula = linear_model(candidate_set)

                new_schema = schema(formula, candidate_set)
                schema_formula = apply_schema(formula, new_schema)

                dummy_reponse, candidate_model_matrix = modelcols(schema_formula, candidate_set)
            else
                @debug "Stopping..."
                break
            end
        else
            model_matrix[best_swap[2][1], :] = candidate_model_matrix[best_swap[2][2], :]
        end

        dispersion_matrix = inv(model_matrix' * model_matrix)

        variance_prediction!(dispersion_matrix,
                             model_matrix,
                             candidate_model_matrix,
                             design_variance,
                             candidate_variance,
                             predictions,
                             deltas)
    end

    current_d_criterion = d_criterion(model_matrix)
    @debug "Final D-Criterion" current_d_criterion

    model_matrix
end

"""
$(TYPEDSIGNATURES)
"""
function measure_optimize_design_refresh(refresh_candidate_set)
    factors = 4
    levels = 2
    experiments = 12
    candidate_set_size = 1000
    iterations = 40

    design = generate_random_design(factors, experiments)
    candidate_set = generate_random_design(factors, candidate_set_size)

    optimized_design = optimize_design(factors = factors,
                                       levels = levels,
                                       experiments = experiments,
                                       design = design,
                                       candidate_set = candidate_set,
                                       iterations = iterations,
                                       refresh_candidate_set = refresh_candidate_set)
end

"""
$(TYPEDSIGNATURES)
"""
function measure_optimize_design(levels)
    factors = 4
    experiments = 6
    iterations = 40
    refresh_candidate_set = false

    candidate_set = expanded_design(factors, levels)
    removed_indices = randperm(levels ^ factors)[1:experiments]

    design = candidate_set[removed_indices, :]
    candidate_set = candidate_set[[i for i in 1:(levels ^ factors) if !(i in removed_indices)], :]

    optimized_design = optimize_design(factors = factors,
                                       levels = levels,
                                       experiments = experiments,
                                       design = design,
                                       candidate_set = candidate_set,
                                       iterations = iterations,
                                       refresh_candidate_set = refresh_candidate_set)
end

seed = 403
Random.seed!(seed)

design = measure_optimize_design_refresh(false)
design_refresh = measure_optimize_design_refresh(true)
design_expanded = measure_optimize_design(2)


text:
module Theories
export Context, UnlabeledContext, LabeledContext,
  Term, Var, Appl,
  Operation, Equation, LabeledEquation, UnlabeledEquation,
  Presentation, UnlabeledPresentation, LabeledPresentation,
  sort_check, check_presentation

using MLStyle
using ..LabeledVectors
import ..LabeledVectors: label, unlabel

# Ideally we would have
#
# struct Context{V}
#   vartypes::V{Int}
# end
#
# const UnlabeledContext = Context{Vector}
# const LabeledContext{S} = Context{LabeledVector{S}}
#
# However Julia does not have support for higher-kinded types

abstract type Context end

struct UnlabeledContext <: Context
  vartypes::Vector{Int}
end

struct LabeledContext{S} <: Context
  vartypes::LabeledVector{S,Int}
end

unlabel(c::LabeledContext) = UnlabeledContext(unlabel(c.vartypes))

label(c::UnlabeledContext, labels::Vector{S}) where {S} = LabeledContext{S}(label(c.vartypes, labels))

@data Term begin
  Var(index::Int)
  Appl(op::Int, args::Vector{Term})
end

struct Operation
  arity::Vector{Int}
  ret::Int
end

abstract type Equation end

struct LabeledEquation{S} <: Equation
  context::LabeledContext{S}
  lhs::Term
  rhs::Term
end

struct UnlabeledEquation <: Equation
  context::UnlabeledContext
  lhs::Term
  rhs::Term
end

unlabel(eq::LabeledEquation) = UnlabeledEquation(unlabel(eq.context), eq.lhs, eq.rhs)

label(eq::UnlabeledEquation, labels::Vector{S}) where {S} =
  LabeledEquation{S}(label(eq.context, labels), eq.lhs, eq.rhs)

abstract type Presentation end

struct UnlabeledPresentation <: Presentation
  sorts::Vector{Tuple{}}
  operations::Vector{Operation}
  equations::Vector{UnlabeledEquation}
end

struct LabeledPresentation{S} <: Presentation
  sorts::LabeledVector{S, Tuple{}}
  operations::LabeledVector{S, Operation}
  equations::LabeledVector{S, LabeledEquation{S}}
end

unlabel(p::LabeledPresentation) =
  UnlabeledPresentation(unlabel(p.sorts), unlabel(p.operations), unlabel.(unlabel(p.equations)))

function label(p::UnlabeledPresentation; sorts::Vector{S}, operations::Vector{S},
               equations::Vector{Tuple{S, Vector{S}}}) where {S}
  LabeledPresentation{S}(
    label(p.sorts, sorts),
    label(p.operations, operations),
    label(map(
      ((eq, ctxlabels)) -> label(eq, ctxlabels),
      zip(p.equations, last.(equations))
    ), first.(equations))
  )
end

function sort_check(p::Presentation, c::Context, t::Term)
  @match t begin
    Var(index) => begin
      @assert 1 <= index <= length(c.vartypes)
      c.vartypes[index]
    end
    Appl(op,args) => begin
      @assert 1 <= op <= length(p.operations)
      @assert p.operations[op].arity == map(arg -> sort_check(p,c,arg), args)
      p.operations[op].ret
    end
  end
end

function check_presentation(p::Presentation)
  for o in p.operations
    @assert 1 <= o.ret <= length(p.sorts)
    for i in o.arity
      @assert 1 <= i <= length(p.sorts)
    end
  end
  for e in p.equations
    @assert sort_check(p, e.context, e.lhs) == sort_check(p, e.context, e.lhs)
  end
  true
end

end


text:
@testset "Regressive product" begin
    @test dual(dual(A) ∧ dual(B)) == A ∨ B
end


text:
module ShowEach

export @showeach, gather_stdout, set_loop_limit, reset_indent

@enum Position Top TopFunction Middle MiddleFunction

const origSTDOUT = STDOUT

const CompactIO = (IOContext(origSTDOUT, limit=true, compact=true, displaysize=(10,20)))

type PrintState
    level::Int
    pos::Vector{Position}
    count::Vector{Int}
    loop_limit::Int
    redirect::Bool
    io::IO    # Does this matter for performance?
end

#PrintState{T<:IO}(level, pos, redirect, io::T) = PrintState{T}(level, pos, redirect, io)

const state = PrintState(0, [], [], 6, false, STDOUT)

inc(pos::Position) = (state.level += 1; push!(state.pos, pos); push!(state.count, 0))
inc() = inc(Top)
dec() = (state.level -= 1; pop!(state.pos); pop!(state.count))

inc_counter() = (state.level > 0 && state.pos[end] != MiddleFunction &&  (state.count[end] += 1))

isloop() = (state.level > 0 && any(state.count .>= state.loop_limit))

function reset_indent()
    state.level = 0
    state.pos = []
    state.count = []
end

gather_stdout(b::Bool) = (reset_indent(); state.redirect = b)
set_loop_limit(n::Int) = (reset_indent(); state.loop_limit = n)

function indent_all_but_top()
    ret = ""
    for (i, pos) in enumerate(state.pos[1:end-1])
        if pos == Middle
            ret = ret*"│  "
        else
            ret = ret*"   "
        end
    end
    ret
end

function indent(str::String)
    ret = indent_all_but_top()
    if state.level == 0 || state.pos[end] in [TopFunction, MiddleFunction]
        return ret*"   "*str
    else
        return ret*"│  "*str
    end
end

function indent()
    ret = indent_all_but_top()
    if state.level > 0
        if state.pos[end] == Top
            ret = ret*"┌─ "
            state.pos[end] = Middle
        elseif state.pos[end] == TopFunction
            ret = ret*"   "
            state.pos[end] = MiddleFunction
        elseif state.pos[end] == MiddleFunction
            ret = ret*"   "
        else
            inc_counter()
            ret = ret*"├─ "
        end
    end
    return ret
end

function show_in_tree(expr::String, value)
    if !isloop()
        value_repr = IOBuffer()
        show(IOContext(value_repr, CompactIO), "text/plain", value)
        l = length(expr) + 3
        output = replace(String(value_repr), "\n", "\n"*indent(" "^l))
        println(state.io, indent(), expr, " = ", output)
        if isloop()
            println(state.io, indent(), " ⋮")
        end
    else
        nothing
    end
end

function show_vector_short(io::IO, v, opn, cls)
    compact, prefix = Base.array_eltype_show_how(v)
    limited = get(io, :limit, false)
    if compact && !haskey(io, :compact)
        io = IOContext(io, :compact => compact)
    end
    print(io, prefix)
    if limited && length(v) > 6
        inds = Base.indices1(v)
        Base.show_delim_array(io, v, opn, ",", "", false, inds[1], inds[1]+2)
        print(io, "  \u2026  ")
        Base.show_delim_array(io, v, "", ",", cls, false, inds[end-2], inds[end])
    else
        Base.show_delim_array(io, v, opn, ",", cls, false)
    end
end

function show_in_tree(expr::String, value::Vector)
    if !isloop()
        value_repr = IOBuffer()
        show_vector_short(IOContext(value_repr, CompactIO), value, "[", "]")
        l = length(expr) + 3
        output = replace(String(value_repr), "\n", "\n"*indent(" "^l))
        println(state.io, indent(), expr, " = ", output)
        if isloop()
            println(state.io, indent(), " ⋮")
        end
    else
        nothing
    end
end

function showexpr(expr, newexpr)
    :(inc(); value = ($newexpr); dec(); show_in_tree($(string(expr)), value); value)
end

function showeach_func(expr::Expr)
    new_inner_body = Expr(:block, :(isloop() || (idt = indent(); isloop() ? println(state.io, indent(), " ⋮") : println(state.io, idt, "in function ", $(string(expr.args[1])), ":"))), :(inc(TopFunction)), (showeach(e) for e in expr.args[2:end])...)
    new_inner_body =  temporary_redirect_stdout(new_inner_body, "function "*string(expr.args[1]))
    new_body = :(try $(new_inner_body) finally dec() end)
    return Expr(expr.head, esc(expr.args[1]), :($new_body))
end

function is_function_definition(expr::Expr)
    expr.head == :function ||
    (
        expr.head == :(=) &&
        (
            Meta.isexpr(expr.args[1], :call) ||
            (
                Meta.isexpr(expr.args[1], :(::)) && Meta.isexpr(expr.args[1].args[1], :call)
            )
        )
    )
end

function showeach(expr::Expr)
    if is_function_definition(expr)
        # function definition -> show body on function call
        return showeach_func(expr)
    elseif expr.head in [:block, :vect, :tuple, :vcat, :hcat, :row, :(:), :(...), :flatten, :return, :if, :generator, :filter]
        # show all arguments recursively
        return Expr(expr.head, (showeach(e) for e in expr.args)...)
    elseif expr.head in [:(=), :(+=), :(-=), :(*=), :(/=), :typed_vcat, :typed_hcat]
        # show all argument but the first recursively
        return Expr(expr.head, esc(expr.args[1]), (showeach(e) for e in expr.args[2:end])...)
    elseif expr.head in [:ref, :comprehension, :(&&), :(||)]
        # show expression and all arguments recursively
        newexpr = Expr(expr.head, (showeach(e) for e in expr.args)...)
        return showexpr(expr, newexpr)
    elseif  expr.head in [:call, :., :typed_comprehension]
        # show expression and all arguments but the first recursively
        newexpr = Expr(expr.head, esc(expr.args[1]), (showeach(e) for e in expr.args[2:end])...)
        return showexpr(expr, newexpr)
    elseif expr.head in [:for, :while]
        # show the kind of loop and all arguments recursively
        newexpr = Expr(expr.head, (showeach(e) for e in expr.args)...)
        return showexpr(string(expr.head) * " ... end", newexpr)
    elseif expr.head == :comparison
        # show expression and every odd argument
        newexpr = Expr(expr.head, (isodd(i) ? showeach(e) : esc(e) for (i,e) in enumerate(expr.args))...)
        return showexpr(expr, newexpr)
    else
        # just return the expression unchanged
        return esc(expr)
    end
end

showeach(e::Symbol) = :(value = $(esc(e)); show_in_tree($(string(e)), value); value)
showeach(e) = esc(e)

function temporary_redirect_stdout(expr, name = nothing)
    if state.redirect
        if name == nothing
            text = "gathered STDOUT:"
        else
            text = "STDOUT from $name:"
        end
        return quote
            if STDOUT == origSTDOUT
                (rd, wr) = redirect_stdout()
                value = begin $(expr) end
                print(" ") # work around blocking of `readavailable` for empty pipe
                close(wr)
                output = chop(String(readavailable(rd)))
                close(rd)
                redirect_stdout($origSTDOUT)
                if length(output) > 0
                    println($text)
                    print(String(output))
                end
                value
            else
                tmpSTDOUT = STDOUT
                (rd, wr) = redirect_stdout()
                value = begin $(expr) end
                print(" ") # work around blocking of `readavailable` for empty pipe
                close(wr)
                output = chop(String(readavailable(rd)))
                close(rd)
                redirect_stdout(tmpSTDOUT)
                if length(output) > 0
                    println($text)
                    print(String(output))
                end
                value
            end
        end
    else
        return expr
    end
end

macro showeach(expr)
    temporary_redirect_stdout(showeach(expr))
end

end # module


text:
# This script create default configuration files that contain the tuning parameters
# for some of the algorithms this package. We assign a default value for each
# compound name key entry in the user-specified
# "compound to parameters filename" JSON file.

# One should manually edit the generated config files according to
# their preferred tuning parameters.

import JSON
import JSON3

name_map_filename = "select_compounds"
load_path = joinpath("/home/roy/Documents/repo/NMRData/input/compound_mapping", "$(name_map_filename).json")
dict_compound_to_filename = JSON.parsefile(load_path)
compound_names = collect( key for (key,value) in dict_compound_to_filename )

##### config file for SH simulation.
function createdefaultcalibrationconfig(compound_names::Vector{String};
    max_cs_shift::Float64 = 0.05, # in ppm.
    α_relative_threshold::Float64 = 0.05,
    Δc_partition_radius::Float64 = 0.3)

    db_dict = Dict()
    for name in compound_names

        db_dict[name] = Dict("maximum chemical shift" => max_cs_shift,
            "relative amplitude threshold" => α_relative_threshold,
            "maximum Δc deviation" => Δc_partition_radius)
    end

    return db_dict
end

save_folder = "/home/roy/Documents/repo/NMRData/input/SH_configs"
save_name = "$(name_map_filename)_SH_configs_default.json"
save_path = joinpath(save_folder, save_name)

# write the default SH configurations to json.
dict_out = createdefaultSHconfig(compound_names)

stringdata = JSON.json(dict_out)

open(save_path, "w") do f
    JSON3.pretty(f, stringdata)
    println(f)
end


text:
try
    using JuliaFormatter
catch error
    if isa(error, LoadError)
        using Pkg
        Pkg.add("JuliaFormatter")
        using JuliaFormatter
    else
        rethrow()
    end
end

import JuliaFormatter
JuliaFormatter.format(".")


text: RpcApi.Statistics(498,1,322013,400,300,String["orgMutationProbs: [100,300,95,0,1,1,1]","orgCloneMutation: 1","ORGANISM_MAX_MUTATIONS_ON_CLONE: 100","orgRainMutationPeriod: 300","ORGANISM_MAX_MUTATION_PERIOD: 10000","ORGANISM_MUTATION_AMOUNT: 2","ORGANISM_MAX_MUTATION_AMOUNT: 100","orgStartAmount: 10","orgStartEnergy: 5000","ORGANISM_MAX_ENERGY: 4294967295","orgEnergySpendPeriod: 3","orgEnergySpendAmount: 1","orgRemoveWeakPeriod: 300","orgRemoveWeakAmount: 10","orgClonePeriod: 5","orgStartColor: 1","codeFuncParamAmount: 2","worldWidth: 400","worldHeight: 300","WORLD_FRAME_DELAY: 0","worldIps: 1","worldMaxOrgs: 500","worldMinOrgs: 50","worldStartEnergyDots: 1000","worldStartEnergyInDot: 65280","worldEnergyCheckPercent: 0.5","worldEnergyCheckPeriod: 5000","worldZoom: 2","backupPeriod: 1","backupAmount: 7","STAT_WIDTH: 650","STAT_HEIGHT: 500","STAT_FRAME_DELAY: 5","conStepEnergySpendPercent: 20","conServerPort: 2000","conFastServerPort: 2001","conServerIp: 127.0.0.1","conLeftServerPort: 0","conLeftServerIp: 127.0.0.1","conRightServerPort: 0","conRightServerIp: 127.0.0.1","conUpServerPort: 0","conUpServerIp: 127.0.0.1","conDownServerPort: 0","conDownServerIp: 127.0.0.1"],RpcApi.SimpleOrganism(0x000000000004e992,:(function (o::Creature.Organism,)

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
include("mul.jl")
include("muladd.jl")
include("inv.jl")
include("add.jl")


text:
"""
    solveBnB!(x::BnBSolver,y::BnBModel)

Solves the branch and bound problem with the input model and solver object.
"""
function solveBnB!(x::BnBSolver,y::BnBModel)

  # initializes counts
  k_int::Int64 = 0
  k_nod::Int64 = length(y.LBD)

  # terminates when max nodes or iteration is reach, or when node stack is empty
  while (x.Term_Check(x,y,k_int))

    # fathom nodes with lower bound greater than global upper bound
    y.LBDg = minimum(y.LBD)
    push!(y.LBDg_hist,minimum(y.LBD))

    # selects node (deletion included)
    nsBox,LBDn,UBDn,id,pos = x.Node_Select(y)

    # prints node in full verbosity mode
    print_node!(x,id,LBDn,nsBox)

    # solves preprocessing/LBD/UBD/postprocessing once to get timing right
    feas_Pre::Bool = true
    feas_Post::Bool = true
    if (k_int == 0)
      #println("pre-check")
      nsBox1 = copy(nsBox)
      yUBDg1 = copy(y.UBDg)
      yLBDg1 = copy(y.LBDg)
      kint1 = copy(k_int)
      pos1 = copy(pos)
      xopt1 = copy(x.opt)
      x.Preprocess(feas_Pre,nsBox1,yUBDg1,kint1,pos1,xopt1,LBDn,UBDn,x,y)
      println("preprocess-check")
      BD_valt, LBD_solt, LBD_feast, temp_objtL = x.Lower_Prob(nsBox1,kint1,pos1,xopt1,yUBDg1)
      println("lower-check")
      UBD_valt, UBD_solt, UBD_feast, temp_objtU = x.Upper_Prob(nsBox1,kint1,pos1,xopt1,yUBDg1)
      println("upper-check")
      x.Postprocess(feas_Post,nsBox1,kint1,pos1,xopt1,
                    temp_objtL,temp_objtU,yLBDg1,yUBDg1)
      println("postprocess-check")
    end
    feas_Pre = true
    feas_Post = true

    # performs prepocessing and times
    tic()
    feas_Pre,nsBox = x.Preprocess(feas_Pre,nsBox,y.UBDg,k_int,pos,x.opt,LBDn,UBDn,x,y)
    push!(y.Pretime,y.Pretime[end]+toq())

    UBD_feas = false
    if (feas_Pre)
      # solves & times lower bounding problem
      tic()
      LBD_val,LBD_sol,LBD_feas,temp_objL = x.Lower_Prob(nsBox,k_int,pos,x.opt,y.UBDg)
      push!(y.LBDgtime,y.LBDgtime[end]+toq())
      y.lbcnt += 1
      print_results!(x,LBD_val,LBD_sol,LBD_feas,true)
      int_info = LBD_val,LBD_sol,LBD_feas,temp_objL,nsBox

      # checks for infeasibility stores solution
      if (LBD_feas)
        println("ran LBD feas")
        println("boolean: $(x.converged(x,y.UBDg,LBD_val))")
        if (~x.converged(x,y.UBDg,LBD_val))

          # solves & times upper bounding problem
          println("ran Upper bound")
          tic()
          UBD_val,UBD_sol,UBD_feas,temp_objU = x.Upper_Prob(nsBox,k_int,pos,x.opt,y.UBDg)
          push!(y.UBDgtime,y.UBDgtime[end]+toq())
          y.ubcnt += 1
          print_results!(x,UBD_val,UBD_sol,UBD_feas,false)

          # fathoms by value dominance
          if (UBD_feas)
            # update to handle equality constraints
            if (y.lastgap>(UBD_val-LBD_val) && (UBD_val - y.UBDg < x.BnB_atol/10.0))
              y.lastgap = UBD_val-LBD_val
              y.soln = UBD_sol
            end
            if (UBD_val < y.UBDg)
              y.feas_fnd = true
              y.first_num = y.lbcnt
              y.UBDg = UBD_val
              y.soln = UBD_sol
              push!(y.UBDg_hist,UBD_val)
            else
              push!(y.UBDg_hist,y.UBDg)
            end
          else
            push!(y.UBDg_hist,y.UBDg_hist[end])
          end

          # performs post processing and times
          tic()
          feas_Post,nsBox_Post = x.Postprocess(feas_Post,nsBox,k_int,pos,x.opt,temp_objL,temp_objU,y.LBDg,y.UBDg)
          push!(y.Posttime,y.Posttime[end]+toq())

          # branch if criteria met
          if (feas_Post)
            if x.Repeat_Chk(x,y,nsBox,nsBox_Post)
              BM_Single!(x,y,LBD_val,UBD_val,nsBox,pos)
            else
              # branch the nodes & stores nodes to stack
              Y1,Y2 = x.Bisect_Func(x,y,nsBox)
              x.Branch_Sto(x,y,LBD_val,UBD_val,Y1,Y2,pos)
            end
          end

        elseif (~x.exhaust && y.feas_fnd)
        end# end converged check
      end # end LBD feasibility check

      # prints relative statistics for the iteration
      fathom!(y)
    else
      LBD_val = -Inf
      LBD_feas = false
      UBD_feas = false
    end

    print_int!(x,k_int,length(y.LBD),id,LBD_val,y.LBDg,y.UBDg,LBD_feas,UBD_feas)

    k_int+=1

  end
  y.soln_val = y.UBDg

  # prints the solution
  if ((x.Verbosity == "Full")||(x.Verbosity == "Normal"))
    print_sol!(x,y,y.ubcnt,y.lbcnt,y.UBDgtime[end],y.LBDgtime[end])
  end
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ a6abf980-2042-11eb-326f-67792f308488
begin
	using Revise

	# dataset
	using DataFrames
	using CSV

	using BenchmarkTools
	using Clustering
	using ISODATA
	using MLDataUtils
	using MLLabelUtils
	using StatsBase
end

# ╔═╡ f277c7dc-4f64-11eb-147d-95a2b10bab9f

# ╔═╡ c4a2ee6c-4f8d-11eb-24f5-e53ed7aa6789
begin
	# load iris dataset
	idf = CSV.read("data/iris.data", DataFrame, header=false)
	iris_raw = convert(Matrix, idf)
	iX, iY = iris_raw[:, 1:end-1], iris_raw[:, end]
	iY = convertlabel(
		MLLabelUtils.LabelEnc.Indices,
		iY,
		labelenc(iY))
	iX = Array{Float64}(transpose(iX))
	iris = cat(iX, transpose(iY), dims=1)
end

# ╔═╡ a1777880-5065-11eb-197d-0dba960318ed
begin
	si = summarize(iris)
	si, iso_param(iris[1:end-1, :], si)
end

# ╔═╡ 7d36be2c-4f70-11eb-3b81-2b78885c3771
begin
	# load glass dataset
	gdf = CSV.read("data/glass.data", DataFrame, header=false)
	glass_raw = convert(Matrix, gdf[:, 2:end])
	gX, gY = glass_raw[:, 1:end-1], glass_raw[:, end]
	gY = convertlabel(
		MLLabelUtils.LabelEnc.Indices,
		gY,
		labelenc(gY))
	for i = 5:7
		inds = findall(x -> x == i, gY)
		gY[inds] .= i - 1
	end
	gX = Array{Float64}(transpose(gX))
	glass = cat(gX, transpose(gY), dims=1)
end

# ╔═╡ a3565518-5065-11eb-1b6f-27696c4519a6
begin
	sg = summarize(glass)
	sg, iso_param(glass[1:end-1, :], sg)
end

# ╔═╡ b2ce85e2-4f8e-11eb-3c52-0fce3a112046
begin
	# load wine dataset
	wdf = CSV.read("data/wine.data", DataFrame, header=false)
	wine_raw = convert(Matrix, wdf)
	wX, wY = wine_raw[:, 2:end], wine_raw[:, 1]
	wine = cat(transpose(wX), transpose(wY), dims=1)
end

# ╔═╡ 666a9e28-4df1-11eb-108b-fd83c3af1a27
struct ExperimentResult
	iter::Float64
	centers::Float64
	adj_idx::Float64
	rand_idx::Float64
	mirkin_idx::Float64
	hubert_idx::Float64
end

# ╔═╡ 7a1ba842-4f6b-11eb-222b-2b41308b3d35
# experiment2: computed isodata params
function experiment2(
	data,
	k,
	num_trials,
	iseeds;
	normalize=false,
	_θn_bias::Float64=0.20,
	_θe_bias::Float64=0.20,
	_θc_bias::Float64=0.20
)
	results = Vector{ExperimentResult}(undef, num_trials)

	for i in 1:num_trials
		# train test split
		train, test = splitobs(shuffleobs(data), at = 0.8)
		train_x, train_y = train[1:end-1, :], round.(Int, train[end, :])
		test_x, test_y = test[1:end-1, :], round.(Int, test[end, :])

		if normalize
			train_z = fit(ZScoreTransform, train_x, dims=2)
			train_x = StatsBase.transform(train_z, train_x)
			test_x = StatsBase.transform(train_z, test_x)
		end

		# calculate isodata parameters
		summaries = summarize(train)
		θn, θe, θc = iso_param(
			train_x,
			summaries,
			_θn_bias=_θn_bias,
			_θe_bias=_θe_bias,
			_θc_bias=_θc_bias)

		# run isodata
		res = base_conv(train_x, k, θn=θn, θe=θe, θc=θc, iseeds=iseeds)

		# fit centers to test data
		fit_y = clusterfit(test_x, res.centers)

		# cluster analysis
		adj_idx, rand_idx, mirkin_idx, hubert_idx = randindex(test_y, fit_y)
		num_centers = size(res.centers, 2)

		# save results
		results[i] = ExperimentResult(
			res.iter, num_centers, adj_idx, rand_idx, mirkin_idx, hubert_idx)
	end

	iters = 0.0
	centers = 0.0
	adjs = 0.0
	rands = 0.0
	mirkins = 0.0
	huberts = 0.0
	for i in 1:num_trials
		iters += results[i].iter
		centers += results[i].centers
		adjs += results[i].adj_idx
		rands += results[i].rand_idx
		mirkins += results[i].mirkin_idx
		huberts += results[i].hubert_idx
	end

	avg_iters = iters/num_trials
	avg_centers = centers/num_trials
	avg_adjs = adjs/num_trials
	avg_rands = rands/num_trials
	avg_mirkins = mirkins/num_trials
	avg_huberts = huberts/num_trials
	return ExperimentResult(
		avg_iters, avg_centers, avg_adjs, avg_rands, avg_mirkins, avg_huberts)
end

# ╔═╡ ebed5ca6-5068-11eb-36f1-51e85780d0f1
"EXPERIMENT 2: ISODATA - IRIS"

# ╔═╡ c43a6cda-4f6b-11eb-3480-910efacac2e6
experiment2(iris, 3, 1000, isodata_s)

# ╔═╡ 06f60a46-4f6c-11eb-141f-6d444b5d56a7
experiment2(iris, 3, 1000, kmppcenters)

# ╔═╡ 4d3eb44e-4f6c-11eb-19e5-55f8d6eb03ad
experiment2(iris, 3, 1000, randomcenters)

# ╔═╡ e20c5c64-5068-11eb-3a40-d5c965657350
"EXPERIMENT 2: ISODATA - GLASS"

# ╔═╡ dd40e1d6-4f71-11eb-026c-6fc896e10fcd
experiment2(glass, 6, 1000, isodata_s, normalize=false)

# ╔═╡ 3452e904-4f72-11eb-288c-4d397ede82ab
experiment2(glass, 6, 1000, kmppcenters, normalize=false)

# ╔═╡ 2f78af8a-4f73-11eb-227e-8733704a70d3
experiment2(glass, 6, 1000, randomcenters, normalize=false)

# ╔═╡ 4b011c6c-4f8f-11eb-0c84-1d844e23378b
"EXPERIMENT 2: ISODATA - WINE"

# ╔═╡ 2efc4230-4f8f-11eb-2b8e-85ee78f4ebaf
experiment2(wine, 3, 1000, isodata_s)

# ╔═╡ 3c6826a0-4f8f-11eb-3912-81f43868d9f3
experiment2(wine, 3, 1000, kmppcenters)

# ╔═╡ 44fbc7fe-4f8f-11eb-1003-bb016d171cc9
experiment2(wine, 3, 1000, randomcenters)

# ╔═╡ 226e09e6-4f93-11eb-2fc9-d15f03f16ddb
# kmeans with new cluster center
function experiment3(
	data,
	k,
	num_trials,
	iseeds
)
	results = Vector{ExperimentResult}(undef, num_trials)

	for i in 1:num_trials
		# train test split
		train, test = splitobs(shuffleobs(data), at = 0.8)
		train_x, train_y = train[1:end-1, :], round.(Int, train[end, :])
		test_x, test_y = test[1:end-1, :], round.(Int, test[end, :])

		# normalize
		train_z = fit(ZScoreTransform, train_x, dims=2)
		train_x = StatsBase.transform(train_z, train_x)
		test_x = StatsBase.transform(train_z, test_x)

		# run isodata
		res = ISODATA.kmeans(train_x, k, iseeds=iseeds)

		# fit centers to test data
		fit_y = clusterfit(test_x, res.centers)

		# cluster analysis
		adj_idx, rand_idx, mirkin_idx, hubert_idx = randindex(test_y, fit_y)
		num_centers = size(res.centers, 2)

		# save results
		results[i] = ExperimentResult(
			res.iter, num_centers, adj_idx, rand_idx, mirkin_idx, hubert_idx)
	end

	iters = 0.0
	centers = 0.0
	adjs = 0.0
	rands = 0.0
	mirkins = 0.0
	huberts = 0.0
	for i in 1:num_trials
		iters += results[i].iter
		centers += results[i].centers
		adjs += results[i].adj_idx
		rands += results[i].rand_idx
		mirkins += results[i].mirkin_idx
		huberts += results[i].hubert_idx
	end

	avg_iters = iters/num_trials
	avg_centers = centers/num_trials
	avg_adjs = adjs/num_trials
	avg_rands = rands/num_trials
	avg_mirkins = mirkins/num_trials
	avg_huberts = huberts/num_trials
	return ExperimentResult(
		avg_iters, avg_centers, avg_adjs, avg_rands, avg_mirkins, avg_huberts)
end

# ╔═╡ f81844c8-5068-11eb-1a4a-c3563bd09221
"EXPERIMENT 3: KMEANS - IRIS"

# ╔═╡ 505daeb0-4f93-11eb-1671-e361fdb821e9
experiment3(iris, 3, 1000, isodata_s)

# ╔═╡ 60697ad2-4f93-11eb-017f-69c9c730506f
experiment3(iris, 3, 1000, kmppcenters)

# ╔═╡ 60fd52de-4f93-11eb-2f43-e797797924af
experiment3(iris, 3, 1000, randomcenters)

# ╔═╡ f9a7d948-5068-11eb-069f-27055b9d6bf7
"EXPERIMENT 3: KMEANS - GLASS"

# ╔═╡ 5601a86c-5068-11eb-2981-65b6287e8ce9
experiment3(glass, 6, 1000, isodata_s)

# ╔═╡ 5e47cb8e-5068-11eb-0459-ff6e67a2fcd4
experiment3(glass, 6, 1000, kmppcenters)

# ╔═╡ 5e332274-5068-11eb-2f5c-1f89f5d6da9e
experiment3(glass, 6, 1000, randomcenters)

# ╔═╡ fad939fe-5068-11eb-26d9-351ad4e57198
"EXPERIMENT 3: KMEANS - WINE"

# ╔═╡ 850b91f4-4f93-11eb-2447-071b5f7d8cd1
experiment3(wine, 3, 1000, isodata_s)

# ╔═╡ 84d58f32-4f93-11eb-2224-478298799589
experiment3(wine, 3, 1000, kmppcenters)

# ╔═╡ 84554f52-4f93-11eb-296a-c567ad3db9c7
experiment3(wine, 3, 1000, randomcenters)

# ╔═╡ Cell order:
# ╠═a6abf980-2042-11eb-326f-67792f308488
# ╟─f277c7dc-4f64-11eb-147d-95a2b10bab9f
# ╠═c4a2ee6c-4f8d-11eb-24f5-e53ed7aa6789
# ╠═a1777880-5065-11eb-197d-0dba960318ed
# ╠═7d36be2c-4f70-11eb-3b81-2b78885c3771
# ╠═a3565518-5065-11eb-1b6f-27696c4519a6
# ╠═b2ce85e2-4f8e-11eb-3c52-0fce3a112046
# ╠═666a9e28-4df1-11eb-108b-fd83c3af1a27
# ╠═7a1ba842-4f6b-11eb-222b-2b41308b3d35
# ╟─ebed5ca6-5068-11eb-36f1-51e85780d0f1
# ╠═c43a6cda-4f6b-11eb-3480-910efacac2e6
# ╠═06f60a46-4f6c-11eb-141f-6d444b5d56a7
# ╠═4d3eb44e-4f6c-11eb-19e5-55f8d6eb03ad
# ╟─e20c5c64-5068-11eb-3a40-d5c965657350
# ╠═dd40e1d6-4f71-11eb-026c-6fc896e10fcd
# ╠═3452e904-4f72-11eb-288c-4d397ede82ab
# ╠═2f78af8a-4f73-11eb-227e-8733704a70d3
# ╟─4b011c6c-4f8f-11eb-0c84-1d844e23378b
# ╠═2efc4230-4f8f-11eb-2b8e-85ee78f4ebaf
# ╠═3c6826a0-4f8f-11eb-3912-81f43868d9f3
# ╠═44fbc7fe-4f8f-11eb-1003-bb016d171cc9
# ╠═226e09e6-4f93-11eb-2fc9-d15f03f16ddb
# ╟─f81844c8-5068-11eb-1a4a-c3563bd09221
# ╠═505daeb0-4f93-11eb-1671-e361fdb821e9
# ╠═60697ad2-4f93-11eb-017f-69c9c730506f
# ╠═60fd52de-4f93-11eb-2f43-e797797924af
# ╟─f9a7d948-5068-11eb-069f-27055b9d6bf7
# ╠═5601a86c-5068-11eb-2981-65b6287e8ce9
# ╠═5e47cb8e-5068-11eb-0459-ff6e67a2fcd4
# ╠═5e332274-5068-11eb-2f5c-1f89f5d6da9e
# ╟─fad939fe-5068-11eb-26d9-351ad4e57198
# ╠═850b91f4-4f93-11eb-2447-071b5f7d8cd1
# ╠═84d58f32-4f93-11eb-2224-478298799589
# ╠═84554f52-4f93-11eb-296a-c567ad3db9c7


text:
## Stability derivative cases
function scale_inputs(freestream :: Freestream, b, c)
    # Building input vector
    x0 = [ freestream.alpha;
           freestream.beta;
           rate_coefficient(freestream.omega, freestream.V, b, c) ]

    # Unscaling non-dimensional rate coefficients
    scale = 2freestream.V .* [1/b, 1/c, 1/b]

    x0, scale
end

function print_case(nf, ff, derivs, comp)
    print_coefficients(nf, ff, comp)
    print_derivatives(derivs, comp)
end

function print_case(data, comp)
    nf, ff, derivs = data[comp]
    print_case(nf, ff, derivs, comp)
end

function solve_stability_case(wing :: Union{Wing, HalfWing}, freestream :: Freestream; rho_ref = 1.225, r_ref = zeros(3), area_ref = projected_area(wing), chord_ref = mean_aerodynamic_chord(wing), span_ref = span(wing), span_num :: Union{Integer, Vector{<: Integer}}, chord_num :: Integer, name = "Wing", viscous = false, x_tr = 0.3, print = false, spacing = ["sine"; fill("cosine", length(span_num) - 1)])
    # Reference values and scaling inputs
    S, b, c = area_ref, span_ref, chord_ref
    x, scale = scale_inputs(freestream, b, c)

    # Closure to generate results with input vector
    function stab(x)
        fs 	 = Freestream(freestream.V, rad2deg(x[1]), rad2deg(x[2]), x[3:end] .* scale)
        data = solve_case(wing, fs;
                          rho_ref   = rho_ref,
                          r_ref     = r_ref,
                          area_ref  = S,
                          span_ref  = b,
                          chord_ref = c,
                          chord_num = chord_num,
                          span_num  = span_num,
                          viscous 	= viscous,
                          x_tr 		= x_tr,
                          spacing   = spacing)

        [ data[1]; data[2] ] # Concatenate nearfield and farfield coefficients for DiffResults value
    end

    # Set up Jacobian system and evaluate
    y = ifelse(viscous, zeros(16), zeros(12))
    result = DiffResults.JacobianResult(y, x)
    result = ForwardDiff.jacobian!(result, stab, x)
    vars   = DiffResults.value(result)
    derivs = DiffResults.jacobian(result)

    # Gather results
    nf 	= ifelse(viscous, vars[1:11], vars[1:9]) 	  	# Nearfield coefficients
    ff 	= ifelse(viscous, vars[12:end], vars[10:end]) 	# Farfield coefficients
    dvs = ifelse(viscous, derivs[[1,4:8...],:], derivs[1:6,:])	# Nearfield stability derivatives, uses total drag for viscous cases just in case of generalisations

    # Print if needed
    if print; print_case(nf, ff, derivs, name) end

    nf, ff, dvs
end

function solve_stability_case(aircraft :: Dict{String, Tuple{Matrix{Panel3D{T}}, Matrix{SVector{3,T}}}}, freestream :: Freestream; rho_ref = 1.225, r_ref = zeros(3), area_ref = 1, chord_ref = 1, span_ref = 1, name = "Aircraft", print = false, print_components = false) where T <: Real
    # Reference values and scaling inputs
    S, b, c = area_ref, span_ref, chord_ref
    x, scale = scale_inputs(freestream, b , c)

    # Closure to generate results with input vector
    function stab(x)
        fs 	 = Freestream(freestream.V, rad2deg(x[1]), rad2deg(x[2]), x[3:end] .* scale)
        data = solve_case(aircraft, fs,
                          rho_ref   = rho_ref,
                          r_ref     = r_ref,
                          area_ref  = S,
                          span_ref  = b,
                          chord_ref = c,
                          name 		= name)

        # Creates array of nearfield and farfield coefficients for each component as a row vector.
        coeffs = reduce(hcat, let comps = data[name]; [ comps[1]; comps[2] ] end for name in names)
    end

    names 	= [ reduce(vcat, keys(aircraft)); name ]
    num_comps = length(names)

    y 		= zeros(12, num_comps)
    result 	= DiffResults.JacobianResult(y, x)
    result 	= ForwardDiff.jacobian!(result, stab, x)

    vars 	= DiffResults.value(result)
    derivs 	= DiffResults.jacobian(result)

    # Dictionary assembly
    ranges  = (1:num_comps) .* 12                     # Painful hacking
    bounds  = zip([ 1; ranges[1:end-1] .+ 1], ranges) # Painful hacking
    data 	= Dict(name => (vars[1:9, i], vars[10:end, i], derivs[first(inds):last(inds)-6,:]) for (i, (name, inds)) in (enumerate ∘ zip)(names, bounds))

    # Printing
    if print;            print_case(data, name)       end
    if print_components; print_case.(Ref(data), names)end

    data
end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Unicode

@testset "replace" begin
    @test replace("\u2202", '*' => '\0', ""=>"") == "\u2202"

    @test replace("foobar", 'o' => '0', ""=>"") == "f00bar"
    @test replace("foobar", 'o' => '0', count=1, ""=>"") == "f0obar"
    @test replace("foobar", 'o' => "", ""=>"") == "fbar"
    @test replace("foobar", 'o' => "", count=1, ""=>"") == "fobar"
    @test replace("foobar", 'f' => 'F', ""=>"") == "Foobar"
    @test replace("foobar", 'r' => 'R', ""=>"") == "foobaR"

    @test replace("foofoofoo", "foo" => "bar", ""=>"") == "barbarbar"
    @test replace("foobarfoo", "foo" => "baz", ""=>"") == "bazbarbaz"
    @test replace("barfoofoo", "foo" => "baz", ""=>"") == "barbazbaz"

    @test replace("", "" => "", ""=>"") == ""
    @test replace("", "" => "x", ""=>"") == "x"
    @test replace("", "x" => "y", ""=>"") == ""

    @test replace("abcd", "" => "^", ""=>"") == "^a^b^c^d^"
    @test replace("abcd", "b" => "^", ""=>"") == "a^cd"
    @test replace("abcd", r"b?" => "^", ""=>"") == "^a^c^d^"
    @test replace("abcd", r"b+" => "^", ""=>"") == "a^cd"
    @test replace("abcd", r"b?c?" => "^", ""=>"") == "^a^d^"
    @test replace("abcd", r"[bc]?" => "^", ""=>"") == "^a^^d^"

    @test replace("foobarfoo", r"(fo|ba)" => "xx", ""=>"") == "xxoxxrxxo"
    @test replace("foobarfoo", r"(foo|ba)" => "bar", ""=>"") == "barbarrbar"

    @test replace("foobar", 'o' => 'ø', ""=>"") == "føøbar"
    @test replace("foobar", 'o' => 'ø', count=1, ""=>"") == "føobar"
    @test replace("føøbar", 'ø' => 'o', ""=>"") == "foobar"
    @test replace("føøbar", 'ø' => 'o', count=1, ""=>"") == "foøbar"
    @test replace("føøbar", 'ø' => 'ö', ""=>"") == "fööbar"
    @test replace("føøbar", 'ø' => 'ö', count=1, ""=>"") == "föøbar"
    @test replace("føøbar", 'ø' => "", ""=>"") == "fbar"
    @test replace("føøbar", 'ø' => "", count=1, ""=>"") == "føbar"
    @test replace("føøbar", 'f' => 'F', ""=>"") == "Føøbar"
    @test replace("ḟøøbar", 'ḟ' => 'F', ""=>"") == "Føøbar"
    @test replace("føøbar", 'f' => 'Ḟ', ""=>"") == "Ḟøøbar"
    @test replace("ḟøøbar", 'ḟ' => 'Ḟ', ""=>"") == "Ḟøøbar"
    @test replace("føøbar", 'r' => 'R', ""=>"") == "føøbaR"
    @test replace("føøbaṙ", 'ṙ' => 'R', ""=>"") == "føøbaR"
    @test replace("føøbar", 'r' => 'Ṙ', ""=>"") == "føøbaṘ"
    @test replace("føøbaṙ", 'ṙ' => 'Ṙ', ""=>"") == "føøbaṘ"

    @test replace("ḟøøḟøøḟøø", "ḟøø" => "bar", ""=>"") == "barbarbar"
    @test replace("ḟøøbarḟøø", "ḟøø" => "baz", ""=>"") == "bazbarbaz"
    @test replace("barḟøøḟøø", "ḟøø" => "baz", ""=>"") == "barbazbaz"

    @test replace("foofoofoo", "foo" => "ƀäṙ", ""=>"") == "ƀäṙƀäṙƀäṙ"
    @test replace("fooƀäṙfoo", "foo" => "baz", ""=>"") == "bazƀäṙbaz"
    @test replace("ƀäṙfoofoo", "foo" => "baz", ""=>"") == "ƀäṙbazbaz"

    @test replace("foofoofoo", "foo" => "bar", ""=>"") == "barbarbar"
    @test replace("foobarfoo", "foo" => "ƀäż", ""=>"") == "ƀäżbarƀäż"
    @test replace("barfoofoo", "foo" => "ƀäż", ""=>"") == "barƀäżƀäż"

    @test replace("ḟøøḟøøḟøø", "ḟøø" => "ƀäṙ", ""=>"") == "ƀäṙƀäṙƀäṙ"
    @test replace("ḟøøƀäṙḟøø", "ḟøø" => "baz", ""=>"") == "bazƀäṙbaz"
    @test replace("ƀäṙḟøøḟøø", "ḟøø" => "baz", ""=>"") == "ƀäṙbazbaz"

    @test replace("ḟøøḟøøḟøø", "ḟøø" => "bar", ""=>"") == "barbarbar"
    @test replace("ḟøøbarḟøø", "ḟøø" => "ƀäż", ""=>"") == "ƀäżbarƀäż"
    @test replace("barḟøøḟøø", "ḟøø" => "ƀäż", ""=>"") == "barƀäżƀäż"

    @test replace("ḟøøḟøøḟøø", "ḟøø" => "ƀäṙ", ""=>"") == "ƀäṙƀäṙƀäṙ"
    @test replace("ḟøøƀäṙḟøø", "ḟøø" => "ƀäż", ""=>"") == "ƀäżƀäṙƀäż"
    @test replace("ƀäṙḟøøḟøø", "ḟøø" => "ƀäż", ""=>"") == "ƀäṙƀäżƀäż"

    @test replace("", "" => "ẍ", ""=>"") == "ẍ"
    @test replace("", "ẍ" => "ÿ", ""=>"") == ""

    @test replace("äƀçđ", "" => "π", ""=>"") == "πäπƀπçπđπ"
    @test replace("äƀçđ", "ƀ" => "π", ""=>"") == "äπçđ"
    @test replace("äƀçđ", r"ƀ?" => "π", ""=>"") == "πäπçπđπ"
    @test replace("äƀçđ", r"ƀ+" => "π", ""=>"") == "äπçđ"
    @test replace("äƀçđ", r"ƀ?ç?" => "π", ""=>"") == "πäπđπ"
    @test replace("äƀçđ", r"[ƀç]?" => "π", ""=>"") == "πäππđπ"

    @test replace("foobarfoo", r"(fo|ba)" => "ẍẍ", ""=>"") == "ẍẍoẍẍrẍẍo"

    @test replace("ḟøøbarḟøø", r"(ḟø|ba)" => "xx", ""=>"") == "xxøxxrxxø"
    @test replace("ḟøøbarḟøø", r"(ḟøø|ba)" => "bar", ""=>"") == "barbarrbar"

    @test replace("fooƀäṙfoo", r"(fo|ƀä)" => "xx", ""=>"") == "xxoxxṙxxo"
    @test replace("fooƀäṙfoo", r"(foo|ƀä)" => "ƀäṙ", ""=>"") == "ƀäṙƀäṙṙƀäṙ"

    @test replace("ḟøøƀäṙḟøø", r"(ḟø|ƀä)" => "xx", ""=>"") == "xxøxxṙxxø"
    @test replace("ḟøøƀäṙḟøø", r"(ḟøø|ƀä)" => "ƀäṙ", ""=>"") == "ƀäṙƀäṙṙƀäṙ"

    @test replace("foo", "oo" => uppercase, ""=>"") == "fOO"

    # Issue 13332
    @test replace("abc", 'b' => 2.1, ""=>"") == "a2.1c"

    # test replace with a count for String and GenericString
    # check that replace is a no-op if count==0
    for s in ["aaa", Test.GenericString("aaa")]
        @test_throws DomainError replace(s, 'a' => "", count = -1, ""=>"")
        @test replace(s, 'a' => 'z', count=0, ""=>"") === s
        @test replace(s, 'a' => 'z', count=1, ""=>"") == "zaa"
        @test replace(s, 'a' => 'z', count=2, ""=>"") == "zza"
        @test replace(s, 'a' => 'z', count=3, ""=>"") == "zzz"
        @test replace(s, 'a' => 'z', count=4, ""=>"") == "zzz"
        @test replace(s, 'a' => 'z', count=typemax(Int), ""=>"") == "zzz"
        @test replace(s, 'a' => 'z', ""=>"")    == "zzz"
    end

    for s in ["abc"]
        @test replace(s) === s
        @test replace(s, 'a' => 'z', ""=>"") === "zbc"
        @test replace(s, 'a' => 'z', 'b' => 'y') == "zyc"
        @test replace(s, 'a' => 'z', 'c' => 'x', "b" => 'y') == "zyx"
        @test replace(s, '1' => 'z', ""=>"") == s
        @test replace(s, 'b' => "BbB", ""=>"", count=1) == "aBbBc"
    end

    for s in ["quick quicker quickest"]
        @test replace(s) === s
        @test replace(s, "quick" => 'a', "quicker" => uppercase, "quickest" => 'z') == "a QUICKER z"
        @test replace(s, "quick"=>"Duck", "quicker"=>"is", "quickest"=>"lame", count=2) == "Duck is quickest"
        @test replace(s, "" => '1', ""=>"") == "1q1u1i1c1k1 1q1u1i1c1k1e1r1 1q1u1i1c1k1e1s1t1"
        @test replace(s, "qu" => "QU", "qu" => "never happens", "ick" => "") == "QU QUer QUest"
        @test replace(s, " " => '_', "r " => "r-") == "quick_quicker-quickest"
        @test replace(s, r"[aeiou]" => "ä", "ui" => "ki", "i" => "I") == "qkick qkickär qkickäst"
        @test replace(s, r"[^ ]+" => "word", "quicker " => "X", count=big"99") == "word Xword"

        @test replace(s, r"(quick)(e)"=>s"\2-\1", "x"=>"X") == "quick e-quickr e-quickst"

        @test replace(s, 'q'=>'Q', 'u'=>'U') == "QUick QUicker QUickest"
        @test replace(s, 'q'=>'Q', r"u"=>'U') == "QUick QUicker QUickest"
        @test replace(s, 'q'=>'Q', equalto('u')=>uppercase) == "QUick QUicker QUickest"
        @test replace(s, 'q'=>'Q', islower=>'-') == "Q---- Q------ Q-------"
        @test replace(s, ['q', 'u']=>'K') == "KKick KKicker KKickest"
        @test replace(s, occursin("uq")=>'K') == "KKick KKicker KKickest"
        @test replace(s, equalto('q')=>"B") == "Buick Buicker Buickest"

        @test replace(s, "qui"=>"A", 'r'=>'R') == "Ack AckeR Ackest"
        @test replace(s, 'r'=>'x', islower=>uppercase) == "QUICK QUICKEx QUICKEST"
        @test replace(s, islower=>uppercase, 'r'=>'x') == "QUICK QUICKER QUICKEST"
        @test replace(s, "q"=>"z", islower=>uppercase, 'r'=>'x') == "zUICK zUICKER zUICKEST"
        @test replace(s, "qui"=>"A", 'r'=>'x', islower=>uppercase) == "ACK ACKEx ACKEST"
        @test replace(s, "qui"=>"A", 'r'=>'x', islower=>uppercase) == "ACK ACKEx ACKEST"
        @test replace(s, r"q"=>"z", islower=>uppercase, 'r'=>'x') == "zUICK zUICKER zUICKEST"
        @test_throws ErrorException("type String has no field match_data") replace(s, "q"=>s"a\1b")
        @test_throws ErrorException("PCRE error: unknown substring") replace(s, r"q"=>s"a\1b")
    end

end


text:
export RationalFunctions,
    RationalBasisFunction

"A basis of functions of the form `1/(z-z_j)` where `z_j` are the poles."
struct RationalFunctions{T} <: Dictionary{T,T}
    poles   :: Vector{T}
end

RationalFunctions(n::Int) = RationalFunctions{Float64}(n)
RationalFunctions{T}(n::Int) where {T} = RationalFunctions(rand(T, n))

size(dict::RationalFunctions) = (length(dict.poles),)

pole(dict::RationalFunctions, idx) = dict.poles[idx]

unsafe_eval_element(dict::RationalFunctions, idx::Int, x) = 1/(x-dict.poles[idx])

support(dict::RationalFunctions{T}) where {T} = DomainSets.FullSpace{T}()

function unsafe_eval_element_derivative(dict::RationalFunctions, idx, x, order)
    @assert order == 1
    -1/(x-dict.poles[idx])^2
end

function similar(dict::RationalFunctions, ::Type{T}, n::Int) where {T}
    @assert n == length(dict)
    RationalFunctions{T}(dict.poles)
end

struct RationalBasisFunction{T} <: Polynomial{T}
    pole    ::  T
end

pole(r::RationalBasisFunction) = r.pole

RationalBasisFunction{T}(r::RationalBasisFunction) where {T} = RationalBasisFunction{T}(r.pole)

show(io::IO, d::RationalBasisFunction{T}) where {T<:Real} =
    print(io, "1/(x-$(pole(r)))  (rational function)")
show(io::IO, d::RationalBasisFunction{T}) where {T<:Complex} =
    print(io, "1/(z-$(pole(r)))  (rational function)")

convert(::Type{TypedFunction{T,T}}, r::RationalBasisFunction) where {T} = RationalBasisFunction{T}(r.pole)

(r::RationalBasisFunction)(x) = 1/(x-pole(r))

basisfunction(dict::RationalFunctions{T}, idx) where {T} = RationalBasisFunction{T}(pole(dict, idx))


text:
"""
    read_tetgen

Read mesh from Tetgen.
"""
function read_tetgen(filename)

    @info "Reading Tetgen FE mesh from " * filename

    # Read points
    points = open(filename * ".node", "r") do io
        npoint = parse(Int, split(readline(io))[1])
        points = zeros(3, npoint)
        for ipoint = 1:npoint
            points[:, ipoint] = parse.(Float64, split(readline(io))[2:4])
        end
        points
    end

    # Read facets and their associated boundaries
    facets, facetmarkers = open(filename * ".face", "r") do io
        nfacet = parse(Int, split(readline(io))[1])
        facets = zeros(Int, 3, nfacet)
        facetmarkers = zeros(Int, nfacet)
        for ifacet = 1:nfacet
            tmp = parse.(Int, split(readline(io)))
            facets[:, ifacet] = tmp[2:4]
            facetmarkers[ifacet] = tmp[5]
        end
        facets, facetmarkers
    end

    # Read elements and their associated compartments
    elements, elementmarkers = open(filename * ".ele", "r") do io
        nelement = parse(Int, split(readline(io))[1])
        elements = zeros(Int, 4, nelement)
        elementmarkers = zeros(Int, nelement)
        for ielement = 1:nelement
            tmp = parse.(Int, split(readline(io)))
            elements[:, ielement] = tmp[2:5]
            elementmarkers[ielement] = tmp[6]
        end
        elements, elementmarkers
    end

    (; points, facets, facetmarkers, elements, elementmarkers)
end


text:
using CategoricalArrays, JSON3, JLD2, Random
using MLDataUtils, StatsBase
using KCenters, KNearestCenters, TextSearch, TextClassification

function normstring(t)
    if occursin("&", t)
        t = replace(t, ">" => ">")
        t = replace(t, "<" => "<")
        t = replace(t, "&" => "&")
    end

    replace(t, r"#(USER|URL|HASHTAG)#" => p->"_" * lowercase(p[2:end-1]))
end

function loadjson(filename, normstring=normstring; textkey="text", labelkey=nothing)
    corpus = []
    labels = []

    for line in readlines(filename)
        d = JSON3.read(line)
        t = d[textkey]
        if t isa AbstractString
            t = normstring(t)
        else
            t = [normstring(t[i]) for i in eachindex(t)]
        end
        push!(corpus, t)

        if labelkey !== nothing
            push!(labels, d[labelkey])
        end
    end

    (corpus=corpus, labels=length(labels) > 0 ? categorical(labels) : nothing, filename=filename, textkey=textkey, labelkey=labelkey)
end

function stratified_split(labels; p=0.7, shuffle=true)
    n = length(labels)
    X = collect(1:n)
    (train_X, train_labels), (test_X, test_labels) = stratifiedobs((X, labels), p=p, shuffle=shuffle)
    train_X, test_X
end

function stratified_folds(labels; folds=3, shuffle=true)
    n = length(labels)
    X = collect(1:n)
    shuffle && shuffle!(X)
    p = tuple([round(1/folds, RoundUp, digits=2) for i in 1:(folds-1)]...)
    L = stratifiedobs((X, labels), p=p)
    F = []
    for i in 1:folds
        P = L[[j for j in 1:folds if i != j]]
        itrain = vcat(getindex.(P, 1)...)
        itest = L[i][1]
        push!(F, (itrain, itest))
    end

    F
end

max_stars_repo_name: sadit/text-classification-scripts

max_issues_repo_name: sadit/text-classification-scripts

max_forks_repo_name: sadit/text-classification-scripts


text:
using Base: Float64
include("include.jl")

PATHCASE = "example"

options = zefiroOptions(PATHCASE=PATHCASE)
sizes = zefiroSizes()
data = zefiroData()
output = zefiroOutput()

load!(options,sizes,data,output)

evaluateEnergy!(options,sizes,data,output)

PCCost!(options,sizes,data,output)
PACost!(options,sizes,data,output)
ICCost!(options,sizes,data,output)

CAPEX!(options,sizes,data,output)

OMCost!(options,sizes,data,output)
DDCost!(options,sizes,data,output)

postprocess!(options,sizes,data,output)
output!(options,sizes,data,output)

max_stars_repo_name: Thiago-NovaesB/zefiro

max_issues_repo_name: Thiago-NovaesB/zefiro

max_forks_repo_name: Thiago-NovaesB/zefiro


text:
using Test, ScientificTypes, MLJScientificTypes, Random
using Tables, CategoricalArrays, CSV, DataFrames, ColorTypes, PersistenceDiagramsBase
using Dates

const Arr  = AbstractArray
const CArr = CategoricalArray
const Cat  = CategoricalValue
const Vec  = AbstractVector

include("type_tests.jl")
include("basic_tests.jl")
include("extra_tests.jl")
include("autotype.jl")
include("extra_coercion.jl")

include("coverage.jl")


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

import FirstOrderLp

include("primal_dual_hybrid_gradient.jl")
include("utils.jl")

# This script runs a experiment on the itest2 instance from the netlib database.
# It assumes the file is located in data/itest2.mps, see the README for details
# about how to obtain this file.
output_csv_path = "/tmp/blinear_current_stats.csv"
output_plots_dir = "/tmp"
print("Creating LP instance\n")
""" Returns a 2D bilinear problem.
"""
lambda_min = 1e-1
simple_lp = FirstOrderLp.linear_programming_problem(
  [-Inf, -Inf],  # variable_lower_bound
  [Inf, Inf],  # variable_upper_bound
  [0.0, 0.0],  # objective_vector
  0.0,                 # objective_constant
  [
    1.0 0.0
    0.0 lambda_min
  ],           # constraint_matrix
  [0.0, 0.0],      # right_hand_side
  2,                     # num_equalities
)

params = PdhgParameters(
  0.8, # step_size (forces the solver to use a provably correct step size)
  1.0, # primal_weight
  100, # printing frequency
  true, # verbose
  500, # iteration limit
  TwoDimensionalSubspace([0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0]),
  0.0, # eps_primal_infeasible
  0.0, # eps_dual_infeasible
  false, # store_iterates
  nothing, # reference_point
)
print("About to start optimizing\n")
solver_output = optimize(params, simple_lp, [1.0, 1.0], [1.0, 1.0])
save_stats_to_csv(solver_output, output_csv_path)
generate_plots_from_stats_csv(output_csv_path, output_plots_dir)


text:
using Enzyme
using Test

@testset "Threads $(Threads.nthreads())" begin
    function tasktest(M, x)
        xr = Ref(x)
        task = Threads.@spawn begin
            @inbounds M[1] = xr[]
        end
        @inbounds M[2] = x
        wait(task)
        nothing
    end

    R = Float64[0., 0.]
    dR = Float64[2., 3.]

    @test 5.0 ≈ Enzyme.autodiff(Reverse, tasktest, Duplicated(R, dR), Active(2.0))[1]
    @test Float64[2.0, 2.0] ≈ R
    @test Float64[0.0, 0.0] ≈ dR

    Enzyme.autodiff(Forward, tasktest, Duplicated(R, dR), Duplicated(2.0, 1.0))
    @test Float64[1.0, 1.0] ≈ dR

    function tasktest2(M, x)
        task = Threads.@spawn begin
           return
        end
        Base.wait(task)
        nothing
    end
    # The empty return previously resulted in an illegal instruction error
    @test 0.0 ≈ Enzyme.autodiff(Reverse, tasktest2, Duplicated(R, dR), Active(2.0))[1]
    @test () === Enzyme.autodiff(Forward, tasktest, Duplicated(R, dR), Duplicated(2.0, 1.0))
end

@testset "Advanced Threads $(Threads.nthreads())" begin
    function foo(y)
        Threads.@threads for i in 1:3
            y[i] *= 2
        end
        nothing
    end

    x = [1.0, 2.0, 3.0]
    dx = [1.0, 1.0, 1.0]
    Enzyme.autodiff(Reverse, foo, Duplicated(x, dx))
    @test 2.0 ≈ x[1]
    @test 4.0 ≈ x[2]
    @test 6.0 ≈ x[3]
    @test 2.0 ≈ dx[1]
    @test 2.0 ≈ dx[2]
    @test 2.0 ≈ dx[3]

    x = [1.0, 2.0, 3.0]
    dx = [1.0, 1.0, 1.0]
    @test_broken Enzyme.autodiff(Forward, foo, Duplicated(x, dx))
end

@testset "Closure-less threads $(Threads.nthreads())" begin
    function bf(i, x)
      x[i] *= x[i]
      nothing
    end

    function psquare0(x)
      Enzyme.pmap(10, bf, x)
    end

    xs = Float64[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    dxs = ones(10)

    Enzyme.autodiff(Reverse, psquare0, Duplicated(xs, dxs))
    @test Float64[2, 4, 6, 8, 10, 12, 14, 16, 18, 20] ≈ dxs

    function psquare1(x)
      Enzyme.@parallel x for i = 1:10
        @inbounds x[i] *= x[i]
      end
    end

    xs = Float64[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    dxs = ones(10)

    Enzyme.autodiff(Reverse, psquare1, Duplicated(xs, dxs))
    @test Float64[2, 4, 6, 8, 10, 12, 14, 16, 18, 20] ≈ dxs

    function psquare2(x, y)
      Enzyme.@parallel x y for i = 1:10
        @inbounds x[i] *= y[i]
      end
    end

    xs = Float64[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    dxs = ones(10)

    Enzyme.autodiff(Reverse, psquare2, Duplicated(xs, dxs), Duplicated(xs, dxs))
    @test Float64[2, 4, 6, 8, 10, 12, 14, 16, 18, 20] ≈ dxs
end

# TODO on 1.8 having `Inactive threads` after `UndefVar` in the main `runtest.jl` leads to a GC verification bug
@testset "Inactive threads" begin
    function thr_inactive(x, y)
        if x
            Threads.@threads for N in 1:5:20
                println("The number of this iteration is $N")
            end
        end
        y
    end
    @test 1.0 ≈ autodiff(Reverse, thr_inactive, false, Active(2.14))[1]
    @test 1.0 ≈ autodiff(Forward, thr_inactive, false, Duplicated(2.14, 1.0))[1]
end