text:
mutable struct SCPHeader
    sender::Vector{UInt8} # 64-byte sender string
    receiver::Vector{UInt8} # 64-byte receiver string
    len::Int32

    function SCPHeader(senderstr::String, receiverstr::String, len::Int32)
        sender_short = string_to_zeroed_null_terminated_array(senderstr, VIRES_HEADER_STRING_LEN)
        receiver_short = string_to_zeroed_null_terminated_array(receiverstr, VIRES_HEADER_STRING_LEN)
        SCPHeader(sender_short,receiver_short,len)
    end

    function SCPHeader(sender_in::Vector{UInt8},receiver_in::Vector{UInt8},len::Int32)
        sender = zeros(UInt8, VIRES_HEADER_STRING_LEN)
        receiver = zeros(UInt8, VIRES_HEADER_STRING_LEN)
        copy!(sender, 1, sender_in, 1, min(length(sender_in), VIRES_HEADER_STRING_LEN))
        copy!(receiver, 1, receiver_in, 1, min(length(receiver_in), VIRES_HEADER_STRING_LEN))
        sender[end] = UInt8(0) # ensure null-terminate
        receiver[end] = UInt8(0)

        new(sender, receiver, len)
    end
end
mutable struct SCPMessage
    header::SCPHeader
    payload::Vector{UInt8} # <- xml string

    function SCPMessage(sender::String, receiver::String, payloadstr::String)
        payload = string_to_zeroed_null_terminated_array(payloadstr)
        header = SCPHeader(sender, receiver, Int32(length(payload)))
        new(header, payload)
    end
    SCPMessage(payloadstr::String) = SCPMessage("SCP", "TaskControl", payloadstr)
    #SCPMessage(payloadstr::String) = SCPMessage("ScVis", "SCP", payloadstr) #Xiaobai

    function SCPMessage(sender::Vector{UInt8},receiver::Vector{UInt8}, payload::Vector{UInt8})
        if payload[end] != NULL
            push!(payload, NULL)
        end
        header = SCPHeader(sender, receiver, Int32(length(payload)))
        new(header, payload)
    end
end

function string_to_zeroed_null_terminated_array(str::String)
    strlen = length(str)
    arr = zeros(UInt8, strlen+1)
    arr[1:strlen] = convert(Vector{UInt8}, str)
    arr[end] = NULL
    arr
end
function string_to_zeroed_null_terminated_array(str::String, array_length::Int)
    arr = zeros(UInt8, array_length)
    len = min(array_length-1, length(str))
    arr[1:len] = convert(Vector{UInt8}, str)[1:len]
    arr
end
function string_from_buffer(arr::Vector{UInt8})
    nulindex = findfirst(arr, NULL)
    if nulindex == 0
        nulindex = length(arr)
    end
    convert(String,arr[1:nulindex-1])
end

function Base.write(io::IO, header::SCPHeader)
    write(io, VIRES_MAGIC_NUMBER)
    write(io, VIRES_VERSION_HEADER)
    write(io, header.sender)
    write(io, header.receiver)
    write(io, header.len)
end
function Base.write(io::IO, message::SCPMessage)
    write(io, message.header)
    write(io, message.payload)
    #print_with_color(:blue, STDOUT, "SENDING\n")
    #print_message(message)
end

function Base.read(io::IO, ::Type{SCPMessage}, already_has_magic_number::Bool=false, verbose::Bool=false)
    if !already_has_magic_number
        scan_for_value(io, VIRES_MAGIC_NUMBER)
    end

    version_number = read(io, UInt16)
    sender = Array{UInt8}(64)
    @assert readbytes!(io, sender) == 64
    receiver = Array{UInt8}(64)
    @assert readbytes!(io, receiver) == 64
    payloadsize = read(io, Int32)
    payload = Array{UInt8}(payloadsize)
    @assert readbytes!(io, payload) == payloadsize
    message = SCPMessage(sender,receiver,payload)
    if verbose
        println(STDOUT, "READING")
        print_message(message)
    end
    message
end
function write_bytehex_to_stdout(f::Function)
    # ex: write_bytehex_to_stdout((buf)->write(buf, 0x08))

    # io -> f(io)
    buf = IOBuffer()
    f(buf)
    println(STDOUT, bytes2hex(takebuf_array(buf)))
end
function print_message(message::SCPMessage)
    # println(STDOUT, "\tSENDER:       ", bytes2hex(message.header.sender))
    # println(STDOUT, "\tSENDER:       ", string_from_buffer(message.header.sender))
    # println(STDOUT, "\tRECEIVER:     ", bytes2hex(message.header.receiver))
    # println(STDOUT, "\tRECEIVER:     ", string_from_buffer(message.header.receiver))
    # println(STDOUT, "\tPAYLOAD SIZE: ", message.header.len)
    println(STDOUT, string_from_buffer(message.header.sender), " -> ", string_from_buffer(message.header.receiver), ": ")
    println(STDOUT, "\tPAYLOAD:      ", string_from_buffer(message.payload))
    # println(STDOUT, "\t              ", bytes2hex(message.payload))
end
function print_message_bytes(message::SCPMessage)
    dat = Array(UInt8, VIRES_HEADER_LEN + message.header.len + 1)
    buf = IOBuffer(dat, true, true)
    write(buf, message)
    println(STDOUT, "STRING BUFFER: ")
    println(STDOUT, bytes2hex(dat))
end
function idle_and_print_messages(io::IO, sleeptime::Float64)
    start_time = time()
    while time() - start_time < sleeptime
        message = read(io, SCPMessage)
        print_with_color(:red, "RECEIVED\n")
        print_message(message)
    end
end

function scan_for_value(io::IO, target::UInt8)
    #=
    Blocks the stream until you get the byte you want
    =#

    while read(io, UInt8) != target
        # do nothing
    end
    nothing
end
function scan_for_value(io::IO, value::UInt16)

    #=
    Blocks the stream until you get the bytes you want
    Note: value[1] must be found first, then value[2], etc.
    =#

    byte_index = 1
    lo = UInt8(value & 0x00FF)
    hi = UInt8((value & 0xFF00) >> 8)

    while byte_index ≤ 2
        val = read(io, UInt8)
        if byte_index == 1
            byte_index += Int(val == lo)
        else
            byte_index = (val == hi) ? 3 : 1
        end
    end
    nothing
end

function wait_for_packet_with_element(io::IO, name::String, elementname::String; timeout::Float64=TIMEOUT_DEFAULT)

    finished = false
    start_time = time()
    while !finished && time() - start_time < timeout
        received = read(io, SCPMessage)
        finished = message_contains_packet_with_element(received, name, elementname)
    end

    finished # if finished is true, we did not time out
end
function wait_for_mirrored_messsage(io::IO, message::SCPMessage; timeout::Float64=TIMEOUT_DEFAULT)

    sent::ETree = message_payload_to_etree(message)

    finished = false
    start_time = time()
    while !finished && time() - start_time < timeout
        received = read(io, SCPMessage)
        finished = (string_from_buffer(received.header.sender) == string_from_buffer(received.header.receiver) == CLIENT_NAME) &&
                sent == message_payload_to_etree(received)
    end

    finished # if finished = true, we did not time out
end
function write_and_wait_for_mirrored_message(io::IO, message::SCPMessage; timeout::Float64=TIMEOUT_DEFAULT)
    write(io, message)
    wait_for_mirrored_messsage(io, message, timeout=timeout)
end
function write_and_wait_for_mirrored_message(io::IO, message::SCPMessage, error_message::String; timeout::Float64=TIMEOUT_DEFAULT)
    write(io, message)
    if !wait_for_mirrored_messsage(io, message, timeout=timeout)
        error(error_message)
    end
    true
end

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
module GOF3R

using s3gof3r_jll

export s3stream, s3getfile, s3upload

"""
    s3stream(bucket, path)

Return an `IO` object streaming `path` from `bucket`.
"""
function s3stream(bucket, path)
    endpoint = AWS_ENDPOINT[]
    s3gof3r_jll.gof3r() do exe
        return open(`$exe get -b $(bucket) -k $(path) --endpoint=$(endpoint)`)
    end
end

"""
    s3getfile(bucket, path, outfile)

Write the file at `bucket:path` to `outfile`.
"""
function s3getfile(bucket, path, outfile)
    endpoint = AWS_ENDPOINT[]
    s3gof3r_jll.gof3r() do exe
        open(`$exe get -b $(bucket) -k $(path) --endpoint=$(endpoint)`) do io
            write(outfile, io)
        end
    end
end

"""
    s3upload(bucket, file, key)

Upload `file` to `bucket:key`.
"""
function s3upload(bucket, file, key)
    endpoint = AWS_ENDPOINT[]
    s3gof3r_jll.gof3r() do exe
        run(`$exe cp $(file) s3://$(bucket)/$(key) --endpoint=$(endpoint)`)
    end
end

"""
    s3stream(f, bucket, path)

Call `f(s3stream(bucket, path))`, ensuring that the stream is closed properly once `f`
has finished executing, including in the case of errors.
"""
function s3stream(f, bucket, path)
    stream = s3stream(bucket, path)
    ret = try
        f(stream)
    catch
        kill(stream)
        rethrow()
    finally
        close(stream.in)
    end
    success(stream) || Base.pipeline_error(stream)
    return ret
end

const AWS_ENDPOINT = Ref{String}("")

function __init__()
    # make sure we have a region set
    if !haskey(ENV, "AWS_REGION")
        ENV["AWS_REGION"] = "us-east-2"
    end
    AWS_ENDPOINT[] = get(ENV, "AWS_ENDPOINT", "s3.$(ENV["AWS_REGION"]).amazonaws.com")
end

end # module

ext: jl


text:
## ----------- Linear GradRule
# Defined by probGrad(h = hMin), which varies by college, and hMax where probGrad hits max (common to all colleges).

"""
	$(SIGNATURES)

Switches for linear graduation rule. It can operate on h levels or gains.

On h levels: Grad prob is linear in h between hMin and hMax. Each college has a calibrated values for prob grad(h = hMin).

On h gains: Grad prob is linear in (h - h0). Each college has a calibrated value for prob grad(h gain = hMin)
"""
Base.@kwdef mutable struct GradRuleSwitchesLinear <: GradRuleSwitches
    tFirstGrad :: TimeInt = 4
    "Number of courses required for graduation"
    minNcForGrad :: ncInt = 40
    "Number of 2 year colleges (cannot graduate from those)"
    n2Year :: ncInt = 1
    "Number of colleges"
    nColleges :: ncInt = 4
    # Min and max grad prob when graduation is feasible.
    gradProbMin :: Double = 0.05
    gradProbMax :: Double = 0.95
    # Grad prob depends on levels or gains
    useHcLevel :: Bool = true
    # Does grad rule vary by college
    byCollege :: Bool = true

    # h level (or h-h0 value) for which min grad prob is reached.
    hMin :: Double = 2.0
    # Are better colleges harder (or the other way around)?
    betterHarder :: Bool = true
end

Base.@kwdef mutable struct GradRuleSetLinear <: GradRuleSet
    @grad_rule_set_common
    # Grad prob at h = hMin (or h-h0 = 0)
    probHminV :: BoundedVector{Double}
    # h level (or h-h0) for which min grad prob is reached.
    hMin :: Double
    # h (or h - h0) at which gradProbMax is reached (the same for all colleges)
    hMax :: Double
    switches :: GradRuleSwitchesLinear
end

## For one college
struct GradRuleLinear <: GradRule
    # Can graduate after at least this many years
    tFirstGrad :: TimeInt
    # Grad prob at h = hMin  or at (h-h0) = hMin
    probHmin :: Double
    # h level (or h-h0) for which min grad prob is reached.
    hMin :: Double
    # Grad prob reaches max at h = hMax (or at h - h0 = hMax)
    hMax :: Double
    gradProbMin :: Double
    gradProbMax :: Double
    # Min no of courses attempted (cumulative)
    nMin :: ncInt
    useHcLevel :: Bool
end

function Base.show(io :: IO, g :: GradRuleLinear)
    print(io, "Linear graduation rule in")
    if use_hc_level(g)
        print(io, " h levels. ");
    else
        print(io, " h gains. ");
    end
    print(io, "hMin = $(h_min(g)), hMax = $(h_max(g))");
end

"""
	$(SIGNATURES)

Table with settings. Columns are explanation and values.
All formatted into strings.
"""
function settings_table(gs :: GradRuleSwitchesLinear)
    nc = n_colleges(gs);
    tFirstGrad = Int(t_first_grad(gs, nc));
    minNcToGrad = Int(min_courses_to_grad(gs, nc; modelUnits = false));
    harderStr = gs.betterHarder  ?  "harder"  :  "easier";
    return [
        "Graduation rule"  "Linear in h";
        "Time to graduation"  "$tFirstGrad";
        "Min courses for graduation"  "$minNcToGrad";
        "Better colleges are"  harderStr
    ]
end

function settings_list(gs :: GradRuleSwitchesLinear, st)
    nc = n_colleges(gs);
    tFirstGrad =  symbol_entry(st, :tGradMin, Int(t_first_grad(gs, nc)));
    minNcToGrad = symbol_entry(st, :nGrad, Int(min_courses_to_grad(gs, nc; modelUnits = false)));
    probGradMin = symbol_entry(st, :probGradMin, grad_prob_min(gs));
    probGradMax = symbol_entry(st, :probGradMax, grad_prob_max(gs));
    return [tFirstGrad, minNcToGrad, probGradMin, probGradMax]
end

use_hc_level(switches :: GradRuleSwitchesLinear) =
    switches.useHcLevel;

h_min(switches :: GradRuleSwitchesLinear) = switches.hMin;
h_min(gs :: GradRuleSetLinear) = h_min(gs.switches);
h_min(g :: GradRuleLinear) = g.hMin;

h_max(gs :: GradRuleSetLinear) = gs.hMax;
h_max(g :: GradRuleLinear) = g.hMax;

grad_colleges(gs :: GradRuleSetLinear) = grad_colleges(gs.switches);

grad_prob_hmin(g :: GradRuleLinear) = g.probHmin;

varies_by_college(switches :: GradRuleSwitchesLinear) =
    switches.byCollege;
grad_rule_same!(switches :: GradRuleSwitchesLinear) =
    switches.byCollege = false;

better_easier!(switches :: GradRuleSwitchesLinear) =
    switches.betterHarder = false;

"""
	$(SIGNATURES)

For all colleges.
Assumes college 1 is a 2 year college. Others are 4 year colleges.
"""
function make_grad_rule_set(objId :: ObjectId,
    grSwitches :: GradRuleSwitchesLinear)

    pProbHmin = init_prob_hmin(objId, grSwitches);
    hMin = h_min(grSwitches);

    if use_hc_level(grSwitches)
        hMax = Double(4.0);
        pHRatio = Param(:hMax, ldescription(:hGradMax), lsymbol(:hGradMax),
            hMax, hMax, hMin + Double(1.0), hMin + Double(6.0), true);
    else
        # Change in h at which grad prob hits max
        lb = hMin + Double(0.2);
        hMax = lb + Double(0.8);
        pHRatio = Param(:hMax, ldescription(:dhGradMax), lsymbol(:dhGradMax),
            hMax, hMax, lb, Double(5.0), true);
    end

    pvec = ParamVector(objId,  [pHRatio]);
    gs = GradRuleSetLinear(objId = objId,
        hMin = hMin, hMax = hMax, probHminV = pProbHmin,
        switches = grSwitches, pvec = pvec);
    @assert validate_grad_set(gs)
    return gs
end

function validate_grad_set(gr :: GradRuleSetLinear)
    isValid = (n_colleges(gr) > 1);
    if any(grad_prob_at_hmin(gr) .>= gr.switches.gradProbMax)
        isValid = false;
        @warn "Grad prob at hMin > gradProbMax in $gr"
    end
    return isValid
end

# Only for colleges that produce graduates.
# Must be less than gradProbMax
function init_prob_hmin(parentId :: ObjectId,
    grSwitches :: GradRuleSwitchesLinear)

    objId = make_child_id(parentId, :probHminV);
    if varies_by_college(grSwitches)
        nc = length(grad_colleges(grSwitches));
    else
        # The same grad prob at hMin for all colleges.
        nc = 1;
    end
    # The Bool argument implies a decreasing vector
    pMax = grSwitches.gradProbMax - 0.05;
    # The Bool argument means that better colleges are harder or easier
    if grSwitches.betterHarder
        slope = :decreasing;
    else
        slope = :increasing;  # should be :nonmonotone ++++++
    end
    p = BoundedVector(objId, ParamVector(objId),
        slope,
        zero(Double), pMax,
        fill(Double(0.5), nc));
    set_pvector!(p; description = ldescription(:probHminV),
        symbol = lsymbol(:probHminV));
    return p
end

function grad_prob_at_hmin(gs :: GradRuleSetLinear, ic :: Integer)
    if can_graduate(gs, ic)
        if varies_by_college(gs)
            # Index into list of colleges that can graduate (also into values of probH1V)
            gIdx = findfirst(grad_colleges(gs) .== ic);
        else
            # There is only one value to retrieve. The same for all colleges.
            gIdx = 1;
        end
        probH1 = ModelParams.values(gs.probHminV, gIdx);
    else
        probH1 = zero(Double);
    end
    return probH1
end

grad_prob_at_hmin(gs :: GradRuleSetLinear) =
    [grad_prob_at_hmin(gs, ic)  for ic = 1 : n_colleges(gs)];

"""
	$(SIGNATURES)

Make grad rule for one college
"""
function make_grad_rule(gs :: GradRuleSetLinear, iCollege :: Integer)
    probHmin = grad_prob_at_hmin(gs, iCollege);
    g = GradRuleLinear(t_first_grad(gs, iCollege),
        probHmin, h_min(gs), h_max(gs), grad_prob_min(gs), grad_prob_max(gs),
        min_courses_to_grad(gs, iCollege),
        use_hc_level(gs));
    @assert validate_gr(g)  "Invalid $g"
    return g
end

# function make_test_grad_rule(switches :: GradRuleSwitchesLinear)
#     return GradRuleLinear(t_first_grad(switches, 3), Double(0.5),
#         h_min(switches), Double(4.0),
#         switches.gradProbMin, switches.gradProbMax, switches.minNcForGrad,
#         use_hc_level(switches))
# end

function validate_gr(g :: GradRuleLinear)
    isValid = true;
    if h_min(g) > h_max(g)
        isValid = false;
        @warn "h_min = $(h_min(g)) > h_max = $(h_max(g))"
    end
    return isValid
end

## Graduation probability at end of t
# One state
function grad_prob(g :: GradRuleLinear, t :: Integer, h :: Double, n :: T1,
    h0 :: Double) where T1 <: Integer

    if can_graduate(g, t)
        if n >= min_courses_to_grad(g, t)
            if use_hc_level(g)
                # At h = 1: g.probH1. Then linear in up to g.hMax.
                dh = h - h_min(g);
                dhMax = h_max(g) - h_min(g);
            else
                # At dh = 0: g.probH1. Then linear in dh up to g.hMax.
                dh = h - h0 - h_min(g);
                dhMax = h_max(g) - h_min(g);
            end
            @assert dhMax > 0.0  "Negative dhMax: $dhMax in $g"
            dProb = g.gradProbMax - grad_prob_hmin(g);
            @assert dProb > 0.01  "Invalid dProb: $dProb"
            gProb = grad_prob_hmin(g) + dh / dhMax * dProb;
            gProb = min(g.gradProbMax, max(g.gradProbMin, gProb));
        else
            gProb = zero(Double);
        end
    else
        gProb = zero(Double);
    end
    return gProb
end

# -----------------


text:
function flatten(weights::Params)::Array{Float32, 1}
    cat([ndims(w) > 1 ? reshape(w, length(w)) : w for w in weights]..., dims=1)
end

function unflatten_set(directions::Array{Float32, 1}, weights::Params)::Params
    unflattened = deepcopy(weights)
    idx = 1
    for uf in unflattened
        uf[:] = reshape(directions[idx:idx + length(uf) - 1], size(uf))
        idx += length(uf)
    end
    @assert idx - 1 == length(directions) "Failed to correctly unflatten, please report this issue."
    unflattened
end

function pca_directions(
    weights::Params, checkpoints::Vector{String},
)::Array{Float32, 2}
    directions_matrix = []
    for checkpoint in checkpoints
        @load checkpoint checkpoint_weights
        directions = create_directions(weights, params(checkpoint_weights)) |> flatten
        push!(directions_matrix, directions)
    end
    directions_matrix = cat(directions_matrix..., dims=2)
    model = fit(PCA, directions_matrix, maxoutdim=2)
    projection(model)
end

function project(
    directions::Array{Float32, 1}, pca_projection::Array{Float32, 2},
    method::Symbol = :cos,
)::Array{Float32, 1}
    if method == :lstsq
        return pca_projection \ directions
    end
    normalizer = [norm(pca_projection[:, i]) for i in 1:size(pca_projection, 2)]
    (pca_projection' * directions) ./ normalizer
end

function trajectory(
    weights::Params, pca_projection::Array{Float32, 2},
    checkpoints::Vector{String}, projection::Symbol = :cos,
)::Array{Float32, 2}
    coordinates = zeros(Float32, length(checkpoints), 2)
    @inbounds for (i, checkpoint) in enumerate(checkpoints)
        @load checkpoint checkpoint_weights
        directions = create_directions(weights, params(checkpoint_weights)) |> flatten
        projection = project(directions, pca_projection, :cos)
        coordinates[i, :] = projection
    end
    coordinates
end

"""
```julia
function create_trajectory(
    weights::Params, checkpoints::Vector{String}, projection::Symbol = :cos,
)::Array{Float32, 2}
```

Create trajectory between model checkpoints.
This allows for path visualization that optimizer took
during training, for example.

# Arguments
- `weights::Params`: Target weights of the model. E.g. from last epoch.
- `checkpoints::Vector{String}`: Checkpoints of the model from previous epochs.
They will be loaded in sorted order and should contain weights
under `checkpoint_weights` key.
- `projection::Symbol`: Projection method. Either `:lstsq` or `:cos`.

!!! note
    Checkpoint files should contain weights under `checkpoint_weights` key.
"""
function create_trajectory(
    weights::Params, checkpoints::Vector{String}, projection::Symbol = :cos,
)::Array{Float32, 2}
    directions = pca_directions(weights, checkpoints)
    trajectory(weights, directions, checkpoints, projection)
end


text:
module Day9

using BenchmarkTools
using LoopVectorization
using ..Utils

function find_risk_level(array)
    n, m = size(array)
    risk_level = 0
    for j in 1:m
        running_uphill = false
        i = 1
        while i <= n
            while i < n && array[i, j] >= array[i+1, j] # running downhill or level, keep going
                i += 1
            end

            # Either we hit a rise or the end of the row.
            # Check if the other directions are satisfied.
            if (j == 1 || array[i, j] < array[i, j-1]) && (j == m || array[i, j] < array[i, j+1])
                risk_level += array[i, j] + 1
            end

            i += 1

            while 1 < i <= n && array[i-1, j] <= array[i, j] # running uphill, keep going
                i += 1
            end
        end
    end

    risk_level
end

function find_basin_sizes_bfs(array)
    n, m = size(array)
    basin_sizes = []
    # Queue of unassigned index pairs
    unassigned = Base.Iterators.product(1:n, 1:m) |> Set
    while !isempty(unassigned)
        next = pop!(unassigned)
        # If it's a 9, then continue to the next.
        k = array[next...]
        if k == 9
            continue
        else
            # Otherwise, it's the start of a new basin
            current_basin_queue = [next]
            current_basin_size = 1
            current_basin = Set(current_basin_queue)
            while !isempty(current_basin_queue)
                current_idx = pop!(current_basin_queue)
                i, j = current_idx
                for neighbor in [(i+1, j), (i-1, j), (i, j-1), (i, j+1)]
                    if neighbor ∈ current_basin
                        continue
                    end
                    k, l = neighbor
                    # Check the neighbor is inbounds
                    if !(1 <= k <= n && 1 <= l <= m)
                        continue
                    end
                    if array[k, l] == 9
                        continue
                    else
                        delete!(unassigned, neighbor)
                        push!(current_basin, neighbor)
                        push!(current_basin_queue, neighbor)
                        current_basin_size += 1
                    end
                end
            end
            push!(basin_sizes, current_basin_size)
        end
    end
    basin_sizes
end

function find_basins_coalesce(array)
    n, m = size(array)
    # Start with an array of tuples of coordinates
    basin_roots = Base.Iterators.product(1:n, 1:m) |> collect

    # Collect each row together
    contiguous_ranges = []
    for j in 1:m
        col_ranges = []
        range_start = 1
        for i in 1:n
            k = array[i, j]
            if k == 9
                if range_start < i
                    push!(col_ranges, range_start:i-1)
                end
                range_start = i+1
                basin_roots[i, j] = (0, 0)
            elseif i > 1 && array[i-1, j] != 9
                basin_roots[i, j] = basin_roots[i-1, j]
            end
        end
        if range_start <= n
            push!(col_ranges, range_start:n)
        end
        push!(contiguous_ranges, col_ranges)
    end

    # display(basin_roots)
    # display(contiguous_ranges)
end

function run()
    test = """2199943210
3987894921
9856789892
8767896789
9899965678
""";
    test_array = parse_array(test)
    @assert find_risk_level(test_array) == 15

    input = read("resources/day9_input.txt", String)
    input_array = parse_array(input)
    @show find_risk_level(input_array)

    top3_prod(a) = sort(a)[end-2:end] |> prod
    @assert find_basin_sizes_bfs(test_array) |> top3_prod == 1134
    # @benchmark (find_basin_sizes_bfs($input_array) |> $top3_prod) samples=100
    @time find_basin_sizes_bfs(input_array) |> top3_prod

    @time find_basins_coalesce(input_array)
end

end


text:
module test_julia_regex

using Test

s = """
hello(HELLO)
world(WORLD)
"""

re = r"""<A HREF="(?<page>[\w.]*)">\w*(?<inner>[\w()]*)"""

matches = collect(eachmatch(re, s))

@test matches[1][:page] == "a.html"
@test matches[1][:inner] == "(HELLO)"
@test matches[2][:page] == "b.html"
@test matches[2][:inner] == "(WORLD)"

end # module test_julia_regex


text:
# ---
# title: 1573. Number of Ways to Split a String
# id: problem1573
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/number-of-ways-to-split-a-string/description/>
# hidden: true
# ---
#
# Given a binary string `s` (a string consisting only of '0's and '1's), we can
# split `s` into 3 **non-empty** strings s1, s2, s3 (s1+ s2+ s3 = s).
#
# Return the number of ways `s` can be split such that the number of characters
# '1' is the same in s1, s2, and s3.
#
# Since the answer may be too large, return it modulo 10^9 + 7.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "10101"
#     Output: 4
#     Explanation: There are four ways to split s in 3 parts where each part contain the same number of letters '1'.
#     "1|010|1"
#     "1|01|01"
#     "10|10|1"
#     "10|1|01"
#
#
# **Example 2:**
#
#
#
#     Input: s = "1001"
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: s = "0000"
#     Output: 3
#     Explanation: There are three ways to split s in 3 parts.
#     "0|0|00"
#     "0|00|0"
#     "00|0|0"
#
#
# **Example 4:**
#
#
#
#     Input: s = "100100010100110"
#     Output: 12
#
#
#
#
# **Constraints:**
#
#   * `3 <= s.length <= 10^5`
#   * `s[i]` is `'0'` or `'1'`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Autogenerated wrapper script for PPL_jll for x86_64-apple-darwin
export libppl, libppl_c, ppl_config, ppl_lcdd, ppl_pips

using GMP_jll
JLLWrappers.@generate_wrapper_header("PPL")
JLLWrappers.@declare_library_product(libppl, "@rpath/libppl.14.dylib")
JLLWrappers.@declare_library_product(libppl_c, "@rpath/libppl_c.4.dylib")
JLLWrappers.@declare_executable_product(ppl_config)
JLLWrappers.@declare_executable_product(ppl_lcdd)
JLLWrappers.@declare_executable_product(ppl_pips)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libppl,
        "lib/libppl.14.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libppl_c,
        "lib/libppl_c.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        ppl_config,
        "bin/ppl-config",
    )

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

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using NBInclude

@nbinclude("GroundState(QOptics).ipynb")

using ITensors

using PyPlot

function HHeisdmrg(N, Jxy, Jz, Nsweep, Bonddim, truncerr)
  sites = siteinds("S=1/2",N)

  ampo = OpSum()
  for j=1:N-1
    ampo += Jz , "Sz",j,"Sz",j+1
    ampo += Jxy, "Sx",j,"Sx",j+1
    ampo += Jxy, "Sy",j,"Sy",j+1
  end
H = MPO(ampo,sites)
psi0 = randomMPS(sites,10)

  sweeps = Sweeps(Nsweep)
  setmaxdim!(sweeps, Bonddim)
  setcutoff!(sweeps, truncerr)

  energy, psi = dmrg(H,psi0, sweeps)

  return energy, psi
end

Jxy=1
Jz=1
truncerr=1E-12
nsweep = 5
Bonddim= 5
bondlist = 1:2:10
e0s = []

Nlist = 7:10;
E0sIT=[]
for l in Nlist
    e0, psi = HHeisdmrg(l, Jxy, Jz, nsweep, Bonddim, truncerr)
    push!(E0sIT,real(e0)/l)
end
E0p=[E0sIT[i][1] for i=1:length(Nlist)];

plot(Nlist,E0pQO,color="blue",linestyle = ":",linewidth = 3)
plot(Nlist,E0p,color="red",linewidth = 3)
scatter(Nlist,E0pQO,facecolor = "blue", edgecolors="blue", s=35)
scatter(Nlist,E0p,facecolor = "red", edgecolors="red", s=35)
legend(["ED Data(QOptics)", "ED Data(ITensors)","E0(QOptics)", "E0(ITensors)"])
xlabel("N")
ylabel("E0/N")
title("QOptics/ITensors-GS Energies(N)")
savefig( "GS(QOIT).pdf", bbox_inches = "tight", pad_inches = 0.1 )

Nlist = 7:10;
E0s=[];
for N in Nlist
    for bonddim in bondlist
        e0, psi = HHeisdmrg(N, Jxy, Jz, nsweep, bonddim, truncerr)
        push!(E0s,real(e0)/N)
    end

end
E0pIT=[E0s[i][1] for i=1:length(Nlist)]

e07=E0s[1:5]
e08=E0s[6:10]
e09=E0s[11:15]
e01=E0s[16:20]

plot(bondlist,e07,color = "black",linewidth = 0.6)
plot(bondlist,e08,color = "green",linewidth = 0.6)
plot(bondlist,e09,color = "blue",linewidth = 0.6)
plot(bondlist,e01,color = "red",linewidth = 0.6)
legend(["E07","E08","E09","E010"])
xlabel("Bond Dim")
ylabel("E0/N")
savefig( "ITensorE0(Bondim).pdf", bbox_inches = "tight", pad_inches = 0.1 )

max_stars_repo_name: BlackBodyRadiation/JupyterTNN

max_issues_repo_name: BlackBodyRadiation/JupyterTNN

max_forks_repo_name: BlackBodyRadiation/JupyterTNN


text:
using Random

struct Flip{T}
    node::T
    population::T
    old_assignment::T
    new_assignment::T
end

function random_flip(graph::GraphIndex{Int64}, plan::Plan{Int64})::Flip{Int64}
    edge_index = rand(plan.cut_edges)
    edge_side = rand(Int64[1, 2])
    @inbounds node = graph.edges[edge_side, edge_index]
    @inbounds node_pop = graph.population[node]
    @inbounds old_assignment = plan.assignment[node]
    @inbounds adj_node = graph.edges[3 - edge_side, edge_index]
    @inbounds new_assignment = plan.assignment[adj_node]
    return Flip{Int64}(node, node_pop, old_assignment, new_assignment)
end

max_stars_repo_name: pjrule/gridwalk

max_issues_repo_name: pjrule/gridwalk

max_forks_repo_name: pjrule/gridwalk


text:
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.

"""
    struct FameRange … end

Range includes a frequency code and two endpoints start and end. This is a type
used by the CHLI library.
"""
struct FameRange
    r_freq::Cint
    r_start::Clonglong
    r_end::Clonglong
end
FameRange() = FameRange(0, 0, 0)
# Base.length(r::FameRange) = (-1 ∈ (r.r_end, r.r_start)) ? 0 : Int64(r.r_end)-Int64(r.r_start)+1

@inline unsafe_read!!(::Val{:precision}, key, name, range, vals) =
    (@fame_call_check(fame_get_precisions,
        (Cint, Cstring, Ptr{FameRange}, Ref{Float64}),
        key, name, range, vals);
    vals)

@inline unsafe_read!!(::Val{:numeric}, key, name, range, vals) =
    (@fame_call_check(fame_get_numerics,
        (Cint, Cstring, Ptr{FameRange}, Ref{Float32}),
        key, name, range, vals);
    vals)

@inline unsafe_read!!(::Val{:boolean}, key, name, range, vals) =
    (@fame_call_check(fame_get_booleans,
        (Cint, Cstring, Ptr{FameRange}, Ref{Int32}),
        key, name, range, vals);
    vals)

@inline function unsafe_read!!(::Val{:string}, key, name, range, vals::Ref{String})
    vals[] = unsafe_read!!(Val(:string), key, name, range, String[""])[1]
end

function unsafe_read!!(::Val{:string}, key, name, range, vals::AbstractVector{String})
    lens = Vector{Cint}(undef, length(vals))
    @fame_call_check(fame_len_strings,
        (Cint, Cstring, Ptr{FameRange}, Ref{Int32}),
        key, name, range, lens)
    for ind in eachindex(lens)
        vals[ind] = repeat("\0", lens[ind])
    end
    @fame_call_check(fame_get_strings,
        (Cint, Cstring, Ptr{FameRange}, Ptr{Ptr{UInt8}}, Ref{Int32}, Ptr{Int32}),
        key, name, range, vals, lens, C_NULL)
    return vals
end

# the case where values are dates and so the type is a frequency
@inline function unsafe_read!!(::Val{TY}, key, name, range, vals) where {TY}
    @assert(haskey(fame_freq, TY), "Cannot read data of type $TY.")
    @fame_call_check(fame_get_dates,
        (Cint, Cstring, Ptr{FameRange}, Ref{FameDate}),
        key, name, range, vals)
    vals
end

# the case of namelist.  It's always a scalar
const HNLALL = Int32(-1)
function unsafe_read!!(::Val{:namelist}, key, name, range, vals)
    len = Ref{Cint}(-1)
    @cfm_call_check(cfmnlen,
        (Cint, Cstring, Cint, Ref{Cint}),
        key, name, HNLALL, len)
    # len[] += 1
    vals[] = repeat("\0", len[])
    @cfm_call_check(cfmgtnl,
        (Cint, Cstring, Cint, Ptr{UInt8}, Cint, Ref{Cint}),
        key, name, HNLALL, vals[], len[], len)
    return vals
end

"""
    do_read!(obj::FameObject, db::FameDatabase)

Read data for the given [`FameObject`](@ref) from the given
[`FameDatabase`](@ref). `obj` is modified in place and returned.

A [`FameObject`](@ref) can be created directly, or it could be
returned by [`quick_info`](@ref) or [`listdb`](@ref).
"""
function do_read! end
export do_read!

function do_read!(obj::FameObject{:scalar,FT,FR}, db::FameDatabase) where {FT,FR}
    # number of observations = 0
    # range = CNULL
    unsafe_read!!(Val(FT), db.key, obj.name, C_NULL, obj.data)
    return obj
end

function do_read!(obj::FameObject{:series,FT,FR}, db::FameDatabase) where {FT,FR}
    if obj.first_index[] == FAME_INDEX_NC || obj.last_index[] == FAME_INDEX_NC
        # empty series
        empty!(obj.data)
        return obj
    end
    range = Ref(FameRange(val_to_int(FR, fame_freq), obj.first_index[], obj.last_index[]))
    resize!(obj.data, obj.last_index[] - obj.first_index[] + 1)
    unsafe_read!!(Val(FT), db.key, obj.name, range, obj.data)
    return obj
end

# const HNMVAL = Int32(0) # Normal value; not missing or magic
# const HNCVAL = Int32(1) # Missing NC - Not Computable
# const HNAVAL = Int32(2) # Missing NA - Not Available
# const HNDVAL = Int32(3) # Missing ND - Not Defined
# const HMGVAL = Int32(4) # Magic value; for internal use only


text:
# Method and method-table pretty-printing

function argtype_decl(n, t) # -> (argname, argtype)
    if isa(n,Expr)
        n = n.args[1]  # handle n::T in arg list
    end
    s = string(n)
    i = search(s,'#')
    if i > 0
        s = s[1:i-1]
    end
    if t === Any && !isempty(s)
        return s, ""
    end
    if t <: Vararg && t !== None && t.parameters[1] === Any
        return string(s, "..."), ""
    end
    return s, string(t)
end

function arg_decl_parts(m::Method)
    tv = m.tvars
    if !isa(tv,Tuple)
        tv = (tv,)
    end
    li = m.func.code
    e = uncompressed_ast(li)
    argnames = e.args[1]
    s = symbol("?")
    decls = [argtype_decl(get(argnames,i,s), m.sig[i]) for i=1:length(m.sig)]
    return tv, decls, li.file, li.line
end

function show(io::IO, m::Method)
    print(io, m.func.code.name)
    tv, decls, file, line = arg_decl_parts(m)
    if !isempty(tv)
        show_delim_array(io, tv, '{', ',', '}', false)
    end
    print(io, "(")
    print_joined(io, [isempty(d[2]) ? d[1] : d[1]*"::"*d[2] for d in decls],
                 ",", ",")
    print(io, ")")
    if line > 0
        print(io, " at ", file, ":", line)
    end
end

function show_method_table(io::IO, mt::MethodTable, max::Int=-1, header::Bool=true)
    name = mt.name
    n = length(mt)
    if header
        m = n==1 ? "method" : "methods"
        print(io, "# $n $m for generic function \"$name\":")
    end
    d = mt.defs
    n = rest = 0
    while !is(d,())
        if max==-1 || n<max || (rest==0 && n==max && d.next === ())
            println(io)
            show(io, d)
            n += 1
        else
            rest += 1
        end
        d = d.next
    end
    if rest > 0
        println(io)
        print(io,"... $rest methods not shown (use methods($name) to see them all)")
    end
end

show(io::IO, mt::MethodTable) = show_method_table(io, mt)

inbase(m::Module) = m == Base ? true : m == Main ? false : inbase(module_parent(m))
fileurl(file) = let f = find_source_file(file); f == nothing ? "" : "file://"*f; end
function url(m::Method)
    M = m.func.code.module
    (m.func.code.file == :null || m.func.code.file == :string) && return ""
    file = string(m.func.code.file)
    line = m.func.code.line
    line <= 0 || ismatch(r"In\[[0-9]+\]", file) && return ""
    if inbase(M)
        return "https://github.com/JuliaLang/julia/tree/$(Base.GIT_VERSION_INFO.commit)/base/$file#L$line"
    else
        try
            d = dirname(file)
            u = Git.readchomp(`config remote.origin.url`, dir=d)
            u = match(Git.GITHUB_REGEX,u).captures[1]
            root = cd(d) do # dir=d confuses --show-toplevel, apparently
                Git.readchomp(`rev-parse --show-toplevel`)
            end
            if beginswith(file, root)
                commit = Git.readchomp(`rev-parse HEAD`, dir=d)
                return "https://github.com/$u/tree/$commit/"*file[length(root)+2:end]*"#L$line"
            else
                return fileurl(file)
            end
        catch
            return fileurl(file)
        end
    end
end

function writemime(io::IO, ::MIME"text/html", m::Method)
    print(io, m.func.code.name)
    tv, decls, file, line = arg_decl_parts(m)
    if !isempty(tv)
        print(io,"")
        show_delim_array(io, tv, '{', ',', '}', false)
        print(io,"")
    end
    print(io, "(")
    print_joined(io, [isempty(d[2]) ? d[1] : d[1]*"::"*d[2]*""
                      for d in decls], ",", ",")
    print(io, ")")
    if line > 0
        u = url(m)
        if isempty(u)
            print(io, " at ", file, ":", line)
        else
            print(io, """ at """,
                  file, ":", line, "")
        end
    end
end

function writemime(io::IO, mime::MIME"text/html", mt::MethodTable)
    name = mt.name
    n = length(mt)
    meths = n==1 ? "method" : "methods"
    print(io, "$n $meths for generic function $name:
")
    d = mt.defs
    while !is(d,())
        print(io, "
 ")
        writemime(io, mime, d)
        d = d.next
    end
    print(io, "
")
end

# pretty-printing of Vector{Method} for output of methodswith:

function writemime(io::IO, mime::MIME"text/html", mt::AbstractVector{Method})
    print(io, summary(mt))
    if !isempty(mt)
        print(io, ":
")
        for d in mt
            print(io, "
 ")
            writemime(io, mime, d)
        end
        print(io, "
")
    end
end

# override usual show method for Vector{Method}: don't abbreviate long lists
writemime(io::IO, mime::MIME"text/plain", mt::AbstractVector{Method}) =
    showarray(io, mt, limit=false)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
export addprocs_scyld, ScyldManager

struct ScyldManager <: ClusterManager
    np::Integer
end

function launch(manager::ScyldManager, params::Dict, instances_arr::Array, c::Condition)
    try
        dir = params[:dir]
        exename = params[:exename]
        exeflags = params[:exeflags]
        np = manager.np

        beomap_cmd = `bpsh -1 beomap --no-local --np $np`
        out,beomap_proc = open(beomap_cmd)
        wait(beomap_proc)
        if !success(beomap_proc)
            error("node availability inaccessible (could not run beomap)")
        end
        nodes = split(chomp(readline(out)),':')
        for (i,node) in enumerate(nodes)
            cmd = `cd $dir '&&' $exename $exeflags $(worker_arg())`
            cmd = detach(`bpsh $node sh -l -c $(Base.shell_escape(cmd))`)
            config = WorkerConfig()

            config.io,_ = open(cmd)
            config.io.line_buffered = true
            config.userdata = Dict{Symbol, Any}()
            config.userdata[:node] = node

            push!(instances_arr, config)
            notify(c)
        end
   catch e
        println("Error launching beomap")
        println(e)
   end
end

function manage(manager::ScyldManager, id::Integer, config::WorkerConfig, op::Symbol)
    if op == :interrupt
        if !isnull(config.ospid)
            node = config.userdata[:node]
            if !success(`bpsh $node kill -2 $(get(config.ospid))`)
                println("Error sending Ctrl-C to julia worker $id on node $node")
            end
        else
            # This state can happen immediately after an addprocs
            println("Worker $id cannot be presently interrupted.")
        end
    end
end

addprocs_scyld(np::Integer) = addprocs(ScyldManager(np))


text:
using BinDeps
using Compat
using Libdl  # For the Windows HACK

@BinDeps.setup

libSDL2 = library_dependency("libSDL2", aliases = ["sdl2", "libsdl2-2.0", "libsdl2-2.0-0", "libSDL","SDL2","libSDL2-2.0.so.0","libSDL2-2.0.so.0.4.0"])
if !Sys.iswindows()
    # HACK: These definitions must come later for windows.
    libSDL2_ttf = library_dependency("libSDL2_ttf", aliases = ["SDL_ttf","SDL2_ttf", "libsdl2-ttf-2.0-0","libSDL2_ttf-2.0.so.0"])
    libSDL2_mixer = library_dependency("libSDL2_mixer", aliases = ["SDL_mixer","SDL2_mixer", "libsdl2-mixer-2.0-0","libSDL2_mixer-2.0.so.0"])
end

if Sys.isapple()
    using Homebrew
    provides(Homebrew.HB, "sdl2", libSDL2, os = :Darwin)
    provides(Homebrew.HB, "SDL2_ttf", libSDL2_ttf, os = :Darwin)
    provides(Homebrew.HB, "SDL2_mixer", libSDL2_mixer, os = :Darwin)
end

if Sys.islinux()
    provides(AptGet, "libsdl2-2.0-0", libSDL2)
    provides(AptGet, "libsdl2-ttf-2.0-0",libSDL2_ttf)
    provides(AptGet, "libsdl2-mixer-2.0-0",libSDL2_mixer)
    provides(Yum, "SDL2", libSDL2)
    provides(Pacman, "sdl2", libSDL2)
    provides(Pacman, "sdl2_ttf", libSDL2_ttf)
    provides(Pacman, "sdl2_mixer", libSDL2_mixer)
end

if Sys.iswindows()
    # HACK: First, install just libSDL2, so that the other libs have access to it.
    provides(Binaries, URI("https://www.libsdl.org/release/SDL2-2.0.7-win32-x64.zip"), libSDL2, unpacked_dir=".", os = :Windows)
    @BinDeps.install Dict(
        "libSDL2" => "libSDL2",
    )
    # Now that SDL2 is installed, we can install the extension libraries.
    libSDL2_ttf = library_dependency("libSDL2_ttf", aliases = ["SDL_ttf","SDL2_ttf"])
    libSDL2_mixer = library_dependency("libSDL2_mixer", aliases = ["SDL_mixer","SDL2_mixer"])
    # HACK: In order to get them to pass validation, we have to manually load
    # SDL2.dll, since SDL2_ttf and SDL2_mixer won't dlopen unless SDL2 is
    # already opened. As far as I can tell, there's no way to tell BinDeps about
    # that kind of hard inter-library dependency.
    Libdl.dlopen(joinpath(@__DIR__, "libSDL2", "SDL2.dll"))
    provides(Binaries, URI("https://www.libsdl.org/projects/SDL_ttf/release/SDL2_ttf-2.0.14-win32-x64.zip"), libSDL2_ttf, unpacked_dir=".", os = :Windows)
    provides(Binaries, URI("https://www.libsdl.org/projects/SDL_mixer/release/SDL2_mixer-2.0.2-win32-x64.zip"), libSDL2_mixer, unpacked_dir=".", os = :Windows)
    # Regenerate the deps.jl file with all three libs.
    @BinDeps.install Dict(
        "libSDL2" => "libSDL2",
        "libSDL2_ttf" => "libSDL2_ttf",
        "libSDL2_mixer" => "libSDL2_mixer",
    )
    # HACK: Finally, hilariously, reorder the @checked_lib lines to make sure
    # SDL2 is loaded first at runtime.
    deps_jl = joinpath(@__DIR__, "deps.jl")
    depslines = readlines(deps_jl)
    liblines = [occursin(r"SDL2", line) for line in depslines]
    depslines[liblines] = sort(depslines[liblines])
    write(deps_jl, join(depslines, "\n"))
else
    @BinDeps.install Dict(
        "libSDL2" => "libSDL2",
        "libSDL2_ttf" => "libSDL2_ttf",
        "libSDL2_mixer" => "libSDL2_mixer",
    )
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using RadMod2D

include("./models2D.jl")

# test parameter
d1 = 0.8
d2 = 1.6
elemsize = 0.005
n = 15
criteria = 1E-5

# analytical
r1 = d1 / 2
r2 = d2 / 2
vfmata = zeros(Float64, 2, 2)
vfmata[1,1] = 0
vfmata[1,2] = 1
vfmata[2,1] = r1 / r2
vfmata[2,2] = 1- vfmata[2,1]

# RadMod2D
m = model_circle_in_circle(0.8, 1.6, elemsize)
vfmat = zeros(Float64, m.nelem, m.nelem)
@time existing_vf!(m, vfmat)
dx, dy = create_tiles(m, n)
@time t_occ = check_tile_occupation(m, dx, dy, n)
@time blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
@time calculating_vf!(m, vfmat, normit = false)
vfmatp = vfmat_to_parts(m, vfmat, normit = false)

# compare results to check if algorithm is working properly
println("analytical results:")
for i = 1:size(vfmata,1)
    println(vfmata[i,:])
end
println("RadMod2D results:")
for i = 1:size(vfmatp,1)
    println(vfmatp[i,:])
end
diff = vfmata - vfmatp
# println("difference:")
# for i = 1:size(diff,1)
#     println(diff[i,:])
# end
for i = 1:4
    if diff[i] > criteria
        error("Radmod2D not working properly")
    end
end

# include("./test/run_test_vf.jl")


text:
module Input

using StaticArrays
using UnPack

using GLFW: GLFW, Joystick as JoystickSlot, DeviceConfigEvent, JoystickPresent,
        GetJoystickAxes, GetJoystickButtons, GetJoystickName

export connected_joysticks, init_joysticks, update_joystick, joystick_callback
export get_axis_data, get_button_state, get_button_change
export ButtonChange, button_unchanged, button_pressed, button_released
export exp_axis_curve
export XBoxController

abstract type AbstractJoystick end

const connected_joysticks = Dict{JoystickSlot, AbstractJoystick}()

###################

@enum ButtonChange begin
    button_unchanged = 0
    button_pressed = 1
    button_released = 2
end
Base.convert(::Type{ButtonChange}, n::Integer) = ButtonChange(n)
Base.zero(::Type{ButtonChange}) = button_unchanged

######################### XBox Controller ########################

const XBoxAxisLabels = (
    :left_analog_x, :left_analog_y, :right_analog_x, :right_analog_y,
    :left_trigger, :right_trigger
)
const XBoxButtonLabels = (
    :button_A, :button_B, :button_X, :button_Y, :left_bumper, :right_bumper,
    :view, :menu, :left_analog, :right_analog,
    :dpad_up, :dpad_right, :dpad_down, :dpad_left
)

const default_axis_mapping = NamedTuple{XBoxAxisLabels}(1:6)
const default_button_mapping = NamedTuple{XBoxButtonLabels}(1:14)

Base.@kwdef struct XBoxAxes
    mapping::NamedTuple{XBoxAxisLabels, NTuple{6,Int}} = default_axis_mapping
    data::MVector{6,Float32} = zeros(MVector{6, Float32})
end

function Base.getindex(axes::XBoxAxes, s::Symbol)
    data = getproperty(axes, :data); mapping = getproperty(axes, :mapping)
    data[mapping[s]]
end

function Base.setindex!(axes::XBoxAxes, v, s::Symbol)
    data = getproperty(axes, :data); mapping = getproperty(axes, :mapping)
    data[mapping[s]] = v
end

Base.@kwdef struct XBoxButtons
    mapping::NamedTuple{XBoxButtonLabels, NTuple{14,Int}} = default_button_mapping
    state::MVector{14,Bool} = zeros(MVector{14,Bool})
    change::MVector{14,ButtonChange} = zeros(MVector{14,ButtonChange})
end

Base.@kwdef struct XBoxController <: AbstractJoystick
    slot::JoystickSlot = GLFW.JOYSTICK_1
    axes::XBoxAxes = XBoxAxes()
    buttons::XBoxButtons = XBoxButtons()
end

XBoxController(slot::JoystickSlot) = XBoxController(; slot)

function update_joystick(joystick::XBoxController)

    @unpack slot, axes, buttons = joystick

    axes.data .= GetJoystickAxes(slot)
    axes[:left_trigger] = 0.5*(1 + axes[:left_trigger])
    axes[:right_trigger] = 0.5*(1 + axes[:right_trigger])

    buttons_state_new = SVector{14,Bool}(GetJoystickButtons(slot))

    for i in 1:length(buttons.state)
        if !buttons.state[i] && buttons_state_new[i]
            buttons.change[i] = button_pressed
        elseif buttons.state[i] && !buttons_state_new[i]
            buttons.change[i] = button_released
        else
            buttons.change[i] = button_unchanged
        end
    end

    buttons.state .= buttons_state_new

end

get_axis_data(joy::XBoxController) = NamedTuple{XBoxAxisLabels}(Tuple(joy.axes.data))
get_button_state(joy::XBoxController) = NamedTuple{XBoxButtonLabels}(Tuple(joy.buttons.state))
get_button_change(joy::XBoxController) = NamedTuple{XBoxButtonLabels}(Tuple(joy.buttons.change))

function get_axis_data(joystick::XBoxController, s::Symbol)
    @unpack data, mapping = joystick.axes
    data[mapping[s]]
end

function get_button_state(joystick::XBoxController, s::Symbol)
    @unpack state, mapping = joystick.buttons
    state[mapping[s]]
end

function get_button_change(joystick::XBoxController, s::Symbol)
    @unpack change, mapping = joystick.buttons
    change[mapping[s]]
end

##################### Joystick configuration ################

#adds any joysticks already connected when GLFW is first imported. from this
#moment, connections and disconnections will need to be caught by the callback
#by calling GLFW.PollEvents()
function init_joysticks()

    for slot in instances(JoystickSlot)
        delete!(connected_joysticks, slot)
        if JoystickPresent(slot)
            add_joystick(slot)
        end
    end

end

function joystick_callback(slot::JoystickSlot, event::DeviceConfigEvent)

    if event === GLFW.CONNECTED
        add_joystick(slot)
    elseif event === GLFW.DISCONNECTED
        remove_joystick(slot)
    end

end

function add_joystick(slot::JoystickSlot)

    if !JoystickPresent(slot)
        println("Error while adding device at slot $slot: not found")
        return
    end

    joystick_model = GetJoystickName(slot)

    if joystick_model === "Xbox Controller"
        joystick = XBoxController(slot)
        println("$joystick_model now active at slot $slot")
    else
        error("Interface for joystick $joystick_model not implemented")
    end

    connected_joysticks[slot] = joystick

end

function remove_joystick(slot::JoystickSlot)

    if JoystickPresent(slot)
        joystick_model = GetJoystickName(slot)
        println("Can't remove $joystick_model from slot $slot: device still connected")
        return
    end

    delete!(connected_joysticks, slot)
    println("$slot is no longer active")
end

####################### Response Curves ####################
            # f_Oc_xc = cls.exp_curve(-state.lstick_y, strength=1, deadzone=0.3),

function exp_axis_curve(x::Real; strength::Real = 0.0, deadzone::Real = 0.0)

    a = strength
    x0 = deadzone

    abs(x) <= 1 || throw(ArgumentError("Input to exponential curve must be within [-1, 1]"))
    (x0 >= 0 && x0 <= 1) || throw(ArgumentError("Exponential curve deadzone must be within [0, 1]"))

    if x > 0
        y = max(0, (x - x0)/(1 - x0)) * exp( a * (abs(x) -1) )
    else
        y = min(0, (x + x0)/(1 - x0)) * exp( a * (abs(x) -1) )
    end
end

end #module


text:
#=

@author : Spencer Lyon <spencer.lyon@nyu.edu>
          Victoria Gregory <victoria.gregory@nyu.edu>

=#
using Plots
pyplot()
using Distributions
using LaTeXStrings

# == Set parameters == #
srand(42)  # reproducible results
n = 250    # Choice of n
k = 10000  # Number of draws of Y_n
dist = Exponential(1./2.)  # Exponential distribution, lambda = 1/2
μ, s = mean(dist), std(dist)

# == Draw underlying RVs. Each row contains a draw of X_1,..,X_n == #
data = rand(dist, k, n)

# == Compute mean of each row, producing k draws of \bar X_n == #
sample_means = mean(data, 2)

# == Generate observations of Y_n == #
Y = sqrt(n) * (sample_means .- μ)

# == Plot == #
xmin, xmax = -3 * s, 3 * s
histogram(Y, nbins=60, alpha=0.5, xlims=(xmin, xmax),
          norm=true, label="")
xgrid = linspace(xmin, xmax, 200)
plot!(xgrid, pdf.(Normal(0.0, s), xgrid), color=:black,
      linewidth=2, label=LaTeXString("\$N(0, \\sigma^2=$(s^2))\$"),
      legendfont=font(12))


text:
function recursion_wolff!(grid::Array{Int, 2},         # Spin grig
                          cluster::BitArray{2},        # Cluster grid
                          i::Int,                      # Row coordinate
                          j::Int;                      # Column coordinate
                          h::Float64            = 0.0, # External field
                          temp::Float64         = 1.0) # Temperature

    cluster[i, j] = true
    for (x, y) in neighbors(grid, i, j)
        if cluster[x,y] == false && grid[x,y] == grid[i, j]
            # Add the spin to cluster
            if rand() < 1 - exp(-2/temp) recursion_wolff!(grid, cluster, x, y, h=h, temp=temp) end
        end
    end
end

# Runs the Wolff algorithm at the grid and returns the final magnetization
function wolff!(grid::Array{Int, 2};
                h::Float64           = 0.0,  # External field
                temp::Float64        = 1.0,  # Temperature
                iters::Integer       = 30,   # Number of iterations
                plot::Bool           = true, # Plot flag
                verbose::Bool        = true) # Verbose flag

    m = Float64[]
    for i in 1:iters
        # Randomly pick a position within the grid
        x   = rand(1:size(grid, 1))
        y   = rand(1:size(grid, 2))

        # Create a cluster grid full of "falses"
        cluster = falses(size(grid))

        # Run the algorithm and calculate the final magnetization
        recursion_wolff!(grid, cluster, x, y, h=h, temp=temp)
        # Calculate ΔE
        ΔE = -2h*clusterspin(grid, cluster)
        # Change spin accordingly
	    if ΔE <= 0 || rand() < exp(-ΔE) flip!(grid, cluster) end
        push!(m, magnetization(grid))
    end

    if verbose println("(T=$temp) Wolff ended with magnetization $(m[end]) after $iters iterations") end
    if plot
        PyPlot.plot(1:length(m), m, "o", color="blue")
        PyPlot.plot(1:length(m), m, "-", color="blue")
        PyPlot.title("Wolff for T=$temp for size $size ($(length(m)) iterations and H=$h)")
        PyPlot.xlabel("Number of Iterations")
        PyPlot.ylabel("Magnetization")
        PyPlot.ylim(0,1.1)
        PyPlot.savefig("wolff_T=$temp.png")
        PyPlot.close()
    end

    return m
end


text:
# Uniforms are OpenGL variables that stay the same for the entirety of a drawcall.
# There are a lot of functions, to upload them, as OpenGL doesn't rely on multiple dispatch.
# here is my approach, to handle all of the uniforms with one function, namely gluniform
# For uniforms, the Vector and Matrix types from ImmutableArrays should be used, as they map the relation almost 1:1

const GLSL_COMPATIBLE_NUMBER_TYPES = (GLfloat, GLint, GLuint, GLdouble)
const NATIVE_TYPES = Union{
    StaticArray, GLSL_COMPATIBLE_NUMBER_TYPES...,
    ZeroIndex{GLint}, ZeroIndex{GLuint},
    GLBuffer, GPUArray, Shader, GLProgram, NativeMesh
}

opengl_prefix(T)  = error("Object $T is not a supported uniform element type")
opengl_postfix(T) = error("Object $T is not a supported uniform element type")

opengl_prefix(x::Type{T}) where {T <: Union{FixedPoint, Float32, Float16}} = ""
opengl_prefix(x::Type{T}) where {T <: Float64} = "d"
opengl_prefix(x::Type{Cint}) = "i"
opengl_prefix(x::Type{T}) where {T <: Union{Cuint, UInt8, UInt16}} = "u"

opengl_postfix(x::Type{Float64}) = "dv"
opengl_postfix(x::Type{Float32}) = "fv"
opengl_postfix(x::Type{Cint})    = "iv"
opengl_postfix(x::Type{Cuint})   = "uiv"

function uniformfunc(typ::DataType, dims::Tuple{Int})
    Symbol(string("glUniform", first(dims), opengl_postfix(typ)))
end
function uniformfunc(typ::DataType, dims::Tuple{Int, Int})
    M, N = dims
    Symbol(string("glUniformMatrix", M == N ? "$M" : "$(M)x$(N)", opengl_postfix(typ)))
end

function gluniform(location::Integer, x::FSA) where FSA <: Union{StaticArray, Colorant}
    xref = [x]
    gluniform(location, xref)
end

_size(p) = size(p)
_size(p::Colorant) = (length(p),)
_size(p::Type{T}) where {T <: Colorant} = (length(p),)
_ndims(p) = ndims(p)
_ndims(p::Type{T}) where {T <: Colorant} = 1

@generated function gluniform(location::Integer, x::Vector{FSA}) where FSA <: Union{StaticArray, Colorant}
    func = uniformfunc(eltype(FSA), _size(FSA))
    callexpr = if _ndims(FSA) == 2
        :($func(location, length(x), GL_FALSE, x))
    else
        :($func(location, length(x), x))
    end
    quote
        $callexpr
    end
end

#Some additional uniform functions, not related to Imutable Arrays
gluniform(location::Integer, target::Integer, t::Texture) = gluniform(GLint(location), GLint(target), t)
gluniform(location::Integer, target::Integer, t::GPUVector) = gluniform(GLint(location), GLint(target), t.buffer)
gluniform(location::Integer, target::Integer, t::Observable) = gluniform(GLint(location), GLint(target), to_value(t))
gluniform(location::Integer, target::Integer, t::TextureBuffer) = gluniform(GLint(location), GLint(target), t.texture)
function gluniform(location::GLint, target::GLint, t::Texture)
    activeTarget = GL_TEXTURE0 + UInt32(target)
    glActiveTexture(activeTarget)
    glBindTexture(t.texturetype, t.id)
    gluniform(location, target)
end
gluniform(location::Integer, x::Enum) = gluniform(GLint(location), GLint(x))

function gluniform(loc::Integer, x::Observable{T}) where T
    gluniform(GLint(loc), to_value(x))
end

gluniform(location::Integer, x::Union{GLubyte, GLushort, GLuint}) = glUniform1ui(GLint(location), x)
gluniform(location::Integer, x::Union{GLbyte, GLshort, GLint, Bool}) = glUniform1i(GLint(location),  x)
gluniform(location::Integer, x::GLfloat) = glUniform1f(GLint(location),  x)
gluniform(location::Integer, x::GLdouble) = glUniform1d(GLint(location),  x)

#Uniform upload functions for julia arrays...
gluniform(location::GLint, x::Vector{Float32}) = glUniform1fv(location,  length(x), x)
gluniform(location::GLint, x::Vector{GLdouble}) = glUniform1dv(location,  length(x), x)
gluniform(location::GLint, x::Vector{GLint}) = glUniform1iv(location,  length(x), x)
gluniform(location::GLint, x::Vector{GLuint}) = glUniform1uiv(location, length(x), x)

glsl_typename(x::T) where {T} = glsl_typename(T)
glsl_typename(t::DataType) = error("Datatype $(t) not supported")
glsl_typename(t::Type{Nothing}) = "Nothing"
glsl_typename(t::Type{GLfloat}) = "float"
glsl_typename(t::Type{GLdouble}) = "double"
glsl_typename(t::Type{GLuint}) = "uint"
glsl_typename(t::Type{GLint}) = "int"
glsl_typename(t::Type{T}) where {T <: Union{StaticVector, Colorant}} = string(opengl_prefix(eltype(T)), "vec", length(T))
glsl_typename(t::Type{TextureBuffer{T}}) where {T} = string(opengl_prefix(eltype(T)), "samplerBuffer")

function glsl_typename(t::Texture{T, D}) where {T, D}
    str = string(opengl_prefix(eltype(T)), "sampler", D, "D")
    t.texturetype == GL_TEXTURE_2D_ARRAY && (str *= "Array")
    str
end

function glsl_typename(t::Type{T}) where T <: SMatrix
    M, N = size(t)
    string(opengl_prefix(eltype(t)), "mat", M==N ? M : string(M, "x", N))
end
toglsltype_string(t::Observable) = toglsltype_string(to_value(t))
toglsltype_string(x::T) where {T<:Union{Real, StaticArray, Texture, Colorant, TextureBuffer, Nothing}} = "uniform $(glsl_typename(x))"
#Handle GLSL structs, which need to be addressed via single fields
function toglsltype_string(x::T) where T
    if isa_gl_struct(x)
        string("uniform ", T.name.name)
    else
        error("can't splice $T into an OpenGL shader. Make sure all fields are of a concrete type and isbits(FieldType)-->true")
    end
end
toglsltype_string(t::Union{GLBuffer{T}, GPUVector{T}}) where {T} = string("in ", glsl_typename(T))
# Gets used to access a
function glsl_variable_access(keystring, t::Texture{T, D}) where {T,D}
    fields = SubString("rgba", 1, length(T))
    if t.texturetype == GL_TEXTURE_BUFFER
        return string("texelFetch(", keystring, "index).", fields, ";")
    end
    return string("getindex(", keystring, "index).", fields, ";")
end
function glsl_variable_access(keystring, ::Union{Real, GLBuffer, GPUVector, StaticArray, Colorant})
    string(keystring, ";")
end
function glsl_variable_access(keystring, s::Observable)
    glsl_variable_access(keystring, to_value(s))
end
function glsl_variable_access(keystring, t::Any)
    error("no glsl variable calculation available for : ", keystring, " of type ", typeof(t))
end

function uniform_name_type(program::GLuint)
    uniformLength = glGetProgramiv(program, GL_ACTIVE_UNIFORMS)
    Dict{Symbol, GLenum}(ntuple(uniformLength) do i # take size and name
        name, typ = glGetActiveUniform(program, i-1)
    end)
end
function attribute_name_type(program::GLuint)
    uniformLength = glGetProgramiv(program, GL_ACTIVE_ATTRIBUTES)
    Dict{Symbol, GLenum}(ntuple(uniformLength) do i
    	name, typ = glGetActiveAttrib(program, i-1)
    end)
end
function istexturesampler(typ::GLenum)
    return (
        typ == GL_SAMPLER_BUFFER || typ == GL_INT_SAMPLER_BUFFER || typ == GL_UNSIGNED_INT_SAMPLER_BUFFER ||
    	typ == GL_IMAGE_2D ||
        typ == GL_SAMPLER_1D || typ == GL_SAMPLER_2D || typ == GL_SAMPLER_3D ||
        typ == GL_UNSIGNED_INT_SAMPLER_1D || typ == GL_UNSIGNED_INT_SAMPLER_2D || typ == GL_UNSIGNED_INT_SAMPLER_3D ||
        typ == GL_INT_SAMPLER_1D || typ == GL_INT_SAMPLER_2D || typ == GL_INT_SAMPLER_3D ||
        typ == GL_SAMPLER_1D_ARRAY || typ == GL_SAMPLER_2D_ARRAY ||
        typ == GL_UNSIGNED_INT_SAMPLER_1D_ARRAY || typ == GL_UNSIGNED_INT_SAMPLER_2D_ARRAY ||
        typ == GL_INT_SAMPLER_1D_ARRAY || typ == GL_INT_SAMPLER_2D_ARRAY
    )
end

gl_promote(x::Type{T}) where {T <: Integer} = Cint
gl_promote(x::Type{Union{Int16, Int8}}) = x

gl_promote(x::Type{T}) where {T <: Unsigned} = Cuint
gl_promote(x::Type{Union{UInt16, UInt8}}) = x

gl_promote(x::Type{T}) where {T <: AbstractFloat} = Float32
gl_promote(x::Type{Float16}) = x

gl_promote(x::Type{T}) where {T <: Normed} = N0f32
gl_promote(x::Type{N0f16}) = x
gl_promote(x::Type{N0f8}) = x

const Color3{T} = Colorant{T, 3}
const Color4{T} = Colorant{T, 4}

gl_promote(x::Type{Bool})                  = GLboolean
gl_promote(x::Type{T}) where {T <: Gray} = Gray{gl_promote(eltype(T))}
gl_promote(x::Type{T}) where {T <: Color3} = RGB{gl_promote(eltype(T))}
gl_promote(x::Type{T}) where {T <: Color4} = RGBA{gl_promote(eltype(T))}
gl_promote(x::Type{T}) where {T <: BGRA} = BGRA{gl_promote(eltype(T))}
gl_promote(x::Type{T}) where {T <: BGR} = BGR{gl_promote(eltype(T))}

gl_promote(x::Type{T}) where {T <: StaticVector} = similar_type(T, gl_promote(eltype(T)))

gl_promote(x::Type{T}) where {T <: GeometryBasics.Mesh} = NativeMesh{T}

gl_convert(x::AbstractVector{Vec3f}) = x

gl_convert(x::T) where {T <: Number} = gl_promote(T)(x)
gl_convert(x::T) where {T <: Colorant} = gl_promote(T)(x)
gl_convert(x::T) where {T <: AbstractMesh} = gl_convert(x)
gl_convert(x::T) where {T <: GeometryBasics.Mesh} = gl_promote(T)(x)
gl_convert(x::Observable{T}) where {T <: GeometryBasics.Mesh} = gl_promote(T)(x)

gl_convert(s::Vector{Matrix{T}}) where {T<:Colorant} = Texture(s)
gl_convert(s::Nothing) = s

isa_gl_struct(x::AbstractArray) = false
isa_gl_struct(x::NATIVE_TYPES) = false
isa_gl_struct(x::Colorant) = false
function isa_gl_struct(x::T) where T
    !isconcretetype(T) && return false
    if T <: Tuple
        return false
    end
    fnames = fieldnames(T)
    !isempty(fnames) && all(name -> isconcretetype(fieldtype(T, name)) && isbits(getfield(x, name)), fnames)
end
function gl_convert_struct(x::T, uniform_name::Symbol) where T
    if isa_gl_struct(x)
        return Dict{Symbol, Any}(map(fieldnames(x)) do name
            (Symbol("$uniform_name.$name") => gl_convert(getfield(x, name)))
        end)
    else
        error("can't convert $x to a OpenGL type. Make sure all fields are of a concrete type and isbits(FieldType)-->true")
    end
end

# native types don't need convert!
gl_convert(a::T) where {T <: NATIVE_TYPES} = a
gl_convert(s::Observable{T}) where {T <: NATIVE_TYPES} = s
gl_convert(s::Observable{T}) where T = const_lift(gl_convert, s)
gl_convert(x::StaticVector{N, T}) where {N, T} = map(gl_promote(T), x)
gl_convert(x::SMatrix{N, M, T}) where {N, M, T} = map(gl_promote(T), x)
gl_convert(a::AbstractVector{<: AbstractFace}) = indexbuffer(s)
gl_convert(t::Type{T}, a::T; kw_args...) where T <: NATIVE_TYPES = a
gl_convert(::Type{<: GPUArray}, a::StaticVector) = gl_convert(a)

function gl_convert(T::Type{<: GPUArray}, a::AbstractArray{X, N}; kw_args...) where {X, N}
    T(convert(AbstractArray{gl_promote(X), N}, a); kw_args...)
end

gl_convert(::Type{<: GLBuffer}, x::GLBuffer; kw_args...) = x
gl_convert(::Type{Texture}, x::Texture) = x
gl_convert(::Type{<: GPUArray}, x::GPUArray) = x

function gl_convert(::Type{T}, a::Vector{Array{X, 2}}; kw_args...) where {T <: Texture, X}
    T(a; kw_args...)
end
gl_convert(::Type{<: GPUArray}, a::Observable{<: StaticVector}) = gl_convert(a)

function gl_convert(::Type{T}, a::Observable{<: AbstractArray{X, N}}; kw_args...) where {T <: GPUArray, X, N}
    TGL = gl_promote(X)
    s = (X == TGL) ? a : lift(x-> convert(Array{TGL, N}, x), a)
    T(s; kw_args...)
end

lift_convert(a::AbstractArray, T, N) = lift(x -> convert(Array{T, N}, x), a)
function lift_convert(a::ShaderAbstractions.Sampler, T, N)
    ShaderAbstractions.Sampler(
        lift(x -> convert(Array{T, N}, x.data), a),
        minfilter = a[].minfilter, magfilter = a[].magfilter,
        x_repeat = a[].repeat[1],
        y_repeat = a[].repeat[min(2, N)],
        z_repeat = a[].repeat[min(3, N)],
        anisotropic = a[].anisotropic, swizzle_mask = a[].swizzle_mask
    )
end

gl_convert(f::Function, a) = f(a)


text:
@doc raw"""
    BiasedSKCE(k)

Biased estimator of the squared kernel calibration error (SKCE) with kernel `k`.

Kernel `k` on the product space of predictions and targets has to be a `Kernel` from the
Julia package
[KernelFunctions.jl](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl)
that can be evaluated for inputs that are tuples of predictions and targets.

# Details

The estimator is biased and guaranteed to be non-negative. Its sample complexity is
``O(n^2)``, where ``n`` is the total number of samples.

Let ``(P_{X_i}, Y_i)_{i=1,\ldots,n}`` be a data set of predictions and corresponding
targets. The estimator is defined as
```math
\frac{1}{n^2} \sum_{i,j=1}^n h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big)
```
where
```math
\begin{aligned}
h_k\big((μ, y), (μ', y')\big) ={}&   k\big((μ, y), (μ', y')\big)
                                   - 𝔼_{Z ∼ μ} k\big((μ, Z), (μ', y')\big) \\
                                 & - 𝔼_{Z' ∼ μ'} k\big((μ, y), (μ', Z')\big)
                                   + 𝔼_{Z ∼ μ, Z' ∼ μ'} k\big((μ, Z), (μ', Z')\big).
\end{aligned}
```

The bias of the estimator is
```math
\frac{1}{n} \Big(\mathbb{E}_{P_X,Y} h_k\big((P_X, Y), (P_X, Y)\big) - \mathrm{SKCE}_k\Big).
```

# References

Widmann, D., Lindsten, F., & Zachariah, D. (2019). [Calibration tests in multi-class
classification: A unifying framework](https://proceedings.neurips.cc/paper/2019/hash/1c336b8080f82bcc2cd2499b4c57261d-Abstract.html).
In: Advances in Neural Information Processing Systems (NeurIPS 2019) (pp. 12257–12267).

Widmann, D., Lindsten, F., & Zachariah, D. (2021). [Calibration tests beyond
classification](https://openreview.net/forum?id=-bxf89v3Nx).

See also: [`UnbiasedSKCE`](@ref), [`BlockUnbiasedSKCE`](@ref)
"""
struct BiasedSKCE{K<:Kernel} <: SKCE
    """Kernel of estimator."""
    kernel::K
end

function (skce::BiasedSKCE)(predictions::AbstractVector, targets::AbstractVector)
    @unpack kernel = skce

    # obtain number of samples
    nsamples = length(predictions)
    nsamples ≥ 1 || error("there must be at least one sample")

    @inbounds begin
        # evaluate kernel function for the first sample
        prediction = predictions[1]
        target = targets[1]
        hij = unsafe_skce_eval(kernel, prediction, target, prediction, target)

        # initialize the calibration error estimate
        estimate = hij / 1

        # for all other pairs of samples
        n = 1
        for i in 2:nsamples
            predictioni = predictions[i]
            targeti = targets[i]

            for j in 1:(i - 1)
                predictionj = predictions[j]
                targetj = targets[j]

                # evaluate the kernel function
                hij = unsafe_skce_eval(kernel, predictioni, targeti, predictionj, targetj)

                # update the estimate (add two terms due to symmetry!)
                n += 2
                estimate += 2 * (hij - estimate) / n
            end

            # evaluate the kernel function
            hij = unsafe_skce_eval(kernel, predictioni, targeti, predictioni, targeti)

            # update the estimate
            n += 1
            estimate += (hij - estimate) / n
        end
    end

    return estimate
end


text:
using Test, DiracNotation, Random

include("array_test.jl")
include("qo_test.jl")


text:
using Test

# A recursive descent parser for a string representation of boolean expressions

mutable struct ExprString
    data::String
    pos::Int
    ExprString(data, pos) = pos != 1 ? error("start must be 1") : new(data, 1)
end

ExprString(s::String) = ExprString(s, 1)

hasmore(e::ExprString) = e.pos <= length(e.data)
isparen(c::Char) = c == '(' || c == ')'
isbinop(c::Char) = c == '&' || c == '|'
isunop(c::Char) = c == '!'
isalphanum(c::Char) = isletter(c) || isnumeric(c)

function skip_when(f, s, start)::Int
    len = length(s)
    if start <= len
        for i in start:len
            if !f(s[i])
                return i
            end
        end
    end
    return len + 1
end

function next_pos(s::String, pos::Int)::Int
    c = s[pos]
    if isparen(c) || isunop(c) || isbinop(c)
        return skip_when(isspace, s, pos + 1)
    elseif isletter(c) #We're inside a variable name, skip alphanumeric
        return skip_when(isalphanum, s, pos + 1)
    else
        throw(ArgumentError("Expected alphanum, parens, \"&\", \"|\", or \"!\", found \"$c\""))
    end
end

function peek(e::ExprString)::String
    c = e.data[e.pos]
    if isparen(c) || isbinop(c)
        return string(c)
    end
    endPos = next_pos(e.data, e.pos)
    strip(e.data[e.pos:endPos-1])
end

function consume(e::ExprString, s::String)
    nextStartPos = next_pos(e.data, e.pos)
    expected = strip(e.data[e.pos:nextStartPos-1])
    if expected == s
        e.pos = nextStartPos
    else
        throw(ArgumentError("Expected $expected, saw $s"))
    end
end

function parse_expr(e::ExprString)::BoolExpr
    term = parse_term(e)
    if hasmore(e)
        tok_string = peek(e)
        if tok_string == "&"
            consume(e, "&")
            return And([term, parse_expr(e)])
        elseif tok_string == "|"
            consume(e, "|")
            return Or([term, parse_expr(e)])
        end
    end
    return term
end

function parse_term(e::ExprString)::BoolExpr
    tok_string = peek(e)
    if tok_string == "!"
        consume(e, "!")
        result = parse_term(e)
        return Not(result)
    elseif tok_string == "("
        consume(e, "(")
        result = parse_expr(e)
        consume(e, ")")
        return result
    elseif true # isAlpha(tok_string)
        consume(e, tok_string)
        return Var(tok_string)
    else
        throw(ArgumentError("Syntax error: expected variable name or \"(\", found \"$tok_string\")"))
    end
end

Base.show(io::IO, exprString::ExprString) = print(io, "ExprString(data=$e.data, pos=$e.pos)")
Base.show(io::IO, boolExpr::BoolExpr) = BoolExprs.show(io, boolExpr)
mkexpr(s) = parse_expr(ExprString(s))

"""
"""
#######

Base.:(==)(x::BoolExpr, y::BoolExpr)::Bool = equiv(x, y)

@testset "Make expressions from strings" begin
    x1 = mkexpr("x1")
    @test x1 == Var("x1")
    x2 = mkexpr("x2")
    @test x2 == Var("x2")
    x3 = mkexpr("x3")
    @test x3 == Var("x3")

    not_x3 = mkexpr("!x3")
    @test not_x3 == Not(Var("x3"))

    and_x1_x2 = mkexpr("x1&x2")
    @test and_x1_x2 == And([x1, x2])

    and_x2_x3 = And([x2, x3])
    @test and_x2_x3 == And([x2, x3])

    or_x1_x2 = mkexpr("x1|x2")
    @test or_x1_x2 == Or([x1, x2])

    or_x2_x3 = mkexpr("x2|x3")
    @test or_x2_x3 == Or([x2, x3])
end
#######

@testset "DNF" begin
    s = "(x1&x2) & (x2&x3)"
    dnf_expr = dnf(mkexpr(s))
    @test dnf_expr == Or([And([Var("x1"), Var("x2"), Var("x2"), Var("x3")])])
    println("dnf($s) = $dnf_expr")
    s = "(x1&x2) | (x2&x3)"
    dnf_expr = dnf(mkexpr(s))
    @test dnf_expr == Or([And([Var("x1"), Var("x2")]), And([Var("x2"), Var("x3")])])
    println("dnf($s) = $dnf_expr")
    s = "(x1|x2) & (x2|x3)"
    dnf_expr = dnf(mkexpr(s))
    @test dnf_expr == Or([And([Var("x1"), Var("x2")]), And([Var("x1"), Var("x3")]), And([Var("x2"), Var("x2")]), And([Var("x2"), Var("x3")])])
    println("dnf($s) = $dnf_expr")
end

@testset "CNF" begin
    s = "(x1&x2) & (x2&x3)"
    cnf_expr = cnf(mkexpr(s))
    @test cnf_expr == And([Var("x1"), Var("x2"), Var("x2"), Var("x3")])
    println("cnf($s) = $cnf_expr")
    s = "(x1&x2) | (x2&x3)"
    cnf_expr = cnf(mkexpr(s))
    @test cnf_expr == And([Or([Var("x1"), Var("x2")]), Or([Var("x1"), Var("x3")]), Or([Var("x2"), Var("x2")]), Or([Var("x2"), Var("x3")])])
    println("cnf($s) = $cnf_expr")
    s = "(x1|x2) & (x2|x3)"
    cnf_expr = cnf(mkexpr(s))
     @test cnf_expr == And([Or([Var("x1"), Var("x2")]), Or([Var("x2"), Var("x3")])])
   println("cnf($s) = $cnf_expr")
end


text:
##########################
# FunctionCell Interface #
##########################
"""
	lda(distances=true [;kwargs])

Constructs an untrained cell that when piped data inside, returns a LDA
model based on the input data and labels. The argument `distances` indicates
whether execution of the model will result in distances (from test observations
to the class-wise mean vectors of the transformed reference space) or the transformed
space of the test objects (distances can be subsequently calculated using this data
and the model).

# Keyword arguments (same as in `MultivariateStats`)
  * `method` can be `:gevd` e.g. generalized eigenvalue decomp or `:whiten` e.g. whiten first (default `:gevd`)
  * `outdim` output dimension e.g. of the transformed space (default `min(<number of variables>, <number_of_classes>-1)`
  * `regcoef` regularization coefficient, improves numerical stability (default `1e-6`)

Read the `MultivariateStats.jl` documentation for more information.
"""
lda(distances::Bool=true; kwargs...) = FunctionCell(lda, (distances,), ModelProperties(), kwtitle(distances ? "LDA (distance)" : "LDA (transform)", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	lda(x, distances=true [;kwargs])

Trains a LDA model that when executed will compute either the distances from each input
observation to the mean vectors of `x` in the transformed space (distances = true)
or alternatively, the LDA transform of the input observations.

"""
# Training
lda(x::T where T<:CellDataL, distances::Bool=true; kwargs...) = lda((getx!(x), gety(x)), distances; kwargs...)
lda(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, distances::Bool=true; kwargs...) = lda((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), distances; kwargs...)
lda(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, distances::Bool=true; kwargs...) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[lda] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Transform labels first
	enc = labelencn(x[2])
	yenc = label2ind.(x[2],enc)

	# Train model
	ldadata = fit(MultivariateStats.MulticlassLDA, length(enc.label), getobs(x[1]), yenc; kwargs...)

	# Build model properties
	modelprops = ModelProperties(ldadata.stats.dim, distances ? length(enc.label) : size(ldadata.proj,2), enc)

	FunctionCell(lda, Model((distances,ldadata), modelprops), kwtitle(distances ? "LDA (distance)" : "LDA (transform)", kwargs))
end

# Execution
lda(x::T where T<:CellData, model::Model{<:Tuple{<:Bool,<:MultivariateStats.MulticlassLDA}}) =
	datacell(lda(getx!(x), model), gety(x))
lda(x::T where T<:AbstractVector, model::Model{<:Tuple{<:Bool,<:MultivariateStats.MulticlassLDA}}) =
	lda(mat(x, LearnBase.ObsDim.Constant{2}()), model)
lda(x::T where T<:AbstractMatrix, model::Model{<:Tuple{<:Bool,<:MultivariateStats.MulticlassLDA}}) = begin

	if model.data[1] == true
		# Return distances between the mean class vectors and each sample in x
		Distances.pairwise(Distances.Euclidean(), model.data[2].pmeans, MultivariateStats.transform(model.data[2], getobs(x)))
	else
		# Return the transformed observations
		MultivariateStats.transform(model.data[2], getobs(x))
	end
end

##########################
# FunctionCell Interface #
##########################
"""
	ldasub(distances=true [;kwargs])

Constructs an untrained cell that when piped data inside, returns a LDA subspace
model based on the input data and labels. The argument `distances` indicates
whether execution of the model will result in distances (from test observations
to the class-wise mean vectors of the transformed reference space) or the transformed
space of the test objects (distances can be subsequently calculated using this data
and the model).

# Keyword arguments (same as in `MultivariateStats`)
  * `normalize` regularization coefficient, improves numerical stability (default `false`)

Read the `MultivariateStats.jl` documentation for more information.
"""
ldasub(distances::Bool=true; kwargs...) = FunctionCell(ldasub, (distances,), ModelProperties(),
					 kwtitle(distances ? "LDA-Subspace (distance)" : "LDA-Subspace (transform)", kwargs);
					 kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	ldasub(x, distances=true [;kwargs])

Trains a LDA subspace model that when executed will compute either the distances from each input
observation to the mean vectors of `x` in the transformed space (distances = true)
or alternatively, the LDA transform of the input observations.

"""
# Training
ldasub(x::T where T<:CellDataL, distances::Bool=true; kwargs...) = ldasub((getx!(x), gety(x)), distances; kwargs...)
ldasub(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, distances::Bool=true; kwargs...) = ldasub((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), distances; kwargs...)
ldasub(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, distances::Bool=true; kwargs...) = begin

	@assert nobs(x[1]) == nobs(x[2])

	# Transform labels first
	enc = labelencn(x[2])
	yenc = label2ind.(x[2],enc)

	# Train model
	ldadata = fit(MultivariateStats.SubspaceLDA, getobs(x[1]), yenc, length(enc.label); kwargs...)

	# Build model properties
	modelprops = ModelProperties(size(ldadata.cmeans,1), distances ? length(enc.label) : size(ldadata.projLDA,2), enc)

	FunctionCell(ldasub, Model((distances,ldadata), modelprops), kwtitle(distances ? "LDA-Subspace (distance)" : "LDA-subspace (transform)", kwargs))
end

# Execution
ldasub(x::T where T<:CellData, model::Model{<:Tuple{<:Bool,<:MultivariateStats.SubspaceLDA}}) =
	datacell(ldasub(getx!(x), model), gety(x))
ldasub(x::T where T<:AbstractVector, model::Model{<:Tuple{<:Bool,<:MultivariateStats.SubspaceLDA}}) =
	ldasub(mat(x, LearnBase.ObsDim.Constant{2}()), model)
ldasub(x::T where T<:AbstractMatrix, model::Model{<:Tuple{<:Bool,<:MultivariateStats.SubspaceLDA}}) = begin

	if model.data[1] == true
		# Return distances between the mean class vectors and each sample in x
		Distances.pairwise(Distances.Euclidean(),
		     MultivariateStats.transform(model.data[2], model.data[2].cmeans),
		     MultivariateStats.transform(model.data[2], getobs(x))
		)
	else
		# Return the transformed observations
		MultivariateStats.transform(model.data[2], getobs(x))
	end
end


text:
module ITensors

function pause()
  println("(paused)")
  readline(stdin)
end

using Random,
      Printf,
      LinearAlgebra,
      StaticArrays # For SmallString

import Base.adjoint,
       Base.conj,
       Base.convert,
       Base.copy,
       Base.copyto!,
       Base.eltype,
       Base.fill!,
       Base.getindex,
       Base.in,
       Base.isapprox,
       Base.isless,
       Base.iterate,
       Base.length,
       Base.push!,
       Base.setindex!,
       Base.eachindex,
       Base.show,
       Base.similar,
       Base.size,
       Base.!=,
       Base.+,
       Base.-,
       Base.*,
       Base./,
       Base.^,
       Base.complex,
       Base.setdiff,  # Since setdiff doesn't
                      # work with IndexSet, overload it
       LinearAlgebra.axpby!,
       LinearAlgebra.axpy!,
       LinearAlgebra.dot,
       LinearAlgebra.norm,
       LinearAlgebra.mul!,
       LinearAlgebra.rmul!,
       LinearAlgebra.normalize!,
       Random.randn!

#TODO: continue work on SmallString, use as Tags
include("smallstring.jl")
include("tagset.jl")
include("index.jl")
include("indexset.jl")
include("storage/tensorstorage.jl")
include("storage/dense.jl")
include("storage/contract.jl")
include("storage/svd.jl")
#export CProps, contract!, compute!, compute_contraction_labels, contract_inds, contract
include("itensor.jl")
include("decomp.jl")
include("iterativesolvers.jl")

###########################################################
# MPS/MPO
#
include("mps/siteset.jl")
include("mps/initstate.jl")
include("mps/mps.jl")
include("mps/mpo.jl")
include("mps/sweeps.jl")
include("mps/projmpo.jl")
include("mps/dmrg.jl")
include("mps/autompo.jl")

###########################################################
# Physics
#
include("physics/lattices.jl")
include("physics/sitesets/spinhalf.jl")
include("physics/sitesets/spinone.jl")
include("physics/sitesets/electron.jl")
include("physics/sitesets/tj.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 01db1fc8-5ffd-11eb-0dae-cda0ff89e63d
begin
	using MLJ, MLJModels
	using Random
	using BenchmarkTools
	using DataFrames
	using LeastSquaresSVM
end

# ╔═╡ b84726f0-6003-11eb-2e0d-6f14682c9588
md"""
# Baseline benchmark for multiclass classification

In this document we will benchmark the multiclass classification problem of the
Least Squares Support Vector Classifier from `LeastSquaresSVM` and the one from
the fantastic and very well-known `libSVM` solver.

The idea is to test both their prediction and training complexities, both
in time and space. To accomplish this, we will be using the `BenchmarkTools`
package.

We will test **two** cases here:

1. Training _with_ hyperparameter tuning;
2. training _without_ tuning.

This is because a well-tuned SVM is _very_ important to predict good results.
We almost always need a well-tuned SVM, so this is an important step.

The dataset we'll be using is the [Iris dataset](https://archive.ics.uci.edu/ml/datasets/Iris/).
This is a very simple and easy dataset, and we will deal with larger datasets
in another document elsewhere.
"""

# ╔═╡ 5d73b9a0-6005-11eb-2ee6-2164f4169fea
md"""
First, some imports and all that good stuff.
"""

# ╔═╡ 143d4bd4-5ffe-11eb-23dc-1751a1913a05
@load SVC pkg=LIBSVM

# ╔═╡ 6a111cfe-6005-11eb-15df-875c96454e44
md"""
## Loading the dataset

We now load the dataset, and we'll convert it to a `DataFrame` to inspect it.
We just want to check that the `scitypes` are correct.
"""

# ╔═╡ 281afffa-5ffd-11eb-24ce-219f746515da
X, y = @load_iris;

# ╔═╡ 6e9dcf16-5ffd-11eb-0523-e1cfd36fd57a
dfIris = DataFrame(X);

# ╔═╡ 82244ce8-5ffd-11eb-003d-23214d962d20
dfIris.y = y;

# ╔═╡ 878df3be-5ffd-11eb-2b19-7f74692d0189
first(dfIris, 3)

# ╔═╡ 8c358782-6005-11eb-2787-8b32b5749b57
md"""
## Training setup

Great, so everything looks good with the dataset.
Let us now setup the dataset for training and testing. We'll do a 60-40% split
on this dataset.
Recall that this dataset only has 150 datasamples. We also always shuffle the
splitting.
"""

# ╔═╡ 2bbc10ce-6000-11eb-2162-398983b5cf48
train, test = partition(eachindex(y), 0.6, shuffle=true);

# ╔═╡ cb51aff4-6005-11eb-1a7e-6b229f087a31
md"""
## Training functions

To correctly benchmark the training steps, we need to setup and create some
functions that collect all the training logic.

Here, we are doing two pairs of functions, one _without_ hyperparameter tuning
and the other _with_ tuning.
You can tell them apart because the ones _with_ hyperparameter tuning have
the suffix `hp`.

Also, to differentiate between both implementations, `svm` corresponds to the
`libsvm` implementation. On the other hand, `lssvm` corresponds to the
`LeastSquaresSVM` implementation.
"""

# ╔═╡ f6f1ffe8-6004-11eb-0d26-3d51bcb24d61
function train_lssvm(X, y, train)
	model = LSSVClassifier()

	mach = machine(model, X, y)
	fit!(mach, rows=train, verbosity=0)

	return mach
end;

# ╔═╡ feaf46c8-6004-11eb-1fe3-7792a33bbc0e
function train_svm(X, y, train)
	model = SVC()

	mach = machine(model, X, y)
	fit!(mach, rows=train, verbosity=0)

	return mach
end;

# ╔═╡ 92b48212-5ffd-11eb-3032-d1c9519703f2
function train_lssvm_hp(X, y, train)
	model = LSSVClassifier()

	r1 = range(model, :σ, lower=1e-2, upper=10.0)
	r2 = range(model, :γ, lower=1, upper=200.0)

	self_tuning_model = TunedModel(
		model=model,
		tuning=Grid(goal=625),
		resampling=CV(nfolds=5),
		range=[r1, r2],
		measure=accuracy
	)

	mach = machine(self_tuning_model, X, y)
	fit!(mach, rows=train)

	return mach
end;

# ╔═╡ 5973126c-5ffe-11eb-2219-85731b9c0786
function train_svm_hp(X, y, train)
	model = SVC()

	r1 = range(model, :cost, lower=1, upper=1000.0)
	r2 = range(model, :gamma, lower=1e-2, upper=10.0)

	self_tuning_model = TunedModel(
		model=model,
		tuning=Grid(goal=625),
		resampling=CV(nfolds=5),
		range=[r1, r2],
		measure=accuracy
	)

	mach = machine(self_tuning_model, X, y)
	fit!(mach, rows=train)

	return mach
end;

# ╔═╡ 2ee014fc-6006-11eb-367a-dbc9f9f97225
md"""
## Dataset prep

With that out of the way, we need to remove the mean from the dataset, as
well as setting its standard deviation to one.
This is referred to as **standardizing** the dataset, and `MLJ` already has
a model for such tasks.

We'll standardize our model before feeding it to our models, because SVMs
are quite picky of non-standardized data.
"""

# ╔═╡ 5e76126c-5ffe-11eb-2b0b-6fe6136909c5
Xstd = MLJ.transform(MLJ.fit!(MLJ.machine(Standardizer(), X)), X);

# ╔═╡ 6f0d1bec-6006-11eb-106b-bd780f26fef9
md"""
## Benchmarking

Let's begin with benchmarking.

### Without tuning

First, let's do a baseline benchmark, _without_ hyperparameter tuning.
Here, we are looking for time and space estimations for both implementations.

We'll start with the `LeastSquaresSVM` implementation.
"""

# ╔═╡ 0ab3cf70-6005-11eb-3f80-03a804a70776
@benchmark train_lssvm($Xstd, $y, $train)

# ╔═╡ ad027e4a-6006-11eb-3ae1-f11f7d5fecdd
md"""
#### Memory allocation

Right out of the bat, we can see that there is quite a lot of allocation.
In the order of KiB it might not seem much, but this will surely be
more promiment in the hyperparameter tuning estimations.

This is expected due to the following reasons:

- This formulation will **always** need **all** the dataset for both training and testing.
- The solver (conjugate gradient, the iterative formulation) has to deal with allocations due to having to create a basis for each iteration.

#### Time estimation

Following the [manual for `BenchmarkTools`](https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/doc/manual.md#which-estimator-should-i-use) we will use the **median** as out metric for
comparing both models. This is because it is not really affected by
noise or some other type of outliers.
"""

# ╔═╡ 85b1dd2c-6008-11eb-0689-9572c7411b30
md"""
With this in mind our median time is of ≈ 400 microseconds, which is pretty okay.
Some of this might come from some compiler optimizations and clever machine
code.

We'll move on to the next implementation to have something to compare it with.
"""

# ╔═╡ 0a75a496-6005-11eb-2135-eb8b6d4d106f
@benchmark train_svm($Xstd, $y, $train)

# ╔═╡ 8ee94646-6008-11eb-089a-7b2ea0a9497b
md"""
#### Memory allocation

Wow, this is an order of maginitude lower compared to the `LeastSquaresSVM`
implementation. But then, again, this is expected.
The `libSVM` solver is **sparse**, meaning that it only selects a couple of data
instances as its _support vectors_.

#### Time estimation

By comparing with the previous _median_ value, this implementation is faster,
and better.
I kind of expected this result due to the fact the the `libSVM` solver is very
_highly_ optimized code, and written almost entirely in C/C++.

And by highly optimized code I mean that they use really clever tricks in various
numerical computations throughout their implementation.

But even so the `LeastSquaresSVM` implementation is very close. Not too shabby.
"""

# ╔═╡ f6c9278a-6009-11eb-300b-1fb5f3a9ad13
md"""
### With tuning. Setup and general info.

The hyperparameter search is done with the [`Grid`](https://alan-turing-institute.github.io/MLJ.jl/dev/tuning_models/#MLJTuning.Grid)
algorithm from `MLJ`, which is a conventional grid seach.

The reason for this choice is simple: we want to be fair.

Each model tunes its hyperparameters using 5-fold cross-validation on the
_training_ set only.

Each model has to tune **two** hyperparameters:

- an _intrinsic_ one; _cost_ for the `libSVM` and _γ_ for the `LeastSquaresSVM`,
- the kernel hyperparameter; both models are using the same RBF kernel.

Each model will search for the best pair of hyperparameters in a
square grid of 625 points. That is, 25 values for each hyperparameter.

It might seem like very few elements, but keep in mind that no parallelization
is used for this search (which normally it is recommended to do).
To make the comparison fair, we are doing this in a single thread, for each
model.

The performance metric will be the classification _accuracy._
"""

# ╔═╡ 5e3c2068-5ffe-11eb-3173-59eb4527939f
t1 = @benchmark train_lssvm_hp($Xstd, $y, $train)

# ╔═╡ da5d6972-602b-11eb-3cce-a50ed90928e0
md"""
#### Memory allocations

Well, this is not a surprise. For each pair of points, a model is trained.
The whole dataset set (the training portion of it) is used each time. For now
it is only comprised of $(length(train)) data instances. But imagine how large
can this grow for datasets with more than just 150 data samples.

This is one of the reasons why the least squares formulation is not the popular
one. Even though it is very simple to implement, it consumes memory like there
is no tomorrow.

#### Time estimation

This is quite a surprise. I seriously did not expect it to perform as well as
it does.
It is fast, efficient, but I think this is mostly due to the nature of the iterative
solver it uses.
If a more direct approach were used, I firmly believe the running time would have
blown up.
"""

# ╔═╡ 5e286ffa-5ffe-11eb-22c3-2bba58505b97
t2 = @benchmark train_svm_hp($Xstd, $y, $train)

# ╔═╡ abe9a28a-602c-11eb-18a9-493b4b05d461
md"""
#### Memory allocations

Again, this is not a surprise. By using just a few data instances, the `libSVM`
implementation is _very_ conservative, memory-wise.

Do note, however, that the difference in memory usage is mind-blowing.
The difference between this implementation and the least squares one is
$(round(t1.memory / t2.memory; digits=4)) smaller.

#### Time estimation

This is really the weird part. I would have expected that this implementation
would be faster.
Maybe it has something to do with the number of samples from the benchmark.
Anyway, the results show that both implementations are equivalent.

This is great news, because it means that no matter which implementation you use,
you will always have your results in a timely manner.
You can choose between both implementations and the prediction time is almost the
same.
"""

# ╔═╡ 00fdae46-601a-11eb-1db6-e73e0f4ead1e
md"""
## Prediction

In this last section I just wanted to show the prediction time for both
implementations.

We are assuming that both models have been tuned correctly.

First, I will wrap the prediction steps in a function, evaluation the
_accuracy_ of prediction for each model.
"""

# ╔═╡ 5e0dc6e6-5ffe-11eb-0b25-8fa52bd3cfc1
function predict_model(y, test, mach)
	results = MLJ.predict(mach, rows=test)
	acc = MLJ.accuracy(results, y[test])

	return acc
end;

# ╔═╡ 61c5201a-602e-11eb-15a7-1f0334dfa695
md"""
I need a tuned least squares SVM, so I do that first.
"""

# ╔═╡ 66054bac-6002-11eb-063c-5525eb2a0ca1
mach1 = train_lssvm_hp(Xstd, y, train);

# ╔═╡ c55d8cb4-6001-11eb-2e6b-655ac8539748
@benchmark predict_model($y, $test, $mach1)

# ╔═╡ 949df2e8-6002-11eb-14cb-e5a269cf2846
acc1 = predict_model(y, test, mach1)

# ╔═╡ 77449dee-602e-11eb-3d23-6baa12c216c1
md"""
### Results for least squares implementation

Good accuracy, as expected. Very high memory consumption, as expected.
The running time is quite high, but nothing to really worry about.
"""

# ╔═╡ 6e868db0-6002-11eb-0cc1-edf2da791ccb
mach2 = train_svm_hp(Xstd, y, train);

# ╔═╡ c5477f70-6001-11eb-3a08-87356dfbbd29
@benchmark predict_model($y, $test, $mach2)

# ╔═╡ c53166f2-6001-11eb-19d9-5d104a8e7107
acc2 = predict_model(y, test, mach2)

# ╔═╡ a75bd45c-602e-11eb-3206-8b097f43be6c
md"""
### Results for `libSVM` implementation

Good accuracy, as expected. Low memory consumption, as expected.

Now, regarding the running time, this is definetely what I was expecting
on the _training_ step. We can clearly see, by comparing the median time between
both benchmarks, that the `libSVM` implementation is almost **twice** as fast
as the least squares implementation.

If it is using less data, that means less computations. But also, recall that this
solver (and code) is very optimized.

## Conclusions

Most of the results presented here represent a baseline benchmark result on a simple
multiclass classification problem.

The results are expected in most cases. Although it is good to see the
`LeastSquaresSVM` doing well on the running time, both in _training_ and in
_prediction._
"""

# ╔═╡ Cell order:
# ╟─b84726f0-6003-11eb-2e0d-6f14682c9588
# ╟─5d73b9a0-6005-11eb-2ee6-2164f4169fea
# ╠═01db1fc8-5ffd-11eb-0dae-cda0ff89e63d
# ╠═143d4bd4-5ffe-11eb-23dc-1751a1913a05
# ╟─6a111cfe-6005-11eb-15df-875c96454e44
# ╠═281afffa-5ffd-11eb-24ce-219f746515da
# ╠═6e9dcf16-5ffd-11eb-0523-e1cfd36fd57a
# ╠═82244ce8-5ffd-11eb-003d-23214d962d20
# ╠═878df3be-5ffd-11eb-2b19-7f74692d0189
# ╟─8c358782-6005-11eb-2787-8b32b5749b57
# ╠═2bbc10ce-6000-11eb-2162-398983b5cf48
# ╟─cb51aff4-6005-11eb-1a7e-6b229f087a31
# ╠═f6f1ffe8-6004-11eb-0d26-3d51bcb24d61
# ╠═feaf46c8-6004-11eb-1fe3-7792a33bbc0e
# ╠═92b48212-5ffd-11eb-3032-d1c9519703f2
# ╠═5973126c-5ffe-11eb-2219-85731b9c0786
# ╟─2ee014fc-6006-11eb-367a-dbc9f9f97225
# ╠═5e76126c-5ffe-11eb-2b0b-6fe6136909c5
# ╟─6f0d1bec-6006-11eb-106b-bd780f26fef9
# ╠═0ab3cf70-6005-11eb-3f80-03a804a70776
# ╟─ad027e4a-6006-11eb-3ae1-f11f7d5fecdd
# ╟─85b1dd2c-6008-11eb-0689-9572c7411b30
# ╠═0a75a496-6005-11eb-2135-eb8b6d4d106f
# ╟─8ee94646-6008-11eb-089a-7b2ea0a9497b
# ╟─f6c9278a-6009-11eb-300b-1fb5f3a9ad13
# ╠═5e3c2068-5ffe-11eb-3173-59eb4527939f
# ╟─da5d6972-602b-11eb-3cce-a50ed90928e0
# ╠═5e286ffa-5ffe-11eb-22c3-2bba58505b97
# ╟─abe9a28a-602c-11eb-18a9-493b4b05d461
# ╟─00fdae46-601a-11eb-1db6-e73e0f4ead1e
# ╠═5e0dc6e6-5ffe-11eb-0b25-8fa52bd3cfc1
# ╟─61c5201a-602e-11eb-15a7-1f0334dfa695
# ╠═66054bac-6002-11eb-063c-5525eb2a0ca1
# ╠═c55d8cb4-6001-11eb-2e6b-655ac8539748
# ╠═949df2e8-6002-11eb-14cb-e5a269cf2846
# ╟─77449dee-602e-11eb-3d23-6baa12c216c1
# ╠═6e868db0-6002-11eb-0cc1-edf2da791ccb
# ╠═c5477f70-6001-11eb-3a08-87356dfbbd29
# ╠═c53166f2-6001-11eb-19d9-5d104a8e7107
# ╟─a75bd45c-602e-11eb-3206-8b097f43be6c

max_stars_repo_name: edwinb-ai/LSSVM-benchmarks

max_issues_repo_name: edwinb-ai/LSSVM-benchmarks

max_forks_repo_name: edwinb-ai/LSSVM-benchmarks