text:
"This is a channel logger, used as the default callback function"
function log_by_channel(ios, devconf, accumulated, dsprev, ds)
    #dic = Dict{Int, String}(value => key for (key, value) in devconf.specbynamedict)
    dic = devconf.specbynamedict
    chna = collect(keys(dic))
    if length(chna)!= 0
        # Channel names sorted by startbit
        startbits = map(chna) do na
            addressfields = dic[na]
            addressfields.bytepos + addressfields.bitpad
        end
        p = sortperm(startbits)
        sortna = chna[p]

        # Current channel values in sequence of sortna
        v = map(sortna) do s
            spec = dic[s]
            extract(ds.values, spec...)
        end
        # Previous channel values in sequence of sortna
        vprev = if length(dsprev.values) > 0
            map(sortna) do s
                spec = dic[s]
                extract(dsprev.values, spec...)
            end
        else
            v
        end
        # Change per channel
        Δv = v - vprev

        # Format values, blank for unchanged values
        sv = map((x , Δ)-> x * Δ==0 ? "" : string(x) , v, Δv)
        str = join(lpad.(sv, 7))

        # Plain headers for files
        strh = join(lpad.(sortna, 7))
        println(ios, strh * "\n" * str)

        # Colorful headers for repl, indicating changes
        iob = IOBuffer()
        ioc = IOContext(iob, :color => true)

        for (na, Δ) in zip(sortna, Δv)
            s = lpad(na, 7)
            if Δ == 0
                printstyled(ioc, s, color= :blue)
            else
                print(ioc, s)
            end
        end
        print(ioc, "\n")
        println(ioc, str)
        println(stderr, String(take!(iob)))
    else
        println(ios, "log_by_channel: No device configuration available for \n\t$(devconf.filename)")
        throw(ErrorException("log_by_channel: No device configuration available for \n\t$(devconf.filename)"))
    end
end

"""
This is a bit logger. It is useable for interpreting the ABI
or 'Array fields', endinanness. In deed, it is easier to inspect
this than to implement a general USB feature report.
Compare the log_by_bit output with the feature reports in
'.julia_hid/DocumentationFo/'
"""
function log_by_bit(ios, devconf, accumulated, dsprev, ds)
    stringbits = string(ds.values)
    stringbitsprev = if length(dsprev.values) != length(ds.values)
        stringbits
    else
        string(dsprev.values)
    end
    stringbitsdiff = map(stringbits, stringbitsprev) do c, cprev
        if c == '\n'
            '\n'
        else
            c == cprev ? " " : "_"
        end
    end |> join
    # align the difference indication to what 'show' outputs
    sbitsdiff = "\t" * split(stringbitsdiff, '\n')[2][2:end]
    println(ios, sbitsdiff)
    println(ios, ds)
    println(stderr, sbitsdiff)
    println(stderr, ds.values)
end

"""
This is a bit logger which accumulates changes in single bits over time.
The number of times each bit has been set is shown.
Otherwise, works as log_by_bit.
"""
function log_by_bit_accumulate_changes(ios, devconf, accumulated, dsprev, ds)
    stringbits = string(ds.values)
    stringbitsprev = if length(dsprev.values) != length(ds.values)
        stringbits
    else
        string(dsprev.values)
    end
    stringbitsdiff = map(stringbits, stringbitsprev) do c, cprev
        if c == '\n'
            '\n'
        else
            c == cprev ? " " : "_"
        end
    end |> join
    # align the difference indication to what 'show' outputs
    sbitsdiff = "\t" * split(stringbitsdiff, '\n')[2][2:end]
    # accumulate changes in a vector, counting changes
    changevec = map(collect(sbitsdiff[2:end])) do c
        if c == '_'
            1
        else
            0
        end
    end
    accumulatedchanges = if length(changevec) != length(accumulated)
        changevec
    else
        changevec + accumulated
    end

    println(ios, sbitsdiff)
    println(ios, ds)
    #print_barplot(io,
    println(stderr, sbitsdiff)
    println(stderr, ds.values)
    print(stderr, "\t")
    print_barplot(stderr, accumulatedchanges)
    accumulatedchanges
end

"Find device state change"
function minus_ds(x::T, y::T) where {T <: DevState}
    if length(x.values) == length(y.values)
        DevState(x.values.data - y.values.data, x.timestamp - y.timestamp, x.proceed || y.proceed)
    else
        @debug "Could not subtract, vecor lengths, x $(length(x.values.data)), y $(length(y.values.data))"
        x
    end
end

"Prints the vector as as a string of chars from minimum(blank) to maximum (full)"
function print_barplot(io, v::Vector{T}) where {T<:Real}
    # 1..12
    bars = [' ', '_','▁', '▂', '▃', '▄', '▅', '▆', '▇', '█', '▔']
    mx = maximum(v)
    @assert minimum(v) >= 0 "horplot_string: Can't handle negative numbers"
    for i = 1:length(v)
        val = mx > 0 ? v[i] / mx : v[i]
        bucket = Int((ceil(10 * val)) + 1)
        print(io, bars[bucket])
    end
    print(io, '\n')
end


text:
module CGA

"""
Generates a binary array of values using a probability vector.
Each single element of the probability vector is the probability of bit having
the value of 1. When the probability vector is [1, 1, 1, ..., 1] then the sampled
vector is [1.0, 1.0, 1.0, ..., 1.0] whereas it is [0.0, 0.0, 0.0, ..., 0.0] when the probability vector
is a vector of zeros. The CGA (compact genetic algorithms) search is started using the
probability vector of [0.5, 0.5, 0.5, ..., 0.5] which produces random vectors of either
zeros or ones.

# Examples
```jldoctest
julia> sample([1, 1, 1, 1, 1])
5-element Array{Bool,1}:
 1
 1
 1
 1
 1
julia> cgasample(ones(10) * 0.5)
10-element Array{Bool,1}:
 1
 1
 1
 1
 0
 0
 0
 1
 1
 0
```
"""
function cgasample(probvector::Array{Float64,1})::Array{Int}
    len = length(probvector)
    ch = Array{Bool}(undef, len)
    rands = rand(Float64, len)
    @inbounds for i in eachindex(ch)
        ch[i] = Int(rands[i] < probvector[i])
    end
    return ch
end

"""
Performs a CGA (Compact Genetic Algorithm) search for minimization of an objective function.
In the example below, the objective function is to minimize sum of bits of a binary vector.
The search method results the optimum vector of [0, 0, ..., 0] where the objective function is zero.

# Examples
```jldoctest
julia> function f(x)
           return sum(x)
       end
f (generic function with 1 method)
julia> cga(chsize = 10, costfunction = f, popsize = 100)
10-element Array{Bool,1}:
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0
```
"""
function cga(; chsize::Int, costfunction::Function, popsize::Int)::Array{Bool,1}
    probvector = ones(Float64, chsize) * 0.5
    mutation = 1.0 / convert(Float64, popsize)
    while !(all(x -> (x <= mutation) || (x >= 1.0 - mutation), probvector))
        ch1 = cgasample(probvector)
        ch2 = cgasample(probvector)
        cost1 = costfunction(ch1)
        cost2 = costfunction(ch2)
        winner = ch1
        loser = ch2
        if (cost2 < cost1)
            winner = ch2
            loser = ch1
        end
        for i = 1:chsize
            if winner[i] != loser[i]
                if isone(winner[i])
                    probvector[i] += mutation
                else
                    probvector[i] -= mutation
                end
            end
        end
    end
    return cgasample(probvector)
end

end # End of module

ext: jl


text:
module adapters_mysql_structured_test

using Test # @test
using Octo.Adapters.MySQL # Schema from SELECT DISTINCT FROM WHERE GROUP BY HAVING ORDER DESC COUNT

struct User
end

u = from(User)
@test_throws Schema.TableNameError to_sql([FROM u])

Schema.model(User, table_name="users")

@test to_sql([FROM u]) == "FROM users"
@test to_sql([SELECT * FROM u]) == "SELECT * FROM users"
@test to_sql([SELECT u.id FROM u]) == "SELECT id FROM users"
@test to_sql([SELECT (u.id, u.name) FROM u]) == "SELECT id, name FROM users"
@test to_sql([WHERE u.id == 2]) == "WHERE id = 2"
@test to_sql([WHERE u.id == 2 AND u.name == "John"]) == "WHERE id = 2 AND name = 'John'"

@test to_sql([SELECT COUNT(*) FROM u]) == "SELECT COUNT(*) FROM users"
@test to_sql([SELECT COUNT(u.name) FROM u]) == "SELECT COUNT(name) FROM users"

@test to_sql([SELECT DISTINCT u.name FROM u]) == "SELECT DISTINCT name FROM users"

u = from(User, :u)
@test to_sql([SELECT COUNT(u.name) FROM u          GROUP BY u.name HAVING COUNT(u.name) > 5 ORDER BY COUNT(u.name) DESC]) ==
             "SELECT COUNT(u.name) FROM users AS u GROUP BY u.name HAVING COUNT(u.name) > 5 ORDER BY COUNT(u.name) DESC"

end # module adapters_mysql_structured_test


text:
using ExcelReaders
using Mimi: EmpiricalDistribution

function load_vector(path, range, header=false)
    tups = collect(load(path, range, header=header))
    name = keys(tups[1])[1]   # name of first item in NamedTuple
    map(obj -> getfield(obj, name), tups)
end

function load_empirical_dist(excel_file::AbstractString,
                             values_range::AbstractString,
                             probs_range::AbstractString="")
    values = load_vector(excel_file, values_range)
    probs = probs_range == "" ? nothing : load_vector(excel_file, probs_range)
    return EmpiricalDistribution(values, probs)
end


text:
## show an svg file
## reload(Pkg.dir("PySide", "examples", "svg.jl"))

using PySide
using Gadfly, Compose

w = Widget()
lyt = VBoxLayout(w)
setLayout(w, lyt)

sv = QtSvg[:QSvgWidget](w.w)
addWidget(lyt, sv)

sl = Slider(w, "Horizontal", 2:10)
push!(lyt, sl)

function draw_plot(n)
    p = plot(sin, 0, n * pi)
    nm = tempname() * ".svg"
    draw(SVG(nm, 5inch, 5inch), p)
    sv[:load](nm)
end

draw_plot(1)

qconnect(sl, :valueChanged, draw_plot) ## value changed passes in value to slot

raise(w)


text:
module TestTypes

using Base.Test
using RLEVectors

@testset begin

# Initialization
expected_run_values = [1,2,3]
expected_run_ends = [2,4,7]
x = RLEVector([1,1,2,2,3,3,3])
y = RLEVector(expected_run_values, expected_run_ends)
@test isa(x,RLEVector)
@test x.runvalues == expected_run_values
@test x.runends == expected_run_ends
x = RLEVector(["a","b","c","d","d","d","e","f"],[3,6,9,12,15,18,21,24])
y = RLEVector(["a","b","c","d","e","f"],[3,6,9,18,21,24])
@test x.runvalues == y.runvalues
@test x.runends == y.runends
x = RLEVector(BitVector(5),collect(1:5))
y = RLEVector([false],[5])
@test x.runvalues == y.runvalues
@test x.runends == y.runends
x = RLEVector([4,5,6],[2,4,6])
y = RLEVector([4,4,5,5,6,6])
@test x.runvalues == y.runvalues
@test x.runends == y.runends

@test_throws ArgumentError RLEVector([1,2,3],[2,4,3]) # Sorted runends
@test IntegerRle([1,2,3],[2,4,6]) == RLEVector([1,2,3],Int32[2,4,6])
@test FloatRle([1,2,3],[2,4,6]) == RLEVector(Float64[1,2,3],Int32[2,4,6])
@test BoolRle([true,false,true],[2,4,6]) == RLEVector([true,false,true],Int32[2,4,6])
@test StringRle(["a","bob","joe"],[2,4,6]) == RLEVector(["a","bob","joe"],Int32[2,4,6])

# Creating
y = RLEVector([1.0,1,2,2,3,3,3])
@test similar(y) == RLEVector([0.0],[length(y)])
@test similar(y,4) == RLEVector(zeros(Real,1),Int[4])

# Conversion
x = RLEVector([4,4,5,5,6,7,8])
@test convert(Vector,x) == [4,4,5,5,6,7,8]
@test convert(Set,x) == Set([4,5,6,7,8])

# Expanding
@test collect(RLEVector([1,1,2,2,3,3])) == [1,1,2,2,3,3]
@test convert(Vector, RLEVector([1,1,2,2,3,3])) == [1,1,2,2,3,3]

## Describing
x = RLEVector(expected_run_values,expected_run_ends)
@test eltype(x) == eltype(expected_run_values)
#@test typeof(show(x)) == Void # At least test that show does not give error

# Getters and setters
x = RLEVector([1,2,3],[2,9,22])
@test rfirst(x) == [1,3,10]
@test rfirst(x,1) == 1
@test rfirst(x,3) == 10
@test rwidth(x) == [2,7,13]
@test rwidth(x,1) == 2
@test rwidth(x,2) == 7
@test rwidth(x,3) == 13
@test rlast(x) == [2,9,22]
@test rvalue(x) == [1,2,3]

# Hashing
x = RLEVector([1,1,2,2,3,3,3])
d = Dict("bob" => x)
typeof(d["bob"]) == typeof(x)

end # testset

end # module


text:
include("rocketcoll.jl")

# need to import show to extend. Without import we are shadowing show
import Base: show

const tab = "   "

function show(io::IO, r::Rocket, depth::Integer = 0)
    println(io, tab^depth, "Rocket(")
    depth += 1

    show(io, r.tank, depth)
    println(io, ",")
    show(io, r.engine, depth)

    if r.payload != nopayload
        println(io, ",")
        # Engheim left out dealing with non-Rocket Payload
        # recurse if payload is rocket
        if typeof(r.payload) == Rocket
            show(io, r.payload, depth)
        else
            # for "final" payloads just print out
            println(io, tab^depth, "Payload(")
            depth += 1
            println(io, tab^depth, r.payload)
            depth -= 1
            print(io, tab^depth, ")")
        end
    end
    println(io)

    depth -= 1
    print(io, tab^depth, ")")
end

function show(io::IO, engine::Engine, depth::Integer = 0)
   # TODO: this function only handles SingleEngine not EngineCluster
   println(io, tab^depth, "Engine(")
   depth += 1
   println(io, tab^depth,
               "thrust = ",
               thrust(engine), ",")
   println(io, tab^depth,
               "Isp    = ",
               Isp(engine), ",")
   println(io, tab^depth,
               "mass   = ",
               mass(engine))
   depth -= 1
   print(io, tab^depth, ")")
end

function show(io::IO, tank::Tank, depth::Integer = 0)
   println(io, tab^depth, "Tank(")
   depth += 1
   println(io, tab^depth,
               "dry_mass   = ",
               tank.dry_mass, ",")
   println(io, tab^depth,
               "propellant = ",
               tank.propellant, ",")
   println(io, tab^depth,
               "total_mass = ",
               tank.total_mass, ",")
   depth -= 1
   print(io, tab^depth, ")")
end

function show(io::IO, ship::SpaceVehicle)
     println(io, "SpaceVehicle(")
     show(io, ship.active_stage, 1)
     println(io)
     println(io, ")")
end

# Example
# julia> show(ship)
# You don't actually need the show you can just type
# julia> ship
# SpaceVehicle(
#    Rocket(
#       Tank(
#          dry_mass   = 4000.0,
#          propellant = 107500.0,
#          total_mass = 111500.0,
#       ),
#       Engine(
#          thrust = 7.605e6,
#          Isp    = 282.0,
#          mass   = 4230.0
#       ),
#       Rocket(
#          Tank(
#             dry_mass   = 4000.0,
#             propellant = 107500.0,
#             total_mass = 111500.0,
#          ),
#          Engine(
#             thrust = 845000.0,
#             Isp    = 282.0,
#             mass   = 470.0
#          ),
#          Payload(
#             SpaceProbe(22000.0)
#          )
#       )
#    )
# )

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
using Pkg
Pkg.activate(@__DIR__)
cd(@__DIR__)
#Pkg.instantiate()

using Weave

function weave_md(filename)
    @info filename
    Weave.weave(filename; doctype = "github", fig_ext= ".png")
    @info "Done"
end

weave_md("find_valid_solutions.jl")


text:
module ThinkDSP

using FFTW
using LaTeXStrings
using LinearAlgebra
using Statistics
using RecipesBase

import Base: +, filter, *
import Latexify: latexify
import WAV: wavplay

include("Signal.jl")
include("Spectrum.jl")
include("Wave.jl")
include("Noise.jl")

include("struct_functions.jl")
include("signal_functions.jl")

export Signal, Wave, Spectrum, Integrated_Spectrum
export BrownianNoise, PinkNoise, GaussianNoise, WhiteNoise
export slope, latexify, triangle, squarewave

end


text:
macro _load_data(const_name, force, load_data)
    quote
        if !$force && !isnothing($const_name[])
            return $const_name[]
        end
        $const_name[] = $load_data
    end |> esc
end


text:
using Test
using printerpolate

data=rand(10,2)

for n in 1:10
    @test round(interpolate(data[n,1],data),digits=6) ==round(data[n,2], digits=6)
end


text:
function plotblock! end

abstract type PlotContext end
struct TextContext <: PlotContext end
struct NDContext{N} <: PlotContext end
struct NDOverlayContext{N} <: PlotContext end
struct NDOnImage{N} <: PlotContext end

struct Plot2D <: PlotContext end

"""
    plotcontext(block)
    plotcontext(blocks)

Return a [`PlotContext`](#) describing how `block` can be plotted.

Note that the data block plotting API is not stable yet, and may change
in the future.
"""
plotcontext(blocks::Tuple) = map(plotcontext, blocks)

plotcontext(::Image{N}) where N = NDContext{N}()
plotblock!(ax, ::Image{2}, data; kwargs...) = plotimage!(ax, data; kwargs...)

plotcontext(::Mask{N}) where N = NDOverlayContext{N}()
plotblock!(ax, block::Mask{2}, data; kwargs...) = plotmask!(ax, data, block.classes; kwargs...)

plotcontext(::Label) = TextContext()
function plotblock!(ax, ::Label, data)
    ax.title[] = ax.title[] * string(data)
end

plotcontext(::LabelMulti) = TextContext()
function plotblock!(ax, ::LabelMulti, data)
    ax.title[] = ax.title[] * string(data)
end

plotcontext(::Keypoints{N}) where N = NDOnImage{N}()
function plotblock!(img, ::Keypoints{2}, data)
    for k in data
        _drawkeypoint!(img, k)
    end
end

# plotsample!

function plotsample!(f, method::BlockMethod, sample)
    contexts = plotcontext(method.blocks)
    plotsample!(f, contexts, method.blocks, sample)
end

function plotsample!(f, ctxs::Tuple{NDContext{2}, TextContext}, blocks, datas)
    f[1, 1] = ax = imageaxis(f)
    plotblock!(ax, blocks[1], datas[1])  # plots an N-D image
    plotblock!(ax, blocks[2], datas[2])  # appends to axis title
end

function plotsample!(f, ctxs::Tuple{NDContext{2}, NDOverlayContext{2}}, blocks, datas)
    f[1, 1] = ax = imageaxis(f)
    plotblock!(ax, blocks[2], datas[2], alpha = 1.)  # overlays transparently over image
    plotblock!(ax, blocks[1], datas[1], alpha = 0.6)  # plots an N-D image
end

function plotsample!(
        f,
        ctxs::NTuple{N, NDContext{2}},
        blocks::NTuple{N},
        datas::NTuple{N}) where N
    for i in 1:N
        f[1, i] = ax = imageaxis(f)
        plotblock!(ax, blocks[i], datas[i])  # plots an image
    end
end

function plotsample!(
        f,
        ::Tuple{NDContext{2}, NDOnImage{2}},
        blocks,
        datas)
    # image with title showing ground truth and prediction
    ax = f[1, 1] = imageaxis(f)
    img = copy(datas[1])
    plotblock!(img, blocks[2], datas[2])
    plotblock!(ax, blocks[1], img)
end

# plotxy! decodes and falls back to `plotsample`

function plotxy!(f, method::BlockMethod, x, y)
    xyblocks = encodedblock(method.encodings, method.blocks)
    input, target = decode(method.encodings, Validation(), xyblocks, (x, y))
    plotsample!(f, method, (input, target))
end

function plotprediction!(f, method::BlockMethod, x, ŷ, y)
    xblock, yblock = encodedblock(method.encodings, method.blocks)
    ŷblock = method.outputblock
    blocks = (xblock, ŷblock, yblock)
    input, target_pred, target = decode(method.encodings, Validation(), blocks, (x, ŷ, y))
    inblocks = decodedblock(method.encodings, blocks, true)
    contexts = plotcontext(inblocks)
    plotprediction!(f, contexts, inblocks, (input, target_pred, target))
end

function plotprediction!(
        f,
        ::Tuple{NDContext{2}, TextContext, TextContext},
        blocks,
        datas)

    ax = f[1, 1] = imageaxis(f)
    plotblock!(ax, blocks[1], datas[1])
    ax.title[] = ax.title[] * "Prediction: "
    plotblock!(ax, blocks[2], datas[2])
    ax.title[] = ax.title[] * " | True: "
    plotblock!(ax, blocks[3], datas[3])
end

function plotprediction!(
        f,
        ::Tuple{NDContext{2}, NDOverlayContext{2}, NDOverlayContext{2}},
        blocks,
        datas)

    # image with title showing ground truth and prediction
    ax1 = f[1, 1] = imageaxis(f)
    ax1.title[] = "True"
    plotblock!(ax1, blocks[3], datas[3])
    plotblock!(ax1, blocks[1], datas[1], alpha=0.4)

    ax2 = f[1, 2] = imageaxis(f)
    ax2.title[] = "Prediction"
    plotblock!(ax2, blocks[2], datas[2])
    plotblock!(ax2, blocks[1], datas[1], alpha=0.4)
end

function plotprediction!(
        f,
        ::Tuple{NDContext{2}, NDOnImage{2}, NDOnImage{2}},
        blocks,
        datas)

    # image with title showing ground truth and prediction
    ax1 = f[1, 1] = imageaxis(f)
    ax1.title[] = "True"
    img1 = copy(datas[1])
    plotblock!(img1, blocks[3], datas[3])
    plotblock!(ax1, blocks[1], img1)

    ax2 = f[1, 2] = imageaxis(f)
    ax2.title[] = "Prediction"
    img2 = copy(datas[1])
    plotblock!(img2, blocks[2], datas[2])
    plotblock!(ax2, blocks[1], img2)
end


text:
using WebIO
using Test

@testset "Scope id deprecations" begin
    @test_deprecated Scope("manual-ids-are-deprecated")
    @test_deprecated Scope(; id="still-deprecated")

    # Make sure that non-deprecated arguments are still handled as expected.
    # This surfaced in https://github.com/JuliaGizmos/Interact.jl/issues/315.
    @test_deprecated !isempty(Scope("foo"; imports=["foo.js"]).imports)
end

@testset "register_renderable function deprecation" begin
    MyTypeName = gensym()
    @eval struct $MyTypeName end
    @eval WebIO.render(x::$MyTypeName) = x.dom
    @test_deprecated @eval WebIO.register_renderable($MyTypeName)
    @test hasmethod(Base.show, (IO, MIME"text/html", @eval($MyTypeName)))
end


text:
using Test

@testset "BinaryTraits Tests" begin
    include("test_super_type.jl")
    include("test_syntax.jl")
    include("test_single_trait.jl")
    include("test_multiple_traits.jl")
    include("test_composite_traits.jl")
    include("test_verbose.jl")
    include("test_parametric_type.jl")
    include("test_cross_module.jl")
    include("test_interfaces.jl")
    include("test_traitfn.jl")
    VERSION >= v"1.1" && include("test_return_types.jl")
end


text:
include("helpers.jl")

# Functions here have similar structure and thus similar arguments:
# details arguments here:
# - df : data frame include all bio-tracers values and tru origins

# increase the number of bio-tracers
function scr_nbio_v2(df, n_bio, n_distr, n_sample, mxcb, n_rep, n_tot_bio,
    n_tot_indiv, n_reg, aug, sig, path, pca = false)

    if pca
        cb_cols = [1:n_bio]
    else
        cb_cols = sample_cols(n_bio, n_tot_bio, mxcb)
    end

    # loop over bio-tracers combinations
    for j in 1:size(cb_cols, 1)
        println("  -- combin = ", j, "--")
        # loop over repetition
        for k in 1:n_rep
            # println("   -> replicate #", k)
            # all samples that will be used
            spl0 = samp_samples(n_distr + n_sample, n_tot_indiv, n_reg)
            # spl1 => sample in spl0 that'll be used in the training set
            spl1 = samp_samples(n_distr, n_distr + n_sample, n_reg)
            df_tmp = df[:, cb_cols[j]]
            # add origin (last column) back
            df_tmp.origin = df.origin
            fl = @sprintf "%sml_nbio_%02d_%03d_%05d.txt" path n_bio j k
            train = df_tmp[spl0[spl1], :]
            valid = df_tmp[setdiff(spl0, spl0[spl1]), :]
            res = simu_ml(train, valid, sig, fl, aug)
        end
    end
end

# increasing the number of bio-tracers
function scr_nbio(df, n_bio_max, n_distr, n_sample, mxcb, n_rep, n_tot_bio,
    n_tot_indiv, n_reg, aug, sig, path, pca = false)
    for i in 1:n_bio_max
        println("------- n_bio = ", i, "------")
        if pca
            cb_cols = [1:n_bio_max]
        else
            cb_cols = sample_cols(i, n_tot_bio, mxcb)
        end
        for j in 1:size(cb_cols, 1)
            println("  -- combin = ", j, "--")
            for k in 1:n_rep

                # println("   -> replicate #", k)
                # all samples that will be used
                spl0 = samp_samples(n_distr + n_sample, n_tot_indiv, n_reg)
                # spl1 => sample in spl0 that'll be used in the training set
                spl1 = samp_samples(n_distr, n_distr + n_sample, n_reg)
                df_tmp = df[:, cb_cols[j]]
                # add origin (last column) back
                df_tmp.origin = df.origin
                fl = @sprintf "%sml_nbio_%02d_%03d_%03d.txt" path i j k
                train = df_tmp[spl0[spl1], :]
                valid = df_tmp[setdiff(spl0, spl0[spl1]), :]
                res = simu_ml(train, valid, sig, fl, aug)
            end
        end
    end
end

# increasing the size of the training set
function scr_ndistr(df, n_bio, n_sample, mxcb, n_rep, n_tot_bio, n_tot_indiv,
    n_reg, aug, sig, id_rep)
    # set seed
    Random.seed!(randperm(1000)[id_rep])
    # loop over training set size
    for i in 5:25
        println("------- n_distr = ", i, "------")
        cb_cols = sample_cols(n_bio, n_tot_bio, mxcb)
        for j in 1:size(cb_cols, 1)
            println("  -- combin = ", j, "--")
            # loop over repetition
            for k in 1:n_rep
                println("   -> replicate #", k)
                # all samples that will be used
                spl0 = samp_samples(i + n_sample, n_tot_indiv, n_reg)
                # spl1 => sample in spl0 that'll be used in the training set
                spl1 = samp_samples(i, i + n_sample, n_reg)
                df_tmp = df[:, cb_cols[j]]
                df_tmp.origin = df.origin
                train = df_tmp[spl0[spl1], :]
                valid = df_tmp[setdiff(spl0, spl0[spl1]), :]
                fl = @sprintf "ml_ndistr_%02d_%02d.txt" n_bio id_rep
                res = simu_ml2(train, valid, sig, fl, i, j, k, aug)
            end
        end
    end
end

# increasing the number of sample (not used anymore as we used nbio instead)
function scr_nsample(df, n_bio, n_distr, n_sample_max, mxcb, n_rep, n_tot_bio,
    n_tot_indiv, n_reg, aug, sig, path, pca = false)
    for i in 1:n_sample_max
        println("------- n_sample = ", i, "------")
        if pca
            cb_cols = [1:n_bio_max]
        else
            cb_cols = sample_cols(n_bio, n_tot_bio, mxcb)
        end
        for j in 1:size(cb_cols, 1)
            println("  -- combin = ", j, "--")
            for k in 1:n_rep
                println("   -> replicate #", k)
                # all samples that will be used
                spl0 = samp_samples(n_distr + i, n_tot_indiv, n_reg)
                # spl1 => sample in spl0 that'll be used in the training set
                spl1 = samp_samples(n_distr, n_distr + i, n_reg)
                df_tmp = df[:, cb_cols[j]]
                df_tmp.origin = df.origin
                fl = @sprintf "%sml_nsamp_%02d_%03d_%03d.txt" path i j k
                res = simu_ml(df_tmp[spl0[spl1], :],
                    df_tmp[setdiff(spl0, spl0[spl1]), :], sig, fl, aug)
            end
        end
    end
end

# increase noise level
function scr_noise(df, n_bio, n_distr, n_sample, mxcb, n_rep, n_tot_bio,
    n_tot_indiv, n_reg, aug, sig, id_rep)
    # set seed
    Random.seed!(randperm(1000)[id_rep])
    # noise sequence
    noise = [exp10(k) for k in -4:.2:1]
    for (i, s) in enumerate(noise)
        println("- noise = ", s, "---------")
        cb_cols = sample_cols(n_bio, n_tot_bio, mxcb)
        for j in 1:size(cb_cols, 1)
            println("  -- combin = ", j, "--")
            for k in 1:n_rep
                println("   -> replicate #", k)
                df_tmp0 = add_noise_df(df, s)
                # all samples that will be used
                spl0 = samp_samples(n_distr + n_sample, n_tot_indiv, n_reg)
                # spl1 => sample in spl0 that'll be used in the training set
                spl1 = samp_samples(n_distr, n_distr + n_sample, n_reg)
                df_tmp = df_tmp0[:, cb_cols[j]]
                df_tmp.origin = df.origin
                fl = @sprintf "ml_noise_%02d_%02d.txt" n_bio id_rep
                res = simu_ml2(df_tmp[spl0[spl1], :],
                    df_tmp[setdiff(spl0, spl0[spl1]), :], sig, fl, s, j, k, aug)
            end
        end
    end
end

# testing the influence if noise addition and data augmentation
function scr_train(df, col, fl)
    sig = [0 0.00001 0.0001 0.001 0.005 0.01 0.05 0.1 .2 0.5 1 2]
    aug = [1 2 5 10 20 50 100 200 500 1000 2000 5000]
    for (i, s) in enumerate(sig)
        for (j, a) in enumerate(aug)
        println("- train: i = ", i, "---- j = ", j)
        res = simu_ml2(
            df[vcat(1:20, 31:50, 61:80), col],
            df[vcat(21:30, 51:60, 81:90), col],
            s, fl, size(col, 1), s, a, a)
        end
    end
end

# df = get_data(false, "../data/data_f_cs.csv")
# scr_train(df, 1:18, "res_train_all.txt")
# scr_train(df, vcat(1:3, 18), "res_train_si.txt")
# scr_train(df, vcat(4:6, 18), "res_train_fa.txt")


text:
using Test, OMEinsum
using OMEinsum: IndexGroup, NestedEinsum, parse_nested, DynamicEinCode, isleaf, collect_ixs
@testset "einsequence" begin
    @test push!(IndexGroup([],1), 'c').inds == IndexGroup(['c'], 1).inds
    @test isempty(IndexGroup([],1))

    @test_throws ArgumentError parse_nested("((ij,jk),km")
    @test_throws ArgumentError parse_nested("((ij,jk),k1)")

    a, b, c = rand(2,2), rand(2,2), rand(2,2)
    abc1 = ein"(ij,jk),km -> im"(a,b,c)
    abc2 = ein"((ij,jk),km) -> im"(a,b,c)
    abc3 = ein"ij,jk,km -> im"(a,b,c)

    @test abc1 ≈ abc2 ≈ abc3
    size_info = Dict('k'=>2)
    a, b, c, d = randn(2), randn(2,2), randn(2), randn(2)
    @test ein"((i,ij),i),j->ik"(a, b, c, d; size_info=size_info) ≈ ein"i,ij,i,j->ik"(a, b, c, d; size_info=size_info)
    @test collect_ixs(ein"((i,ij),i),j->ik") == collect_ixs(ein"i,ij,i,j->ik") == collect_ixs(DynamicEinCode(ein"i,ij,i,j->ik")) == [['i'], ['i','j'], ['i'], ['j']]
end

@testset "macro" begin
    b, c, d = rand(2,2), rand(2,2,2), rand(2,2,2,2)
    @ein a[i,j] := b[i,k] * c[k,k,l] * d[l,m,m,j]
    @ein a2[i,j] := b[i,k] * (c[k,k,l] * d[l,m,m,j])
    @test a ≈ a2
end

@testset "flatten" begin
    ne = ein"(ij,j),k->"
    @test OMEinsum.flatten(ne) === ein"ij,j,k->"
    todynamic(ne::NestedEinsum) = isleaf(ne) ? NestedEinsum{DynamicEinCode{Char}}(ne.tensorindex) : NestedEinsum(todynamic.(ne.args), DynamicEinCode(ne.eins))
    ne2 = todynamic(ne)
    @test OMEinsum.flatten(ne2) isa DynamicEinCode && OMEinsum.flatten(ne2) == DynamicEinCode(ein"ij,j,k->")
end


text:
# NOTE: the implementation in openCV is very optimal, I have to study that one and hopefully
# get to implement it.

connectivity4() = ( (-1,0),(0,1),(1,0),(0,-1) )
connectivity8() = ( (-1,0,0),(1,0,0),(0,-1,0),(0,1,0),(0,0,-1),(0,0,1) )
default_connect(N) = ( N == 2 ) ? connectivity4() : connectivity8();

""" one at a time connected component labelling in 2D """
function OAATccl( mask::Array{UInt8,2}; connectivity=connectivity4(), ref=UInt8(1) )

    h, w  = size(mask);
    mask_ = copy(mask);

    components = Array{Array{Tuple{Int64,Int64},1},1}(undef,0);
    comp_idx   =  0;
    stack      = [];

    for col in 1:w, row in 1:h
        # starting the process
        if mask_[ row, col ] == ref
            comp_idx += 1;
            push!( components, Array{Tuple{Int64,Int64},1}(undef,0) )
            push!( stack, (row,col) );
        end

        # Iterating
        while length(stack) > 0
            idx = pop!( stack )
            mask_[ idx... ] = UInt8(!Bool(ref))
            push!( components[comp_idx], idx );

            for off in connectivity
                r, c = min.((h,w), max.(1, idx .+ off))
                if mask_[r,c] == ref
                    push!( stack, ( r, c ) );
                end
            end
        end
    end
    return components;
end

function OAATcclrev( mask::Array{UInt8,2}; connectivity=connectivity4() )
	return OAATccl( mask, connectivity=connectivity, ref=UInt8(0) )
end

""" Boolean image, mask """
function OAATccl( mask::Array{Bool,2}; connectivity=connectivity4() )

    h, w  = size(mask);
    mask_ = copy(mask);

    components = Array{Array{Tuple{Int64,Int64},1},1}(undef,0);
    comp_idx   =  0;
    stack      = [];

    for col in 1:w, row in 1:h

        # starting the process
        if mask_[ row, col ]
            comp_idx += 1;
            push!( components, Array{Tuple{Int64,Int64},1}(undef,0) )
            push!( stack, (row,col) );
            mask_[ row, col ] = false
        end

        # Iterating
        while length(stack) > 0
            idx = pop!( stack )
            push!( components[comp_idx], idx );

            for off in connectivity
                r, c = min.((h,w), max.(1, idx .+ off))
                if mask_[r,c]
                    push!( stack, ( r, c ) );
                    mask_[r,c] = false
                end
            end
        end
    end
    return components;
end

function OAATccl_rev( mask::Array{Bool,2}; connectivity=connectivity4() )

    h, w  = size(mask);
    mask_ = copy(mask);

    components = Array{Array{Tuple{Int64,Int64},1},1}(undef,0);
    comp_idx   =  0;
    stack      = [];

    for col in 1:w, row in 1:h

        # starting the process
        if !mask_[ row, col ]
            comp_idx += 1;
            push!( components, Array{Tuple{Int64,Int64},1}(undef,0) )
            push!( stack, (row,col) );
            mask_[ row, col ] = true
        end

        # Iterating
        while length(stack) > 0
            idx = pop!( stack )
            push!( components[comp_idx], idx );

            for off in connectivity
                r, c = min.((h,w), max.(1, idx .+ off))
                if !mask_[r,c]
                    push!( stack, ( r, c ) );
                    mask_[r,c] = true
                end
            end
        end
    end
    return components;
end

function getBorders_Centroid( ccl, mask )
    borders  = Array{typeof(ccl[1]),1}( undef,0 );
    centroid = zeros( Float32, length( ccl[1] ) );
    offs     = default_connect( length(ccl[1] ) );

    for coords in ccl
        centroid .+= coords;

        if any( coords .== 1 ) || any( coords .== size(mask) )
            push!( borders, coords )
        else
            for off in offs
                Bool( mask[ (coords .+ off)... ] ) || ( push!( borders, coords ); break; )
            end
        end

    end
    return borders, centroid ./ length(ccl)
end

# Needs borders and Centroid
function circularities( ccls, mask )

    circularities = Array{Float32,1}(undef,length(ccls));
    for idx in 1:length(ccls)
        ccl = ccls[idx]
        bords, centroid = getBorders_Centroid( ccl, mask )
        rad = 0.0
        for coords in bords
            rad += sum( abs.( coords .- centroid ) );
        end
        rad = rad / length(bords)
        circ = 0.0
        for coords in bords
            circ +=  sum( abs.( abs.( coords .- centroid ) .- rad ) );
        end
        circularities[idx] = circ/length(bords)
    end

    return circularities;
end

# Needs borders and Centroid
function circularity( ccl, mask )

    bords, centroid = getBorders_Centroid( ccl, mask )
    rad = 0.0
    for coords in bords
        rad += sum( abs.( coords .- centroid ) );
    end
    rad = rad / length(bords)
    circ = 0.0
    for coords in bords
        circ +=  sum( abs.( abs.( coords .- centroid ) .- rad ) );
    end

    return circ/length(bords)
end

""" one at a time connected component labelling in OAAT 3D """

function OAATccl( mask::Array{Bool,3}; ref=true, connectivity=connectivity8() )

    h, w, d = size(mask);
    mask_   = copy(mask);

    components = Array{Array{Tuple{Int64,Int64,Int64},1},1}(undef,0);
    comp_idx   =  0;
    stack      = [];

	for zet in 1:d, col in 1:w, row in 1:h

        if mask_[ row, col, zet ] == ref
            comp_idx += 1;
            push!( components, Array{Tuple{Int64,Int64,Int64},1}(undef,0) )
            push!( stack, (row,col,zet) );
            mask_[ row,col,zet ] = !ref
        end

        # Iterating
        while length(stack) > 0
            idx = pop!( stack )
            push!( components[comp_idx], idx );

            for off in connectivity
                r, c, z = min.( (h,w,d), max.( 1, idx .+ off ) )
                if mask_[r,c,z] == ref
                    push!( stack, ( r, c, z ) );
                    mask_[ r,c,z ] = !ref
                end
            end
        end

    end

    return components;
end

function OAATccl( mask::Array{UInt8,3}; ref=UInt8(1), connectivity=connectivity8() )

    h, w, d = size(mask);
    mask_   = copy(mask);

    components = Array{Array{Tuple{Int64,Int64,Int64},1},1}(undef,0);
    comp_idx   =  0;
    stack      = [];

	for zet in 1:d, col in 1:w, row in 1:h

        if mask_[ row, col, zet ] == ref
            comp_idx += 1;
            push!( components, Array{Tuple{Int64,Int64,Int64},1}(undef,0) )
            push!( stack, (row,col,zet) );
        end

        # Iterating
        while length(stack) > 0
            idx = pop!( stack )
            mask_[ idx... ] = UInt8(!Bool(ref))
            push!( components[comp_idx], idx );

            for off in connectivity
                r, c, z = min.( (h,w,d), max.( 1, idx .+ off ) )
                if mask_[r,c,z] == ref
                    push!( stack, ( r, c, z ) );
                end
            end
        end

    end

    return components;
end

function OAATcclrev( mask::Array{UInt8,3}; connectivity=connectivity8() )

	return OAATccl( mask, connectivity=connectivity, ref=UInt8(0) )
end

""" Utility functions """

function CCLlimits( ccl )
	miny = typemax(eltype(ccl[1]))
	minx = typemax(eltype(ccl[1]))
	maxy = 0
	maxx = 0
	@simd for coords in ccl
		miny = ( coords[1] < miny ) ? coords[1] : miny
		minx = ( coords[2] < minx ) ? coords[2] : minx
		maxy = ( coords[1] > maxy ) ? coords[1] : maxy
		maxx = ( coords[2] > maxx ) ? coords[2] : maxx
	end
	return miny, minx, maxy, maxx
end

function CCLcentroid( ccl::Array{Tuple{Int64,Int64,Int64},1} )
	cy, cx, cz = zeros( Float32, 3 )
	for coords in ccl
		cy += coords[1]
		cx += coords[2]
		cz += coords[3]
	end
	return (cy,cx,cz)./length(ccl)
end

function CCLcentroid( ccl::Array{Tuple{Int64,Int64},1} )
	cy, cx = zeros( Float32, 2 )
	for coords in ccl
		cy += coords[1]
		cx += coords[2]
	end
	return (cy,cx)./length(ccl)
end

function infoAboutPositiveCCLS( input::Array{<:Real,N}; areaTH=0 ) where {N}

	mean  = 0.0
	@inbounds @simd for idx in 1:length(input)
		mean += input[idx]
	end
	mean  = mean/length(input)
	signs = Array{UInt8,N}(undef,size(input))
	@inbounds @simd for idx in 1:length(input)
		signs[idx] = UInt8( sign(input[idx]-mean) > 0 )
	end

	CCLS    = OAATccl( signs )
    cont    = 0
	maxArea = 0
    for ccl in CCLS
        cont   += length(ccl) > areaTH
		maxArea = ( length(ccl) ) > maxArea ? length(ccl) : maxArea
    end
	return cont, maxArea
end

function countNegativeCCLS( input::Array{<:Real,N}; areaTH=0 ) where {N}

	mean  = 0.0
	@inbounds @simd for idx in 1:length(input)
		mean += input[idx]
	end
	mean = mean/length(input)
	signs = Array{UInt8,N}(undef,size(input))
	@inbounds @simd for idx in 1:length(input)
		signs[idx] = UInt8( sign(input[idx]-mean) < 0 )
	end

	CCLS = OAATccl( signs )
    cont = 0
    for ccl in CCLS
        cont += length(ccl) > areaTH
    end
	return cont
end

# Extra
function OAATcclID( mask::Array{UInt16,N}, ID; connectivity=default_connect(N) ) where {N}
	isID = Array{UInt8,N}(undef,size(mask))
	@inbounds @simd for idx in 1:length(mask)
		isID[idx] == ( mask[idx] == ID )
	end
	return OAATccl( mask, connectivity=connectivity )
end

function OAATccl( mask::Array{UInt16,3};
                  connectivity=( (-1,0,0), (1,0,0),
                                 (0,-1,0), (0,1,0),
                                 (0,0,-1), (0,0,1) ) )

    h, w, d = size(mask);
    mask_   = copy(mask);
	num     = length( unique( mask[:] ) );

    components = [ Array{Tuple{Int64,Int64,Int64},1}(undef,0) for x in 1:num-1 ]

	for zet in 1:d, col in 1:w, row in 1:h

        if mask_[ row, col, zet ] > 0

            push!( components[ mask_[row,col,zet] ], (row,col,zet) )
			mask_[row,col,zet] = 0
        end
	end

    return components;
end

max_stars_repo_name: BrucePBon/ImageAnalysis

max_issues_repo_name: BrucePBon/ImageAnalysis

max_forks_repo_name: BrucePBon/ImageAnalysis


text:
tsname="joSincInterp"
@testset "$tsname" begin
    xin = range(0.0, 10.0, length=201)
    xout = range(0.0, 10.0, length=401)
    # Difference in dx
    scale = 2
    # Input coarser
    Is = joSincInterp(xin, xout)
    @test isadjoint(Is)[1]
    v = sin.(pi .* xin)
    vi = Is*v
    # Makes sure we do not have unwanted zeros
    @test vi[2] != 0
    # Verifies some known values
    @test isapprox(maximum(vi), maximum(v))
    @test isapprox(v[11], vi[21])

    # Output coarser
    Is = joSincInterp(xout, xin)
    @test isadjoint(Is)[1]
    v = sin.(pi .* xout)
    vi = Is*v
    # Verifies some known values
    @test isapprox(maximum(vi), maximum(v))
    @test isapprox(v[21], vi[11])
end


text:
using Test
using PolyPairPots, JuLIP, LinearAlgebra, Test
using JuLIP.Testing, JuLIP.MLIPs

randr() = 1.0 + rand()
randcoeffs(B) = rand(length(B)) .* (1:length(B)).^(-2)

##

@testset "PolyPairPots.jl" begin

include("test_basics.jl")
include("test_repulsion.jl")

end


text:
module FileTrees

using FilePathsBase
using AbstractTrees
import AbstractTrees: children
import Base: parent, getindex

export FileTree, File, name, Path, path, maketree, children, rename, setvalue, setparent

import FilePathsBase: /, Path, @p_str
export @p_str, get, values, nodes, files, dirs

include("datastructure.jl")

import Base: cp, mv, rm, touch, mkpath

include("tree-ops.jl")

import Glob: GlobMatch, @glob_str

export @glob_str
export mapsubtrees

include("patterns.jl")

export mapvalues, reducevalues, NoValue, hasvalue # load, save

include("values.jl")

import Dagger
import Dagger: compute, delayed, Chunk, Thunk

export lazy, exec, compute

include("parallelism.jl")

end # module


text:
using CImGui.CSyntax
using CImGui.CSyntax.CStatic

#TODO what to do with contexts, do we actually need to save them? Should they all be in the canvas?
@component @with_kw mutable struct GuiInfo
    show_gui::Bool = false
end

@component struct GuiText
    text::String
end

@component @with_kw struct GuiFuncs
    funcs::Vector{Function} = Function[]
end

struct GuiRenderer <: AbstractRenderSystem end

Overseer.requested_components(::GuiRenderer) = (GuiFuncs, GuiText)

function Overseer.prepare(::GuiRenderer, dio::Diorama)
    return isempty(dio[GuiFuncs]) &&
           (dio[Entity(1)] = GuiFuncs(); dio[Entity(1)] = GuiInfo())
end

function Overseer.update(::GuiRenderer, m::AbstractLedger)
    camera = singleton(m, Camera3D)
    gui_info = singleton(m, GuiInfo)
    keyboard = singleton(m, Keyboard)

    if pressed(keyboard) && keyboard.button == GLFW.KEY_G
        gui_info.show_gui = keyboard.modifiers != 1
    end

    if gui_info.show_gui
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        Gui.NewFrame()
        # if Gui.Checkbox("Show GUI")
        if !isempty(m[GuiText])
            Gui.Begin("User Text")
            camera.locked = Gui.IsWindowFocused() || camera.locked
            Gui.SetWindowFontScale(2.0f0)
            for e in m[GuiText]
                Gui.Text(e.text)
            end
            Gui.End()
        end

        # Submitted Gui Funcs
        fs = singleton(m, GuiFuncs).funcs
        for (i, f) in enumerate(fs)
            f()
        end
        empty!(fs)
        camera.locked = Gui.IsAnyItemActive() || camera.locked
        Gui.Render()
        ImGui_ImplOpenGL3_RenderDrawData(Gui.GetDrawData())
    end

    # Components Debug
    # Gui.Begin("Components")
    #    Gui.SetWindowFontScale(2.0f0)
    # for c in filter(x -> !isempty(x), renderer.data.components)
    #     if Gui.TreeNode(replace("$(eltype(c))", "Glimpse." => ""))
    #         eids = @valid_entities_in(c)
    #         if length(eids) < 5
    #             for e in eids
    #                 if Gui.TreeNode("$e")
    #                     Gui.Text("$(c[e])")
    #                     Gui.TreePop()
    #                 end
    #             end
    #         end
    #         Gui.TreePop()
    #     end
    # end
    # Gui.End()
    ###
end


text:
################################################################################
#                            MATH SYMBOL METHODS
################################################################################
# Support additional math symbols beyond what JuMP does
function _infopt_math_symbol(::Type{JuMP.REPLMode}, name::Symbol)::String
    if name == :intersect
        return Sys.iswindows() ? "and" : "∩"
    elseif name == :prop
        return "~"
    elseif name == :partial
        return Sys.iswindows() ? "d" : "∂"
    elseif name == :expect
        return Sys.iswindows() ? "E" : "𝔼"
    elseif name == :integral
        return Sys.iswindows() ? "integral" : "∫"
    else
        return JuMP._math_symbol(JuMP.REPLMode, name)
    end
end

# Support additional math symbols beyond what JuMP does
function _infopt_math_symbol(::Type{JuMP.IJuliaMode}, name::Symbol)::String
    if name == :intersect
        return "\\cap"
    elseif name == :prop
        return "\\sim"
    elseif name == :partial
        return "\\partial"
    elseif name == :open_rng
        return "\\left["
    elseif name == :close_rng
        return "\\right]"
    elseif name == :expect
        return "\\mathbb{E}"
    elseif name == :integral
        return "\\int"
    else
        return JuMP._math_symbol(JuMP.IJuliaMode, name)
    end
end

################################################################################
#                            INFINITE SET METHODS
################################################################################
# Return "s" if n is greater than one
_plural(n) = (isone(n) ? "" : "s")

## Return the string of an infinite set
# IntervalSet
function set_string(print_mode, set::IntervalSet)::String
    return string("[", JuMP._string_round(JuMP.lower_bound(set)), ", ",
                  JuMP._string_round(JuMP.upper_bound(set)), "]")
end

# DistributionSet
function set_string(print_mode, set::DistributionSet)::String
    return string(set.distribution)
end

# CollectionSet
function set_string(print_mode, set::CollectionSet)::String
    sets = collection_sets(set)
    num_sets = length(sets)
    cs_string = string("CollectionSet with ", num_sets, " set",
                       _plural(num_sets), ":")
    for s in sets
        cs_string *= string("\n ", set_string(print_mode, s))
    end
    return cs_string
end

# Fallback
function set_string(print_mode, set::S) where {S <: AbstractInfiniteSet}
    # hack fix because just calling `string(set)` --> StackOverflow Error
    values = map(f -> getfield(set, f), fieldnames(S))
    return string(S, "(", join(values, ", "), ")")
end

## Return in set strings of infinite sets
# Extend to return of in set string for interval sets
function JuMP.in_set_string(print_mode, set::IntervalSet)::String
    if JuMP.lower_bound(set) != JuMP.upper_bound(set)
        return string(_infopt_math_symbol(print_mode, :in), " ",
                      set_string(print_mode, set))
    else
        return string(_infopt_math_symbol(print_mode, :eq), " ",
                      JuMP._string_round(JuMP.lower_bound(set)))
    end
end

# Extend to return of in set string for distribution sets
function JuMP.in_set_string(print_mode, set::DistributionSet)::String
    dist = set.distribution
    name = string(typeof(dist))
    bracket_index = findfirst(isequal('{'), name)
    if bracket_index !== nothing
        name = name[1:bracket_index-1]
    end
    if !isempty(size(dist))
        dims = size(dist)
        name *= string("(dim", _plural(length(dims)), ": (", join(dims, ", "), "))")
    end
    return string(_infopt_math_symbol(print_mode, :prop), " ", name)
end

# Extend to return of in set string of other sets
function JuMP.in_set_string(print_mode, set::AbstractInfiniteSet)::String
    return string(_infopt_math_symbol(print_mode, :in), " ",
                  set_string(print_mode, set))
end

## Extend in set string to consider parameter bounds
# IntervalSet
function JuMP.in_set_string(print_mode,
    pref::GeneralVariableRef,
    set::IntervalSet,
    bounds::ParameterBounds{GeneralVariableRef})::String
    # determine if in bounds
    in_bounds = haskey(bounds, pref)
    # make the string
    interval = in_bounds ? bounds[pref] : set
    return JuMP.in_set_string(print_mode, interval)
end

# InfiniteScalarSet
function JuMP.in_set_string(print_mode,
    pref::GeneralVariableRef,
    set::InfiniteScalarSet,
    bounds::ParameterBounds{GeneralVariableRef})::String
    # determine if in bounds
    if haskey(bounds, pref)
        bound_set = bounds[pref]
        if JuMP.lower_bound(bound_set) == JuMP.upper_bound(bound_set)
            return JuMP.in_set_string(print_mode, bound_set)
        else
            return  string(JuMP.in_set_string(print_mode, set), " ",
                           _infopt_math_symbol(print_mode, :intersect),
                           " ", set_string(print_mode, bound_set))
        end
    else
        return JuMP.in_set_string(print_mode, set)
    end
end

################################################################################
#                           MEASURE STRING METHODS
################################################################################
## Convert measure data into a useable string for measure printing
# 1-D DiscreteMeasureData/FunctionalDiscreteMeasureData
function measure_data_string(print_mode,
    data::Union{DiscreteMeasureData{GeneralVariableRef},
                FunctionalDiscreteMeasureData{GeneralVariableRef}}
    )::String
    pref = parameter_refs(data)
    lb = JuMP.lower_bound(data)
    ub = JuMP.upper_bound(data)
    nan_bound = isnan(lb) || isnan(ub)
    if nan_bound
        return JuMP.function_string(print_mode, pref)
    else
        set = IntervalSet(lb, ub)
        return string(JuMP.function_string(print_mode, pref), " ",
                      JuMP.in_set_string(print_mode, set))
    end
end

# Multi-D DiscreteMeasureData/FunctionalDiscreteMeasureData
function measure_data_string(print_mode,
    data::Union{DiscreteMeasureData{Vector{GeneralVariableRef}},
                FunctionalDiscreteMeasureData{Vector{GeneralVariableRef}}}
    )::String
    prefs = parameter_refs(data)
    lbs = JuMP.lower_bound(data)
    ubs = JuMP.upper_bound(data)
    has_bounds = !isnan(first(lbs)) && !isnan(first(ubs))
    homo_bounds = has_bounds && _allequal(lbs) && _allequal(ubs)
    names = map(p -> _remove_name_index(p), prefs)
    homo_names = _allequal(names)
    num_prefs = length(prefs)
    if homo_names && homo_bounds
        set = IntervalSet(first(lbs), first(ubs))
        return string(first(names), " ", JuMP.in_set_string(print_mode, set),
                      "^", num_prefs)
    elseif has_bounds
        str_list = [JuMP.function_string(print_mode, prefs[i]) * " " *
                    JuMP.in_set_string(print_mode, IntervalSet(lbs[i], ubs[i]))
                    for i in eachindex(prefs)]
        return _make_str_value(str_list)[2:end-1]
    elseif homo_names
        return first(names)
    else
        return _make_str_value(prefs)
    end
end

# extract the most compact parameter name possible
function _get_root_parameter_name(data::AbstractMeasureData)::String
    prefs = parameter_refs(data)
    names = map(p -> _remove_name_index(p), prefs)
    if _allequal(names)
        return first(names)
    else
        return _make_str_value(prefs)
    end
end

# Fallback for measure_data_string
function measure_data_string(print_mode, data::AbstractMeasureData)::String
    return _get_root_parameter_name(data)
end

# Make strings to represent measures in REPLMode
function variable_string(m::Type{JuMP.REPLMode}, mref::MeasureRef)::String
    data = measure_data(mref)
    data_str = measure_data_string(m, data)
    func_str = JuMP.function_string(m, measure_function(mref))
    name = JuMP.name(mref)
    if name == "integral"
        name = _infopt_math_symbol(m, :integral)
    elseif name == "expect"
        name = _infopt_math_symbol(m, :expect)
    end
    return string(name, "{", data_str, "}[", func_str, "]")
end

# Make strings to represent measures in IJuliaMode
function variable_string(m::Type{JuMP.IJuliaMode}, mref::MeasureRef)::String
    data = measure_data(mref)
    data_str = measure_data_string(m, data)
    func_str = JuMP.function_string(m, measure_function(mref))
    name = JuMP.name(mref)
    if name == "integral"
        param_name = _get_root_parameter_name(data)
        end_str = length(param_name) == 1 ? "d$param_name" : "d($param_name)"
        return string("\\int_{", data_str, "}", func_str, end_str)
    else
        name = name == "expect" ? "\\mathbb{E}" : string("\\text{", name, "}")
        return string(name, "_{", data_str, "}",
                  InfiniteOpt._infopt_math_symbol(m, :open_rng), func_str,
                  InfiniteOpt._infopt_math_symbol(m, :close_rng))
    end
end

################################################################################
#                          VARIABLE STRING METHODS
################################################################################
## helper function for getting the variable names
# REPLMode
function _get_base_name(::Type{JuMP.REPLMode}, vref)::String
    var_name = JuMP.name(vref)
    if !isempty(var_name)
        return var_name
    else
        return "noname"
    end
end

# IJuliaMode
function _get_base_name(::Type{JuMP.IJuliaMode}, vref)::String
    var_name = JuMP.name(vref)
    if !isempty(var_name)
        # TODO: This is wrong if variable name constains extra "]"
        return replace(replace(var_name, "[" => "_{", count = 1), "]" => "}")
    else
        return "noname"
    end
end

# Helper method for infinite variable construction
function _add_on_parameter_refs(base_name::String, prefs::VectorTuple)::String
    param_name_tuple = "("
    for i in 1:size(prefs, 1)
        element_prefs = prefs[i, :]
        type = _index_type(first(element_prefs))
        if type == DependentParameterIndex
            param_name = _remove_name_index(first(element_prefs))
        elseif length(element_prefs) == 1
            param_name = JuMP.name(first(element_prefs))
        else
            # TODO this isn't quite right with a subset of an independent container
            names = map(p -> _remove_name_index(p), element_prefs)
            if _allequal(names)
                param_name = first(names)
            else
                param_name = string("[", join(element_prefs, ", "), "]")
            end
        end
        if i != size(prefs, 1)
            param_name_tuple *= string(param_name, ", ")
        else
            param_name_tuple *= string(param_name, ")")
        end
    end
    return string(base_name, param_name_tuple)
end

# Make a string for InfiniteVariableRef
function variable_string(print_mode, vref::InfiniteVariableRef)::String
    base_name = _get_base_name(print_mode, vref)
    if !haskey(_data_dictionary(vref), JuMP.index(vref))
        return base_name
    else
        prefs = raw_parameter_refs(vref)
        return _add_on_parameter_refs(base_name, prefs)
    end
end

# Make a string for ParameterFunctionRef
function variable_string(print_mode, vref::ParameterFunctionRef)::String
    base_name = _get_base_name(print_mode, vref)
    if !haskey(_data_dictionary(vref), JuMP.index(vref))
        return base_name
    else
        prefs = raw_parameter_refs(vref)
        return _add_on_parameter_refs(base_name, prefs)
    end
end

## Make helper function for making derivative operators
# REPL
function _deriv_operator(m::Type{JuMP.REPLMode}, pref)::String
    return string(_infopt_math_symbol(m, :partial), "/",
                  _infopt_math_symbol(m, :partial), variable_string(m, pref))
end

# IJulia
function _deriv_operator(m::Type{JuMP.IJuliaMode}, pref)::String
    return string("\\frac{", _infopt_math_symbol(m, :partial), "}{",
                  _infopt_math_symbol(m, :partial), " ",
                  variable_string(m, pref), "}")
end

# TODO implement more intelligent naming for nested derivatives (i.e., use exponents)
# TODO account for container naming when variable macro is used (maybe deal with this at the macro end)
# Make a string for DerivativeRef
function variable_string(print_mode, dref::DerivativeRef)::String
    base_name = _get_base_name(print_mode, dref)
    if !haskey(_data_dictionary(dref), JuMP.index(dref))
        return base_name
    elseif base_name != "noname"
        prefs = raw_parameter_refs(dref)
        return _add_on_parameter_refs(base_name, prefs)
    else
        vref = dispatch_variable_ref(derivative_argument(dref))
        pref = operator_parameter(dref)
        return string(_deriv_operator(print_mode, pref),
                      _infopt_math_symbol(print_mode, :open_rng),
                      variable_string(print_mode, vref),
                      _infopt_math_symbol(print_mode, :close_rng))
    end
end

## Return the parameter value as an appropriate string
# Number
function _make_str_value(value)::String
    return JuMP._string_round(value)
end

# Array{<:Number}
function _make_str_value(values::Array)::String
    if length(values) == 1
        return _make_str_value(first(values))
    end
    if length(values) <= 4
        str_value = "["
        for i in eachindex(values)
            if i != length(values)
                str_value *= JuMP._string_round(values[i]) * ", "
            else
                str_value *= JuMP._string_round(values[i]) * "]"
            end
        end
        return str_value
    else
        return string("[", JuMP._string_round(first(values)), ", ..., ",
                      JuMP._string_round(last(values)), "]")
    end
end

# Make a string for PointVariableRef
# TODO improve so numerator of derivative contains the point
function variable_string(print_mode, vref::PointVariableRef)::String
    if !haskey(_data_dictionary(vref), JuMP.index(vref)) || !isempty(JuMP.name(vref))
        return _get_base_name(print_mode, vref)
    else
        ivref = dispatch_variable_ref(infinite_variable_ref(vref))
        if ivref isa InfiniteVariableRef || !isempty(JuMP.name(ivref))
            base_name = _get_base_name(print_mode, ivref)
        else
            base_name = variable_string(print_mode, ivref) # we have a derivative
        end
        prefs = raw_parameter_refs(ivref)
        values = raw_parameter_values(vref)
        name = string(base_name, "(")
        for i in 1:size(prefs, 1)
            if i != size(prefs, 1)
                name *= string(_make_str_value(values[prefs.ranges[i]]), ", ")
            else
                name *= string(_make_str_value(values[prefs.ranges[i]]), ")")
            end
        end
        return name
    end
end

# Make a string for ReducedVariableRef
# TODO improve so numerator of derivative contains the parameter tuple
function variable_string(print_mode, vref::ReducedVariableRef)::String
    if !haskey(_data_dictionary(vref), JuMP.index(vref)) || !isempty(JuMP.name(vref))
        return _get_base_name(print_mode, vref)
    else
        ivref = dispatch_variable_ref(infinite_variable_ref(vref))
        if ivref isa InfiniteVariableRef || !isempty(JuMP.name(ivref))
            base_name = _get_base_name(print_mode, ivref)
        else
            base_name = variable_string(print_mode, ivref) # we have a derivative
        end
        prefs = raw_parameter_refs(ivref)
        eval_supps = eval_supports(vref)
        raw_list = [i in keys(eval_supps) ? eval_supps[i] : prefs[i]
                    for i in eachindex(prefs)]
        param_name_tuple = "("
        for i in 1:size(prefs, 1)
            value = raw_list[prefs.ranges[i]]
            if i != size(prefs, 1)
                param_name_tuple *= string(_make_str_value(value), ", ")
            else
                param_name_tuple *= string(_make_str_value(value), ")")
            end
        end
        return string(base_name, param_name_tuple)
    end
end

# Fallback
function variable_string(print_mode, vref::JuMP.AbstractVariableRef)::String
    return _get_base_name(print_mode, vref)
end

# Extend function string for DispatchVariableRefs (REPL)
function JuMP.function_string(::Type{JuMP.REPLMode},
                              vref::DispatchVariableRef)::String
    return variable_string(JuMP.REPLMode, vref)
end

# Extend function string for DispatchVariableRefs (IJulia)
function JuMP.function_string(::Type{JuMP.IJuliaMode},
                              vref::DispatchVariableRef)::String
    return variable_string(JuMP.IJuliaMode, vref)
end

# Extend function string for GeneralVariableRefs (REPL)
function JuMP.function_string(::Type{JuMP.REPLMode},
                              vref::GeneralVariableRef)::String
    return variable_string(JuMP.REPLMode, dispatch_variable_ref(vref))
end

# Extend function string for GeneralVariableRefs (IJulia)
function JuMP.function_string(::Type{JuMP.IJuliaMode},
                              vref::GeneralVariableRef)::String
    return variable_string(JuMP.IJuliaMode, dispatch_variable_ref(vref))
end

################################################################################
#                         CONSTRAINT STRING METHODS
################################################################################
# Return parameter bound list as a string
function bound_string(print_mode, bounds::ParameterBounds{GeneralVariableRef}
                      )::String
    string_list = ""
    for (pref, set) in bounds
        string_list *= string(JuMP.function_string(print_mode, pref), " ",
                              JuMP.in_set_string(print_mode, set), ", ")
    end
    return string_list[1:end-2]
end

## Return the parameter domain string given the object index
# IndependentParameter
function _param_domain_string(print_mode, model::InfiniteModel,
                              index::IndependentParameterIndex,
                              bounds::ParameterBounds{GeneralVariableRef}
                              )::String
    pref = dispatch_variable_ref(model, index)
    set = infinite_set(pref)
    gvref = GeneralVariableRef(model, MOIUC.key_to_index(index), typeof(index))
    return string(JuMP.function_string(print_mode, pref), " ",
                  JuMP.in_set_string(print_mode, gvref, set, bounds))
end

# DependentParameters
function _param_domain_string(print_mode, model::InfiniteModel,
                              index::DependentParametersIndex,
                              bounds::ParameterBounds{GeneralVariableRef}
                              )::String
    # parse the infinite set
    first_gvref = GeneralVariableRef(model, MOIUC.key_to_index(index),
                                     DependentParameterIndex, 1)
    set = _parameter_set(dispatch_variable_ref(first_gvref))
    # iterate over the individual scalar sets if is a collection set
    if set isa CollectionSet
        domain_str = ""
        for i in eachindex(collection_sets(set))
            sset = collection_sets(set)[i]
            gvref = GeneralVariableRef(model, MOIUC.key_to_index(index),
                                       DependentParameterIndex, i)
            domain_str *= string(JuMP.function_string(print_mode, gvref), " ",
                        JuMP.in_set_string(print_mode, gvref, sset, bounds), ", ")
        end
        domain_str = domain_str[1:end-2]
    else
        # determine if bounds contain equalities and filter to the bounds of interest
        is_eq = haskey(bounds, first_gvref) && JuMP.lower_bound(bounds[first_gvref]) == JuMP.upper_bound(bounds[first_gvref])
        filtered_bounds = filter(e -> e[1].raw_index == MOIUC.key_to_index(index) &&
                                 e[1].index_type == DependentParameterIndex, bounds)
        # build the domain string
        if is_eq
            domain_str = bound_string(print_mode, filtered_bounds)
        else
            domain_str = string(_remove_name_index(first_gvref), " ",
                            JuMP.in_set_string(print_mode, set))
            if !isempty(filtered_bounds)
                domain_str *= string(" ", _infopt_math_symbol(print_mode, :intersect),
                                     " (", bound_string(print_mode, filtered_bounds), ")")
            end
        end
    end
    return domain_str
end

# Return string of an infinite constraint
function JuMP.constraint_string(print_mode, cref::InfOptConstraintRef;
                                in_math_mode = false)::String
    # get the function and set strings
    func_str = JuMP.function_string(print_mode, _core_constraint_object(cref))
    in_set_str = JuMP.in_set_string(print_mode, _core_constraint_object(cref))
    # check if constraint if finite
    obj_nums = _object_numbers(cref)
    if isempty(obj_nums)
        bound_str = ""
    else
        # get the parameter bounds if there are any
        if has_parameter_bounds(cref)
            bounds = parameter_bounds(cref)
        else
            bounds = ParameterBounds()
        end
        # prepare the parameter domains
        model = JuMP.owner_model(cref)
        bound_str = string(", ", JuMP._math_symbol(print_mode, :for_all), " ")
        for index in _param_object_indices(model)[_object_numbers(cref)]
            bound_str *= string(_param_domain_string(print_mode, model, index, bounds),
                                ", ")
        end
        bound_str = bound_str[1:end-2]
    end
    # form the constraint string
    if print_mode == JuMP.REPLMode
        lines = split(func_str, '\n')
        lines[1 + div(length(lines), 2)] *= " " * in_set_str * bound_str
        constr_str = join(lines, '\n')
    else
        constr_str = string(func_str, " ", in_set_str, bound_str)
    end
    # format for IJulia
    if print_mode == JuMP.IJuliaMode && !in_math_mode
        constr_str = JuMP._wrap_in_inline_math_mode(constr_str)
    end
    # add name if it has one
    name = JuMP.name(cref)
    if isempty(name) || (print_mode == JuMP.IJuliaMode && in_math_mode)
        return constr_str
    else
        return string(name, " : ", constr_str)
    end
end

# return list of constraint
function JuMP.constraints_string(print_mode, model::InfiniteModel)::Vector{String}
    # allocate the string vector
    strings = Vector{String}(undef, JuMP.num_constraints(model))
    # produce a string for each constraint
    counter = 1
    for (index, object) in model.constraints
        cref = _make_constraint_ref(model, index)
        strings[counter] = JuMP.constraint_string(print_mode, cref,
                                                  in_math_mode = true)
        counter += 1
    end
    return strings
end

################################################################################
#                           OBJECTIVE STRING METHODS
################################################################################
# Return the objective string
function JuMP.objective_function_string(print_mode, model::InfiniteModel)
    return JuMP.function_string(print_mode, JuMP.objective_function(model))
end

################################################################################
#                           DATATYPE SHOWING METHODS
################################################################################
# Show infinite sets in REPLMode
function Base.show(io::IO, set::AbstractInfiniteSet)
    print(io, set_string(JuMP.REPLMode, set))
    return
end

# Show infinite sets in IJuliaMode
function Base.show(io::IO, ::MIME"text/latex", set::AbstractInfiniteSet)
    print(io, set_string(JuMP.IJuliaMode, set))
end

# TODO use ... when necessary --> Subdomain bounds: t in [0, 1], x[1] == 0, x[2] == 0, ... x[6] == 0, a in [2, 3]
# Show ParameterBounds in REPLMode
function Base.show(io::IO, bounds::ParameterBounds)
    print(io, "Subdomain bounds (", length(bounds), "): ",
          bound_string(JuMP.REPLMode, bounds))
    return
end

# Show ParameterBounds in IJuliaMode
function Base.show(io::IO, ::MIME"text/latex", bounds::ParameterBounds)
    print(io, "Subdomain bounds (", length(bounds), "): ",
          bound_string(JuMP.IJuliaMode, bounds))
end

# Show constraint in REPLMode
function Base.show(io::IO, ref::InfOptConstraintRef)
    print(io, JuMP.constraint_string(JuMP.REPLMode, ref))
    return
end

# Show constraint in IJuliaMode
function Base.show(io::IO, ::MIME"text/latex", ref::InfOptConstraintRef)
    print(io, JuMP.constraint_string(JuMP.IJuliaMode, ref))
end

# Show the backend information associated with the optimizer model
function JuMP.show_backend_summary(io::IO, model::InfiniteModel)
    println(io, "Optimizer model backend information: ")
    JuMP.show_backend_summary(io, optimizer_model(model))
    return
end

# Show the objective function type
function JuMP.show_objective_function_summary(io::IO, model::InfiniteModel)
    println(io, "Objective function type: ",
            JuMP.objective_function_type(model))
    return
end

# Return constraint summary in JuMP like manner
function JuMP.show_constraints_summary(io::IO, model::InfiniteModel)
    for (F, S) in JuMP.list_of_constraint_types(model)
        n_constraints = JuMP.num_constraints(model, F, S)
        println(io, "`$F`-in-`$S`: $n_constraints constraint",
                _plural(n_constraints))
    end
    return
end

# Have the infinitemodel object show similar to Model
function Base.show(io::IO, model::InfiniteModel)
    # indicate is an infinite model
    println(io, "An InfiniteOpt Model")
    # show objective ssnse info
    sense = JuMP.objective_sense(model)
    if sense == MOI.MAX_SENSE
        print(io, "Maximization")
    elseif sense == MOI.MIN_SENSE
        print(io, "Minimization")
    else
        print(io, "Feasibility")
    end
    println(io, " problem with:")
    # show finite parameter info
    num_finite_params = num_parameters(model, FiniteParameter)
    println(io, "Finite Parameter", _plural(num_finite_params), ": ",
            num_finite_params)
    # show infinite parameter info
    num_infinite_params = num_parameters(model, InfiniteParameter)
    println(io, "Infinite Parameter", _plural(num_infinite_params), ": ",
            num_infinite_params)
    # show variable info
    num_vars = JuMP.num_variables(model)
    println(io, "Variable", _plural(num_vars), ": ", num_vars)
    # show the derivative info
    num_derivs = num_derivatives(model)
    println(io, "Derivative", _plural(num_derivs), ": ", num_derivs)
    # show measure info
    num_meas = num_measures(model)
    println(io, "Measure", _plural(num_meas), ": ", num_meas)
    # show objective function info
    if sense != MOI.FEASIBILITY_SENSE
        JuMP.show_objective_function_summary(io, model)
    end
    # show constraint info
    JuMP.show_constraints_summary(io, model)
    # show other info
    names_in_scope = sort!(collect(keys(JuMP.object_dictionary(model))))
    if !isempty(names_in_scope)
        print(io, "Names registered in the model: ",
              join(string.(names_in_scope), ", "), "\n")
    end
    JuMP.show_backend_summary(io, model)
end


text:
using Hecke

using Test

using LinearAlgebra

for scope in Hecke.ASSERT_SCOPE
  set_assert_level(scope, 3)
end

#for scope in Hecke.VERBOSE_SCOPE
#  set_verbose_level(scope, 3)
#end

k, a = quadratic_field(5)
@assert fmpz(1) - a == -(a - 1)
@assert 1 - a == -(a - 1)

@time include("NumField.jl")
@time include("AlgAss.jl")
@time include("AlgAssAbsOrd.jl")
@time include("AlgAssRelOrd.jl")
@time include("EllCrv.jl")
@time include("GrpAb.jl")
@time include("Grp.jl")
@time include("LinearAlgebra.jl")
@time include("Map.jl")
@time include("Misc.jl")
@time include("NfAbs.jl")
@time include("NfOrd.jl")
@time include("NfRel.jl")
@time include("RCF.jl")
@time include("Examples.jl")
@time include("Sparse.jl")
@time include("QuadForm.jl")
@time include("LocalField.jl")


text:
@testset "canonical domains" begin
    @test canonicaldomain(2..3) isa ChebyshevInterval
    @test canonicaldomain(2..3) == ChebyshevInterval{Float64}()
    @test mapto(UnitInterval(), UnitInterval()) isa IdentityMap
    @test mapto(UnitInterval(), ChebyshevInterval()) isa AffineMap
end


text:
module Infernal

include("cmd.jl")

end


text:
"""
    project_reach(Rsets, vars, options)

Projection of a reachability analysis result in 2D.

### Input

- `Rsets`   -- reachable states representation
- `vars`    -- variables to plot; two-dimensional index vector
- `options` -- options

### Notes

The `vars` argument is required even if the optional argument
`projection_matrix` is passed, because we also determine whether time is used as
a dimension from this variable.
"""
function project_reach(
        Rsets::Vector{<:AbstractReachSet{<:LazySets.LazySet{N}}},
        vars::Vector{Int64},
        options::AbstractOptions)::Vector{<:AbstractReachSet} where {N<:Real}
    # parse input
    @assert length(vars) == 2 "we only support projection to two dimensions"
    xaxis = vars[1]
    yaxis = vars[2]
    got_time = (xaxis == 0)
    if xaxis < 0
        throw(DomainError("value $xaxis for x variable not allowed"))
    elseif yaxis <= 0
        throw(DomainError("value $yaxis for y variable not allowed"))
    end
    n = options[:n]
    projection_matrix_high_dimensional = options[:projection_matrix]
    transformation_matrix = options[:transformation_matrix]

    # allocate output and define overapproximation function
    ε = options[:ε_proj]
    if ε < Inf
        oa = x -> overapproximate(x, HPolygon, ε)
        RsetsProj = Vector{ReachSet{HPolygon{N}}}(undef, length(Rsets))
    else
        set_type = options[:set_type_proj]
        oa = x -> overapproximate(x, set_type)
        RsetsProj = Vector{ReachSet{<:set_type{N}}}(undef, length(Rsets))
    end

    function _project(rs)
        X = set(rs)
        if projection_matrix_high_dimensional == nothing
            # use a simple projection to state variables
            if got_time
                # projection to a single state variable
                projection_matrix = sparse([1], [yaxis], [1.0], 1, n)
            else
                # projection to two state variables
                projection_matrix =
                    sparse([1, 2], [xaxis, yaxis], [1.0, 1.0], 2, n)
            end
        else
            @assert size(projection_matrix_high_dimensional, 1) == 1 "currently " *
                "we only support one-dimensional projection matrices"
            if got_time
                # create a 1-row matrix
                projection_matrix =
                    sparse(fill(1, n), 1:n, projection_matrix_high_dimensional[1, :], 1, n)
            else
                # create a 2-row matrix
                projection_matrix =
                    sparse(fill(2, n), 1:n, projection_matrix_high_dimensional[1, :], 2, n)
                projection_matrix[1, xaxis] = 1.0
            end
        end

        # apply optional transformation to projection matrix
        if transformation_matrix != nothing
            if got_time
                # add another dimension for time: block matrix [S 0; 0 1]
                transformation_matrix =
                    sparse(cat([1, 2], transformation_matrix, [1]))
            end
            projection_matrix = projection_matrix * transformation_matrix
        end

        # project set
        rs = project(rs, projection_matrix)
        projected = set(rs)

        # add time dimension
        t0 = time_start(rs)
        t1 = time_end(rs)
        if got_time
            time_interval = Interval(t0, t1)
            projected = CartesianProduct(time_interval, projected)
        end
        return ReachSet(oa(projected), t0, t1)
    end

    return [_project(rs) for rs in Rsets]
end

"""
    project(Rsets, options)

Projects a sequence of sets according to the settings defined in the options.

### Input

- `Rsets`   -- solution of a reachability problem
- `options` -- options structure

### Notes

A projection matrix can be given in the options structure, or passed as a
dictionary entry.
"""
function project(Rsets::Vector{<:AbstractReachSet}, options::AbstractOptions)
    return project_reach(Rsets, options[:plot_vars], options)
end

function project(flowpipe::Flowpipe, options::AbstractOptions)
    return Flowpipe(project(flowpipe.reachsets, options))
end

function project(reach_sol::ReachSolution)
    proj_flowpipes = Vector{Flowpipe}(undef, length(reach_sol.flowpipes))
    for (i, flowpipe) in enumerate(reach_sol.flowpipes)
        proj_flowpipes[i] = project(flowpipe, reach_sol.options)
    end
    return ReachSolution(proj_flowpipes, reach_sol.options)
end

project(Rsets::Vector{<:AbstractReachSet}, options::Pair{Symbol,<:Any}...) =
    project(Rsets, Options(Dict{Symbol,Any}(options)))


text:
"""
	histogram(cmd0::String="", arg1=nothing; kwargs...)

Examines the first data column to calculate histogram parameters based on the bin-width provided.
Alternatively, show histograms of GMTimage & GMTgrid objects directly. In the GMTimage case, the
options 'auto=true' or 'thresholds=(0, 0.1)' will find the histogram bounds of UInt16 images
convenient for contrast enhancement (histogram stretch). The values represent the percentage of
countings used to estimate the boundings. The option 'zoom=true' will set 'auto=true' and show
histogram only on the region of interest.

Full option list at [`pshistogram`]($(GMTdoc)histogram.html)

Parameters
----------

- $(GMT.opt_J)
- **A** | **horizontal** :: [Type => Bool]

    Plot the histogram horizontally from x = 0 [Default is vertically from y = 0].
    ($(GMTdoc)histogram.html#a)
- $(GMT.opt_B)
- **C** | **color** | **cmap** :: [Type => Str | GMTcpt]

    Give a CPT. The mid x-value for each bar is used to look-up the bar color.
    ($(GMTdoc)histogram.html#c)
- **D** | **annot** | **annotate** | **counts** :: [Type => Str | Tuple]

    Annotate each bar with the count it represents.
    ($(GMTdoc)histogram.html#d)
- **E** | **width** :: [Type => Bool]			`Arg = width[+ooffset]`

    Use an alternative histogram bar width than the default set via T, and optionally shift all bars by an offset.
    ($(GMTdoc)histogram.html#e)
- **F** | **center** :: [Type => Bool]

    Center bin on each value. [Default is left edge].
    ($(GMTdoc)histogram.html#f)
- **G** | **fill** :: [Type => Number | Str]

    Select filling of bars [if no G, L or C set G=100].
    ($(GMTdoc)histogram.html#g)
- **I** | **inquire** | **bins** :: [Type => Bool | :O | :o | bins=(all=true,) | bins=(no_zero=true,) ]

    Inquire about min/max x and y after binning OR output the binned array.
    ($(GMTdoc)histogram.html#i)
- **L** | **out_range** :: [Type => Str]			`Arg = l|h|b`

    Handling of extreme values that fall outside the range set by **T**.
    ($(GMTdoc)histogram.html#l)
- **N** | **distribution** :: [Type => Str]

    Draw the equivalent normal distribution; append desired pen [0.5p,black].
    ($(GMTdoc)histogram.html#n)
- $(GMT.opt_P)
- **Q** | **cumulative** :: [Type => Bool | "r"]

    Draw a cumulative histogram. Append r to instead compute the reverse cumulative histogram.
    ($(GMTdoc)histogram.html#q)
- **R** | **region** :: [Type => Str]

    Specifies the ‘region’ of interest in (r,azimuth) space. r0 is 0, r1 is max length in units.
    ($(GMTdoc)histogram.html#r)
- **S** | **stairs** :: [Type => Str | number]

    Draws a stairs-step diagram which does not include the internal bars of the default histogram.
    ($(GMTdoc)histogram.html#s)
- **T** | **range** | **bin** :: [Type => Str]			`Arg = [min/max/]inc[+n] | file|list]`

    Defines the range of the new CPT by giving the lowest and highest z-value and interval.
    ($(GMTdoc)histogram.html#t)
- **W** | **pen** :: [Type => Str | Tuple]

	Set pen attributes for sector outline or rose plot. [Default is no outline].
	($(GMTdoc)histogram.html#w)
- **Z** | **kind** :: [Type => Number | Str]

    Choose between 6 types of histograms.
    ($(GMTdoc)histogram.html#z)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function histogram(cmd0::String="", arg1=nothing; first=true, kwargs...)

	arg2 = nothing		# May be needed if GMTcpt type is sent in via C
	N_args = (arg1 === nothing) ? 0 : 1

	gmt_proggy = (IamModern[1]) ? "histogram "  : "pshistogram "
	length(kwargs) == 0 && return monolitic(gmt_proggy, cmd0, arg1, arg2)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode

	cmd = ""
	opt_Z = add_opt(d, "", 'Z', [:Z :kind], (counts="0", count="0", freq="1", log_count="2", log_freq="3",
	                                         log10_count="4", log10_freq="5", weights="+w"), true, "")
	opt_T = parse_opt_range(d, "", "")		# [:T :range :inc :bin]
	(isa(arg1, GMTimage) || isa(arg1, GMTgrid)) && occursin("/", opt_T) && error("here 'inc' must be a scalar")

	# If inquire, no plotting so do it and return
	opt_I = add_opt(d, "", "I", [:I :inquire :bins], (all="O", no_zero="o"))
	if (opt_I != "")
		cmd *= opt_I
		if ((r = dbg_print_cmd(d, cmd)) !== nothing)  return r  end
		cmd, arg1, = read_data(d, cmd0, cmd, arg1, " ")
		if (isa(arg1, GMTimage))		# If it's an image with no bin option, default to bin=1
			arg1, cmd = loc_histo(arg1, cmd, opt_T, opt_Z)
		else
			cmd *= opt_Z
			if (opt_T != "")  cmd *= " -T" * opt_T  end
		end
		return gmt(gmt_proggy * cmd, arg1)
	end

	K, O = set_KO(first)		# Set the K O dance

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, cmd, "histogram", O, " -JX12c/12c")
	cmd = parse_common_opts(d, cmd, [:UVXY :JZ :c :e :p :t :params], first)[1]
	cmd = parse_these_opts(cmd, d, [[:A :horizontal], [:F :center], [:Q :cumulative], [:S :stairs]])
	cmd = add_opt_fill(cmd, d, [:G :fill], 'G')
	cmd = add_opt(d, cmd, 'D', [:D :annot :annotate :counts], (beneath="_+b", font="+f", offset="+o", vertical="_+r"))
	cmd = parse_INW_coast(d, [[:N :distribution :normal]], cmd, "N")
	(show_kwargs[1]) && print_kwarg_opts(symbs, "NamedTuple | Tuple | Dict | String")
	(isa(arg1, GMTimage) || isa(arg1, GMTgrid)) && occursin("/", opt_T) && error("here 'inc' must be a scalar")

	if (GMTver >= 6.2)  cmd = add_opt(d, cmd, 'E', [:E :width], (width="", off="+o", offset="+o"))  end

	# If file name sent in, read it and compute a tight -R if this was not provided
	if (opt_R == "")  opt_R = " "  end			# So it doesn't try to find the -R in next call
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R)
	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_args, arg1, arg2)

	cmd  = add_opt(d, cmd, 'L', [:L :out_range], (first="l", last="h", both="b"))
	cmd *= add_opt_pen(d, [:W :pen], "W", true)     	# TRUE to also seek (lw|lt,lc,ls)
	if (!occursin("-G", cmd) && !occursin("-C", cmd) && !occursin("-S", cmd))
		cmd *= " -G150"
	elseif (occursin("-S", cmd) && !occursin("-W", cmd))
		cmd *= " -W0.3p"
	end

	limit_L = nothing
	if (isa(arg1, GMTimage))		# If it's an image with no bin option, default to bin=1
		do_clip = (isa(arg1[1], UInt16) && (val = find_in_dict(d, [:full_histo])[1]) === nothing) ? true : false
		do_auto = ((val_auto = find_in_dict(d, [:auto :threshols])[1]) !== nothing) ? true : false	# Automatic bounds detetion
		do_zoom = ((find_in_dict(d, [:zoom])[1]) !== nothing) ? true : false	# Automatic zoom to interesting region
		if (do_zoom && !do_auto)  val_auto = nothing  end	# I.e. 'zoom' sets also the auto mode
		hst, cmd = loc_histo(arg1, cmd, opt_T, opt_Z)
		if (do_clip && (all(hst[3:10,2] .== 0)))  hst[1,2] = 0; hst[2,2] = 0  end
		if (do_auto || do_zoom)
			limit_L, limit_R = find_histo_limits(arg1, val_auto, 200)
			if (do_zoom)
				mm = extrema(hst, dims=1)		# 1×2 Array{Tuple{UInt16,UInt16},2}
				x_max = min(limit_R * 1.15, hst[end,1])		# 15% to the right but not fall the cliff
				opt_R_ = " -R$(limit_L*0.85)/$x_max/0/$(mm[2][2] * 1.1) "
				(opt_R != " ") && @warn("'zoom' option overrides the requested region limits and sets its own")
				cmd = replace(cmd, opt_R => opt_R_, count=1)
				opt_R = opt_R_					# It will be needed further down in vlines
			end
		end
		arg1 = hst		# We want to send the histogram, not the GMTimage
	elseif (isa(arg1, GMTgrid))
		if (opt_T != "")
			inc = parse(Float64, opt_T) + eps()		# + EPS to avoid the extra last bin at right with 1 count only
			n_bins = Int(ceil((arg1.range[6] - arg1.range[5]) / inc))
		else
			n_bins = Int(ceil(sqrt(length(arg1.z))))
			inc = (arg1.range[6] - arg1.range[5]) / n_bins + eps()
		end
		(!isa(inc, Real) || inc <= 0) && error("Bin width must be a > 0 number and no min/max")
		hst = zeros(n_bins,2)
		@inbounds for k = 1:length(arg1)  hst[Int(floor((arg1[k]-arg1.range[5])/inc)+1), 2] += 1  end
		[@inbounds hst[k,1] = arg1.range[5] + inc*(k-1) for k = 1:n_bins]

		cmd = (opt_Z == "") ? cmd * " -Z0" : cmd * opt_Z
		if (!occursin("+w", cmd))  cmd *= "+w"  end		# Pretending to be weighted is crutial for the trick
		cmd *= " -T$(arg1.range[5])/$(arg1.range[6])/$inc"
		arg1 = hst		# We want to send the histogram, not the GMTgrid
	else
		if (opt_T != "")  opt_T = " -T" * opt_T  end	# It lacked the -T so that it could be used in loc_histo()
		cmd *= opt_T * opt_Z
	end

	# The following looks a bit messy but it's needed to auto plotting verical lines with the limits
	show_ = false;		fmt_ = "ps";		savefig_ = nothing
	if (limit_L !== nothing)
		if (haskey(d, :show))  show_ = (d[:show] != 0)  end				# Backup the :show val
		d[:show] = false
		if (haskey(d, :fmt))  fmt_ = d[:fmt]; delete!(d, :fmt)  end		# Backup the :show val
		if ((val = find_in_dict(d, [:savefig :figname :name])[1]) !== nothing) savefig_ = val  end
	end

	out2 = nothing;		Vd_ = 0				# Backup values
	if (haskey(d, :Vd))  Vd_ = d[:Vd]  end

	# Plot the histogram
	out1 = finish_PS_module(d, gmt_proggy * cmd, "", K, O, true, arg1, arg2)

	# And if wished, plot the two vertical lines with the limits annotated in them
	if (limit_L !== nothing)
		if (opt_R == " ")					# Set a region for the vlines
			mm = extrema(hst, dims=1)
			opt_R = " -R$(mm[1][1])/$(mm[1][2])/0/$(mm[2][2])"
		end
		vlines!([limit_L], pen="0.5p,blue,dashed", decorated=(quoted=true, n_labels=1, const_label="$limit_L", font=9, pen=(0.5,:red)), R=opt_R[4:end], Vd=Vd_)
		out2 = vlines!([limit_R], pen="0.5p,blue,dashed", decorated=(quoted=true, n_labels=1, const_label="$limit_R", font=9, pen=(0.5,:red)), R=opt_R[4:end], fmt=fmt_, savefig=savefig_, show=show_, Vd=Vd_)
	end
	out = (out1 !== nothing && out2 !== nothing) ? [out1;out2] : ((out1 !== nothing) ? out1 : out2)

end

# ---------------------------------------------------------------------------------------------------
function find_histo_limits(In, thresholds=nothing, width=200)
	# Find the histogram limits of a UInt16 GMTimage that allow to better stretch the histogram
	# THRESHOLDS is an optional Tuple input containing the left and right thresholds, in percentage,
	# between which the histogram values will be retained. Defaults are (0,0.5) percent
	# WIDTH is bin width used to obtain a rough histogram that is used to compute the limits.
	if (isa(In, Array{UInt16,3}) || isa(In, Array{UInt8,3}))
		L1 = find_histo_limits(view(In,:,:,1), thresholds, width)
		L2 = find_histo_limits(view(In,:,:,2), thresholds, width)
		L3 = find_histo_limits(view(In,:,:,3), thresholds, width)
		return (L1[1], L1[2], L2[1], L2[2], L3[1], L3[2])
	end
	hst = loc_histo(In, "", string(width), "")[1]
	if (size(hst,1) >= 5)
		all(hst[2:5,2] .== 0) && (hst[1,2] = 0)	# Here we always check for high counts in zero bin
	end
	max_ = maximum(hst, dims=1)[2]
	(max_ == 0) && error("This histogram had nothing but countings ONLY in first bin. No point to proceed.")
	thresh_l = 0.0;		thresh_r = 0.005
	if (isa(thresholds, Tuple) && length(thresholds) == 2)
		thresh_l, thresh_r = thresholds[:] ./ 100
	end
	thresh_l *= max_
	thresh_r *= max_
	kl = 1;		kr = size(hst,1)
	while (hst[kl,2] == 0 || hst[kl,2] < thresh_l)  kl += 1  end
	while (hst[kr,2] == 0 || hst[kr,2] < thresh_r)  kr -= 1  end
	return Int(hst[kl,1]), Int(min(hst[kr,1] + width, hst[end,1]))
end

# ---------------------------------------------------------------------------------------------------
function loc_histo(in, cmd::String="", opt_T::String="", opt_Z::String="")
	# Very simple function to compute histograms of images (integers)
	# We put the countings in a Mx2 arrray to trick GMT (pshistogram) to think it's recieving a weighted input.
	(!isa(in[1], UInt16) && !isa(in[1], UInt8)) && error("Only UInt8 or UInt16 image types allowed here")

	inc = (opt_T != "") ? Float64(Meta.parse(opt_T)) : 1.0
	(!isa(inc, Real) || inc <= 0) && error("Bin width must be a > 0 number and no min/max")

	n_bins = (isa(in[1], UInt8)) ? 256 : Int(ceil((maximum(in) + 1) / inc))	# For UInt8 use the full [0 255] range
	hst = zeros(n_bins,2)
	pshst_wall!(in, hst, inc, n_bins)

	cmd = (opt_Z == "") ? cmd * " -Z0" : cmd * opt_Z
	if (!occursin("+w", cmd))  cmd *= "+w"  end		# Pretending to be weighted is crutial for the trick

	return hst, cmd * " -T0/$(n_bins * inc)/$inc"
end

function pshst_wall!(in, hst, inc, n_bins::Int)
	# Function barrier for type instability. With the body of this in calling fun the 'inc' var
	# introduces a mysterious type instability and execution times multiply by 3.
	if (inc == 1)
		@inbounds for k = 1:length(in)  hst[in[k]+1, 2] += 1  end
	else
		@inbounds for k = 1:length(in)  hst[Int(floor(in[k]/inc)+1), 2] += 1  end
	end
	[@inbounds hst[k,1] = inc*(k-1) for k = 1:n_bins]
	return nothing
end

# ---------------------------------------------------------------------------------------------------
histogram!(cmd0::String="", arg1=nothing; kw...) = histogram(cmd0, arg1; first=false, kw...)
histogram(arg1; kw...)  = histogram("", arg1; first=true, kw...)
histogram!(arg1; kw...) = histogram("", arg1; first=false, kw...)

const pshistogram  = histogram			# Alias
const pshistogram! = histogram!			# Alias