text:
@kwdef mutable struct SAddHandler <: AbstractHandler
    id = 0
end

@register(SAddHandler, +, begin
    [SBit, SBit] => [SBit]
    [SBit, Vector{SBit}] => [Vector{SBit}]
    [Vector{SBit}, SBit] => [Vector{SBit}]
    [Vector{SBit}, Vector{SBit}] => [Vector{SBit}]
    [SBit, Matrix{SBit}] => [Matrix{SBit}]
    [Matrix{SBit}, SBit] => [Matrix{SBit}]
    [Matrix{SBit}, Matrix{SBit}] => [Matrix{SBit}]
end)

function (handler::SAddHandler)(netlist::Netlist,
                                inputs::Vector{Variable},
                                outputs::Vector{Variable})
    # update netlist with inputs
    setsigned!(netlist, getname(inputs[1]), true)
    setsigned!(netlist, getname(inputs[2]), true)

    # compute output size
    lname, rname, outsize = handlebroadcast(inputs[1].name, inputs[2].name,
                                            getsize(netlist, getname(inputs[1])),
                                            getsize(netlist, getname(inputs[2])))

    # update netlist with output
    setsigned!(netlist, getname(outputs[1]), true)

    outstring = """
        $stdcomment
        // BEGIN add$(handler.id)
        stoch_add_mat #(
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) add$(handler.id)_pp (
                .CLK(CLK),
                .nRST(nRST),
                .A($(lname("_p"))),
                .B($(rname("_p"))),
                .Y($(outputs[1].name)_p)
            );
        stoch_add_mat #(
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) add$(handler.id)_mm (
                .CLK(CLK),
                .nRST(nRST),
                .A($(lname("_m"))),
                .B($(rname("_m"))),
                .Y($(outputs[1].name)_m)
            );
        // END add$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end


text:
##=================####==============##
# MZH Agents Interactions
##=================####==============##

"""
    get_next_state!(agent, params)
Finds `agent.new_state` according to its contacts and meetings
"""
function get_next_state!(agents, ag_id, params)

    agent = agents[ag_id]

    # IF AGENT IS SUCEPTIBLE AND IS *NOT* AT HOME...
    if agent.state == "S" && agent.at_home == false

        if agent.num_meets <= agent.degree_t
            meets = [agents[id] for id in sample(agent.contacts_t, agent.num_meets, replace=false)]
        else
            meets = [agents[id] for id in agent.contacts_t]
        end

        # println(agent.id, "|meets:", [x.id for x in meets])

        for c_agent in meets # LOOP THROUGH AGENTS MET THAT TIMESTEP
            coll_state = agent.state * c_agent.state
            # SUCEPTIBLE <- INFECTED
            if coll_state == "SI" && c_agent.at_home == false
                # println(coll_state)
                if rand() < params.attack_rate
                    agent.new_state = "I"
                    agent.infection_t = params.now_t
                    break
                end
            end
        end
    elseif agent.state == "I" # INFECTED AGENT TO RECOVER

        time_from_infection = params.now_t - agent.infection_t
        if time_from_infection == agent.recovery_t
            agent.new_state = "R"
            # println(agent.id, "|", time_from_infection, "|RECOVERED!")
        end

    end
end

##=================####==============##

ext: jl


text:
@add_problem affine function affine_negate_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable()
    p = minimize(-x, [x <= 0]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = atol rtol = rtol
        @test evaluate(-x) ≈ 0 atol = atol rtol = rtol
    end
end

@add_problem affine function affine_kron_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = ComplexVariable(3, 3)
    y = [1.0 2.0; 3.0 4.0]
    if test
        @test size(kron(x, y)) == (6, 6)
        @test size(kron(y, x)) == (6, 6)
    end
end

@add_problem affine function affine_multiply_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(1)
    p = minimize(2.0 * x, [x >= 2, x <= 4]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol = atol rtol = rtol
        @test (evaluate(2.0x))[1] ≈ 4 atol = atol rtol = rtol
    end

    x = Variable(2)
    A = 1.5 * eye(2)
    p = minimize([2 2] * x, [A * x >= [1.1; 1.1]]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2.93333 atol = atol rtol = rtol
        @test (evaluate([2 2] * x))[1] ≈ 2.93333 atol = atol rtol = rtol
        @test vec(evaluate(A * x)) ≈ [1.1; 1.1] atol = atol rtol = rtol
    end

    y = Variable(1)
    x = Variable(3)
    z = [1.0, 2.0, 3.0] * y
    k = -y * [1.0, 2.0, 3.0]
    c = [y <= 3.0, y >= 0.0, x >= ones(3), k <= x, x <= z]
    o = 3 * y
    p = Problem{T}(:minimize, o, c)
    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 3 atol = atol rtol = rtol
    end

    p = Problem{T}(:minimize, o, c...)
    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 3 atol = atol rtol = rtol
    end

    # Check #274
    x = ComplexVariable(2, 2)
    p = minimize(
        real([1.0im, 0.0]' * x * [1.0im, 0.0]),
        [x == [1.0 0.0; 0.0 1.0]];
        numeric_type = T,
    )

    handle_problem!(p)
    if test
        @test p.optval ≈ 1.0 atol = atol rtol = rtol
    end
end

@add_problem affine function affine_dot_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2)
    p = minimize(dot([2.0; 2.0], x), x >= [1.1; 1.1]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4.4 atol = atol rtol = rtol
        @test (evaluate(dot([2.0; 2.0], x)))[1] ≈ 4.4 atol = atol rtol = rtol
    end
end

@add_problem affine function affine_dot_atom_for_matrix_variables(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2, 2)
    p = minimize(dot(fill(2.0, (2, 2)), x), x >= 1.1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 8.8 atol = atol rtol = rtol
        @test (evaluate(dot(fill(2.0, (2, 2)), x)))[1] ≈ 8.8 atol = atol rtol =
            rtol
    end
end

@add_problem affine function affine_add_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(1)
    y = Variable(1)
    p = minimize(x + y, [x >= 3, y >= 2]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 5 atol = atol rtol = rtol
        @test evaluate(x + y) ≈ 5 atol = atol rtol = rtol
    end

    x = Variable(1)
    p = minimize(x, [eye(2) + x >= eye(2)]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = atol rtol = rtol
        @test evaluate(eye(2) + x) ≈ eye(2) atol = atol rtol = rtol
    end

    y = Variable()
    p = minimize(y - 5, y >= -1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ -6 atol = atol rtol = rtol
        @test evaluate(y - 5) ≈ -6 atol = atol rtol = rtol
    end
end

@add_problem affine function affine_transpose_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2)
    c = ones(2, 1)
    p = minimize(x' * c, x >= 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol = atol rtol = rtol
        @test (evaluate(x' * c))[1] ≈ 2 atol = atol rtol = rtol
    end

    X = Variable(2, 2)
    c = ones(2, 1)
    p = minimize(c' * X' * c, [X >= ones(2, 2)]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol = atol rtol = rtol
        @test (evaluate(c' * X' * c))[1] ≈ 4 atol = atol rtol = rtol
    end

    rows = 2
    cols = 3
    r = rand(rows, cols)
    r_2 = rand(cols, rows)
    x = Variable(rows, cols)
    c = ones(1, cols)
    d = ones(rows, 1)
    p = minimize(
        c * x' * d + d' * x * c' + (c * x''''' * d)',
        [x' >= r_2, x >= r, x''' >= r_2, x'' >= r];
        numeric_type = T,
    )
    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    s = sum(max.(r, r_2')) * 3
    if test
        @test p.optval ≈ s atol = atol rtol = rtol
        @test (evaluate(
            c * x' * d + d' * x * c' + (c * ((((x')')')')' * d)',
        ))[1] ≈ s atol = atol rtol = rtol
    end
end

@add_problem affine function affine_index_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2)
    p = minimize(x[1] + x[2], [x >= 1]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol = atol rtol = rtol
        @test (evaluate(x[1] + x[2]))[1] ≈ 2 atol = atol rtol = rtol
    end

    x = Variable(3)
    I = [true true false]
    p = minimize(sum(x[I]), [x >= 1]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol = atol rtol = rtol
        @test (evaluate(sum(x[I])))[1] ≈ 2 atol = atol rtol = rtol
    end

    rows = 6
    cols = 8
    n = 2
    X = Variable(rows, cols)
    A = randn(rows, cols)
    c = rand(1, n)
    p = minimize(c * X[1:n, 5:5+n-1]' * c', X >= A; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    s = c * A[1:n, 5:5+n-1]' * c'
    if test
        @test p.optval ≈ s[1] atol = atol rtol = rtol
        @test evaluate(c * (X[1:n, 5:(5+n)-1])' * c') ≈ s atol = atol rtol =
            rtol
    end
end

@add_problem affine function affine_sum_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2, 2)
    p = minimize(sum(x), x >= 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol = atol rtol = rtol
        @test evaluate(sum(x)) ≈ 4 atol = atol rtol = rtol
    end

    x = Variable(2, 2)
    p = minimize(sum(x) - 2 * x[1, 1], x >= 1, x[1, 1] <= 2; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol = atol rtol = rtol
        @test (evaluate(sum(x) - 2 * x[1, 1]))[1] ≈ 1 atol = atol rtol = rtol
    end

    x = Variable(10)
    a = rand(10, 1)
    p = maximize(sum(x[2:6]), x <= a; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ sum(a[2:6]) atol = atol rtol = rtol
        @test evaluate(sum(x[2:6])) ≈ sum(a[2:6]) atol = atol rtol = rtol
    end
end

@add_problem affine function affine_diag_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2, 2)
    p = minimize(sum(diag(x, 1)), x >= 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol = atol rtol = rtol
        @test evaluate(sum(diag(x, 1))) ≈ 1 atol = atol rtol = rtol
    end

    x = Variable(4, 4)
    p = minimize(sum(diag(x)), x >= 2; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 8 atol = atol rtol = rtol
        @test evaluate(sum(diag(x))) ≈ 8 atol = atol rtol = rtol
    end
end

@add_problem affine function affine_trace_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2, 2)
    p = minimize(tr(x), x >= 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol = atol rtol = rtol
        @test evaluate(tr(x)) ≈ 2 atol = atol rtol = rtol
    end
end

@add_problem affine function affine_dot_multiply_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(3)
    p = maximize(sum(dot(*)(x, [1, 2, 3])), x <= 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 6 atol = atol rtol = rtol
        @test evaluate(sum((dot(*))(x, [1, 2, 3]))) ≈ 6 atol = atol rtol = rtol
    end

    x = Variable(3, 3)
    p = maximize(sum(dot(*)(x, eye(3))), x <= 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 3 atol = atol rtol = rtol
        @test evaluate(sum((dot(*))(x, eye(3)))) ≈ 3 atol = atol rtol = rtol
    end

    x = Variable(5, 5)
    p = minimize(x[1, 1], dot(*)(3, x) >= 3; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol = atol rtol = rtol
        @test (evaluate(x[1, 1]))[1] ≈ 1 atol = atol rtol = rtol
    end

    x = Variable(3, 1)
    p = minimize(sum(dot(*)(ones(3, 3), x)), x >= 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 9 atol = atol rtol = rtol
        @test (evaluate(x[1, 1]))[1] ≈ 1 atol = atol rtol = rtol
    end

    x = Variable(1, 3)
    p = minimize(sum(dot(*)(ones(3, 3), x)), x >= 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 9 atol = atol rtol = rtol
        @test (evaluate(x[1, 1]))[1] ≈ 1 atol = atol rtol = rtol
    end

    x = Variable(1, 3, Positive())
    p = maximize(sum(dot(/)(x, [1 2 3])), x <= 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 11 / 6 atol = atol rtol = rtol
        @test evaluate(sum((dot(/))(x, [1 2 3]))) ≈ 11 / 6 atol = atol rtol =
            rtol
    end

    # Broadcast fusion works
    x = Variable(5, 5)
    a = 2.0 .* x .* ones(Int, 5)
    if test
        @test a isa Convex.DotMultiplyAtom
    end
end

@add_problem affine function affine_reshape_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    A = rand(2, 3)
    X = Variable(3, 2)
    c = rand()
    p = minimize(sum(reshape(X, 2, 3) + A), X >= c; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ sum(A .+ c) atol = atol rtol = rtol
        @test evaluate(sum(reshape(X, 2, 3) + A)) ≈ sum(A .+ c) atol = atol rtol =
            rtol
    end

    b = rand(6)
    p = minimize(sum(vec(X) + b), X >= c; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ sum(b .+ c) atol = atol rtol = rtol
        @test evaluate(sum(vec(X) + b)) ≈ sum(b .+ c) atol = atol rtol = rtol
    end

    x = Variable(4, 4)
    c = ones(16, 1)
    reshaped = reshape(x, 16, 1)
    a = collect(1:16)
    p = minimize(c' * reshaped, reshaped >= a; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    # TODO: why is accuracy lower here?
    if test
        @test p.optval ≈ 136 atol = 10atol atol = atol rtol = rtol
        @test (evaluate(c' * reshaped))[1] ≈ 136 atol = 10atol atol = atol rtol =
            rtol
    end
end

@add_problem affine function affine_hcat_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(4, 4)
    y = Variable(4, 6)
    p = maximize(
        sum(x) + sum([y fill(4.0, 4)]),
        [x y fill(2.0, (4, 2))] <= 2;
        numeric_type = T,
    )

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 96 atol = atol rtol = rtol
        @test evaluate(sum(x) + sum([y fill(4.0, 4)])) ≈ 96 atol = atol rtol =
            rtol
        @test evaluate([x y fill(2.0, (4, 2))]) ≈ fill(2.0, (4, 12)) atol = atol rtol =
            rtol
    end
end

@add_problem affine function affine_vcat_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(4, 4)
    y = Variable(4, 6)

    # TODO: fix dimension mismatch [y 4*eye(4); x -ones(4, 6)]
    p = maximize(
        sum(x) + sum([y 4*eye(4); x -ones(4, 6)]),
        [x; y'] <= 2;
        numeric_type = T,
    )

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    # TODO: why is accuracy lower here?
    if test
        @test p.optval ≈ 104 atol = 10atol atol = atol rtol = rtol
        @test evaluate(sum(x) + sum([y 4*eye(4); x -(ones(4, 6))])) ≈ 104 atol =
            10atol atol = atol rtol = rtol
        @test evaluate([x; y']) ≈ 2 * ones(10, 4) atol = atol rtol = rtol
    end
end

@add_problem affine function affine_single_hcat_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(4, 4)
    p = maximize(tr(hcat(x)), hcat(x) <= 2; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.optval ≈ 8 atol = atol rtol = rtol
    end
end

@add_problem affine function affine_single_vcat_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(4, 4)
    p = maximize(tr(vcat(x)), vcat(x) <= 2; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.optval ≈ 8 atol = atol rtol = rtol
    end
end

@add_problem affine function affine_Diagonal_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2, 2)
    if test
        @test_throws ArgumentError Diagonal(x)
    end

    x = Variable(4)
    p = minimize(sum(Diagonal(x)), x == [1, 2, 3, 4]; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 10 atol = atol rtol = rtol
        @test all(abs.(evaluate(Diagonal(x)) - Diagonal([1, 2, 3, 4])) .<= atol)
    end

    x = Variable(3)
    c = [1; 2; 3]
    p = minimize(c' * Diagonal(x) * c, x >= 1, sum(x) == 10; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 21 atol = atol rtol = rtol
    end

    x = Variable(3)
    p = minimize(sum(x), x >= 1, Diagonal(x)[1, 2] == 1; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.status != MOI.OPTIMAL
    end
end

@add_problem affine function affine_conv_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(3)
    h = [1, -1]
    p = minimize(sum(conv(h, x)) + sum(x), x >= 1, x <= 2; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 3 atol = atol rtol = rtol
        @test evaluate(sum(conv(h, x))) ≈ 0 atol = atol rtol = rtol
    end

    x = Variable(3)
    h = [1, -1]
    p = minimize(sum(conv(x, h)) + sum(x), x >= 1, x <= 2; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 3 atol = atol rtol = rtol
        @test evaluate(sum(conv(h, x))) ≈ 0 atol = atol rtol = rtol
    end
end

@add_problem affine function affine_satisfy_problems(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable()
    p = satisfy(x >= 0; numeric_type = T)

    add_constraints!(p, x >= 1)
    add_constraints!(p, [x >= -1, x <= 4])
    handle_problem!(p)
    if test
        @test p.status == MOI.OPTIMAL
    end

    p = satisfy([x >= 0, x >= 1, x <= 2]; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.status == MOI.OPTIMAL
    end

    p = maximize(1, [x >= 1, x <= 2]; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.status == MOI.OPTIMAL
    end

    constr = x >= 0
    constr += x >= 1
    constr += x <= 10
    constr2 = x >= 0
    constr2 += [x >= 2, x <= 3] + constr
    p = satisfy(constr; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.status == MOI.OPTIMAL
    end
end

@add_problem affine function affine_dualvalue(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable()
    p = minimize(x, x >= 0; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.constraints[1].dual ≈ 1 atol = atol rtol = rtol
    end

    x = Variable()
    p = maximize(x, x <= 0; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.constraints[1].dual ≈ 1 atol = atol rtol = rtol
    end

    x = Variable()
    p = minimize(x, x >= 0, x == 2; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.constraints[1].dual ≈ 0 atol = atol rtol = rtol
        @test abs.(p.constraints[2].dual) ≈ 1 atol = atol rtol = rtol
    end

    x = Variable(2)
    A = 1.5 * eye(2)
    p = minimize(dot([2.0; 2.0], x), [A * x >= [1.1; 1.1]]; numeric_type = T)

    handle_problem!(p)
    if test
        dual = [4 / 3; 4 / 3]
        @test all(abs.(p.constraints[1].dual - dual) .<= atol)
    end
end

@add_problem affine function affine_Partial_transpose(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    dims = [2, 3, 4]
    d = prod(dims)
    A = rand(ComplexF64, 2, 2)
    B = rand(ComplexF64, 3, 3)
    C = rand(ComplexF64, 4, 4)
    M = kron(A, B, C)
    Mt1 = kron(transpose(A), B, C)
    Mt2 = kron(A, transpose(B), C)
    Mt3 = kron(A, B, transpose(C))

    Rt1 = ComplexVariable(d, d)
    Rt2 = ComplexVariable(d, d)
    Rt3 = ComplexVariable(d, d)
    S = rand(ComplexF64, d, d)
    handle_problem!(
        satisfy(partialtranspose(Rt1, 1, dims) == S; numeric_type = T),
    )

    handle_problem!(
        satisfy(partialtranspose(Rt2, 2, dims) == S; numeric_type = T),
    )

    handle_problem!(
        satisfy(partialtranspose(Rt3, 3, dims) == S; numeric_type = T),
    )

    if test
        @test partialtranspose(M, 1, dims) ≈ Mt1 atol = atol rtol = rtol
        @test partialtranspose(M, 2, dims) ≈ Mt2 atol = atol rtol = rtol
        @test partialtranspose(M, 3, dims) ≈ Mt3 atol = atol rtol = rtol
        @test partialtranspose(S, 1, dims) ≈ evaluate(Rt1) atol = atol rtol =
            rtol
        @test partialtranspose(S, 2, dims) ≈ evaluate(Rt2) atol = atol rtol =
            rtol
        @test partialtranspose(S, 3, dims) ≈ evaluate(Rt3) atol = atol rtol =
            rtol
    end

    if test
        @test_throws ArgumentError partialtrace(rand(6, 6), 3, [2, 3])
        @test_throws ArgumentError partialtrace(rand(6, 6), 1, [2, 4])
        @test_throws ArgumentError partialtrace(rand(3, 4), 1, [2, 3])
    end
end

@add_problem affine function affine_permuteddims_matrix(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    #this function is used in the partial transpose
    for n in (2, 3, 4, 5)
        dims = ntuple(i -> rand(2:5), n)
        d = prod(dims)
        v = rand(d)
        p = randperm(n)
        out1 = vec(permutedims(reshape(v, dims), p))
        out2 = Convex.permutedims_matrix(dims, p) * v
        if test
            @test out1 ≈ out2 atol = atol rtol = rtol
        end
    end
end


text:
module TestRequest
import HTTP: Layer, request, Response

abstract type TestLayer{Next <: Layer} <: Layer{Next} end
abstract type LastLayer{Next <: Layer} <: Layer{Next} end
export TestLayer, LastLayer, request

function request(::Type{TestLayer{Next}}, io::IO, req, body; kw...)::Response where Next
		return request(Next, io, req, body; kw...)
end

const FLAG = Ref(false)
function request(::Type{LastLayer{Next}}, resp)::Response where Next
    FLAG[] = true
		return request(Next, resp)
end

end


text:
################################################################################
#
#   SQUARE OCTAGON LATTICE
#
################################################################################

# Implementation
# - implementation 1 (Standard)
# - labels <: Any
function getUnitcellShastrySutherland(
            unitcell_type  :: Type{U},
            implementation :: Val{1}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2}, U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # Bravais lattice vectors
        Vector{Float64}[
            Float64[2, 0],
            Float64[0, 2]
        ],
        # Sites
        S[
            newSite(S, Float64[0,0], getDefaultLabelN(LS,1)),
            newSite(S, Float64[0,1], getDefaultLabelN(LS,2)),
            newSite(S, Float64[1,0], getDefaultLabelN(LS,3)),
            newSite(S, Float64[1,1], getDefaultLabelN(LS,4))
        ],
        # Bonds
        B[
            newBond(B, 1, 2, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 3, 4, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 1, 4, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 1, 3, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 2, 4, getDefaultLabel(LB), ( 0, 0)),

            newBond(B, 2, 1, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 4, 3, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 4, 1, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 3, 1, getDefaultLabel(LB), ( 0, 0)),
            newBond(B, 4, 2, getDefaultLabel(LB), ( 0, 0)),

            newBond(B, 2, 1, getDefaultLabel(LB), ( 0, 1)),
            newBond(B, 1, 2, getDefaultLabel(LB), ( 0,-1)),

            newBond(B, 4, 3, getDefaultLabel(LB), ( 0, 1)),
            newBond(B, 3, 4, getDefaultLabel(LB), ( 0,-1)),

            newBond(B, 3, 1, getDefaultLabel(LB), ( 1, 0)),
            newBond(B, 1, 3, getDefaultLabel(LB), (-1, 0)),

            newBond(B, 4, 2, getDefaultLabel(LB), ( 1, 0)),
            newBond(B, 2, 4, getDefaultLabel(LB), (-1, 0)),

            newBond(B, 2, 3, getDefaultLabel(LB), (-1, 1)),
            newBond(B, 3, 2, getDefaultLabel(LB), ( 1,-1))
        ]
    )
end

# Implementation
# - implementation 4 (Kitaev)
# - labels <: Any
function getUnitcellShastrySutherland(
            unitcell_type  :: Type{U},
            implementation :: Val{4}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2}, U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # Bravais lattice vectors
        Vector{Float64}[
            Float64[2, 0],
            Float64[0, 2]
        ],
        # Sites
        S[
            newSite(S, Float64[0,0], getDefaultLabelN(LS,1)),
            newSite(S, Float64[0,1], getDefaultLabelN(LS,2)),
            newSite(S, Float64[1,0], getDefaultLabelN(LS,3)),
            newSite(S, Float64[1,1], getDefaultLabelN(LS,4))
        ],
        # Bonds
        B[
            newBond(B, 1, 2, getDefaultLabelY(LB), ( 0, 0)),
            newBond(B, 3, 4, getDefaultLabelY(LB), ( 0, 0)),
            newBond(B, 1, 4, getDefaultLabelZ(LB), ( 0, 0)),
            newBond(B, 1, 3, getDefaultLabelY(LB), ( 0, 0)),
            newBond(B, 2, 4, getDefaultLabelY(LB), ( 0, 0)),

            newBond(B, 2, 1, getDefaultLabelY(LB), ( 0, 0)),
            newBond(B, 4, 3, getDefaultLabelY(LB), ( 0, 0)),
            newBond(B, 4, 1, getDefaultLabelZ(LB), ( 0, 0)),
            newBond(B, 3, 1, getDefaultLabelY(LB), ( 0, 0)),
            newBond(B, 4, 2, getDefaultLabelY(LB), ( 0, 0)),

            newBond(B, 2, 1, getDefaultLabelX(LB), ( 0, 1)),
            newBond(B, 1, 2, getDefaultLabelX(LB), ( 0,-1)),

            newBond(B, 4, 3, getDefaultLabelX(LB), ( 0, 1)),
            newBond(B, 3, 4, getDefaultLabelX(LB), ( 0,-1)),

            newBond(B, 3, 1, getDefaultLabelX(LB), ( 1, 0)),
            newBond(B, 1, 3, getDefaultLabelX(LB), (-1, 0)),

            newBond(B, 4, 2, getDefaultLabelX(LB), ( 1, 0)),
            newBond(B, 2, 4, getDefaultLabelX(LB), (-1, 0)),

            newBond(B, 2, 3, getDefaultLabelZ(LB), (-1, 1)),
            newBond(B, 3, 2, getDefaultLabelZ(LB), ( 1,-1))
        ]
    )
end

# Implementation
# - implementation -4 (Kitaev with signs)
# - labels <: Any
function getUnitcellShastrySutherland(
            unitcell_type  :: Type{U},
            implementation :: Val{-4}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2}, U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # Bravais lattice vectors
        Vector{Float64}[
            Float64[2, 0],
            Float64[0, 2]
        ],
        # Sites
        S[
            newSite(S, Float64[0,0], getDefaultLabelN(LS,1)),
            newSite(S, Float64[0,1], getDefaultLabelN(LS,2)),
            newSite(S, Float64[1,0], getDefaultLabelN(LS,3)),
            newSite(S, Float64[1,1], getDefaultLabelN(LS,4))
        ],
        # Bonds
        B[
            newBond(B, 1, 2, getDefaultLabelY(LB), ( 0, 0)),
            newBond(B, 3, 4, getDefaultLabelY(LB), ( 0, 0)),
            newBond(B, 1, 4, getDefaultLabelZ(LB), ( 0, 0)),
            newBond(B, 1, 3, getDefaultLabelY(LB), ( 0, 0)),
            newBond(B, 2, 4, getDefaultLabelY(LB), ( 0, 0)),

            newBond(B, 2, 1, flipLabel(getDefaultLabelY(LB)), ( 0, 0)),
            newBond(B, 4, 3, flipLabel(getDefaultLabelY(LB)), ( 0, 0)),
            newBond(B, 4, 1, flipLabel(getDefaultLabelZ(LB)), ( 0, 0)),
            newBond(B, 3, 1, flipLabel(getDefaultLabelY(LB)), ( 0, 0)),
            newBond(B, 4, 2, flipLabel(getDefaultLabelY(LB)), ( 0, 0)),

            newBond(B, 2, 1, getDefaultLabelX(LB), ( 0, 1)),
            newBond(B, 1, 2, flipLabel(getDefaultLabelX(LB)), ( 0,-1)),

            newBond(B, 4, 3, getDefaultLabelX(LB), ( 0, 1)),
            newBond(B, 3, 4, flipLabel(getDefaultLabelX(LB)), ( 0,-1)),

            newBond(B, 3, 1, getDefaultLabelX(LB), ( 1, 0)),
            newBond(B, 1, 3, flipLabel(getDefaultLabelX(LB)), (-1, 0)),

            newBond(B, 4, 2, getDefaultLabelX(LB), ( 1, 0)),
            newBond(B, 2, 4, flipLabel(getDefaultLabelX(LB)), (-1, 0)),

            newBond(B, 2, 3, getDefaultLabelZ(LB), (-1, 1)),
            newBond(B, 3, 2, flipLabel(getDefaultLabelZ(LB)), ( 1,-1))
        ]
    )
end


text:
function generateGridData()

  return u
end


text:
includet("OpCodes.jl")
using .OpCodes

data = readline("data/Dec9_data.txt")
getoutput(compute!(IntComputer(data), 1))

# question 2
getoutput(compute!(IntComputer(data), 2))

# test cases
ic = IntComputer("109,1,204,-1,1001,100,1,100,1008,100,16,101,1006,101,0,99")
getoutput(compute!(ic))

ic = IntComputer("1102,34915192,34915192,7,4,7,99,0")
getoutput(compute!(ic))


text:
module NumberUnions

export IEEEFloat,
       MachInt, MachUInt, MachFloat,
       SysInt, SysUInt, SysFloat,
       Integer128, Integer64, Integer32, Integer8,
       IntFloat64, IntFloat32, IntFloat16,
       bytes2Int, bytes2UInt, bytes2Float

import Base: IEEEFloat

const MachInt   = Union{   Int64,   Int32 }
const MachUInt  = Union{  UInt64,  UInt32 }
const MachFloat = Union{ Float64, Float32 }

const IntFloat64  = Union{ Int64, Float64 }
const IntFloat32  = Union{ Int32, Float32 }
const IntFloat16  = Union{ Int16, Float16 }

const SysInt   = Union{ Int128,   Int64,   Int32,  Int16,  Int8 }
const SysUInt  = Union{ UInt128, UInt64,  UInt32, UInt16, UInt8 }
const SysFloat = Union{ Float64, Float32, Float16 }

const Integer128 = Union{ Int128, UInt128 }
const Integer64  = Union{ Int64, UInt64 }
const Integer32  = Union{ Int32, UInt32 }
const Integer16  = Union{ Int16, UInt16 }
const Integer8   = Union{ Int8, UInt8 }

"""
   lookup primitive numeric {Int, UInt, Float} type using sizeof(type)

```julia

indexbysize(nbytes) = 1 + (nbytes >> 1) - (nbytes >> 3) - ((nbytes >> 4) << 1);

floatsizes = (sizeof(Float16), sizeof(Float32), sizeof(Float64))
(2, 4, 8)

indices_for_float_types =  [indexbysize(nbytes) for nbytes in floatsizes];

SysFloatsBySize[ indices_for_float_types ] == [Float16, Float32, Float64]
true

intsizes = (sizeof(Int8), sizeof(Int16), sizeof(Int32), sizeof(Int64), sizeof(Int128))
(1, 2, 4, 8, 16)

indices_for_int_types =  [indexbysize(nbytes) for nbytes in intsizes];

SysIntsBySize[ indicies_for_int_types ] == [Int8,    Int16,   Int32,   Int64,   Int128]
true
```
""" bytes2Int, bytes2UInt, bytes2Float

const SysIntSizes   = [ Int8,    Int16,   Int32,   Int64,   Int128  ]
const SysUIntSizes  = [ UInt8,   UInt16,  UInt32,  UInt64,  UInt128 ]
const SysFloatSizes = [ Float16, Float16, Float32, Float64, Float64 ]

@inline indexbysize(nbytes) = 1 + (nbytes >> 1) - (nbytes >> 3) - ((nbytes >> 4) << 1)

bytes2Int(nbytes) = SysIntSizes[ indexbysize(nbytes) ]
bytes2UInt(nbytes) = SysUIntSizes[ indexbysize(nbytes) ]
bytes2Float(nbytes) = SysFloatSizes[ indexbysize(nbytes) ]

end # module


text: println("There are no tests defined AWSSDK.jl. See tests in AWSCore.jl")

max_stars_repo_name: JuliaCloud/AWSSDK

max_issues_repo_name: JuliaCloud/AWSSDK

max_forks_repo_name: JuliaCloud/AWSSDK


text:
"""
Opaque Decoder struct
"""
mutable struct OpusDecoder  # mutable so that finalizer can be applied
    v::Ptr{Cvoid}
    fs::Int32
    channels::Cint

    """
    Create new OpusDecoder object with given samplerate and channels
    """
    function OpusDecoder(samplerate, channels)
        errorptr = Ref{Cint}(0);
        # Create new decoder object with the given samplerate and channel
        decptr = ccall((:opus_decoder_create,libopus), Ptr{Cvoid}, (Int32, Cint, Ref{Cint}), samplerate, channels, errorptr)
        err = errorptr[]
        dec = new(decptr, samplerate, channels)
        if err != OPUS_OK
            error("opus_decoder_create() failed: $(OPUS_ERROR_MESSAGE_STRS[err])")
        end

        # Register finalizer to cleanup this decoder
        finalizer(dec) do x
            ccall((:opus_decoder_destroy,libopus),Cvoid,(Ptr{Cvoid},),x.v)
        end
        return dec
    end
end

function get_nb_samples(data, fs)
    num_samples = ccall((:opus_packet_get_nb_samples, libopus), Cint, (Ptr{UInt8}, Int32, Int32), data, length(data), fs)
    if num_samples < 0
        error("opus_packet_get_nb_samples() failed: $(OPUS_ERROR_MESSAGE_STRS[num_samples])")
    end
    return num_samples
end

function decode_packet(dec::OpusDecoder, packet::Vector{UInt8})
    packet_samples = get_nb_samples(packet, dec.fs)
    output = Vector{Float32}(undef, packet_samples*dec.channels)

    #print("opus_decode_float: ")
    num_samples = ccall((:opus_decode_float,libopus), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Int32, Ptr{Float32}, Cint, Cint),
                        dec.v, packet, length(packet), output, packet_samples*dec.channels, 0)
    if num_samples < 0
        error("opus_decode_float() failed: $(OPUS_ERROR_MESSAGE_STRS[num_samples])")
    end
    return output
end

"""
Packets go in, floating-point audio stream comes out
"""
function decode_all_packets(dec::OpusDecoder, packets::Vector{Vector{UInt8}})
    # Skip past header packets
    start_idx = 1
    while is_header_packet(packets[start_idx])
        start_idx += 1
    end

    # Figure out the length of each packet
    packet_lens = map(packet -> get_nb_samples(packet, dec.fs)*dec.channels, packets[start_idx:end])
    packet_lens = vcat(zeros(Int32, start_idx - 1), packet_lens)

    # Allocate output now that we know the length of all our audio
    output = Vector{Float32}(undef, sum(packet_lens))

    # Decode each packet into its corresponding chunk of output
    out_idx = 1
    for idx in start_idx:length(packets)
        output[out_idx:out_idx + packet_lens[idx] - 1] = decode_packet(dec, packets[idx])
        out_idx += packet_lens[idx]
    end

    # Return our hard-earned goodness, reshaping if we're stereo!
    if dec.channels > 1
        return reshape(output, (Int64(dec.channels), div(length(output), dec.channels)))'
    end
    return output''
end

"""
Returns (audio, fs) unless the ogg file has no opus streams
"""
function load(fio::IO)
    audio = nothing
    packets = Ogg.load(fio)
    for serial in keys(packets)
        opus_head = OpusHead()
        opus_tags = OpusTags()
        try
            # Find the first stream that is Opus and decode it
            opus_head = OpusHead(packets[serial][1])
            opus_tags = OpusTags(packets[serial][2])
        # TODO: throw a more specific exception. Catching all exceptions here
        # makes things hard to debug if there's a problem (or Julia syntax changes)
        catch
            continue
        end
        dec = Opus.OpusDecoder(48000, opus_head.channels)
        audio = decode_all_packets(dec, packets[serial])
        break
    end

    if audio == nothing
        error("Could not find any Opus streams!")
    end
    return audio, 48000
end

function load(file_path::Union{File{format"OPUS"},AbstractString})
    open(file_path) do fio
        return load(fio)
    end
end


text:
using Primes

@testset "Rational reconstruction" begin
    ms = [107, 199, 509, 2^31 - 1]
    as = [QQ(0), QQ(1), QQ(1, 2), QQ(3, 2), QQ(4, 5), QQ(1, 7), QQ(6)]

    for m in ms
        for a in as
            gf = GF(m)
            ac = lift( numerator(a) // gf(denominator(a)) )
            ar1 = Groebner.rational_reconstruction(ac, BigInt(m))
            ar2 = Groebner.rational_reconstruction(Int(ac), m)
            @test ar1 == ar2 == a
        end
    end

    # random tests
    ms = Primes.primes(2^20, 2^21)[1:3]
    as = [QQ(rand(1:100), rand(1:100)) for _ in 1:10]
    for m in ms
        for a in as
            gf = GF(m)
            ac = lift( numerator(a) // gf(denominator(a)) )
            ar1 = Groebner.rational_reconstruction(ac, BigInt(m))
            ar2 = Groebner.rational_reconstruction(Int(ac), m)
            @test ar1 == ar2 == a
        end
    end

end


text:
#===============================================================================

        Linear transformation of the underlying process, disturbed
        by Gaussian noise. The main struct introduced here is
        `LinearGsnObs`.

===============================================================================#
@doc raw"""
    struct LinearGsnObs{Tag,D,T,FPT,S,R,K} <: Observation{D,T}
        L::S
        μ::T
        Σ::R
        obs::T
        t::Float64
        full_obs::Bool
        θ::Vector{K}
    end

Observation of the underlying process that is of the form:
```math
LX+ξ,\quad\mbox{where }ξ∼N(μ,Σ),
```
and $L\in\mathbb{R}^{D×d}$, $Σ\in\mathbb{R}^{D×D}$ and $μ\in\mathbb{R}^D$, with
$d>0$ the dimension of the underlying process $X$. `FPT` stores information
about first-passage times. `full_obs` is an indicator for whether it is a full
observation of the process (as it grants the use of Markov Property). `θ` is a
container that may contain parameters that enter `L`, `μ` and `Σ`, whereas `Tag`
is a disambiguation flag used at compile time or for multiple dispatch to
differentiate between different observation types.

    function LinearGsnObs(
        t,
        obs::T;
        L=_default_L(obs), # defaults to I
        Σ=_default_Σ(obs), # defaults to ϵI, with ϵ=1e-11
        μ=_default_μ(obs), # defaults to 0 vector
        fpt=NoFirstPassageTimes(),
        full_obs=false,
        Tag=0,
        θ=[],
    ) where T

Base constructor of LinearGsnObs. `t` is the time of the observation and `obs`
is the actual observation. The remaining, named parameters are self-explanatory.
"""
struct LinearGsnObs{Tag,D,T,FPT,TL,Tμ,TΣ,K} <: Observation{D,T}
    L::TL
    μ::Tμ
    Σ::TΣ
    obs::T
    t::Float64
    full_obs::Bool
    θ::Vector{K}

    function LinearGsnObs(
            t, obs::T, L::TL, Σ::TΣ, μ::Tμ, fpt::FPT, full_obs, θ::Vector{K},
            Tag,
        ) where {T,TL,TΣ,Tμ,FPT,K}
        D = length(obs)
        new{Tag,D,T,FPT,TL,Tμ,TΣ,K}(L, μ, Σ, obs, t, full_obs, θ)
    end
end

_default_L(obs) = _default_L(obs, ismutable(obs))
_default_L(obs::T, ::Val{true}) where T = Diagonal(
    repeat([one(eltype(T))], length(obs))
)
_default_L(obs::T, ::Val{false}) where T = SDiagonal{length(obs)}(I)

_default_Σ(obs, ϵ=1e-11) = _default_Σ(obs, ismutable(obs), ϵ)
_default_Σ(obs::T, ::Val{true}, ϵ=1e-11) where T = ϵ*Diagonal(
    repeat([one(eltype(T))], length(obs))
)
_default_Σ(obs::T, ::Val{false}, ϵ=1e-11) where T = ϵ*SDiagonal{length(obs)}(I)

_default_μ(obs) = _default_μ(obs, ismutable(obs))
_default_μ(obs::T, ::Val{true}) where T = zeros(eltype(T),length(obs))
_default_μ(obs::T, ::Val{false}) where T = zero(T)

function LinearGsnObs(
        t,
        obs::T;
        L=_default_L(obs),
        Σ=_default_Σ(obs),
        μ=_default_μ(obs),
        fpt=NoFirstPassageTimes(),
        full_obs=false,
        Tag=0,
        θ=[],
    ) where T
    @assert length(θ) == 0 || Tag > 0
    LinearGsnObs(t, obs, L, Σ, μ, fpt, full_obs, θ, Tag)
end

"""
    fpt_info(::LinearGsnObs{D,T,Tag,FPT}) where {D,T,Tag,FPT}

Return information about first-passage times
"""
fpt_info(::LinearGsnObs{Tag,D,T,FPT}) where {Tag,D,T,FPT} = FPT

get_tag(::LinearGsnObs{Tag}) where Tag = Tag

function Base.summary(io::IO, obs::LinearGsnObs)
    _L, _μ, _Σ, v = L(obs), μ(obs), Σ(obs), ν(obs)
    t, full_obs = obs.t, obs.full_obs
    θ = parameters(obs)
    parametric_type = get_tag(obs) > 0
    L_size = size(_L)
    println(io, repeat("⏤", 40 ))
    println(
        io,
        "|Observation `v = Lx+ξ`, where `L` is a $L_size-matrix, `x` is a ",
        "state of the stochastic process and `ξ`∼N(μ,Σ)."
    )
    println(io, "|...")
    println(io, "|| ν: $v (observation),\n|| \t→ typeof(ν): ", typeof(v),",")
    println(io, "|| made at time $t.")
    println(io, "|...")
    println(io, "|L: $_L,\n|\t→ typeof(L): ", typeof(_L))
    println(io, "|μ: $_μ,\n|\t→ typeof(μ): ", typeof(_μ))
    println(io, "|Σ: $_Σ,\n|\t→ typeof(Σ): ", typeof(_Σ))
    println(io, "|...")
    println(io, "|This is ", (full_obs ? "" : "NOT ") , "an exact observation.")
    println(io, "|...")
    println(
        io,
        "|It ",
        parametric_type ?
        "depends on additional parameters, which are set to: $θ." :
        "does not depend on any additional parameters."
    )
    println(io, "|...")
    _summary_fpt(io, fpt_info(obs); prepend="|")
    println(io, repeat("⋆ ", 10))
end

Base.summary(o::Observation) = summary(Base.stdout, o)

"""
    L(o::LinearGsnObs)

Return matrix L from the observation scheme ν = Lx+ξ, where ξ∼N(μ,Σ)
"""
L(o::LinearGsnObs) = o.L

"""
    μ(o::LinearGsnObs)

Return vector μ from the observation scheme ν = Lx+ξ, where ξ∼N(μ,Σ)
"""
μ(o::LinearGsnObs) = o.μ

"""
    Σ(o::LinearGsnObs)

Return matrix Σ from the observation scheme ν = Lx+ξ, where ξ∼N(μ,Σ)
"""
Σ(o::LinearGsnObs) = o.Σ

"""
    Λ(o::LinearGsnObs)

Return matrix Λ:=Σ⁻¹ from the observation scheme ν = Lx+ξ, where ξ∼N(μ,Σ)
"""
Λ(o::LinearGsnObs) = inv(o.Σ)

"""
    ν(o::Observation)

Return the observation
"""
ν(o::Observation) = o.obs

"""
    obs(o::Observation)

Alias to ν. Return the observation.
"""
const obs = ν

clone(o::LinearGsnObs{0}, args...) = o
update_params!(o::LinearGsnObs{0}, new_params...) = o

function update_params!(o::LinearGsnObs, new_params::AbstractArray)
    update_params!(o, new_params...)
end
function update_params!(o::LinearGsnObs, new_params...)
    @assert length(new_params) == length(o.θ)
    o.θ .= new_params
    o
end

var_parameter_pos(o::Observation) = tuple()
parameters(o::Observation) = tuple(o.θ...)

@doc raw"""
    Base.rand([rng::Random.AbstractRNG], o::LinearGsnObs, X)

Sample an observation according to
```math
V ∼ L X + ξ,\quad \mbox{where } ξ∼N(μ,Σ),
```
with $L$, $μ$ and $Σ$ defined in `o`.
"""
function Base.rand(rng::Random.AbstractRNG, o::LinearGsnObs, x)
    o.full_obs && return x
    L(o) * x + rand(Gaussian(μ(o), Σ(o)))
end

Base.rand(o::LinearGsnObs, x) = rand(Random.GLOBAL_RNG, o, x)


text:
@info "Test iterator pattern..."

for i = AddOne(0, 10)
    println(i)
end

println()

a = [5,4,3,2,1]
for i = ArrayIterator(a, 1, 3)
    println(i)
end

println()

for i = ArrayIterator(a, 2, 3)
    println(i)
end


text:
abstract type AbstractElemRestriction end

"""
    ElemRestrictionNone()

Returns the singleton object corresponding to libCEED's `CEED_ELEMRESTRICTION_NONE`
"""
struct ElemRestrictionNone <: AbstractElemRestriction end
Base.getindex(::ElemRestrictionNone) = C.CEED_ELEMRESTRICTION_NONE[]

"""
    ElemRestriction

Wraps a `CeedElemRestriction` object, representing the restriction from local vectors to
elements. An `ElemRestriction` object can be created using [`create_elem_restriction`](@ref)
or [`create_elem_restriction_strided`](@ref).
"""
mutable struct ElemRestriction <: AbstractElemRestriction
    ref::RefValue{C.CeedElemRestriction}
    function ElemRestriction(ref)
        obj = new(ref)
        finalizer(obj) do x
            # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
            destroy(x)
        end
        return obj
    end
end
destroy(r::ElemRestriction) = C.CeedElemRestrictionDestroy(r.ref) # COV_EXCL_LINE
Base.getindex(r::ElemRestriction) = r.ref[]
Base.show(io::IO, ::MIME"text/plain", e::ElemRestriction) =
    ceed_show(io, e, C.CeedElemRestrictionView)

@doc raw"""
    create_elem_restriction(
        ceed::Ceed,
        nelem,
        elemsize,
        ncomp,
        compstride,
        lsize,
        offsets::AbstractArray{CeedInt},
        mtype::MemType=MEM_HOST,
        cmode::CopyMode=COPY_VALUES,
    )

Create a `CeedElemRestriction`.

!!! warning "Zero-based indexing"
    In the below notation, we are using **0-based indexing**. libCEED expects the offset
    indices to be 0-based.

# Arguments:
- `ceed`:       The [`Ceed`](@ref) object
- `nelem`:      Number of elements described in the `offsets` array
- `elemsize`:   Size (number of "nodes") per element
- `ncomp`:      Number of field components per interpolation node (1 for scalar fields)
- `compstride`: Stride between components for the same L-vector "node". Data for node $i$,
                component $j$, element $k$ can be found in the L-vector at index `offsets[i
                + k*elemsize] + j*compstride`.
- `lsize`:      The size of the L-vector. This vector may be larger than the elements and
                fields given by this restriction.
- `offsets`:    Array of shape `(elemsize, nelem)`. Column $i$ holds the ordered list of the
                offsets (into the input [`CeedVector`](@ref)) for the unknowns corresponding
                to element $i$, where $0 \leq i < \textit{nelem}$. All offsets must be in
                the range $[0, \textit{lsize} - 1]$.
- `mtype`:      Memory type of the `offsets` array, see [`MemType`](@ref)
- `cmode`:      Copy mode for the `offsets` array, see [`CopyMode`](@ref)
"""
function create_elem_restriction(
    c::Ceed,
    nelem,
    elemsize,
    ncomp,
    compstride,
    lsize,
    offsets::AbstractArray{CeedInt};
    mtype::MemType=MEM_HOST,
    cmode::CopyMode=COPY_VALUES,
)
    ref = Ref{C.CeedElemRestriction}()
    C.CeedElemRestrictionCreate(
        c[],
        nelem,
        elemsize,
        ncomp,
        compstride,
        lsize,
        mtype,
        cmode,
        offsets,
        ref,
    )
    ElemRestriction(ref)
end

@doc raw"""
    create_elem_restriction_strided(ceed::Ceed, nelem, elemsize, ncomp, lsize, strides)

Create a strided `CeedElemRestriction`.

!!! warning "Zero-based indexing"
    In the below notation, we are using **0-based indexing**. libCEED expects the offset
    indices to be 0-based.

# Arguments:
- `ceed`:     The [`Ceed`](@ref) object
- `nelem`:    Number of elements described by the restriction
- `elemsize`: Size (number of "nodes") per element
- `ncomp`:    Number of field components per interpolation node (1 for scalar fields)
- `lsize`:    The size of the L-vector. This vector may be larger than the elements and
              fields given by this restriction.
- `strides`:  Array for strides between [nodes, components, elements]. Data for node $i$,
              component $j$, element $k$ can be found in the L-vector at index `i*strides[0]
              + j*strides[1] + k*strides[2]`. [`STRIDES_BACKEND`](@ref) may be used with
              vectors created by a Ceed backend.
"""
function create_elem_restriction_strided(c::Ceed, nelem, elemsize, ncomp, lsize, strides)
    ref = Ref{C.CeedElemRestriction}()
    C.CeedElemRestrictionCreateStrided(c[], nelem, elemsize, ncomp, lsize, strides, ref)
    ElemRestriction(ref)
end

"""
    apply!(
        r::ElemRestriction,
        u::CeedVector,
        ru::CeedVector;
        tmode=NOTRANSPOSE,
        request=RequestImmediate(),
    )

Use the [`ElemRestriction`](@ref) to convert from L-vector to an E-vector (or apply the
tranpose operation). The input [`CeedVector`](@ref) is `u` and the result stored in `ru`.

If `tmode` is `TRANSPOSE`, then the result is added to `ru`. If `tmode` is `NOTRANSPOSE`,
then `ru` is overwritten with the result.
"""
function apply!(
    r::ElemRestriction,
    u::CeedVector,
    ru::CeedVector;
    tmode=NOTRANSPOSE,
    request=RequestImmediate(),
)
    C.CeedElemRestrictionApply(r[], tmode, u[], ru[], request[])
end

"""
    apply(r::ElemRestriction, u::AbstractVector; tmode=NOTRANSPOSE)

Use the [`ElemRestriction`](@ref) to convert from L-vector to an E-vector (or apply the
tranpose operation). The input is given by `u`, and the result is returned as an array of
type `Vector{CeedScalar}`.
"""
function apply(r::ElemRestriction, u::AbstractVector; tmode=NOTRANSPOSE)
    ceed_ref = Ref{C.Ceed}()
    ccall(
        (:CeedElemRestrictionGetCeed, C.libceed),
        Cint,
        (C.CeedElemRestriction, Ptr{C.Ceed}),
        r[],
        ceed_ref,
    )
    c = Ceed(ceed_ref)
    uv = CeedVector(c, u)
    if tmode == NOTRANSPOSE
        ruv = create_evector(r)
    else
        ruv = create_lvector(r)
        ruv[] = 0.0
    end
    apply!(r, uv, ruv; tmode=tmode)
    Vector(ruv)
end

"""
    create_evector(r::ElemRestriction)

Return a new [`CeedVector`](@ref) E-vector.
"""
function create_evector(r::ElemRestriction)
    ref = Ref{C.CeedVector}()
    C.CeedElemRestrictionCreateVector(r[], C_NULL, ref)
    CeedVector(ref)
end

"""
    create_lvector(r::ElemRestriction)

Return a new [`CeedVector`](@ref) L-vector.
"""
function create_lvector(r::ElemRestriction)
    ref = Ref{C.CeedVector}()
    C.CeedElemRestrictionCreateVector(r[], ref, C_NULL)
    CeedVector(ref)
end

"""
    create_vectors(r::ElemRestriction)

Return an (L-vector, E-vector) pair.
"""
function create_vectors(r::ElemRestriction)
    l_ref = Ref{C.CeedVector}()
    e_ref = Ref{C.CeedVector}()
    C.CeedElemRestrictionCreateVector(r[], l_ref, e_ref)
    CeedVector(l_ref), CeedVector(e_ref)
end

"""
    getcompstride(r::ElemRestriction)

Get the L-vector component stride.
"""
function getcompstride(r::ElemRestriction)
    lsize = Ref{CeedInt}()
    C.CeedElemRestrictionGetCompStride(r[], lsize)
    lsize[]
end

"""
    getnumelements(r::ElemRestriction)

Get the total number of elements in the range of an [`ElemRestriction`](@ref).
"""
function getnumelements(r::ElemRestriction)
    result = Ref{CeedInt}()
    C.CeedElemRestrictionGetNumElements(r[], result)
    result[]
end

"""
    getelementsize(r::ElemRestriction)

Get the size of elements in the given [`ElemRestriction`](@ref).
"""
function getelementsize(r::ElemRestriction)
    result = Ref{CeedInt}()
    C.CeedElemRestrictionGetElementSize(r[], result)
    result[]
end

"""
    getlvectorsize(r::ElemRestriction)

Get the size of an L-vector for the given [`ElemRestriction`](@ref).
"""
function getlvectorsize(r::ElemRestriction)
    result = Ref{CeedInt}()
    C.CeedElemRestrictionGetLVectorSize(r[], result)
    result[]
end

"""
    getnumcomponents(r::ElemRestriction)

Get the number of components in the elements of an [`ElemRestriction`](@ref).
"""
function getnumcomponents(r::ElemRestriction)
    result = Ref{CeedInt}()
    C.CeedElemRestrictionGetNumComponents(r[], result)
    result[]
end

"""
    getmultiplicity!(r::ElemRestriction, v::AbstractCeedVector)

Get the multiplicity of nodes in an [`ElemRestriction`](@ref). The [`CeedVector`](@ref) `v`
should be an L-vector (i.e. `length(v) == getlvectorsize(r)`, see [`create_lvector`](@ref)).
"""
function getmultiplicity!(r::ElemRestriction, v::AbstractCeedVector)
    @assert length(v) == getlvectorsize(r)
    C.CeedElemRestrictionGetMultiplicity(r[], v[])
end

"""
    getmultiplicity(r::ElemRestriction)

Convenience function to get the multiplicity of nodes in the [`ElemRestriction`](@ref),
where the result is returned in a newly allocated Julia `Vector{CeedScalar}` (see also
[`getmultiplicity!`](@ref)).
"""
function getmultiplicity(r::ElemRestriction)
    v = create_lvector(r)
    getmultiplicity!(r, v)
    Vector(v)
end

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
rawdata = []
open("data.txt") do file
    for line in eachline(file)
        push!(rawdata,parse(Int,line))
    end
end

function printSuccess(a,b,c)
    println("Found it")
    println("a=", a)
    println("b=", b)
    println("c=", c)
    answer = a * b * c
    println("answer=", answer)
end

sorteddata = sort(rawdata)

for z in 200:-1:1
    x = 1
    y = x+1
    a = sorteddata[x]
    b = sorteddata[y]
    c = sorteddata[z]

    while (a + b + c != 2020)
        println("x=", x," y=", y, " z=", z, " a=", a, " b=", b, " c=", c)
        if (a + b + c > 2020)
            break
        end
        if (a + b + c < 2020)
            for x in 1:100
                if (a + b + c < 2020)
                    if (a + b + c < 2020)
                        y += 1
                        b = sorteddata[y]
                    end
                    if (a + b + c > 2020)
                        break
                    end
                    if a + b + c == 2020
                        printSuccess(a,b,c)
                        exit(0)
                    end
                end
                if (a + b + c < 2020)
                    x += 1
                    a = sorteddata[y]
                end
                if (a + b + c > 2020)
                    break
                end
                if a + b + c == 2020
                    printSuccess(a,b,c)
                    exit(0)
                end
                x += 1
                a = sorteddata[x]
            end
        end

        if (a + b + c == 2020)
            printSuccess(a,b,c)
            exit(0)
        end
    end
end


text:
mutable struct GenericModel{T} <: AbstractModel{T}
    parameters::ParameterVector{T}
    steady_state::ParameterVector{T}
    keys::Dict{Symbol, Int}
    endogenous_states::Dict{Symbol,Int}             # these fields used to create matrices in the
    exogenous_shocks::Dict{Symbol,Int}              # measurement and equilibrium condition equations.
    observables::Dict{Symbol,Int}                   #
    pseudo_observables::Dict{Symbol,Int}            #
    spec::String
    subspec::String
    settings::Dict{Symbol, Setting}
    test_settings::Dict{Symbol, Setting}
    rng::MersenneTwister
    testing::Bool
    observable_mappings::Dict{Symbol, Observable}
    pseudo_observable_mappings::Dict{Symbol, PseudoObservable}
end

function GenericModel(spec::String = "", subspec::String = ""; custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(), testing = false)
    settings = Dict{Symbol, Setting}()
    test_settings = Dict{Symbol, Setting}()
    rng = MersenneTwister(0)

    m = GenericModel{Float64}(
        Vector{AbstractParameter{Float64}}(), Vector{AbstractParameter{Float64}}(),
        Dict{Symbol, Int}(),
        Dict{Symbol, Int}(), Dict{Symbol, Int}(), Dict{Symbol, Int}(), Dict{Symbol, Int}(),
        spec, subspec, settings, test_settings, rng, testing,
        Dict{Symbol, Observable}(), Dict{Symbol, PseudoObservable}())

    model_settings!(m)
    default_test_settings!(m)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    return m
end

description(m::GenericModel) = "A generic model: You can put parameters, states, shocks, observables, pseudo-observables, spec, subpsec, settings inside of me."

function model_settings!(m::GenericModel)

end


text:
###################
###Miscellaneous###
###################
"""
    getLength(x::Int64)

Compute the physical length of a density matrix from the length x of its vector representation.
"""
@inline function getLength(x::Int64)
    return (63-leading_zeros(x))>>1
end

"""
    triming!(A; tol::Float64=2e-15)

Put to 0 all the elements of A smaller in absolute value than tol.
"""
@inline function triming!(A; tol::Float64=2e-15)
    for x in 1:length(A)
        if abs(A[x])<tol
            A[x]=zero(typeof(A[x]))
        end
    end
end
"""
    triming(x; tol::Float64=2e-15)

Return 0 if x is smaller than tol, x else
"""
@inline function triming(x; tol::Float64=2e-15)
    if abs(x)<tol
        return zero(typeof(x))
    else
        return x
    end
end
#######################
###Entropy functions###
#######################
"""
    entropy(eigen)

From a list of svs, return the vN entropy. sv should be real and between 0 and 1
"""
@inline function entropy(eigen)
    entr=0.
    @simd for x in eigen
        #println(x)
        if 1.>x>0.
            #entr-=real(x)*log(real(x))
            entr-=x*log(x)
        end
    end
    return entr
end

"""
    computeEntropy(rho)

Compute the vN entanglement entropy of the density matrix rho, stored as a vector
"""
@inline function computeEntropy(rho)  #Compute vonNeumannEE
    len=getLength(length(rho));
    return entropy(svdvals(reshape(rho, (2^len, 2^len))))
end

"""
    mutualInformation(rho, start, nb1, nbtot)

Compute the mutual Information between start:start+nb1-1 and start+nb1:start+nbtot-1
"""
function mutualInformation(rho, start, nb1, nbtot, Traces)
    if nbtot==nb1
        return 0
    end
    len=getLength(length(rho));
    return computeEntropy(buildReducedDensityMatrix(rho, start+nb1, nbtot-nb1, Traces))+computeEntropy(buildReducedDensityMatrix(rho, start, nb1, Traces))-computeEntropy(buildReducedDensityMatrix(rho, start, nbtot, Traces))
end

"""
    myTrace(rho)

Compute the trace of a DM while removing the imaginary part of diag term
"""
@inline function myTrace(rho)
    tr=0.0
    len=getLength(length(rho));
    stepp=2^len+1; offset=1
    for j=1:2^len
        @inbounds rho[offset]=rho[offset].re
        @inbounds tr+=rho[offset].re
        offset+=stepp
    end
    return tr
end

"""
    normalizeDM!(rho)

Normalize a DM stocked as a vector by its trace
"""
function normalizeDM!(rho)
    a=myTrace(rho)
    @. rho.=rho./a
end

max_stars_repo_name: LoicHerviou/MLQuantumMarginal

max_issues_repo_name: LoicHerviou/MLQuantumMarginal

max_forks_repo_name: LoicHerviou/MLQuantumMarginal


text:
using AbstractPlotting
using AbstractPlotting.MakieLayout
using GLMakie

using LinearAlgebra
using Roots
using ForwardDiff
D(f) = x -> ForwardDiff.derivative(f, float(x))

descr = """
An old optimization problem is to find the shortest distance between two
points when the rate of travel differs due to some medium (darker means slower).
In this example, the relative rate can be adjusted (with the slider), and the
various points  (by clicking on and dragging a point). From there, the
user can adjust the crossing point to identify the time.
"""

## From http://juliaplots.org/MakieReferenceImages/gallery//edit_polygon/index.html:
function add_move!(scene, points, pplot)
    idx = Ref(0); dragstart = Ref(false); startpos = Base.RefValue(Point2f0(0))
    on(events(scene).mousedrag) do drag
        if ispressed(scene, Mouse.left)
            if drag == Mouse.down
                plot, _idx = mouse_selection(scene)
                if plot == pplot
                    idx[] = _idx; dragstart[] = true
                    startpos[] = to_world(scene, Point2f0(scene.events.mouseposition[]))
                end
            elseif drag == Mouse.pressed && dragstart[] && checkbounds(Bool, points[], idx[])
                pos = to_world(scene, Point2f0(scene.events.mouseposition[]))

                # very wierd, but we work with components
                # not vector
                z = zero(eltype(pos))
                x,y = pos

                ptidx = idx[]

                x = clamp(x, 0, 5)

                if ptidx == 1
                    y = clamp(y, z, 5)
                elseif ptidx == 2
                    y = z
                elseif ptidx == 3
                    y = clamp(y, -5, z)
                end

                points[][idx[]] = [x,y]
                points[] = points[]
            end
        else
            dragstart[] = false
        end
        return
    end
end

upperpoints = Point2f0[(0,0), (5, 0), (5, 5), (0,5)]
lowerpoints = (Point2f0[(0,0), (5, 0), (5, -5), (0,-5)])

points = Node(Point2f0[(1, 4), (3, 0), (4,-4.0)])

# where we lay our scene:
scene, layout = layoutscene()
layout.halign = :left
layout.valign = :top

p = layout[1:3, 1] = LScene(scene)
rowsize!(layout, 1, Auto(1))
colsize!(layout, 1, Auto(1))

flayout = GridLayout()
layout[1,2] = flayout

flayout[1,1] = LText(scene, chomp(descr))

details = flayout[2, 1] = LText(scene, "...")

λᵣ = flayout[3,1] =  LText(scene, "λ = v₁/v₂ = 1")
λ = flayout[4,1] = LSlider(scene, range = -3:0.2:3, startvalue = 0.0)

tm = lift(λ.value, points) do λ, pts
    x0,y0 = pts[1]
    x, y = pts[2]
    x1, y1 = pts[3]

    v1 = 1
    v2 = v1/2.0^λ

    t(x) = sqrt((x-x0)^2 + y0^2)/v1 + sqrt((x1-x)^2 + y1^2)/v2
    val = t(x)

    details.text[] = "Time: $(round(val, digits=2)) units"

    val

end

a = lift(λ.value, points) do λ, pts
    x0,y0 = pts[1]
    x1, y1 = pts[3]
    v1 = 1
    v2 = v1/2.0^λ

    t(x) = sqrt((x-x0)^2 + y0^2)/v1 + sqrt((x1-x)^2 + y1^2)/v2
    x′ = fzero(D(t), x0, x1)
    t(x′)
end

λcolor = lift(λ.value) do val
    # val = v1/v2 ∈ [1/8, 8]
    n = floor(Int, 50 - val/3 * 25)
    Symbol("gray" * string(n))
end

linecolor = lift(a) do val
    abs(val - tm[]) <= 1e-2 ? :green : :white
end

on(λ.value) do val
    v = round(2.0^(val), digits=2)
    txt = "λ = v₁/v₂ = $v"
    λᵣ.text[] = txt

end

poly!(p.scene, upperpoints, color = :gray50) # neutral color
poly!(p.scene, lowerpoints, color = λcolor) # color depends on λ

lines!(p.scene, Point2f0[(0,0), (5, 0)], color=:black, strokewidth=5, strokecolor=:black, raw=true)
lines!(p.scene, points, color = linecolor, strokewidth = 10, markersize = 0.05, strokecolor = :black, raw = true)
scatter!(p.scene, points, color = linecolor, strokewidth = 10, markersize = 0.05, strokecolor = :black, raw = true)
xlims!(p.scene, (0,5))
ylims!(p.scene, (-5,5))

add_move!(p.scene, points, p.scene[end])

scene


text:
abstract type Field end

Base.:/(x::F, y::F) where F <: Field = x * inv(y)

Base.one(x::F) where F <: Field = one(F)
Base.zero(x::F) where F <: Field = zero(F)

square(x::Field) = x*x
Base.inv(x::Field) = x^(order(x) - 1)

Base.literal_pow(::typeof(^), x::F, ::Val{2}) where F <: Field = square(x)
Base.literal_pow(::typeof(^), x::F, ::Val{0}) where F <: Field = one(F)

abstract type BinaryField <: Field end

Base.:-(x::F, y::F) where F <: BinaryField = x + y

<|(::Type{F}, x) where F <: BinaryField = convert(F, x)

function Base.convert(::Type{F}, x::Integer) where F <: BinaryField
    if x == 0
        return zero(F)
    elseif x == 1
        return one(F)
    else
        error("Conversion of $x not possible")
    end
end

function Base.:^(g::F, a::Integer) where F <: BinaryField

    if a == 0
        return one(F)
    end

    e = tobits(a)

    r = length(e) - 1

    x = g

    for i in 2:r+1

        x = square(x)

        if e[i] == 1
            x = g*x
        end
    end

    return x
end

function (::Type{F})(x::Integer) where F <: BinaryField
    if x == 0
        return zero(F)
    elseif x==1
        return one(F)
    else
        error("Not supported")
    end
end

function Base.convert(::Type{F}, x) where F <: BinaryField

    N = bitlength(F)

    bits = tobits(x)

    return convert(F, bits[end - N + 1:end])
end

Base.convert(::Type{F}, x::BinaryField) where F <: BinaryField = error("Can't be done blindly.")
Base.convert(::Type{F}, x::F) where F <: BinaryField = x

tobits(::Type{F}, x::BinaryField) where F <: BinaryField = convert(BitVector, x)

Base.show(io::IO, x::BinaryField) = print(io, join(i ? "1" : "0" for i in tobits(x)))

Base.isless(x::F, y::F) where F <: BinaryField = tobits(x) < tobits(y)

abstract type PrimeField <: Field end

<|(::Type{F}, x) where F <: PrimeField = F(x)

function modulus end
function value end

Base.show(io::IO, x::PrimeField) = print(io, value(x))

modinv(s, q) = mod(gcdx(s, q)[2], q)

Base.inv(x::F) where F <: PrimeField = F(modinv(value(x), modulus(x)))

Base.zero(::Type{F}) where F <: PrimeField = F(0)
Base.one(::Type{F}) where F <: PrimeField = F(1)

Base.:+(x::F, y::F) where F <: PrimeField = F(mod(value(x) + value(y), modulus(F)))
Base.:-(x::F, y::F) where F <: PrimeField = x + F(modulus(F) - value(y))
Base.:-(x::F) where F <: PrimeField = F(modulus(F) - value(x))

Base.:*(x::F, y::F) where F <: PrimeField = F(mod(value(x) * value(y), modulus(F)))
Base.:*(x::Integer, y::F) where F <: PrimeField = F(x) * y
Base.:*(x::F, y::Integer) where F <: PrimeField = x * F(y)

Base.:^(x::F, n::Integer) where F <: PrimeField = F(powermod(value(x), n, modulus(F)))

Base.:(==)(x::F, y::F) where F <: PrimeField = value(x) == value(y)

Base.isless(x::F, y::F) where F <: PrimeField = value(x) < value(y)


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

using Markdown
using InteractiveUtils

# ╔═╡ cac0bc38-ead9-4f99-a939-fc2645874efd
md"""
# Types in Julia
"""

# ╔═╡ ecb2b8d9-ceda-4008-8137-b29bd7ef3ba9
# Uses of types.jl

# ╔═╡ b4f0668c-8b25-4588-9176-52bc7bfe797b
md"""
In this notebook we will start exploring types by implementing a type to implement algorithmic differentiation.
"""

# ╔═╡ 0ece0b55-ede7-4db1-9b48-61e1d4534761
md"""
### What is a type?
"""

# ╔═╡ 6dcc01fe-c551-43d2-930b-ff1c93d15077
md"""
A *type* can be thought of as a label that is associated with data stored in memory; this label tells Julia how to interpret the data. For example:
"""

# ╔═╡ 0d352bcc-dfc1-429f-ad6c-40dc8c4cd038
begin
	x = 1
	y = 1.3
end

# ╔═╡ 81863f06-8dac-48fd-a6e1-cf5df30c63f6
sizeof(x), sizeof(y)

# ╔═╡ 3de40205-6b5d-4d99-b886-1a1dc6a07955
typeof(x), typeof(y)

# ╔═╡ d4c280cf-910c-4bbf-8dca-b01b14711296
# rinterpret the bit-representation of x (an int) as another type, e.g. Float64
z = reinterpret(Float64, x)

# ╔═╡ 7082e444-6268-4721-a278-2511a22daf68
bitstring(x)

# ╔═╡ 43538919-4349-43f1-a196-759fe65d3be3
bitstring(z)

# ╔═╡ e0c1df0c-ff62-4982-8f3d-3ccbb44b8d2f
md"""
Both variables are stored in 8 bytes (64 bits), but one is interpreted as an integer and the other as a floating-point number.  Similarly, a pair of two numbers may be intepreted as a complex number, or an interval, or a dual number, or...; although the same information may be stored (two numbers), we want each of these different *kinds* or *types* of objects to be treated differently.
"""

# ╔═╡ e6281275-0b87-47f4-a081-bb366ced3466
md"""
## Algorithmic differentiation
"""

# ╔═╡ 6056ff4a-858f-447c-85cc-0f3f6e0e53ca
md"""
In the previous notebook we used `ForwardDiff.jl` to automatically differentiate a function. Here we will see how to implement a simple version of this.
"""

# ╔═╡ db17c624-91b6-44b1-bd10-4d72f4cb9776
md"""
The idea is to approximate a (nice enough) function $f$ near a point $a$ by a Taylor series of order 1, i.e. a straight line passing through $(a, f(a))$, with slope equal to the derivative $f'(a)$:

$f(x) \simeq f(a) + \epsilon f'(a)$,

where $\epsilon := x - a$.

We now use this to derive the standard rules for the derivative of a sum and product:

$$f(x) + g(x) \simeq [f(a) + g(a)] + \epsilon [f'(a) + g'(a)]$$

$$f(x) \cdot g(x) \simeq [f(a) \cdot g(a)] + \epsilon [f(a) g'(a) + g(a) f'(a)],$$

where we suppose that $\epsilon$ is small enough that $\epsilon^2 = 0$, or alternatively just "take the linear part".
"""

# ╔═╡ b3b15bf5-115d-4831-9115-36c88e9b899c
md"""
### Defining a composite type
"""

# ╔═╡ 70843c81-6abe-40b2-8d0b-8d779fccedc8
md"""
We see that by using just two pieces of information, namely the value $f(a)$ and the derivative $f'(a)$, we can represent a function $f$ near a given point $a$.

The pair $(f(a), f'(a))$ is often called a **dual number**. We see that is has certain **behaviours** under arithmetic operations. Whenever we have a new behaviour, a *new type is lurking*!

We group the two values into a **composite type**. We can think of a composite type as specifying the structure of a box containing several pieces of information (data) inside. Defining a composite type with two **fields** (pieces of information) has the following syntax:
"""

# ╔═╡ 82c79506-68ea-41dd-888d-840371c08430
struct MyType
    a
    b::Int
end

# ╔═╡ 291315e0-ad0e-4a57-b100-2baa5af04c04
md"""
Here we have additionally specified that the information stored in the field `b` must be of type `Int` using the **type annotation operator**, `::`.
"""

# ╔═╡ f2421fbf-5cfc-4332-be60-b61fa4e02d8c
md"""
Creating an object of that type is accomplished as follows:
"""

# ╔═╡ 710612bd-4856-4e51-b2d6-2bdc972f196d
x1 = MyType(3, 4)

# ╔═╡ 70eb5167-6288-4f84-b140-21b4b329cb45
md"""
We can extract information as follows:
"""

# ╔═╡ 9fb42f3c-c7ff-4f31-b5e6-9e96de16bac5
x1.a

# ╔═╡ 07808ac8-c910-4e92-ac34-c9ff225f4be8
typeof(x1.a)

# ╔═╡ 8777c6a7-90e3-4dc7-ab61-3d12531f2ea9
md"""
#### Exercise 1

1. Define a composite type `Dual` with fields `value` and `deriv` of type `Float64`.

2. Create two `Dual` numbers `x` and `y`.

3. What happens if you try to add `x` and `y` together?

4. Make a function `add` that adds `x` and `y` and returns a new `Dual` number, following the rules we found above.
"""

# ╔═╡ 5488c55f-8d0a-4c9b-8f8b-40ed4fe14e67
struct Dual
    value::Float64
    deriv::Float64
end

# ╔═╡ 4f5e5eea-13df-477e-baea-c37ab6074aea
x_D = Dual(3, 4)

# ╔═╡ da932afe-990f-4ca7-aa81-9d85e1ca3bac
methods(Dual)

# ╔═╡ fac29cb2-dc57-4a06-9728-567169ac3768
convert(Float64, 3)

# ╔═╡ 4118613c-e58a-4f63-81d1-0e197e66cf98
# this raises an error
convert(Float64, "hello")

# ╔═╡ ac8fa4de-d8ef-44f2-b709-aceb1e0ecb04
y_D = Dual(5, 6)

# ╔═╡ 92aa5678-478c-4cb6-9187-3285534aac8c
x_D + y_D

# ╔═╡ 04a21e57-a45c-4e4b-977a-42c6ee228b4a
newadd(x::Dual, y::Dual) = Dual(x.value + y.value, x.deriv + y.deriv)

# ╔═╡ 915bfb27-7234-43ea-a052-093f89d7bdcd
newadd(x_D, y_D)

# ╔═╡ 9b213fb9-382a-4c5b-93fe-22252d29d0c6
md"""
## Implementing arithmetic for a type
"""

# ╔═╡ 8eed2b5c-dde2-437e-8383-4ed206ce72de
md"""
We would like to be able to use `+` and `*` for our new `Dual` type, rather than typing `add(x, y)`. To do so, we need to do the following
"""

# ╔═╡ 62a97a4b-dfa6-404d-958b-b7b4ab7d25ca
# +(x::Dual, y::Dual) = Dual(x.value + y.value, x.deriv + y.deriv)

# ╔═╡ 9b0d361a-8c30-4ec0-a183-ac3890895f0e
# Base.+(x::Dual, y::Dual) = Dual(x.value + y.value, x.deriv + y.deriv)

# ╔═╡ 9b8f41b5-9ade-4f34-8d09-cb548bd6c242
# or
# +(x::Dual, y::Dual) = add(x, y)

# ╔═╡ 5b29f183-bda6-4ffb-bb78-5199980039d3
# or
# function +(x::Dual, y::Dual)
#    return add(x, y)
# end

# ╔═╡ fc624260-ba8a-4057-8a1f-2659e2343654
md"""
In Julia, `+` and `*` are just functions. They are defined in `Base` (a module containing basic function definitions) and must be `import`ed before being **extended**. They consist of many different **methods** (versions):
"""

# ╔═╡ abcbc0a1-33e7-430f-b0a8-e8b65c8666b2
@which +(3)

# ╔═╡ 17498a94-5db9-4c5e-9268-e588fddeae54
-(3)

# ╔═╡ c847f0c2-1eac-4900-b80c-7f89846b2924
methods(+)

# ╔═╡ 6caea9da-3edc-4309-bad8-f348693960b2
x

# ╔═╡ 17cce3ab-1b4f-402d-bcc3-c7905360a441
y

# ╔═╡ 526a204c-6a72-4cde-900a-568667e45c64
x + y

# ╔═╡ 86aa82b1-caae-4124-b48a-14497420f471
md"""
We can add more methods that work on our own types. (We are not allowed to modify their behaviour on combinations of types that to not contain our user-defined types; doing so is known as "type piracy" and can affect other people's code in unexpected ways.)
"""

# ╔═╡ e359dd62-65f0-4033-98c2-aff400ed0bb7
md"""
#### Exercise 2

1. Import the `+` and `*` functions from `Base` and implement them for the `Dual` type.
They should return a new `Dual` object.

2. Check that the number of methods has changed.

3. Use `@which x + y` to check that Julia knows which method to use when adding two `Dual`s.

4. Can you define `x + a` for a `Dual` number `x` and a real number `a`? What happens
"""

# ╔═╡ 86141559-f6f4-43a6-8291-550adf31b25b

# ╔═╡ cee6af83-0079-4acd-9e14-cd7b7eb6947f
begin
	value(f::Dual) = f.value
	deriv(f::Dual) = f.deriv
end

# ╔═╡ 4e222c6a-054b-41e6-bef4-b0726344ee94
let
	∂(f::Dual) = f.deriv
	*(f::Dual, g::Dual) = Dual(value(f) * value(g),
    	                     value(f) * ∂(g) + ∂(f) * value(g))
	42 # return this to avoid conflict
end

# ╔═╡ a1b1b945-ba2e-4668-8698-4df02660d379
md"""
(Try using `@code_llvm`)
"""

# ╔═╡ 8857e0aa-da38-4ae6-8ef9-efa66dd51cfc
begin
	f = Dual(3.0, 4.0)
	g = Dual(5.0, 6.0)
end

# ╔═╡ ad671591-a625-4de2-9b18-84c7773cb4b8
f, g

# ╔═╡ 053e90a7-5fa6-4063-9634-63d48ee16690
f + g

# ╔═╡ d521f0e3-1636-4e8b-b8f3-639a90780690
f * g

# ╔═╡ b7b81d9d-354d-4ec0-b50c-302dbcb8f0c4
md"""
Amazingly, we now have enough to be able to differentiate simple Julia functions involving only `+` and `*`. Define
"""

# ╔═╡ d67fa8a9-f43e-4eca-bb0e-1603c9f7da5d
begin
	a = 3.0
	xx = Dual(a, 1.0)  # "the identity function x ↦ x, with derivative 1"
end

# ╔═╡ 3bc3d5a3-9b9f-4d07-9d59-f24d5c005038
md"""
We initialize the derivative as 1.0 when we make a `Dual`. If we use `x` then we automatically differentiate!
"""

# ╔═╡ fd34be78-af76-4d46-bdbd-19fd19155785
md"""
#### Exercise 3
"""

# ╔═╡ caccd695-118c-46a8-b59f-376e4cbc83ce
md"""
1. Define `a = 3.0` and `xx = Dual(a, 1.0)`.

    (i) Compute `xx + xx`. The result should have the value $2a$ and the derivative $2$ -- write a test that it does so.

    (ii) Do the same for `xx * xx`.

2. Define the function `f(x) = x * x + x`. Compute `f(xx)` and check that it gives the correct value and derivative!

3. Does this work for the function `f(x) = x^2 + x`?  What do you need to do?

4. What happens for `f(x) = x^2 + 2x`? What do you need to do?

5. What should you do for `f(x) = sin(x) + x`?
"""

# ╔═╡ 6418d0a4-8dd3-4816-8e5b-092cb769ebb6
a = 3.0
xx = Dual(a, 1.0)

# ╔═╡ 8ebb42e9-11ec-434d-87c9-4259c26f6526
xx + xx

# ╔═╡ 84b6a870-4515-4e6c-b306-d6b6d2a0fb8d
xx * xx

# ╔═╡ a0865cb0-7077-44ff-a130-eb73cbbc3179
xx * xx + xx  # the function x ↦ x^2 + x,  derivative 2x + 1

# ╔═╡ f1fc7637-4804-4fc1-81ac-14874e66f9d6
ff(x) = x*x*x + x*x + x*x  # x^3 + 2x^2

# ╔═╡ ded7adc0-619f-4735-bd93-9c26c4c21858
ff(xx)

# ╔═╡ 1284374b-ab26-46b8-a74e-119655f5112e
ff(a)

# ╔═╡ 4e06b4f6-9753-4847-8c6c-120208520dc5
n = 4
fff(x) = x^n

# ╔═╡ 420c0ef0-3f7e-4044-be30-c2873636c6e8
dff(x) = 3x^2 + 4x

# ╔═╡ f426bed2-05f1-4598-9817-bb90262911ef
dff(a)

# ╔═╡ 09bfa122-ca5b-4b4a-a9eb-09b1b5482a76
fff(x) = x^2

# ╔═╡ d08a434c-c816-4f43-9c93-2d0eda57723c
fff(xx)

# ╔═╡ 131373df-8e58-49b0-be5f-3cd8d8deb61d
fff(xx)

# ╔═╡ 76a7a87b-0283-4a85-9d88-0eec06ef689f
fff4(x) = x^4

# ╔═╡ d608bc38-5653-4ad5-b6e4-6b2236299e4c
fff4(xx)

# ╔═╡ bf697142-0eaf-4717-9434-eaa89946beb4
md"""
`Base.literal_pow`
"""

# ╔═╡ 108212df-c86f-4b75-8576-65f3964687e9
md"""
#### Exercise 4
"""

# ╔═╡ eb94cc27-ac9e-4202-acdc-248809eebdb9
md"""
1. Define a function `differentiate` that differentiates a function `f` at a point `a` using `Dual` numbers, by following the above pattern. (It should return just the derivative at the given point.)
"""

# ╔═╡ d7f6c80a-4922-42be-9bce-e493e603c373

# ╔═╡ 70ec3ada-f7ae-4526-be35-5aa8cd4533ca
md"""
This is the basis of ("forward-mode") automatic differentiation. The `ForwardDiff.jl` method contains a sophisticated implementation of this method.
"""

# ╔═╡ 5011cf98-128e-4c94-bae0-109016427c50
md"""
### Parametric types
"""

# ╔═╡ d57937d8-3f1e-4f85-bf6b-afc8cfe3e99d
md"""
For simplicity, in the above we fixed the fields in the `Dual` type to be of type `Float64`. By doing so we are actually *losing power*. Instead we should let Julia "fill in" the types.

To do so, we specify that we want to use a **type parameter** `T`. We can think of this as a "special kind of variable" that can only take on certain kinds of values. We specify this with the following syntax:
"""

# ╔═╡ df1c033c-8a99-4b45-bc9c-3287a3a959ec
struct MyType2{T}
    a::T
    b::T
end

# ╔═╡ b0902a74-01e0-456a-a027-a5030d494d7f
md"""
[Note that we have not reused the name `MyType` since Julia *does not allow types to be redefined in a different way*.]
"""

# ╔═╡ 9226e5a3-4445-4420-9b54-be31468e6b07
md"""
Here we are specifying that both fields `a` and `b` must share the same type `T`, but we have not restricted what values `T` can take. When we create an object, Julia will *infer* (work out) the type:
"""

# ╔═╡ ccc258fa-c928-4175-a48e-bd72b7d22f31
x_M = MyType2(3, 4)

# ╔═╡ c9f3dde8-7989-45e8-9831-894d82e57e90
y_M = MyType2(3.1, 4.2)

# ╔═╡ 2d3876d6-5fcf-4bfb-8dba-4bdf252c73fc
z_M = MyType2(1, 5.3)

# ╔═╡ a0d7a33a-e067-4e69-9c70-14005f41a2c6
struct MyType3{S,T}
    a::S
    b::T
end

# ╔═╡ 812149d8-0040-4d67-b2a2-f1304af4b225
x_M2 = MyType3(1, 5.3)

# ╔═╡ 75fa60fe-6583-4a6a-b6b0-06fb81006218
md"""
Note that `x` and `y` have *different* types.
"""

# ╔═╡ 4146cf8f-ab2f-457d-97ef-e6d39e93246e
md"""
We can define functions acting on parametric types without necessarily talking about the type parameter:
"""

# ╔═╡ 3e61a354-ca21-410b-b422-92ea362b8c10
md"""
#### Exercise 5

1. Define a function that takes an object of type `MyType2`, *without* mentioning the type parameter, and returns the sum of the two fields.

   What happens when you apply this function to `x` and `y`?

2. Define a type `Dual2` with a type parameter `T` and the same functions `+` and `*` as before.

3. Define the function `f(x) = x * x + x`. What happens if you pass in `Dual` numbers with different type parameters?
"""

# ╔═╡ 907f65aa-e1a1-4838-a1a0-cdaa0a1cfa45
struct DualReal
    value::Real
    deriv::Real
end

# ╔═╡ fc5839ce-935a-46aa-bbe7-01f203b0fdb5
DualReal(big(1.0), big"1.4")

# ╔═╡ c445bab5-c012-4890-89bc-e7d6950fa7b4
DualReal(3.1, big"1.4")

# ╔═╡ 9008cc1d-e7c2-4b24-ad08-968c085baa36
struct DualUntyped
    value
    deriv
end

# ╔═╡ d5cb6844-6f9d-422e-a241-96450aff70f5
DualUntyped("hello", "David")

# ╔═╡ 7b2f03b4-d3dd-4fa6-a569-c878fc9f6557
struct DualBig
    value::BigFloat
    deriv::BigFloat
end

# ╔═╡ 97097783-2948-4e9c-b769-8085c7943227
begin
	struct Dual3{T<:Real}
	    a::T
    	b::T
	end

	Dual3(a::Real, b::Real) = Dual3(promote(a, b))

	# Dual3(a::Real, b::Real) = Dual3(promote(a, b)...)
end

# ╔═╡ 4bdb5bfc-82c2-4de5-982a-4a69872749d1
Dual3("a", "b")

# ╔═╡ 9c05cca8-76d9-4abd-8774-28e7f41d3b09
Dual3(big"1.0", big"1.4")

# ╔═╡ 05db96f7-3eeb-47aa-a4dc-dd67b0f02a2b
Dual3(1, 3.1)

# ╔═╡ b490d647-85f5-4b82-ab29-2d18ce1c7db4
methods(Dual3)

# ╔═╡ d981b89c-4833-4e38-9f80-246ed5f09445
Dual3(3, 4.5)

# ╔═╡ 9bfe7470-afec-4a24-9f96-43272e2f9ecd
promote(1, 3.4)

# ╔═╡ 61cc131f-2d02-4de1-b79a-9105d6decc15
@which promote_type(Int, Float64)

# ╔═╡ a5c1c463-2393-499d-852c-16e69549819a
Complex(3, 4.5)

# ╔═╡ 749ba068-6deb-4a7d-a650-412e5fbdb670
@which Complex(3, 4.5)

# ╔═╡ 226f5af6-702f-431c-a269-5dd4cafa2699
v = Dual(3, 4)

# ╔═╡ c3831c37-8e09-4846-963d-27e1d4a824ac
v isa Number

# ╔═╡ fd8bf8f8-68dd-4925-aa9b-8eca1937c577
subtypes(Real)

# ╔═╡ e3d64165-170e-4ee0-a9bc-ef759ce6afd7
subtypes(AbstractFloat)

# ╔═╡ 713dc49e-6038-40b7-a917-299277f7ea7a
subtypes(Float64)

# ╔═╡ 68d54b7c-4003-49e8-b2f9-666b882a5737
isconcretetype(Float64)

# ╔═╡ 7ae57ebb-d6bc-4b9b-913a-f9c4bb11d380
q = 3.1

# ╔═╡ 99dd8229-8658-4efe-9133-3eb8cd00b6b3
typeof(q)

# ╔═╡ ef28ca54-38dd-4e8e-a5e8-897c802f81e1
q isa AbstractFloat

# ╔═╡ a31cfb26-9914-441c-aea8-6f341f2635fb
gg(x::AbstractFloat) = sqrt(x)

# ╔═╡ b05dcbc4-0149-48a3-8533-1bd6515d0b32
gg(x::T) where {T <: AbstractFloat} = sqrt(x)

# ╔═╡ f44efbd4-7ce2-44db-ac0e-eeb4aeafef70
gg(π)

# ╔═╡ 7df9a91e-01d4-471f-9ad5-55ac56c245d3
struct Hello{T}
    a::T
    b::T
    c::T
end

# ╔═╡ ef1c3d05-7020-4339-b431-63ad63b244cc
add(x::Hello, y::Hello) = x.a + y.a, x.b + y.b, ...

# ╔═╡ 186148a4-698e-460b-82dd-7b203698de45
w = Hello(3, 4, 45)

# ╔═╡ b9ac03c2-d93e-4372-a3cb-4568baf015bb
myfields = fieldnames(typeof(w))

# ╔═╡ 522a7a5c-8400-4e1b-9ebd-6d76b33cb87f
getfield(w, :a)

# ╔═╡ d70384a3-51fc-4f71-a3d6-d14b03719495
x

# ╔═╡ ab382c7f-933e-4611-af23-fe31967ff8cc
[getfield(x, name) for name in myfields]

# ╔═╡ 8df1bf08-5ee2-4cfe-a273-62c71dc9b8fb
abstract type AbstractDual{S<:Real,T<:Real}

struct Dual5{S<:Real, T<:Real} <: AbstractDual{S,T}
    a::S
    b::T
end

# ╔═╡ babed906-7888-431b-8b11-d6e5889073fa
getfield.(myfields)

# ╔═╡ d7e3e2b2-6f3f-4761-9141-924ae1448343

# ╔═╡ ce4bd2c6-6269-47e3-8f43-50969a19487a
x::Dual ^ n::Integer = ..

# ╔═╡ 555aa973-a5ce-41f9-aab7-261a697758cc
*(f::Dual, g::Dual) = Dual(f.value * g.value,
                         f.value * g.deriv + f.deriv * g.value)

# Alternatively and equivalently
# *(f::Dual, g::Dual) = Dual(value(f) * value(g),
#                          value(f) * deriv(g) + deriv(f) * value(g))

# ╔═╡ eea59dfb-acc7-4d56-82b2-59f0ed40afc7
import Base: +, *

# ╔═╡ b0c794fb-0d77-4667-9d47-0bcdfaa8fecb
begin
	import Base: ^
	^(x::Dual, n::Integer) = Base.power_by_squaring(x, n)
end

# ╔═╡ Cell order:
# ╟─cac0bc38-ead9-4f99-a939-fc2645874efd
# ╠═ecb2b8d9-ceda-4008-8137-b29bd7ef3ba9
# ╟─b4f0668c-8b25-4588-9176-52bc7bfe797b
# ╟─0ece0b55-ede7-4db1-9b48-61e1d4534761
# ╟─6dcc01fe-c551-43d2-930b-ff1c93d15077
# ╠═0d352bcc-dfc1-429f-ad6c-40dc8c4cd038
# ╠═81863f06-8dac-48fd-a6e1-cf5df30c63f6
# ╠═3de40205-6b5d-4d99-b886-1a1dc6a07955
# ╠═d4c280cf-910c-4bbf-8dca-b01b14711296
# ╠═7082e444-6268-4721-a278-2511a22daf68
# ╠═43538919-4349-43f1-a196-759fe65d3be3
# ╟─e0c1df0c-ff62-4982-8f3d-3ccbb44b8d2f
# ╟─e6281275-0b87-47f4-a081-bb366ced3466
# ╟─6056ff4a-858f-447c-85cc-0f3f6e0e53ca
# ╟─db17c624-91b6-44b1-bd10-4d72f4cb9776
# ╟─b3b15bf5-115d-4831-9115-36c88e9b899c
# ╟─70843c81-6abe-40b2-8d0b-8d779fccedc8
# ╠═82c79506-68ea-41dd-888d-840371c08430
# ╟─291315e0-ad0e-4a57-b100-2baa5af04c04
# ╟─f2421fbf-5cfc-4332-be60-b61fa4e02d8c
# ╠═710612bd-4856-4e51-b2d6-2bdc972f196d
# ╟─70eb5167-6288-4f84-b140-21b4b329cb45
# ╠═9fb42f3c-c7ff-4f31-b5e6-9e96de16bac5
# ╠═07808ac8-c910-4e92-ac34-c9ff225f4be8
# ╟─8777c6a7-90e3-4dc7-ab61-3d12531f2ea9
# ╠═5488c55f-8d0a-4c9b-8f8b-40ed4fe14e67
# ╠═4f5e5eea-13df-477e-baea-c37ab6074aea
# ╠═da932afe-990f-4ca7-aa81-9d85e1ca3bac
# ╠═fac29cb2-dc57-4a06-9728-567169ac3768
# ╠═4118613c-e58a-4f63-81d1-0e197e66cf98
# ╠═ac8fa4de-d8ef-44f2-b709-aceb1e0ecb04
# ╠═92aa5678-478c-4cb6-9187-3285534aac8c
# ╠═04a21e57-a45c-4e4b-977a-42c6ee228b4a
# ╠═915bfb27-7234-43ea-a052-093f89d7bdcd
# ╟─9b213fb9-382a-4c5b-93fe-22252d29d0c6
# ╟─8eed2b5c-dde2-437e-8383-4ed206ce72de
# ╠═eea59dfb-acc7-4d56-82b2-59f0ed40afc7
# ╠═62a97a4b-dfa6-404d-958b-b7b4ab7d25ca
# ╠═9b0d361a-8c30-4ec0-a183-ac3890895f0e
# ╠═9b8f41b5-9ade-4f34-8d09-cb548bd6c242
# ╠═5b29f183-bda6-4ffb-bb78-5199980039d3
# ╟─fc624260-ba8a-4057-8a1f-2659e2343654
# ╠═abcbc0a1-33e7-430f-b0a8-e8b65c8666b2
# ╠═17498a94-5db9-4c5e-9268-e588fddeae54
# ╠═c847f0c2-1eac-4900-b80c-7f89846b2924
# ╠═6caea9da-3edc-4309-bad8-f348693960b2
# ╠═17cce3ab-1b4f-402d-bcc3-c7905360a441
# ╠═526a204c-6a72-4cde-900a-568667e45c64
# ╟─86aa82b1-caae-4124-b48a-14497420f471
# ╟─e359dd62-65f0-4033-98c2-aff400ed0bb7
# ╠═86141559-f6f4-43a6-8291-550adf31b25b
# ╠═cee6af83-0079-4acd-9e14-cd7b7eb6947f
# ╠═555aa973-a5ce-41f9-aab7-261a697758cc
# ╠═4e222c6a-054b-41e6-bef4-b0726344ee94
# ╟─a1b1b945-ba2e-4668-8698-4df02660d379
# ╠═8857e0aa-da38-4ae6-8ef9-efa66dd51cfc
# ╠═ad671591-a625-4de2-9b18-84c7773cb4b8
# ╠═053e90a7-5fa6-4063-9634-63d48ee16690
# ╠═d521f0e3-1636-4e8b-b8f3-639a90780690
# ╟─b7b81d9d-354d-4ec0-b50c-302dbcb8f0c4
# ╠═d67fa8a9-f43e-4eca-bb0e-1603c9f7da5d
# ╟─3bc3d5a3-9b9f-4d07-9d59-f24d5c005038
# ╟─fd34be78-af76-4d46-bdbd-19fd19155785
# ╟─caccd695-118c-46a8-b59f-376e4cbc83ce
# ╠═6418d0a4-8dd3-4816-8e5b-092cb769ebb6
# ╠═8ebb42e9-11ec-434d-87c9-4259c26f6526
# ╠═84b6a870-4515-4e6c-b306-d6b6d2a0fb8d
# ╠═a0865cb0-7077-44ff-a130-eb73cbbc3179
# ╠═f1fc7637-4804-4fc1-81ac-14874e66f9d6
# ╠═ded7adc0-619f-4735-bd93-9c26c4c21858
# ╠═1284374b-ab26-46b8-a74e-119655f5112e
# ╠═420c0ef0-3f7e-4044-be30-c2873636c6e8
# ╠═f426bed2-05f1-4598-9817-bb90262911ef
# ╠═09bfa122-ca5b-4b4a-a9eb-09b1b5482a76
# ╠═d08a434c-c816-4f43-9c93-2d0eda57723c
# ╠═b0c794fb-0d77-4667-9d47-0bcdfaa8fecb
# ╠═4e06b4f6-9753-4847-8c6c-120208520dc5
# ╠═131373df-8e58-49b0-be5f-3cd8d8deb61d
# ╠═76a7a87b-0283-4a85-9d88-0eec06ef689f
# ╠═d608bc38-5653-4ad5-b6e4-6b2236299e4c
# ╠═ce4bd2c6-6269-47e3-8f43-50969a19487a
# ╟─bf697142-0eaf-4717-9434-eaa89946beb4
# ╟─108212df-c86f-4b75-8576-65f3964687e9
# ╟─eb94cc27-ac9e-4202-acdc-248809eebdb9
# ╠═d7f6c80a-4922-42be-9bce-e493e603c373
# ╟─70ec3ada-f7ae-4526-be35-5aa8cd4533ca
# ╟─5011cf98-128e-4c94-bae0-109016427c50
# ╟─d57937d8-3f1e-4f85-bf6b-afc8cfe3e99d
# ╠═df1c033c-8a99-4b45-bc9c-3287a3a959ec
# ╟─b0902a74-01e0-456a-a027-a5030d494d7f
# ╟─9226e5a3-4445-4420-9b54-be31468e6b07
# ╠═ccc258fa-c928-4175-a48e-bd72b7d22f31
# ╠═c9f3dde8-7989-45e8-9831-894d82e57e90
# ╠═2d3876d6-5fcf-4bfb-8dba-4bdf252c73fc
# ╠═a0d7a33a-e067-4e69-9c70-14005f41a2c6
# ╠═812149d8-0040-4d67-b2a2-f1304af4b225
# ╟─75fa60fe-6583-4a6a-b6b0-06fb81006218
# ╟─4146cf8f-ab2f-457d-97ef-e6d39e93246e
# ╟─3e61a354-ca21-410b-b422-92ea362b8c10
# ╠═907f65aa-e1a1-4838-a1a0-cdaa0a1cfa45
# ╠═fc5839ce-935a-46aa-bbe7-01f203b0fdb5
# ╠═c445bab5-c012-4890-89bc-e7d6950fa7b4
# ╠═9008cc1d-e7c2-4b24-ad08-968c085baa36
# ╠═d5cb6844-6f9d-422e-a241-96450aff70f5
# ╠═7b2f03b4-d3dd-4fa6-a569-c878fc9f6557
# ╠═97097783-2948-4e9c-b769-8085c7943227
# ╠═4bdb5bfc-82c2-4de5-982a-4a69872749d1
# ╠═9c05cca8-76d9-4abd-8774-28e7f41d3b09
# ╠═05db96f7-3eeb-47aa-a4dc-dd67b0f02a2b
# ╠═b490d647-85f5-4b82-ab29-2d18ce1c7db4
# ╠═d981b89c-4833-4e38-9f80-246ed5f09445
# ╠═9bfe7470-afec-4a24-9f96-43272e2f9ecd
# ╠═61cc131f-2d02-4de1-b79a-9105d6decc15
# ╠═a5c1c463-2393-499d-852c-16e69549819a
# ╠═749ba068-6deb-4a7d-a650-412e5fbdb670
# ╠═226f5af6-702f-431c-a269-5dd4cafa2699
# ╠═c3831c37-8e09-4846-963d-27e1d4a824ac
# ╠═fd8bf8f8-68dd-4925-aa9b-8eca1937c577
# ╠═e3d64165-170e-4ee0-a9bc-ef759ce6afd7
# ╠═713dc49e-6038-40b7-a917-299277f7ea7a
# ╠═68d54b7c-4003-49e8-b2f9-666b882a5737
# ╠═7ae57ebb-d6bc-4b9b-913a-f9c4bb11d380
# ╠═99dd8229-8658-4efe-9133-3eb8cd00b6b3
# ╠═ef28ca54-38dd-4e8e-a5e8-897c802f81e1
# ╠═a31cfb26-9914-441c-aea8-6f341f2635fb
# ╠═f44efbd4-7ce2-44db-ac0e-eeb4aeafef70
# ╠═b05dcbc4-0149-48a3-8533-1bd6515d0b32
# ╠═7df9a91e-01d4-471f-9ad5-55ac56c245d3
# ╠═ef1c3d05-7020-4339-b431-63ad63b244cc
# ╠═186148a4-698e-460b-82dd-7b203698de45
# ╠═b9ac03c2-d93e-4372-a3cb-4568baf015bb
# ╠═522a7a5c-8400-4e1b-9ebd-6d76b33cb87f
# ╠═d70384a3-51fc-4f71-a3d6-d14b03719495
# ╠═ab382c7f-933e-4611-af23-fe31967ff8cc
# ╠═8df1bf08-5ee2-4cfe-a273-62c71dc9b8fb
# ╠═babed906-7888-431b-8b11-d6e5889073fa
# ╠═d7e3e2b2-6f3f-4761-9141-924ae1448343

max_stars_repo_path: 2. Uses of types.jl

max_issues_repo_path: 2. Uses of types.jl

max_forks_repo_path: 2. Uses of types.jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ContainerStateRunning <: SwaggerModel
    startedAt::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: startedAt

    function IoK8sApiCoreV1ContainerStateRunning(;startedAt=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ContainerStateRunning, Symbol("startedAt"), startedAt)
        setfield!(o, Symbol("startedAt"), startedAt)
        o
    end
end # type IoK8sApiCoreV1ContainerStateRunning

const _property_map_IoK8sApiCoreV1ContainerStateRunning = Dict{Symbol,Symbol}(Symbol("startedAt")=>Symbol("startedAt"))
const _property_types_IoK8sApiCoreV1ContainerStateRunning = Dict{Symbol,String}(Symbol("startedAt")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApiCoreV1ContainerStateRunning }) = collect(keys(_property_map_IoK8sApiCoreV1ContainerStateRunning))
Swagger.property_type(::Type{ IoK8sApiCoreV1ContainerStateRunning }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ContainerStateRunning[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ContainerStateRunning }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ContainerStateRunning[property_name]

function check_required(o::IoK8sApiCoreV1ContainerStateRunning)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ContainerStateRunning }, name::Symbol, val)
end


text:
using SemiLagrangian:
    splititr,
    splitvec,
    transposition,
    totuple,
    tovector,
    tupleshape,
    getextarray,
    CircEdge,
    AbstractInterpolation,
    interpolatemod!

@testset "split util" begin
    data = [
        (3, 15, [1:5, 6:10, 11:15])
        (3, 11, [1:4, 5:8, 9:11], 13)
        (5, 24, [1:5, 6:10, 11:15, 16:20, 21:24])
    ]
    for d in data
        @test d[3] == splititr(d[1], d[2])
    end
    data2 = [
        (47, 3, 15, [47:51, 52:56, 57:61])
        (34, 4, 5, [34:35, 36:36, 37:37, 38:38])
        (23, 7, 13, [23:24, 25:26, 27:28, 29:30, 31:32, 33:34, 35:35])
    ]
    for d in data2
        @test d[4] == splitvec(d[2], collect(d[1]:d[1]+d[3]-1))
    end
end
@testset "test tools" begin

    v = collect(1:53)
    t = splitvec(5, v)
    @test t[1] == collect(1:11)
    @test t[2] == collect(12:22)
    @test t[3] == collect(23:33)
    @test t[4] == collect(34:43)
    @test t[5] == collect(44:53)

    @test transposition(1, 2, 5) == [2, 1, 3, 4, 5]
    @test transposition(4, 2, 7) == [1, 4, 3, 2, 5, 6, 7]

end

@testset "transposition" begin
    @test transposition(1, 2, 5) == [2, 1, 3, 4, 5]
    @test transposition(1, 2, 5) == [2, 1, 3, 4, 5]
    @test transposition(3, 7, 7) == [1, 2, 7, 4, 5, 6, 3]
    @test transposition(7, 3, 7) == [1, 2, 7, 4, 5, 6, 3]
end

@testset "totuple tovector" begin
    @test (3, 4, 1) == totuple([3, 4, 1])
    @test [5, 1, 9, 2] == tovector((5, 1, 9, 2))
end

@testset "tupleshape" begin
    @test (1, 1, 71, 1) == tupleshape(3, 4, 71)
    @test (53, 1) == tupleshape(1, 2, 53)
    @test (1, 1, 67) == tupleshape(3, 3, 67)
end

function test_extarray(T, sz, decbeg, decend)
    tabor = rand(T, sz)
    tabext = getextarray(tabor, decbeg, decend)
    for ind in CartesianIndices(tabext)
        @test tabext[ind] ==
              tabor[CartesianIndex(mod.(ind.I .- decbeg .- 1, sz) .+ 1)]
    end
    @test size(tabext) == sz .+ decbeg .+ decend
end
@testset "ExtArray test" begin
    test_extarray(Float64, (20, 10), (3, 5), (2, 7))
    test_extarray(Float64, (4, 10, 7), (2, 3, 5), (1, 2, 4))
    test_extarray(Float64, (20,), (3,), (7,))
end


text:
######################################################################
# backend.jl
#
# Currently uses PyCall and PyNEC as the backend
######################################################################

using PyCall
@pyimport PyNEC as pn

function geo_card(geo::PyObject, w::Wire, tag_id::Int)
    geo[:wire](tag_id, w.n_segments, w.a..., w.b..., w.rad, w.rdel, w.rrad)
end

function ex_card(context::PyObject, vs::VoltageSource, s)

    # loop over wires to find matching segment
    tag_id = 0
    for (sei, se) in enumerate(s)
        se == vs.w && (tag_id = sei)
    end

    context[:ex_card](0,            # 0 for voltage source
                      tag_id,       # I2 - tag number of source segment
                      vs.m,         # I3 - m
                      0,            # I4 -
                      real(vs.v),   # F1 - Real part of voltage in volts
                      imag(vs.v),   # F2 - Imaginary part of voltage
                      0,            # F3 -
                      0,            # F4 -
                      0,            # F5 -
                      0             # F6 -
                     )
end

function ex_card(context::PyObject, lw::LinearWave, s)
    context[:ex_card](1,            # I1 - 1 for linear wave
                      lw.n_thetas,  # I2 - number of theta angles
                      lw.n_phis,    # I3 - number of phi angles
                      0,            # I4 - not sure what this is
                      lw.theta,     # F1
                      lw.phi,       # F2
                      lw.eta,       # F3
                      lw.d_theta,   # F4
                      lw.d_phi,     # F5
                      0             # F6 - min/maj axis
                     )
end

function RadiationPattern(s, ex::Excitation, fmhz::Real;
                          thetas::RealOrRange=0:180,
                          phis::RealOrRange=0:360
                         )
    # create context
    context = pn.nec_context()

    # put structure in
    geo = context[:get_geometry]()
    for (sei, se) in enumerate(s)
        geo_card(geo, se, sei)
    end
    context[:geometry_complete](0)

    # ground
    # TODO: obviously, don't just do it this way
    # nullifies previous parameters and sets free-space condition
    context[:gn_card](-1, 0, 0, 0, 0, 0, 0, 0)

    # excitation
    ex_card(context, ex, s)

    # frequency
    #context[:fr_card](0, 2, 2400, 100.0e6)
    context[:fr_card](0, 1, fmhz, 0)

    # rp_card
    context[:rp_card](0,                    # I1 - 0 = normal mode
                      length(thetas),       # I2 - num theta values
                      length(phis),         # I3 - num phi values
                      0,                # I4-X - output format
                      0,                # I4-N - normalization
                      0,                # I4-D - powergain or directive gain
                      0,                # I4-A - Averaging?
                      istart(thetas),       # F1 - initial theta angle
                      istart(phis),         # F2 - initial phi angle
                      step(thetas),         # F3 - theta increment
                      step(phis),           # F4 - phi increment
                      0.0,              # F5 - Radial distance
                      0.0               # F6 - Gain normalization factor
                     )

    # get radiation pattern
    rp = context[:get_radiation_pattern](0)
    return RadiationPattern(rp)
end

function RadiationPattern(rp::PyObject)
    gains_db = rp[:get_gain]()
    gains = 10.0 .^ (gains_db / 10.0)
    thetas = rp[:get_theta_angles]()
    phis = rp[:get_phi_angles]()
    return RadiationPattern(gains_db, gains, thetas, phis)
end


text:
using SpecialFunctions: loggamma
using Printf: @printf

T = BigFloat
case = 1
if case == 1
    α_vals = T[ n/10 for n = 1:2:9 ]
    x_vals = T[ n*5 for n = 1:8 ]
    β = 13 / parse(T, "17")
elseif case == 2
    α_vals = T[ (994+n)/1000 for n = 1:5 ]
    x_vals = T[ n*5 for n = 5:10 ]
    β = one(T)
end
@printf("\nβ = %10.6f\n", β)

tol = parse(T, "10")^(-16)
include("assess.jl")

@printf("\nMinimum N such that Nth term is less than %0.2e\n\n", tol)
@printf("  x    ")
for α in α_vals
    @printf("   α = %5.3f  ", α)
end
@printf("\n\n")
for x in x_vals
    @printf("%5.2f | ", x)
    for α in α_vals
        Nmin, term = assess(α, β, x, tol)
        if term < tol
            @printf("%3d %8.2e |", Nmin, term)
        else
            @printf("%3s %8s |", "**", "**")
        end
    end
    @printf("\n")
end
@printf("\n\nOptimal truncation at N = x^(1/α):\n")
@printf("\n  x    ")
for α in α_vals
    @printf("   α = %5.3f  ", α)
end
@printf("\n\n")
for x in x_vals
    @printf("%5.2f | ", x)
    for α in α_vals
        Nopt = ceil(Integer, x^(1/α))
        if Nopt > 10_000_000_000
            @printf(" %11s |", "**")
        else
            @printf(" %11d |", Nopt)
        end
    end
    @printf("\n")
end


text:
using BasisFunctions, DomainSets, Test
@testset "Spaces " begin
s = GenericFunctionSpace{Float64,ComplexF64}()
@test typeof(zero(s)(rand()))==ComplexF64
@test typeof(one(s)(rand()))==ComplexF64
@test zero(s)(rand())==0
@test one(s)(rand())==1

s = MeasureSpace(FourierWeight())
@test measure(s) ==FourierWeight()
@test domain(s) == UnitInterval()
@test space(FourierWeight()) == s
@test FourierSpace() == s
@test ChebyshevTSpace() == MeasureSpace(ChebyshevTWeight{Float64}())
@test Inf ∈ domain(L2())

s1 = GenericFunctionSpace{Float64,ComplexF64}()
s2 = MeasureSpace(FourierWeight())
@test components(tensorproduct(s1,s2)) == (s1,s2)
end


text:
# This file includes and example file for HInterp2D interpolation

using FractalTools
using GeometryBasics
using Makie

# Generate data
f(x, y) = [
    x^2 + y^2 + 1,
    x^2 - y^2
    ]
ngon = Triangle(
    Point(BigFloat(0.), BigFloat(0.)),
    Point(BigFloat(1.), BigFloat(0.)),
    Point(BigFloat(0.5), BigFloat(1.)))
npts = 100
pts = getdata(f, ngon, npts)
interp = interpolate(pts, HInterp2D(0.01 * ones(2,2)))

function err(x, y)
    fval = f(x, y)
    ival = interp(x, y)
    map(item -> abs.(item), fval - ival) ./ map(item -> abs.(item), fval) * 100
end

tpts = getdata(ngon, npts)

fig, ax, plt = trisurf(tpts, (x, y) -> f(x, y)[1], meshcolor3=:red)
trisurf!(ax, tpts, (x, y) -> interp(x, y)[1], meshcolor3=:blue)
fig, ax, plt = trisurf(tpts, (x, y) -> err(x, y)[1], meshcolor3=:red)