text:
# Generate a random test problem
rng = MersenneTwister(1234)
B = 3 # number of blocks

block_model = BlockNLPModel()
T = rand(rng, 1:5, B)
for i = 1:B
    model = Model()
    @variable(model, x >= 0, start = rand(rng))
    @objective(model, Min, (x - T[i])^2)
    @constraint(model, x <= 1)
    nlp = MathOptNLPModel(model)
    add_block(block_model, nlp)
end
A = ones(Float64, B)
b = rand(rng, 50*B:100*B) / 100

links = Dict(1 => [A[1]])
for i = 2:B
    links[i] = sparse([A[i]])
end
add_links(block_model, 1, links, b)

@testset "Testing counters" begin
    @test block_model.problem_size.block_counter ≈ B
    @test block_model.problem_size.link_counter ≈ 1
    @test block_model.problem_size.var_counter ≈ B
    @test block_model.problem_size.con_counter ≈ B
end

@testset "Testing BlockNLPModel.jl " begin
    @test length(block_model.blocks) ≈ B
    @test length(block_model.linking_constraints) ≈ 1
    for i = 1:B
        block = block_model.blocks[i]
        @test block.meta.nvar ≈ 1
        @test block.meta.ncon ≈ 1
        @test block.var_idx == i:i
        @test block.con_idx == i:i
        @test block.linking_constraint_id == [1]
    end
    @test block_model.linking_constraints[1].link_map == Dict(1 => [2, 3, 1])
    @test get_linking_matrix(block_model) == sparse([1, 1, 1], [1, 2, 3], [1.0, 1.0, 1.0])
    for i = 1:B
        @test get_linking_matrix_blocks(block_model)[i] == sparse([1], [1], [1.0])
    end
    @test get_rhs_vector(block_model) == [b]
    @test n_constraints(block_model) == 4
end

@testset "Testing FullSpaceModel.jl" begin
    f(x) = sum((x[i] - T[i])^2 for i = 1:B)
    ∇f(x) = [2 * (x[i] - T[i]) for i = 1:B]
    g(x) = vcat(x, get_linking_matrix(block_model) * x)
    J(x) = [1.0 0 0; 0 1.0 0; 0 0 1.0; 1 1 1]
    H(x, y) = [2.0 0 0; 0 2.0 0; 0 0 2.0]

    fm = FullSpaceModel(block_model)
    @test fm.meta.nvar ≈ B
    @test fm.meta.ncon ≈ B + 1
    @test fm.meta.nnzh ≈ B
    @test fm.meta.nnzj ≈ 6

    x = rand(rng, 0:100, B) / 100
    y = rand(rng, 0:100, B + 1) / 100
    @test obj(fm, x) ≈ f(x)
    @test grad(fm, x) ≈ ∇f(x)
    @test cons(fm, x) ≈ g(x)
    @test jac(fm, x) ≈ J(x)
    @test hess(fm, x, y) ≈ tril(H(x, y))

    # Solve the full-space problem with ipopt
    stats = ipopt(fm, print_level = 0)
    @test stats.solution ≈ [0.41, 1, 1]

    # Solve the full-space problem with MadNLP
    stats = madnlp(fm, print_level = MadNLP.WARN)
    @test stats.solution ≈ [0.4099999816058684, 1.000000008750505, 1.000000009643627]
end

@testset "Testing dualized_block.jl" begin
    A = get_linking_matrix_blocks(block_model)
    y = ones(Float64, 1)
    f(x, i) = (x[1] - T[i])^2 + dot(y, A[i], x)
    ∇f(x, i) = 2 * (x[1] - T[i]) .+ A[i]' * y

    for i = 1:B
        x = rand(rng, 0:100, 1) / 100
        dualized_block = DualizedNLPBlockModel(block_model.blocks[i].problem_block, y, A[i])
        @test dualized_block.meta.nvar ≈ 1
        @test dualized_block.meta.ncon ≈ 1
        @test obj(dualized_block, x) ≈ f(x, i)
        @test grad(dualized_block, x) ≈ ∇f(x, i)

        # Solve the block using ipopt
        result = ipopt(dualized_block; print_level = 0)
        @test result.solution[1] ≈ 1.0
    end
end

@testset "Testing augmented_block.jl" begin
    A = get_linking_matrix(block_model)
    b = get_rhs_vector(block_model)
    sol = rand(rng, 0:100, B) / 100
    ρ = 0.1
    y = ones(Float64, 1)

    f(x, i) =
        (x[i] - T[i])^2 +
        dot(y, A[:, block_model.blocks[i].var_idx], x[i]) +
        (ρ / 2) * norm(A * x - b)^2
    ∇f(x, i) =
        2 * (x[i] - T[i]) .+ (
            A[:, block_model.blocks[i].var_idx]' * y .+
            ρ .* A[:, block_model.blocks[i].var_idx]' * (A * x - b)
        )
    H(x, i) =
        2.0 .+
        ρ * A[:, block_model.blocks[i].var_idx]' * A[:, block_model.blocks[i].var_idx]

    for i = 1:B
        x = rand(rng, 0:100, 1) / 100
        current_sol = deepcopy(sol)
        current_sol[i] = x[1]
        augmented_block = AugmentedNLPBlockModel(block_model.blocks[i], y, ρ, A, b, sol)
        @test augmented_block.meta.nvar ≈ 1
        @test augmented_block.meta.ncon ≈ 1
        @test obj(augmented_block, x) ≈ f(current_sol, i)
        @test grad(augmented_block, x) ≈ ∇f(current_sol, i)
        @test hess(augmented_block, x, [1.0]) ≈ H(x, i)

        # Solve the block using ipopt
        result = ipopt(augmented_block; print_level = 0)
        @test result.solution[1] ≈ 1.0
    end
end

@testset "Testing prox_augmented_block.jl" begin
    A = get_linking_matrix(block_model)
    b = get_rhs_vector(block_model)
    sol = rand(rng, 0:100, B) / 100
    ρ = 0.1
    y = ones(Float64, 1)
    P = sparse(0.1 * I, 1, 1)

    f(x, i) =
        (x[i] - T[i])^2 +
        dot(y, A[:, block_model.blocks[i].var_idx], x[i]) +
        (ρ / 2) * norm(A * x - b)^2 +
        1 / 2 * dot((x[i] - sol[i]), P, (x[i] - sol[i]))
    ∇f(x, i) =
        2 * (x[i] - T[i]) .+ (
            A[:, block_model.blocks[i].var_idx]' * y .+
            ρ .* A[:, block_model.blocks[i].var_idx]' * (A * x - b)
        ) .+ P * (x[i] - sol[i])
    H(x, i) =
        2.0 .+
        ρ * A[:, block_model.blocks[i].var_idx]' * A[:, block_model.blocks[i].var_idx] .+ P

    for i = 1:B
        x = rand(rng, 0:100, 1) / 100
        current_sol = deepcopy(sol)
        current_sol[i] = x[1]
        prox_augmented_block =
            ProxAugmentedNLPBlockModel(block_model.blocks[i], y, ρ, A, b, sol, P)
        @test prox_augmented_block.meta.nvar ≈ 1
        @test prox_augmented_block.meta.ncon ≈ 1
        @test obj(prox_augmented_block, x) ≈ f(current_sol, i)
        @test grad(prox_augmented_block, x) ≈ ∇f(current_sol, i)
        @test hess(prox_augmented_block, x, [1.0]) ≈ H(x, i)

        # Solve the block using ipopt
        result = ipopt(prox_augmented_block; print_level = 0)
        @test result.solution[1] ≈ 1.0
    end
end


text:
using Documenter
using Hankel
import LinearAlgebra: mul!, ldiv!

makedocs(
    sitename="Hankel.jl",
    pages = [
        "Documentation" => "index.md",
        "Derivations" => "derivations.md",
        "Automatic Differentiation" => "autodiff.md"
    ]
)

deploydocs(
    repo = "github.com/chrisbrahms/Hankel.jl.git",
)

ext: jl

max_stars_repo_name: chrisbrahms/Hankel

max_issues_repo_name: chrisbrahms/Hankel

max_forks_repo_name: chrisbrahms/Hankel


text:
function _default(mesh::TOrSignal{M}, s::Style, data::Dict) where M <: GeometryBasics.Mesh
    return @gen_defaults! data begin
        shading = true
        main = mesh
        vertex_color = Vec4f0(0)
        texturecoordinates = Vec2f0(0)
        image = nothing => Texture
        shader = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "standard.vert", "standard.frag",
            view = Dict("light_calc" => light_calc(shading))
        )
    end
end


text:
##### Precipitation tendencies

#####
##### First order fluxes
#####

function flux(::Rain, ::PrecipitationFlux, atmos, args)
    @unpack state, aux = args
    FT = eltype(state)
    u = state.ρu / state.ρ
    q_rai = state.precipitation.ρq_rai / state.ρ

    v_term_rai::FT = FT(0)
    param_set = parameter_set(atmos)
    if q_rai > FT(0)
        v_term_rai =
            terminal_velocity(param_set, CM1M.RainType(), state.ρ, q_rai)
    end

    k̂ = vertical_unit_vector(atmos, aux)
    return state.precipitation.ρq_rai * (u - k̂ * v_term_rai)
end

function flux(::Snow, ::PrecipitationFlux, atmos, args)
    @unpack state, aux = args
    FT = eltype(state)
    u = state.ρu / state.ρ
    q_sno = state.precipitation.ρq_sno / state.ρ
    param_set = parameter_set(atmos)
    v_term_sno::FT = FT(0)
    if q_sno > FT(0)
        v_term_sno =
            terminal_velocity(param_set, CM1M.SnowType(), state.ρ, q_sno)
    end

    k̂ = vertical_unit_vector(atmos, aux)
    return state.precipitation.ρq_sno * (u - k̂ * v_term_sno)
end

#####
##### Second order fluxes
#####

function flux(::Rain, ::Diffusion, atmos, args)
    @unpack state, diffusive = args
    @unpack D_t = args.precomputed.turbulence
    d_q_rai = (-D_t) .* diffusive.precipitation.∇q_rai
    return d_q_rai * state.ρ
end

function flux(::Snow, ::Diffusion, atmos, args)
    @unpack state, diffusive = args
    @unpack D_t = args.precomputed.turbulence
    d_q_sno = (-D_t) .* diffusive.precipitation.∇q_sno
    return d_q_sno * state.ρ
end

#####
##### Sources
#####

function source(::Rain, s::WarmRain_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_qr
end

function source(::Rain, s::RainSnow_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_qr
end

function source(::Snow, s::RainSnow_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_qs
end


text:
module RandomWalk

export randomwalk1d, randomwalk2d, rwtrap1d_lifetime, rwtrap1d_lifetime_exp, dlaline!

"""
1D Random Walk

# Arguments

-`time`: the number of steps of the random walk
-`p`: the probability of going right
-`walks`: the number of walks, defaults to `100*time`

# Returns

a vector containing the coordinates of each point and a vector counting in how many times
the walker landed on each point ath the end of the walk (only non-zero points are returned)
"""
function randomwalk1d(time::Integer, p::Float64; walks=100*time)
    path = zeros(2 * time + 1)
    coordinates = collect(Integer, -time:time)

    for _ in 1:walks
        pos = time + 1
        for _ in 1:time
            if rand() <= p
                pos += 1
            else
                pos -= 1
            end
        end
        path[pos] += 1
    end

    nonzero = findall(x -> x != 0, path)
    return coordinates[nonzero], path[nonzero]
end

"""
2D Random Walk

# Arguments

-`time::Integer`: the number of steps of the random walk
-`walks::Integer`: the number of walks, defaults to the value of `time`

# Returns

two tuples of two integers giving the x and y limits of the points and
an integer matrix counting the number of times the walker landed on each point.
"""
function randomwalk2d(time::Integer; walks::Integer=100*time)
    x, y = collect(Integer, -time:time), collect(Integer, -time:time)
    grid = zeros(Integer, 2 * time + 1, 2 * time + 1)
    steps = ((0, 1), (1, 0), (0, -1), (-1, 0))

    for _ in 1:walks
        xpos, ypos = time + 1, time + 1
        for _ in 1:time
            step = steps[rand(1:4)]
            xpos += step[1]
            ypos += step[2]
        end
        grid[xpos, ypos] += 1
    end

    return randomwalk2dcrop(x, y, grid)
end

"""
crop the output of `randomwalk2d`
"""
@inline function randomwalk2dcrop(x::Vector{<:Integer}, y::Vector{<:Integer},
        grid::Matrix{<:Integer})
    height, width = size(grid)
    bottom = left = 1
    top = height
    right = width

    while all(==(0), grid[top, :])
        top -= 1
    end
    while all(==(0), grid[bottom, :])
        bottom += 1
    end
    while all(==(0), grid[:, left])
        left += 1
    end
    while all(==(0), grid[:, right])
        right -= 1
    end

    return x[left:right], y[bottom:top], grid[bottom:top, left:right]
end

"""
Calculate the average lifetime of a 1D random walker inside trapping boundaries

# Arguments

-`len::Integer`: The length of the path. The coordinates are in the range [-len÷2, len÷2]
-`initpos::Integer`: the initial position of the walker
-`walks::Integer`: the number of walks, defaults to the value of `time`
-`maxtime::Integer`: The maximum number of steps before the walker is abandoned;
    in case the number of steps exceeds `maxtime`, the lifetime is taken to be `maxtime`.
    defaults to `100*len`.

# Returns

The average lifetime, which is the average time passed before the walker gets trapped
in the boundaries of the path
"""
function rwtrap1d_lifetime(len::Integer, initpos::Integer, walks::Integer;
        maxtime::Integer=100*len)
    lifetime_sum = 0

    for _ in 1:walks
        pos = initpos
        time = 0

        while time < maxtime
            time += 1
            if rand() <= 0.5
                pos += 1
                if pos > len ÷ 2
                    break
                end
            else
                pos -= 1
                if pos < -len ÷ 2
                    break
                end
            end
        end

        lifetime_sum += time
    end

    lifetime_avg = lifetime_sum / walks
    return lifetime_avg
end

"""
Calculate the expected lifetime of a 1D random walker inside trapping boundaries

# Arguments

-`len::Integer`: The length of the path. The coordinates are in the range [-len÷2, len÷2]
-`initpos::Integer`: the position of the boundaries (traps)
-`maxtime::Integer`: The maximum number of steps before the walker is abandoned;
    in case the number of steps exceeds `maxtime`, the lifetime is taken to be `maxtime`.

# Returns

The average lifetime, which is the average time passed before the walker gets trapped
in the boundaries of the path
"""
function rwtrap1d_lifetime_exp(len::Integer, initpos::Integer; maxtime::Integer=100*len)
    lifetime_avg = 0

    probs = zeros(len + 3)
    probs[initpos + len ÷ 2 + 2] = 1

    for _ in 1:maxtime
        newprobs = zeros(len + 3)
        for i in 2:len+2
            newprobs[i - 1] += probs[i] * 0.5
            newprobs[i + 1] += probs[i] * 0.5
        end
        probs = newprobs

        lifetime_avg += sum(probs[2:end-1])
    end

    return lifetime_avg
end

"""
Applies diffusion-limited aggregation on a 2D surface, assuming the "seed" is a line

the line is at the bottom and new particles appear from the top of the pile

# Arguments

- `surface::BitMatrix`: The simulation surface. Particles occupy points that are true.
- `particles::Integer`: the number of particles (walkers) simulated
- `escapethresh::Integer`: if the particle gets this much higher from the pile, it will
    "escape". It will be removed from the simulation and a new particle will be generated
    in its place. Defaults to 5.

# Returns
the `surface`
"""
function dlaline!(surface::BitMatrix, particles::Integer; escapethresh::Integer=5)
    height, width = size(surface)
    maxheight = height
    while true ∉ surface[maxheight - 1, :]
        maxheight -= 1
    end
    steps = ((0, 1), (1, 0), (0, -1), (-1, 0))

    while particles > 0
        ypos, xpos = maxheight, rand(1:width)

        while true
            xprev, yprev = xpos, ypos
            step = steps[rand(1:4)]
            ypos += step[1]
            xpos += step[2]
            # periodic boundary conditions (if xpos = width + 1, then it will be shifted
            # left by width, and if xpos = 0, it will be shifted right by width)
            xpos += (width + 1 - 2xpos) ÷ (width + 1) * width

            if ypos > min(maxheight + escapethresh, height)
                break
            elseif surface[ypos, xpos]
                surface[yprev, xprev] = true
                particles -= 1
                if yprev > maxheight
                    maxheight = yprev
                end
                break
            end
        end
    end
    return surface # for convenience
end

end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


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

# simple keyword args case
kwf1(ones; tens=0, hundreds=0) = ones + 10*tens + 100*hundreds
@test kwf1(2) == 2
@test kwf1(2, hundreds=6) == 602
@test kwf1(2, tens=6) == 62
@test kwf1(1, hundreds=2, tens=7) == 271
@test kwf1(3, tens=7, hundreds=2) == 273

@test_throws MethodError kwf1()             # no method, too few args
@test_throws MethodError kwf1(1, z=0)       # unsupported keyword
@test_throws MethodError kwf1(1, 2)         # no method, too many positional args

# keyword args plus varargs
kwf2(x, rest...; y=1) = (x, y, rest)
@test isequal(kwf2(0), (0, 1, ()))
@test isequal(kwf2(0,1,2), (0, 1, (1,2)))
@test isequal(kwf2(0,1,2,y=88), (0, 88, (1,2)))
@test isequal(kwf2(0,y=88,1,2), (0, 88, (1,2)))
@test_throws MethodError kwf2(0, z=1)
@test_throws MethodError kwf2(y=1)

# Test for #13919
test13919(x::Vararg{Int}; key=100) = (x, key)
@test test13919(1, 1)[1] === (1, 1)
@test test13919(1, 1)[2] === 100
@test test13919(1, 1, key=10)[1] === (1, 1)
@test test13919(1, 1, key=10)[2] === 10

# keyword arg with declared type
kwf3(x; y::Float64 = 1.0) = x + y
@test kwf3(2) == 3.0
@test kwf3(2, y=3.0) == 5.0
@test_throws TypeError kwf3(2, y=3)  # wrong type keyword

# function with only keyword args
kwf4(;a=1,b=2) = (a,b)
@test isequal(kwf4(), (1,2))
@test isequal(kwf4(b=10), (1,10))

# in-order evaluation of keyword args
kwf9(;read=true,write=!read) = (read,write)
@test kwf9() === (true,false)
@test kwf9(read=false) === (false,true)
@test kwf9(write=true) === (true,true)

# rest keywords
kwdelegator(ones;kw...) = kwf1(ones;kw...)
@test kwdelegator(4,hundreds=8) == 804

# optional positional args
opaf1(a,b=1,c=2,d=3) = (a,b,c,d)
@test isequal(opaf1(0), (0,1,2,3))
@test isequal(opaf1(0,2), (0,2,2,3))
@test isequal(opaf1(0,2,4), (0,2,4,3))
@test isequal(opaf1(0,2,4,6), (0,2,4,6))
@test_throws MethodError opaf1()
@test_throws MethodError opaf1(0,1,2,3,4)

# optional positional plus varargs
opaf2(a=1,rest...) = (a,rest)
@test isequal(opaf2(), (1,()))
@test isequal(opaf2(2), (2,()))
@test isequal(opaf2(2,3), (2,(3,)))

# optional positional plus keyword args
opkwf1(a=0,b=1;k=2) = (a,b,k)
@test isequal(opkwf1(), (0,1,2))
@test isequal(opkwf1(10), (10,1,2))
@test isequal(opkwf1(10,20), (10,20,2))
@test_throws MethodError opkwf1(10,20,30)
@test isequal(opkwf1(10,20,k=8), (10,20,8))
@test isequal(opkwf1(11;k=8),    (11, 1,8))
@test isequal(opkwf1(k=8),       ( 0, 1,8))

# dictionaries as keywords
@test kwf1(4; Dict(:hundreds=>9, :tens=>5)...) == 954

# with inner function
let
    function kwf_maker()
        f(;k=0) = 2k+1
    end
    kwf5 = kwf_maker()
    @test kwf5() == 1
    @test kwf5(k=2) == 5
    @test_throws MethodError kwf5(1)
end

# with every feature!
extravagant_args(x,y=0,rest...;color="blue",kw...) =
    (x,y,rest,color,kwf1(6;tens=8,kw...))

@test isequal(extravagant_args(1), (1,0,(),"blue",86))
@test isequal(extravagant_args(1;hundreds=7), (1,0,(),"blue",786))
@test isequal(extravagant_args(1,2,3;Dict(:color=>"red", :hundreds=>3)...),
              (1,2,(3,),"red",386))

# passing empty kw container to function with no kwargs
@test sin(1.0) == sin(1.0; Dict()...)
# issue #18845
@test (@eval sin(1.0; $([]...))) == sin(1.0)
f18845() = 2
@test f18845(;) == 2
@test f18845(; []...) == 2
@test (@eval f18845(; $([]...))) == 2

# passing junk kw container
@test_throws BoundsError extravagant_args(1; Any[[]]...)

# keyword args with static parmeters
kwf6{T}(x; k::T=1) = T
@test kwf6(1) === Int
@test kwf6(1;k=2.5) === Float64

kwf7{T}(x::T; k::T=1) = T
@test kwf7(2) === Int
@test kwf7(1.5;k=2.5) === Float64
@test_throws MethodError kwf7(1.5)
@test_throws TypeError kwf7(1.5;k=2)

# try to confuse it with quoted symbol
kwf8{T}(x::MIME{:T};k::T=0) = 0
@test kwf8(MIME{:T}()) === 0

# issue #4538
macro TEST4538()
    quote
        function $(esc(:test4538))(x=1)
            return x
        end
    end
end
@TEST4538
@test test4538() == 1
@test test4538(2) == 2

macro TEST4538_2()
    quote
        function $(esc(:test4538_2))(;x=1)
            return x
        end
    end
end
@TEST4538_2
@test test4538_2() == 1
@test_throws MethodError test4538_2(2)
@test test4538_2(x=2) == 2

# that, but in a module
module Foo4538
macro TEST()
    quote
        function $(esc(:test4538_foo_2))(;x=1)
            return x
        end
    end
end
end
@Foo4538.TEST()
@test test4538_foo_2() == 1
@test test4538_foo_2(x=2) == 2

f4538_3(;x=1) = x
macro TEST4538_3()
    quote
        x = 2
        f4538_3(x=3)
    end
end
@test (@TEST4538_3) == 3

# issue #4801
mutable struct T4801{X}
    T4801{X}(;k=0) where X = new()
end
@test isa(T4801{Any}(k=0), T4801{Any})

# issue #4974
function Foo4974(f; kwargs...)
    (f(), kwargs)
end

function test4974(;kwargs...)
    Foo4974(;kwargs...) do
        2
    end
end

@test test4974(a=1) == (2, [(:a, 1)])

# issue #7704, computed keywords
@test kwf1(1; (:tens, 2)) == 21
let
    p = (:hundreds, 3)
    q = (:tens, 1)
    @test kwf1(0; p, q) == 310
    @test kwf1(0; q, hundreds=4) == 410
end

# with anonymous functions, issue #2773
let f = (x;a=1,b=2)->(x, a, b)
    @test f(0) === (0, 1, 2)
    @test f(1,a=10,b=20) === (1,10,20)
    @test f(0,b=88) === (0, 1, 88)
    @test_throws MethodError f(0,z=1)
end
@test ((a=2)->10a)(3) == 30
@test ((a=2)->10a)() == 20
@test ((a=1,b=2)->(a,b))() == (1,2)
@test ((a=1,b=2)->(a,b))(5) == (5,2)
@test ((a=1,b=2)->(a,b))(5,6) == (5,6)

# issue #9948
f9948, getx9948 = let
    x = 3
    h(;x=x) = x
    getx() = x
    h, getx
end
@test_throws UndefVarError f9948()
@test getx9948() == 3
@test f9948(x=5) == 5
@test_throws UndefVarError f9948()
@test getx9948() == 3

# issue #17785 - handle all sources of kwargs left-to-right
g17785(; a=1, b=2) = (a, b)
let opts = (:a=>3, :b=>4)
    @test g17785(; a=5, opts...) == (3, 4)
    @test g17785(; opts..., a=5) == (5, 4)
    @test g17785(; opts..., a=5, b=6) == (5, 6)
    @test g17785(; b=0, opts..., a=5) == (5, 4)
end

# pr #18396, kwargs before Base is defined
@eval Core.Inference begin
    f18396(;kwargs...) = g18396(;kwargs...)
    g18396(;x=1,y=2) = x+y
end
@test Core.Inference.f18396() == 3

# issue #7045, `invoke` with keyword args
f7045(x::Float64; y=true) = y ? 1 : invoke(f7045,Tuple{Real},x,y=y)
f7045(x::Real; y=true) = y ? 2 : 3
@test f7045(1) === 2
@test f7045(1.0) === 1
@test f7045(1, y=false) === 3
@test f7045(1.0, y=false) === 3

# issue #20804
struct T20804{T}
    y::T
end
(f::T20804)(;x=10) = f.y + x
let x = T20804(4)
    @test x() == 14
    @test x(x=8) == 12
end

# issue #21147
function f21147(f::Tuple{A}; kwargs...) where {B,A<:Tuple{B}}
    return B
end
@test f21147(((1,),)) === Int
@test f21147(((1,),), k = 2) === Int
function g21147(f::Tuple{A}, k = 2) where {B,A<:Tuple{B}}
    return B
end
@test g21147(((1,),)) === Int
@test g21147(((1,),), 2) === Int

max_stars_repo_name: liveinjs/julia

max_issues_repo_name: liveinjs/julia

max_forks_repo_name: liveinjs/julia


text:
module ExtendedVariantModeFlagToggledExtendedVariantTrigger

using ..Ahorn, Maple

# obsolete
@mapdef Trigger "ExtendedVariantMode/FlagToggledExtendedVariantTrigger" FlagToggledExtendedVariantTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	variantChange::String="Gravity", enable::Bool=true, newValue::Integer=10, revertOnLeave::Bool=false, revertOnDeath::Bool=true, withTeleport::Bool=false, flag::String="flag_name", inverted::Bool=false)

end


text:
using Colors
using CoordinateTransformations
using GeometryBasics
using MeshCat
using Rotations

"""
    Hopper
    	model inspired by "Dynamically Stable Legged Locomotion"
		s = (x, z, t, r)
			x - lateral position
			z - vertical position
			t - body orientation
			r - leg length
"""

# Kinematics
kinematics(q) = [q[1] + q[4] * sin(q[3]), q[2] - q[4] * cos(q[3])]
# kinematics(q) = [q[3], q[4]]

# Visualization
function visualize!(vis, q, Δt; fix_camera = true)

	# foot radius
    r_foot = 0.05

	# leg width
    r_leg = 0.5 * r_foot

	# set background to white
	setvisible!(vis["/Background"], true)
    setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setvisible!(vis["/Axes"], false)

	# create body
    setobject!(vis["body"], Sphere(Point3f0(0),
        convert(Float32, 0.1)),
        MeshPhongMaterial(color = RGBA(0, 1, 0, 1.0)))

	# create foot
    setobject!(vis["foot"], Sphere(Point3f0(0),
        convert(Float32, r_foot)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	# create leg
    n_leg = 100
    for i = 1:n_leg
        setobject!(vis["leg$i"], Sphere(Point3f0(0),
            convert(Float32, r_leg)),
            MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
    end

    p_leg = [zeros(3) for i = 1:n_leg]

	# animation
    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(q)
        p_body = [q[t][1], 0.0, q[t][2]]
        p_foot = [kinematics(q[t])[1], 0.0, kinematics(q[t])[2]]

        q_tmp = Array(copy(q[t]))
        r_range = range(0, stop = q[t][4], length = n_leg)
        for i = 1:n_leg
            q_tmp[4] = r_range[i]
            p_leg[i] = [kinematics(q_tmp)[1], 0.0, kinematics(q_tmp)[2]]
        end
        q_tmp[4] = q[t][4]
        p_foot = [kinematics(q_tmp)[1], 0.0, kinematics(q_tmp)[2]]

        z_shift = [0.0; 0.0; r_foot]

        MeshCat.atframe(anim, t) do
            settransform!(vis["body"], Translation(p_body + z_shift))
            settransform!(vis["foot"], Translation(p_foot + z_shift))

            for i = 1:n_leg
                settransform!(vis["leg$i"], Translation(p_leg[i] + z_shift))
            end
        end
    end

	if fix_camera
		settransform!(vis["/Cameras/default"],
			compose(Translation(0.0, 0.5, -1.0), LinearMap(RotZ(-pi / 2.0))))
	end

    MeshCat.setanimation!(vis, anim)
end


text:
# include("../src/jl")
# include("arr_fcn.jl")

import SparseArrays.spzeros

using Distributed
import SparseArrays.dropzeros!
# using JLD
using JLD2
# using FileIO
using Logging
BoolOrNothing = Union{Bool, Nothing}

# @everywhere import SparseArrays.dropzeros!
# @everywhere using JLD
# @everywhere using Logging
# @everywhere using Distributed
# @everywhere BoolOrNothing = Union{Bool, Nothing}

# @everywhere _boundary3_storage = Dict()
_boundary3_storage = Dict();
_larmodel_storage = Dict();
# _global_boundary3_storage = @spawn _boundary3_storage;
_param_boundary_allow_memory = true
_param_boundary_allow_read_files = true
_param_boundary_allow_write_files = true
_param_larmodel_allow_memory = true
_param_larmodel_allow_read_files = false
_param_larmodel_allow_write_files = false

# using arr_fcn

function set_param(;
    boundary_allow_memory::BoolOrNothing=nothing,
    boundary_allow_read_files::BoolOrNothing=nothing,
    boundary_allow_write_files::BoolOrNothing=nothing,
    larmodel_allow_memory::BoolOrNothing=nothing,
    larmodel_allow_read_files::BoolOrNothing=nothing,
    larmodel_allow_write_files::BoolOrNothing=nothing
    )
    global _param_boundary_allow_read_files
    global _param_boundary_allow_write_files
    global _param_boundary_allow_memory
    global _param_larmodel_allow_read_files
    global _param_larmodel_allow_write_files
    global _param_larmodel_allow_memory
    if boundary_allow_read_files != nothing
        _param_boundary_allow_read_files = boundary_allow_read_files
    end
    if boundary_allow_write_files != nothing
        _param_boundary_allow_write_files = boundary_allow_write_files
    end
    if boundary_allow_memory != nothing
        _param_boundary_allow_memory = boundary_allow_memory
    end
    if larmodel_allow_read_files != nothing
        _param_larmodel_allow_read_files = larmodel_allow_read_files
    end
    if larmodel_allow_write_files != nothing
        _param_boundary_allow_write_files = larmodel_allow_write_files
    end
    if larmodel_allow_memory != nothing
        _param_larmodel_allow_memory = larmodel_allow_memory
    end
end

function reset(;
    boundary_storage::BoolOrNothing=nothing,

    )
    if boundary_storage != nothing
        # _boundary3_storage = @spawn Dict()
        _boundary3_storage = Dict()
    end
end

function _create_name_for_boundary(block_size::Array, prefix::String="boundary_matrix_")
    len = length(block_size)
    if len == 3
        fn = prefix * string(block_size[1]) * "x" * string(block_size[2]) * "x"  * string(block_size[3]) * ".jld2"
    else
         error("Function not defined for this dimension")

        fn = ""
    end
    return fn
end

function get_larmodel(block_size::Array)
    @debug "getting larmodel"
    if _param_larmodel_allow_memory & haskey(_larmodel_storage, block_size)
        larmodel = _larmodel_storage[block_size]
    else
        fn = _create_name_for_boundary(block_size::Array, "larmodel_")

        # fn = _create_name_for_boundary(block_size::Array)
        if _param_larmodel_allow_read_files & isfile(fn)
            # larmodel = JLD.load(fn)["larmodel"]
            larmodel = FileIO.load(fn, "larmodel")
            # print("R")
            # @debug "R"
        else
            # println("==== caluculate larmodel")
            larmodel::Lar.LARmodel = Lar.cuboidGrid(block_size, true)
            __fix_cuboidGrid_FV!(larmodel)
            # println("==== boundary calculated")
            if _param_larmodel_allow_write_files
                # JLD.save(fn, "larmodel", larmodel)
                FileIO.save(fn, "larmodel", larmodel)
                # print("W")
                # @debug "W"
            else
                # print("C")
                @debug "C"
            end
        end
        _larmodel_storage[block_size] = larmodel
    end
    return larmodel
end

function get_boundary3(block_size::Array, return_larmodel=true)
    # println("== get_boundary3 function called ", typeof(_boundary3_storage), " ", keys(_boundary3_storage))
    global _global_boundary3_storage
    # _boundary3_storage = fetch(_global_boundary3_storage)
    # println("== fetch finished")
    if _param_boundary_allow_memory & haskey(_boundary3_storage, block_size)
        bMatrix = _boundary3_storage[block_size]
        # println("==== boundary from memory")
        @debug "."
        if return_larmodel
            larmodel = get_larmodel(block_size)
        end

    else
        fn = _create_name_for_boundary(block_size::Array)
        if _param_boundary_allow_read_files & isfile(fn)
            # bMatrix = JLD.load(fn)["boundary_matrix"]
            # bMatrix = FileIO.load(fn, "boundary_matrix")
            @debug "reading boundary matrix from file"
            @JLD2.load fn I J V
            bMatrix = SparseArrays.sparse(I,J,V)
            # print("R")
            @debug "R"
            if return_larmodel
                larmodel = get_larmodel(block_size)
            end
        else
            @debug("Caluculate boundary")
            bMatrix, larmodel = calculate_boundary3(block_size)
            # println("==== boundary calculated")
            if _param_boundary_allow_write_files
                I, J, V = findnz(bMatrix)
                # JLD.save(fn, "boundary_matrix", bMatrix)
                # FileIO.save(fn, "boundary_matrix", bMatrix)
                @JLD2.save fn I J V
                # print("W")
                @debug "W"
            else
                # print("C")
                @debug "C"
            end
        end
        _boundary3_storage[block_size] = bMatrix
        # _boundary3_storage = @spawn local_boundary3_storage;
    end
    # println("== get_boundary3 function end")
    if return_larmodel
        return bMatrix, larmodel
    else
        return bMatrix
    end
#
end

# @everywhere begin
    """
    In cuboidGrid the order for node id in faces is wrong.
    Here is the fix by swapping last two array elements.
    """
    function __fix_cuboidGrid_FV!(larmodel::Lar.LARmodel)
        FV = larmodel[2][3]
        for oneF in FV
            tmp = oneF[4]
            oneF[4] = oneF[3]
            oneF[3] = tmp
        end
    end

    function calculate_boundary3(block_size)
    # function get_boundary3(block_size)
        if typeof(block_size) == Tuple{Int64,Int64,Int64}
            block_size = [block_size[1], block_size[2], block_size[3]]
        end
    #     V, CVill = Lar.cuboidGrid([block_size[1], block_size[2], block_size[3]])

        # A lot of work can be done by this:
        lmodel = get_larmodel(block_size)
        # lmodel::Lar.LARmodel = Lar.cuboidGrid(block_size, true)
        V, (VV, EV, FV, CV) = lmodel
    #     model =
        # __fix_cuboidGrid_FV!(lmodel)

        CVchar = Lar.characteristicMatrix(CV)
        FVchar = Lar.characteristicMatrix(FV)
        b3 = CVchar * FVchar'

        # b3 = sparse_filter!(b3, 4, 1, 0)
        b3 = sparse_filter!(b3, 4, 1, 0)
        dropzeros!(b3)

    #     return b3

    #     model
    #     println("get_boundary3: 3")
        return b3, lmodel #(V, (VV, EV, FV, CV))
    end
# end


text:
using PrettyPlots
using Base.Test

# write your own tests here

# plot an obstacle feild and test the position plot
@test 1 == 1

# obstacle plotting tests
#using NLOptControl

#n=NLOpt();
#c=defineCase();
#obstaclePlot(n,c);

#@test obstaclePlot


text:
struct LArray{T,N,D<:AbstractArray{T,N},Syms} <: DenseArray{T,N}
  __x::D
  LArray{Syms}(__x) where Syms = new{eltype(__x),ndims(__x),typeof(__x),Syms}(__x)
  LArray{T,N,D,Syms}(__x) where {T,N,D,Syms} = new{T,N,D,Syms}(__x)
end

#####################################
# NamedTuple compatibility
#####################################
## LArray to named tuple
function Base.convert(::Type{NamedTuple}, x::LArray{T,N,D,Syms}) where {T,N,D,Syms}
    tup = NTuple{length(Syms),T}(x)
    NamedTuple{Syms,typeof(tup)}(tup)
end
Base.keys(x::LArray{T,N,D,Syms}) where {T,N,D,Syms} = Syms

## Named tuple to LArray
#=
    1. `LArray((2,2), (a=1, b=2, c=3, d=4))` (need to specify size)
    2. `LArray((2,2); a=1, b=2, c=3, d=4)` : alternative form using kwargs
    3. `LVector((a=1, b=2))` : can infer size for vectors
    4. `LVector(a=1, b=2)` : alternative form using kwargs
=#
function LArray(size::NTuple{S,Int}, tup::NamedTuple{Syms,Tup}) where {S,Syms,Tup}
    __x = reshape(collect(tup), size)
    LArray{Syms}(__x)
end
LArray(size::NTuple{S,Int}; kwargs...) where {S} = LArray(size, kwargs.data)
LVector(tup::NamedTuple) = LArray((length(tup),), tup)
LVector(;kwargs...) = LVector(kwargs.data)

## pairs iterator
Base.pairs(x::LArray{T,N,D,Syms}) where {T,N,D,Syms} =
    # (label => getproperty(x, label) for label in Syms) # not type stable?
    (Syms[i] => x[i] for i in 1:length(Syms))

#####################################
# Array Interface
#####################################
Base.size(x::LArray) = size(getfield(x,:__x))
Base.@propagate_inbounds Base.getindex(x::LArray,i...) = getfield(x,:__x)[i...]
Base.@propagate_inbounds function Base.setindex!(x::LArray, y, i...)
  getfield(x,:__x)[i...] = y
  return x
end

Base.propertynames(::LArray{T,N,D,Syms}) where {T,N,D,Syms} = Syms
symnames(::Type{LArray{T,N,D,Syms}}) where {T,N,D,Syms} = Syms

Base.@propagate_inbounds function Base.getproperty(x::LArray,s::Symbol)
    if s == :__x
        return getfield(x,:__x)
    end
    return getindex(x,Val(s))
end

Base.@propagate_inbounds function Base.setproperty!(x::LArray,s::Symbol,y)
    if s == :__x
        return setfield!(x,:__x,y)
    end
    setindex!(x,y,Val(s))
end

Base.@propagate_inbounds Base.getindex(x::LArray, s::Symbol) = getindex(x, Val(s))
Base.@propagate_inbounds Base.getindex(x::LArray, s::Val) = __getindex(x, s)
Base.@propagate_inbounds Base.setindex!(x::LArray, v, s::Symbol) = setindex!(x, v, Val(s))

@generated function Base.setindex!(x::LArray,y,::Val{s}) where s
  syms = symnames(x)
  if syms isa NamedTuple
    idxs = syms[s]
    return quote
      Base.@_propagate_inbounds_meta
      setindex!(getfield(x, :__x), y, $idxs)
      return x
    end
  else # Tuple
    idx = findfirst(y -> y==s, symnames(x))
    return quote
      Base.@_propagate_inbounds_meta
      setindex!(getfield(x, :__x), y, $idx)
      return x
    end
  end
end

Base.@propagate_inbounds Base.getindex(x::LArray,s::AbstractArray{Symbol,1}) = [getindex(x,si) for si in s]

function Base.similar(x::LArray{T,K,D,Syms},::Type{S},dims::NTuple{N,Int}) where {T,K,D,Syms,S,N}
    tmp = similar(x.__x,S,dims)
    LArray{S,N,typeof(tmp),Syms}(tmp)
end

# Allow copying LArray of uninitialized data, as with regular Array
Base.copy(x::LArray) = typeof(x)(copy(getfield(x,:__x)))
Base.deepcopy(x::LArray) = typeof(x)(deepcopy(getfield(x,:__x)))
Base.copyto!(x::LArray,y::LArray) = copyto!(getfield(x,:__x),getfield(y,:__x))

# enable the usage of LAPACK
Base.unsafe_convert(::Type{Ptr{T}}, a::LArray{T,N,D,S}) where {T,N,D,S} = Base.unsafe_convert(Ptr{T}, getfield(a,:__x))

Base.convert(::Type{T},x) where {T<:LArray} = T(x)
Base.convert(::Type{T},x::T) where {T<:LArray} = x
Base.convert(::Type{<:Array},x::LArray) = convert(Array,getfield(x,:__x))

ArrayInterface.restructure(x::LArray{T,N,D,Syms},y::LArray{T2,N2,D2,Syms}) where {T,N,D,T2,N2,D2,Syms} = reshape(y,size(x)...)

#####################################
# Broadcast
#####################################
struct LAStyle{T,N,L} <: Broadcast.AbstractArrayStyle{N} end
LAStyle{T,N,L}(x::Val{i}) where {T,N,L,i} = LAStyle{T,N,L}()
Base.BroadcastStyle(::Type{LArray{T,N,D,L}}) where {T,N,D,L} = LAStyle{T,N,L}()
Base.BroadcastStyle(::LabelledArrays.LAStyle{T,N,L}, ::LabelledArrays.LAStyle{E,N,L}) where{T,E,N,L} =
    LAStyle{promote_type(T,E),N,L}()

function Base.similar(bc::Broadcast.Broadcasted{LAStyle{T,N,L}}, ::Type{ElType}) where {T,N,L,ElType}
    tmp = similar(Array{ElType},axes(bc))
    if axes(bc) != axes(Tuple(L))
        return tmp
    else
        return LArray{ElType,N,typeof(tmp),L}(tmp)
    end
end

# Broadcasting checks for aliasing with Base.dataids but the fallback
# for AbstractArrays is very slow. Instead, we just call dataids on the
# wrapped buffer
Base.dataids(A::LArray) = Base.dataids(A.__x)

"""
    @LArray Eltype Size Names
    @LArray Values Names

Creates an `LArray` with names determined from the `Names`
vector and values determined from the `Values` array. Otherwise, and eltype
and size are used to make an LArray with undefined values.

For example:

    a = @LArray Float64 (2,2) (:a,:b,:c,:d)
    b = @LArray [1,2,3] (:a,:b,:c)
    c = @LArray [1,2,3] (a=1:2, b=2:3)
    d = @LArray [1 2; 3 4] (a=(2, :), b=2:3)
"""
macro LArray(vals,syms)
  vals = esc(vals)
  syms = esc(syms)
  return quote
      LArray{$syms}($vals)
  end
end

macro LArray(type,size,syms)
  type = esc(type)
  size = esc(size)
  syms = esc(syms)
  return quote
      LArray{$syms}(Array{$type}(undef,$size...))
  end
end

"""
    @LVector Type Names

Creates an `LArray` of dimension 1 with eltype and undefined values.
Length is via the number of names given.

For example:

    b = @LVector [1,2,3] (:a,:b,:c)
"""
macro LVector(type,syms)
  type = esc(type)
  syms = esc(syms)
  return quote
      LArray{$syms}(Vector{$type}(undef,length($syms)))
  end
end

#the following gives errror: TypeError: in <:, expected Type, got TypeVar
#symbols(::LArray{T,N,D<:AbstractArray{T,N},Syms}) where {T,N,D,Syms} = Syms

"""
    symbols(::LArray{T,N,D,Syms})

Returns the labels of the `LArray` .

For example:

    z = @LVector Float64 (:a,:b,:c,:d)
    symbols(z)  # NTuple{4,Symbol} == (:a, :b, :c, :d)
"""
symbols(::LArray{T,N,D,Syms}) where {T,N,D,Syms} = Syms isa NamedTuple ? keys(Syms) : Syms

# copy constructors

"""
LVector(v1::Union{SLArray,LArray}; kwargs...)

Creates a 1D copy of v1 with corresponding items in kwargs replaced.

For example:

    z = LVector(a=1, b=2, c=3);
    z2 = LVector(z; c=30)
"""
function LVector(v1::Union{SLArray,LArray}; kwargs...)
  t2 = merge(convert(NamedTuple, v1), kwargs.data)
  LVector(t2)
end

"""
LVector(v1::Union{SLArray,LArray}; kwargs...)

Creates a copy of v1 with corresponding items in kwargs replaced.

For example:

    ABCD = @SLArray (2,2) (:a,:b,:c,:d);
    B = ABCD(1,2,3,4);
    B2 = LArray(B; c=30 )
"""
function LArray(v1::Union{SLArray,LArray}; kwargs...)
  t2 = merge(convert(NamedTuple, v1), kwargs.data)
  LArray(size(v1),t2)
end

# moved vom slarray.js to here because LArray need to be known
"""
SLVector(v1::SLArray; kwargs...)

Creates a 1D copy of v1 with corresponding items in kwargs replaced.

For example:

    z = SLVector(a=1, b=2, c=3);
    z2 = SLVector(z; c=30)
"""
function SLVector(v1::Union{SLArray,LArray}; kwargs...)
  t2 = merge(convert(NamedTuple, v1), kwargs.data)
  SLVector(t2)
end

"""
SLVector(v1::SLArray; kwargs...)

Creates a copy of v1 with corresponding items in kwargs replaced.

For example:

    ABCD = @SLArray (2,2) (:a,:b,:c,:d);
    B = ABCD(1,2,3,4);
    B2 = SLArray(B; c=30 )
"""
function SLArray(v1::Union{SLArray{S,T,N,L,Syms},LArray{T,N,D,Syms}}; kwargs...) where {S,T,N,L,Syms,D}
  t2 = merge(convert(NamedTuple, v1), kwargs.data)
  SLArray{S}(t2)
end


text:
abstract type SimulationModel end

struct MassMatrixModel <: SimulationModel end
struct ImplicitModel <: SimulationModel end


text:
#
# This file is part of the Guards.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

function Actors.onmessage(A::_ACT, ::Val{:guard}, msg::Init)
    A.init = msg.x          # save the init Func
    A.sta = A.init()        # exec the init function
    become(deepcopy, A.sta) # set the behavior
end

function Actors.onmessage(A::_ACT, ::Val{:guard}, msg::Call)
    if isempty(msg.x)
        send(msg.from, Response(A.bhv(), A.self))
    else
        A.res = first(msg.x)(A.sta, msg.x[2:end]...)
        send(msg.from, Response(deepcopy(A.res), A.self))
    end
end

function Actors.onmessage(A::_ACT, ::Val{:guard}, msg::Cast)
    A.res = first(msg.x)(A.sta, msg.x[2:end]...)
end


text:
using DataFrames
using Distributions
using NearestNeighbors
using StaticArrays
using StatsBase

struct ShapePrior{N}
    std_values::MeanVec{N};
    std_value_stds::MeanVec{N};
    n_samples::Int;
end

distributions(prior::ShapePrior) = Normal.(prior.std_values, prior.std_value_stds)
rand(prior::ShapePrior) = rand.(distributions(prior))

var_posterior(prior::ShapePrior{N}, eigen_values::T where T <: Union{MeanVec{N}, StaticArray{Tuple{N},Float64,1}}; n_samples::TR where TR <: Real) where N =
    var_posterior.(prior.std_values, prior.std_value_stds, eigen_values; n_samples=n_samples, prior_n_samples=prior.n_samples)

# f(dx) = sign(dx) * sqrt(|dx/std|) * std
function var_posterior(prior_std::Float64, prior_std_std::Float64, eigen_value::Float64; n_samples::TR1 where TR1 <: Real, prior_n_samples::TR2 where TR2 <: Real)
    d_std = (sqrt(eigen_value) - prior_std)
    std_adj = (prior_std + sign(d_std) * (sqrt(abs(d_std) / prior_std_std + 1) - 1) * prior_std_std)

    return ((prior_n_samples * prior_std + n_samples * std_adj) / (prior_n_samples + n_samples)) ^ 2
end

function sample_var(d::Normal)
    @assert d.μ > 0
    while true
        v = rand(d)
        if v > 0
            return v.^2
        end
    end
end
sample_var(prior::ShapePrior) = sample_var.(distributions(prior))

mutable struct Component{N}
    position_params::MvNormalF{N};
    composition_params::CategoricalSmoothed;
    n_samples::Int;
    prior_probability::Float64;
    confidence::Float64;

    shape_prior::Union{Nothing, ShapePrior{N}};
    n_molecules_per_segment::Dict{Int, Int};

    guid::Int;

    Component(position_params::MvNormalF{L}, composition_params::CategoricalSmoothed;
              n_samples::Int=0, confidence::Float64=1.0, shape_prior::Union{Nothing, ShapePrior{L}}=nothing, guid::Int=-1) where L =
        new{L}(position_params, composition_params, n_samples, 1.0, confidence, shape_prior, Dict{Int, Int}(), guid)
end

# # This function would be useful when processing QC score as another kind of confidences
# function maximize!(c::Component{N} where N, pos_data::T1 where T1 <: AbstractMatrix{Float64}, comp_data::T2 where T2 <: Union{AbstractVector{Int}, AbstractVector{Union{Int, Missing}}},
#         conf_data::T3 where T3 <: AbstractVector{Float64}; nuclei_probs::Union{<:AbstractVector{Float64}, Nothing}=nothing, min_nuclei_frac::Float64=0.1)
#     c.n_samples = size(pos_data, 2) # TODO: need to replace it with confidences, but for that I need to re-write all prior segmentation code to work with confidences as well. Also may need to adjust some other parts
#     maximize!(c.composition_params, comp_data, conf_data);
#     if nuclei_probs === nothing
#         maximize!(c.position_params, pos_data, conf_data);
#     else
#         maximize!(c.position_params, pos_data, conf_data .* nuclei_probs);
#         if length(nuclei_probs) > 1
#             c.confidence = quantile(nuclei_probs, 1 - min_nuclei_frac)
#         end
#     end

#     if c.shape_prior !== nothing
#         try
#             adjust_cov_by_prior!(c.position_params.Σ, c.shape_prior; n_samples=sum(conf_data))
#         catch
#             @show c.position_params.Σ
#             @show c.shape_prior
#             @show sum(conf_data)
#             rethrow()
#         end
#     end

#     return c
# end

function maximize!(c::Component{N} where N, pos_data::T1 where T1 <: AbstractMatrix{Float64}, comp_data::T2 where T2 <: Union{AbstractVector{Int}, AbstractVector{Union{Int, Missing}}};
        nuclei_probs::Union{<:AbstractVector{Float64}, Nothing}=nothing, min_nuclei_frac::Float64=0.1)
    c.n_samples = size(pos_data, 2)
    maximize!(c.composition_params, comp_data);
    if nuclei_probs === nothing
        maximize!(c.position_params, pos_data);
    else
        maximize!(c.position_params, pos_data, nuclei_probs);
        if length(nuclei_probs) > 1
            c.confidence = quantile(nuclei_probs, 1 - min_nuclei_frac)
        end
    end

    if c.shape_prior !== nothing
        adjust_cov_by_prior!(c.position_params.Σ, c.shape_prior; n_samples=c.n_samples)
    end

    return c
end

pdf(comp::Component{N}, x::SVector{N, Float64}, gene::Missing; use_smoothing::Bool=true) where N =
    comp.prior_probability * comp.confidence * pdf(comp.position_params, x)

pdf(comp::Component{N}, x::SVector{N, Float64}, gene::Int64; use_smoothing::Bool=true) where N =
    comp.prior_probability * comp.confidence * pdf(comp.position_params, x) * pdf(comp.composition_params, gene; use_smoothing=use_smoothing)

function adjust_cov_by_prior!(Σ::CovMat{N}, prior::ShapePrior{N}; n_samples::TR where TR <: Real) where N
    if (Σ[2, 1] / max(Σ[1, 1], Σ[2, 2])) < 1e-5 # temporary fix untill https://github.com/JuliaArrays/StaticArrays.jl/pull/694 is merged
        Σ[1, 2] = Σ[2, 1] = 0.0
    end

    fact = eigen(Σ)
    eigen_values_posterior = var_posterior(prior, fmax.(fact.values, 0.0); n_samples=n_samples)
    Σ .= fact.vectors * SMatrix{N, N, Float64}(diagm(0 => eigen_values_posterior)) * inv(fact.vectors)
    Σ .= fmax.(Σ, Σ')

    return adjust_cov_matrix!(Σ)
end

function normal_posterior!(μ::MeanVec{N}, Σ::CovMat{N}, μ_prior::MeanVec{N}, Σ_prior::CovMat{N}; n::TR1 where TR1 <: Real, n_prior::TR2 where TR2 <: Real) where N
    Σ = (n .* Σ .+ n_prior .* Σ_prior .+ (n * n_prior / (n + n_prior)) .* ((μ .- μ_prior) * (μ .- μ_prior)')) ./ (n + n_prior)
    μ .= (n .* μ .+ n_prior .* μ_prior) ./ (n + n_prior)
    adjust_cov_matrix!(Σ)
    return nothing
end


text:
module MatFuncs

# package code goes here

import Base: getindex, size

export MatFunc

struct MatFunc{F, T} <: AbstractMatrix{T}
    func::F
    dims::Tuple{Int,Int}
end

size(fm::MatFunc) = fm.dims
getindex(fm::MatFunc, i::Int, j::Int) = fm.func(i,j)

MatFunc(T::DataType, func::F, dims::Tuple{Int,Int}) where {F<:Function} = MatFunc{F,T}(func,dims)

end # module


text:
module OnlineLearning

using NumericExtensions, NumericFuns, Devectorize, ArrayViews

import NumericExtensions.evaluate

export Learner

abstract Learner
#should implement an update! and a predict! method
# update!{T<:FloatingPoint}(obj::Learner{T}, x::Matrix{T}, y::Vector{T})
#and possibly additional kw args

typealias DSMat{Tv, Ti <: Integer} Union(DenseMatrix{Tv}, SparseMatrixCSC{Tv, Ti})

predict(obj::Learner, x::DSMat; kwargs...) =
	predict!(obj, Array(eltype(x), size(x,1)), x; kwargs...)

include("extra.jl")
include("sgd.jl")
include("glm.jl")
include("svm.jl")
include("utils.jl")

end # module


text:
using NeXLSpectrum
using Test

@testset "XRF Spectrum" begin
    fn1 = joinpath(@__DIR__, "Other", "20170510_D1_Mo_50kv_100mA_acryl 0.spx")
    sp1a = loadspectrum(fn1)
    #sp1b = FileIO.load(fn1)
    @test sp1a[:RealTime] == 310.610
    @test sp1a[:LiveTime] == 300.294
    @test sp1a[:BrukerThroughput] == 130000
    @test sp1a[:DetectorModel] == "Bruker XFlash 430"
    @test sp1a[:DetectorSerialNumber] == "11881_0239"
    @test sp1a[:DetectorThickness] == 0.1 * 0.45 # cm
    @test sp1a[:DeadLayerThickness] == 0.1 * 0.029 # ???
    @test sp1a[:BeamEnergy] == 1000.0 * 50.0
    @test sp1a[:XRFTubeAnode] == n"Mo"
    @test sp1a[:ProbeCurrent] == 1000.0 * 99
    @test sp1a[:XRFTubeIncidentAngle] == deg2rad(84.0)
    @test sp1a[:XRFTubeTakeOffAngle] == deg2rad(6.0)
    @test sp1a[:XRFExcitationAngle] == deg2rad(50.0)
    @test sp1a[:Elevation] == deg2rad(50.0)
    @test sp1a[:TakeOffAngle] == deg2rad(50.0)
    @test sp1a[:XRFExcitationPathLength] == 0.1 * 10.0
    @test sp1a[:XRFDetectionPathLength] == 0.1 * 20.0
    @test sp1a[:DetectorSolidAngle] == 0.0065
    @test sp1a[:ChamberPressure] == 20.0
    @test sp1a[:ChamberAtmosphere] == "Air"
    @test sp1a[:XRFSampleTilt] == 0.0
    f1 = sp1a[:XRFTubeWindow]
    @test isequal(f1.material, pure(n"Be"))
    @test isequal(f1.thickness, 200.0e-4)
end

@testset "SEM/EDS" begin
    fn2 = joinpath(@__DIR__, "Other", "Albite_ChMxd_20kV0p5nA130kHzTC_5p9kHzOCR_4ks.spx")
    sp2a = loadspectrum(fn2)
    #sp2b = FileIO.load(fn2)
    @test sp2a[:RealTime] == 4028.702
    @test sp2a[:LiveTime] == 4000.005
    @test sp2a[:BrukerThroughput] == 130000
    @test sp2a[:DetectorModel] == "Bruker XFlash 6|30"
    @test sp2a[:DetectorSerialNumber] == "13773"
    @test sp2a[:DetectorThickness] == 0.1 * 0.45
    @test sp2a[:DeadLayerThickness] == 0.1 * 0.029
    @test sp2a[:Elevation] == deg2rad(35.0)
    win = sp2a[:Window]
    @test length(win) == 5
    @test isapprox(win[1], Film(pure(n"B"), 1.3E1 * 1.0e-7))
    @test isapprox(win[2], Film(pure(n"C"), 1.45E1 * 1.0e-7))
    @test isapprox(win[3], Film(pure(n"N"), 4.5 * 1.0e-7))
    @test isapprox(win[4], Film(pure(n"O"), 8.5 * 1.0e-7))
    @test isapprox(win[5], Film(pure(n"Al"), 3.5 * 1.0e-7))
    @test sp2a[:BeamEnergy] == 20.0e3

    @test n"C" in sp2a[:Elements]
    @test n"O" in sp2a[:Elements]
    @test n"Na" in sp2a[:Elements]
    @test n"Al" in sp2a[:Elements]
    @test n"Si" in sp2a[:Elements]
    @test sp2a[:TakeOffAngle] == deg2rad(35.0)
    @test sp2a[:WorkingDistance] == 1.1
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
export Transformation, LinearTransformation, NonlinearTransformation

abstract type Transformation{T,D} <: AbstractSequence{T,D} end

abstract type LinearTransformation{T,D} <: Transformation{T,D} end

abstract type NonlinearTransformation{T,D} <: Transformation{T,D} end

for trans in (:Cumprod,:Cumsum)
    @eval begin
        export $trans
        mutable struct $trans{T,D} <: LinearTransformation{T,D}
            sequence::Vector{T}
            generator::D
        end
        $trans(generator) = $trans(eltype(generator)[],generator)
    end
end

for trans in (:Difference,:Shift)
    @eval begin
        export $trans, order
        mutable struct $trans{T,D,O} <: LinearTransformation{T,D}
            sequence::Vector{T}
            generator::D
        end
        $trans(sequence,generator,O::Int) = (@assert O ≥ 0; $trans{eltype(sequence),typeof(generator),O}(sequence,generator))
        $trans(generator,O::Int) = $trans(eltype(generator)[],generator,O)
        $trans(generator) = $trans(generator,1)
        order(::$trans{T,D,O}) where {T,D,O} = O
        order(::Type{$trans{T,D,O}}) where {T,D,O} = O
        *(T1::Type{$trans},T2::$trans) = $trans(sequence(T2),generator(T2),order(T2)+1)
    end
end

export Δ, E
Δ{T,D} = Difference{T,D}
E{T,D} = Shift{T,D}

for (op,trans) in ((:(Base.cumprod),:Cumprod),(:(Base.cumsum),:Cumsum),(:(Base.diff),:Δ))
    @eval begin
        $op(S::AbstractSequence) = $trans($op(sequence(S)),S)
    end
end

for trans in (:Aitken,:Euler)
    @eval begin
        export $trans
        mutable struct $trans{T,D} <: NonlinearTransformation{T,D}
            sequence::Vector{T}
            generator::D
        end
        $trans(generator) = (D = eltype(generator); T = typeof(one(D)/one(D));$trans(T[],generator))
    end
end

for trans in (:ϵ,:ρ)
    @eval begin
        export $trans
        mutable struct $trans{T,D} <: NonlinearTransformation{T,D}
            sequence::Vector{T}
            generator::D
            auxiliary::Vector{T}
        end
        $trans(generator) = (D = eltype(generator); T = typeof(one(D)/one(D));$trans(T[],generator,T[]))
        auxiliary(T::$trans) = T.auxiliary
    end
end

Shanks{T,D} = ϵ{T,D}

for trans in (:Drummond,:Levin,:Weniger)
    @eval begin
        export $trans
        mutable struct $trans{T,D,Symbol} <: NonlinearTransformation{T,D}
            sequence::Vector{T}
            generator::D
            numerator::Vector{T}
            denominator::Vector{T}
        end
        $trans(generator,ω::Symbol) = (D = eltype(generator);T = typeof(one(D)/one(D)); $trans{T,typeof(generator),ω}(T[],generator, T[], T[]))
        numerator(T::$trans) = T.numerator
        denominator(T::$trans) = T.denominator
    end
end

for trans in (:Drummond,:Levin)
    @eval begin
        remainder(S::$trans{T,D,:u},k::Int) where {T,D} = (s = generator(S); k > 1 ? (k+1)*(s[k]-s[k-1]) : 2s[1])
        remainder(S::$trans{T,D,:t},k::Int) where {T,D} = (s = generator(S); k > 1 ? s[k]-s[k-1] : s[1])
        remainder(S::$trans{T,D,:d},k::Int) where {T,D} = (s = generator(S); s[k+1]-s[k])
        remainder(S::$trans{T,D,:v},k::Int) where {T,D} = (s = generator(S); k > 1 ? (s[k]-s[k-1])*(s[k+1]-s[k])/(2s[k]-s[k-1]-s[k+1]) : s[1]*(s[2]-s[1])/(2s[1]-s[2]))
    end
end

remainder(S::Weniger{T,D,:y},k::Int) where {T,D} = (s = generator(S); k > 1 ? (k+1)*(s[k]-s[k-1]) : 2s[1])
remainder(S::Weniger{T,D,:τ},k::Int) where {T,D} = (s = generator(S); k > 1 ? s[k]-s[k-1] : s[1])
remainder(S::Weniger{T,D,:δ},k::Int) where {T,D} = (s = generator(S); s[k+1]-s[k])
remainder(S::Weniger{T,D,:φ},k::Int) where {T,D} = (s = generator(S); k > 1 ? (s[k]-s[k-1])*(s[k+1]-s[k])/(2s[k]-s[k-1]-s[k+1]) : s[1]*(s[2]-s[1])/(2s[1]-s[2]))

Drummond(generator) = Drummond(generator,:u)
Levin(generator) = Levin(generator,:u)
Weniger(generator) = Weniger(generator,:y)

max_stars_repo_name: MikaelSlevinsky/SequenceTransformations

max_issues_repo_name: MikaelSlevinsky/SequenceTransformations

max_forks_repo_name: MikaelSlevinsky/SequenceTransformations


text:
module Day1

using DelimitedFiles

count(vals) = sum(diff(vals) .> 0)

windowed_count(vals; k=3) = sum(diff(window_sums(vals, k=k)) .> 0)

test = [199
        200
        208
        210
        200
        207
        240
        269
        260
        263];

@assert count(test) == 7
@assert windowed_count(test) == 7

input = readdlm("resources/day1_input", '\t', Int, '\n') |> vec;

display(count(input))
display(windowed_count(input))

end


text:
@inline telemetry_alloc(::NoTelemetry, args...) = nothing

#####
##### Telemetry
#####

@enum TelemetryAction::UInt begin
    # Allocation
    AllocFast
    AllocSlow

    # Dellocations
    DeallocFast
    DeallocSlow

    # Transitions
    ReadTransition
    WriteTransition
    ReleaseTransition
    TransitionUnknown

    # Movement
    PrimaryInLocal
    PrimaryInRemote

    # Cleanup
    UnsafeFreed
    GarbageCollected
end

struct TelemetryRecord
    accesstime::UInt
    # Index into a global backtrace buffer.
    backtrace_key::Int
    action::TelemetryAction
end

struct AllocationRecord
    id::UInt
    size::Int
    basepointer::UInt
end

const BacktraceType = Vector{Union{Ptr{Nothing},Base.InterpreterIP}}

struct Telemetry
    # Map allocated IDs to their sequence of actions.
    logs::Dict{UInt, Vector{TelemetryRecord}}
    objects::Dict{UInt, AllocationRecord}
    backtraces::Dict{BacktraceType,Int}
    lock::Base.Threads.SpinLock
end

function Telemetry()
    return Telemetry(
        Dict{UInt, Vector{TelemetryRecord}}(),
        Dict{UInt, AllocationRecord}(),
        Dict{BacktraceType,Int}(),
        Base.Threads.SpinLock(),
    )
end

#####
##### Accessors
#####

ids(telemetry::Telemetry) = sort(collect(keys(telemetry.logs)))
getlog(telemetry::Telemetry, id) = telemetry.logs[id]

function getkey(telemetry::Telemetry)
    bt = backtrace()
    backtraces = telemetry.backtraces
    return get!(backtraces, bt, length(backtraces))
end

function Base.stacktrace(telemetry::Telemetry, record::TelemetryRecord)
    key = record.backtrace_key
    if key > 0
        return stacktrace(telemetry.backtraces[key])
    end
    return nothing
end

function Base.stacktrace(telemetry::Telemetry, key::Integer)
    bt = get(telemetry.backtraces, key, nothing)
    if bt !== nothing
        return stacktrace(bt)
    end
    return nothing
end

#####
##### Hooks
#####

function telemetry_alloc(telemetry::Telemetry, bytes, id, ptr, location)
    # Get the current time.
    now = time_ns()
    action = (location == Local) ? AllocFast : AllocSlow
    # Find where we are on the stack.
    # If this is a new location, mark it as such.
    allocation_record = AllocationRecord(id, bytes, UInt(ptr))
    @spinlock telemetry.lock begin
        telemetry.objects[id] = allocation_record
        key = getkey(telemetry)
        v = get!(() -> TelemetryRecord[], telemetry.logs, id)
        push!(v, TelemetryRecord(now, key, action))
    end
    return nothing
end

function telemetry_dealloc(telemetry::Telemetry, id, location)
    now = time_ns()
    action = (location == Local) ? DeallocFast : DeallocSlow
    @spinlock telemetry.lock begin
        key = getkey(telemetry)
        v = telemetry.logs[id]
        push!(v, TelemetryRecord(now, key, action))
    end
end

function telemetry_unsafefree(telemetry::Telemetry, id)
    now = time_ns()
    @spinlock telemetry.lock begin
        key = getkey(telemetry)
        v = telemetry.logs[id]
        push!(v, TelemetryRecord(now, key, UnsafeFreed))
    end
end

function telemetry_gc(telemetry::Telemetry, id)
    now = time_ns()
    @spinlock telemetry.lock begin
        log = telemetry.logs[id]
        push!(log, TelemetryRecord(now, 0, GarbageCollected))
    end
    return nothing
end

function telemetry_change(telemetry::Telemetry, id, to)
    # Determine the state change
    now = time_ns()
    state = TransitionUnknown
    if to == :ReadOnly
        state = ReadTransition
    elseif to == :ReadWrite
        state = WriteTransition
    elseif to == :NotBusy
        state = ReleaseTransition
    end

    # Find where we are on the stack.
    # If this is a new location, mark it as such.
    @spinlock telemetry.lock begin
        key = getkey(telemetry)
        log = telemetry.logs[id]
        push!(log, TelemetryRecord(now, key, state))
    end
    return nothing
end

function telemetry_primary(telemetry::Telemetry, id, location)
    action = location == Local ? PrimaryInLocal : PrimaryInRemote
    record = TelemetryRecord(time_ns(), 0, action)
    @spinlock telemetry.lock begin
        push!(telemetry.logs[id], record)
    end
    return nothing
end

#####
##### Analysis Passes
#####

# Find the backtrace key for all allocation sites.
function allocation_ids(telemetry::Telemetry)
    # Keep track of both the keys as well as the earliest access time so we can
    # reconstruct the call graph in some kind of ordered manner.
    #
    # Though to be honest, this will probably be equivalent to just straight up sorting
    # the keys ...
    keys = Dict{Int,UInt}()
    for log in values(telemetry.logs)
        record = first(log)
        @assert in(record.action, (AllocFast, AllocSlow))
        key = record.backtrace_key
        accesstime = record.accesstime

        keys[key] = min(accesstime, get!(keys, key, accesstime))
    end
    return sort(collect(Base.keys(keys)); by = x -> keys[x])
end

function stacktraces_for(telemetry::Telemetry, keys)
    return stacktrace.(Ref(telemetry), keys)
end

function estimate_lifetime(library::Vector{TelemetryRecord})
    # If this hasn't been GC'd yet, than return 0 as a sentinel value.
    last(library).action != GarbageCollected && return 0

    # Otherwise, the lifetime will be the difference between the first and last access
    # times. We need to be a bit careful if there wasn't a state transition for some
    # reason. In this case, the length of the record should be 2.
    #
    # The best we can do is say that the object was alive from when it was allocated
    # to when it was garbage collected.
    if length(library) == 2
        return last(library).accesstime - first(library).accesstime
    else
        return library[end-1].accesstime - first(library).accesstime
    end
end

function create_timeline(actions::Dict{K,<:Any}, objects) where {K}
    timeline = Vector{NamedTuple{(:accesstime, :id, :size, :action),Tuple{Int,Int,Int,Symbol}}}()
    for (id_str, logs) in actions
        sz = objects[id_str].size
        id = id_str
        for log in logs
            nt = (
                accesstime = log.accesstime,
                id = id,
                size = sz,
                action = log.action
            )
            push!(timeline, nt)
        end
    end
    sort!(timeline; by = x -> x.accesstime)
    return timeline
end

function heap_usage(timeline::AbstractVector{<:NamedTuple}; name = "Fast")
    alloc_str = Symbol("Alloc$name")
    dealloc_str = Symbol("Dealloc$name")

    usage = [(time = 0, utilization = 0)]
    for nt in timeline
        if (nt.action == alloc_str)
            utilization = usage[end].utilization + nt.size
            push!(usage, (time = nt.accesstime, utilization))
        elseif (nt.action == dealloc_str)
            utilization = usage[end].utilization - nt.size
            push!(usage, (time = nt.accesstime, utilization))
        end
    end
    popfirst!(usage)
    return usage
end

#####
##### Save Telemetry
#####

const GLOBAL_FILES = [
    "./boot.jl",
    "./client.jl",
    "./essentials.jl",
    "REPL/src/REPL.jl",
    "REPL[",
]

function prettify(stack::Vector{Base.StackTraces.StackFrame})
    # Remove all CachedArray files from the front of the trace.
    while true
        frame = first(stack)
        if occursin("CachedArrays/src", string(frame.file))
            popfirst!(stack)
        else
            break
        end
    end

    # Remove all REPL related stuff from the end of the trace.
    while any(x -> occursin(x, string(last(stack).file)), GLOBAL_FILES)
        pop!(stack)
    end
    return stack
end

# Custom JSON serialization for Telemetry
struct TelemetrySerialization <: JSON.CommonSerialization end
struct NoBacktraces{T}
    telemetry::T
end

function JSON.lower(x::Telemetry)
    # Reorder backtraces
    reverse_backtrace = Dict(v => k for (k,v) in x.backtraces)
    backtraces = [get(reverse_backtrace, i, "null") for i in Base.OneTo(length(reverse_backtrace))]

    return Dict(
        :logs => x.logs,
        :objects => x.objects,
        :backtraces => backtraces,
    )
end

function JSON.lower(x::NoBacktraces)
    telemetry = x.telemetry
    return Dict(
        :logs => telemetry.logs,
        :objects => telemetry.objects,
    )
end

const SC = JSON.StructuralContext

function JSON.show_json(io::SC, ::TelemetrySerialization, x::BacktraceType)
    # Step 1: Convert from a vector of pointers to a full on stacktrace.
    trace = prettify(stacktrace(x))
    return JSON.show_json(io, TelemetrySerialization(), trace)
end

function JSON.show_json(io::SC, ::TelemetrySerialization, x::Base.StackTraces.StackFrame)
    # Don't worry about macro expansions.
    if x.func == Symbol("macro expansion")
        return JSON.show_json(io, JSON.StandardSerialization(), "omitted")
    end

    if startswith(String(x.func), "top-level") || startswith(String(x.file), "REPL")
        return JSON.show_json(io, JSON.StandardSerialization(), "omitted")
    end

    dict = Dict(
        :func => x.func,
        :file => x.file,
        :line => x.line,
        :linfo => x.linfo
    )
    return JSON.show_json(io, TelemetrySerialization(), dict)
end

# Make type printing a little less heinous
tostring(x) = string(x)
function tostring(x::DataType)
    if isempty(x.parameters)
        return string(x)
    else
        inner = join(tostring.(x.parameters), ", ")
        return "$(x.name.module).$(x.name.name){$inner}"
    end
end

function tostring(::Type{NamedTuple{names,T}}) where {names,T}
    return "NamedTuple{$names, $(tostring(T))}"
end

tostring(x::Core.TypeofBottom) = "Union{}"
tostring(::Type{<:CachedArray{T,N}}) where {T,N} = "CachedArray{$T,$N}"
tostring(::Type{<:CacheManager}) = "CacheManager"

# Top level save functions
function JSON.show_json(io::SC, ::TelemetrySerialization, x::Core.MethodInstance)
    return JSON.show_json(io, JSON.StandardSerialization(), tostring.(x.specTypes.parameters))
end

function save_trace(file::AbstractString, x::Union{Telemetry, NoBacktraces})
    open(io -> save_trace(io, x), file; write = true)
end

function save_trace(io::IO, x::Union{Telemetry, NoBacktraces})
    return JSON.show_json(io, TelemetrySerialization(), x; indent = 4)
end

# #####
# ##### For mistakes
# #####
#
# just_comma(str) = match(r"^\s*,\s*$", str) !== nothing
# just_space(str) = all(isspace, str)
#
# function clean_json(dst::AbstractString, src::AbstractString)
#     ispath(dst) && rm(dst)
#     open(dst; write = true, create = true) do io
#         # Load up things into an IOBuffer
#         buf = IOBuffer()
#         last = ""
#         for ln in eachline(src)
#             last = clean_json(buf, ln, last)
#         end
#         println(buf, last)
#         seekstart(buf)
#
#         # Stream out remainders
#         last = ""
#         for ln in eachline(buf)
#             if match(r"^\s*{\s*", ln) !== nothing && endswith(last, "}")
#                 last = "$last,"
#             end
#             println(io, last)
#             last = ln
#         end
#         println(io, last)
#     end
#     return nothing
# end
#
# function clean_json(io::IO, str, last)
#     if just_comma(str)
#         last = replace(last, "," => "")
#     end
#     isempty(last) || println(io, last)
#
#     if just_comma(str) || just_space(str) || isempty(str)
#         return ""
#     end
#     return str
# end


text:
@doc """
	sparsegrid( D::Int, order::Int, f::Function=gausshermite; sym::Bool=true )

Computation of sparse grid nodes and the associated weights

- `D` : Dimension of integrant
- `k` : Order of quadrature rule
- `f` : Function generating 1D nodes and weights -- in that order -- for an integer input
- `sym` : Boolean variable determining if the nodes should be symmetrized

If the nodes are supposed to be symmetric (as those in the Gauss-Hermite rule),
they should be so in order to correctly identify multiply occuring nodes in the
union of sparse sets
"""->
function sparsegrid( D::Integer, order::Integer, f::Function=gausshermite; sym::Bool=true )
	# Final nodes and weights in D dimensions
	nodes = Array{Float64}(D, 0)
	weights = Array{Float64}(0)

	# Compute univariate nodes and weights
	nodes1D = Vector{Vector{Float64}}(order)
	weights1D = similar(nodes1D)

	for k in 1:order
		nodes1D[k], weights1D[k] = f( k )

		if sym
			symmetrize!( nodes1D[k] )
		end
	end

	nodes, weights = sparsegrid( D, nodes1D, weights1D )

	return nodes, weights
end

# Computation of sparse grid nodes and the associated weights
# from collection of one-dimensional nodes and weights

function sparsegrid( D::Integer, nodes1D::Vector{Vector{Float64}}, weights1D::Vector{Vector{Float64}} )
	order = length( nodes1D )

	# Final nodes and weights in D dimensions
	nodes = Array{Float64}(D, 0)
	weights = Array{Float64}(0)

	mink = max(0, order-D)
	maxk = order - 1

	# Temporary nodes and weights for 1 dimension
	N = Vector{Vector{Float64}}(D)
	W = similar(N)

	for k in mink:maxk
		alpha = listNdq(D, D+k)
		nalpha = size(alpha, 2)

		for n in 1:nalpha
			# The nodes and weights for this alpha mixture
			for d in 1:D
				N[d] = nodes1D[ alpha[d,n] ]
				W[d] = weights1D[ alpha[d,n] ]
			end

			# Compute all the possible combinations of D-dimensional nodes
			combN = combvec(N)

			# Compute the associated weights
			cw = combvec(W)
			combW = (-1)^(maxk-k) * binomial(D-1, D+k-order) * prod(cw, 1)

			# Save nodes and weights
			nodes = hcat(nodes, combN)
			weights = vcat(weights, vec(combW))
		end
	end

	# Remove redundant nodes
	unodes, uweights = uniquenodes(nodes, weights)

	return unodes, uweights
end

# ------------------------------------------------------------
# To correctly reduce "overlapping" nodes the middle node in an
# uneven number must be exactly zero

function symmetrize!( nodes::Vector{Float64} )
	N = length(nodes)

	if isodd(N)
		midpoint = div( N-1, 2 )
		nodes[ midpoint+1 ] = 0.0
	else
		midpoint = div( N, 2 )
	end

	for n in 1:midpoint
		nodes[n] = -nodes[N+1-n]
	end
end

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

@doc """
	listNdq( D::Int, q::Int )

Find elements in the set `N^D_q = { i \in N^D : sum(i) = q }`.

The algorithm and the formula for computing the number of elements in this set is found in the thesis mentioned in the README
"""->
function listNdq( D::Integer, q::Integer )
	if q < D
		error("listNdq: q must be larger than D")
	end

	M = binomial(q-1, D-1)
	L = ones(Int, D, M)

	k = ones(Int, D)
	maxk = q - D + 1
	khat = maxk*ones(Int, D)

	p = 1
	count = 0

	while k[D] < maxk
		k[p] += 1

		if k[p] > khat[p]
			k[p] = 1
			p += 1
		else
			for i in 1:p-1
				khat[i] = khat[p] - k[p] + 1
			end

			k[1] = khat[1]
			p = 1

			count += 1
			L[:,count] = k
		end
	end

	return L
end

@doc """
	combvec( vecs::Array{Any} ) -> Matrix

Counterpart of Matlab's combvec:
Creates all combinations of vectors in `vecs`, an array of vectors.
"""->
function combvec{T}( vecs::Vector{Vector{T}} )
	# Construct all Cartesian combinations of elements in vec as tuples
	P = product( vecs... )

	D = length(vecs)
	N = length(P)::Int64
	y = Array{T}(D,N)

	n = 0
	for p in P
		y[:,n+=1] = [p...]
	end

	return y
end

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

# Copy of Base.sortcols that also returns the permutation indices
# TODO: Update when 0.5/0.6 is out
function sortcolsidx(A::AbstractMatrix; kws...)
    r = 1:size(A,1)
    cols = [ view(A,r,i) for i=1:size(A,2) ]
    p = sortperm(cols; kws..., order=Base.Order.Lexicographic)

    return A[:,p], p
end

@doc """
	uniquenodes(nodes::Matrix, weights::Vector)

Find unique nodes and sum the weights of the identical nodes
"""->
function uniquenodes(nodes::AbstractMatrix, weights::AbstractVector)
	# Sort nodes and weights according to nodes
	sortnodes, P = sortcolsidx(nodes)
	weights = weights[P]

	D, N = size(nodes)

	keep = [1]
	lastkeep = 1

	for n in 2:N
		if sortnodes[:,n] == sortnodes[:,n-1]
			weights[lastkeep] += weights[n]
		else
			lastkeep = n
			push!(keep, n)
		end
	end

	# Keep unique nodes
	unodes = sortnodes[:, keep]
	uweights = weights[keep]

	return unodes, uweights
end

# ------------------------------------------------------------
# Compute tensor product grid

@doc """
	tensorgrid( N::Vector, W::Vector, D::Int )

	Compute tensor grid of `N` nodes and corresponding weights `W` for `D` dimensions.
"""->
function tensorgrid( N::Vector, W::Vector, D::Integer )
	NN = repeat( Any[N], outer=[D; 1] )
	WW = repeat( Any[W], outer=[D; 1] )

	tensorN = combvec( NN )
	tensorW = vec(prod( combvec(WW), 1 ))

	return tensorN, tensorW
end

function tensorgrid( D::Integer, order::Integer, f::Function=gausshermite )
	N, W = f(order)

	tensorN, tensorW = tensorgrid(N, W, D)

	return tensorN, tensorW
end


text:
module YaDB

import Base: length, join #, show

export Table, Row, GRow, WhereClause, HavingClause, D_SF,
  id, length, insert, update, delete, coltype,
  select, where, limit, group_by, order_by, #, show
  U_IN, U_SN

const Row = Dict{Symbol, Any}
const GRow = Union{Row, Vector{Pair{Symbol, Any}}, Vector{Any}}  ## Generic Row
const WhereClause = Function
const HavingClause = Function
const S_N = Union{Symbol, Nothing}
# const PK = Pair{S_N, DataType}
const D_SF = Dict{Symbol, Function}
const U_IN = Union{Int, Nothing}
const U_SN = Union{String, Nothing}
const UDT = Union{Union, DataType, Type}

mutable struct Table
  columns::Vector{Symbol}
  types::Vector{UDT}
  _pkey::S_N
  rows::Vector{Row}

  function Table(col_types::Vector{Pair{Symbol, DT}};
               pkey::Pair=(:id => Int)) where DT <: UDT
    # vector of pairs keeps implicit order
    @assert length(col_types) ≥ 1
    #
    cols = map(((k, _)=p) -> k, col_types)
    types = map(((_, v)=p) -> v, col_types)
    cols, types = check_pk(cols, types, pkey)
    new(cols, types, pkey.first, Vector{Row}[])
  end

  function Table(cols::Vector{Symbol}, types::Vector{DT};
               pkey::Pair=(:id => Int)) where DT <: UDT
    cols, types = check_pk(cols, types, pkey)
    new(cols, types, pkey.first, Vector{Row}[])
  end
end

## ================================
## Public API
## ================================

id(self::Table) = self._pkey

length(self::Table) = length(self.rows)

function insert(self::Table, row::Vector{Pair{Symbol, DT}}) where DT <: Any
  """
  row is something like [:name => "FooBar", :num_friends => 3, id(Users) => 3]
  """
  res = filter(p -> p.first == id(self), row)
  has_pkey = length(res) > 0
  ##
  ## NULL values permitted
  ##
  #length(row) ≠ length(self.types) && has_pkey &&
  #  throw(ArgumentError("Mismatch with expected number of columns"))
  ##
  check_value_dtype(self, row)
  ## find the pair <pkey, value> pkey is given by id(self)
  row_ids = filter(p -> p.first == id(self), row)
  if length(row_ids) == 0
    pk_id =  gen_pkey_value(self)
    row = [row..., id(self) => pk_id] # push!(row, id(self) => pk_id)
  else
    id_val = row_ids[1][2]  ## value assoc with this pkey
    row_already_inserted(self, id_val) && (return nothing)
  end
  push!(self.rows, Dict(row...))
end

function insert(self::Table, row::Vector{Any})
  """
  row is something like: [9, "Kate", 2 ] - just pure values
  """
  length(row) ≠ length(self.types) && id(self) ∈ map(p -> p[1], row) &&
    throw(ArgumentError("Mismatch with expected number of columns / $(length(row)) vs $(length(self.types))"))
  #
  check_value_dtype(self, row)
  ## row[1] is the value assoc. with pkey
  row_already_inserted(self, row[1]) && (return nothing)
  push!(self.rows, Dict(zip(self.columns, row)))
end

function insert(self::Table, row::Row)
  """
  row is a  Dict, like: insert(Users, Dict(:name => "Ayumi", :num_friends => 5))
  """
  length(row) ≠ length(self.types) && haskey(row, id(self)) &&
    throw(ArgumentError("Mismatch with expected number of columns"))
  check_value_dtype(self, row)
  if haskey(row, id(self))
    row_already_inserted(self, row[id(self)]) && (return nothing)
  else
    row[id(self)] = gen_pkey_value(self)
  end
  push!(self.rows, row)
end

insert(self::Table, rows::Vector{T}) where T <: GRow = insert.(Ref(self), rows)

function update(self::Table, updates::Row, pred::WhereClause=row -> true)
  ## 1 - Make sure key/columns and values are consistent with column types
  for (col, val) ∈ updates
    col ∉ self.columns && throw(ArgumentError("invalid column $(col)"))
    col_type = coltype(self, col)

    !(typeof(val) <: col_type) && val !== nothing &&
      throw(ArgumentError("Expected $(col_type), got $(typeof(val))"))
  end

  ## 2 - OK, update
  rows_to_update = filter(((ix, r)=t) -> pred(r), collect(enumerate(self.rows)))
  for (ix, row) ∈ rows_to_update
    self.rows[ix] = Row(row..., updates...)
  end
  nothing
end

function delete(self::Table, pred::WhereClause=row -> true)|
  self.rows = [
    row for row ∈ self.rows if !pred(row)
  ]
  nothing
end

function select(self::Table;
    keep_cols=Vector{Symbol}[], add_cols=Dict{Symbol, Function}())::Table
  ##
  length(keep_cols) == 0 && (keep_cols = self.columns)

  ## New column names and types
  new_cols = [keep_cols..., collect(keys(add_cols))...]
  keep_types = [coltype(self, col) for col in keep_cols]

  ## collect the rows for result table
  n_rows = Vector{Any}[]
  add_types = Any[]
  for (ix, row) ∈ enumerate(self.rows)
    n_row = Any[row[col] for col ∈ keep_cols]
    ## as we process the first row, we can get the return type...
    ## ...of each function defined in add_cols
    for (_col_name, fn) ∈ add_cols
      r = fn(row)
      ix == 1 && push!(add_types, typeof(r))
      push!(n_row, r)
    end
    push!(n_rows, n_row)
  end

  ## Create result table
  new_types = [keep_types..., add_types...]
  @assert(length(new_cols) == length(new_types),
    "length(new_cols) $(length(new_cols)) == length(new_types)")

  n_table = create_res_table(self, new_cols, new_types;
                             pred=s -> id(s) ∈ keep_cols)
  insert(n_table, n_rows)
  n_table
end

function where(self::Table, pred::WhereClause=row -> true)::Table
  """
  Only the rows that satisfy pred are returned
  """
  n_table = create_res_table(self, self.columns, self.types;
                             pred=s -> id(s) !== nothing)

  n_rows = Vector{Any}[]
  for row ∈ self.rows
    pred(row) && (push!(n_rows, [row[col] for col ∈ n_table.columns]))
  end
  insert(n_table, n_rows)
  n_table
end

function limit(self::Table, num_rows::Int=5)::Table
  """
  Only the first num_rows are returned
  """
  @assert 1 ≤ num_rows ≤ length(self.rows) "1 ≤ $(num_rows) ≤ $(length(self.rows))"

  n_table = create_res_table(self, self.columns, self.types;
                             pred=s -> id(s) !== nothing)

  ## NOTE: mark vector as GRow type
  # rows = GRow[Any[v] for r ∈ self.rows[1:num_rows] for (_, v) ∈ r]
  rows = Vector{Any}[]
  for row ∈ self.rows[1:num_rows]
    push!(rows, Any[v for (_, v) ∈ row])
  end
  insert(n_table, rows)
  n_table
end

function group_by(self::Table; group_by_cols::Vector{Symbol}, agg::D_SF,
  having=HavingClause=row_gp -> true)::Table
  grouped_rows = Dict{}()

  ## 1 - Populate groups
  for row ∈ self.rows
    key = Tuple(row[col] for col ∈ group_by_cols)
    ary = get(grouped_rows, key, Row[])
    push!(ary, row)
    grouped_rows[key] = ary
  end

  ## 2 - populate rows
  n_rows = Vector{Any}[]
  agg_types = Any[]
  for (ix, (key, rows)) ∈ enumerate(grouped_rows)
    agg_row = Any[]
    if having(rows)
      agg_row = Any[agg_row..., key...]  ## Keep Any[]
      for (_, agg_fn) ∈ agg
        r = agg_fn(rows)
        ix == 1 && push!(agg_types, typeof(r))
        push!(agg_row, r)
      end
    end
    push!(n_rows, agg_row)
  end

  ## 3 - res. table consists of group_by columns and aggregates
  new_cols = [group_by_cols..., collect(keys(agg))...] |> unique
  gp_by_types = [coltype(self, col) for col ∈ group_by_cols]
  new_types = [gp_by_types..., agg_types...]

  n_table = create_res_table(self, new_cols, new_types;
                             pred=s -> id(s) ∈ new_cols)
  insert(n_table, n_rows)
  n_table
end

function order_by(self::Table, order::Function)::Table
  n_table = select(self)
  sort!(n_table.rows, by=order)
  n_table
end

function join(self::Table, otable::Table; left_join=false)::Table
  join_on_cols = [c for c ∈ self.columns if c ∈ otable.columns]
  add_cols = [c for c ∈ otable.columns if c ∉ join_on_cols && c != id(otable)]
  new_cols = [self.columns..., add_cols...]
  new_types = UDT[self.types..., coltype.(Ref(otable), add_cols)...]

  n_table = Table(new_cols, new_types; pkey=nothing => Nothing)

  n_rows = Vector{Any}[]
  for row ∈ self.rows
    is_join = orow -> all(orow[c] == row[c] for c ∈ join_on_cols)

    o_rows = where(otable, is_join).rows
    for o_row ∈ o_rows
      push!(n_rows, Any[[row[c] for c ∈ self.columns]...,
          [o_row[c] for c ∈ add_cols]...])
    end

    if left_join && length(o_rows) == 0
      push!(n_rows, Any[[row[c] for c ∈ self.columns]...,
                        [nothing for _ ∈ add_cols]...])
    end
  end
  insert(n_table, n_rows)
  n_table
end

# function show(io::IO, self::Table)
#   s = "num. records: $(length(self)):\n"
#   for rd ∈ self.rows
#     s₁ = []
#     for (k, v) ∈ rd
#       push!(s₁, "$(k): $(v)")
#     end
#     s = string(s, "  <", join(s₁, ", "), ">\n")
#   end
#   print(io, "$(s)")
# end

function coltype(self::Table, colname::Symbol)::UDT
  ix = findfirst(col -> col == colname, self.columns)
  ix === nothing && throw(ArgumentError("column $(colname) inexistent"))
  self.types[ix]
end

## ================================
## Internals
## ================================

function create_res_table(self::Table, new_cols, new_types;
                          pred::Function=_a -> true)
  ## keep pkey or not...
  pred(self) ?
    Table(new_cols, new_types; pkey=id(self) => coltype(self, id(self))) :
    Table(new_cols, new_types; pkey=nothing => Nothing)
end

function check_pk(cols, types, pkey) # ::Pair
  if pkey.first !== nothing
    @assert pkey[2] <: Integer "Expecting an Integer, got: $(pkey[2])"
  end

  if pkey.first !== nothing && pkey.first ∉ cols
    cols  = [pkey[1], cols...]
    types = [pkey[2], types...]
  end
  (cols, types)
end

function check_value_dtype(self::Table, values::Vector{DT}) where DT <: Any
                           # values::Vector{Any})
  dtypes = self.types
  for (v, dt) ∈ zip(values, dtypes)
    !(typeof(v) <: dt) && v !== nothing &&
      throw("Expected type for <$(v)>: $(dt), got: $(typeof(v))")
  end
end

function check_value_dtype(self::Table,
                           row::Union{Row, Vector{Pair{Symbol, DT}}}) where DT <: Any
                           # row::Union{Row, Vector{Pair{Symbol, Any}}, Vector{Pair{Symbol, Union{Nothing, Any}}}})
  (cols, types) = self.columns, self.types
  col_type = Dict(zip(cols, types))
  for (k, v) ∈ row
    (haskey(col_type, k) && typeof(v) <: col_type[k]) ||
      throw("Expected type/2: $(col_type[k]), got: $(typeof(v))")
  end
end

function gen_pkey_value(self::Table)
  """
  Assuming pkey is instance of a numeric type
  """
  length(self.rows) == 0 && (return 1)
  (map(r -> r[id(self)], self.rows) |> maximum)  + 1
end

function row_already_inserted(self::Table, id_val::Any)::Bool
  id(self) === nothing && return false  ## no pkey => ignore check
  #
  res = filter(r -> r[id(self)] == id_val, self.rows)
  length(res) > 0    ## row already inserted id > 0
end

end  # Module

max_stars_repo_path: Data Science From Scratch/src/YaDB.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Data Science From Scratch/src/YaDB.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Data Science From Scratch/src/YaDB.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
using Compat.Test, MathOptInterface

const MOI  = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

@MOIU.model ModelForCachingOptimizer (MOI.ZeroOne, MOI.Integer) (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval) () () (MOI.SingleVariable,) (MOI.ScalarAffineFunction,) () ()

@testset "Continuous linear problems" begin
    optimizer = MOIU.CachingOptimizer(ModelForCachingOptimizer{Float64}(), Cbc.Optimizer())
    config = MOIT.TestConfig(atol=1e-4, rtol=1e-4, duals = false, infeas_certificates = false)
    MOIT.contlineartest(optimizer, config, [
    "linear1",  ## asks for ConstraintPrimal
    "linear2",  ## asks for ConstraintPrimal
    "linear7",  ## uses vector of constraints
    "linear10", ## asks for ConstraintPrimal
    "linear13", ## asks for ConstraintPrimal
    "linear14", ## asks for ConstraintPrimal
    "linear11", ## asks for ConstraintPrimal
    "linear15"  ## uses vector of constraints
    ])
end

@testset "Integer linear tests" begin

    optimizer = MOIU.CachingOptimizer(ModelForCachingOptimizer{Float64}(), Cbc.Optimizer())
    config = MOIT.TestConfig(atol=1e-4, rtol=1e-4, duals = false, infeas_certificates = false)
    # int1 excluded because asks for ConstraintPrimal
    # int2 excluded because uses vector of constraints
    MOIT.intlineartest(optimizer, config, ["int1", "int2"])

end

@testset "ModelLike tests" begin
    optimizer = MOIU.CachingOptimizer(ModelForCachingOptimizer{Float64}(), Cbc.Optimizer())

    # MOIT.nametest(optimizer) ## uses names/strings
    # @testset "validtest" begin ## at some moment inside the test it asks for the number of constraints of a type that is not supported by my model and caching optimizer returns an error. I cannot return the number of constraints of a specific type
    #     MOIT.validtest(optimizer)
    # end
    # @testset "emptytest" begin ## vector of constraints
    #     MOIT.emptytest(optimizer)
    # end
    @testset "orderedindicestest" begin
        MOIT.orderedindicestest(optimizer)
    end
    # @testset "canaddconstrainttest" begin ## do not pass due to vector of variables
    #     MOIT.canaddconstrainttest(optimizer, Float64, Complex{Float64})
    # end
    # @testset "copytest" begin ## do not pass due to vector of variables
    #     MOIT.copytest(optimizer.optimizer, optimizer)
    # end
end

@testset "Unit Tests" begin
    config = MOIT.TestConfig()
    optimizer = MOIU.CachingOptimizer(ModelForCachingOptimizer{Float64}(), Cbc.Optimizer())

    MOIT.basic_constraint_tests(optimizer, config)

    MOIT.unittest(optimizer, config, [
        "solve_with_lowerbound", ## cannot get with strings
        "solve_blank_obj",
        "solve_qcp_edge_cases",
        "solve_qp_edge_cases",
        "solve_affine_interval", ## cannot get with strings
        "solve_affine_greaterthan", ## cannot get with strings
        "solve_with_upperbound",  ## cannot get with strings
        "solve_singlevariable_obj", ## cannot get with strings
        "solve_affine_equalto",  ## cannot get with strings
        "solve_affine_lessthan",  ## cannot get with strings
        "solve_constant_obj",  ## cannot get with strings
        "solve_affine_deletion_edge_cases", ## do not support vector of constraints
        ## TODO: fix new tests of objective edge cases
        "solve_duplicate_terms_obj",
        "solve_objbound_edge_cases"
    ])
end


text:
function isnorm_fac_elem(K::NfRel{nf_elem}, a::nf_elem)
  Ka, mKa, mkK = collapse_top_layer(K)
  Kas, KasToKa = simplify(Ka)
  Ka = Kas
  #mKa = KasToKa * mKa
  mkK = mkK * inv(KasToKa)

  ZKa = lll(maximal_order(Ka))
  C, mC = class_group(ZKa, use_aut = true)

  S = collect(keys(factor(mkK(a)*ZKa)))

  c = get_attribute(ZKa, :ClassGrpCtx)
  FB = c.FB.ideals::Vector{ideal_type(order_type(Ka))}
  i = length(FB)
  q, mq = quo(C, elem_type(C)[preimage(mC, I) for I = S], false)
  while length(q) > 1
    while FB[i] in S || iszero(mq(preimage(mC, FB[i])))
      i -= 1
    end
    push!(S, FB[i])
    q, mmq = quo(q, [mq(preimage(mC, FB[i]))], false)
    mq = mq*mmq
  end

  s = Set(ideal_type(order_type(AnticNumberField))[minimum(mkK, I) for I = S])
  #make S relative Galois closed:
  PS = IdealSet(ZKa)
  S = reduce(vcat, Vector{ideal_type(ZKa)}[collect(keys(factor(PS(mkK, p)))) for p = s], init = Vector{ideal_type(ZKa)}())

  local U::GrpAbFinGen

  if length(S) == 0
    U, mU = unit_group_fac_elem(ZKa)
  else
    U, mU = sunit_group_fac_elem(collect(S))
  end

  class_group(parent(a))

  local u::GrpAbFinGen

  if length(s) == 0
    u, mu = unit_group_fac_elem(maximal_order(parent(a)))
  else
    u, mu = sunit_group_fac_elem(collect(s))
  end
  No = hom(U, u, elem_type(u)[preimage(mu, norm(mkK, mU(g))) for g = gens(U)])
  aa = preimage(mu, FacElem(a))::GrpAbFinGenElem
  fl, so = haspreimage(No, aa)
  fl || return fl, FacElem(one(K))
  return true, FacElem(K, Dict{elem_type(K), fmpz}([image(KasToKa * mKa, k) => v for (k,v) = (mU(so)::FacElem{elem_type(Ka), typeof(Ka)})]))
end

function isnorm(K::NfRel{nf_elem}, a::nf_elem)
  fl, s = isnorm_fac_elem(K, a)
  return fl, evaluate(s)
end

function norm_equation(K::NfRel{nf_elem}, a::nf_elem)
  fl, s = isnorm(K, a)
  fl || error("no solution")
  return s
end


text:
export AbstractNodeParams,
    Node,
    Edge,
    NodeGraph,
    INPUT,
    OUTPUT

const INPUT = 0
const OUTPUT = 1

"""
    abstract type AbstractNodeParams end

Abstract type that should be implemented to create new nodes.
"""
abstract type AbstractNodeParams end

"""
    mutable struct Node

Represents a node in the graph. Each node has the following fields:
- `id::Int`: Unique identifier of the node
- `name::String`: Name of the node
- `x::Float64`, `y::Float64`: Position of the node
- `params::AbstractNodeParams`: Parameters of the node

Each node has a finite number of input and output connectors, which can be
used to join them to other nodes through [`Edge`](@ref)s. The number of
input and output connectors is found through the [`input`](@ref) and
[`output`](@ref) functions corresponding to the type-symbol of `params`.
"""
mutable struct Node{P<:AbstractNodeParams}
    id::Int
    name::String
    x::Float64
    y::Float64
    params::P
end

"""
    struct Edge

Represents an edge between two nodes in the graph. Each edge has the following
fields:
- `from::Int`: Source node ID
- `from_type::Int`: 0 if the edge is from an input connector, 1 if
  from an output connector
- `from_conn::Int`: Source connector index
- `to::Int`: Destination node ID
- `to_type::Int`: 0 if the edge is from an input connector, 1 if
  from an output connector
- `to_conn::Int`: Destination connector index
"""
struct Edge
    from::Int
    from_type::Int
    from_conn::Int
    to::Int
    to_type::Int
    to_conn::Int
end

"""
    mutable struct NodeGraph

Struct containing vectors of the [`Node`](@ref)s and [`Edge`](@ref)s of
the graph.
"""
mutable struct NodeGraph
    nodes::Vector{Node}
    edges::Vector{Edge}
end

NodeGraph() = NodeGraph(Node[], Edge[])

function NodeGraph(path::String)
    load(path, "ng")
end

StructTypes.StructType(::Type{NodeGraph}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Edge}) = StructTypes.Struct()
StructTypes.StructType(::Type{N}) where {N<:Node} = StructTypes.DictType()
StructTypes.keyvaluepairs(n::Node) = [
    :node_id => n.id,
    :node_name => n.name,
    :pos => (x = n.x, y = n.y),
    :inputs => inputs(typeof(n.params)),
    :outputs => outputs(typeof(n.params)),
]
StructTypes.StructType(::Type{T}) where {T<:AbstractNodeParams} = StructTypes.DictType()
StructTypes.keyvaluepairs(t::T) where {T<:AbstractNodeParams} = [
    :type => type_symbol(t),
    (p => to_editor_format(getproperty(t, p)) for p in propertynames(t))...
]


text:
using HamNNNNclassifier
using Base.Test

# write your own tests here
@test train() == println("training...")
@test printdatafiledescription("/Users/henryolders/hammingnn/other\ datasets/TestAdjustInputs8.tsv") == ""


text:
# Adapted from 'Ordinary Differential Equations' test/units_tests.jl
using Test
using MechGlueDiffEqBase
using MechanicalUnits: @import_expand, dimension, NoDims, ∙
using OrdinaryDiffEq
@import_expand(N, s, m, km, kg)

#using OrdinaryDiffEq, Unitfu
algs = [Euler(),Midpoint(),Heun(),Ralston(),RK4(),SSPRK104(),SSPRK22(),SSPRK33(),
        SSPRK43(),SSPRK432(),BS3(),BS5(),DP5(),DP8(),Feagin10(),Feagin12(),
        Feagin14(),TanYam7(),Tsit5(),TsitPap8(),Vern6(),Vern7(),Vern8(),Vern9()]

@testset "Unitful time with unitless state" begin
    u0 = 30.0
    tspan = (0.0, 10.0)s
    prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
    prob2 = ODEProblem((u,t,p)    -> (-0.2∙s⁻¹ * u[1]), u0, tspan)
    prob3 = ODEProblem((u,t,p)    -> [-0.2∙s⁻¹* u[1]], [u0],tspan)
    for prob in [prob1, prob2, prob3]
        @test solve(prob, Tsit5()).retcode === :Success
    end
end

@testset "Scalar units" begin
    f(y,p,t) = 0.5*y / 3.0s
    u0 = 1.0∙N
    prob = ODEProblem(f,u0,(0.0,1.0)s)

    for alg in algs
        @show alg
        @test solve(prob,alg,dt=1∙s/10).retcode === :Success
    end
end

@testset "2D units" begin
    f(dy,y,p,t) = (dy .= 0.5.*y ./ 3.0s)
    u0 = [1.0 2.0
          3.0 1.0]N
    prob = ODEProblem(f,u0,(0.0, 1.0)s)

    for alg in algs
        @show alg
        sol = solve(prob,alg,dt=1∙s/10)
    end
    sol = solve(prob,ExplicitRK())
end

@testset "Without ArrayPartition - broken" begin
    # coordinate: u = [position, momentum]
    # parameters: p = [mass, force constanst]
    function f_harmonic!(du,u,p,t)
        du[1] = u[2]/p[1]
        du[2] = -p[2]*u[1]
    end

    mass = 1.0∙kg
    k = 1.0∙N/m
    p = [mass, k]

    u0 = [1.0m, 0.0∙kg∙m/s] # initial values (position, momentum)
    tspan = (0.0∙s, 10.0∙s)
    prob = ODEProblem(f_harmonic!, u0, tspan, p)
    @test_broken solve(prob, Tsit5()).retcode == :Success
end