text:
#= ######################################################
inside functions
###################################################### =#

function get_vectors_for_tilewalk(p1::Point2D{T}, p2::Point2D{T}) where T<:AbstractFloat
    # calculate necessary data for tilewalk
    vec = p2 - p1
    vec_l = norm(vec)
    vec_norm = normit(vec)
    tol = 1E-8
    if vec_norm.x > tol
        dirx = 1
    elseif vec_norm.x < ((-1) * tol)
        dirx = -1
    else
        dirx = 0
    end
    if vec_norm.y > tol
        diry = 1
    elseif vec_norm.y < ((-1) * tol)
        diry = -1
    else
        diry = 0
    end
    return Index2D(dirx, diry), vec_l, vec_norm
end

function get_start_tile(p1::Point2D{T1}, dir::Index2D{T2}, n::T2, dx::T1, dy::T1) where {T1<:AbstractFloat, T2<:Integer}
    # get tile of starting point
    tx_step = 1
    ty_step = 1
    tol = 1E-8
    if dir.x > 0
        tx_step += floor(Integer,(p1.x + dir.x * tol) / dx)
    elseif dir.x < 0
        tx_step += floor(Integer,(p1.x + dir.x * tol) / dx)
    else
        # wenn exakt auf einer tile Linie -> 1D walk
        # entscheidet Vorzeichen vor tol wo der Bucket liegt
        # wichtig wenn linie am Rand
        if p1.x > (dx * n - tol)
            tx_step += floor(Integer,(p1.x - tol) / dx)
        else
            tx_step += floor(Integer,(p1.x + tol) / dx)
        end
    end
    if dir.y > 0
        ty_step += floor(Integer,(p1.y + dir.y * tol) / dy)
    elseif dir.y < 0
        ty_step += floor(Integer,(p1.y + dir.y * tol) / dy)
    else
        # wenn exakt auf einer tile Linie -> 1D walk
        # entscheidet Vorzeichen vor tol wo der Bucket liegt
        # wichtig wenn linie am Rand
        if p1.y > (dy * n - tol)
            ty_step += floor(Integer,(p1.y - tol) / dy)
        else
            ty_step += floor(Integer,(p1.y + tol) / dy)
        end
    end
    return Index2D(tx_step, ty_step)
end

function Sx(m)
    val = sqrt(1 + m^2)
    return val
end

function Sy(m)
    val = sqrt(1 + (m^(-1))^2)
    return val
end

function get_lengths_to_next_tiles(p1::Point2D{T1}, p2::Point2D{T1}, dir::Index2D{T2}, tile::Index2D{T2}, dx::T1, dy::T1)::Tuple{T1, T1} where {T1<:AbstractFloat, T2<:Integer}
    # get lengths in x and y direction to next tiles
    if dir.x > 0
        dx_step = dx * (tile.x)
    elseif dir.x < 0
        dx_step = dx * (tile.x - 1)
    else
        dx_step = NaN
    end
    if dir.y > 0
        dy_step = dy * (tile.y)
    elseif dir.y < 0
        dy_step = dy * (tile.y - 1)
    else
        dy_step = NaN
    end
    # println("dx: ", dx_step, "  dy: ", dy_step)
    lx = dir.x * (dx_step - p1.x) * Sx((p2.y-p1.y) / (p2.x-p1.x))
    ly = dir.y * (dy_step - p1.y) * Sy((p2.y-p1.y) / (p2.x-p1.x))
    # correction for 1D case and changing NaN into Inf
    if isnan(lx)
        lx = Inf
    end
    if isnan(ly)
        ly = Inf
    end
    # println("lx: ",lx,"  ly: ",ly)
    return lx, ly
end

function get_next_tile(tile::Index2D{T2}, dir::Index2D{T2}, lx::T1, ly::T1)::Tuple{Index2D{T2},T1,Symbol} where {T1<:AbstractFloat, T2<:Integer}
    # get next tile
    # bei einem 2D Fall exakt durch die Knoten entscheidet
    # hier < und <= welche Seite das next tile liegt
    # default von mir: < geht immer in nächste y richtung
    if abs(lx) < abs(ly)
        # println("moving in x direction")
        tile_new = tile + Index2D(dir.x, 0)
        lcurrent = lx
        color = :cyan
    else
        # println("moving in y direction")
        tile_new = tile + Index2D(0, dir.y)
        lcurrent = ly
        color = :coral
    end
    return tile_new, lcurrent, color
end

function create_randomly_occupied_tiles(n::T; density = 0.1) where T<:Integer
    # create random occ tiles
    # only for testing
    t_occ = rand(n,n)
    t_occ[t_occ .> (1-density)] .= 1
    t_occ[t_occ .<= (1-density)] .= 0
    t_occ = convert.(Integer,t_occ)
    t_occ_n = Matrix{Union{Vector{T},Missing}}(missing,n,n)
    # t_occ_n[t_occ .== 1] .= Vector{T}(undef,2)
    # t_occ = rand([0, 1], n, n)
    # t_occ = zeros(Integer, n, n)
    for i = 1:n, j = 1:n
        if t_occ[i,j] == 1
            t_occ_n[i,j] = [1,1]
        end
    end
    return t_occ_n
end

function get_max_steps(n::T) where T<:Integer
    m1 = round(Int64,sqrt(n*n+n*n)*2)
    # m2 = convert(typeof(n),m1)
    return m1
end

#= ######################################################
final functions
###################################################### =#

function tilewalk_with_check_for_occ_tiles(fig, ax, p1::Point2D{T1}, p2::Point2D{T1}, dx::T1, dy::T1, n::T2, t_occ::Matrix{Union{Vector{T2},Missing}}) where {T1<:AbstractFloat, T2<:Integer}
    # do tilewalk between two points and check for occupied tiles
    dir, ltot, lvec = get_vectors_for_tilewalk(p1, p2)
    tile = get_start_tile(p1, dir, n, dx, dy)
    println("start tile number: ", tile.x, ", ", tile.y)
    max_steps = get_max_steps(n)
    if !ismissing(t_occ[tile.x, tile.y])
        println("------> occupied ")
        start_walk = false
    else
        # print start tile in color
        poly!(ax, Rect((tile.x-1) * dx, (tile.y-1) * dy, dx, dy), color = (:blue, 0.3))
        start_walk = true
    end
    if start_walk
        # do walk
        for step = 1:max_steps
            lx, ly = get_lengths_to_next_tiles(p1, p2, dir, tile, dx, dy)
            # check if end reached
            tol2 = 1E-8 # necessary because of not exact numbers
            if lx >= (ltot - tol2) && ly >= (ltot - tol2)
                break
            end
            tile, lcurrent, pcolor = get_next_tile(tile, dir, lx, ly)
            # println("lcurrent: ",lcurrent)
            # plot next point
            p = p1 + lvec * lcurrent
            scatter!(ax, Point2f0(p.x, p.y), color = pcolor, markersize = 15)
            println("next tile number: ", tile.x, ", ", tile.y)
            # check this tile
            if !ismissing(t_occ[tile.x, tile.y])
                println("------> occupied ")
                # plot hitting point
                p = p1 + lvec * lcurrent
                scatter!(ax, Point2f0(p.x, p.y), color = :blue, markersize = 15)
                break
            end
            # print tile in color
            poly!(ax, Rect((tile.x-1) * dx, (tile.y-1) * dy, dx, dy), color = (:blue, 0.3))
        end
    end
end

function tilewalk_with_return(fig, ax, p1::Point2D{T1}, p2::Point2D{T1}, dx::T1, dy::T1, n::T2) where {T1<:AbstractFloat, T2<:Integer}
    # do tilewalk between two points and check for occupied tiles
    max_steps = get_max_steps(n)
    tile_list = Vector{Index2D{T2}}(undef,max_steps)
    dir, ltot, lvec = get_vectors_for_tilewalk(p1, p2)
    tile = get_start_tile(p1, dir, n, dx, dy)
    println("start tile number: ", tile.x, ", ", tile.y)
    hit = 1
    tile_list[hit] = Index2D(tile.x, tile.y)
    # print start tile in color
    poly!(ax, Rect((tile.x-1) * dx, (tile.y-1) * dy, dx, dy), color = (:blue, 0.3))
    # do walk
    for step = 1:max_steps
        lx, ly = get_lengths_to_next_tiles(p1, p2, dir, tile, dx, dy)
        # check if end reached
        tol2 = 1E-8 # necessary because of not exact numbers
        if lx >= (ltot - tol2) && ly >= (ltot - tol2)
            break
        end
        tile, lcurrent, pcolor = get_next_tile(tile, dir, lx, ly)
        println("next tile number: ", tile.x, ", ", tile.y)
        hit += 1
        tile_list[hit] = Index2D(tile.x, tile.y)
        # print tile in color
        poly!(ax, Rect((tile.x-1) * dx, (tile.y-1) * dy, dx, dy), color = (:blue, 0.3))
    end
    return tile_list[1:hit]
end

function tilewalk_with_return!(tile_list::Vector{Index2D{T2}}, p1::Point2D{T1}, p2::Point2D{T1}, dx::T1, dy::T1, n::T2)::T2 where {T1<:AbstractFloat, T2<:Integer}
    # do tilewalk between two points and check for occupied tiles
    max_steps = get_max_steps(n)
    dir, ltot, lvec = get_vectors_for_tilewalk(p1, p2)
    tile = get_start_tile(p1, dir, n, dx, dy)
    hit = 1
    tile_list[hit] = tile
    for step = 1:max_steps # do walk
        if step == max_steps
            @warn "reached max allowed tiles in tilewalk"
        end
        lx, ly = get_lengths_to_next_tiles(p1, p2, dir, tile, dx, dy)
        # check if end reached
        tol2 = 1E-8 # necessary because of not exact numbers
        if lx >= (ltot - tol2) && ly >= (ltot - tol2)
            break
        end
        tile, nothing, nothing = get_next_tile(tile, dir, lx, ly)
        hit += 1
        tile_list[hit] = tile
    end
    return hit
end


text:
# CSI
# ===
#
# Specification:
#   1: https://github.com/samtools/hts-specs/blob/master/CSIv1.pdf
#   2: https://github.com/samtools/hts-specs/blob/master/CSIv2.pdf
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

struct CSIBin
    # distinct bin
    bin::UInt32

    # virtual file offset
    loffset::UInt64

    # number of records stored in the bin
    n_rec::UInt64

    # chunks
    chunks::Vector{Chunk}
end

struct CSI
    # CSI version (0x01 or 0x02)
    version::UInt8

    # number of bits for the minimal interval
    minshift::Int32

    # depth of the binning index
    depth::Int32

    # auxilary data
    aux::Vector{UInt8}

    # bins
    bins::Vector{Vector{CSIBin}}

    # number of unmapped reads
    n_no_coor::Nullable{Int}
end

function Base.show(io::IO, csi::CSI)
    print(io, summary(csi), "(<version=$(csi.version),#seqs=$(length(csi.bins))>)")
end

"""
    CSI(input::AbstractString)

Load a CSI index from `input`.
"""
function CSI(input::AbstractString)
    return open(x->CSI(x), input)
end

"""
    CSI(input::IO)

Load a CSI index from `input`.
"""
function CSI(input_::IO)
    input = CodecBGZF.BGZFDecompressorStream(input_)

    # check magic bytes
    C = read(input, UInt8)
    S = read(input, UInt8)
    I = read(input, UInt8)
    v = read(input, UInt8)
    if C != UInt32('C') || S != UInt32('S') || I != UInt8('I') || v ∉ (0x01, 0x02)
        error("invalid csi magic bytes")
    end

    minshift = read(input, Int32)
    depth = read(input, Int32)
    l_aux = read(input, Int32)
    aux = read(input, l_aux)
    n_ref = read(input, Int32)
    bins = Vector{CSIBin}[]
    for _ in 1:n_ref
        n_bin = read(input, Int32)
        refbins = CSIBin[]
        for _ in 1:n_bin
            bin = read(input, UInt32)
            loffset = read(input, UInt64)
            if v == 0x02
                n_rec = read(input, UInt64)
            else
                n_rec = typemax(UInt64)
            end
            n_chunk = read(input, Int32)
            chunks = Chunk[]
            for _ in 1:n_chunk
                chunk_beg = read(input, UInt64)
                chunk_end = read(input, UInt64)
                push!(chunks, (chunk_beg, chunk_end))
            end
            push!(refbins, CSIBin(bin, loffset, n_rec, chunks))
        end
        push!(bins, refbins)
    end

    if !eof(input)
        n_no_coor = Nullable{Int}(read(input, UInt64))
    else
        n_no_coor = Nullable{Int}()
    end

    return CSI(v, minshift, depth, aux, bins, n_no_coor)
end

function overlapchunks(csi::CSI, seqid::Integer, interval::UnitRange{<:Integer})
    if !(1 ≤ seqid ≤ endof(csi.bins))
        throw(ArgumentError("sequence id $(seqid) is out of range"))
    elseif isempty(interval)
        return Chunk[]
    end
    chunks = Chunk[]
    bins = reg2bins(interval, csi.minshift, csi.depth)
    binindex = csi.bins[seqid]
    for bin in bins
        for x in binindex
            if x.bin == bin
                for chunk in x.chunks
                    if chunk[2] > x.loffset
                        push!(chunks, chunk)
                    end
                end
            end
        end
    end
    sort!(chunks)
    reduce!(chunks)
    return chunks
end

# Calculate bins overlapping an interval (inclusive, one-based).
reg2bins(interval::UnitRange{<:Integer}, minshift::Integer, depth::Integer) =
    reg2bins(convert(UnitRange{Int}, interval),
             convert(Int, minshift),
             convert(Int, depth))
function reg2bins(interval::UnitRange{Int}, minshift::Int, depth::Int)
    bins = Int[]
    s = minshift + depth * 3
    bin_start = 0
    for d in 0:depth
        for k in ((first(interval)-1)>>s):((last(interval)-1)>>s)
            push!(bins, bin_start + k)
        end
        s -= 3
        bin_start = 8 * bin_start + 1
    end
    return bins
end

ext: jl


text:
axis_spec =
[
    (:_type, AbstractString, "x"),
    (:scale, AbstractString, "x"),
    (:orient, AbstractString, nothing),
    (:title, AbstractString, nothing),
    (:titleOffset, Number, nothing),
    (:format, AbstractString, nothing),
    (:formatType, AbstractString, nothing),
    (:ticks, Number, nothing),
    (:values, Vector, nothing),
    (:subdivide, Number, nothing),
    (:tickPadding, Number, nothing),
    (:tickSize, Number, nothing),
    (:tickSizeMajor, Number, nothing),
    (:tickSizeMinor, Number, nothing),
    (:tickSizeEnd, Number, nothing),
    (:offset, Union{Number, Dict{Any, Any}}, nothing),
    (:layer, AbstractString, "front"),
    (:grid, Bool, false),
    (:properties, Union{Dict{Any,Any}, VegaMarkPropertySet}, nothing)
]

primitivefactory(:VegaAxis, axis_spec)

max_stars_repo_name: jpthor/vega


text:
# functions that classify regression results
isFERegressionResult(r::FixedEffectModel) = has_fe(r)
isIVRegressionResult(r::FixedEffectModel) = has_iv(r)
isOLSRegressionResult(r::FixedEffectModel) = !has_iv(r)

# FE and IV regression not supported in GLM.jl
isFERegressionResult(r::TableRegressionModel) = false
isIVRegressionResult(r::TableRegressionModel) = false
isOLSRegressionResult(r::TableRegressionModel) = isa(r.model, LinearModel)

# get display labels for terms
function name(t::InteractionTerm)
    n = name(t.terms[1])
    for term in t.terms[2:end]
        n = "$(n) * $(name(term))"
    end
    n
end
function name(t::FunctionTerm)
    string(t.args_parsed[1].sym)
end


text:
using Pkg
Pkg.activate("./")
Pkg.instantiate()
using Mandelbrot
using Plots
using BenchmarkTools

#%% code to make some of the example images
# starting point
cmap = Mandelbrot.cycle_cmap(:inferno, 5)
xmin = -2.2
xmax = 0.8
ymin = -1.2
ymax = 1.2

fractal0_data = FractalData(
    xmin,
    xmax,
    ymin,
    ymax,
    width = Mandelbrot.w_4k,
    height = Mandelbrot.h_4k,
    colormap = cmap,
    maxIter = 1500,
    scale_function = x -> x,
)

cmap1 = Mandelbrot.cycle_cmap(:inferno, 5)
xmin1 = -1.744453831814658538530
xmax1 = -1.744449945239591698236
ymin1 = 0.022017835126305555133
ymax1 = 0.022020017997233506531

fractal1_data = FractalData(
    xmin1,
    xmax1,
    ymin1,
    ymax1,
    width = Mandelbrot.w_4k,
    height = Mandelbrot.h_4k,
    colormap = cmap1,
    maxIter = 1500,
    scale_function = x -> x,
)

cmap2 = cgrad(:ice)
scale2 = log
maxIter2 = 5000
xmin2 = 0.308405876790033128474
xmax2 = 0.308405910247503605302
ymin2 = 0.025554220954294027410
ymax2 = 0.025554245987221578418

fractal2_data = FractalData(
    xmin2,
    xmax2,
    ymin2,
    ymax2,
    width = Mandelbrot.w_4k,
    height = Mandelbrot.h_4k,
    colormap = cmap2,
    maxIter = maxIter2,
    scale_function = scale2,
)

# cmap3b = Mandelbrot.pumpkin(4)
cmap3c = Mandelbrot.cycle_cmap(:inferno, 50)
scale3 = x -> x^-5
maxIter3 = 5000
xmin3 = 0.307567454839614329536
xmax3 = 0.307567454903142214608
ymin3 = 0.023304267108419154581
ymax3 = 0.023304267156089095573

fractal3_data = FractalData(
    xmin3,
    xmax3,
    ymin3,
    ymax3,
    width = Mandelbrot.w_4k,
    height = Mandelbrot.h_4k,
    colormap = cmap3c,
    maxIter = maxIter3,
    scale_function = scale3,
)

fractal3b_data = FractalData(
    xmin3,
    xmax3,
    ymin3,
    ymax3,
    width = 2 * Mandelbrot.w_4k,
    height = 2 * Mandelbrot.h_4k,
    colormap = cmap3c,
    maxIter = maxIter3,
    scale_function = scale3,
)

cmap4 = Mandelbrot.alien_space(50)
scale4 = x -> x^-5
maxIter4 = 50000
xmin4 = 0.2503006273651145643691
xmax4 = 0.2503006273651201870891
ymin4 = 0.0000077612880963380370
ymax4 = 0.0000077612881005550770

xmin4b = BigFloat("0.2503006273651145643691")
xmax4b = BigFloat("0.2503006273651201870891")
ymin4b = BigFloat("0.0000077612880963380370")
ymax4b = BigFloat("0.0000077612881005550770")

fractal4b_data = FractalData(
    xmin4b,
    xmax4b,
    ymin4b,
    ymax4b,
    width = Mandelbrot.w_HD,
    height = Mandelbrot.h_HD,
    colormap = cmap4,
    maxIter = maxIter4,
    scale_function = scale4,
)

#%% number 1
computeMandelbrot!(fractal1_data)

display_fractal(fractal1_data,
    # filename = "images/mandelbrot1.png"
)
#%%
display_fractal(
    fractal1_data,
    scale = log,
    # filename = "mandelbrot-fractal/images/mandelbrot1b.png"
)

#%% let's change the colormap and the scale function
fractal1_data.colormap = Mandelbrot.fire_and_ice(2)
fractal1_data.background_color = :black
scale = x-> x^-1

display_fractal(fractal1_data, scale = scale,
 # filename="images/mandelbrot1c",
 offset=0)

#%% number 2
computeMandelbrot!(fractal2_data)

display_fractal(fractal2_data,
    # filename = "mandelbrot-fractal/images/mandelbrot2.png"
)

#%% number 3
computeMandelbrot!(fractal3_data)

display_fractal(fractal3_data,
    # filename = "mandelbrot-fractal/images/mandelbrot3d.png"
)
#%% number 3 using GPU
@btime computeMandelbrot!(fractal3_data, use_GPU=true, verbose=false)
@btime computeMandelbrot!(fractal3_data, use_GPU=false, verbose=false)

display_fractal(fractal3_data,
    # filename = "mandelbrot-fractal/images/mandelbrot3d.png"
)

#%% number 4
computeMandelbrot!(fractal4_data)

display_fractal(fractal4_data,
    # filename = "mandelbrot-fractal/images/mandelbrot4.png"
)
#%% navigation

fractal0_data.maxIter = 50
preview_fractal(fractal0_data, scale = :linear)
Mandelbrot.move_center!(fractal0_data, -41, 0)

fractal0_data.maxIter = 500 #increas maximum number of iterations
Mandelbrot.zoom!(fractal0_data, 100)
Mandelbrot.move_center!(fractal0_data, -35, 0)
Mandelbrot.zoom!(fractal0_data, 10)
Mandelbrot.move_center!(fractal0_data, -30, 0)
fractal0_data.maxIter = 1000
preview_fractal(fractal0_data, scale = x->1/log10(x)) #nice fractal!

coords = Mandelbrot.get_coords(fractal0_data) #let's save the coordinates for future use

Mandelbrot.set_coords(fractal0_data, coords...)

computeMandelbrot!(fractal0_data) # compute 4k resolution image
display_fractal(fractal0_data, scale = x->1/log10(x),
    filename = "images/mandelbrot_movement.png"
) #plot and save the fractal

#%% compute an animation up to point 3

Mandelbrot.create_animation(
    (xmin3, xmax3, ymin3, ymax3),
    n_frames = 500,
    scale = log10,
    colormap = Mandelbrot.fire_and_ice(),
)


text:
using Debugger: Debugger, @enter, execute_command, RunDebugger, @make_frame
import JuliaInterpreter: JuliaInterpreter, Frame, @lookup, Compiled, pc_expr

using Test
using .Meta: isexpr
import REPL

if !isdefined(Main, :isnothing)
    isnothing(x) = x === nothing
end

Debugger.set_highlight(Debugger.HIGHLIGHT_OFF)

#@testset "Main tests" begin
    include("utils.jl")
    include("misc.jl")
    include("stepping.jl")
    include("interpret.jl")
    include("evaling.jl")
    include("ui.jl")
#end


text:
@params struct AugLag2 <: AbstractOptimizer
    primaloptimizer
    dualoptimizer
end
function AugLag2(;
    primaloptimizer = Optim.ConjugateGradient(linesearch=Optim.LineSearches.BackTracking(iterations = 10)),
    dualoptimizer = Optim.GradientDescent(linesearch=Optim.LineSearches.BackTracking(iterations = 10)),
)
    return AugLag2(primaloptimizer, dualoptimizer)
end

@params struct AugLag2Options
    primaloptions
    dualoptions
    maxiter
    tol
    quadfactor
end
function AugLag2Options(alg::AugLag2;
    primaloptions = alg.primaloptimizer isa MMA02 || alg.primaloptimizer isa MMA87 ? MMAOptions(maxiter = 100, tol = Tolerance(kkt = 1e-4)) : Optim.Options(outer_iterations = 10, iterations = 10),
    dualoptions = Optim.Options(outer_iterations = 10, iterations = 10),
    maxiter = 10,
    tol = Tolerance(),
    quadfactor = 10,
)
    return AugLag2Options(primaloptions, dualoptions, maxiter, tol, quadfactor)
end

function Solution(lagmodel::AugLag2Model)
    prevx = copy(getmin(lagmodel))
    x = copy(prevx)
    prevf = Inf
    f = Inf
    λ = copy(getlinweights(lagmodel))
    g = copy(λ)
    convstate = ConvergenceState()
    return Solution(prevx, x, getlinweights(lagmodel), prevf, f, nothing, g, nothing, convstate)
end

@params mutable struct AugLag2Workspace <: Workspace
    model::VecModel
    lagmodel::AugLag2Model
    x0::AbstractVector
    optimizer::AugLag2
    options::AugLag2Options

    solution::Solution
    convcriteria::ConvergenceCriteria
    callback::Function

    trace::Trace
    outer_iter::Int
    iter::Int
	fcalls::Int
end
function AugLag2Workspace(
    model::VecModel,
    optimizer::AugLag2,
    x0::AbstractVector;
    options::AugLag2Options = AugLag2Options(optimizer),
    convcriteria::ConvergenceCriteria = KKTCriteria(),
    plot_trace::Bool = false,
    show_plot::Bool = plot_trace,
    save_plot = nothing,
    callback::Function = plot_trace ? LazyPlottingCallback(; show_plot = show_plot, save_plot = save_plot) : NoCallback(),
    kwargs...,
)
    T = eltype(x0)
    lagmodel = AugLag2Model(model; kwargs...)

    # Convergence
    solution = Solution(lagmodel)
    #assess_convergence!(solution, model, options.tol, convcriteria)

    # Trace
    trace = Trace([])

    # Iteraton counter
    fcalls, outer_iter, iter = 1, 0, 0

    return AugLag2Workspace(
        model,
        lagmodel,
        x0,
        optimizer,
        options,

        solution,
        convcriteria,
        callback,

        trace,
        outer_iter,
        iter,
        fcalls,
    )
end

Workspace(model::VecModel, alg::AugLag2, x0::AbstractVector; kwargs...) = AugLag2Workspace(model, alg, x0; kwargs...)

function optimize!(workspace::AugLag2Workspace)
    @unpack lagmodel, solution, options, convcriteria = workspace
    @unpack callback, optimizer, trace = workspace
    @unpack x0, outer_iter, iter, fcalls = workspace
    @unpack primaloptimizer, dualoptimizer = optimizer
    @unpack maxiter, primaloptions, dualoptions, quadfactor = options
    @unpack prevx, x, λ = solution

    xl, xu = getmin(lagmodel), getmax(lagmodel)

    T = eltype(x)
    ni, nj = length(λ), length(x)

    # Original model
    model = lagmodel |> getparent

    auglag = getobjective(lagmodel)

    cb = (tr; kwargs...) -> begin
        solution = deepcopy(solution)
        solution.prevx .= solution.x
        solution.x .= getx(lagmodel)
        solution.λ .= getλ(lagmodel)
        solution.prevf = solution.f
        solution.f = getorigobjval(lagmodel)
        solution.g .= getorigconstrval(lagmodel)
        assess_convergence!(solution, lagmodel, options.tol, convcriteria)
        callback(solution)
        return hasconverged(solution)
	end

    primaloptimizerfunction(λ) = begin
        setlinweights!(auglag, λ)
        # Primal objective to be minimized
        # Calculates the objective and its gradient
        # Optim-compatible objective
        if primaloptimizer isa MMA02 || primaloptimizer isa MMA87
            primalmodel = Model()
            addvar!(primalmodel, xl, xu)
            set_objective!(primalmodel, getprimalobjective(auglag))
            result = optimize(primalmodel, primaloptimizer, clamp.(getx(lagmodel), xl .+ 1e-3, xu .- 1e-3), options = primaloptions, callback = cb, convcriteria = KKTCriteria())
            fcalls += result.fcalls
        else
            primalobj = getoptimobj(getprimalobjective(auglag), true)
            result = Optim.optimize(
                Optim.only_fg!(primalobj),
                xl,
                xu,
                clamp.(getx(lagmodel), xl .+ 1e-3, xu .- 1e-3),
                Optim.Fminbox(primaloptimizer),
                @set(primaloptions.callback = cb),
            )
            fcalls += result.f_calls
        end
        if getx(lagmodel) != result.minimizer
            primalobj(result.minimizer)
            setx!(lagmodel, result.minimizer)
        end
        if debugging[]
            @show getx(lagmodel)
        end
        return getx(lagmodel), getorigobjval(lagmodel), getorigconstrval(lagmodel)
    end

    # Dual objective to be minimized - original objective will be maximized
    # Calculates negative the objective and its gradient
    # Optim-compatible objective
    dualobj = getoptimobj(getdualobjective(auglag, primaloptimizerfunction), false)

    # Lower and upper bounds on the dual variables
    λl = zeros(ni) .+ 1e-10
    λu = fill(Inf, ni)

    # Solve the dual problem by minimizing negative the dual objective value
    for i in 1:maxiter
        setquadweight!(lagmodel, min(getquadweight(lagmodel) * quadfactor, 1e10))
        if debugging[]
            @show getquadweight(lagmodel)
        end
        λresult = Optim.optimize(
            Optim.only_fg!(dualobj),
            λl,
            λu,
            max.(copy(λ), 1e-3),
            Optim.Fminbox(dualoptimizer),
            dualoptions, #@set(dualoptions.callback = cb),
        )
        λ .= λresult.minimizer
        if debugging[]
            @show λ
            @show getx(lagmodel)
        end
        if λ != getλ(lagmodel)
            dualobj(λ)
            setλ!(lagmodel, λ)
        end
    end
    if debugging[]
        #@show getx(lagmodel)
    end

    @pack! workspace = iter, fcalls
    solution.x .= getx(lagmodel)
    solution.λ .= getλ(lagmodel)
    solution.f = getorigobjval(lagmodel)
    solution.g .= getorigconstrval(lagmodel)
    callback(solution, update = true)

    results = GenericResult(
        optimizer,
        x0,
        solution.x,
        solution.f,
        iter,
        iter == options.maxiter,
        options.tol,
        solution.convstate,
        fcalls,
    )
    return results
end

function getoptimobj(obj, minimize = true)
    optimobj(z) = optimobj(1.0, nothing, z)
    function optimobj(F, G, z)
        if G !== nothing
            val, grad = value_gradient(obj, z)
            if minimize
                G[:] .= grad
            else
                G[:] .= .-grad
            end
            if F !== nothing
                if minimize
                    return val
                else
                    return -val
                end
            end
        end
        # No gradient necessary, just return the log joint.
        if F !== nothing
            if minimize
                return obj(z)
            else
                return -obj(z)
            end
        end
        return nothing
    end
    return optimobj
end


text:
using ParserCombinator: Parsers.GML

function _gml_read_one_graph(gs)
    nodes = [x[:id] for x in gs[:node]]
    g = Graph()
    sds = [(Int(x[:source]), Int(x[:target])) for x in gs[:edge]]
    for (s,d) in (sds)
        addedge!(g, s, d)
    end
    return g
end

function loadgml(gname::AbstractString)
    p = GML.parse_dict(readall(gname))
    for gs in p[:graph]
		return _gml_read_one_graph(gs)
    end
    error("Graph $gname not found")
end

function savegml(gname::AbstractString, g::Graph, c)
    io = open(gname, "w")
    println(io, "graph")
    println(io, "[")
    println(io, "directed 0")
    i = 0
    nodemap = Dict{Int,Int}()
    for u in keys(g)
        nodemap[u] = i
        println(io,"\tnode")
        println(io,"\t[")
        println(io,"\t\tid $i")
        println(io,"\t\tlabel $u")
        length(c) > 0 && println(io,"\t\tvalue $(c[u])")
        println(io,"\t]")
        i += 1
    end
    for u in keys(g)
        for v in g[u]
            if u < v
                println(io,"\tedge")
                println(io,"\t[")
                println(io,"\t\tsource $(nodemap[u])")
                println(io,"\t\ttarget $(nodemap[v])")
                println(io,"\t]")
            end
        end
    end
    println(io, "]")
    close(io)
    return 1
end


text:
####
#Debugging
####
#The following functions are primarily meant for debugging.
#Not used in computations.

#Uses a discretization in the primal space
#h is a 1D function defining moment constraint for the original problem
#S, M, mu are the original parameters
function _vopp_ub_primal(S, M, mu, h; N=1000)
	#translate everything to the Vc space
	c = mu * (1 - M)
	hbar(t) = h((mu - c) * t + c)
	Sc = vopp.comp_Sc(S, M)

	#discretize the space
	@assert isfinite(Sc) "Currently only supports finite S"

	grid = range(0, stop=Sc, length=N)  #grid includes 0 and Sc
	m = Model(solver=GurobiSolver(OutputFlag=0))
	@variable(m, probs[1:N] >= 0)
	@constraint(m, sum(probs) == 1)  #normalization
	@constraint(m, sum([probs[i] * grid[i] for i = 1:N]) == 1)  #mean

	@constraint(m, sum( [hbar(grid[i]) * probs[i] for i = 1:N] ) == 0)  #moment
	@variable(m, y)

	for i = 1:N
		@constraint(m, y >= grid[i] * sum(probs[i:N]))  #pricing inequality
	end
	@objective(m, Min, y)
	solve(m)
	1/getvalue(y)
end

#Uses a discretization in the primal space of rectangular distributions
#S, M, mu, mode are the original parameters of NON-standardized dist
function _vopp_lb_primal_unimodal_MAD(mu, S, M, D, mode; N=1000, pj=-1)
	#convert to standardized distribution
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	h(t) = M * abs(t - 1)/2 - D

	#H(t) = 1/abs(mode_c - t) * int_<t, mode_c> h if t != mode_c
	function H(t)
		TOL = 1e-8
	    #Dirac case
	    if abs(t - mode_c) <= TOL
	        return h(mode_c)
	    end

	    if 1 >= min(t, mode_c) && 1 <= max(t, mode_c)
	    	return quadgk(h, t, 1, mode_c)[1] / (mode_c - t)
	    end

	    return quadgk(h, t, mode_c)[1]/ (mode_c - t)
	end

	function G(p, m, t)
		if abs(m - t) <= 1e-8
			return m >= p ? 1. : 0.
		elseif p > max(m, t)
			return 0.
		elseif p < min(m, t)
			return 1.
		end
		return (max(m, t) - p) / abs(m - t)
	end

	#ensure the dirac is in mixture
	t_grid = collect(range(0, stop=Sc, length=N))
	if !(mode_c in t_grid)
		push!(t_grid, mode_c)
		sort!(t_grid)
	end
	N = length(t_grid)

	#maximize pricing at p
	r_mp = -Inf
	best_price = -Inf
	dist = zero(t_grid)

	if pj < 0
		p_grid = t_grid[:]
	else
		p_grid = [pj]
	end

	for p in p_grid
		m = Model(solver=GurobiSolver(OutputFlag=0))
		@variable(m, probMs[1:N] >= 0)
		@constraint(m, sum(probMs) == 1)  #normalization

		#mean
		@constraint(m, dot(t_grid .+ mode_c, probMs) == 2)

		#moment
		@constraint(m, dot(map(H, t_grid), probMs) == 0)

		@objective(m, Max, p * dot(map(t-> G(p, mode_c, t), t_grid),
									probMs)
					)

		status = solve(m)
		if status != :Optimal
			println("Status: $status for price $p  mode $mode_c")
		elseif getobjectivevalue(m) > r_mp
			r_mp = getobjectivevalue(m)
			best_price = p
			dist[:] = getvalue(probMs)
		end
	end #end loop over prices
	return 1/r_mp, best_price, t_grid, dist
end

#NOTE:  Worst-Case distribution is a dirac at 0 (which is mode) and a uniform <0, S>
#As S goes to infinity, you can achieve a D of 2.
#returns maximum coef of dev, mode and also the mixture distribution
#recall all mixtures written as <t, m> uniforms
#Inputs should be normalized
##VG Possibly Deprecate?
function max_D_unimodal(S; N=100)
	t_grid = range(0, stop=S, length=N)

	h(t) = abs(t-1)
	function H(t, m)
		TOL = 1e-8
		if abs(t-m) <= TOL
			return h(m)
		end
		if 1 <= min(t, m) || 1 >= max(t, m)
			return quadgk(s-> h(s), t, m)[1] / (m-t)
		else
			return quadgk(s-> h(s), t, 1, m)[1] / (m-t)
		end
		#unreachable
	end

	max_D = -Inf
	max_mode = -Inf
	dist = zero(t_grid)
	for mode in range(0, stop=2, length=N)
		m = Model(solver=GurobiSolver(OutputFlag=0))
		@variable(m, probMs[1:N] >= 0)
		@constraint(m, sum(probMs) == 1)  #normalization

		#mean
		@constraint(m, dot(t_grid .+ mode, probMs) == 2)

		#max coef of dev
		@objective(m, Max, .5 * dot(map(t->H(t, mode), t_grid), probMs))
		status = solve(m)
		if status != :Optimal
			println("Status: $status for mode $mode")
		elseif getobjectivevalue(m) > max_D
			max_D = getobjectivevalue(m)
			max_mode = mode
			dist[:] = getvalue(probMs)
		end
	end #loop over mode

	return max_D, max_mode, dist, t_grid
end

#Uses a discretization in the primal space of rectangular distributions
#S, M, mu, mode are the original parameters of NON-standardized dist
function _vopp_lb_primal_unimodal_CV(mu, S, M, C, mode; N=1000, price=-1)
	#convert to standardized distribution
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	h(t) = M^2 * (t-1)^2 - C^2

	#H(t) = 1/abs(mode_c - t) * int_<t, mode_c> h if t != mode_c
	H(t) = abs(mode_c - t) <= 1e-8 ? h(mode_c) : quadgk(h, t, mode_c)[1]/ (mode_c - t)

	function G(p, m, t)
		if abs(m - t) <= 1e-8
			return m >= p ? 1. : 0.
		elseif p > max(m, t)
			return 0.
		elseif p < min(m, t)
			return 1.
		end
		return (max(m, t) - p) / abs(m - t)
	end

	#ensure the dirac is in mixture
	t_grid = collect(range(0, stop=Sc, length=N))
	if !(mode_c in t_grid)
		push!(t_grid, mode_c)
		sort!(t_grid)
	end
	N = length(t_grid)

	#maximize pricing at p
	r_mp = -Inf
	best_price = -Inf
	dist = zero(t_grid)

	if price > 0
		p_grid = [price]
	else
		p_grid = t_grid
	end

	for p in p_grid
		m = Model(solver=GurobiSolver(OutputFlag=0))
		@variable(m, probMs[1:N] >= 0)
		@constraint(m, sum(probMs) == 1)  #normalization

		#mean
		@constraint(m, dot(t_grid .+ mode_c, probMs) == 2)

		#moment
		@constraint(m, dot(map(H, t_grid), probMs) == 0)

		@objective(m, Max, p * dot(map(t-> G(p, mode_c, t), t_grid),
									probMs)
					)

		status = solve(m)
		if status != :Optimal
			println("Status: $status for price $p  mode $mode_c")
		elseif getobjectivevalue(m) > r_mp
			r_mp = getobjectivevalue(m)
			best_price = p
			dist[:] = getvalue(probMs)
		end
	end #end loop over prices
	return 1/r_mp, best_price, dist
end

#conjecture that maximal CV occurs for uniform
function max_cv_unimodal_guess(mu, S, M, mode)
	#convert to standardized distribution
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	if Sc < 2  #Unif[2-Sc, Sc]
		return Sc / sqrt(12)
	else
		return 2 / sqrt(12)
	end
end


text:
function fuse_batchnorm(d::Dense, bn::BatchNorm)
    d.σ != identity &&
        throw(ArgumentError("Can't fuse Dense layer with activation $(d.σ)."))
    scale = safedivide(bn.γ, sqrt.(bn.σ²))
    W = scale .* d.weight
    b = scale .* (d.bias - bn.μ) + bn.β
    return Dense(W, b, bn.λ)
end

function fuse_batchnorm(c::Conv, bn::BatchNorm)
    c.σ != identity && throw(ArgumentError("Can't fuse Conv layer with activation $(c.σ)."))
    scale = safedivide(bn.γ, sqrt.(bn.σ²))
    W = c.weight .* reshape(scale, 1, 1, 1, :)
    b = scale .* (c.bias - bn.μ) + bn.β
    return Conv(W, b, bn.λ)
end

"""
    try_fusing(model, i)

Attempt to fuse pair of model layers at indices `i` and `i+1`.
Returns fused model and `true` if layers were fused, unmodified model and `false` otherwise.
"""
function try_fusing(model, i)
    l1 = model[i]
    l2 = model[i + 1]
    if l1 isa Union{Dense,Conv} && l2 isa BatchNorm && activation(l1) == identity
        if i == length(model) - 1
            model = Chain(model[1:(i - 1)]..., fuse_batchnorm(l1, l2))
        end
        model = Chain(model[1:(i - 1)]..., fuse_batchnorm(l1, l2), model[(i + 2):end]...)
        return model, true
    end
    return model, false
end

"""
    canonize(model)

Canonize model by flattening it and fusing BatchNorm layers into preceding Dense and Conv
layers with linear activation functions.
"""
function canonize(model::Chain)
    model = flatten_model(model)
    i = 1
    while i < length(model)
        model, fused = try_fusing(model, i)
        !fused && (i += 1)
    end
    return model
end


text:
"""
	TransferOperatorRectangularBinning

A transfer operator estimated over rectangular binning of a set of points, using the
`TransferOperatorEstimatorRectangularBinning` estimator.

The return type for the `TransferOperatorEstimatorRectangularBinning` estimator.

## Fields
-  **`transfermatrix`**. The estimated transfer matrix.

"""
struct TransferOperatorRectangularBinning <: AbstractTransferOperator
	transfermatrix::AbstractArray{Float64, 2}
end

export TransferOperatorRectangularBinning


text:
module VarianceComponentTypeTest

using VarianceComponentModels, MathProgBase, Ipopt, Base.Test

srand(123)

# generate data from a d-variate response variane component model
n = 100   # no. observations
d = 2     # no. categories
m = 2     # no. variance components
p = 2     # no. covariates
X = randn(n, p)
B = ones(p, d)
Σ = ntuple(x -> zeros(d, d), m)
for i in 1:m
  Σi = randn(d, d)
  copy!(Σ[i], Σi' * Σi)
end
## make the first variance component 0 matrix
#fill!(Σ[1], 0.0)
V = ntuple(x -> zeros(n, n), m)
for i = 1:m-1
  Vi = randn(n, 50)
  copy!(V[i], Vi * Vi')
end
copy!(V[m], eye(n))
# form Ω
Ω = zeros(n*d, n*d)
for i = 1:m
  Ω += kron(Σ[i], V[i])
end
Ωchol = cholfact(Ω)
Y = X * B + reshape(Ωchol[:L] * randn(n*d), n, d)

info("VarianceComponentModel type")
#@code_warntype VarianceComponentModel(B, Σ)
@inferred VarianceComponentModel(B, Σ)
#@code_warntype VarianceComponentModel(Σ)
@inferred VarianceComponentModel(Σ)
vcmodel = VarianceComponentModel(Σ)
@test isempty(vcmodel.B)

info("TwoVarCompModelRotate type")
F = eigfact(Σ[1], Σ[2])
#@code_warntype TwoVarCompModelRotate(B * F.vectors, F.values, F.vectors, logdet(Σ[2]))
@inferred TwoVarCompModelRotate(B * F.vectors, F.values, F.vectors, logdet(Σ[2]))

info("VarianceComponentVariate type")
#@code_warntype VarianceComponentVariate(Y, V)
@inferred VarianceComponentVariate(Y, V)
vcdata = VarianceComponentVariate(Y, V)
@test isempty(vcdata.X)

info("Forming VarianceComponentModel from data")
vcdata = VarianceComponentVariate(Y, X, V)
#@code_warntype VarianceComponentVariate(Y, X, V)
@inferred VarianceComponentVariate(Y, X, V)
vcdata = VarianceComponentVariate(Y, X, V)
#@code_warntype VarianceComponentModel(vcdata)
@inferred VarianceComponentModel(vcdata)
vcmodel = VarianceComponentModel(vcdata)
@test vecnorm(vcmodel.B) ≈ 0.0
@test vecnorm(vcmodel.Σ[1] - eye(d)) ≈ 0.0
@test vecnorm(vcmodel.Σ[2] - eye(d)) ≈ 0.0

info("Pre-compute eigen-decomposition and rotate data")
#@code_warntype TwoVarCompVariateRotate(vcdata)
@inferred TwoVarCompVariateRotate(vcdata)
vcdatarot = TwoVarCompVariateRotate(vcdata)
#@code_warntype TwoVarCompModelRotate(vcmodel)
@inferred TwoVarCompModelRotate(vcmodel)
vcmodelrot = TwoVarCompModelRotate(vcmodel)

info("Test generalized eigen-decomposition")
vcdatacopy = deepcopy(vcdata)
scale!(vcdatacopy.V[2], 2.0) # V[2] = 2I
vcdatacopyrot = TwoVarCompVariateRotate(vcdatacopy)
@test vcdatacopyrot.logdetV2 ≈ n * log(2.0) # det(V[2]) = 2^n

info("VarianceComponentModel from TwoVarCompVariateRotate")
@inferred VarianceComponentModel(vcdatarot)
vcmodelfromrot = VarianceComponentModel(vcdatarot)
@test size(vcmodelfromrot.B, 1) == p
@test size(vcmodelfromrot.B, 2) == d
@test length(vcmodelfromrot.Σ) == m

info("VarianceComponentAuxData from VarianceComponentVariate")
vcobsaux1 = VarianceComponentAuxData(vcdata)
vcobsaux2 = VarianceComponentAuxData(vcdatarot)
@test size(vcobsaux1.res) == size(vcobsaux2.res)
@test size(vcobsaux1.Xwork) == size(vcobsaux2.Xwork)
@test size(vcobsaux1.ywork) == size(vcobsaux2.ywork)
@test size(vcobsaux1.obswt) == size(vcobsaux2.obswt)

info("Query functions")
@test eltype(vcmodel) == eltype(B)
@test eltype(vcdata) == eltype(Y)
@test eltype(vcmodelrot) == eltype(vcmodel)
@test eltype(vcdatarot) == eltype(vcdata)
@test length(vcmodel) == d
@test length(vcdata) == d
@test length(vcmodelrot) == d
@test length(vcdatarot) == d
@test size(vcdata) == (n, d)
@test size(vcdatarot) == (n, d)
@test nvarcomps(vcmodel) == m
@test nvarcomps(vcdata) == m
@test nvarcomps(vcmodelrot) == m
@test nvarcomps(vcdatarot) == m
@test nmeanparams(vcmodel) == p * d
@test nmeanparams(vcmodelrot) == p * d
@test nmeanparams(vcdata) == p * d
@test nmeanparams(vcdatarot) == p * d
@test nvarparams(vcmodel) == m * binomial(d + 1, 2)
@test nvarparams(vcmodelrot) == m * binomial(d + 1, 2)
@test nparams(vcmodel) == p * d + m * binomial(d + 1, 2)
@test nparams(vcmodelrot) == p * d + m * binomial(d + 1, 2)

info("Mean, covariance, and residual of model")
#@code_warntype mean(vcmodel, vcdata)
@inferred mean(vcmodel, vcdata)
μ = mean(vcmodel, vcdata)
@test size(μ, 1) == n
@test size(μ, 2) == d
#@code_warntype cov(vcmodel, vcdata)
@inferred cov(vcmodel, vcdata)
Ω = cov(vcmodel, vcdata)
@test size(Ω, 1) == n * d
@test size(Ω, 2) == n * d
#@code_warntype residual(vcmodel, vcdata)
@inferred residual(vcmodel, vcdata)
res = residual(vcmodel, vcdata)
@test size(res, 1) == n
@test size(res, 2) == d
@inferred residual(vcmodelrot, vcdatarot)
resrot = residual(vcmodelrot, vcdatarot)
@test size(resrot, 1) == n
@test size(resrot, 2) == d

end # module VarianceComponentTypeTest


text:
# single channel EEG

using Neuroimaging, DataDeps

data_path = joinpath(
    datadep"ExampleSSR",
    "Neuroimaging.jl-example-data-master",
    "neuroimaingSSR.bdf",
)

s = read_EEG(data_path)

plot(s, "F6")

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
function save_noise!(W::NoiseFunction)

end

Base.show(io::IO, A::NoiseFunction) =
           invoke(show, Tuple{typeof(io), Any}, io, A)
Base.show(io::IO, ::MIME"text/plain", A::NoiseFunction) = show(io, A)

function interpolate!(W::NoiseFunction,t)
  W.W(t)
end

function interpolate!(out1,out2,W::NoiseFunction,t)
  W.W(out1,out2,t)
end

function calculate_step!(W::NoiseFunction,dt)
  if isinplace(W)
    W(W.dW,W.dZ,W.curt+dt)
    W.dW .-= W.curW
    if W.Z != nothing
      W.dZ .-= W.curZ
    end
  else
    new_W, new_Z = W(W.curt+dt)
    W.dW = new_W - W.curW
    if W.Z != nothing
      W.dZ = new_Z - W.curZ
    end
  end
  W.dt = dt
end

function accept_step!(W::NoiseFunction,dt,setup_next=true)
  if isinplace(W)
    W.curW .+= W.dW
  else
    W.curW += W.dW
  end
  W.curt += W.dt
  if W.Z != nothing
    if isinplace(W)
      W.curZ .+= W.dZ
    else
      W.curZ += W.dZ
    end
  end

  W.dt = dt #dtpropose
  if setup_next
    calculate_step!(W,dt)
  end
end

function reject_step!(W::NoiseFunction,dtnew)
  calculate_step!(W,dtnew)
end

function setup_next_step!(W::NoiseFunction)
  calculate_step!(W,W.dt)
end


text:
function Obj_crit_mix(model::EoSModel,F,z,V,T)
    f(x) = eos(model,V,T,x)
    H(x) = ForwardDiff.hessian(f,x)/8.134/T
    L(x) = det(H(x))
    dL(x) = ForwardDiff.gradient(L,x)
    M(x) = [H(x)[1:end-1,:];transpose(dL(x))]
    F[1] = L(z)
    F[2] = det(M(z))
    return F
end

function crit_mix(model::EoSModel,z;v0=nothing)
    if v0 === nothing
        v0 = x0_crit_mix(model,z)
    end
    f! = (F,x) -> Obj_crit_mix(model, F, z, exp10(x[1]), x[2])
    r  = Solvers.nlsolve(f!,v0[1:end],LineSearch(Newton()))
    sol = Solvers.x_sol(r)
    T_c = sol[2]
    V_c = exp10(sol[1])
    p_c = pressure(model, V_c, T_c, z)
    return (T_c, p_c, V_c)
end

function x0_crit_mix(model::EoSModel,z)
    pure = split_model(model)
    crit = crit_pure.(pure)
    n_c  = length(pure)
    V_c  = sum(z[i]*crit[i][3] for i in 1:n_c)
    T_c  = prod(crit[i][1]^z[i] for i ∈ 1:n_c)
    return [log10(V_c),T_c]
end


text: const Nullable{T} = Union{Nothing, T}


text:
export JacobianFreeNewtonKrylovSolver, JacobianAction

"""
mutable struct JacobianAction{FT, AT}
    rhs!
    ϵ::FT
    Q::AT
    Qdq::AT
    Fq::AT
    Fqdq::AT
end

Solve for Frhs = F(q), the Jacobian action is computed

    ∂F(Q)      F(Q + eΔQ) - F(Q)
    ---- ΔQ ≈ -------------------
     ∂Q                e

...
# Arguments
- `rhs!`           : nonlinear operator F(Q)
- `ϵ::FT`          : ϵ used for finite difference, e = e(Q, ϵ)
- `Q::AT`          : cache for Q
- `Qdq::AT`        : container for Q + ϵΔQ
- `Fq::AT`         : cache for F(Q)
- `Fqdq::AT`       : container for F(Q + ϵΔQ)
...
"""
mutable struct JacobianAction{FT, AT}
    rhs!
    ϵ::FT
    Q::AT
    Qdq::AT
    Fq::AT
    Fqdq::AT
end

function JacobianAction(rhs!, Q, ϵ)
    return JacobianAction(
        rhs!,
        ϵ,
        similar(Q),
        similar(Q),
        similar(Q),
        similar(Q),
    )
end

"""
Approximates the action of the Jacobian of a nonlinear
form on a vector `ΔQ` using the difference quotient:

      ∂F(Q)      F(Q + e ΔQ) - F(Q)
JΔQ = ---- ΔQ ≈ -------------------
       ∂Q                e

Compute  JΔQ with cached Q and F(Q), and the direction  dQ
"""
function (op::JacobianAction)(JΔQ, dQ, args...)
    rhs! = op.rhs!
    Q = op.Q
    Qdq = op.Qdq
    ϵ = op.ϵ
    Fq = op.Fq
    Fqdq = op.Fqdq

    FT = eltype(dQ)
    n = length(dQ)
    normdQ = norm(dQ, weighted_norm)

    if normdQ > ϵ
        factor = FT(1 / (n * normdQ))
    else
        # initial newton step, ΔQ = 0
        factor = FT(1 / n)
    end

    β = √ϵ
    e = factor * β * norm(Q, 1, false) + β

    Qdq .= Q .+ e .* dQ

    rhs!(Fqdq, Qdq, args...)

    JΔQ .= (Fqdq .- Fq) ./ e

end

"""
update cached Q and F(Q) before each Newton iteration
"""
function update_Q!(op::JacobianAction, Q, args...)
    op.Q .= Q
    Fq = op.Fq

    op.rhs!(Fq, Q, args...)
end

"""
Solve for Frhs = F(Q), by finite difference

    ∂F(Q)      F(Q + eΔQ) - F(Q)
    ---- ΔQ ≈ -------------------
     ∂Q                e

     Q^n+1 = Q^n - dF/dQ(Q^{n})⁻¹ (F(Q^n) - Frhs)

     set ΔQ = F(Q^n) - Frhs
"""
mutable struct JacobianFreeNewtonKrylovSolver{FT, AT} <: AbstractNonlinearSolver
    # small number used for finite difference
    ϵ::FT
    # tolerances for convergence
    tol::FT
    # Max number of Newton iterations
    M::Int
    # Linear solver for the Jacobian system
    linearsolver
    # container for unknows ΔQ, which is updated for the linear solver
    ΔQ::AT
    # contrainer for F(Q)
    residual::AT
end

"""
JacobianFreeNewtonKrylovSolver constructor
"""
function JacobianFreeNewtonKrylovSolver(
    Q,
    linearsolver;
    ϵ = 1.e-8,
    tol = 1.e-6,
    M = 30,
)
    FT = eltype(Q)
    residual = similar(Q)
    ΔQ = similar(Q)
    return JacobianFreeNewtonKrylovSolver(
        FT(ϵ),
        FT(tol),
        M,
        linearsolver,
        ΔQ,
        residual,
    )
end

"""
JacobianFreeNewtonKrylovSolver initialize the residual
"""
function initialize!(
    rhs!,
    Q,
    Qrhs,
    solver::JacobianFreeNewtonKrylovSolver,
    args...,
)
    # where R = Qrhs - F(Q)
    R = solver.residual
    # Computes F(Q) and stores in R
    rhs!(R, Q, args...)
    # Computes R = R - Qrhs
    R .-= Qrhs
    return norm(R, weighted_norm)
end

"""
Solve for Frhs = F(Q), by finite difference

Q^n+1 = Q^n - dF/dQ(Q^{n})⁻¹ (F(Q^n) - Frhs)

set ΔQ = F(Q^n) - Frhs

...
# Arguments
- `rhs!`:  functor rhs!(Q) =  F(Q)
- `jvp!`:  Jacobian action jvp!(ΔQ)  = dF/dQ(Q) ⋅ ΔQ
- `preconditioner`: approximation of dF/dQ(Q)
- `Q` : Q^n
- `Qrhs` : Frhs
- `solver`: linear solver
...
"""
function donewtoniteration!(
    rhs!,
    jvp!,
    preconditioner::AbstractPreconditioner,
    Q,
    Qrhs,
    solver::JacobianFreeNewtonKrylovSolver,
    args...,
)

    FT = eltype(Q)
    ΔQ = solver.ΔQ
    ΔQ .= FT(0.0)

    # R(Q) == 0, R = F(Q) - Qrhs, where F = rhs!
    # Compute right-hand side for Jacobian system:
    # J(Q)ΔQ = -R
    # where R = Qrhs - F(Q), which is computed at the end of last step or in the initialize function
    R = solver.residual

    # R = F(Q^n) - Frhs
    # ΔQ = dF/dQ(Q^{n})⁻¹ (Frhs - F(Q^n)) = -dF/dQ(Q^{n})⁻¹ R
    iters =
        linearsolve!(jvp!, preconditioner, solver.linearsolver, ΔQ, -R, args...)

    # Newton correction Q^{n+1} = Q^n + dF/dQ(Q^{n})⁻¹ (Frhs - F(Q^n))
    Q .+= ΔQ

    # Compute residual norm and residual for next step
    rhs!(R, Q, args...)
    R .-= Qrhs
    resnorm = norm(R, weighted_norm)

    return resnorm, iters
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
include("help.jl")
function mark_kross_x(r::Robot)
    for side in ((Nord,Ost),(Sud,Ost),(Sud,West),(Nord,West))
        while (isborder(r,side[1])==false && isborder(r,side[2])==false)
            move!(r,side[1])
            move!(r,side[2])
            putmarker!(r)
        end
        while (ismarker(r)==true)
            move!(r, inverse(side[1]))
            move!(r, inverse(side[2]))
        end
    end
    putmarker!(r)
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
module eve
using Reexport
include("SDE.jl")
include("MarketDatas.jl")
include("CREST.jl")
include("EveCentral.jl")
include("LocalData.jl")
include("PriceTools.jl")

@reexport using .SDE
@reexport using .MarketDatas
@reexport using .CREST
@reexport using .EveCentral
@reexport using .LocalData
@reexport using .PriceTools
end # module


text:
module Mjolnir

using IRTools, MacroTools
using IRTools.All
using IRTools: block

export @trace

include("context.jl")
include("cleanup.jl")
include("infer.jl")
include("trace.jl")
include("macros.jl")

include("lib/base.jl")
include("lib/numeric.jl")
include("lib/array.jl")
include("lib/struct.jl")

end # module


text:
using Test
using CartesianMesh
#using Revise
using CutCell
include("useful_routines.jl")

mesh = UniformMesh([0.,0.],[1.,1.],[2,2])
@test allequal(CutCell.elements_per_mesh_side(mesh),[2,2])
cellmaps = CutCell.cell_maps(mesh)
cm1 = CutCell.CellMap([0.,0.],[0.5,0.5])
cm2 = CutCell.CellMap([0.,0.5],[0.5,1.])
cm3 = CutCell.CellMap([0.5,0.],[1.,0.5])
cm4 = CutCell.CellMap([0.5,0.5],[1.,1.])
@test allequal(cellmaps,[cm1,cm2,cm3,cm4])

coords = CutCell.nodal_coordinates(mesh.x0,mesh.widths,mesh.nelements,[5,5])
xrange = range(0.0,stop=1.,length=5)
ycoords = repeat(xrange,5)
xcoords = repeat(xrange,inner=5)
@test allapprox(coords,vcat(xcoords',ycoords'))

connectivity = CutCell.nodal_connectivity([5,5],3,9,mesh.nelements)
c1 = [1,2,3,6,7,8,11,12,13]
c2 = [3,4,5,8,9,10,13,14,15]
c3 = [11,12,13,16,17,18,21,22,23]
c4 = [13,14,15,18,19,20,23,24,25]
testconn = hcat(c1,c2,c3,c4)
@test allequal(connectivity,testconn)

cellconnectivity = CutCell.cell_connectivity(mesh)
testcellconn = [0  1  0  3
                3  4  0  0
                2  0  4  0
                0  0  1  2]
@test allequal(testcellconn,cellconnectivity)

mesh = UniformMesh([0.,0.],[1.,1.],[3,2])
connectivity = CutCell.nodal_connectivity([10,7],4,16,mesh.nelements)
testnodeconn = vcat(43:46,50:53,57:60,64:67)
@test allequal(connectivity[:,5],testnodeconn)

cellconnectivity = CutCell.cell_connectivity(mesh)
c1 = [0,3,2,0]
c2 = [1,4,0,0]
c3 = [0,5,4,1]
c4 = [3,6,0,2]
c5 = [0,0,6,3]
c6 = [5,0,0,4]
testconn = hcat(c1,c2,c3,c4,c5,c6)
@test allequal(cellconnectivity,testconn)

mesh = UniformMesh([0.,0.],[1.,1.],[3,2])
femesh = CutCell.Mesh(mesh,16)
testnodeconn = vcat(43:46,50:53,57:60,64:67)
@test allequal(testnodeconn,femesh.nodalconnectivity[:,5])
@test CutCell.nodes_per_element(femesh) == 16
@test CutCell.nodes_per_mesh_side(femesh) == [10,7]
@test CutCell.number_of_nodes(femesh) == 70

@test allequal(CutCell.bottom_boundary_node_ids(femesh),1:7:64)
@test allequal(CutCell.right_boundary_node_ids(femesh),64:70)
@test allequal(CutCell.top_boundary_node_ids(femesh),7:7:70)
@test allequal(CutCell.left_boundary_node_ids(femesh),1:7)

mesh = CutCell.Mesh([0.,0.],[4.,3.],[4,3],9)
isboundarycell = CutCell.is_boundary_cell(mesh)
testboundarycell = ones(Int,12)
testboundarycell[5] = testboundarycell[8] = 0
@test allequal(isboundarycell,testboundarycell)

isboundarycell = CutCell.is_boundary_cell(mesh.cellconnectivity)
@test allequal(isboundarycell,testboundarycell)

mesh = CutCell.Mesh([2.,1.],[2.,1.],[3,2],4)
@test CutCell.cell_id(mesh,[3.,1.25]) == 3
@test CutCell.cell_id(mesh,[10/3+eps(),1.5+eps()]) == 6


text:
############################################################################################
"""
$(SIGNATURES)
Function that checks for number of particles to achieve target variance of log target estimate.

# Examples
```julia
```

"""
function estimate_Nparticles(
    _rng::Random.AbstractRNG,
    pf::ParticleFilterConstructor,
    objective::Objective,
    variance::AbstractFloat;
    Nchains::Int64 = Threads.nthreads(),
    margin::Float64 = 0.25,
    itermax::Int64 = 100,
    mincoverage::Float64 = 0.01,
    printoutput = true
)
    ArgCheck.@argcheck Nchains > 0
    ArgCheck.@argcheck margin > 0
    ArgCheck.@argcheck itermax > 0
## Assign all kernels
    algorithms = map(iter -> pf(_rng, objective.model, objective.data, objective.temperature, SampleDefault()), Base.OneTo(Nchains))
    models = map(iter -> deepcopy(objective.model), Base.OneTo(Nchains))
## Preallocate buffer for variance and mean of ℓobjective estimate
    ℓobjective_variance = zeros(itermax)
    ℓobjective = zeros(Nchains, itermax)
    Ncoverage = zeros(itermax)
    Ncoverage[1] = pf.default.coverage
    iter = 1
## Run first test for estimates
    for chain in Base.OneTo(Nchains)
        _, diagnostics = propose!(_rng, algorithms[chain], models[chain], objective.data, objective.temperature, UpdateTrue())
        ℓobjective[chain, iter] = diagnostics.base.ℓobjective
    end
    ℓobjective_variance[iter] = Statistics.var(view(ℓobjective, :, iter))
    while iter < itermax
        iter += 1
    ## Compute variance and mean of ℓobjective estimate
        Threads.@threads for chain in Base.OneTo(Nchains)
            _, diagnostics = propose!(_rng, algorithms[chain], models[chain], objective.data, objective.temperature, UpdateTrue())
            ℓobjective[chain, iter] = diagnostics.base.ℓobjective
        end
        ℓobjective_variance[iter] = var(view(ℓobjective, :, iter))
    ## Check if conditions are fullfilled
        _margin = ℓobjective_variance[iter]*margin
        boundaries = (ℓobjective_variance[iter] - _margin, ℓobjective_variance[iter] + _margin)
        if boundaries[begin] < variance < boundaries[end]
            break
        end
    ## Adjust number of particles
        if ℓobjective_variance[iter] > variance
            #Note: Double particles at the beginning and then slowly decrease change over iterations
            Ncoverage[iter] = Ncoverage[iter-1] + Ncoverage[iter-1] * 2.0/iter
        else
            Ncoverage[iter] = max(mincoverage, Ncoverage[iter-1] - Ncoverage[iter-1] * 2.0/iter)
        end
        for chain in Base.OneTo(Nchains)
            algorithms[chain].tune.chains.coverage = Ncoverage[iter]
        end
    ## Print current run diagnostics
        if printoutput
            println(
                "(Target variance ", variance,
                ", Iteration ", iter,
                ") - Coverage set to ", round(Ncoverage[iter]; digits=3),
                " for variance of ", round(ℓobjective_variance[iter]; digits=2), "."
            )
        end
    end
    ## Return Nparticles
    return (
        Coverage = Ncoverage,
        ℓobjective_variance = ℓobjective_variance,
        ℓobjective = ℓobjective,
        itermax = iter,
        variance = variance
        )
end

############################################################################################
#export
export estimate_Nparticles


text:
"""
    Multivector2{T <: Real} <: AbstractMultivector{T}

An immutable 4-tuple of real numbers that represents a member
of a 2-dimensional multivector algebra.

Each `Multivector2` has the form
```math
    p+qB
```
where ``p`` and ``q`` are 1-dimensional multivectors, or
```math
    a+bA+cB+dAB
```
where ``a``, ``b``, ``c``, and ``d`` are real (and of the same type),
and ``A * A = 0``, ``B * B = 0``, and ``AB = A * B = -B * A``.
Here ``*`` is the wedge product.
"""
struct Multivector2{T <: Real} <: AbstractMultivector{T}
    l::Multivector1{T}
    r::Multivector1{T}

    Multivector2{U}(l::Multivector1{U}, r::Multivector1{U}) where {U <: Real} = new(l, r)
end

function Multivector2(x::Multivector1{T}, y::Multivector1{T}) where {T <: Real}
    Multivector2{T}(x, y)
end

function Multivector2(z::Multivector1)
    Multivector2(z, zero(z))
end

function Multivector2(a::T, b::T, c::T, d::T) where {T <: Real}
    Multivector2{T}(Multivector1{T}(a, b), Multivector1{T}(c, d))
end

function Multivector2(a::T, b::T, c::T) where {T <: Real}
    Multivector2{T}(Multivector1{T}(a, b), Multivector1{T}(c, zero(T)))
end

function Multivector2(a::T, b::T) where {T <: Real}
    Multivector2{T}(Multivector1{T}(a, b), zero(Multivector1{T}))
end

function Multivector2(a::T) where {T <: Real}
    Multivector2{T}(Multivector1{T}(a, zero(T)), zero(Multivector1{T}))
end

function Multivector2(a::Real, b::Real, c::Real, d::Real)
    Multivector2(promote(a, b, c, d)...)
end

function Multivector2(a::Real, b::Real, c::Real)
    Multivector2(promote(a, b, c)...)
end

function Multivector2(a::Real, b::Real)
    Multivector2(promote(a, b)...)
end

function Multivector2(x::Multivector1, y::Multivector1)
    Multivector2(promote(x.l, x.r, y.l, y.r)...)
end

function show(io::IO, z::Multivector2)
    print(io, "[1: ")
    print(io, z.l.l)
    print(io, ", A: ")
    print(io, z.l.r)
    print(io, ", B: ")
    print(io, z.r.l)
    print(io, ", AB: ")
    print(io, z.r.r)
    print(io, "]")
end

function real(z::Multivector2)
    real(z.l)
end

function unreal(z::Multivector2)
    vcat(unreal(z.l), asarray(z.r))
end

function zero(z::Multivector2{T}) where {T <: Real}
    Multivector2{T}(zero(z.l), zero(z.r))
end

function zero(::Type{Multivector2{T}}) where {T <: Real}
    Multivector2(zero(T), zero(T), zero(T), zero(T))
end

function one(z::Multivector2{T}) where {T <: Real}
    Multivector2{T}(one(z.l), zero(z.r))
end

function one(::Type{Multivector2{T}}) where {T <: Real}
    Multivector2(one(T), zero(T), zero(T), zero(T))
end

"""
    conj{T <: Real}(z::Multivector2{T})

The `Multivector2` conjugate. If ``z=a+bA+cB+dAB``, then `conj(z)` gives
```math
    a-bA-cB-dAB
```
This operation is an involution.
"""
function conj(z::Multivector2{T}) where {T <: Real}
    Multivector2{T}(conj(z.l), -z.r)
end

"""
    cloak{T <: Real}(z::Multivector2{T})

The cloak conjugate changes the sign of even blades. If ``z=a+bA+cB+dAB``, then `cloak(z)` gives
```math
    -a+bA+cB-dAB
```
This operation is equivalent to `-dagger(z)` and thus is also an involution.
"""
function cloak(z::Multivector2{T}) where {T <: Real}
    Multivector2{T}(cloak(z.l), dagger(z.r))
end

"""
    dagger{T <: Real}(z::Multivector2{T})

The dagger conjugate changes the sign of odd blades. If ``z=a+bA+cB+dAB``, then `dagger(z)` gives
```math
    a-bA-cB+dAB
```
This operation is an involution.
"""
function dagger(z::Multivector2{T}) where {T <: Real}
    Multivector2{T}(dagger(z.l), cloak(z.r))
end

"""
    star{T <: Real}(z::Multivector2{T})

Returns the Hodge star conjugate. If ``z=a+bA+cB+dAB``, then `star(z)` gives
```math
    d-cA+bB+aAB
```
This operation is not an involution, but `star(star(z)) = dagger(z)`.
"""
function star(z::Multivector2{T}) where {T <: Real}
    Multivector2{T}(dagger(star(z.r)), star(z.l))
end

"""
    selfstar(z::Multivector2)

The self-star-conjugate part.
This operation is idempotent.
"""
function selfstar(z::Multivector2)
    (z + star(z) + dagger(z) + star(dagger(z))) / 4
end

"""
    antiselfstar(z::Multivector2)

The anti-self-star-conjugate part.
This operation is idempotent.
"""
function antiselfstar(z::Multivector2)
    (z - star(z) + dagger(z) - star(dagger(z))) / 4
end

function (+)(x::Multivector2, y::Multivector2)
    Multivector2(x.l + y.l, x.r + y.r)
end

function (+)(x::Multivector2, y::Multivector1)
    Multivector2(x.l + y, x.r)
end

function (+)(x::Multivector1, y::Multivector2)
    Multivector2(x + y.l, y.r)
end

function (+)(z::Multivector2, a::Real)
    Multivector2(z.l + a, z.r)
end

function (+)(a::Real, z::Multivector2)
    Multivector2(z.l + a, z.r)
end

function (-)(z::Multivector2{T}) where {T <: Real}
    Multivector2{T}(-z.l, -z.r)
end

function (-)(x::Multivector2, y::Multivector2)
    Multivector2(x.l - y.l, x.r - y.r)
end

function (-)(x::Multivector2, y::Multivector1)
    Multivector2(x.l - y, x.r)
end

function (-)(x::Multivector1, y::Multivector2)
    Multivector2(x - y.l, -y.r)
end

function (-)(z::Multivector2, a::Real)
    Multivector2(z.l - a, z.r)
end

function (-)(a::Real, z::Multivector2)
    Multivector2(a - z.l, -z.r)
end

"""
    (*)(x::Multivector2, y::Multivector2)

Wedge product of 2-dimensional multivectors.
This operation is non-commutative but associative.
"""
function (*)(x::Multivector2, y::Multivector2)
    Multivector2(x.l * y.l,  (y.r * x.l) + (x.r * conj(y.l)))
end

function (*)(x::Multivector2, y::Multivector1)
    Multivector2(x.l * y, x.r * conj(y))
end

function (*)(x::Multivector1, y::Multivector2)
    Multivector2(x * y.l,  y.r * x)
end

"""
    (*)(z::Multivector2, a::Real)
    (*)(a::Real, z::Multivector2)

Scaling and/or reflection of a `Multivector2` by a real number.
"""
function (*)(z::Multivector2, a::Real)
    Multivector2(z.l * a, z.r * a)
end

function (*)(a::Real, z::Multivector2)
    Multivector2(a * z.l, a * z.r)
end

function (/)(x::Multivector2, y::Multivector2)
    if iszerodivisor(y)
        error(ZeroDivisorDenominator)
    end

    x * inv(y)
end

function (\)(y::Multivector2, x::Multivector2)
    if iszerodivisor(y)
        error(ZeroDivisorDenominator)
    end

    inv(y) * x
end

function (/)(z::Multivector2, a::Real)
    Multivector2(z.l / a, z.r / a)
end

function (\)(a::Real, z::Multivector2)
    Multivector2(a \ z.l, a \ z.r)
end

function random(::Type{Multivector2{T}}) where T <: Real
    Multivector2{T}(random(Multivector1{T}), random(Multivector1{T}))
end


text:
using Gcode
using Base.Test

def setUp(self):
    self.outfile = TemporaryFile()
    self.g = G(outfile=self.outfile, print_lines=False,
               aerotech_include=False)

def tearDown(self):
    self.g.teardown()
    del self.outfile
    del self.g

def test_init(self):
    self.assertEqual(self.g.movement_mode, 'relative')

def test_set_home(self):
    g = self.g
    g.set_home()
    expected = 'G92'
    self.assert_output(expected)
    g.set_home(x=10, y=20, A=5)
    expected += """
G92 X10.000000 Y20.000000 A5.000000"""
    self.assert_output(expected)
    self.assert_position({'A': 5.0, 'x': 10.0, 'y': 20.0, 'z': 0})
    g.set_home(y=0)
    self.assert_position({'A': 5.0, 'x': 10.0, 'y': 0.0, 'z': 0})

def test_reset_home(self):
    self.g.reset_home()
    expected = 'G92.1'
    self.assert_output(expected)

def test_relative(self):
    self.g.relative()
    self.assertEqual(self.g.movement_mode, 'relative')
    self.assert_output('G91')

def test_absolute(self):
    self.g.absolute()
    self.assertEqual(self.g.movement_mode, 'absolute')
    self.assert_output('G90')

def test_feed(self):
    self.g.feed(10)
    self.assert_output('F10')

def test_dwell(self):
    self.g.dwell(10)
    self.assert_output('G4 P10')

def test_setup(self):
    self.outfile = TemporaryFile()
    self.g = G(outfile=self.outfile, print_lines=False)
    expected = open(os.path.join(HERE, '../header.txt')).read()
    self.assert_output(expected)

def test_home(self):
    self.g.home()
    expected = """
G90
G1 X0.000000 Y0.000000
G91
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 0})

def test_move(self):
    g = self.g
    g.move(10, 10)
    self.assert_position({'x': 10.0, 'y': 10.0, 'z': 0})
    g.move(10, 10, A=50)
    self.assert_position({'x': 20.0, 'y': 20.0, 'A': 50, 'z': 0})
    expected = """
G1 X10.000000 Y10.000000
G1 X10.000000 Y10.000000 A50.000000
"""
    self.assert_output(expected)

def test_abs_move(self):
    self.g.abs_move(10, 10)
    expected = """
G90
G1 X10.000000 Y10.000000
G91
"""
    self.assert_output(expected)
    self.assert_position({'x': 10, 'y': 10, 'z': 0})

    self.g.abs_move(5, 5)
    expected += """
G90
G1 X5.000000 Y5.000000
G91
"""
    self.assert_output(expected)
    self.assert_position({'x': 5, 'y': 5, 'z': 0})

    self.g.abs_move(15, 0, D=5)
    expected += """
G90
G1 X15.000000 Y0.000000 D5.000000
G91
"""
    self.assert_output(expected)
    self.assert_position({'x': 15, 'y': 0, 'D': 5, 'z': 0})

def test_arc(self):
    with self.assertRaises(RuntimeError):
        self.g.arc()

    self.g.arc(x=10, y=0)
    expected = """
G17
G2 Y0 X10 R5.0
"""
    self.assert_output(expected)
    self.assert_position({'x': 10, 'y': 0, 'z': 0})

    self.g.arc(x=5, A=0, direction='CCW', radius=5)
    expected += """
G16 X Y A
G18
G3 A0 X5 R5
"""
    self.assert_output(expected)
    self.assert_position({'x': 15, 'y': 0, 'A': 0, 'z': 0})

    self.g.arc(x=0, y=10, helix_dim='D', helix_len=10)
    expected += """
G16 X Y D
G17
G2 Y10 X0 R5.0 G1 D10
"""
    self.assert_output(expected)
    self.assert_position({'x': 15, 'y': 10, 'A': 0, 'D': 10, 'z': 0})

    with self.assertRaises(RuntimeError):
        self.g.arc(x=10, y=10, radius=1)

def test_abs_arc(self):
    self.g.abs_arc(x=0, y=10)
    expected = """
G90
G17
G2 Y10 X0 R5.0
G91
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 10, 'z': 0})

    self.g.abs_arc(x=0, y=10)
    expected += """
G90
G17
G2 Y10 X0 R0.0
G91
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 10, 'z': 0})

def test_rect(self):
    self.g.rect(10, 5)
    expected = """
G1 Y5.000000
G1 X10.000000
G1 Y-5.000000
G1 X-10.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 0})

    self.g.rect(10, 5, start='UL')
    expected += """
G1 X10.000000
G1 Y-5.000000
G1 X-10.000000
G1 Y5.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 0})

    self.g.rect(10, 5, start='UR')
    expected += """
G1 Y-5.000000
G1 X-10.000000
G1 Y5.000000
G1 X10.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 0})

    self.g.rect(10, 5, start='LR')
    expected += """
G1 X-10.000000
G1 Y5.000000
G1 X10.000000
G1 Y-5.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 0})

    self.g.rect(10, 5, start='LL', direction='CCW')
    expected += """
G1 X10.000000
G1 Y5.000000
G1 X-10.000000
G1 Y-5.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 0})

    self.g.rect(10, 5, start='UL', direction='CCW')
    expected += """
G1 Y-5.000000
G1 X10.000000
G1 Y5.000000
G1 X-10.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 0})

    self.g.rect(10, 5, start='UR', direction='CCW')
    expected += """
G1 X-10.000000
G1 Y-5.000000
G1 X10.000000
G1 Y5.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 0})

    self.g.rect(10, 5, start='LR', direction='CCW')
    expected += """
G1 Y5.000000
G1 X-10.000000
G1 Y-5.000000
G1 X10.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 0})

def test_meander(self):
    self.g.meander(2, 2, 1)
    expected = """
G91
G1 X2.000000
G1 Y1.000000
G1 X-2.000000
G1 Y1.000000
G1 X2.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 2, 'y': 2, 'z': 0})

    self.g.meander(2, 2, 1.1)
    expected += """
;WARNING! meander spacing updated from 1.1 to 1.0
G91
G1 X2.000000
G1 Y1.000000
G1 X-2.000000
G1 Y1.000000
G1 X2.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 4, 'y': 4, 'z': 0})

    self.g.meander(2, 2, 1, start='UL')
    expected += """
G91
G1 X2.000000
G1 Y-1.000000
G1 X-2.000000
G1 Y-1.000000
G1 X2.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 6, 'y': 2, 'z': 0})

    self.g.meander(2, 2, 1, start='UR')
    expected += """
G91
G1 X-2.000000
G1 Y-1.000000
G1 X2.000000
G1 Y-1.000000
G1 X-2.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 4, 'y': 0, 'z': 0})

    self.g.meander(2, 2, 1, start='LR')
    expected += """
G91
G1 X-2.000000
G1 Y1.000000
G1 X2.000000
G1 Y1.000000
G1 X-2.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 2, 'y': 2, 'z': 0})

    self.g.meander(2, 2, 1, start='LR', orientation='y')
    expected += """
G91
G1 Y2.000000
G1 X-1.000000
G1 Y-2.000000
G1 X-1.000000
G1 Y2.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 4, 'z': 0})

    self.g.meander(3, 2, 1, start='LR', orientation='y')
    expected += """
G91
G1 Y2.000000
G1 X-1.000000
G1 Y-2.000000
G1 X-1.000000
G1 Y2.000000
G1 X-1.000000
G1 Y-2.000000
"""
    self.assert_output(expected)
    self.assert_position({'x': -3, 'y': 4, 'z': 0})

def test_clip(self):
    self.g.clip()
    expected = """
G16 X Y Z
G18
G3 X0 Z4 R2.0
"""
    self.assert_output(expected)
    self.assert_position({'y': 0, 'x': 0, 'z': 4})

    self.g.clip(axis='A', direction='-y', height=10)
    expected += """
G16 X Y A
G19
G2 Y0 A10 R5.0
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 4, 'A': 10})

    self.g.clip(axis='A', direction='-y', height=-10)
    expected += """
G16 X Y A
G19
G3 Y0 A-10 R5.0
"""
    self.assert_output(expected)
    self.assert_position({'x': 0, 'y': 0, 'z': 4, 'A': 0})

### helper functions #####################################################

def assert_output(self, expected):
    if isinstance(expected, basestring):
        expected = expected.split('\n')
    expected = [x.strip() for x in expected if x.strip()]
    self.outfile.seek(0)
    lines = self.outfile.readlines()
    lines = [x.strip() for x in lines if x.strip()]
    self.assertListEqual(lines, expected)

def assert_position(self, expected):
    self.assertEqual(self.g.current_position, expected)


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    WeightedSampler(size, [weights]; replace=false)

A method for weighted sampling from spatial objects that produces
samples of given `size` based on `weights` with or without replacement
depending on the option `replace`. By default weights are uniform.
"""
struct WeightedSampler{W<:Union{Vector{<:Real},Nothing}} <: AbstractSampler
  size::Int
  weights::W
  replace::Bool
end

WeightedSampler(size::Int, weights::W=nothing;
                replace=false) where {W<:Union{Vector{<:Real},Nothing}} =
  WeightedSampler(size, weights, replace)

WeightedSampler(size::Int, weights::AbstractWeights; replace=false) =
  WeightedSampler(size, collect(weights), replace)

function sample(object::AbstractSpatialObject, sampler::WeightedSampler)
  n = npoints(object)
  s = sampler.size
  r = sampler.replace
  if s > n && r == false
    @error "invalid sample size for sampling without replacement"
  end
  w = isnothing(sampler.weights) ? fill(1/n, n) : sampler.weights
  @assert length(w) == n "invalid number of weights for object"

  view(object, sample(1:n, Weights(w), s, replace=r))
end


text:
# List of all valid NYC zipcodes

const NYC_ZIPCODES = [
    10453, 10457, 10460, # Central Bronx
    10458, 10467, 10468, # Bronx Park and Fordham
    10451, 10452, 10456, # High Bridge and Morrisania
    10454, 10455, 10459, 10474, # Hunts Point and Mott Haven
    10463, 10471, # Kingsbridge and Riverdale
    10466, 10469, 10470, 10475, # Northeast Bronx
    10461, 10462, 10464, 10465, 10472, 10473, # Southeast Bronx

    11212, 11213, 11216, 11233, 11238, # Central Brooklyn
    11209, 11214, 11228, # Southwest Brooklyn
    11204, 11218, 11219, 11230, # Borough Park
    11234, 11236, 11239, # Canarsie and Flatlands
    11223, 11224, 11229, 11235, # Southern Brooklyn
    11210, 11205, 11215, 11217, 11231, # Northwest Brooklyn
    11203, 11210, 11225, 11226, # Flatbush
    11207, 11208, # East New York and New Lots
    11211, 11222, # Greenpoint
    11220, 11232, # Sunset Park
    11206, 11221, 11237, # Bushwick and Williamsburg

    10026, 10027, 10030, 10037, 10039, # Central Harlem
    10001, 10011, 10018, 10019, 10020, 10036, # Chelsea and Clinton
    10029, 10035, # East Harlem
    10010, 10016, 10017, 10022, # Gramercy Park and Murray Hill
    10012, 10013, 10014, # Greenwich Village and Soho
    10004, 10005, 10006, 10007, 10038, 10280, # Lower Manhattan
    10002, 10003, 10009, # Lower East Side
    10021, 10028, 10044, 10065, 10075, 10128, # Upper East Side
    10023, 10024, 10025, # Upper West Side
    10031, 10032, 10033, 10034, 10040, # Inwood and Washington Heights

    11361, 11362, 11363, 11364, # Northeast Queens
    11354, 11355, 11356, 11357, 11358, 11359, 11360, # North Queens
    11365, 11366, 11367, # Central Queens
    11412, 11423, 11432, 11433, 11434, 11435, 11436, # Jamaica
    11101, 11102, 11103, 11104, 11105, 11106, # Northwest Queens
    11374, 11375, 11379, 11385, # West Central Queens
    11691, 11692, 11693, 11694, 11695, 11697, # Rockaways
    11004, 11005, 11411, 11413, 11422, 11426, 11427, 11428, 11429, # Southeast Queens
    11414, 11415, 11416, 11417, 11418, 11419, 11420, 11421, # Southwest Queens
    11368, 11369, 11370, 11372, 11373, 11377, 11378, # West Queens

    10302, 10303, 10310, # Port Richmond
    10306, 10307, 10308, 10309, 10312, # South Shore
    10301, 10304, 10305, # Stapleton and St. George
    10314 # Mid-Island
]

const POSITIONS = Dict{String, String}(
    "President of the United States" => "President",
    "Vice-President of the United States" => "Vice President",
    "United States Senate" => "Senator",
    "United States House of Representatives NY-03" => "NY-03 Representative",
    "United States House of Representatives NY-05" => "NY-05 Representative",
    "United States House of Representatives NY-06" => "NY-06 Representative",
    "United States House of Representatives NY-07" => "NY-07 Representative",
    "United States House of Representatives NY-08" => "NY-08 Representative",
    "United States House of Representatives NY-09" => "NY-09 Representative",
    "United States House of Representatives NY-10" => "NY-10 Representative",
    "United States House of Representatives NY-11" => "NY-11 Representative",
    "United States House of Representatives NY-12" => "NY-12 Representative",
    "United States House of Representatives NY-13" => "NY-13 Representative",
    "United States House of Representatives NY-14" => "NY-14 Representative",
    "United States House of Representatives NY-15" => "NY-15 Representative",
    "United States House of Representatives NY-16" => "NY-16 Representative",
    "District Attorney - Bronx" => "Bronx District Attorney",
    "District Attorney - Manhattan" => "Manhattan District Attorney",
    "District Attorney - Queens" => "Queens District Attorney",
    "District Attorney - Richmond" => "Richmond District Attorney"
)


text:
module OntarioCOVID19

using HTTP, CSV, DataFrames, Dates, PlotlyBase
export plot_all_plotly, plot_vaccination_data, plotlyplot

# show plots in static HTML
function plotlyplot(p::Plot)
    println("""
    ~~~

    ~~~
    """)
end

function get_summary_data()
    url = "https://data.ontario.ca/dataset/f4f86e54-872d-43f8-8a86-3892fd3cb5e6/resource/ed270bb8-340b-41f9-a7c6-e8ef587e6d11/download/covidtesting.csv"
    bytes = Vector{UInt8}(replace(String(HTTP.get(url).body),".0"=>""))
    return CSV.File(bytes; normalizenames=true) |> DataFrame
end

function get_vaccine_data()
    url = "https://data.ontario.ca/dataset/752ce2b7-c15a-4965-a3dc-397bf405e7cc/resource/8a89caa9-511c-4568-af89-7f2174b4378c/download/vaccine_doses.csv"
    df = CSV.File(HTTP.get(url).body; normalizenames=true) |> DataFrame
    return df
end

function plot_vaccination_data()
    # vaccination data
    df = OntarioCOVID19.get_vaccine_data()::DataFrame
    df = dropmissing(df, :total_doses_administered, disallowmissing=true)::DataFrame

    date = df.report_date::Vector{Date}
    prev_day_admin = coalesce.(df.previous_day_total_doses_administered,0)::Vector{Int}
    total_admin = coalesce.(df.total_doses_administered,0)::Vector{Int}
    total_complete = coalesce.(df.total_individuals_fully_vaccinated,0)::Vector{Int}

    # insert a data row for 2020-12-29 by using previous_day_doses_administered
    insert!(date, 2, date[2]-Day(1))
    insert!(total_admin, 2, total_admin[2] - prev_day_admin[2])
    insert!(prev_day_admin, 2, 0)
    insert!(total_complete, 2, 0)

    # insert missing rows between 24th & 29th
    splice!(date, 2:1, date[2].-Day.(4:-1:1))
    splice!(total_admin, 2:1, fill(total_admin[1],4))
    splice!(prev_day_admin, 2:1, fill(0,4))
    splice!(total_complete, 2:1, fill(0,4))

    # manipulations to get daily data
    daily_dose = diff([0;total_admin])
    daily_date = date .- Day(1)
    daily_complete = diff([0;total_complete])

    # plot
    layout_options = (xaxis_title="Date", yaxis_rangemode="tozero", legend_xanchor="left", legend_x=0.01, paper_bgcolor="rgba(255,255,255,0)", plot_bgcolor="rgba(255,255,255,0)", margin_l=50, margin_r=50, margin_t=50, margin_b=50)

    ont_pop = 14_733_119
    people_vaccinated = total_admin .- total_complete
    dtick = get_dtick(maximum(people_vaccinated))

    t1 = bar(x=date,y=total_complete,name="Fully Vaccinated")
    t2 = bar(x=date,y=total_admin.-total_complete.*2,name="Partially Vaccinated")
    t3 = scatter(x=date,y=people_vaccinated./ont_pop,yaxis="y2",name="",opacity=0,showlegend=false)
    p1 = Plot([t1,t2,t3], Layout(yaxis_title="People Vaccinated",barmode="stack";layout_options...,
        yaxis_dtick=dtick,yaxis2_dtick=dtick/ont_pop,yaxis2_rangemode="tozero",yaxis2_tickformat=",.1%",yaxis2_side="right",yaxis2_overlaying="y",yaxis2_scaleanchor="y",yaxis2_scaleratio=ont_pop,yaxis2_showgrid=false))

    t1 = bar(x=daily_date,y=daily_complete,name="Second Doses")
    t2 = bar(x=daily_date,y=daily_dose.-daily_complete,name="First Doses")
    t3 = scatter(x=daily_date[7:end],y=moving_average(daily_dose,7),name="7-Day Average")
    p2 = Plot([t1,t2,t3], Layout(yaxis_title="Doses",barmode="stack";layout_options...))

    return p1, p2
end

function plot_all_plotly()
    # download data
    df = get_summary_data()::DataFrame
    df = dropmissing(df, :Total_Cases, disallowmissing=true)::DataFrame

    # cases
    date = df.Reported_Date::Vector{Date}
    total_positive = df.Total_Cases::Vector{Int}
    daily_positive = diff([0;total_positive])
    active = coalesce.(df.Confirmed_Positive,0)::Vector{Int}
    recovered = coalesce.(df.Resolved,0)::Vector{Int}
    deceased = coalesce.(df.Deaths,0)::Vector{Int}
    ventilator = coalesce.(df.Num_of_patients_in_ICU_on_a_ventilator_testing_positive,0)::Vector{Int} .+
        coalesce.(df.Num_of_patients_in_ICU_on_a_ventilator_testing_negative,0)
    icu = coalesce.(df.Number_of_patients_in_ICU_testing_positive_for_COVID_19,0)::Vector{Int} .+
        coalesce.(df.Number_of_patients_in_ICU_testing_negative_for_COVID_19,0)::Vector{Int}
    hospitalized = coalesce.(df.Number_of_patients_hospitalized_with_COVID_19,0)::Vector{Int}
    tests = coalesce.(df.Total_patients_approved_for_testing_as_of_Reporting_Date,0)::Vector{Int}
    investigating = coalesce.(df.Under_Investigation,0)::Vector{Int}
    presumed_pos = coalesce.(df.Presumptive_Positive,0)::Vector{Int}
    presumed_neg = coalesce.(df.Presumptive_Negative,0)::Vector{Int}

    investigating .+= presumed_pos .+ presumed_neg
    negative = max.(tests .- total_positive .- investigating, 0) # don't use confirmed neg col since it's not updated

    daily_total = diff([0; total_positive .+ negative .+ investigating])

    positivity = max.(0,daily_positive./daily_total.*100)

    # plot
    layout_options = (xaxis_range=[date[18],date[end]+Day(1)], xaxis_title="Date", yaxis_rangemode="tozero", legend_xanchor="left", legend_x=0.01, paper_bgcolor="rgba(255,255,255,0)", plot_bgcolor="rgba(255,255,255,0)", margin_l=50, margin_r=50, margin_t=50, margin_b=50)

    t1 = bar(x=date,y=active,name="Active")
    t2 = bar(x=date,y=recovered,name="Recovered")
    t3 = bar(x=date,y=deceased,name="Deceased")
    p1 = Plot([t1,t2,t3], Layout(yaxis_title="Cases",barmode="stack";layout_options...))

    t1 = bar(x=date,y=daily_positive,name="Daily")
    t2 = scatter(x=date[7:end],y=moving_average(daily_positive,7),name="7-Day Average")
    p2 = Plot([t1,t2], Layout(yaxis_title="New Cases",barmode="stack";layout_options...))

    t1 = bar(x=date,y=max.(hospitalized.-icu,0),name="Hospitalized, not in ICU")
    t2 = bar(x=date,y=icu.-ventilator,name="ICU without ventilator")
    t3 = bar(x=date,y=ventilator,name="On a ventilator")
    p3 = Plot([t1,t2,t3], Layout(yaxis_title="Hospitalized Cases",barmode="stack";layout_options...))

    t1 = bar(x=date,y=negative,name="Negative")
    t2 = bar(x=date,y=investigating,name="Pending")
    t3 = bar(x=date,y=total_positive,name="Positive")
    p4 = Plot([t1,t2,t3], Layout(yaxis_title="Tests",barmode="stack";layout_options...))

    t1 = bar(x=date,y=positivity,name="Daily")
    t2 = scatter(x=date[7:end],y=moving_average(positivity,7),name="7-Day Average")
    p5 = Plot([t1,t2], Layout(yaxis_title="Positivity Rate (%)";layout_options...))

    t1 = bar(x=date,y=daily_total,name="Daily")
    t2 = scatter(x=date[7:end],y=moving_average(daily_total,7),name="7-Day Average")
    p6 = Plot([t1,t2], Layout(yaxis_title="Daily Tests";layout_options...))

    idx = icu .!= 0
    frac_icu = icu[idx]./hospitalized[idx].*100
    frac_vent = ventilator[idx]./icu[idx].*100
    X = [fill!(similar(frac_icu),1) 1:length(frac_icu)]
    A1 = X*linreg(X,frac_icu)
    A2 = X*linreg(X,frac_vent)

    t1 = scatter(x=date[idx],y=frac_icu,mode="markers",name="Fraction of hospitalized patients in ICU")
    t2 = scatter(x=date[idx],y=frac_vent,mode="markers",name="Fraction of patients in ICU on a ventilator")
    t3 = scatter(x=date[idx],y=A1,name="")
    t4 = scatter(x=date[idx],y=A2,name="")
    p7 = Plot([t1,t2,t3,t4], Layout(yaxis_title="Percentage (%)",yaxis_range=[0,100];layout_options...))

    return p1, p2, p3, p4, p5, p6, p7
end

# backward-looking moving average
moving_average(vs,n) = [sum(@view vs[i-n+1:i])/n for i in n:length(vs)]

# linear regression
linreg(x::Matrix, y::Array) = x\y

function get_dtick(ymax)
    @assert ymax > 1 # method only works for values larger than 1
    units = [1,2,5,10]
    dtick = round(Int,ymax)
    tdiff = ymax/1
    for nticks in 5:9
        Δ = ymax/nticks
        pow = floor(Int, log10(Δ))
        m = round(Int, Δ/10^pow)
        u = units[argmin(abs.(m.-units))]
        d = u*10^pow
        if abs(d-Δ) < tdiff
            tdiff = abs(d-Δ)
            dtick = d
        end
    end
    return dtick
end

end # module


text:
module Lockblock

using ..Ahorn, Maple

sprites = Dict{String, String}(
    "wood" => "objects/door/lockdoor00",
    "temple_a" => "objects/door/lockdoorTempleA00",
    "temple_b" => "objects/door/lockdoorTempleB00",
    "moon" => "objects/door/moonDoor11"
)

const placements = Ahorn.PlacementDict(
    "Locked Door ($(Ahorn.humanizeVariableName(sprite)))" => Ahorn.EntityPlacement(
        Maple.LockBlock,
        "point",
        Dict{String, Any}(
            "sprite" => sprite
        )
    ) for (sprite, texture) in sprites
)

Ahorn.editingOptions(entity::Maple.LockBlock) = Dict{String, Any}(
    "sprite" => sort(collect(keys(sprites)))
)

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

    return Ahorn.Rectangle(x, y, 32, 32)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.LockBlock, room::Maple.Room)
    sprite = get(entity.data, "sprite", "wood")

    if haskey(sprites, sprite)
        Ahorn.drawSprite(ctx, sprites[sprite], 16, 16)
    end
end

end


text:
import GreyDecision.GreyNumbers.GreyNumber
import GreyDecision.GreyNumbers.kernel

@testset "Grey Numbers" begin

    @testset "Kernel" begin
        @test kernel(GreyNumber(0.0, 10.0)) == 5.0
    end

    @testset "Constructor with single argument" begin
        @test GreyNumber(1.0) == GreyNumber{Float64}(1.0, 1.0)
        @test GreyNumber(0.0) == GreyNumber{Float64}(0.0, 0.0)
        @test GreyNumber(1)   == GreyNumber{Int64}(1, 1)
        @test GreyNumber(0)   == GreyNumber{Int64}(0, 0)
    end

    @testset "Data structure - Zero Constructor" begin
        g = GreyNumber()
        @test g.a == 0.0
        @test g.b == 0.0
    end

    @testset "Data structure - Constructor" begin
        g = GreyNumber(2.0, 5.0)
        @test g.a == 2.0
        @test g.b == 5.0

        h = GreyNumber(5.0, 2.0)
        @test g.a == 2.0
        @test g.b == 5.0
    end

    @testset "Whitenization" begin
        g = GreyNumber(10.0, 20.0)
        w = GreyNumbers.whitenizate(g)
        @test w == 15.0

        w2 = GreyNumbers.whitenizate(g, t = 1.0)
        @test w2 == 10.0

        w3 = GreyNumbers.whitenizate(g, t = 0.0)
        @test w3 == 20.0
    end

    @testset "+" begin
        @testset "Grey + Grey" begin
            g1 = GreyNumber(1, 2)
            g2 = GreyNumber(5, 10)
            greysum = g1 + g2
            greysum2 = g2 + g1
            @test greysum isa GreyNumber
            @test greysum2 isa GreyNumber
            @test greysum == greysum2
            @test greysum.a == 6
            @test greysum.b == 12
        end
        @testset "Grey + scalar" begin
            g1 = GreyNumber(1, 2)
            k = 3
            result  = g1 + k
            result2 = k  + g1
            @test result isa GreyNumber
            @test result2 isa GreyNumber
            @test result2 == result
            @test result.a == 4
            @test result.b == 5
        end
    end

    @testset "Additive inverse (negatition)" begin
        g = GreyNumber(10, 20)
        ginv = -g
        @test ginv isa GreyNumber
        @test ginv.a <= ginv.b
        @test ginv.a == -20
        @test ginv.b == -10
    end

    @testset "Substraction" begin
        @testset "Grey - Grey" begin
            g1 = GreyNumber(10, 20)
            g2 = GreyNumber(1, 2)
            gdiff1 = g1 - g2
            gdiff2 = g2 - g1
            @test gdiff1 isa GreyNumber
            @test gdiff2 isa GreyNumber
            @test gdiff1.a <= gdiff1.b
            @test gdiff2.a <= gdiff2.b
            @test gdiff1.a == 8
            @test gdiff1.b == 19
            @test gdiff2.a == -19
            @test gdiff2.b == -8
        end
        @testset "Grey - scalar" begin
            g = GreyNumber(10, 20)
            k = 3
            result  = g - k
            @test result isa GreyNumber
            @test result.a == 7
            @test result.b == 17
        end
        @testset "Scalar - Grey" begin
            g1 = GreyNumber(100.0, 100.0)
            k = 30.0
            g2 = GreyNumber(k, k)
            @test g1 - k  == g1 - g2
            @test k - g1  == g2 - g1
        end
    end

    @testset "Multiplication" begin
        g1 = GreyNumber(3, 4)
        g2 = GreyNumber(5, 10)
        mult = g1 * g2
        mult2 = g2 * g1
        @test mult isa GreyNumber
        @test mult2 isa GreyNumber
        @test mult == mult2
        @test mult.a == 15
        @test mult.b == 40
    end

    @testset "Multiplication by scalar" begin
        @testset "Scalar on the right" begin
            g1 = GreyNumber(10, 20)
            result = g1 * 10
            @test result isa GreyNumber
            @test result.a == 100
            @test result.b == 200
        end
        @testset "Scalar on the left" begin
            g1 = GreyNumber(10, 20)
            result = 10 * g1
            @test result isa GreyNumber
            @test result.a == 100
            @test result.b == 200
        end
        @testset "Negative scalar" begin
            g1 = GreyNumber(10, 20)
            result = -10 * g1
            @test result isa GreyNumber
            @test result.a <= result.b
            @test result.a == -200
            @test result.b == -100
        end
    end

    @testset "Reciprocal (inverse)" begin
        g = GreyNumber(10, 20)
        repic = inv(g)
        @test repic isa GreyNumber
        @test repic.a <= repic.b
        @test repic.a == 1 / 20
        @test repic.b == 1 / 10
    end

    @testset "Division" begin
        @testset "Division of two numbers" begin
            g1 = GreyNumber(3, 4)
            g2 = GreyNumber(5, 10)
            division = g1 / g2
            @test division isa GreyNumber
            @test division.a <= division.b
            @test division.a == 0.3
            @test division.b == 0.8
        end

        @testset "Division by scalar" begin
            @testset "Positive scalar" begin
                g1 = GreyNumber(30, 40)
                result = g1 / 10
                @test result isa GreyNumber
                @test result.a <= result.b
                @test result.a == 3
                @test result.b == 4
            end

            @testset "Negative scalar" begin
                g1 = GreyNumber(30, 40)
                result = g1 / -10
                @test result isa GreyNumber
                @test result.a <= result.b
                @test result.a == -4
                @test result.b == -3
            end
        end
    end

    @testset "Power" begin
        @testset "All positive" begin
            g = GreyNumber(2, 4)
            powered = g^5
            @test powered isa GreyNumber
            @test powered.a == 32
            @test powered.b == 1024
        end

        @testset "Some negative" begin
            g = GreyNumber(-5, -4)
            powered = g^3
            @test powered isa GreyNumber
            @test powered.a <= powered.b
            @test powered.a == -5 * 5 * 5
            @test powered.b == -4 * 4 * 4
        end
    end

    @testset "sqrt" begin
        g = GreyNumber(36, 100)
        gs = sqrt(g)
        @test gs isa GreyNumber
        @test gs.a == 6.0
        @test gs.b == 10.0
    end

    @testset "cbrt" begin
        k = 5.0
        p = 3.0
        g = GreyNumber(k * k * k, p * p * p)
        result = cbrt(g)

        @test result isa GreyNumber
        @test result.a == p
        @test result.b == k
    end

    @testset "abs" begin
        @testset "all positive" begin
            g = GreyNumber(4, 5)
            gabs = abs(g)
            @test gabs isa GreyNumber
            @test gabs == g
        end

        @testset "all negative" begin
            g = GreyNumber(-4, -5)
            gabs = abs(g)
            @test gabs isa GreyNumber
            @test gabs.a == 4
            @test gabs.b == 5
        end
    end

    @testset "< and >" begin
        @test GreyNumber(4, 5) < GreyNumber(7, 8)
        @test GreyNumber(4, 5) < GreyNumber(7, 10)
        @test GreyNumber(4, 5) > GreyNumber(1, 2)
        @test GreyNumber(5, 40) > GreyNumber(4, 5)
        @test GreyNumber(1, 2) < GreyNumber(1,3)
        @test GreyNumber(1, 5) > GreyNumber(1,3)
    end

    @testset "<= and >=" begin
        @test GreyNumber(4, 5) <= GreyNumber(4, 8)
        @test GreyNumber(4, 5) <= GreyNumber(4, 10)
        @test GreyNumber(4, 5) >= GreyNumber(1, 5)
        @test GreyNumber(5, 40) >= GreyNumber(4, 39)
        @test GreyNumber(1, 2) <= GreyNumber(1,3)
        @test GreyNumber(1, 20) >= GreyNumber(1,3)
    end

    @testset "==" begin
        @test GreyNumber(4, 5) == GreyNumber(4, 5)
        @test GreyNumber(4, 5) != GreyNumber(4, 2)
    end

    @testset "min" begin
        nums = [GreyNumber(1, 20), GreyNumber(2, 20), GreyNumber(3, 20)]
        @test minimum(nums) == GreyNumber(1, 20)
    end

    @testset "max" begin
        nums = [GreyNumber(1, 20), GreyNumber(1, 40), GreyNumber(1, 60)]
        @test maximum(nums) == GreyNumber(1, 60)
    end

    @testset "isinteger" begin
        @test isinteger(GreyNumber(5, 6))
    end

    @testset "isreal" begin
        @test isreal(GreyNumber(5.0, 6.0))
    end

    @testset "isinf" begin
        @test isinf(GreyNumber(0.0, Inf64))
    end

    @testset "iszero" begin
        @test iszero(GreyNumber())
        @test iszero(GreyNumber(0, 0))
        @test iszero(GreyNumber(0.0, 0.0))
    end

    @testset "isone" begin
        @test isone(GreyNumber(1, 1))
        @test isone(GreyNumber(1.0, 1.0))
        @test !isone(GreyNumber())
    end

    @testset "iseven" begin
        # Base.iseven is defined for integers
        @test_throws AssertionError iseven(GreyNumber(2.0, 4.0))
        @test iseven(GreyNumber(2, 4))
    end

    @testset "isodd" begin
        # Base.isodd is defined for integers
        @test_throws AssertionError isodd(GreyNumber(1.0, 3.0))
        @test isodd(GreyNumber(1, 3))
    end

    @testset "elttype" begin
        @test eltype(GreyNumber(0.0, 0.0)) == Float64
        @test eltype(GreyNumber(0, 0)) == Int64
    end

    @testset "sum" begin
        @test sum(GreyNumber()) == 0
        @test sum(GreyNumber(0, 0)) == 0
        @test sum(GreyNumber(0.0, 0.0)) == 0.0
        @test sum(GreyNumber(3.0, 4.0)) == 7.0
    end

    @testset "getindex (g[i])" begin
        @testset "Integers" begin
            g = GreyNumber(7, 8)
            @test g[1] isa Int64
            @test g[2] isa Int64
            @test g[1] == 7
            @test g[2] == 8
        end

        @testset "Floats" begin
            g = GreyNumber(7.0, 8.0)
            @test g[1] isa Float64
            @test g[2] isa Float64
            @test g[1] == 7.0
            @test g[2] == 8.0
        end
    end

    @testset "first" begin
        @testset "Integer" begin
            g = GreyNumber(7, 8)
            @test g[1] isa Int64
            @test first(g) == 7
        end

        @testset "Float" begin
            g = GreyNumber(7.0, 8.0)
            @test g[1] isa Float64
            @test first(g) == 7.0
        end
    end

    @testset "last" begin
        @testset "Integer" begin
            g = GreyNumber(7, 9)
            @test last(g) == 9
        end

        @testset "Floats" begin
            g = GreyNumber(7.0, 9.0)
            @test last(g) == 9.0
        end
    end

end # End of testset GreyNumbers


text:
using Distributions
dists = [
    DiscreteUniform(10,20),
    Binomial(10,0.5),
    Geometric(0.5),
    NegativeBinomial(10,0.5),
    Hypergeometric(30, 40, 10),
    Poisson(5.5)]

println("Distribution \t\t\t\t\t\t Parameters \t Support")
reshape([dists ;  params.(dists) ;
		((d)->(minimum(d),maximum(d))).(dists) ],
		length(dists),3)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
    wnorm(x, smooth=0.0; by="invprior")

`x` is each category's number of samples. `by` option provides two methods:\n
+ `prior`: let biger be more biger, e.g. [1,4] -> [0.25,1.0]\n
+ `invprior`: let smaller be biger, e.g. [1,4] -> [1.0,0.25]\n
to return each category's weight. `help?> prior` and `help?> invprior` for details.
As `p(X|S) = P(S|X)P(X)/P(S) ∝  P(S|X)/P(S)`, `P(S)` can be considered as weight
when we're dealing with imbalanced category samples with a neural network. `P(S)`
can be estimated by the proportion of each category. i.e. ⤦\n
           (# samples from S)
    P(S) = ────────────────── ∝  (# samples from S)
            (# all samples)
"""
function wnorm(x, smooth=0.0; by="invprior")
    @assert smooth>=0.0 "smooth >= 0 , but got $smooth"
    by=="invprior" && return invprior(x, smooth)
    by=="prior"    && return prior(x, smooth)
end

"""
    prior(x, smooth=0.0) -> z

let biger be more biger! ⤦

    y = (x .+ smooth*maximum(x))
    z = y ./ sum(y)
    z = z ./ maximum(z)
# Example
    julia> x = [1,100,1000];
    julia> prior(x,0.5)'
    0.192618  0.230681  0.576701
"""
function prior(x, smooth=0.0)
    n = maximum(x) * smooth
    y = x .+ n
    s = sum(y)
    z = y ./ sum(y)
    return z ./ maximum(z)
end

"""
    invprior(x, smooth=0.0) -> z

let smaller be biger! ⤦

    y = 1 ./ (x .+ smooth*maximum(x))
    z = y ./ sum(y)
    z = z ./ maximum(z)
# Example
    julia> x = [1,100,1000];
    julia> invprior(x,0.5)'
     1.0  0.835  0.334
"""
function invprior(x, smooth=0.0)
    n = maximum(x) * smooth
    y = 1 ./ (x .+ n)
    s = sum(y)
    z = y ./ sum(y)
    return z ./ maximum(z)
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta