text:
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using Psychometrics
using JLD2
using Statistics
using DelimitedFiles
using FileIO
using StatsBase
using Distributions

function online_calibrate()
replications = 1:1
I_operational = 1_000
I_field = 25
I_total = I_field + I_operational
N = 3_000

test_operational  = 25
test_field = 5
test_length = test_field + test_operational

iter_mcmc_latent = 2_000
iter_mcmc_item = 4_000
#after how many responses update item parameter estimates
batch_size = 100
N_T = 500

items = load(string("examples/online-calibration/true values/true_items.jld2"), "items")
examinees = load(string("examples/online-calibration/true values/true_examinees.jld2"), "examinees")[1:N]

for rep in replications

    I_operational = 1000
    I_field = 25
    I_total = I_field + I_operational
    N = 3_000

    test_operational  = 25
    test_field = 5
    test_length = test_field + test_operational

    iter_mcmc_latent = 2_000
    iter_mcmc_item = 2_000
    #after how many responses update item parameter estimates
    batch_size = 50
    N_T = 500
    items = load(string("examples/online-calibration/true values/true_items.jld2"), "items")
    examinees = load(string("examples/online-calibration/true values/true_examinees.jld2"), "examinees")

    retired_items = 0
    retired_items_vector= fill(0, N)

    # RESPONSES

    responses_per_examinee = [Response[] for e in examinees];

    # INITIAL VALUES

    ## ITEMS
    # first I_total-I_field items have true values of parameters (very well estimated)
    items_est_operational = items[1:I_operational];
    a_operational = map(i -> i.parameters.a, items_est_operational)
    b_operational = map(i -> i.parameters.b, items_est_operational)
    a_est_prior = TruncatedNormal(mean(a_operational), std(a_operational), 0,  Inf);
    println(params(a_est_prior))
    a_est_bounds = [1e-5, 5.0];
    b_est_prior = Normal(mean(b_operational), std(b_operational));
    b_est_bounds = [-6.0, 6.0];

    items_est_field = [Item(i+I_operational, string("item_", i+I_operational), ["math"], Parameters2PL(Product([a_est_prior, b_est_prior]), a_est_bounds, b_est_bounds)) for i = 1 : I_field];
    map( i -> begin
        i.parameters.calibrated = false
    end, items_est_field);
    # starting values
    map( i -> begin
        i.parameters.a = Float64(mean(a_operational))
    end, items_est_field);
    map( i -> begin
        i.parameters.b = mean(b_operational)
    end, items_est_field);
    # starting expected information matrices
    map( i -> begin

    i.parameters.expected_information = Psychometrics._item_expected_information(i.parameters, Latent1D(0.0))
    end, items_est_field)
    # append not calibrated items
    items_est = vcat(items_est_operational, items_est_field);
    #responses vectors for not calibrated items
    responses_not_calibrated = Response[]

    ## EXAMINEES
    latent_est_prior = Normal(0, 3);
    latent_est_bounds = [-6.0, 6.0];
    examinees_est = [Examinee(e, string("examinee_",e), Latent1D(latent_est_prior, latent_est_bounds)) for e = 1 : N];

    #set starting value taking a random value in the interval -0.5:0.5
    map(e ->
            begin
                e.latent.val = Distributions.rand(-0.3:0.3)
            end,
    examinees_est)
    items_idx_per_examinee = Vector{Vector{Int64}}(undef, N);
    examinees_est_theta = [[zero(Float64) for i=1:test_length] for n=1:N];
    items_est_parameters = [[[0.0, 0.0] for n=1:ceil(N_T/batch_size)] for i=1:I_field]
    #@save string("examples/online-calibration/true values/priors.jld2") a_est_prior b_est_prior latent_est_prior
    #@save string("examples/online-calibration/true values/bounds.jld2") a_est_bounds b_est_bounds latent_est_bounds

    # START ONLINE-CALIBRATION

    n=1
    while retired_items < I_field #size(filter(i -> i.parameters.calibrated == false, items_est),1) > 0
        examinee = copy(examinees_est[n])
        examinee_true = copy(examinees[n])
        idx_n = examinee.idx
        items_idx_n = Int64[]
        responses_n = Response[]
        println("n: ", n)
        #println("true: ", examinees[n].latent.val)
        available_items_idx = map(i -> i.idx, filter(i2 -> i2.parameters.calibrated, items_est))
        for i in 1:test_operational
            # find best item idx
            next_item_idx = find_best_item(examinee, items_est[available_items_idx]);
            push!(items_idx_n, next_item_idx);
            available_items_idx = setdiff(available_items_idx, next_item_idx)
            # answer to the item
            resp = answer(examinee_true, items[next_item_idx])
            push!(responses_n, resp)
            Psychometrics.assess_examinee_pg!(examinee, items_est[items_idx_n], responses_n; mcmc_iter = iter_mcmc_latent, sampling = false, already_sorted = true)
            #println("est_BIAS: ", examinee.latent.val - examinees[n].latent.val)
            # store theta estimate
            examinees_est_theta[n][i] = examinee.latent.val
        end
        # find available items to calibrate
        available_items_idx = map(i -> i.idx, filter(i2 -> !i2.parameters.calibrated, items_est))
        for i in 1:test_field
            if size(available_items_idx,1) > 0
                # find best item to be calibrated
                next_item_idx = find_best_item(examinee, items_est[available_items_idx]; method = "D-gain");
                push!(items_idx_n, next_item_idx);
                available_items_idx = setdiff(available_items_idx, next_item_idx)
                # answer to the item
                resp = answer(examinee_true, items[next_item_idx])
                push!(responses_n, resp)

                # set the theta prior as theta posterior
                #examinee.latent.prior = examinee.latent.posterior
                examinees_est[n] = copy(examinee)
                # store theta estimate
                #println("examinee.latent.val: ", examinee.latent.val)
                examinees_est_theta[n][i+(test_operational)] = examinee.latent.val
                #println("est_BIAS: ", examinees_est[n].latent.val - examinees[n].latent.val)
                push!(responses_not_calibrated, resp)
                item = copy(items_est[next_item_idx])
                #item.parameters.expected_information += item_expected_information(item.parameters, examinee.latent)
                resp_item = get_responses_by_item_idx(item.idx, responses_not_calibrated)
                if mod(size(resp_item, 1), batch_size) == 0
                    examinees_est_item = examinees_est[map( r -> r.examinee_idx, resp_item)]
                    Psychometrics.calibrate_item_pg!(item, examinees_est_item, resp_item; mcmc_iter = iter_mcmc_item, sampling = false, already_sorted = true)
                    #save chain
                    if rep ==1
                        chain = item.parameters.chain
                        @save string("examples/online-calibration/rep_", rep, "_item_", item.idx, "_resp_", size(resp_item, 1), "_chain.jld2") chain
                    end
                    items_est_parameters[next_item_idx - I_operational][Int64(size(resp_item,1)/batch_size)] = [item.parameters.a, item.parameters.b]
                    if next_item_idx==1001
                    println("est_BIAS a : ", sqrt((item.parameters.a - items[next_item_idx].parameters.a)^2))
                    println("est_BIAS b : ", sqrt((item.parameters.b - items[next_item_idx].parameters.b)^2))
                    end
                    if size(resp_item, 1) >= N_T
                        item.parameters.calibrated = true
                        retired_items +=1
                    end
                    #update expected information
                    item.parameters.expected_information = mapreduce(e -> Psychometrics._item_expected_information(item.parameters, e.latent), +, examinees_est_item)
                else
                    #update expected information
                    item.parameters.expected_information += Psychometrics._item_expected_information(item.parameters, examinee.latent)
                end
                items_est[next_item_idx] = item
            end
        end
        examinees_est[n] = copy(examinee)
        #println("est_BIAS: ", examinee.latent.val - examinees[n].latent.val)
        responses_per_examinee[n] = copy(responses_n)
        items_idx_per_examinee[n] = copy(items_idx_n)
        retired_items_vector[n] = copy(retired_items)
        n += 1
    end
    n -= 1
    @save string("examples/online-calibration/how_many_examinees_rep_",rep,".jld2") n
    @save string("examples/online-calibration/retired_items_vector_rep_",rep,".jld2") retired_items_vector
    @save string("examples/online-calibration/examinees_est_theta_rep_",rep,".jld2") examinees_est_theta
    @save string("examples/online-calibration/items_est_parameters_rep_",rep,".jld2") items_est_parameters

    @save string("examples/online-calibration/items_est_rep_",rep,".jld2") items_est
    @save string("examples/online-calibration/examinees_est_rep_",rep,".jld2") examinees_est

    rmse_a =  sqrt(mean(map( (i, i_est)-> (i.parameters.a - i_est.parameters.a)^2, items[I_operational:I_total], items_est[I_operational:I_total])))
    rmse_b = sqrt(mean(map( (i, i_est)-> (i.parameters.b - i_est.parameters.b)^2, items[I_operational:I_total], items_est[I_operational:I_total])))
    rmse_theta = sqrt(mean(map( (e, e_est)-> (e.latent.val- e_est.latent.val)^2, examinees[1:(n-1)], examinees_est[1:(n-1)])))
    println( "avg RMSE a = ", rmse_a)
    println( "avg RMSE b = ", rmse_b)
    println( "avg RMSE tehta = ", rmse_theta)
    items_est_field = items_est[(I_operational+1):I_total]
    @save string("examples/online-calibration/rep_", rep, "_items_est_field.jld2") items_est_field
    @save string("examples/online-calibration/rep_", rep, "_responses.jld2") responses_per_examinee
    @save string("examples/online-calibration/rep_", rep, "_items_est.jld2") items_est
    @save string("examples/online-calibration/rep_", rep, "_responses.jld2") responses_per_examinee
end
# theta_sets = collect(-6.0:0.5:6.0)
# rmse_theta = zeros(Float64,(size(theta_sets,1) ))
# gap = theta_sets[2] - theta_sets[1]
# for t in 1:size(theta_sets,1)
#     println(t)
#         e_t_idx = map( e -> e.idx, filter( e2-> (e2.latent.val >= theta_sets[t] - (gap/2)) && (e2.latent.val < theta_sets[t] +(gap/2)), examinees[1:n-1]))
#         if size(e_t_idx,1)>1
#             rmse_theta[t] = sqrt(mean(map((e,e_est) -> (e.latent.val-e_est.latent.val)^2, examinees[e_t_idx], examinees_est[e_t_idx])))
#         else
#             rmse_theta[t] = 0.0
#         end
# end
# plot(theta_sets, rmse_theta)

# rep=100

# a_s = zeros(I_total,rep)
# b_s = zeros(I_total,rep)
# for rep = 1:rep
#     @load string("examples/online-calibration/rep_", rep, "_items_est.jld2") items_est
#     a = map(i -> i.parameters.a, items_est)
#     a_s[:, rep] = copy(a)
#     b = map(i -> i.parameters.b, items_est)
#     b_s[:, rep] = copy(b)
# end
a_s = zeros(1020,rep)
b_s = zeros(1020,rep)
for rep = 1:rep
    @load string("examples/online-calibration/rep_", rep, "_items_est.jld2") items_est
    a = map(i -> i.parameters.a, items_est)
    a_s[:, rep] = copy(a)
    b = map(i -> i.parameters.b, items_est)
    b_s[:, rep] = copy(b)
end

# a_true = map(i -> i.parameters.a, items)
# a_true = hcat([a_true for rep=1:rep]...)
# b_true = map(i -> i.parameters.b, items)
# b_true = hcat([b_true for rep=1:rep]...)
# diff_a = (a_s[1001:1020,:] - a_true[1001:1020,:]).^2
# diff_b = (b_s[1001:1020,:] - b_true[1001:1020,:]).^2

# println(sqrt(sum(diff_a)/I_field/rep))
# println(sqrt(sum(diff_b)/I_field/rep))
return examinees, examinees_est, items, items_est, responses
end
examinees, examinees_est, items, items_est, responses = online_calibrate();


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function solution(model::JuMP.Model)::OrderedDict
    instance, T = model[:instance], model[:instance].time
    function timeseries(vars, collection)
        return OrderedDict(
            b.name => [round(value(vars[b.name, t]), digits = 5) for t in 1:T]
            for b in collection
        )
    end
    function production_cost(g)
        return [
            value(model[:is_on][g.name, t]) * g.min_power_cost[t] + sum(
                Float64[
                    value(model[:segprod][g.name, t, k]) *
                    g.cost_segments[k].cost[t] for
                    k in 1:length(g.cost_segments)
                ],
            ) for t in 1:T
        ]
    end
    function production(g)
        return [
            value(model[:is_on][g.name, t]) * g.min_power[t] + sum(
                Float64[
                    value(model[:segprod][g.name, t, k]) for
                    k in 1:length(g.cost_segments)
                ],
            ) for t in 1:T
        ]
    end
    function startup_cost(g)
        S = length(g.startup_categories)
        return [
            sum(
                g.startup_categories[s].cost *
                value(model[:startup][g.name, t, s]) for s in 1:S
            ) for t in 1:T
        ]
    end
    sol = OrderedDict()
    sol["Production (MW)"] =
        OrderedDict(g.name => production(g) for g in instance.units)
    sol["Production cost (\$)"] =
        OrderedDict(g.name => production_cost(g) for g in instance.units)
    sol["Startup cost (\$)"] =
        OrderedDict(g.name => startup_cost(g) for g in instance.units)
    sol["Is on"] = timeseries(model[:is_on], instance.units)
    sol["Switch on"] = timeseries(model[:switch_on], instance.units)
    sol["Switch off"] = timeseries(model[:switch_off], instance.units)
    sol["Reserve (MW)"] = timeseries(model[:reserve], instance.units)
    sol["Reserve shortfall (MW)"] = OrderedDict(
        t =>
            (instance.shortfall_penalty[t] >= 0) ?
            round(value(model[:reserve_shortfall][t]), digits = 5) : 0.0 for
        t in 1:instance.time
    )
    sol["Net injection (MW)"] =
        timeseries(model[:net_injection], instance.buses)
    sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses)
    if !isempty(instance.lines)
        sol["Line overflow (MW)"] = timeseries(model[:overflow], instance.lines)
    end
    if !isempty(instance.price_sensitive_loads)
        sol["Price-sensitive loads (MW)"] =
            timeseries(model[:loads], instance.price_sensitive_loads)
    end
    return sol
end

ext: jl


text:
mutable struct Class_point
    x::Array{Float64,1}
    y::Array{Float64,1}
    s::Array{Float64,1}
    mu::Float64
    primal_scale::Float64

    function Class_point(x::Array{Float64,1},y::Array{Float64,1},s::Array{Float64,1},mu::Float64)
        @assert(length(s) == length(y))
        return new(x, y, s, mu, 1.0)
    end
end

#=function shrink_direction!(dir::Class_point, by::Float64)
    dir.x = dir.x * by
    dir.y = dir.y * by
    dir.s = dir.s * by
    dir.mu = dir.mu * by;
    dir.mu = dir.mu * by;
end=#

function zero_point(dim::Int64,ncon::Int64)
    return Class_point(zeros(dim),zeros(ncon),zeros(ncon),0.0)
end

#=
type Class_point
    x::Array{Float64,1}
    ncon::Int64

    i_l::Array{Int64,1} # variable to constraints
    y_l::Array{Float64,1}
    s_l::Array{Float64,1}

    i_u::Array{Int64,1}
    y_u::Array{Float64,1}
    s_u::Array{Float64,1}

    mu::Float64

    function Class_point(nvar::Int64, i_l::Array{Int64,1}, i_u::Array{Int64,1})
        this = new();
        this.x = zeros(nvar)
        this.i_l = i_l;
        this.i_u = i_u;

        this.y_l = zeros(length(i_l))
        this.s_l = zeros(length(i_l))

        this.s_u = zeros(length(i_u))
        this.y_u = zeros(length(i_u))

        this.mu = 0.0;

        return this
    end
end

function get_y_divide_s(p::Class_point)
    y_divide_s_l = point.y_l ./ point.s_l
    y_divide_s_u = point.y_u ./ point.s_u

    result = zeros(p.ncon)
    result[p.i_l] += y_divide_s_l
    result[p.i_u] += y_divide_s_u

    return result
end

function get_net_y(p::Class_point)
    y = zeros(p.ncon)
    y[p.i_l] += p.y_l
    y[p.i_u] += p.y_u

    return y
end

function zero_point(point::Class_point)
    new_point = deepcopy(point)
    point.x = zeros(length(point.x))

    nlb = length(point.i_l)
    point.s_l = zeros(nlb)
    point.y_l = zeros(nlb)

    nub = length(point.i_u)
    point.s_u = zeros(nub)
    point.y_u = zeros(nub)

    point.mu = 0.0

    return point
end=#


text:
# Linear regression (using iterative update)

using Regression
import Regression: solve, Options

# Data

d = 5
k = 3
n = 1000

W = randn(k, d+1)
X = randn(d, n)
Y = (W[:, 1:d] * X .+ W[:,d+1]) + 0.2 * randn(k, n)

# Solve

ret = solve(linearreg(X, Y; bias=1.0);
            solver=GD(),
            reg=SqrL2Reg(1.0e-3),
            options=Options(verbosity=:iter, grtol=1.0e-6 * n))

println()

# Print results

W_est = ret.sol
@printf("relative error = %.4e\n", sumabs2(W_est - W) / sumabs2(W))


text:
module StarJumpBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Star Jump Block" => Ahorn.EntityPlacement(
        Maple.StarJumpBlock,
        "rectangle"
    ),
    "Star Climb Controller" => Ahorn.EntityPlacement(
        Maple.StarClimbController,
    ),
)

Ahorn.minimumSize(entity::Maple.StarJumpBlock) = 8, 8
Ahorn.resizable(entity::Maple.StarJumpBlock) = true, true

Ahorn.selection(entity::Maple.StarJumpBlock) = Ahorn.getEntityRectangle(entity)

function Ahorn.selection(entity::Maple.StarClimbController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

function getStarjumpRectangles(room::Maple.Room)
    entities = filter(e -> e.name == "starJumpBlock", room.entities)
    rects = Ahorn.Rectangle[
        Ahorn.Rectangle(
            Int(get(e.data, "x", 0)),
            Int(get(e.data, "y", 0)),
            Int(get(e.data, "width", 8)),
            Int(get(e.data, "height", 8))
        ) for e in entities
    ]

    return rects
end

# Is there a star block we should connect to at the offset?
function noAdjacent(entity::Maple.StarJumpBlock, ox::Integer, oy::Integer, rects::Array{Ahorn.Rectangle, 1})
    x, y = Ahorn.position(entity)

    rect = Ahorn.Rectangle(x + ox + 4, y + oy + 4, 1, 1)

    return !any(Ahorn.checkCollision.(rects, Ref(rect)))
end

fillColor = (255, 255, 255, 255) ./ 255
corners = String[
    "objects/starjumpBlock/corner00.png",
    "objects/starjumpBlock/corner01.png",
    "objects/starjumpBlock/corner02.png",
    "objects/starjumpBlock/corner03.png"
]
edgeHs = String[
    "objects/starjumpBlock/edgeH00.png",
    "objects/starjumpBlock/edgeH01.png",
    "objects/starjumpBlock/edgeH02.png",
    "objects/starjumpBlock/edgeH03.png",
]
edgeVs = String[
    "objects/starjumpBlock/edgeV00.png",
    "objects/starjumpBlock/edgeV01.png",
    "objects/starjumpBlock/edgeV02.png",
    "objects/starjumpBlock/edgeV03.png",
]
leftRailings = String[
    "objects/starjumpBlock/leftrailing00.png"
    "objects/starjumpBlock/leftrailing01.png"
    "objects/starjumpBlock/leftrailing02.png"
    "objects/starjumpBlock/leftrailing03.png"
    "objects/starjumpBlock/leftrailing04.png"
    "objects/starjumpBlock/leftrailing05.png"
    "objects/starjumpBlock/leftrailing06.png"
]
rightRailings = String[
    "objects/starjumpBlock/rightrailing00.png"
    "objects/starjumpBlock/rightrailing01.png"
    "objects/starjumpBlock/rightrailing02.png"
    "objects/starjumpBlock/rightrailing03.png"
    "objects/starjumpBlock/rightrailing04.png"
    "objects/starjumpBlock/rightrailing05.png"
    "objects/starjumpBlock/rightrailing06.png"
]
railings = String[
    "objects/starjumpBlock/railing00.png"
    "objects/starjumpBlock/railing01.png"
    "objects/starjumpBlock/railing02.png"
    "objects/starjumpBlock/railing03.png"
    "objects/starjumpBlock/railing04.png"
    "objects/starjumpBlock/railing05.png"
    "objects/starjumpBlock/railing06.png"
]

function renderStarjumpBlock(ctx::Ahorn.Cairo.CairoContext, entity::Maple.StarJumpBlock, room::Maple.Room)
    starJumpRectangles = getStarjumpRectangles(room)
    rng = Ahorn.getSimpleEntityRng(entity)

    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    # Horizontal Border
    w = 8
    while w < width - 8
        if noAdjacent(entity, w, -8, starJumpRectangles)
            edge = rand(rng, edgeHs)
            Ahorn.drawSprite(ctx, edge, w + 4, 4)

            rail = railings[mod1(floor(Int, (x + w) / 8), length(railings))]
            Ahorn.drawSprite(ctx, rail, w, -8, jx=0.0, jy=0.0)
        end

        if noAdjacent(entity, w, height, starJumpRectangles)
            texture = rand(rng, edgeHs)
            Ahorn.drawSprite(ctx, texture, w + 4, height - 4, sy=-1)
        end

        w += 8
    end

    # Vertical Border
    h = 8
    while h < height - 8
        if noAdjacent(entity, -8, h, starJumpRectangles)
            texture = rand(rng, edgeVs)
            Ahorn.drawSprite(ctx, texture, 4, h + 4, sx=-1)
        end

        if noAdjacent(entity, width, h, starJumpRectangles)
            texture = rand(rng, edgeVs)
            Ahorn.drawSprite(ctx, texture, width - 4, h + 4)
        end

        h += 8
    end

    # Top Left Corner
    if noAdjacent(entity, -8, 0, starJumpRectangles) && noAdjacent(entity, 0, -8, starJumpRectangles)
        corner = rand(rng, corners)
        Ahorn.drawSprite(ctx, corner, 4, 4, sx=-1)

        rail = leftRailings[mod1(floor(Int, (x + w) / 8), length(leftRailings))]
        Ahorn.drawSprite(ctx, rail, 0, -8, jx=0.0, jy=0.0)

    elseif noAdjacent(entity, -8, 0, starJumpRectangles)
        corner = rand(rng, edgeVs)
        Ahorn.drawSprite(ctx, corner, 4, 4, sx=-1)

    elseif noAdjacent(entity, 0, -8, starJumpRectangles)
        corner = rand(rng, edgeHs)
        Ahorn.drawSprite(ctx, corner, 4, 4, sx=-1)

        rail = leftRailings[mod1(floor(Int, (x + w) / 8), length(leftRailings))]
        Ahorn.drawSprite(ctx, rail, 0, -8, jx=0.0, jy=0.0)
    end

    # Top Right Corner
    if noAdjacent(entity, width, 0, starJumpRectangles) && noAdjacent(entity, width - 8, -8, starJumpRectangles)
        corner = rand(rng, corners)
        Ahorn.drawSprite(ctx, corner, width - 4, 4)

        rail = rightRailings[mod1(floor(Int, (x + w) / 8), length(rightRailings))]
        Ahorn.drawSprite(ctx, rail, width - 8, -8, jx=0.0, jy=0.0)

    elseif noAdjacent(entity, width, 0, starJumpRectangles)
        corner = rand(rng, edgeVs)
        Ahorn.drawSprite(ctx, corner, width - 4, 4)

    elseif noAdjacent(entity, width - 8, -8, starJumpRectangles)
        corner = rand(rng, edgeHs)
        Ahorn.drawSprite(ctx, corner, width - 4, 4)

        rail = rightRailings[mod1(floor(Int, (x + w) / 8), length(rightRailings))]
        Ahorn.drawSprite(ctx, rail, width - 8, -8, jx=0.0, jy=0.0)
    end

    # Bottom Left Corner
    if noAdjacent(entity, -8, height - 8, starJumpRectangles) && noAdjacent(entity, 0, height, starJumpRectangles)
        corner = rand(rng, corners)
        Ahorn.drawSprite(ctx, corner, 4, height - 4, sx=-1, sy=-1)

    elseif noAdjacent(entity, -8, height - 8, starJumpRectangles)
        corner = rand(rng, edgeVs)
        Ahorn.drawSprite(ctx, corner, 4, height - 4, sx=-1, sy=-1)

    elseif noAdjacent(entity, 0, height, starJumpRectangles)
        corner = rand(rng, edgeHs)
        Ahorn.drawSprite(ctx, corner, 4, height - 4, sx=-1, sy=-1)
    end

    # Bottom Right Corner
    if noAdjacent(entity, width, height - 8, starJumpRectangles) && noAdjacent(entity, width - 8, height, starJumpRectangles)
        corner = rand(rng, corners)
        Ahorn.drawSprite(ctx, corner, width - 4, height - 4, sy=-1)

    elseif noAdjacent(entity, width, height - 8, starJumpRectangles)
        corner = rand(rng, edgeVs)
        Ahorn.drawSprite(ctx, corner, width - 4, height - 4, sy=-1)

    elseif noAdjacent(entity, width - 8, height, starJumpRectangles)
        corner = rand(rng, edgeHs)
        Ahorn.drawSprite(ctx, corner, width - 4, height - 4, sy=-1)
    end
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.StarJumpBlock, room::Maple.Room) = renderStarjumpBlock(ctx, entity, room)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.StarClimbController, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.northernLights, -12, -12)

end

max_stars_repo_name: zyowo/Ahorn

max_issues_repo_name: zyowo/Ahorn

max_forks_repo_name: zyowo/Ahorn


text:
# MetidaCu
# Copyright © 2019-2020 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>
using Metida
using MetidaCu
using Test, CSV, DataFrames, StatsModels, CategoricalArrays


text:
export BrickletJoystickPosition
struct BrickletJoystickPosition
    x::Integer
    y::Integer
end

export BrickletJoystickAnalogValue
struct BrickletJoystickAnalogValue
    x::Integer
    y::Integer
end

export BrickletJoystickPositionCallbackThreshold
struct BrickletJoystickPositionCallbackThreshold
    option::Char
    min_x::Integer
    max_x::Integer
    min_y::Integer
    max_y::Integer
end

export BrickletJoystickAnalogValueCallbackThreshold
struct BrickletJoystickAnalogValueCallbackThreshold
    option::Char
    min_x::Integer
    max_x::Integer
    min_y::Integer
    max_y::Integer
end

export BrickletJoystickIdentity
struct BrickletJoystickIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletJoystick
"""
2-axis joystick with push-button
"""
mutable struct BrickletJoystick <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletJoystick(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_joystick")
        deviceInternal = package.BrickletJoystick(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_position
"""
    $(SIGNATURES)

Returns the position of the joystick. The middle position of the joystick is x=0, y=0.
The returned values are averaged and calibrated (see :func:`Calibrate`).

If you want to get the position periodically, it is recommended to use the
:cb:`Position` callback and set the period with
:func:`Set Position Callback Period`.
"""
function get_position(device::BrickletJoystick)
    return device.deviceInternal.get_position()
end

export is_pressed
"""
    $(SIGNATURES)

Returns *true* if the button is pressed and *false* otherwise.

It is recommended to use the :cb:`Pressed` and :cb:`Released` callbacks
to handle the button.
"""
function is_pressed(device::BrickletJoystick)
    return device.deviceInternal.is_pressed()
end

export get_analog_value
"""
    $(SIGNATURES)

Returns the values as read by a 12-bit analog-to-digital converter.

.. note::
 The values returned by :func:`Get Position` are averaged over several samples
 to yield less noise, while :func:`Get Analog Value` gives back raw
 unfiltered analog values. The only reason to use :func:`Get Analog Value` is,
 if you need the full resolution of the analog-to-digital converter.

If you want the analog values periodically, it is recommended to use the
:cb:`Analog Value` callback and set the period with
:func:`Set Analog Value Callback Period`.
"""
function get_analog_value(device::BrickletJoystick)
    return device.deviceInternal.get_analog_value()
end

export calibrate
"""
    $(SIGNATURES)

Calibrates the middle position of the joystick. If your Joystick Bricklet
does not return x=0 and y=0 in the middle position, call this function
while the joystick is standing still in the middle position.

The resulting calibration will be saved on the EEPROM of the Joystick
Bricklet, thus you only have to calibrate it once.
"""
function calibrate(device::BrickletJoystick)
    device.deviceInternal.calibrate()
end

export set_position_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Position` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Position` callback is only triggered if the position has changed since the
last triggering.
"""
function set_position_callback_period(device::BrickletJoystick, period)
    device.deviceInternal.set_position_callback_period(period)
end

export get_position_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Position Callback Period`.
"""
function get_position_callback_period(device::BrickletJoystick)
    return device.deviceInternal.get_position_callback_period()
end

export set_analog_value_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Analog Value` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Analog Value` callback is only triggered if the analog values have
changed since the last triggering.
"""
function set_analog_value_callback_period(device::BrickletJoystick, period)
    device.deviceInternal.set_analog_value_callback_period(period)
end

export get_analog_value_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Analog Value Callback Period`.
"""
function get_analog_value_callback_period(device::BrickletJoystick)
    return device.deviceInternal.get_analog_value_callback_period()
end

export set_position_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Position Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the position is *outside* the min and max values"
 "'i'",    "Callback is triggered when the position is *inside* the min and max values"
 "'<'",    "Callback is triggered when the position is smaller than the min values (max is ignored)"
 "'>'",    "Callback is triggered when the position is greater than the min values (max is ignored)"
"""
function set_position_callback_threshold(device::BrickletJoystick, option, min_x, max_x, min_y, max_y)
    device.deviceInternal.set_position_callback_threshold(option, min_x, max_x, min_y, max_y)
end

export get_position_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Position Callback Threshold`.
"""
function get_position_callback_threshold(device::BrickletJoystick)
    return device.deviceInternal.get_position_callback_threshold()
end

export set_analog_value_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Analog Value Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the analog values are *outside* the min and max values"
 "'i'",    "Callback is triggered when the analog values are *inside* the min and max values"
 "'<'",    "Callback is triggered when the analog values are smaller than the min values (max is ignored)"
 "'>'",    "Callback is triggered when the analog values are greater than the min values (max is ignored)"
"""
function set_analog_value_callback_threshold(device::BrickletJoystick, option, min_x, max_x, min_y, max_y)
    device.deviceInternal.set_analog_value_callback_threshold(option, min_x, max_x, min_y, max_y)
end

export get_analog_value_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Analog Value Callback Threshold`.
"""
function get_analog_value_callback_threshold(device::BrickletJoystick)
    return device.deviceInternal.get_analog_value_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Position Reached`,
* :cb:`Analog Value Reached`

are triggered, if the thresholds

* :func:`Set Position Callback Threshold`,
* :func:`Set Analog Value Callback Threshold`

keep being reached.
"""
function set_debounce_period(device::BrickletJoystick, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletJoystick)
    return device.deviceInternal.get_debounce_period()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletJoystick)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletJoystick, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
module JuliaToolsForMagTIP
# pls add https://github.com/CGRG-lab/DataFrameTools.jl.git, https://github.com/CGRG-lab/FileTools.jl.git

# Write your package code here.
using Dates
include("datenum.jl")
export datenum, toordinal, chkdatetime

include("isgeomagdatavalid.jl")
export isgeomagdatavalid

end


text:
using CrossSections
using Test

using Unitful: MeV
using Random
rng = MersenneTwister()
proc = ComptonKleinNishina()

for E in [0.1MeV, 1MeV, 10MeV, 100MeV]
    for _ in 1:10
        ret = sample(rng, proc, E)
        @show E
        @show ret
    end
end


text:
module LambdaCases
using MLStyle
using MLStyle.Sugars
using MLStyle.AbstractPatterns: init_cfg

export gen_lambda, @λ
@nospecialize
"""
Code generation for `@λ`.
The first argument must be something like
- `a -> b`
- `begin a -> b; (c -> d)... end`
"""
function gen_lambda(cases, source :: LineNumberNode, mod :: Module)
    TARGET = gensym("λ")
    function make_pair_expr(case, stmts)
        let block = Expr(:block, stmts...)
            :($case => $block)
        end
    end

    @switch cases begin
        @case :($a -> $(bs...))     ||
              :($a => $b) && let bs = [b] end

            pair = make_pair_expr(a, bs)
            cbl = Expr(:block, source, pair)
            match_expr = gen_match(TARGET, cbl, source, mod)

        @case let stmts=[] end && Expr(:block,
            Many[
                a :: LineNumberNode && Do[push!(stmts, a)] ||
                Or[:($a => $b) && let bs=[b] end, :($a -> $(bs...))] &&
                Do[push!(stmts, make_pair_expr(a, bs))]
            ]...
        )

            cbl = Expr(:block, source, stmts...)
            match_expr = gen_match(TARGET, cbl, source, mod)
    end

    Expr(:function,
        Expr(:call, TARGET, TARGET),
        Expr(:block, source, init_cfg(match_expr))
    )
end

"""
Lambda cases.

e.g.

```julia
    xs = [(1, 2), (1, 3), (1, 4)]
    map((@λ (1, x) => x), xs)
    # => [2, 3, 4]

    (2, 3) |> @λ begin
        1 => 2
        2 => 7
        (a, b) => a + b
    end
    # => 5
```
"""
macro λ(cases)
    gen_lambda(cases, __source__, __module__) |> esc
end
@specialize

end


text:
#-----------------------------------------------------------------------------------------------
# 1st run - packages installation

# ERROR: LoadError: ArgumentError: Package DataFrames not found in current path:
# - Run `import Pkg; Pkg.add("DataFrames")` to install the DataFrames package.
# import Pkg; Pkg.add("DataFrames")
# import Pkg;
# Pkg.add("CSV")
# Pkg.add("StatsBase")

#-----------------------------------------------------------------------------------------------
using DataFrames, CSV, Statistics, StatsBase
data = CSV.read("../../../externals/Core.Math.Samples/data/Pejcic_318.csv", copycols = true)

# @show data.ATV
# println()

for col in eachcol(data,true);
    if  (
        eltype(col[2]) <: Real)
        println(mean(col[2])
        );
    end;
end;

# println(names(data),"\n")

# println("Sample Mean: ", mean(data[:ATT]))
# println("Harmonic <= Geometric <= Arithmetic ",(harmmean(data), geomean(data), mean(data)))
# println("Sample Variance: ",var(data))
# println("Sample Standard Deviation: ",std(data))
# println("Minimum: ", minimum(data))
# println("Maximum: ", maximum(data))
# println("Median: ", median(data))
# println("95th percentile: ", percentile(data, 95))
# println("0.95 quantile: ", quantile(data, 0.95))
# println("Interquartile range: ", iqr(data),"\n")


text:
include("BaseEnumeration.jl")
include("Simplexe/main.jl")
global basesolver = BaseEnumeration
function changeSolver!(::Type{T}) where T <: AbstractSolver
    global basesolver = T
end
function solve!(lp::AbstractLP; verbose::Bool = false)
    solver = basesolver(lp)
    solver(lp, verbose = verbose)
    push!(lp.solvers, solver)
    return lp.status
end

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin


text:
@testset "linear problems - partial dualization" begin
    @testset "lp1_test - ignore x_2" begin
        #=
        primal
            min -4x_2 - 1
        s.t.
            x_1 >= 3         :y_2
            x_1 + 2x_2 <= 3  :y_3

        ignore x_2 during dualization

        dual
            obj ignored
        s.t.
            y_2 >= 0
            y_3 <= 0
            y_2 + y_3 == 0    :x_1
        =#
        primal_model = lp1_test()
        dual = Dualization.dualize(
            primal_model,
            variable_parameters = VI[VI(2)],
            ignore_objective = true,
        )
        dual_model = dual.dual_model
        primal_dual_map = dual.primal_dual_map

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 2
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set(
            [
                (SVF, MOI.GreaterThan{Float64})
                (SVF, MOI.LessThan{Float64})
                (SAF{Float64}, MOI.EqualTo{Float64})
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SVF,MOI.GreaterThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SVF,MOI.LessThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 1
        # @show obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())

        eq_con1_fun = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(1),
        )
        eq_con1_set = MOI.get(
            dual_model,
            MOI.ConstraintSet(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(1),
        )
        @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con1_fun) == 0.0
        @test MOI.constant(eq_con1_set) == 0.0

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{SAF{Float64},MOI.LessThan{Float64}}(1)] ==
              [VI(1)]

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(1)] ==
              CI{SAF{Float64},MOI.EqualTo{Float64}}(1)
    end

    @testset "lp7_test - x_1 ignored" begin
        #=
        primal
            min -4x1 -3x2 -1
        s.t.
            2x1 + x2 - 3 <= 0  :y_3
            x1 + 2x2 - 3 <= 0  :y_4
            x1 >= 1            :y_1
            x2 >= 0            :y_2

            ignore x_1 during dualization

        dual
            obj ignored
        s.a.
            y_2 + y_3 + 2y_4 = -3 :x_2
            y_1 >= 0
            y_2 >= 0
            y_3 <= 0
            y_4 <= 0
        =#
        primal_model = lp7_test()
        dual = Dualization.dualize(
            primal_model,
            variable_parameters = VI[VI(1)],
            ignore_objective = true,
        )
        dual_model = dual.dual_model
        primal_dual_map = dual.primal_dual_map

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 4
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set(
            [
                (SVF, MOI.GreaterThan{Float64})
                (SAF{Float64}, MOI.EqualTo{Float64})
                (VVF, MOI.Nonpositives)
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SVF,MOI.GreaterThan{Float64}}(),
        ) == 2
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VVF,MOI.Nonpositives}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 1
        # @show obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())

        eq_con2_fun = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(2),
        )
        eq_con2_set = MOI.get(
            dual_model,
            MOI.ConstraintSet(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(2),
        )
        @test MOI.coefficient.(eq_con2_fun.terms) == [1.0; 2.0; 1.0]
        @test MOI.constant.(eq_con2_fun) == 0.0
        @test MOI.constant(eq_con2_set) == -3.0

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{VAF{Float64},MOI.Nonpositives}(1)] ==
              [VI(1); VI(2)]
        @test primal_con_dual_var[CI{SVF,MOI.GreaterThan{Float64}}(1)] ==
              [VI(3)]
        @test primal_con_dual_var[CI{SVF,MOI.GreaterThan{Float64}}(2)] ==
              [VI(4)]

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(2)] ==
              CI{SAF{Float64},MOI.EqualTo{Float64}}(2)
    end

    @testset "lp12_test - x_1 and x_3 ignored" begin
        #=
        primal
            min 4x_3 + 5
        s.t.
            x_1 + 2x_2 + x_3 <= 20 :y_3
            x_1 <= 1               :y_1
            x_2 <= 3               :y_2

        ignoring x_1 and x_3

        dual
            obj ignored
        s.t.
            #  y_1 + y_3 == 0  :x_1
            y_2 + 2y_3 == 0 :x_2
            #  y_3 == 4.0      :x_3
            y_1 <= 0
            y_2 <= 0
            y_3 <= 0
        =#
        primal_model = lp12_test()
        dual = Dualization.dualize(
            primal_model,
            variable_parameters = VI[VI(1), VI(3)],
            ignore_objective = true,
        )
        dual_model = dual.dual_model
        primal_dual_map = dual.primal_dual_map

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 3
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set(
            [
                (SVF, MOI.LessThan{Float64})
                (SAF{Float64}, MOI.EqualTo{Float64})
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SVF,MOI.LessThan{Float64}}(),
        ) == 3
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 1

        eq_con2_fun = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(1),
        )
        eq_con2_set = MOI.get(
            dual_model,
            MOI.ConstraintSet(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(1),
        )
        @test MOI.coefficient.(eq_con2_fun.terms) == [2.0; 1.0]
        @test MOI.constant.(eq_con2_fun) == 0.0
        @test MOI.constant(eq_con2_set) == 0.0

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{SVF,MOI.LessThan{Float64}}(2)] == [VI(3)]
        @test primal_con_dual_var[CI{SVF,MOI.LessThan{Float64}}(1)] == [VI(2)]
        @test primal_con_dual_var[CI{SAF{Float64},MOI.LessThan{Float64}}(1)] ==
              [VI(1)]

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(2)] ==
              CI{SAF{Float64},MOI.EqualTo{Float64}}(1)
    end
end


text:
bitstype 256 Bits256

immutable BigFloatB <: AbstractFloat
    prec::Clong
    sign::Cint
    exp::Clong
    signif::Bits256
end

precision(x::BigFloatB) = x.prec

unsafe_convert(::Type{Ptr{UInt64}}, bref::Ptr{BigFloatX.Bits256}) = unsafe_convert(Ptr{Limb},unsafe_convert(Ptr{Void},bref))
unsafe_convert(::Type{Ptr{Limb}}, bref::Ref{Bits256}) = unsafe_convert(Ptr{Limb},unsafe_convert(Ptr{Void},bref))

function new_mpfr_bref()
    prec = precision(BigFloat)
    bref = Ref{Bits256}()
    Mpfr_t(prec,0,0,unsafe_convert(Ptr{Limb}, bref)), bref
end

BigFloatB(m::Mpfr_t, s::Ref{Bits256}) =
    BigFloatB(m.prec, m.sign, m.exp, s[])

convert(::Type{Mpfr_t}, x::BigFloatB) =
    Mpfr_t(x.prec, x.sign, x.exp, unsafe_convert(Ptr{Limb},Ref(x.signif)))

for (fJ, fC) in ((:si,:Clong), (:ui,:Culong), (:d,:Float64))
    @eval begin
        function convert(::Type{BigFloatB}, x::($fC))
            m, bref = new_mpfr_bref()
            ccall(($(string(:mpfr_set_,fJ)), :libmpfr), Int32, (Ptr{Mpfr_t}, ($fC), Int32), &m, x, ROUNDING_MODE[end])
            return BigFloatB(m, bref)
        end
    end
end

# Basic arithmetic without promotion
for (fJ, fC) in ((:+,:add), (:*,:mul))
    @eval begin
        # BigFloatB
        function ($fJ)(x::BigFloatB, y::BigFloatB)
            m, bref = new_mpfr_bref()
            ccall(($(string(:mpfr_,fC)),:libmpfr), Int32, (Ptr{Mpfr_t}, Ptr{Mpfr_t}, Ptr{Mpfr_t}, Int32), &m, &x, &y, ROUNDING_MODE[end])
            return BigFloatB(m, bref)
        end

    end
end

zero(::Type{BigFloatB}) = convert(BigFloatB,0)

function string(x::BigFloatB)
    # In general, the number of decimal places needed to read back the number exactly
    # is, excluding the most significant, ceil(log(10, 2^precision(x)))
    k = ceil(Int32, precision(x) * 0.3010299956639812)
    lng = k + Int32(8) # Add space for the sign, the most significand digit, the dot and the exponent
    buf = Array(UInt8, lng + 1)
    # format strings are guaranteed to contain no NUL, so we don't use Cstring
    lng = ccall((:mpfr_snprintf,:libmpfr), Int32, (Ptr{UInt8}, Culong, Ptr{UInt8}, Ptr{Mpfr_t}), buf, lng + 1, "%.Re", &x)
    if lng < k + 5 # print at least k decimal places
        lng = ccall((:mpfr_sprintf,:libmpfr), Int32, (Ptr{UInt8}, Ptr{UInt8}, Ptr{Mpfr_t}), buf, "%.$(k)Re", &x)
    elseif lng > k + 8
        buf = Array(UInt8, lng + 1)
        lng = ccall((:mpfr_snprintf,:libmpfr), Int32, (Ptr{UInt8}, Culong, Ptr{UInt8}, Ptr{Mpfr_t}), buf, lng + 1, "%.Re", &x)
    end
    return bytestring(pointer(buf), (1 <= x < 10 || -10 < x <= -1 || x == 0) ? lng - 4 : lng)
end


text:
"""
$(SIGNATURES)

Compute the convex hull of `points` and store the result in `hull` using the
[Jarvis march (gift wrapping) algorithm](https://en.wikipedia.org/wiki/Gift_wrapping_algorithm).
This algorithm has ``O(nh)`` complexity, where ``n`` is the number of points and ``h`` is
the number of vertices of the convex hull.
"""
function jarvis_march!(hull::ConvexHull, points::AbstractVector{<:PointLike}; atol=eps(eltype(hull)))
    # Adapted from https://www.algorithm-archive.org/contents/jarvis_march/jarvis_march.html.
    op = orientation_comparator(hull)
    n = length(points)
    vertices = hull.vertices
    @inbounds begin
        # Preallocate
        resize!(vertices, n)
        if n <= 2
            vertices .= points
        else
            # Find an initial hull vertex using lexicographic ordering.
            current = first(points)
            for i in 2 : n
                p = points[i]
                if Tuple(p) < Tuple(current)
                    current = p
                end
            end

            i = 1
            while true
                # Add point
                vertices[i] = current

                # Next point is the one with extremal internal angle.
                next = first(points)
                δnext = next - current
                for j in 2 : n
                    p = points[j]
                    δ = p - current
                    c = cross2(δnext, δ)

                    # Note the last clause here, which ensures strong convexity in the presence of
                    # collinear points by accepting `p` if it's farther away from `current` than
                    # `next`.
                    if next == current || (op(0, c) && abs(c) > atol) || (abs(c) <= atol && δ ⋅ δ > δnext ⋅ δnext)
                        next = p
                        δnext = δ
                    end
                end
                current = next
                current == first(vertices) && break
                i += 1
                if i > n
                    # println(IOContext(stdout, :compact=>false), points)
                    # println(IOContext(stdout, :compact=>false), vertices)
                    error("should never get here.")
                end
            end

            # Shrink to computed number of vertices.
            resize!(vertices, i)
        end
    end
    return hull
end


text:
# Proximal methods for LASSO

using Regression
import Regression: Options, solve

## prepare for simulation

d = 100
num_zeros = 60
n = 5000

w = randn(d)
iz = sort(randperm(d)[1:num_zeros])
w[iz] = 0.0
@assert countnz(w) == d - num_zeros

X = randn(d, n)
y = X'w + 0.1 * randn(n)

## solve w using proximal gradient descent (i.e. ISTA)

ret = solve(linearreg(X, y);
            reg=L1Reg(0.01 * n),
            solver=ProxGD(),
            options=Options(verbosity=:iter, armijo=0.2, ftol=1.0e-8 * n))

println()

## examine solution

w_e = ret.sol
@printf("relative error   = %.4e\n", sumabs2(w_e - w) / sumabs2(w))

numz_w = d - countnz(w)
numz_we = d - countnz(w_e)
@printf("w:   num_zeros = %d\n", numz_w)
@printf("w_e: num_zeros = %d\n", numz_we)
@printf("overlap of zeros = %.1f%%\n",
    100.0 * countnz((w .== 0) & (w_e .== 0)) / max(numz_w, numz_we))


text:
module BadelineBoss

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Badeline Boss" => Ahorn.EntityPlacement(
        Maple.BadelineBoss,
    )
)

Ahorn.editingOptions(entity::Maple.BadelineBoss) = Dict{String, Any}(
    "patternIndex" => Maple.badeline_boss_shooting_patterns
)

Ahorn.nodeLimits(entity::Maple.BadelineBoss) = 0, -1

sprite = "characters/badelineBoss/charge00.png"

function Ahorn.selection(entity::Maple.BadelineBoss)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.getSpriteRectangle(sprite, nx, ny))
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.BadelineBoss)
    px, py = Ahorn.position(entity)

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)

        Ahorn.drawArrow(ctx, px, py, nx, ny, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny)

        px, py = nx, ny
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.BadelineBoss, room::Maple.Room)
    x, y = Ahorn.position(entity)
    Ahorn.drawSprite(ctx, sprite, x, y)
end

end


text:
using Conda

Conda.add("pymatgen", channel = "conda-forge")


text:
module Unitary

# unitary U
# tangent vectors Δ = U*A with A' = -A

using TensorKit
import TensorKit: similarstoragetype, SectorDict
using ..TensorKitManifolds: projectantihermitian!, projectisometric!, PolarNewton
import ..TensorKitManifolds: base, checkbase, inner, retract, transport, transport!

struct UnitaryTangent{T<:AbstractTensorMap, TA<:AbstractTensorMap}
    W::T
    A::TA
    function UnitaryTangent(W::AbstractTensorMap{S,N₁,N₂},
                            A::AbstractTensorMap{S,N₂,N₂}) where {S,N₁,N₂}
        T = typeof(W)
        TA = typeof(A)
        return new{T,TA}(W,A)
    end
end
Base.copy(Δ::UnitaryTangent) = UnitaryTangent(Δ.W, copy(Δ.A))
Base.getindex(Δ::UnitaryTangent) = Δ.W * Δ.A
base(Δ::UnitaryTangent) = Δ.W
checkbase(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) = Δ₁.W == Δ₂.W ? Δ₁.W :
    throw(ArgumentError("tangent vectors with different base points"))

# Basic vector space behaviour
Base.:+(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) =
    UnitaryTangent(checkbase(Δ₁,Δ₂), Δ₁.A + Δ₂.A)
Base.:-(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) =
    UnitaryTangent(checkbase(Δ₁,Δ₂), Δ₁.A - Δ₂.A)
Base.:-(Δ::UnitaryTangent) = (-1)*Δ

Base.:*(Δ::UnitaryTangent, α::Real) = rmul!(copy(Δ), α)
Base.:*(α::Real, Δ::UnitaryTangent) = lmul!(α, copy(Δ))
Base.:/(Δ::UnitaryTangent, α::Real) = rmul!(copy(Δ), inv(α))
Base.:\(α::Real, Δ::UnitaryTangent) = lmul!(inv(α), copy(Δ))

Base.zero(Δ::UnitaryTangent) = UnitaryTangent(Δ.W, zero(Δ.A))

function TensorKit.rmul!(Δ::UnitaryTangent, α::Real)
    rmul!(Δ.A, α)
    return Δ
end
function TensorKit.lmul!(α::Real, Δ::UnitaryTangent)
    lmul!(α, Δ.A)
    return Δ
end
function TensorKit.axpy!(α::Real, Δx::UnitaryTangent, Δy::UnitaryTangent)
    checkbase(Δx, Δy)
    axpy!(α, Δx.A, Δy.A)
    return Δy
end
function TensorKit.axpby!(α::Real, Δx::UnitaryTangent, β::Real, Δy::UnitaryTangent)
    checkbase(Δx, Δy)
    axpby!(α, Δx.A, β, Δy.A)
    return Δy
end

function TensorKit.dot(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent)
    checkbase(Δ₁, Δ₂)
    return dot(Δ₁.A, Δ₂.A)
end
TensorKit.norm(Δ::UnitaryTangent, p::Real = 2) = norm(Δ.A, p)

# tangent space methods
function inner(W::AbstractTensorMap, Δ₁::UnitaryTangent, Δ₂::UnitaryTangent;
                metric = :euclidean)
    @assert metric == :euclidean
    Δ₁ === Δ₂ ? norm(Δ₁)^2 : real(dot(Δ₁,Δ₂))
end
function project!(X::AbstractTensorMap, W::AbstractTensorMap; metric = :euclidean)
    @assert metric == :euclidean
    P = W'*X
    A = projectantihermitian!(P)
    return UnitaryTangent(W, A)
end
project(X, W; metric = :euclidean) = project!(copy(X), W; metric = :euclidean)

# geodesic retraction, coincides with Stiefel retraction (which is not geodesic for p < n)
function retract(W::AbstractTensorMap, Δ::UnitaryTangent, α; alg = nothing)
    W == base(Δ) || throw(ArgumentError("not a valid tangent vector at base point"))
    E = exp(α*Δ.A)
    W′ = projectisometric!(W*E; alg = SDD())
    A′ = Δ.A
    return W′, UnitaryTangent(W′, A′)
end

# isometric vector transports compatible with above retraction (also with differential of retraction)
function transport!(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α::Real, W′;
                    alg = :stiefel)
    if alg == :parallel
        return transport_parallel!(Θ, W, Δ, α, W′)
    elseif alg == :stiefel
        return transport_stiefel!(Θ, W, Δ, α, W′)
    else
        throw(ArgumentError("unknown algorithm: `alg = $metric`"))
    end
end
function transport(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α::Real, W′;
                    alg = :stiefel)
    return transport!(copy(Θ), W, Δ, α, W′; alg = alg)
end

# transport_parallel correspondings to the torsion-free Levi-Civita connection
# transport_stiefel is compatible to Stiefel.transport and corresponds to a non-torsion-free connection
function transport_parallel!(Θ::UnitaryTangent,
                                W::AbstractTensorMap,
                                Δ::UnitaryTangent, α, W′)
    W == checkbase(Δ,Θ) || throw(ArgumentError("not a valid tangent vector at base point"))
    E = exp((α/2)*Δ.A)
    A′ = projectantihermitian!(E'*Θ.A*E) # exra projection for stability
    return UnitaryTangent(W′, A′)
end
transport_parallel(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α, W′) =
    transport_parallel!(copy(Θ), W, Δ, α, W′)

function transport_stiefel!(Θ::UnitaryTangent,
                                W::AbstractTensorMap,
                                Δ::UnitaryTangent, α, W′)
    W == checkbase(Δ,Θ) || throw(ArgumentError("not a valid tangent vector at base point"))
    A′ = Θ.A
    return UnitaryTangent(W′, A′)
end
transport_stiefel(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α, W′) =
    transport_stiefel!(copy(Θ), W, Δ, α, W′)

end


text:
using DashBootstrapComponents

badges = html_span([
    dbc_badge(
        "Primary",
        href = "#",
        color = "primary",
        className = "me-1 text-decoration-none",
    ),
    dbc_badge(
        "Secondary",
        href = "#",
        color = "secondary",
        className = "me-1 text-decoration-none",
    ),
    dbc_badge(
        "Success",
        href = "#",
        color = "success",
        className = "me-1 text-decoration-none",
    ),
    dbc_badge(
        "Warning",
        href = "#",
        color = "warning",
        className = "me-1 text-decoration-none",
    ),
    dbc_badge(
        "Danger",
        href = "#",
        color = "danger",
        className = "me-1 text-decoration-none",
    ),
    dbc_badge("Info", href = "#", color = "info", className = "me-1 text-decoration-none"),
    dbc_badge(
        "Light",
        href = "#",
        color = "light",
        text_color = "dark",
        className = "me-1 text-decoration-none",
    ),
    dbc_badge("Dark", href = "#", color = "dark", className = "text-decoration-none"),
],);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module ParetoFront

export update_pareto!,  get_pareto

# We have two n-vectors x and y.  min_idxs is a vector of indices at which
# x[i] < y[i], if x is to dominate y. Likewise, max_idxs is a vector of
# indices at which x[i] > y[i], if x is to dominate y
function domvals(x, y, min_idxs, max_idxs)
    n_better = 0
    n_worse  = 0

    for i in min_idxs
        n_better += x[i] < y[i]
        n_worse  += x[i] > y[i]
    end

    for i in max_idxs
        n_better += x[i] > y[i]
        n_worse  += x[i] < y[i]
    end

    return n_better, n_worse
end

# Conditionally update the pareto set (a set of n-vectors) with a new candidate
# vector x. The vector min_idxs contains the indices of the vector x whose
# values we want to find the minimum and analogously for the vector max_idxs.
# Note that if an index does not appear in either min_idxs or max_idxs, it will
# not be used in calculating the Pareto front
function update_pareto!(pareto, x, min_idxs, max_idxs)
    remove_pts = Set()
    for p in pareto
        n_better, n_worse = domvals(x, p, min_idxs, max_idxs)

        # Current candidate is dominated by point j in pareto so ignore it
        # To avoid including duplicate points, change "n_worse > 0" to "n_worse >= 0"
        n_worse > 0 && n_better == 0 && return pareto

        # Current candidate dominates point p in pareto, so remove p from pareto
        n_better > 0 && n_worse == 0 && push!(remove_pts, p)
    end

    setdiff!(pareto, remove_pts)
    push!(pareto, x)
    return pareto
end

# Use get_pareto when you already have a matrix (whose rows are points)
# out of which you want to extract the Pareto front
function get_pareto(V::Matrix{T}; min_idxs, max_idxs) where {T}
    pareto = Set{Vector{T}}()
    @inbounds for x in eachrow(V)
        update_pareto!(pareto, x, min_idxs, max_idxs)
    end
    return pareto
end

end


text:
module BenchCat

using BenchmarkTools
using Transducers

itr = (y for x in 1:1000 for y in 1:x);

suite = BenchmarkGroup()

suite["xf"] = @benchmarkable foldl(+, eduction($itr))
suite["base"] = @benchmarkable sum($itr)

end  # module
BenchCat.suite


text:
include("1.jl")
include("2.jl")
function updowntwicestep(r::Robot,side1::HorizonSide)
    godownost(r,Sud,West)
    putmarker!(r)
    putmarkhole(r,side1)
    while !isborder(r,Nord)
        move!(r,Nord)
        putmarkhole(r,inverse(side1))
        backing(r, side1)
        move!(r,inverse(side1))
    end
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
"""
    flux_boundary_maxwell!(fw, bc, w, inK, γ, dt, rot)
    flux_boundary_maxwell!(fw, bc, w, inK, γ, dt, len, rot)
    flux_boundary_maxwell!(fw, fh, fb, bc, h, b, u, ω, inK, dt, rot)
    flux_boundary_maxwell!(fw, ff, bc, f, u, v, ω, dt, len, rot)
    flux_boundary_maxwell!(fw, fh, fb, bc, h, b, u, v, ω, inK, dt, len, rot)
    flux_boundary_maxwell!(fw, ff, bc, f, u, v, w, ω, dt, area, rot)

Maxwell's diffusive boundary flux

- @args: distribution functions and their slopes at left/right sides of interface
- @args: particle velocity quadrature points and weights
- @args: time step

"""
function flux_boundary_maxwell!(
    fw::AbstractVector{T1},
    bc::AbstractVector{T2},
    w::AbstractVector{T3},
    inK,
    γ,
    dt,
    rot, # 1 / -1
) where {
    T1<:AbstractFloat,
    T2<:Real,
    T3<:Real,
} # 1D continuum

    @assert length(bc) == 3

    primL, primR = ifelse(rot == 1, (bc, conserve_prim(w, γ)), (conserve_prim(w, γ), bc))

    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    Muv_L = moments_conserve(MuL1, Mxi1, 1, 0)
    Muv_R = moments_conserve(MuR2, Mxi2, 1, 0)

    ρ = ifelse(rot == 1, -primR[1] * Muv_R[1] / Muv_L[1], -primL[1] * Muv_L[1] / Muv_R[1])

    @. fw = ρ * (Muv_L + Muv_R) * dt

    return nothing

end

#--- 1D2F1V ---#
function flux_boundary_maxwell!(
    fw::AbstractVector{T1},
    fh::T2,
    fb::T2,
    bc::AbstractVector{T3},
    h::T4,
    b::T4,
    u::T5,
    ω::T5,
    inK,
    dt,
    rot = 1,
) where {
    T1<:AbstractFloat,
    T2<:AbstractVector{<:AbstractFloat},
    T3<:Real,
    T4<:AbstractVector{<:AbstractFloat},
    T5<:AbstractVector{<:AbstractFloat},
} # 1D2F1V

    @assert length(bc) == 3

    δ = heaviside.(u .* rot)
    SF = sum(ω .* u .* h .* (1.0 .- δ))
    SG = (bc[end] / π)^0.5 * sum(ω .* u .* exp.(-bc[end] .* (u .- bc[2]) .^ 2) .* δ)
    prim = [-SF / SG; bc[2:end]]

    H = maxwellian(u, prim)
    B = H .* inK ./ (2.0 * prim[end])

    hWall = H .* δ .+ h .* (1.0 .- δ)
    bWall = B .* δ .+ b .* (1.0 .- δ)

    fw[1] = discrete_moments(hWall, u, ω, 1) * dt
    fw[2] = discrete_moments(hWall, u, ω, 2) * dt
    fw[3] =
        (
            0.5 * discrete_moments(hWall .* u .^ 2, u, ω, 1) +
            0.5 * discrete_moments(bWall, u, ω, 1)
        ) * dt

    @. fh = u * hWall * dt
    @. fb = u * bWall * dt

    return nothing

end

#--- 2D Continuum ---#
function flux_boundary_maxwell!(
    fw::AbstractVector{T1},
    bc::AbstractVector{T2},
    w::AbstractVector{T3},
    inK,
    γ,
    dt,
    len,
    rot,
) where {
    T1<:Real,
    T2<:Real,
    T3<:Real,
}

    @assert length(bc) == 4

    primL, primR = ifelse(rot == 1, (bc, conserve_prim(w, γ)), (conserve_prim(w, γ), bc))

    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    Muv_L = moments_conserve(MuL1, Mv1, Mxi1, 1, 0, 0)
    Muv_R = moments_conserve(MuR2, Mv2, Mxi2, 1, 0, 0)

    ρ = ifelse(rot == 1, -primR[1] * Muv_R[1] / Muv_L[1], -primL[1] * Muv_L[1] / Muv_R[1])

    @. fw = ρ * (Muv_L + Muv_R) * dt * len

    return nothing

end

# ------------------------------------------------------------
# 2D1F2V
# ------------------------------------------------------------
function flux_boundary_maxwell!(
    fw::AbstractVector{T1},
    fh::AbstractMatrix{T2},
    bc::AbstractVector{T3},
    h::AbstractMatrix{T4},
    u::T5,
    v::T5,
    ω::T5,
    dt,
    len,
    rot = 1,
) where {
    T1<:AbstractFloat,
    T2<:AbstractFloat,
    T3<:Real,
    T4<:AbstractFloat,
    T5<:AbstractMatrix{<:AbstractFloat},
}

    @assert length(bc) == 4

    δ = heaviside.(u .* rot)

    SF = sum(ω .* u .* h .* (1.0 .- δ))
    SG =
        (bc[end] / π) *
        sum(ω .* u .* exp.(-bc[end] .* ((u .- bc[2]) .^ 2 .+ (v .- bc[3]) .^ 2)) .* δ)
    prim = [-SF / SG; bc[2:end]]

    H = maxwellian(u, v, prim)

    hWall = H .* δ .+ h .* (1.0 .- δ)

    fw[1] = discrete_moments(hWall, u, ω, 1) * len * dt
    fw[2] = discrete_moments(hWall, u, ω, 2) * len * dt
    fw[3] = discrete_moments(hWall .* u, v, ω, 1) * len * dt
    fw[4] = 0.5 * discrete_moments(hWall .* (u .^ 2 .+ v .^ 2), u, ω, 1) * len * dt

    @. fh = u * hWall * len * dt

    return nothing

end

# ------------------------------------------------------------
# 2D2F2V
# ------------------------------------------------------------
function flux_boundary_maxwell!(
    fw::AbstractVector{T1},
    fh::T2,
    fb::T2,
    bc::AbstractVector{T3},
    h::T4,
    b::T4,
    u::T5,
    v::T5,
    ω::T5,
    inK,
    dt,
    len,
    rot = 1,
) where {
    T1<:AbstractFloat,
    T2<:AbstractMatrix{<:AbstractFloat},
    T3<:Real,
    T4<:AbstractMatrix{<:AbstractFloat},
    T5<:AbstractMatrix{<:AbstractFloat},
}

    @assert length(bc) == 4

    δ = heaviside.(u .* rot)

    SF = sum(ω .* u .* h .* (1.0 .- δ))
    SG =
        (bc[end] / π) *
        sum(ω .* u .* exp.(-bc[end] .* ((u .- bc[2]) .^ 2 .+ (v .- bc[3]) .^ 2)) .* δ)
    prim = [-SF / SG; bc[2:end]]

    H = maxwellian(u, v, prim)
    B = H .* inK ./ (2.0 * prim[end])

    hWall = H .* δ .+ h .* (1.0 .- δ)
    bWall = B .* δ .+ b .* (1.0 .- δ)

    fw[1] = discrete_moments(hWall, u, ω, 1) * len * dt
    fw[2] = discrete_moments(hWall, u, ω, 2) * len * dt
    fw[3] = discrete_moments(hWall .* u, v, ω, 1) * len * dt
    fw[4] =
        (
            0.5 * discrete_moments(hWall .* (u .^ 2 .+ v .^ 2), u, ω, 1) +
            0.5 * discrete_moments(bWall, u, ω, 1)
        ) *
        len *
        dt

    @. fh = u * hWall * len * dt
    @. fb = u * bWall * len * dt

    return nothing

end

# ------------------------------------------------------------
# 1F3V
# ------------------------------------------------------------
function flux_boundary_maxwell!(
    fw::AbstractVector{T1},
    ff::AbstractArray{T2,3},
    bc::AbstractVector{T3},
    f::AbstractArray{T4,3},
    u::T5,
    v::T5,
    w::T5,
    ω::T5,
    dt,
    area,
    rot = 1,
) where {
    T1<:AbstractFloat,
    T2<:AbstractFloat,
    T3<:Real,
    T4<:AbstractFloat,
    T5<:AbstractArray{<:AbstractFloat,3},
}

    @assert length(bc) == 5

    δ = heaviside.(u .* rot)
    SF = sum(ω .* u .* f .* (1.0 .- δ))
    SG =
        (bc[end] / π)^1.5 * sum(
            ω .* u .*
            exp.(
                -bc[end] .* ((u .- bc[2]) .^ 2 .+ (v .- bc[3]) .^ 2 .+ (w .- bc[4]) .^ 2),
            ) .* δ,
        )
    prim = [-SF / SG; bc[2:end]]

    M = maxwellian(u, v, w, prim)
    fWall = M .* δ .+ f .* (1.0 .- δ)

    fw[1] = discrete_moments(fWall, u, ω, 1) * area * dt
    fw[2] = discrete_moments(fWall, u, ω, 2) * area * dt
    fw[3] = discrete_moments(fWall .* u, v, ω, 1) * area * dt
    fw[4] = discrete_moments(fWall .* u, w, ω, 1) * area * dt
    fw[5] =
        (0.5 * discrete_moments(fWall .* (u .^ 2 .+ v .^ 2 + w .^ 2), u, ω, 1)) * area * dt

    @. ff = u * fWall * area * dt

    return nothing

end

"""
    flux_boundary_specular!(fw, ff, f, u, ω, dt)
    flux_boundary_specular!(fw, fh, fb, h, b, u, ω, dt)

Specular reflection boundary flux

- @args: fluxes of conservative variables and distribution functions
- @args: distribution functions
- @args: velocity quadrature points and weights
- @args: time step

"""
function flux_boundary_specular!(
    fw::AbstractVector{T1},
    ff::AbstractVector{T2},
    f::AbstractVector{T3},
    u::T4,
    ω::T4,
    dt,
) where {
    T1<:Real,
    T2<:AbstractFloat,
    T3<:AbstractFloat,
    T4<:AbstractVector{<:AbstractFloat},
} # 1D1F1V

    fWall = similar(f)
    for i in eachindex(f)
        fWall[i] = f[end-i+1]
    end

    fw[1] = discrete_moments(fWall, u, ω, 1) * dt
    fw[2] = discrete_moments(fWall, u, ω, 2) * dt
    fw[3] = 0.5 * discrete_moments(fWall .* u .^ 2, u, ω, 1) * dt

    @. ff = u * fWall * dt

    return nothing

end

# ------------------------------------------------------------
# 1D2F1V
# ------------------------------------------------------------
function flux_boundary_specular!(
    fw::AbstractVector{T1},
    fh::T2,
    fb::T2,
    h::T3,
    b::T3,
    u::T4,
    ω::T4,
    dt,
) where {
    T1<:Real,
    T2<:AbstractVector{<:AbstractFloat},
    T3<:AbstractVector{<:AbstractFloat},
    T4<:AbstractVector{<:AbstractFloat},
}

    hWall = similar(h)
    bWall = similar(b)
    for i in eachindex(h)
        hWall[i] = h[end-i+1]
        bWall[i] = b[end-i+1]
    end

    fw[1] = discrete_moments(hWall, u, ω, 1) * dt
    fw[2] = discrete_moments(hWall, u, ω, 2) * dt
    fw[3] =
        (
            0.5 * discrete_moments(hWall .* u .^ 2, u, ω, 1) +
            0.5 * discrete_moments(bWall, u, ω, 1)
        ) * dt

    @. fh = u * hWall * dt
    @. fb = u * bWall * dt

    return nothing

end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

# contributed by Jarret Revels
# based on the Javascript program

using Printf

A(i,j) = 1.0 / ((i+j)*(i+j+1.0)/2.0+i+1.0)

function Au(u,w)
    n = length(u)
    for i = 1:n, j = 1:n
        j == 1 && (w[i] = 0)
        w[i] += A(i-1,j-1) * u[j]
    end
end

function Atu(w,v)
    n = length(w)
    for i = 1:n, j = 1:n
        j == 1 && (v[i] = 0)
        v[i] += A(j-1,i-1) * w[j]
    end
end

function perf_spectralnorm(n::Int=100)
    u = ones(Float64,n)
    v = zeros(Float64,n)
    w = zeros(Float64,n)
    vv = vBv = 0
    for i = 1:10
        Au(u,w)
        Atu(w,v)
        Au(v,w)
        Atu(w,u)
    end
    for i = 1:n
        vBv += u[i]*v[i]
        vv += v[i]*v[i]
    end
    return sqrt(vBv/vv)
end

n = parse(Int,ARGS[1])
@printf("%.9f\n", perf_spectralnorm(n))


text:
export UD_English
module UD_English
using ..MLDatasets

const defdir = joinpath(Pkg.dir("MLDatasets"), "datasets/ud_english")

traindata(dir=defdir) = readdata(dir, "en-ud-train.conllu")

devdata(dir=defdir) = readdata(dir, "en-ud-dev.conllu")

testdata(dir=defdir) = readdata(dir, "en-ud-test.conllu")

function readdata(dir, filename)
    mkpath(dir)
    url = "https://raw.githubusercontent.com/UniversalDependencies/UD_English/master"
    path = joinpath(dir, filename)
    isfile(path) || download("$(url)/$(filename)", path)
    CoNLL.read(path)
end

end


text:
module TestArrow

import Arrow
import MemoryAccesses
using ..Example: schedule_timings, two_arrays
using Test

function check_table(record1)
    io = IOBuffer()
    Arrow.write(io, [(sample = record1,)])
    seekstart(io)
    record2 = Arrow.Table(io).sample[1]

    table1 = MemoryAccesses.table(record1)
    table2 = MemoryAccesses.table(record2)
    @test table2 == table1
end

function test_table()
    @testset "$f" for f in [schedule_timings, two_arrays]
        check_table(f(100))
    end
end

end  # module


text:
# The AHL21 integrator WITHOUT derivatives.
"""

Carries out AHL21 mapping with compensated summation, WITHOUT derivatives
"""
function ahl21!(s::State{T},h::T) where T<:AbstractFloat
    h2 = 0.5*h;
    h6 = h/6
    kickfast!(s,h6)
    drift!(s,h2)
    drift_kepler!(s,h2)
    phic!(s,h)
    phisalpha!(s,h,T(2.0))
    kepler_drift!(s,h2)
    drift!(s,h2)
    kickfast!(s,h6)
    return
end

"""

Drifts all particles with compensated summation.
"""
function drift!(s::State{T},h::T) where {T <: Real}
    @inbounds for i=1:s.n, j=1:NDIM
        s.x[j,i],s.xerror[j,i] = comp_sum(s.x[j,i],s.xerror[j,i],h*s.v[j,i])
    end
    return
end

"""

Computes "fast" kicks for pairs of bodies in lieu of -drift+Kepler with compensated summation
"""
function kickfast!(s::State{T},h::T) where {T <: Real}
    s.rij .= zero(T)
    @inbounds for i=1:s.n-1
        for j = i+1:s.n
            if s.pair[i,j]
                for k=1:3
                    s.rij[k] = s.x[k,i] - s.x[k,j]
                    #r2 += s.rij[k]*s.rij[k]
                end
                r2inv = 1.0/dot_fast(s.rij)
                r3inv = r2inv*sqrt(r2inv)
                fac2 = h*GNEWT*r3inv
                for k=1:3
                    fac = fac2*s.rij[k]
                    s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i],-s.m[j]*fac)
                    s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j], s.m[i]*fac)
                end
            end
        end
    end
    return
end

"""

Computes correction for pairs which are kicked.
"""
function phic!(s::State{T},h::T) where {T <: Real}
    s.a .= zero(T)
    s.rij .= zero(T)
    s.aij .= zero(T)
    @inbounds for i=1:s.n-1, j = i+1:s.n
        if s.pair[i,j] # kick group
            for k=1:3
                s.rij[k] = s.x[k,i] - s.x[k,j]
                #r2 += s.rij[k]^2
            end
            r2inv = inv(dot_fast(s.rij))
            r3inv = r2inv*sqrt(r2inv)
            for k=1:3
                fac = GNEWT*s.rij[k]*r3inv
                facv = fac*2*h/3
                s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i],-s.m[j]*facv)
                s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j],s.m[i]*facv)
                s.a[k,i] -= s.m[j]*fac
                s.a[k,j] += s.m[i]*fac
            end
        end
    end
    coeff = h^3/36*GNEWT
    @inbounds for i=1:s.n-1 ,j=i+1:s.n
        if s.pair[i,j] # kick group
            for k=1:3
                s.aij[k] = s.a[k,i] - s.a[k,j]
                s.rij[k] = s.x[k,i] - s.x[k,j]
            end
            r2 = dot_fast(s.rij)
            r1 = sqrt(r2)
            ardot = dot_fast(s.aij,s.rij)
            fac1 = coeff/(r2*r2*r1)
            fac2 = 3*ardot
            for k=1:3
                fac = fac1*(s.rij[k]*fac2-r2*s.aij[k])
                s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i],s.m[j]*fac)
                s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j],-s.m[i]*fac)
            end
        end
    end
    return
end

"""

Computes the 4th-order correction with compensated summation.
"""
function phisalpha!(s::State{T},h::T,alpha::T) where {T <: Real}
    s.a .= zero(T)
    s.rij .= zero(T)
    s.aij .= zero(T)
    coeff = alpha*h^3/96*2*GNEWT

    fac = zero(T); fac1 = zero(T); fac2 = zero(T); r1 = zero(T); r2 = zero(T); r3 = zero(T)
    @inbounds for i=1:s.n-1
        for j = i+1:s.n
            if ~s.pair[i,j] # correction for Kepler pairs
                for k=1:3
                    s.rij[k] = s.x[k,i] - s.x[k,j]
                end
                r2 = dot_fast(s.rij)
                r3 = r2*sqrt(r2)
                fac2 = GNEWT/r3
                for k=1:3
                    fac = fac2*s.rij[k]
                    s.a[k,i] -= s.m[j]*fac
                    s.a[k,j] += s.m[i]*fac
                end
            end
        end
    end
    # Next, compute \tilde g_i acceleration vector (this is rewritten
    # slightly to avoid reference to \tilde a_i):
    @inbounds for i=1:s.n-1
        for j=i+1:s.n
            if ~s.pair[i,j] # correction for Kepler pairs
                for k=1:3
                    s.aij[k] = s.a[k,i] - s.a[k,j]
                    s.rij[k] = s.x[k,i] - s.x[k,j]
                end
                r2 = dot_fast(s.rij)
                r1 = sqrt(r2)
                ardot = dot_fast(s.aij,s.rij)
                # fac1 = coeff/r1^5
                fac1 = coeff/(r2*r2*r1)
                fac2 = (2*GNEWT*(s.m[i]+s.m[j])/r1 + 3*ardot)
                for k=1:3
                    fac = fac1*(s.rij[k]*fac2-r2*s.aij[k])
                    s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i], s.m[j]*fac)
                    s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j],-s.m[i]*fac)
                end
            end
        end
    end
    return
end

"""

Take a combined -drift+kepler step for each pair of bodies.
"""
function drift_kepler!(s::State{T}, h::T) where T<:Real
    n = s.n
    @inbounds for i in 1:n-1, j in i+1:n
        if ~s.pair[i,j]
            kepler_driftij_gamma!(s,i,j,h,true)
        end
    end
    return
end

"""

Take a combined kepler-drift step for each pair of bodies
"""
function kepler_drift!(s::State{T}, h::T) where T<:Real
    n = s.n
    @inbounds for i in n-1:-1:1, j in n:-1:i+1
        if ~s.pair[i,j]
            kepler_driftij_gamma!(s,i,j,h,false)
        end
    end
    return
end

"""

Carries out a Kepler step and reverse drift for bodies i & j with compensated summation.
"""
function kepler_driftij_gamma!(s::State{T},i::Int64,j::Int64,h::T,drift_first::Bool) where {T <: Real}
    @inbounds for k=1:NDIM
        s.x0[k] = s.x[k,i] - s.x[k,j] # x0 = positions of body i relative to j
        s.v0[k] = s.v[k,i] - s.v[k,j] # v0 = velocities of body i relative to j
    end
    gm = GNEWT*(s.m[i]+s.m[j])
    if gm == 0; return; end
    # Compute differences in x & v over time step:
    jac_delxv_gamma!(s,gm,h,drift_first,grad=false)
    mijinv =1.0/(s.m[i] + s.m[j])
    mi = s.m[i]*mijinv # Normalize the masses
    mj = s.m[j]*mijinv
    @inbounds for k=1:3
        # Add kepler-drift differences, weighted by masses, to start of step:
        s.x[k,i],s.xerror[k,i] = comp_sum(s.x[k,i],s.xerror[k,i], mj*s.delxv[k])
        s.x[k,j],s.xerror[k,j] = comp_sum(s.x[k,j],s.xerror[k,j],-mi*s.delxv[k])
        s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i], mj*s.delxv[3+k])
        s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j],-mi*s.delxv[3+k])
    end
    return
end


text:
module QueryOperators

using DataStructures
using IteratorInterfaceExtensions
using TableShowUtils
import DataValues

export Grouping, key

include("operators.jl")
include("NamedTupleUtilities.jl")

include("enumerable/enumerable.jl")
include("enumerable/enumerable_groupby.jl")
include("enumerable/enumerable_join.jl")
include("enumerable/enumerable_groupjoin.jl")
include("enumerable/enumerable_orderby.jl")
include("enumerable/enumerable_map.jl")
include("enumerable/enumerable_gather.jl")
include("enumerable/enumerable_filter.jl")
include("enumerable/enumerable_mapmany.jl")
include("enumerable/enumerable_defaultifempty.jl")
include("enumerable/enumerable_count.jl")
include("enumerable/enumerable_take.jl")
include("enumerable/enumerable_drop.jl")
include("enumerable/enumerable_unique.jl")
include("enumerable/show.jl")

include("source_iterable.jl")

end # module


text:
#######################################################################
# A simple formatting example using XlsxWriter.
# This program demonstrates the indentation cell format.
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("cell_indentation.xlsx")

	ws = add_worksheet!(wb)

	indent1 = add_format!(wb, Dict("indent"=>1))
	indent2 = add_format!(wb, Dict("indent"=>2))

	set_column!(ws, "A:A", 40)

	write!(ws, "A1", "This text is indented 1 level", fmt=indent1)
	write!(ws, "A2", "This text is indented 2 levels", fmt=indent2)

	close(wb)
	isfile("cell_indentation.xlsx")
end

test()


text:
module PlotReferenceImages
# A module is required for `pathof(PlotReferenceImages)` to work.

using Requires

# ------------------------------------------------------------------------------
# Plots docs and tests reference images

function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plotdocs.jl")
end

end # module


text:
import Random.rand!

"""
    _rescale!(mat)

Changes the matrix `mat` so that it is between `0` and `1`.
"""
function _rescale!(mat)
    mat .-= NaNMath.minimum(mat)
    mat ./= NaNMath.maximum(mat)
end

"""
    mask!(array::AbstractArray{<:AbstractFloat}, maskarray::AbstractArray{<:AbstractBool})

Modifies `array` so that the positions at which `maskarray` is `false` are
replaced by `NaN`.
"""
function mask!(array::AbstractArray{<:Float64}, maskarray::AbstractArray{<:Bool})
    (size(array) == size(maskarray)) || throw(DimensionMismatch("The dimensions of array, $(size(array)), and maskarray, $(size(maskarray)), must match. "))
    array[.!maskarray] .= NaN
    array
end

"""
    rand(alg, dims::Tuple{Vararg{Int64,2}}; mask=nothing) where {T <: Integer}

Creates a landscape of size `dims` (a tuple of two integers) following the model
defined by `alg`. The `mask` argument accepts a matrix of boolean values, and is
passed to `mask!` if it is not `nothing`.
"""
function Base.rand(alg::T, dims::Tuple{Int64,Int64}; mask=nothing) where {T <: NeutralLandscapeMaker}
    ret = Matrix{Float64}(undef, dims...)
    rand!(ret, alg; mask=mask)
end
Base.rand(alg::T, dims::Integer...; mask=nothing) where {T <: NeutralLandscapeMaker} = rand(alg, dims; mask = mask)

"""
    rand!(mat, alg) where {IT <: Integer}

Fill the matrix `mat` with a landscape created following the model defined by
`alg`. The `mask` argument accepts a matrix of boolean values, and is passed to
`mask!` if it is not `nothing`.
"""
function rand!(mat::AbstractArray{<:AbstractFloat,2} where N, alg::T; mask=nothing) where {T <: NeutralLandscapeMaker}
    _landscape!(mat, alg)
    isnothing(mask) || mask!(mat, mask)
    _rescale!(mat)
end


text:
abstract type lossFun end

export lossFun

### binaryCrossentropy

"""
    return the average cross entropy loss over vector of labels and predictions

    input:
        a := (?1, c,m) matrix of predicted values, where c is the number of classes
        y := (?1, c,m) matrix of predicted values, where c is the number of classes

        Note: in case the number of classes is one (1) it is okay to have
              a scaler values for a and y

    output:
        J := scaler value of the cross entropy loss
"""
abstract type binaryCrossentropy <: lossFun end

function binaryCrossentropy(a, y)

    aNew = prevnextfloat.(a)
    J = .-(y .* log.(aNew) .+ (1 .- y) .* log.(1 .- aNew))
    return J
end #binaryCrossentropy

export binaryCrossentropy

"""
    compute the drivative of cross-entropy loss function to the input of the
    layer dZ
"""
function dbinaryCrossentropy(a, y)
    dJ = a .- y
    return dJ
end #dbinaryCrossentropy

export dbinaryCrossentropy

### categoricalCrossentropy

abstract type categoricalCrossentropy <: lossFun end

function categoricalCrossentropy(a, y)
    aNew = prevnextfloat.(a)
    J = .-(y .* log.(aNew))
    return J
end

function dcategoricalCrossentropy(a, y)
    dJ = a .- y
    return dJ
end

export categoricalCrossentropy, dcategoricalCrossentropy

"""
    return previous float if x == 1 and nextfloat if x == 0
"""
prevnextfloat(x) = x==0 ? nextfloat(x) : x==1 ? prevfloat(x) : x

export prevnextfloat

### cost function

"""
    function cost(
        loss::Type{categoricalCrossentropy},
        A::AbstractArray{T1,N},
        Y::AbstractArray{T2,N},
    ) where {T1, T2, N}

Compute the cost for `categoricalCrossentropy` loss function

"""
function cost(
    loss::Type{categoricalCrossentropy},
    A::AbstractArray{T1,N},
    Y::AbstractArray{T2,N},
) where {T1, T2, N}

    c, m = size(A)[N-1:N]
    costs = sum(loss(Float64.(A),Y)) / m
    return Float64.(costs)

end #function cost

"""
    function cost(
        loss::Type{binaryCrossentropy},
        A::AbstractArray{T1,N},
        Y::AbstractArray{T2,N},
    ) where {T1, T2, N}

Compute the cost for `binaryCrossentropy` loss function

"""
function cost(
    loss::Type{binaryCrossentropy},
    A::AbstractArray{T1,N},
    Y::AbstractArray{T2,N},
) where {T1, T2, N}

    c, m = size(A)[N-1:N]
    costs = sum(loss(Float64.(A),Y)) / (c*m)
    return Float64.(costs)

end #function cost


text:
__precompile__(true)

module StringDistances

##############################################################################
##
## Export
##
##############################################################################
import Base: eltype, length, iterate, ==, hash, isless, convert, show
import Distances: evaluate, Hamming, hamming, PreMetric, SemiMetric
import IterTools: chain
export
evaluate,
compare,
Hamming,
Levenshtein,
DamerauLevenshtein,
Jaro,
QGram,
Cosine,
Jaccard,
SorensenDice,
Overlap,
RatcliffObershelp,
Winkler,
Partial,
TokenSort,
TokenSet,
TokenMax

##############################################################################
##
## include
##
##############################################################################
include("utils.jl")
include("distances/edit.jl")
include("distances/qgram.jl")
include("distances/RatcliffObershelp.jl")
include("compare.jl")

end


text:
@testset "collatz" begin

# A little puzzle inspired by the collatz conjecture, reach one by using the collatz conjecture moves (3n + 1 can be used even if the number is even)
neighbours(n) = n % 2 == 0 ? [n ÷ 2, 3n + 1] : [3n + 1]

# The best case is if the number is a power of 2, then the moves are just the log2, otherwise this is a very optimistic heuristic (it doesn't overestimate)
heuristic(n, goal) = floor(Int64, log2(n))

@testset "12" begin
  start = 12
  res = astar(neighbours, start, 1; heuristic)
  @test res.status == :success
  @test res.path == [12, 6, 3, 10, 5, 16, 8, 4, 2, 1]
  @test res.cost == 9
end

@testset "12 with maxcost" begin
  start = 12
  res = astar(neighbours, start, 1; heuristic, maxcost=3)
  @test res.status == :nopath
  @test res.path == [12, 6, 3]
  @test res.cost == 2
end

@testset "12 with cost function" begin
  start = 12

  # dividing costs 1, multiplying costs 4
  cost(state, neighbour) = state / neighbour == 2 ? 1 : 4

  res = astar(neighbours, start, 1; heuristic, cost)
  @test res.status == :success
  @test res.path == [12, 6, 3, 10, 5, 16, 8, 4, 2, 1]
  @test res.cost == 15
end

@testset "12 with maxdepth" begin
  start = 12
  res = astar(neighbours, start, 1; heuristic, maxdepth=3)
  @test res.status == :nopath
  @test res.path == [12, 6, 3]
  @test res.cost == 2
end

end


text:
#=
day11:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-11
=#

using AoC, Test, Memoize
using Random: shuffle

NEIGHBOURS = setdiff([CartesianIndex(i,j) for j in -1:1 for i in -1:1], [CartesianIndex(0,0)])
NOOP = (x...;y...)->nothing

function flash!(M, I, onUpdate=NOOP)
    onUpdate(M)
    for J in shuffle(NEIGHBOURS)
        if checkbounds(Bool, M, I+J) && M[I+J] <= 9
            M[I+J] += 1
            if M[I+J] > 9
                flash!(M, I+J, onUpdate)
            end
        end
    end
    return M
end

@test flash!([[1,2] [3,10]], CartesianIndex(2,2)) == [[2,3] [4,10]]
@test flash!([[1,2] [8,10]], CartesianIndex(2,2)) == [[2,3] [9,10]]
@test flash!([[1,2] [9,10]], CartesianIndex(2,2)) == [[3,4] [10,10]]

function round!(M, onUpdate=NOOP)
    for I in shuffle(CartesianIndices(M))
        if M[I] <= 9
            M[I] += 1
            onUpdate(M)
            if M[I] > 9
                flash!(M, I, onUpdate)
            end
        end
    end
    onUpdate(M; key=true)

    flashes = 0
    for I in eachindex(M)
        if M[I] > 9
            flashes += 1
            M[I] = 0
        end
    end

    return flashes
end

test1 = [[1,1,1,1,1] [1,9,9,9,1] [1,9,1,9,1] [1,9,9,9,1] [1,1,1,1,1]]
@test round!(test1) == 9
@test test1 == [[3,4,5,4,3] [4,0,0,0,4] [5,0,0,0,5] [4,0,0,0,4] [3,4,5,4,3]]

function rounds!(M, n, onUpdate=NOOP)
    flashes = 0
    for i in 1:n
        flashes += round!(M, onUpdate)
    end
    return flashes
end
@test rounds!(exampleIntMap(11,1), 1) == 0
@test rounds!(exampleIntMap(11,1), 2) == 35
@test rounds!(exampleIntMap(11,1), 3) == 35+45
@test rounds!(exampleIntMap(11,1), 10) == 204

function findBrightFlash(M, onUpdate=NOOP)
    round = 0
    while true
        round += 1
        flashes = round!(M, onUpdate)
        if flashes == length(M)
            return round
        end
    end
end

part1(M, onUpdate=NOOP) = rounds!(M, 100, onUpdate)
@test part1(exampleIntMap(11,1)) == 1656

part2(M,onUpdate=NOOP) = findBrightFlash(M, onUpdate)
@test part2(exampleIntMap(11,1)) == 195

intMap(11) |> m -> @time part1(m) |> show
intMap(11) |> m -> @time part2(m) |> show

using FileIO, ColorTypes, FixedPointNumbers

@memoize circle(r) = filter(I->I[1]^2 + I[2]^2<=r^2, [(i,j) for j in -r:r for i in -r:r]) .|> CartesianIndex

function writeImage(f, M, magnification, radius, flashRadius=radius)

    animation = zeros(RGB{N0f8}, (size(M) .* magnification .+ 1)..., 1000)
    frames = 0
    updates = 0

    octopus = circle(radius)
    flash = circle(flashRadius)

    function mapUpdate(M; key=false)
        updates += 1
        if !key && (updates % 20) != 0
            return
        end
        frames += 1
        if size(animation, 3) < frames
            newAnim = similar(animation, size(animation,1),size(animation,2),size(animation,3)*2)
            copyto!(newAnim, CartesianIndices(animation), animation, CartesianIndices(animation))
            animation = newAnim
        end
        frame = view(animation, :,:,frames)
        if flashRadius != radius
            for I in CartesianIndices(M)
                J = CartesianIndex(Tuple(I) .* magnification .+ 1 .- (magnification ÷ 2))
                if M[I] > 9
                    for K in flash
                        if checkbounds(Bool, frame, J+K)
                            frame[J+K] = 0.3
                        end
                    end
                end
            end
        end
        for I in CartesianIndices(M)
            J = CartesianIndex(Tuple(I) .* magnification .+ 1 .- (magnification ÷ 2))
            if M[I] <= 9
                for K in octopus
                    frame[J+K] = M[I]/20
                end
            else
                for K in octopus
                    frame[J+K] = 1.0
                end
            end
        end
    end

    f(M, mapUpdate)

    save("target/day11.gif", view(animation,:,:,1:frames))
end

writeImage((M,onUpdate)->rounds!(M,300,onUpdate), intMap(11), 20, 5, 15)