text:
using Test

using CSV
using DataFrames
using FlightMechanicsUtils
using OrdinaryDiffEq
using StaticArrays

using FlightMechanicsSimulator

# ---------- PROPAGATE COORDINATED TURN ----------
# Check that this is a trimmed condition
# Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
# and simulation: dynamics, controls design, and autonomous systems. John Wiley
# & Sons.
# Example 3.6-2 (page 191)

# Trimmed conditions: page 192
# C     X(1)  -> vt (m/s)
# C     X(2)  -> alpha (rad)
# C     X(3)  -> beta (rad)
# C     X(4)  -> phi (rad)
# C     X(5)  -> theta (rad)
# C     X(6)  -> psi (rad)
# C     X(7)  -> P (rad/s)
# C     X(8)  -> Q (rad/s)
# C     X(9)  -> R (rad/s)
# C     X(10) -> North (m)
# C     X(11) -> East (m)
# C     X(12) -> Altitude (m)
# C     X(13) -> Pow
x_stev = [
    502.0 * FT2M,
    0.2392628,
    0.0005061803,
    1.366289,
    0.05000808,
    0.2340769,
    -0.01499617,
    0.2933811,
    0.06084932,
    0.0 * FT2M,
    0.0 * FT2M,
    0.0 * FT2M,
    64.12363,
]
controls_stev = [0.8349601, -1.481766, 0.09553108, -0.4118124]
xcg = 0.35
dssd, outputs = f(
    time,
    SixDOFAeroEuler(x_stev),
    controls_stev,
    F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
    F16StevensAtmosphere(x_stev[12]),
    LHDownGravity(FlightMechanicsSimulator.F16Stevens.GD*FT2M),
)

x_dot = get_xdot(dssd)

# Linear acceleration
@test isapprox(x_dot[1:3], zeros(3), atol = 5e-4)
# Angular acceleration
@test isapprox(x_dot[7:9], zeros(3), atol = 1e-5)
# Cost function
cost =
    x_dot[1]^2 +
    100 * (x_dot[2]^2 + x_dot[3]^2) +
    10 * (x_dot[7]^2 + x_dot[8]^2 + x_dot[9]^2)

@test isapprox(cost, 0.0, atol = 1e-5)

# Desired turn rate
@test isapprox(x_dot[6], 0.3, atol = 1e-6)
# No vertical changes
@test isapprox(x_dot[12], 0.0, atol = 1e-4)

# RETRIM to refine flying condition
dssd, controls_trim, outputs_trim, cost = trim(
    SixDOFAeroEuler(x_stev),
    controls_stev,
    F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
    F16StevensAtmosphere(x_stev[12]),
    LHDownGravity(FlightMechanicsSimulator.F16Stevens.GD*FT2M),
    0.0,
    0.3,
)

x_trim = get_x(dssd)
x_dot_trim = get_xdot(dssd)

# Linear acceleration
@test isapprox(x_dot_trim[1:3], zeros(3), atol=1e-12)
# Angular acceleration
@test isapprox(x_dot_trim[7:9], zeros(3), atol=1e-12)
# Desired turn rate
@test isapprox(x_dot_trim[6], 0.3, atol = 1e-6)
# No vertical changes
@test isapprox(x_dot_trim[12], 0.0, atol = 1e-4)

dt = 0.01  # s
t0 = 0.0  # s
t1 = 180.0  # s

x = x_trim

results = simulate(
    t0,
    t1,
    SixDOFAeroEuler(x),
    controls_trim,
    F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
    F16StevensAtmosphere,
    LHDownGravity(FlightMechanicsSimulator.F16Stevens.GD*FT2M);
    solver=RK4(),
    solve_args=Dict(:reltol=>1e-10, :saveat=>dt),
    )

# Check X, Y against Stevens
# Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
# and simulation: dynamics, controls design, and autonomous systems. John Wiley
# & Sons.
# Example 3.6-6 (page 198)
xy_trajectory_data = [
    # time (s)  X (ft)  Y(ft)
    0.00e0      0.00e0    0.00e0
    1.00e1      2.36e2    3.33e3
    2.00e1     -4.68e2    6.65e1
    3.00e1      6.90e2    3.20e3
    4.00e1     -8.97e2    2.61e2
    5.00e1      1.09e3    2.94e3
    6.00e1     -1.26e3    5.68e2
    7.00e1      1.40e3    2.59e3
    8.00e1     -1.51e3    9.62e2
    9.00e1      1.60e3    2.16e3
    1.00e2     -1.66e3    1.41e3
    1.10e2      1.67e3    1.70e3
    1.20e2     -1.66e3    1.89e3
    1.30e2      1.61e3    1.22e3
    1.40e2     -1.53e3    2.35e3
    1.50e2      1.42e3    7.87e2
    1.60e2     -1.28e3    2.76e3
    1.70e2      1.11e3    4.22e2
    1.80e2     -9.21e2    3.07e3
 ]

 for case in eachrow(xy_trajectory_data)
    r = results[findall(in(case[1]), results.time), :]
    @test isapprox(case[2], r[1, :x] * M2FT, atol=20)
    @test isapprox(case[3], r[1, :y] * M2FT, atol=20)
 end

# Select last time step state
 x = results[end, :]

# Check that TAS, α, β, θ, ϕ, p, q, r, alt, pow remain constant
@test isapprox(x[:tas], x_trim[1])  # TAS
@test isapprox(x[:α], x_trim[2])  # α
@test isapprox(x[:β], x_trim[3])  # β
@test isapprox(x[:ϕ], x_trim[4])  # ϕ
@test isapprox(x[:θ], x_trim[5])  # θ
@test isapprox(x[:p], x_trim[7])  # p
@test isapprox(x[:q], x_trim[8])  # q
@test isapprox(x[:r], x_trim[9])  # r
@test isapprox(x[:z], x_trim[12], atol=1e-12)  # alt
@test isapprox(x[:pow], x_trim[13])  # pow


text:
using DataFrames, TidyStanza
using Test

include("relocate-test.jl")

df = DataFrame(x = repeat(1:3, inner = 2, outer = 2),
       a = repeat(4:6, inner = 2, outer = 2),
       b = repeat(7:9, inner = 2, outer = 2),
       val1 = ["ce_val1_1","cf_val1_1","ce_val1_2","cf_val1_2","ce_val1_3","cf_val1_3","de_val1_1","df_val1_1","de_val1_2","df_val1_2","de_val1_3","df_val1_3"],
       val2 = ["ce_val2_1","cf_val2_1","ce_val2_2","cf_val2_2","ce_val2_3","cf_val2_3","de_val2_1","df_val2_1","de_val2_2","df_val2_2","de_val2_3","df_val2_3"],
       cname1 = repeat(["c", "d"], inner = 6),
       cname2 = repeat(["e", "f"], 6)
)

wider_df1 = DataFrame(
    x = repeat(1:3, inner=2),
    a =  repeat(4:6, inner=2),
    b =  repeat(7:9, inner=2),
    cname2 = repeat(["e", "f"], 3),
    val1_c = ["ce_val1_1","cf_val1_1", "ce_val1_2", "cf_val1_2", "ce_val1_3", "cf_val1_3"],
    val1_d = ["de_val1_1","df_val1_1", "de_val1_2", "df_val1_2", "de_val1_3", "df_val1_3"],
    val2_c = ["ce_val2_1","cf_val2_1", "ce_val2_2", "cf_val2_2", "ce_val2_3", "cf_val2_3"],
    val2_d = ["de_val2_1","df_val2_1", "de_val2_2", "df_val2_2", "de_val2_3", "df_val2_3"]
)

wider_df2 = DataFrame(
    x = repeat(1:3, inner = 4),
    a = repeat(4:6, inner = 4),
    b = repeat(7:9, inner = 4),
    val2 = ["ce_val2_1","cf_val2_1","de_val2_1","df_val2_1","ce_val2_2","cf_val2_2","de_val2_2","df_val2_2","ce_val2_3","cf_val2_3","de_val2_3","df_val2_3"],
    val1_c_e = vcat("ce_val1_1", repeat([missing],3), "ce_val1_2",repeat([missing],3),  "ce_val1_3", repeat([missing],3) ),
    val1_c_f = vcat(repeat([missing],1), "cf_val1_1", repeat([missing],3), "cf_val1_2",repeat([missing],3),  "cf_val1_3", repeat([missing],2) ),
    val1_d_e = vcat(repeat([missing],2), "de_val1_1", repeat([missing],3), "de_val1_2",repeat([missing],3),  "de_val1_3", repeat([missing],1) ),
    val1_d_f = vcat(repeat([missing],3), "df_val1_1", repeat([missing],3), "df_val1_2",repeat([missing],3),  "df_val1_3")
)

wider_df3 = DataFrame(
    x = 1:3,
    a = 4:6,
    b = 7:9,
    val1_c_e = ["ce_val1_1", "ce_val1_2", "ce_val1_3"],
    val1_c_f = ["cf_val1_1", "cf_val1_2", "cf_val1_3"],
    val1_d_e = ["de_val1_1", "de_val1_2", "de_val1_3"],
    val1_d_f = ["df_val1_1", "df_val1_2", "df_val1_3"],
    val2_c_e = ["ce_val2_1", "ce_val2_2", "ce_val2_3"],
    val2_c_f = ["cf_val2_1", "cf_val2_2", "cf_val2_3"],
    val2_d_e = ["de_val2_1", "de_val2_2", "de_val2_3"],
    val2_d_f = ["df_val2_1", "df_val2_2", "df_val2_3"]
)

df_noid = DataFrame(t = [:a, :b, :c, :a, :b, :c], x = 1:6, y = 11:16)
df_noid_res = DataFrame(
    x_a = [1 , missing, missing, 4, missing, missing],
    x_b = [missing, 2, missing, missing, 5, missing],
    x_c = [missing, missing, 3, missing, missing, 6],
    y_a = [11 , missing, missing, 14, missing, missing],
    y_b = [missing, 12, missing, missing, 15, missing],
    y_c = [missing, missing, 13, missing, missing, 16],
)

@testset "TidyStanza.jl" begin
    # Write your tests here.
    @test pivot_wider(df, names_from = :cname1, values_from = [:val1,:val2]) == wider_df1
    @test isequal(pivot_wider(df, names_from = [:cname1, :cname2], values_from = :val1),wider_df2)
    @test pivot_wider(df, names_from = [:cname1,:cname2], values_from = [:val1,:val2]) == wider_df3
    @test pivot_wider(df, names_from = "cname1", values_from = ["val1","val2"]) == wider_df1
    @test isequal(pivot_wider(df, names_from = ["cname1", "cname2"], values_from = "val1"),wider_df2)
    @test pivot_wider(df, names_from = ["cname1", "cname2"], values_from = ["val1","val2"]) == wider_df3
    @test isequal(pivot_wider(df_noid, names_from = [:t], values_from = [:x, :y]), df_noid_res)
end

ext: jl


text:
# using GivensBoson
using LinearAlgebra
include("./GivensBoson.jl")
using Main.GivensBoson

Nrange = 4:40:300
result = zeros(size(Nrange)[1],4)
for id in CartesianIndices(Nrange)
    N = Nrange[id]
    ϵ = 1E-11
    A = 1
    λ = 1
    bc = "PBC"
    begin
        Q = zeros(N,N)
        for id in CartesianIndices((1:N,1:N))
            if id[1]-id[2]==1
                Q[id] = 1.0*A
            elseif id[2] - id[1] ==1
                Q[id] = -1.0*A
            end
        end

        if bc == "PBC"
            Q[N,1] = -1.0*A
            Q[1,N] =  1.0*A
        end
        Q = Q/2.0
        Q
    end
    ϵ = 1E-10
    Lambda = diagm([λ for i = 1:N])
    M = [Lambda Q;transpose(Q) Lambda]
    A = M
    origin_A = copy(A)
    # S,V = recanonicalize(A,"AntiSymmetry")
    S,V = given_eigen_solver(A,zeromode=true)

    η = diagm(vcat([1.0 for i=1:N],[-1.0 for i=1:N]))
    result[id,1]=N
    if N%4 == 0
        t = transpose(V)*η*V
        print("Number of zero mode is: ",sum(abs.(t-η)),'\n')
        for i = 1:2N
            if abs(t[i,i])<1E-11
                t[i,i] = i>N ? -1.0 : 1.0
            end
        end
    end
    # display(t)
    result[id,2]=norm(t-η)
    result[id,3]=norm(S-diagm(diag(S)))
    result[id,4]=norm(transpose(V)*origin_A*V-S)
    # if N == 12
        # display(V)
        # display(transpose(V)*origin_A*V-S)
    # end
end
vscodedisplay(result)
# display(result)


text:
mutable struct Cell <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    status::Bool
end

"""
``` julia
game_of_life(;
    rules::Tuple = (2, 3, 3, 3),
    dims = (100, 100),
    Moore = true
)
```
Same as in [Conway's game of life](@ref).
"""
function game_of_life(;
    rules::Tuple = (2, 3, 3, 3),
    dims = (100, 100),
    Moore = true
)
    space = GridSpace(dims, moore = Moore)
    properties = Dict(:rules => rules)
    model = ABM(Cell, space; properties = properties)
    node_idx = 1
    for x in 1:dims[1]
        for y in 1:dims[2]
            add_agent_pos!(Cell(node_idx, (x, y), false), model)
            node_idx += 1
        end
    end
    return model, dummystep, game_of_life_model_step!
end

function game_of_life_model_step!(model)
    new_status = fill(false, nagents(model))
    for (agid, ag) in model.agents
        nlive = nlive_neighbors(ag, model)
        if ag.status == true && (nlive ≤ model.rules[4] && nlive ≥ model.rules[1])
            new_status[agid] = true
        elseif ag.status == false && (nlive ≥ model.rules[3] && nlive ≤ model.rules[4])
            new_status[agid] = true
        end
    end

    for k in keys(model.agents)
        model.agents[k].status = new_status[k]
    end
end

function nlive_neighbors(ag, model)
    neighbors_coords = node_neighbors(ag, model)
    nlive = 0
    for nc in neighbors_coords
        nag = model.agents[Agents.coord2vertex((nc[2], nc[1]), model)]
        if nag.status == true
            nlive += 1
        end
    end
    return nlive
end


text:
using Sylvia
using Test

using Sylvia: commuteswith, isfalse, istrue, sym

@sym obj
@sym [Number] a b c d
@sym [Bool] x y z
@sym [Matrix{Float64}] A B
@sym [Vector{Float64}] v w

@testset "identities" begin
    @test obj == sym(:obj)
    @test a == sym(Number, :a)
    @test a == S"a::Number" == @!(a::Number) == @! :a::Number
    @test a + b == S"(a + b)::Number" == @!((a + b)::Number) == @! :a::Number + :b::Number
    @test S"Float64" == @! Float64
    @test S"AbstractMatrix" == @! AbstractMatrix

    @test iszero(zero(a))
    @test iszero(zero(typeof(a)))
    @test isone(one(a))
    @test isone(one(typeof(a)))
    @test iszero(zero(S"Float64"))
    @test isone(one(S"Float64"))
    @test a[1] == getindex(a, 1)
    @test a[:b] == getindex(a, S":(:b)")
    @test a.b == getproperty(a, :b)
    @test a.b == getproperty(a, S":(:b)")
    @test a(1) == a(S"1") == Sylvia.apply(a, 1)
end

@testset "conversion" begin
    @sym [Float64] q r
    @test convert(Tuple{Sylvia.Sym{Float64}, Sylvia.Sym{Float64}}, (:q, :r)) == (q, r)
    @test convert(Tuple{Sylvia.Sym, Sylvia.Sym}, (:a, :b)) == (a, b)
end

@testset "promotion" begin
    for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
        @test @!(a::T) + @!(b::T) isa Sylvia.Sym{T}
        @test @!(a::T) - @!(b::T) isa Sylvia.Sym{T}
        @test @!(a::T) * @!(b::T) isa Sylvia.Sym{T}
        @test @!(a::T) / @!(b::T) isa Sylvia.Sym{T}
        @test cos(@! a::T) isa Sylvia.Sym{T}
        @test sin(@! a::T) isa Sylvia.Sym{T}
        @test exp(@! a::T) isa Sylvia.Sym{T}
        @test log(@! a::T) isa Sylvia.Sym{T}
        @test (@! a::T)^2 isa Sylvia.Sym{T}
        @test (@! a::T)^-1 isa Sylvia.Sym{T}
    end
end

@testset "wild" begin
    @sym [Wild] wild_default [Wild{Number}] wild_number [Wild{Float64}] wild_float
    @test tagof(wild_default) == Sylvia.DEFAULT_TAG
    @test tagof(wild_number) == Number
    @test tagof(wild_float) == Float64
    @test !Sylvia.iswild(a)
    @test @! unwrap Sylvia.iswild(a::Wild{Number})
    @test Sylvia.iswild(wild_default)
    @test Sylvia.iswild(wild_number)
    @test Sylvia.iswild(wild_float)
    @test Sylvia.ismatch(Sylvia.match(a, wild_number))
    @test Sylvia.ismatch(Sylvia.match(wild_float, wild_default))
    @test !Sylvia.ismatch(Sylvia.match(wild_default, wild_float))
    @test Sylvia.ismatch(Sylvia.match(wild_number, wild_default))
    @test Sylvia.ismatch(Sylvia.match(@!(a::Float64), wild_float))
    @test !Sylvia.ismatch(Sylvia.match(a, wild_float))
end

@testset "mock" begin
    @sym [Mock] mock_default [Mock{Real}] mock_real [Mock{Float64}] mock_float
    function f(x::Number)
        y = x + 1
        return 2y
    end
    g(x::Number) = x

    @test Sylvia.ismocking(mock_default) && Sylvia.ismocking(mock_real) && Sylvia.ismocking(mock_float)

    @test_throws MethodError f(mock_default)
    @test_throws MethodError f(mock_real)
    @test_throws MethodError f(mock_float)
    @test_throws MethodError g(mock_default)
    @test_throws MethodError g(mock_real)
    @test_throws MethodError g(mock_float)

    Sylvia.@register f 1
    Sylvia.@register g 1

    @test f(mock_default) == Sylvia._apply(f, mock_default)
    @test f(mock_real) == 2(mock_real + 1)
    @test f(mock_float) == Sylvia._apply(f, mock_float)
    @test g(mock_default) == Sylvia._apply(g, mock_default)
    @test g(mock_real) == mock_real
    @test g(mock_float) == Sylvia._apply(g, mock_float)
end

@testset "default rules" begin
    @test commuteswith(+, a, b)
    @test commuteswith(+, A, B)
    @test commuteswith(+, A, b)
    @test commuteswith(+, a, B)

    @test commuteswith(*, a, b)
    @test commuteswith(*, a, B)
    @test commuteswith(*, A, b)
    @test commuteswith(*, A, B) isa Sylvia.Sym

    @test commuteswith(&, x, y)
    @test commuteswith(|, x, y)
end

@testset "dynamic rules" begin
    @scope begin
        @! set iszero(a) --> true
        @! set a in b --> true
        @scope begin
            @! set a in b --> false
            @test iszero(a)
            @test !(a in b)
            @! unset a in b
            @! set iszero(a) --> false
            @test a in b
            @test !iszero(a)
            @! clear!
            @test iszero(a)
        end
        @test iszero(a)
        @test a in b
    end
end

@testset "expr" begin
    @test Sylvia.expr(a) === :a
    @test Sylvia.expr(a + b) == :(a + b)
    @test Sylvia.expr((a, b, c, d)) == :((a, b, c, d))
    @test Sylvia.expr([a, b, c, d]) == :([a, b, c, d])
    @test Sylvia.getsymbols((a, b, c, d)) == [:a, :b, :c, :d]
    @test Sylvia.getsymbols([a, b, c, d]) == [:a, :b, :c, :d]
    @test Sylvia.getops((a + b, b * c, sin(d))) == [+, *, sin]
end

@testset "combine" begin
    @test a[] == Sylvia.combine(:ref, a)
    @test a[1] == Sylvia.combine(:ref, a, S"1")
    @test a[1, 2] == Sylvia.combine(:ref, a, S"1", S"2")
    @test a' == Sylvia.combine(Symbol("'"), a)
end

@testset "tags" begin
    @test Sylvia.tagof(S"1") === Int
    @test Sylvia.tagof(S"1.0") === Float64
    @test Sylvia.issubtag(Number, Any)
    @test Sylvia.issubtag(Real, Number)
    @test Sylvia.issubtag(:a, Any)
    @test Sylvia.issubtag(1, Any)
    @test Sylvia.issubtag(:a, :a)
    @test !Sylvia.issubtag(:a, Number)
    @test !Sylvia.issubtag(:a, :b)
end

@testset "simplification" begin
    @test gather(a + b + c + d + a + b + c + d) == 2a + 2b + 2c + 2d
    @scope begin
        @! set iszero(a) --> true
        @! set isone(b) --> true
        @test gather(a + b*c + c) == 2c
    end

    @scope begin
        @! set istrue(x) --> true
        @! set isfalse(y) --> true

        @test gather(x & y) == y
        @test gather(x | y) == x
        @test gather(x | y | z) == x
        @test gather((x | y) & z) == z
    end
end

@testset "substitution" begin
    @test substitute(a + b, a => b) == b + b
    @test substitute(a + b, a => A) == A + b
    @test substitute(a + b, a => b, b => a) == b + a
    @test substitute(a + b, a => b, strict=true) == b + b
    @test substitute(a + b, a => x, strict=true) == x + b
    @test_throws ErrorException substitute(a + b, a => A, strict=true)
end

@testset "circular queries" begin
    @! set :f(:x::Wild) --> :x::Wild
    @! set :f(:x::Wild + :y::Wild) --> :f(:f(:x::Wild) + :f(:y::Wild))
    @! set :hop1(:x::Wild) --> :hop2(:x::Wild)
    @! set :hop2(:x::Wild) --> :hop3(:x::Wild)
    @! set :hop3(:x::Wild) --> :hop1(:x::Wild)

    @test (@! :f(a + b)) == (@! raw :f(a + b))
    @test (@! :hop1(a)) == (@! raw :hop1(a))
    @test (@! :hop2(a)) == (@! raw :hop2(a))
    @test (@! :hop3(a)) == (@! raw :hop3(a))
end

@testset "arrays" begin
    @test all(Vector(A, 2) .== [A[1], A[2]])
    @test all(Matrix(A, 2, 2) .== [A[1,1] A[1,2]; A[2,1] A[2,2]])

    @! set length(v) --> 4

    @test collect(v) == Vector(v, 4)
    for (i, vi) in enumerate(v)
        @test v[i] == vi
    end
end

@testset "@!" begin
    function g end
    @test @!(g(a)) == S"g(a)"
    @! set g(a) --> a + 1
    @test a + 1 == @! resolve g(a)
end

# This function definition has to live at the top level
function h end

@testset "function generation" begin
    X = randn(2, 2)
    x = a + b + c
    f = @! eval (a, b, c) -> x
    @! eval :g(A) = A^2
    @! eval function h(a, b, c)
        x
    end

    @test f(1, 2, 3) == 6
    @test h(1, 2, 3) == 6
    @test g(X) ≈ X^2
end

@testset "type generation" begin
    @! eval struct :MyType <: :Number
        a
        b
    end
    instance = MyType(1, 2)

    @test MyType <: Number
    @test fieldnames(MyType) == (:a, :b)
    @test fieldtype(MyType, :a) == Number
    @test fieldtype(MyType, :b) == Number
    @test instance.a == 1 && instance.b == 2
end


text:
using Test

include("fibonacci.jl")

const fib = [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352, 24157817, 39088169, 63245986, 102334155, 165580141, 267914296, 433494437, 701408733, 1134903170, 1836311903, 2971215073, 4807526976, 7778742049, 12586269025]

@testset "Fib is defined" begin
    @test isdefined(@__MODULE__, :Fib)
    @test typeof(Fib(100)) == Fib
end

@testset "Can be iterated" begin
    @test [a for a in Fib(50)] == fib
end

@testset "Fib does not contain fields to store state internally" begin
    # Backport fieldtypes for Julia 1.0
    if VERSION < v"1.1"
        @eval fieldtypes(T::Type) = ntuple(i -> fieldtype(T, i), fieldcount(T))
    end

    # Ensure only one numeric value (n) can be stored
    # This should prevent Fib types that use a tuple or struct to store more state
    @test fieldcount(Fib) == 1
    @test fieldtypes(Fib)[1] <: Number
end

@testset "Can be collected" begin
    @test collect(Fib(50)) == fib
    @test Base.IteratorSize(Fib) == Base.HasLength()
end

@testset "Has an eltype" begin
    @test typeof(collect(Fib(50))) == Vector{Int}
    @test Base.IteratorEltype(Fib) == Base.HasEltype()
end


text:
using Pkg
Pkg.activate(@__DIR__)

using Random

# for bug with display
ENV["GKSwstype"] = "100"

example_files = filter(readdir(@__DIR__, join=true)) do f
    endswith(f, ".jl") && !occursin("large", f) && !occursin("result", f) && f != "activate.jl"
end

example_shuffle = randperm(length(example_files))

for file in example_files[example_shuffle[1:2]]
    @info "running example $file"
    run(`julia $file`)
end


text:
# This file was generated, do not modify it. # hide
r_depth = range(dtr_model, :max_depth, lower=2, upper=20)

tm = TunedModel(model=dtr_model, ranges=[r_depth], tuning=Grid(resolution=10),
                resampling=CV(nfolds=5, rng=1254), measure=rms)
mtm = machine(tm, X, y)

fit!(mtm, rows=train)

ypred = predict(mtm, rows=test)
round(rms(ypred, y[test]), sigdigits=3)


text:
# Force transform of input range for the following formats
const VIO_GRAY_SCALE_TYPES = Set((AV_PIX_FMT_GRAY8, AV_PIX_FMT_GRAY10LE))

const VIO_AVCOL_PRI_TO_SWS_CS = Dict{Cint, Cint}(
    AVCOL_PRI_BT709     => SWS_CS_ITU709,
    AVCOL_PRI_SMPTE170M => SWS_CS_ITU601,
    AVCOL_PRI_SMPTE240M => SWS_CS_SMPTE240M,
    AVCOL_PRI_BT2020    => SWS_CS_BT2020,
    AVCOL_PRI_BT470M    => SWS_CS_FCC
)

mutable struct SwsTransform
    sws_context::SwsContextPtr
    srcframe::AVFramePtr
    dstframe::AVFramePtr
end

function SwsTransform(src_w, src_h, src_pix_fmt, src_primaries, src_color_range,
                      dst_w, dst_h, dst_pix_fmt, dst_primaries, dst_color_range,
                      sws_color_options, sws_scale_options)
    sws_context = SwsContextPtr()
    sws_options = Dict(:srcw       => string(src_w),
                       :srch       => string(src_h),
                       :src_format => string(src_pix_fmt),
                       :dstw       => string(dst_w),
                       :dsth       => string(dst_h),
                       :dst_format => string(dst_pix_fmt))
    for (key, val) in pairs(sws_scale_options)
        sws_options[key] = string(val)
    end
    if sws_options[:srcw] != string(src_w) ||
        sws_options[:dstw] != string(dst_w) ||
        sws_options[:srch] != string(src_h) ||
        sws_options[:dsth] != string(dst_h)
        throw(ArgumentError("Changing pixel dimensions with sws_scale_options not yet supported"))
    end
    set_class_options(sws_context, sws_options)
    ret = sws_init_context(sws_context, C_NULL, C_NULL)
    ret < 0 && error("Could not initialize sws context")

    inv_table = _vio_primaries_to_sws_table(src_primaries)
    table = _vio_primaries_to_sws_table(dst_primaries)
    ret = sws_update_color_details(sws_context;
                                   inv_table = inv_table,
                                   src_range = src_color_range,
                                   table = table,
                                   dst_range = dst_color_range,
                                   sws_color_options...)
    ret || error("Could not set sws color details")

    srcframe = AVFramePtr()
    dstframe = AVFramePtr()
    SwsTransform(sws_context, srcframe, dstframe)
end

function SwsTransform(src_w, src_h, src_pix_fmt, src_primaries, src_color_range,
             dst_pix_fmt, dst_primaries, dst_color_range, sws_color_options,
                      sws_scale_options)
    SwsTransform(src_w, src_h, src_pix_fmt, src_primaries, src_color_range,
                 src_w, src_h, dst_pix_fmt, dst_primaries, dst_color_range,
                 sws_color_options, sws_scale_options)
end

mutable struct GrayTransform
    srcframe::AVFramePtr
    src_depth::Int # Does not include padding
    dstframe::AVFramePtr
    dst_depth::Int # Does not include padding
end

function GrayTransform()
    srcframe = AVFramePtr()
    dstframe = AVFramePtr()
    GrayTransform(srcframe, 0, dstframe, 0)
end

const GraphType = Union{AVFramePtr, GrayTransform, SwsTransform}

graph_input_frame(s::SwsTransform) = s.srcframe
graph_output_frame(s::SwsTransform) = s.dstframe

graph_input_frame(s::AVFramePtr) = s
graph_output_frame(s::AVFramePtr) = s

graph_input_frame(s::GrayTransform) = s.srcframe
graph_output_frame(s::GrayTransform) = s.dstframe

null_graph(s::AVFramePtr) = AVFramePtr(C_NULL)
null_graph(s::SwsTransform) =
    SwsTransform(SwsContextPtr(C_NULL), AVFramePtr(C_NULL), AVFramePtr(C_NULL))
null_graph(s::GrayTransform) = GrayTransform(AVFramePtr(C_NULL), 0,
                                             AVFramePtr(C_NULL), 0)

function exec!(s::SwsTransform)
    src_ptr = unsafe_convert(Ptr{AVFrame}, s.srcframe)
    @preserve s begin
        # s is preserved if r is preserved, hopefully? A little confusing.
        src_data_ptr = field_ptr(src_ptr, :data)
        src_linesize_ptr = field_ptr(src_ptr, :linesize)
        dst_ptr = unsafe_convert(Ptr{AVFrame}, s.dstframe)
        dst_data_ptr = field_ptr(dst_ptr, :data)
        dst_linesize_ptr = field_ptr(dst_ptr, :linesize)

        sigatomic_begin()
        ret = sws_scale(s.sws_context, src_data_ptr, src_linesize_ptr,
                        0, s.dstframe.height, dst_data_ptr,
                        dst_linesize_ptr)
        sigatomic_end()
        ret <= 0 && error("Could not scale frame")
    end
    nothing
end

function scale_gray_frames!(f, dstframe, ::Type{T}, ::Type{S}, srcframe
                            ) where {T<:Unsigned, S<:Unsigned}
    height = srcframe.height
    width = srcframe.width
    ok = height == dstframe.height && width == dstframe.width
    ok || throw(ArgumentError("Frames and not the same size"))
    srcls = srcframe.linesize[1]
    dstls = dstframe.linesize[1]
    srcdp = srcframe.data[1]
    dstdp = dstframe.data[1]
    for h in 0:height-1
        src_line_offset = srcls * h
        dst_line_offset = dstls * h
        src_line_ptr = convert(Ptr{S}, srcdp + src_line_offset)
        dst_line_ptr = convert(Ptr{T}, dstdp + dst_line_offset)
        for w in 1:width
            input_el = unsafe_load(src_line_ptr, w)
            unsafe_store!(dst_line_ptr, f(input_el), w)
        end
    end
end

function gray_range(av_colrange, bitdepth)
    if av_colrange == AVCOL_RANGE_JPEG
        r = (0, 2 ^ bitdepth - 1)
    else
        b = 16 * 2 ^ (bitdepth - 8)
        r = (b, b + 219 * 2 ^ (bitdepth - 8))
    end
    r
end

function make_scale_function(src_t, src_bounds, dst_t, dst_bounds)
    src_span = src_bounds[2] - src_bounds[1]
    dst_span = dst_bounds[2] - dst_bounds[1]
    f = x -> round(dst_t, (dst_span * UInt32(x - src_bounds[1])) //
                   src_span + dst_bounds[1])
end

function make_scale_function(s::GrayTransform)
    src_bounds = gray_range(s.srcframe.color_range, s.src_depth)
    dst_bounds = gray_range(s.dstframe.color_range, s.dst_depth)
    dst_t = strip_interpretation(VIO_PIX_FMT_DEF_ELTYPE_LU[s.dstframe.format])
    src_t = strip_interpretation(VIO_PIX_FMT_DEF_ELTYPE_LU[s.srcframe.format])
    make_scale_function(src_t, src_bounds, dst_t, dst_bounds), src_t, dst_t
end

function exec!(s::GrayTransform)
    f, src_t, dst_t = make_scale_function(s)
    scale_gray_frames!(f, s.dstframe, dst_t, src_t, s.srcframe)
end

exec!(f::AVFramePtr) = nothing

function _vio_sws_getCoefficients(cs_num)
    table_ptr = sws_getCoefficients(cs_num)
    unsafe_wrap(Array, table_ptr, 4)
end

function _vio_primaries_to_sws_table(avcol_pri)
    sws_cs = get(VIO_AVCOL_PRI_TO_SWS_CS, avcol_pri, SWS_CS_DEFAULT)
    _vio_sws_getCoefficients(sws_cs)
end

function sws_get_color_details(sws_context)
    inv_table_dp = Ref(Ptr{Cint}())
    srcRange = Ref{Cint}()
    table_dp = Ref(Ptr{Cint}())
    dstRange = Ref{Cint}()
    brightness = Ref{Cint}()
    contrast = Ref{Cint}()
    saturation = Ref{Cint}()
    ret = sws_getColorspaceDetails(sws_context, inv_table_dp, srcRange,
                                   table_dp, dstRange, brightness, contrast,
                                   saturation)
    ret < 0 && return nothing
    inv_table = Vector{Cint}(undef, 4)
    table = similar(inv_table)
    unsafe_copyto!(pointer(inv_table), inv_table_dp[], 4)
    unsafe_copyto!(pointer(table), table_dp[], 4)
    return inv_table, srcRange[], table, dstRange[], brightness[], contrast[], saturation[]
end

function sws_set_color_details(sws_context, inv_table, src_range, table,
                               dst_range, brightness, contrast, saturation)
    length(inv_table) == length(table) == 4 || throw(ArgumentError("tables wrong size"))
    ret = sws_setColorspaceDetails(sws_context, inv_table, src_range, table,
                                   dst_range, brightness, contrast, saturation)
    ret != -1
end

function sws_update_color_details(sws_context; inv_table = nothing,
                                  src_range = nothing, table = nothing,
                                  dst_range = nothing, brightness = nothing,
                                  contrast = nothing, saturation = nothing)
    outs = sws_get_color_details(sws_context)
    outs === nothing && return false
    def_inv_table, def_src_range, def_table, def_dst_range, def_brightness,
    def_contrast, def_saturation = outs
    new_inv_table  = inv_table  === nothing ? def_inv_table  : inv_table
    new_src_range  = src_range  === nothing ? def_src_range  : src_range
    new_table      = table      === nothing ? def_table      : table
    new_dst_range  = dst_range  === nothing ? def_dst_range  : dst_range
    new_brightness = brightness === nothing ? def_brightness : brightness
    new_contrast   = contrast   === nothing ? def_contrast   : contrast
    new_saturation = saturation === nothing ? def_saturation : saturation
    sws_set_color_details(sws_context, new_inv_table, new_src_range, new_table,
                          new_dst_range, new_brightness, new_contrast,
                          new_saturation)
end

function pix_fmt_to_bits_per_pixel(pix_fmt)
    fmt_desc = av_pix_fmt_desc_get(pix_fmt)
    check_ptr_valid(fmt_desc, false) || error("Unknown pixel format $dst_pix_fmt")
    padded_bits = av_get_padded_bits_per_pixel(fmt_desc)
    bits = av_get_bits_per_pixel(fmt_desc)
    bits, padded_bits
end


text:
using Plots
using StaticArrays
using DifferentialEquations
pyplot()

function simulation(kp::T, kv::T, Ts::Y, t0::Y, tend::Y, r1::Y , r2::Y) where {T<:AbstractMatrix, Y<: AbstractFloat}

    function myrobot(du, u, p, t)
        m = SVector{2}([23.902, 1.285])
        l= SVector{2}([.45, .45])
        r= SVector{2}([0.091, 0.048])
        I= SVector{2}([1.266, 0.093])
        g = 9.81
        θ = SVector{2}(u[1:2])
        dθ = SVector{2}(u[3:4])
        const cos1 = cos(θ[1])
        const cos2 = cos(θ[2])
        const sin1 = sin(θ[1])
        const sin2 = sin(θ[2])
        const cos12 = cos(θ[1]+θ[2])

        #Matriz de inércia
        m11 = m[1] * r[1]^2 + m[2]* (l[1]^2 + r[2]^2 + 2*l[1] * r[2] * cos2)+ I[1] + I[2]
        m12 = m[2] * (r[2]^2 + l[1] *r[2] * cos2) + I[2]
        m22 = m[2] * r[2]^2 + I[2]
        M = SMatrix{2,2}([m11 m12; m12 m22])

        #Matriz C
        h = -m[2] * l[1] * r[2] * sin2
        c11 = h * dθ[2]
        c12 = h * (dθ[1] + dθ[2])
        c21 = -h * dθ[1]
        c22 = 0
        C = SMatrix{2,2}([c11 c12; c21 c22])

        #Matriz gravidade
        g2 = m[2] * r[2] * g * cos12
        g1 = (m[1] * r[1] + m[2] * l[1]) * g * cos1 + g2
        G = SVector{2}([g1; g2])

        xr = SVector{2}([r1, r2])
        e = xr - θ
        tau = kp*e - kv*dθ
        du[1:2] = dθ
        du[3:4] = inv(M)*(tau - C * dθ - G)
    end

    tspan = (t0,tend)
    prob2 = ODEProblem(myrobot,[0., 0., 0., 0.], tspan)
    sol2 = solve(prob2, Tsit5(), saveat = Ts,force_dtmin=true, maxiters = 1e10)

    organize(2,sol2)
end

Ts = 0.08 # Intervalo entre leituras da saída
tend = 20.0 # tempo de simulação
t0 = 0.0 # instante inicial
r1 =  0.6#0.6#1.6# referência junta 1
r2 =  2.2#1.8#2.0# referência junta 2

kp_vec = [280.,55.]
kv_vec = [100.,25.]
kp = SMatrix{2,2}(diagm(kp_vec))
kv = SMatrix{2,2}(diagm(kv_vec))
x, v, t, a, ta, j, tj = simulation(kp, kv, Ts, t0, tend, r1, r2)

p1 = plot(t,x[1], label = "desejado 1")
plot!([r1],seriestype= :hline, label = "referência")
p2 = plot(t,x[2], label = "desejado 2")
plot!([r2],seriestype= :hline, label = "referência")
println("posição final1 = $(x[1][end])| posição final2 = $(x[2][end])|")
plot(p1,p2)


text:
using PAPA
using Test

@testset "PAPA.jl" begin
    # Write your own tests here.
end


text:
formatters = Dict{ Compat.ASCIIString, Function }()

function sprintf1( fmt::Compat.ASCIIString, x )
    global formatters
    f = generate_formatter( fmt )
    f( x )
end

function generate_formatter( fmt::Compat.ASCIIString )
    global formatters
    if haskey( formatters, fmt )
        return formatters[fmt]
    end
    func = @compat Symbol("sprintf_", replace(base64encode(fmt), "=", "!"))

    if !contains( fmt, "'" )
        test = Base.Printf.parse( fmt )
        if length( test ) != 1 || !( typeof( test[1] ) <: Tuple )
            error( "Only one AND undecorated format string is allowed")
        end

        code = quote
            function $func( x )
                @sprintf( $fmt, x )
            end
        end
    else
        conversion = fmt[end]
        if !in( conversion, "sduifF" )
            error( "thousand separator not defined for " * string( conversion ) * " conversion")
        end
        fmtactual = replace( fmt, "'", "", 1 )
        test = Base.Printf.parse( fmtactual )
        if length( test ) != 1 || !( typeof( test[1] ) <: Tuple )
            error( "Only one AND undecorated format string is allowed")
        end
        if in( conversion, "sfF" )
            code = quote
                function $func{T<:Real}( x::T )
                    s = @sprintf( $fmtactual, x )
                    # commas are added to only the numerator
                    if T <: Rational && endswith( $fmtactual, "s" )
                        spos = findfirst( s, '/' )
                        s = addcommas( s[1:spos-1] ) * s[spos:end]
                    else
                        dpos = findfirst( s, '.' )
                        if dpos != 0
                            s = addcommas( s[1:dpos-1] ) * s[ dpos:end ]
                        else # find the rightmost digit
                            for i in length( s ):-1:1
                                if isdigit( s[i] )
                                    s = addcommas( s[1:i] ) * s[i+1:end]
                                    break
                                end
                            end
                        end
                    end
                    s
                end
            end
        else
            code = quote
                function $func( x )
                    s = @sprintf( $fmtactual, x )
                    for i in length( s ):-1:1
                        if isdigit( s[i] )
                            s = addcommas( s[1:i] ) * s[i+1:end]
                            break
                        end
                    end
                    s
                end
            end
        end
    end
    f = eval( code )
    formatters[ fmt ] = f
    f
end

function addcommas( s::Compat.ASCIIString )
    len = length(s)
    t = ""
    for i in 1:3:len
        subs = s[max(1,len-i-1):len-i+1]
        if i == 1
            t = subs
        else
            if match( r"[0-9]", subs ) != nothing
                t = subs * "," * t
            else
                t = subs * t
            end
        end
    end
    return t
end

function generate_format_string(;
        width::Int=-1,
        precision::Int= -1,
        leftjustified::Bool=false,
        zeropadding::Bool=false,
        commas::Bool=false,
        signed::Bool=false,
        positivespace::Bool=false,
        alternative::Bool=false,
        conversion::Compat.ASCIIString="f" #aAdecEfFiosxX
        )
    s = "%"
    if commas
        s *= "'"
    end
    if alternative && in( conversion[1], "aAeEfFoxX" )
        s *= "#"
    end
    if zeropadding && !leftjustified && width != -1
        s *= "0"
    end

    if signed
        s *= "+"
    elseif positivespace
        s *= " "
    end

    if width != -1
        if leftjustified
            s *= "-" * string( width )
        else
            s *= string( width )
        end
    end
    if precision != -1
        s *= "." * string( precision )
    end
    s * conversion
end

function format{T<:Real}( x::T;
        width::Int=-1,
        precision::Int= -1,
        leftjustified::Bool=false,
        zeropadding::Bool=false, # when right-justified, use 0 instead of space to fill
        commas::Bool=false,
        signed::Bool=false, # +/- prefix
        positivespace::Bool=false,
        stripzeros::Bool=(precision== -1),
        parens::Bool=false, # use (1.00) instead of -1.00. Used in finance
        alternative::Bool=false, # usually for hex
        mixedfraction::Bool=false,
        mixedfractionsep::AbstractString="_",
        fractionsep::AbstractString="/", # num / den
        fractionwidth::Int = 0,
        tryden::Int = 0, # if 2 or higher, try to use this denominator, without losing precision
        suffix::AbstractString="", # useful for units/%
        autoscale::Symbol=:none, # :metric, :binary or :finance
        conversion::Compat.ASCIIString=""
        )
    checkwidth = commas
    if conversion == ""
        if T <: AbstractFloat || T <: Rational && precision != -1
            actualconv = "f"
        elseif T <: Unsigned
            actualconv = "x"
        elseif T <: Integer
            actualconv = "d"
        else
            conversion = "s"
            actualconv = "s"
        end
    else
        actualconv = conversion
    end
    if signed && commas
        error( "You cannot use signed (+/-) AND commas at the same time")
    end
    if T <: Rational && conversion == "s"
        stripzeros = false
    end
    if ( T <: AbstractFloat && actualconv == "f" || T <: Integer ) && autoscale != :none
        actualconv = "f"
        if autoscale == :metric
            scales = [
                (1e24, "Y" ),
                (1e21, "Z" ),
                (1e18, "E" ),
                (1e15, "P" ),
                (1e12, "T" ),
                (1e9,  "G"),
                (1e6,  "M"),
                (1e3,  "k") ]
            if abs(x) > 1
                for (mag, sym) in scales
                    if abs(x) >= mag
                        x /= mag
                        suffix = sym * suffix
                        break
                    end
                end
            elseif T <: AbstractFloat
                smallscales = [
                    ( 1e-12, "p" ),
                    ( 1e-9,  "n" ),
                    ( 1e-6,  "μ" ),
                    ( 1e-3,  "m" ) ]
                for (mag,sym) in smallscales
                    if abs(x) < mag*10
                        x /= mag
                        suffix = sym * suffix
                        break
                    end
                end
            end
        else
            if autoscale == :binary
                scales = [
                    (1024.0 ^8,  "Yi" ),
                    (1024.0 ^7,  "Zi" ),
                    (1024.0 ^6,  "Ei" ),
                    (1024^5,  "Pi" ),
                    (1024^4,  "Ti" ),
                    (1024^3,  "Gi"),
                    (1024^2,  "Mi"),
                    (1024,    "Ki")
                ]
            else # :finance
                scales = [
                    (1e12, "t" ),
                    (1e9,  "b"),
                    (1e6,  "m"),
                    (1e3,  "k") ]
            end
            for (mag, sym) in scales
                if abs(x) >= mag
                    x /= mag
                    suffix = sym * suffix
                    break
                end
            end
        end
    end

    nonneg = x >= 0
    fractional = 0
    if T <: Rational && mixedfraction
        actualconv = "d"
        actualx = trunc( Int, x )
        fractional = abs(x) - abs(actualx)
    else
        if parens && !in( actualconv[1], "xX" )
            actualx = abs(x)
        else
            actualx = x
        end
    end
    s = sprintf1( generate_format_string( width=width,
        precision=precision,
        leftjustified=leftjustified,
        zeropadding=zeropadding,
        commas=commas,
        signed=signed,
        positivespace=positivespace,
        alternative=alternative,
        conversion=actualconv
    ),actualx)

    if T <:Rational && conversion == "s"
        if mixedfraction && fractional != 0
            num = fractional.num
            den = fractional.den
            if tryden >= 2 && mod( tryden, den ) == 0
                num *= div(tryden,den)
                den = tryden
            end
            fs = string( num ) * fractionsep * string(den)
            if length(fs) < fractionwidth
                fs = repeat( "0", fractionwidth - length(fs) ) * fs
            end
            s = rstrip(s)
            if actualx != 0
                s = rstrip(s) * mixedfractionsep * fs
            else
                if !nonneg
                    s = "-" * fs
                else
                    s = fs
                end
            end
            checkwidth = true
        elseif !mixedfraction
            s = replace( s, "//", fractionsep )
            checkwidth = true
        end
    elseif stripzeros && in( actualconv[1], "fFeEs" )
        dpos = findfirst( s, '.')
        if in( actualconv[1], "eEs" )
            if in( actualconv[1], "es" )
                epos = findfirst( s, 'e' )
            else
                epos = findfirst( s, 'E' )
            end
            if epos == 0
                rpos = length( s )
            else
                rpos = epos-1
            end
        else
            rpos = length(s)
        end
        # rpos at this point is the rightmost possible char to start
        # stripping
        stripfrom = rpos+1
        for i = rpos:-1:dpos+1
            if s[i] == '0'
                stripfrom = i
            elseif s[i] ==' '
                continue
            else
                break
            end
        end
        if stripfrom <= rpos
            if stripfrom == dpos+1 # everything after decimal is 0, so strip the decimal too
                stripfrom = dpos
            end
            s = s[1:stripfrom-1] * s[rpos+1:end]
            checkwidth = true
        end
    end

    s *= suffix

    if parens && !in( actualconv[1], "xX" )
        # if zero or positive, we still need 1 white space on the right
        if nonneg
            s = " " * strip(s) * " "
        else
            s = "(" * strip(s) * ")"
        end

        checkwidth = true
    end

    if checkwidth && width != -1
        if length(s) > width
            s = replace( s, " ", "", length(s)-width )
            if length(s) > width && endswith( s, " " )
                s = reverse( replace( reverse(s), " ", "", length(s)-width ) )
            end
            if length(s) > width
                s = replace( s, ",", "", length(s)-width )
            end
        elseif length(s) < width
            if leftjustified
                s = s * repeat( " ", width - length(s) )
            else
                s = repeat( " ", width - length(s) ) * s
            end
        end
    end

    s
end


text:
"""
    Population{T<:Soma,
               NT<:AbstractArray{T, 1},
               WT<:AbstractMatrix{<:Real},
               ST<:AbstractArray{<:AbstractSynapse, 2},
               LT<:AbstractLearner} <: AbstractArray{T, 1}

A population of neurons is an array of [`Soma`](@ref)s,
  a weighted matrix of [synapses](@ref Synapse Models), and a [learner](@ref Learning).

Fields:
- `somas::AbstractArray{<:Soma, 1}`: a vector of somas
- `weights::AbstractMatrix{<:Real}`: a weight matrix
- `synapses::AbstractArray{<:AbstractSynapse, 2}`: a matrix of synapses
- `learner::AbstractLearner`: a learning mechanism
"""
struct Population{T<:Soma, NT<:AbstractArray{T, 1}, WT<:AbstractMatrix{<:Real}, ST<:AbstractArray{<:AbstractSynapse, 2}, LT<:AbstractLearner} <: AbstractArray{T, 1}
    somas::NT
    weights::WT
    synapses::ST
    learner::LT
end

"""
    size(pop::Population)

Return the number of neurons in a population.
"""
Base.size(pop::Population) = length(pop.somas)

Base.IndexStyle(::Type{<:Population}) = IndexLinear()
Base.getindex(pop::Population, i::Int) = Neuron(pop.synapses[:, i], pop.somas[i])
function Base.setindex!(pop::Population, neuron::Neuron, i::Int)
    pop.synapses[:, i] = neuron.synapses
    pop.somas[i] = neuron.soma
end

Base.show(io::IO, pop::Population{T, <:Any, <:Any, ST, LT}) where {T, ST, LT} =
    print(io, "Population{$(nameof(eltype(pop.somas.body))), $(nameof(eltype(ST))), $(nameof(LT))}($(size(pop)))")
Base.show(io::IO, ::MIME"text/plain", pop::Population) = show(io, pop)

_instantiate(x::AbstractArray, I...) = x[I...]
_instantiate(x, I...) = x()

"""
    Population(weights::AbstractMatrix{<:Real};
               cell = LIF, synapse = Synapse.Delta,
               threshold = Threshold.Ideal, learner = George())
Create a population by specifying the `weights`
  and optionally the cell type, synapse type, threshold type, and learner.

# Keyword Fields:
- `cell::AbstractCell`: a constructor or function that creates a cell body, or a vector of pre-constructed cells
- `synapse::AbstractSynapse`: a constructor or function that creates a synapse, or a matrix of pre-constructed synapses
- `threshold::AbstractThreshold`: a constructor or function that creates a threshold, or a vector of pre-constructed thresholds
- `learner::AbstractLearner`: a learner object
"""
function Population(weights::AbstractMatrix{<:Real}; cell = LIF, synapse = Synapse.Delta, threshold = Threshold.Ideal, learner = George())
    n = _checkweights(weights)
    synapses = StructArray(_instantiate(synapse, i, j) for i in 1:n, j in 1:n)
    somas = StructArray(Soma(_instantiate(cell, i), _instantiate(threshold, i)) for i in 1:n;
                        unwrap = t -> t <: AbstractCell || t <: AbstractThreshold)

    Population(somas, weights, synapses, learner)
end

"""
    neurons(pop::Population)

Return an array of neurons within the population.
"""
neurons(pop::Population) = pop[1:end]

"""
    synapses(pop::Population)

Return an array of edges representing the synapses within the population.
"""
synapses(pop::Population) = pop.synapses

function _processspikes!(pop::Population, spikes; dt::Real = 1.0)
    # record spikes with learner
    record!(pop.learner, pop.weights, spikes; dt = dt)

    # excite post-synaptic neurons
    map((row, s) -> (s > 0) && excite!(row, s + 1), eachrow(pop.synapses), spikes)

    # apply refactory period to synapses
    map((col, s) -> (s > 0) && spike!(col, s), eachcol(pop.synapses), spikes)
end

"""
    evaluate!(pop::Population, t::Integer; dt::Real = 1.0, dense = false, inputs = nothing)
    (::Population)(t::Integer; dt::Real = 1.0, dense = false)

Evaluate a population of neurons at time step `t`.
Return a vector of time stamps (`t` if the neuron spiked and zero otherwise).
"""
function evaluate!(pop::Population, t::Integer; dt::Real = 1.0, dense = false, inputs = nothing)
    # evalute inputs
    !isnothing(inputs) && excite!(pop.somas, [input(t; dt = dt) for input in inputs])

    # evaluate synapses and excite neuron bodies w/ current
    current = vec(reduce(+, pop.weights .* evaluate!(pop.synapses, t; dt = dt); dims = 1))
    excite!(pop.somas, current)

    # evaluate somas
    spikes = evaluate!(pop.somas, t; dt = dt)

    # process spike events
    _processspikes!(pop, spikes; dt = dt)

    return spikes
end
(pop::Population)(t::Integer; kwargs...) = evaluate!(pop, t; kwargs...)

"""
    update!(pop::Population, t::Integer; dt::Real = 1.0)

Update synaptic weights within population according to `pop.learner`.
"""
update!(pop::Population, t::Integer; dt::Real = 1.0) = update!(pop.learner, pop.weights, t; dt = dt)

"""
    reset!(pop::Population)

Reset `pop.synapses` and `pop.somas`.
"""
function reset!(pop::Population)
    reset!(pop.synapses)
    reset!(pop.somas)
end

function _recordspikes!(dict::Dict{Int, Array{Int, 1}}, spikes)
    for (id, spiketime) in enumerate(spikes)
        if spiketime > 0
            record = get!(dict, id, Int[])
            push!(record, spiketime)
        end
    end
end

"""
    simulate!(pop::Population, dt::Real = 1.0)

Simulate a population of neurons. Optionally specify a learner. The `prespike` and
`postspike` functions will be called immediately after either event occurs.

Fields:
- `pop::Population`: the population to simulate
- `T::Integer`: number of time steps to simulate
- `dt::Real`: the simulation time step
- `cb::Function`: a callback function that is called after event evaluation (expects `(neuron_id, t)` as input)
- `dense::Bool`: set to `true` to evaluate every time step even in the absence of events
"""
function simulate!(pop::Population, T::Integer; dt::Real = 1.0, cb = () -> (), dense = false, inputs = nothing)
    spiketimes = Dict([(i, Int[]) for i in 1:size(pop)])

    for t = 1:T
        # evaluate callback
        cb()

        # evaluate population once
        spikes = evaluate!(pop, t; dt = dt, dense = dense, inputs = inputs)

        # record spike time
        _recordspikes!(spiketimes, spikes)

        # update weights
        update!(pop, t; dt = dt)
    end

    return spiketimes
end


text:
include("censored.jl")
include("cross_validate.jl")
include("precision_at_k.jl")
include("simple_glrms.jl")
#include("fit_rdataset.jl")


text:
# ============================================================
# Flux for Solution Structures
# ============================================================

"""
Evolution of particle transport

* kinetic flux vector splitting (KFVS)
* kinetic central-upwind (KCU)
* unified gas-kinetic scheme (UGKS)

"""
function evolve!(
    KS::SolverSet,
    uq::AbstractUQ,
    sol::Solution1D1F,
    flux::Flux1D1F,
    dt::AbstractFloat;
    mode = :kfvs::Symbol,
)

    if mode == :kfvs
        @inbounds Threads.@threads for i in eachindex(flux.fw)
            for j in axes(sol.w[1], 2) # over gPC coefficients or quadrature points
                fw = @view flux.fw[i][:, j]
                ff = @view flux.ff[i][:, j]

                KitBase.flux_kfvs!(
                    fw,
                    ff,
                    sol.f[i-1][:, j] .+ 0.5 .* KS.pSpace.dx[i-1] .* sol.sf[i-1][:, j],
                    sol.f[i][:, j] .- 0.5 .* KS.pSpace.dx[i] .* sol.sf[i][:, j],
                    KS.vSpace.u,
                    KS.vSpace.weights,
                    dt,
                    sol.sf[i-1][:, j],
                    sol.sf[i][:, j],
                )
            end
        end
    else
        throw("flux mode isn't available")
    end

end

function evolve!(
    KS::SolverSet,
    uq::AbstractUQ,
    sol::Solution1D2F,
    flux::Flux1D2F,
    dt::AbstractFloat;
    mode = :kfvs::Symbol,
)

    if mode == :kfvs
        @inbounds Threads.@threads for i in eachindex(flux.fw)
            for j in axes(sol.w[1], 2) # over gPC coefficients or quadrature points
                fw = @view flux.fw[i][:, j]
                fh = @view flux.fh[i][:, j]
                fb = @view flux.fb[i][:, j]

                KitBase.flux_kfvs!(
                    fw,
                    fh,
                    fb,
                    sol.h[i-1][:, j] .+ 0.5 .* KS.pSpace.dx[i-1] .* sol.sh[i-1][:, j],
                    sol.b[i-1][:, j] .+ 0.5 .* KS.pSpace.dx[i-1] .* sol.sb[i-1][:, j],
                    sol.h[i][:, j] .- 0.5 .* KS.pSpace.dx[i] .* sol.sh[i][:, j],
                    sol.b[i][:, j] .- 0.5 .* KS.pSpace.dx[i] .* sol.sb[i][:, j],
                    KS.vSpace.u,
                    KS.vSpace.weights,
                    dt,
                    sol.sh[i-1][:, j],
                    sol.sb[i-1][:, j],
                    sol.sh[i][:, j],
                    sol.sb[i][:, j],
                )
            end
        end
    else
        throw("flux mode isn't available")
    end

end

#--- 2D case ---#
function evolve!(
    KS::SolverSet,
    uq::AbstractUQ,
    sol::Solution2D2F,
    flux::Flux2D2F,
    dt::AbstractFloat;
    mode = :kfvs::Symbol,
)

    if mode == :kfvs

        @inbounds for j = 1:KS.pSpace.ny
            for i = 1:KS.pSpace.nx+1
                un = KS.vSpace.u .* flux.n1[i, j][1] .+ KS.vSpace.v .* flux.n1[i, j][2]
                ut = KS.vSpace.v .* flux.n1[i, j][1] .- KS.vSpace.u .* flux.n1[i, j][2]

                for k in axes(sol.w[1, 1], 2)
                    fw1 = @view flux.fw1[i, j][:, k]
                    fh1 = @view flux.fh1[i, j][:, :, k]
                    fb1 = @view flux.fb1[i, j][:, :, k]

                    flux_kfvs!(
                        fw1,
                        fh1,
                        fb1,
                        sol.h[i-1, j][:, :, k] .+
                        0.5 .* KS.pSpace.dx[i-1, j] .* sol.sh[i-1, j][:, :, k, 1],
                        sol.b[i-1, j][:, :, k] .+
                        0.5 .* KS.pSpace.dx[i-1, j] .* sol.sb[i-1, j][:, :, k, 1],
                        sol.h[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dx[i, j] .* sol.sh[i, j][:, :, k, 1],
                        sol.b[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dx[i, j] .* sol.sb[i, j][:, :, k, 1],
                        un,
                        ut,
                        KS.vSpace.weights,
                        dt,
                        0.5 * (KS.pSpace.dy[i-1, j] + KS.pSpace.dy[i, j]),
                        sol.sh[i-1, j][:, :, k, 1],
                        sol.sb[i-1, j][:, :, k, 1],
                        sol.sh[i, j][:, :, k, 1],
                        sol.sb[i, j][:, :, k, 1],
                    )
                    flux.fw1[i, j][:, k] .=
                        global_frame(fw1, flux.n1[i, j][1], flux.n1[i, j][2])
                end
            end
        end

        @inbounds for j = 1:KS.pSpace.ny+1
            for i = 1:KS.pSpace.nx
                vn = KS.vSpace.u .* flux.n2[i, j][1] .+ KS.vSpace.v .* flux.n2[i, j][2]
                vt = KS.vSpace.v .* flux.n2[i, j][1] .- KS.vSpace.u .* flux.n2[i, j][2]

                for k in axes(sol.w[1, 1], 2)
                    fw2 = @view flux.fw2[i, j][:, k]
                    fh2 = @view flux.fh2[i, j][:, :, k]
                    fb2 = @view flux.fb2[i, j][:, :, k]

                    flux_kfvs!(
                        fw2,
                        fh2,
                        fb2,
                        sol.h[i, j-1][:, :, k] .+
                        0.5 .* KS.pSpace.dy[i, j-1] .* sol.sh[i, j-1][:, :, k, 2],
                        sol.b[i, j-1][:, :, k] .+
                        0.5 .* KS.pSpace.dy[i, j-1] .* sol.sb[i, j-1][:, :, k, 2],
                        sol.h[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dy[i, j] .* sol.sh[i, j][:, :, k, 2],
                        sol.b[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dy[i, j] .* sol.sb[i, j][:, :, k, 2],
                        vn,
                        vt,
                        KS.vSpace.weights,
                        dt,
                        0.5 * (KS.pSpace.dx[i, j-1] + KS.pSpace.dx[i, j]),
                        sol.sh[i, j-1][:, :, k, 2],
                        sol.sb[i, j-1][:, :, k, 2],
                        sol.sh[i, j][:, :, k, 2],
                        sol.sb[i, j][:, :, k, 2],
                    )
                    flux.fw2[i, j][:, k] .=
                        global_frame(fw2, flux.n2[i, j][1], flux.n2[i, j][2])
                end
            end
        end

    elseif mode == :kcu

        @inbounds for j = 1:KS.pSpace.ny
            for i = 1:KS.pSpace.nx+1
                un = KS.vSpace.u .* flux.n1[i, j][1] .+ KS.vSpace.v .* flux.n1[i, j][2]
                ut = KS.vSpace.v .* flux.n1[i, j][1] .- KS.vSpace.u .* flux.n1[i, j][2]

                for k in axes(sol.w[1, 1], 2)
                    fw1 = @view flux.fw1[i, j][:, k]
                    fh1 = @view flux.fh1[i, j][:, :, k]
                    fb1 = @view flux.fb1[i, j][:, :, k]

                    flux_kcu!(
                        fw1,
                        fh1,
                        fb1,
                        local_frame(
                            sol.w[i-1, j][:, k] .+
                            0.5 .* KS.pSpace.dx[i-1, j] .* sol.sw[i-1, j][:, k, 1],
                            flux.n1[i, j][1],
                            flux.n1[i, j][2],
                        ),
                        sol.h[i-1, j][:, :, k] .+
                        0.5 .* KS.pSpace.dx[i-1, j] .* sol.sh[i-1, j][:, :, k, 1],
                        sol.b[i-1, j][:, :, k] .+
                        0.5 .* KS.pSpace.dx[i-1, j] .* sol.sb[i-1, j][:, :, k, 1],
                        local_frame(
                            sol.w[i, j][:, k] .-
                            0.5 .* KS.pSpace.dx[i, j] .* sol.sw[i, j][:, k, 1],
                            flux.n1[i, j][1],
                            flux.n1[i, j][2],
                        ),
                        sol.h[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dx[i, j] .* sol.sh[i, j][:, :, k, 1],
                        sol.b[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dx[i, j] .* sol.sb[i, j][:, :, k, 1],
                        un,
                        ut,
                        KS.vSpace.weights,
                        KS.gas.K,
                        KS.gas.γ,
                        KS.gas.μᵣ,
                        KS.gas.ω,
                        KS.gas.Pr,
                        dt,
                        0.5 * (KS.pSpace.dy[i-1, j] + KS.pSpace.dy[i, j]),
                    )
                    flux.fw1[i, j][:, k] .=
                        global_frame(fw1, flux.n1[i, j][1], flux.n1[i, j][2])
                end
            end
        end

        @inbounds for j = 1:KS.pSpace.ny+1
            for i = 1:KS.pSpace.nx
                vn = KS.vSpace.u .* flux.n2[i, j][1] .+ KS.vSpace.v .* flux.n2[i, j][2]
                vt = KS.vSpace.v .* flux.n2[i, j][1] .- KS.vSpace.u .* flux.n2[i, j][2]

                for k in axes(sol.w[1, 1], 2)
                    fw2 = @view flux.fw2[i, j][:, k]
                    fh2 = @view flux.fh2[i, j][:, :, k]
                    fb2 = @view flux.fb2[i, j][:, :, k]

                    flux_kcu!(
                        fw2,
                        fh2,
                        fb2,
                        local_frame(
                            sol.w[i, j-1][:, k] .+
                            0.5 .* KS.pSpace.dy[i, j-1] .* sol.sw[i, j-1][:, k, 2],
                            flux.n2[i, j][1],
                            flux.n2[i, j][2],
                        ),
                        sol.h[i, j-1][:, :, k] .+
                        0.5 .* KS.pSpace.dy[i, j-1] .* sol.sh[i, j-1][:, :, k, 2],
                        sol.b[i, j-1][:, :, k] .+
                        0.5 .* KS.pSpace.dy[i, j-1] .* sol.sb[i, j-1][:, :, k, 2],
                        local_frame(
                            sol.w[i, j][:, k] .-
                            0.5 .* KS.pSpace.dy[i, j] .* sol.sw[i, j][:, k, 2],
                            flux.n2[i, j][1],
                            flux.n2[i, j][2],
                        ),
                        sol.h[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dy[i, j] .* sol.sh[i, j][:, :, k, 2],
                        sol.b[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dy[i, j] .* sol.sb[i, j][:, :, k, 2],
                        vn,
                        vt,
                        KS.vSpace.weights,
                        KS.gas.K,
                        KS.gas.γ,
                        KS.gas.μᵣ,
                        KS.gas.ω,
                        KS.gas.Pr,
                        dt,
                        0.5 * (KS.pSpace.dx[i, j-1] + KS.pSpace.dx[i, j]),
                    )
                    flux.fw2[i, j][:, k] .=
                        global_frame(fw2, flux.n2[i, j][1], flux.n2[i, j][2])
                end
            end
        end

    elseif mode == :ugks

        @inbounds for j = 1:KS.pSpace.ny
            for i = 1:KS.pSpace.nx+1
                un = KS.vSpace.u .* flux.n1[i, j][1] .+ KS.vSpace.v .* flux.n1[i, j][2]
                ut = KS.vSpace.v .* flux.n1[i, j][1] .- KS.vSpace.u .* flux.n1[i, j][2]

                for k in axes(sol.w[1, 1], 2)
                    fw1 = @view flux.fw1[i, j][:, k]
                    fh1 = @view flux.fh1[i, j][:, :, k]
                    fb1 = @view flux.fb1[i, j][:, :, k]

                    flux_ugks!(
                        fw1,
                        fh1,
                        fb1,
                        local_frame(
                            sol.w[i-1, j][:, k] .+
                            0.5 .* KS.pSpace.dx[i-1, j] .* sol.sw[i-1, j][:, k, 1],
                            flux.n1[i, j][1],
                            flux.n1[i, j][2],
                        ),
                        sol.h[i-1, j][:, :, k] .+
                        0.5 .* KS.pSpace.dx[i-1, j] .* sol.sh[i-1, j][:, :, k, 1],
                        sol.b[i-1, j][:, :, k] .+
                        0.5 .* KS.pSpace.dx[i-1, j] .* sol.sb[i-1, j][:, :, k, 1],
                        local_frame(
                            sol.w[i, j][:, k] .-
                            0.5 .* KS.pSpace.dx[i, j] .* sol.sw[i, j][:, k, 1],
                            flux.n1[i, j][1],
                            flux.n1[i, j][2],
                        ),
                        sol.h[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dx[i, j] .* sol.sh[i, j][:, :, k, 1],
                        sol.b[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dx[i, j] .* sol.sb[i, j][:, :, k, 1],
                        un,
                        ut,
                        KS.vSpace.weights,
                        KS.gas.K,
                        KS.gas.γ,
                        KS.gas.μᵣ,
                        KS.gas.ω,
                        KS.gas.Pr,
                        dt,
                        0.5 * KS.pSpace.dx[i-1, j],
                        0.5 * KS.pSpace.dx[i, j],
                        0.5 * (KS.pSpace.dy[i-1, j] + KS.pSpace.dy[i, j]),
                        sol.sh[i-1, j][:, :, k, 1],
                        sol.sb[i-1, j][:, :, k, 1],
                        sol.sh[i, j][:, :, k, 1],
                        sol.sb[i, j][:, :, k, 1],
                    )
                    flux.fw1[i, j][:, k] .=
                        global_frame(fw1, flux.n1[i, j][1], flux.n1[i, j][2])
                end
            end
        end

        @inbounds for j = 1:KS.pSpace.ny+1
            for i = 1:KS.pSpace.nx
                vn = KS.vSpace.u .* flux.n2[i, j][1] .+ KS.vSpace.v .* flux.n2[i, j][2]
                vt = KS.vSpace.v .* flux.n2[i, j][1] .- KS.vSpace.u .* flux.n2[i, j][2]

                for k in axes(sol.w[1, 1], 2)
                    fw2 = @view flux.fw2[i, j][:, k]
                    fh2 = @view flux.fh2[i, j][:, :, k]
                    fb2 = @view flux.fb2[i, j][:, :, k]

                    flux_ugks!(
                        fw2,
                        fh2,
                        fb2,
                        local_frame(
                            sol.w[i, j-1][:, k] .+
                            0.5 .* KS.pSpace.dy[i, j-1] .* sol.sw[i, j-1][:, k, 2],
                            flux.n2[i, j][1],
                            flux.n2[i, j][2],
                        ),
                        sol.h[i, j-1][:, :, k] .+
                        0.5 .* KS.pSpace.dy[i, j-1] .* sol.sh[i, j-1][:, :, k, 2],
                        sol.b[i, j-1][:, :, k] .+
                        0.5 .* KS.pSpace.dy[i, j-1] .* sol.sb[i, j-1][:, :, k, 2],
                        local_frame(
                            sol.w[i, j][:, k] .-
                            0.5 .* KS.pSpace.dy[i, j] .* sol.sw[i, j][:, k, 2],
                            flux.n2[i, j][1],
                            flux.n2[i, j][2],
                        ),
                        sol.h[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dy[i, j] .* sol.sh[i, j][:, :, k, 2],
                        sol.b[i, j][:, :, k] .-
                        0.5 .* KS.pSpace.dy[i, j] .* sol.sb[i, j][:, :, k, 2],
                        vn,
                        vt,
                        KS.vSpace.weights,
                        KS.gas.K,
                        KS.gas.γ,
                        KS.gas.μᵣ,
                        KS.gas.ω,
                        KS.gas.Pr,
                        dt,
                        0.5 * KS.pSpace.dy[i, j-1],
                        0.5 * KS.pSpace.dy[i, j],
                        0.5 * (KS.pSpace.dx[i, j-1] + KS.pSpace.dx[i, j]),
                        sol.sh[i, j-1][:, :, k, 2],
                        sol.sb[i, j-1][:, :, k, 2],
                        sol.sh[i, j][:, :, k, 2],
                        sol.sb[i, j][:, :, k, 2],
                    )
                    flux.fw2[i, j][:, k] .=
                        global_frame(fw2, flux.n2[i, j][1], flux.n2[i, j][2])
                end
            end
        end

    else

        throw("flux mode isn't available")

    end # if

end

"""
Maxwell's diffusive boundary flux

"""
function evolve_boundary!(
    bc::Array,
    KS::SolverSet,
    sol::Solution2D2F,
    flux::Flux2D2F,
    dt::AbstractFloat;
    mode = :maxwell::Symbol,
)

    if mode == :maxwell

        @inbounds for j = 1:KS.pSpace.ny
            un = KS.vSpace.u .* flux.n1[1, j][1] .+ KS.vSpace.v .* flux.n1[1, j][2]
            ut = KS.vSpace.v .* flux.n1[1, j][1] .- KS.vSpace.u .* flux.n1[1, j][2]

            for k in axes(sol.w[1, 1], 2)
                bcL = local_frame(bc[1][:, k], flux.n1[1, j][1], flux.n1[1, j][2])
                bcR = local_frame(bc[2][:, k], flux.n1[1, j][1], flux.n1[1, j][2])

                fw1 = @view flux.fw1[1, j][:, k]
                fh1 = @view flux.fh1[1, j][:, :, k]
                fb1 = @view flux.fb1[1, j][:, :, k]
                flux_boundary_maxwell!(
                    fw1,
                    fh1,
                    fb1,
                    bcL, # left
                    sol.h[1, j][:, :, k],
                    sol.b[1, j][:, :, k],
                    un,
                    ut,
                    KS.vSpace.weights,
                    KS.gas.K,
                    dt,
                    KS.pSpace.dy[1, j],
                    1,
                )
                flux.fw1[1, j][:, k] .=
                    global_frame(flux.fw1[1, j][:, k], flux.n1[1, j][1], flux.n1[1, j][2])

                fw2 = @view flux.fw1[end, j][:, k]
                fh2 = @view flux.fh1[end, j][:, :, k]
                fb2 = @view flux.fb1[end, j][:, :, k]
                flux_boundary_maxwell!(
                    fw2,
                    fh2,
                    fb2,
                    bcR, # right
                    sol.h[KS.pSpace.nx, j][:, :, k],
                    sol.b[KS.pSpace.nx, j][:, :, k],
                    un,
                    ut,
                    KS.vSpace.weights,
                    KS.gas.K,
                    dt,
                    KS.pSpace.dy[KS.pSpace.nx, j],
                    -1,
                )
                flux.fw1[end, j][:, k] .= global_frame(
                    flux.fw1[end, j][:, k],
                    flux.n1[end, j][1],
                    flux.n1[end, j][2],
                )

            end
        end

        @inbounds for i = 1:KS.pSpace.nx
            vn = KS.vSpace.u .* flux.n2[i, 1][1] .+ KS.vSpace.v .* flux.n2[i, 1][2]
            vt = KS.vSpace.v .* flux.n2[i, 1][1] .- KS.vSpace.u .* flux.n2[i, 1][2]

            for k in axes(sol.w[1, 1], 2)
                bcU = local_frame(bc[3][:, k], flux.n2[i, 1][1], flux.n2[i, 1][2])
                bcD = local_frame(bc[4][:, k], flux.n2[i, 1][1], flux.n2[i, 1][2])

                fw3 = @view flux.fw2[i, 1][:, k]
                fh3 = @view flux.fh2[i, 1][:, :, k]
                fb3 = @view flux.fb2[i, 1][:, :, k]
                flux_boundary_maxwell!(
                    fw3,
                    fh3,
                    fb3,
                    bcD, # down
                    sol.h[i, 1][:, :, k],
                    sol.b[i, 1][:, :, k],
                    vn,
                    vt,
                    KS.vSpace.weights,
                    KS.gas.K,
                    dt,
                    KS.pSpace.dx[i, 1],
                    1,
                )
                flux.fw2[i, 1][:, k] .=
                    global_frame(flux.fw2[i, 1][:, k], flux.n2[i, 1][1], flux.n2[i, 1][2])

                fw4 = @view flux.fw2[i, end][:, k]
                fh4 = @view flux.fh2[i, end][:, :, k]
                fb4 = @view flux.fb2[i, end][:, :, k]
                flux_boundary_maxwell!(
                    fw4,
                    fh4,
                    fb4,
                    bcU, # down
                    sol.h[i, KS.pSpace.ny][:, :, k],
                    sol.b[i, KS.pSpace.ny][:, :, k],
                    vn,
                    vt,
                    KS.vSpace.weights,
                    KS.gas.K,
                    dt,
                    KS.pSpace.dx[i, KS.pSpace.ny],
                    -1,
                )
                flux.fw2[i, end][:, k] .= global_frame(
                    flux.fw2[i, end][:, k],
                    flux.n2[i, end][1],
                    flux.n2[i, end][2],
                )
            end
        end

    end # if

end


text:
function ELBO(model::GP, X, y, pr_mean, kernel)
    setpr_mean!(model.f, pr_mean)
    setkernel!(model.f, kernel)
    state = compute_kernel_matrices(model, (;), X, true)
    return objective(model, state, y)
end

function ELBO(model::AbstractGPModel, X, y, pr_means, kernels, state)
    setpr_means!(model, pr_means)
    setkernels!(model, kernels)
    state = compute_kernel_matrices(model, state, X, true)
    return objective(model, state, y)
end

function ELBO(model::AbstractGPModel, X, y, pr_means, kernels, Zs, state)
    setpr_means!(model, pr_means)
    setkernels!(model, kernels)
    setZs!(model, Zs)
    state = compute_kernel_matrices(model, state, X, true)
    return objective(model, state, y)
end

# External ELBO call on internal and new data
@traitfn function ELBO(model::TGP) where {TGP <: AbstractGPModel; IsFull{TGP}}
    return ELBO(model, input(model.data), output(model.data))
end

function ELBO(model::AbstractGPModel, X::AbstractMatrix, y::AbstractArray; obsdim=1)
    return ELBO(model, KernelFunctions.vec_of_vecs(X; obsdim), y)
end

function ELBO(model::AbstractGPModel, X::AbstractVector, y::AbstractArray)
    y = treat_labels!(y, likelihood(model))
    state = compute_kernel_matrices(model, (;), X, true)
    if inference(model) isa AnalyticVI
        local_vars = init_local_vars(likelihood(model), length(X))
        local_vars = local_updates!(
            local_vars,
            likelihood(model),
            y,
            mean_f(model, state.kernel_matrices),
            var_f(model, state.kernel_matrices),
        )
        state = merge(state, (; local_vars))
    end
    return ELBO(model, state, y)
end

function ELBO(
    model::OnlineSVGP, state::NamedTuple, X::AbstractMatrix, y::AbstractArray; obsdim=1
)
    return ELBO(model, state, KernelFunctions.vec_of_vecs(X; obsdim), y)
end

function ELBO(model::OnlineSVGP, state::NamedTuple, X::AbstractVector, y::AbstractArray)
    y = treat_labels!(y, likelihood(model))
    state = compute_kernel_matrices(model, state, X, true)
    if inference(model) isa AnalyticVI
        local_vars = init_local_vars(likelihood(model), length(X))
        local_vars = local_updates!(
            local_vars,
            likelihood(model),
            y,
            mean_f(model, state.kernel_matrices),
            var_f(model, state.kernel_matrices),
        )
        state = merge(state, (; local_vars))
    end
    return ELBO(model, state, y)
end


text:
export bibliography,
       file
using BibTeX
using BibTeXFormat
using Mustache
read_file(path) = read(open(path), String)
function bibliography(cites::Vector{String})
    bibliography = Bibliography(read_file(joinpath(dirname(@__FILE__), "SpectralClustering.bib")))
    formatted_entries = format_entries(UNSRTAlphaStyle,bibliography)
    return write_to_string( HTMLBackend() ,formatted_entries,  cites)
end
function file(name::String)
    file = read_file(joinpath(dirname(@__FILE__),"src", "js",name))
    return file
end

function file(name::String, data::Dict{String, String})
    file = read_file(joinpath(dirname(@__FILE__),"src", "js",name))
    return render(file, data)
end


text:
using TypeDBClient
g = TypeDBClient

client = g.CoreClient("127.0.0.1",1729)
sess = g.CoreSession(client, context[:db_name], g.Proto.Session_Type.SCHEMA, request_timeout = Inf)
trans = g.transaction(sess, g.Proto.Transaction_Type.WRITE)

g.put(g.ConceptManager(trans), AttributeType, "is-alive", VALUE_TYPE.BOOLEAN )
g.put(g.ConceptManager(trans), EntityType, "person")

function cm_call(f::Function, session::Session)
    transaction = Transaction(session, session.sessionID, g.Proto.Transaction_Type.READ, TypeDBOptions())
    concept_manager = ConceptManager(transaction)
    res = f(concept_manager)
    close(transaction)
    return res
end

function thing_call(f::Function, session::Session, name::AbstractString)
    cm_call(session) do concept_manager
        thing_type_res = TypeDBClient.execute(
            concept_manager, TypeDBClient.ConceptManagerRequestBuilder.get_thing_type_req(name)
        )
        thing_type = TypeDBClient.instantiate(thing_type_res.get_thing_type_res.thing_type)
        remote_thing = TypeDBClient.as_remote(thing_type, concept_manager.transaction)
        f(remote_thing)
    end
end

res = thing_call(thing -> get_owns(thing), sess, "person")

g.commit(trans)


text:
# This file was generated, do not modify it. # hide
pipe = @pipeline MyPipe(X -> coerce(X, :age=>Continuous),
                       hot = OneHotEncoder(),
                       knn = KNNRegressor(K=3),
                       target = UnivariateStandardizer());


text:
# =========================================================================
#    BASIC DATA GENERATION
# =========================================================================

const vertex_data = Ref{Any}(nothing)
V() = sort!(collect(keys(vertex_data[])))
V(i) = vertex_data[][i]

const vehicle_data = Ref{Any}(nothing)
K() = sort!(collect(keys(vehicle_data[])))
K(k) = vehicle_data[][k]

const period_data = Ref{Any}(nothing)
T() = period_data[]
T(t) = period_data[][t]

const demand_data = Ref{Any}(nothing)
d() = demand_data[]
d(i,t) = demand_data[][i,t]

const distance_data = Ref{Any}(nothing)
dist() = distance_data[]
dist(i,j) = distance_data[][i,j]

function extract!(path::String;f::String) #extract from excel
    xf = XLSX.readxlsx(path) #READ WORKSHEET
    data = Dict{Symbol,DataFrame}() #DATAFRAME DICT

    for sheets in XLSX.sheetnames(xf) #TURN SHEETS INTO DATAFRAME
        df = DataFrame(XLSX.gettable(xf[sheets])...) #TRANSFORM TABLE INTO DATAFRAME
        data[Symbol(sheets)] = df #DEFINE THE NAME FROM THE WORKSHEET
    end

    V = process_vertex(data[:vertices])
    K = process_vehicle(V,data[:vehicles])

    T = collect( #range from starting month for duration
        range(
            last(data[:periods].start),
            length = last(data[:periods].T),
            step = 1
        )
    )

    dist = JuMP.Containers.DenseAxisArray{Float64}(undef, keys(V), keys(V))
    dist .= 999999999
    for i in keys(V), j in keys(V)
        if i != j
            if f == "haversine"
                dist[i,j] = haversine([V[i].x,V[i].y],[V[j].x,V[j].y],6378.137)
            else #if f == "euclidean"
                dist[i,j] = euclidean([V[i].x,V[i].y],[V[j].x,V[j].y])
            end
        end
    end

    #=d = JuMP.Containers.DenseAxisArray(
        Array{Float64}(data[:non_negative_demands][:,string.(T)]), #dataset
        Array{String}(data[:non_negative_demands].point), #dims 1
        T #dims 2
    )=#

    d = JuMP.Containers.DenseAxisArray(
        Array{Float64}(data[:demands][:,string.(T)]), #dataset
        Array{String}(data[:demands].point), #dims 1
        T #dims 2
    )

    vertex_data[] = V
    vehicle_data[] = K
    period_data[] = T
    demand_data[] = d
    distance_data[] = dist

    return V,K,T,d,dist
end

function process_vertex(df::DataFrame)
    V = Dict{String,vtx}() #INITIATE VERTICES
    for v in eachrow(df) #ITERATE OVER DATA
        V[v.name] = vtx(
            v.x, v.y, v.MAX, v.MIN, v.START, v.h
        )
    end

    return V
end

function process_vehicle(V::Dict,df::DataFrame)
    K = Dict{String,veh}() #INITIATE VEHICLES
    idx_asal = collect(keys(V))
    idx_tujuan = collect(keys(V))
    list_moda = unique(df.Moda)

    for m in list_moda
        trayek = filter(p -> p.Moda == m, df)
        unique!(trayek) #trayek yang unik

        Q = JuMP.Containers.DenseAxisArray{Int64}(undef,idx_asal,idx_tujuan)
        f = JuMP.Containers.DenseAxisArray{Int64}(undef,idx_asal,idx_tujuan)
        g = JuMP.Containers.DenseAxisArray{Int64}(undef,idx_asal,idx_tujuan)
        lim = JuMP.Containers.DenseAxisArray{Int64}(undef,idx_asal,idx_tujuan)

        Q .= 0
        f .= 0
        g .= 0
        lim .= 0

        for r in eachrow(trayek)
            Q[r.Asal,r.Tujuan] = round(r.kapasitas)
            f[r.Asal,r.Tujuan] = round(r.peti_cost)
            g[r.Asal,r.Tujuan] = round(r.trip_cost)
            lim[r.Asal,r.Tujuan] = round(r.limit)
        end

        K[m] = veh(Q,f,g,lim)
    end

    return K
end

const default_optimizer = Ref{Any}(nothing)
set_optimizer!(O) = default_optimizer[] = O
get_optimizer() = default_optimizer[]
reset_optimizer() = default_optimizer[] = nothing


text:
_varconstructor(name::Symbol) = :(Variable{$(esc(Expr(:quote, name)))}())
_makevar(name::Symbol) = :($(esc(name)) = $(_varconstructor(name)))

function _split(expr::Expr)
    @capture(expr, var_Symbol[lb_Int:ub_Int]) || error("Unsupported syntax. Expected x[a:b] for literal integers a and b")
    var, lb, ub
end

function _makevar(expr::Expr)
    var, lb, ub = _split(expr)
    :($(esc(var)) = $(Expr(:tuple, [_varconstructor(Symbol("$var[$i]")) for i in lb:ub]...)))
end

_return_name(s::Symbol) = esc(s)
function _return_name(expr::Expr)
    var, lb, ub = _split(expr)
    esc(var)
end

"""
Construct polynomial variables and bind them to local variables of the same
name. Usage:

Create a single variable named x:

    @polyvar(x)

Create several variables at the same time:

    @polyvar(x, y, z)

Create a tuple of variables x = (x1, x2, x3, x4):

    @polyvar(x[1:4])

Create variables x, y, a = (a1, a2, a3):

    @polyvar(x, y, a[1:3])

You can also assign the results of the macro to a tuple:

    vars = @polyvar(x[1:5])
    @assert typeof(vars[1]) == Variable{:x1}

"""
macro polyvar(names...)
    Expr(:block, _makevar.(names)..., Expr(:tuple, _return_name.(names)...))
end

macro polyvar(name)
    Expr(:block, _makevar(name), _return_name(name))
end


text:
export GradientConfig, GradientDiffResult, config,
    gradient!, evaluate, evaluate!, jacobian, jacobian!, value

const Index = UInt16
const Delimiter = UInt16
const Exponent = UInt16

mutable struct PolyConfig{T}
    monomials_delimiters::Vector{Delimiter}
    monomials::Vector{Index}
    grad_monomials_delimiters::Vector{Vector{Tuple{Delimiter, Exponent}}}
    grad_monomials::Vector{Vector{Index}}
    reduced_exponents_delimiters::Vector{Delimiter}
    reduced_exponents_map::Vector{Tuple{Index, Exponent}}
    reduced_values::Vector{T}
end

function PolyConfig(g::Polynomial{T}, reduced_exponents::Matrix{UInt16}, big_lookup::Matrix{<:Integer}, ::Type{S}) where {T, S}
    exponents = g.exponents
    m, n = size(reduced_exponents)

    reduced_exponents_delimiters = Vector{UInt16}(undef, n)
    reduced_exponents_map = Vector{NTuple{2, UInt16}}()
    for j=1:n
        nindices = 0
        for i=1:m
            if reduced_exponents[i, j] > 0
                nindices += 1
                push!(reduced_exponents_map, (convert(UInt16, i), big_lookup[i, j]))
            end
        end
        reduced_exponents_delimiters[j] = nindices
    end

    # monomials_full = BitMatrix(exponents)
    monomials_delimiters = Vector{UInt16}(undef, n)
    monomials = Vector{UInt16}()
    for j=1:n
        nindices = 0
        for i=1:m
            if exponents[i, j] > 0
                nindices += 1
                push!(monomials, convert(UInt16, i))
            end
        end
        monomials_delimiters[j] = nindices
    end

    grad_monomials_delimiters = Vector{Vector{Tuple{Delimiter, Exponent}}}(undef, m)
    grad_monomials = Vector{Vector{Index}}(undef, m)
    for varindex = 1:m
        imonomials_delimiters = Vector{Tuple{Delimiter, Exponent}}(undef, n)
        imonomials = Vector{Index}()
        for j=1:n
            nindices = 0
            k = exponents[varindex, j]
            if k > 0
                for i=1:m
                    if exponents[i, j] > 0 && i != varindex
                        nindices += 1
                        push!(imonomials, i)
                    end
                end
            end
            imonomials_delimiters[j] = (nindices, k)
        end
        grad_monomials[varindex] = imonomials
        grad_monomials_delimiters[varindex] = imonomials_delimiters
    end

    PolyConfig(
        monomials_delimiters,
        monomials,
        grad_monomials_delimiters,
        grad_monomials,
        reduced_exponents_delimiters,
        reduced_exponents_map,
        zeros(typeof(one(T) * one(S) + one(T) * one(S)), n))
end

function fillreduced_values!(
    cfg::PolyConfig{T},
    g::Polynomial,
    x::AbstractVector,
    diffs_values::AbstractMatrix{T}) where T
    v = cfg.reduced_values
    reds = cfg.reduced_exponents_delimiters
    rem = cfg.reduced_exponents_map
    dv = diffs_values
    cfs = g.coefficients

    N = length(rem)
    n = length(reds)
    k = 1
    j = 1
    nextj_at = j + reds[j]
    @inbounds res = convert(T, cfs[j])
    while k ≤ N || j < n
        togo = nextj_at - k
        if togo == 0
            @inbounds v[j] = res
            j += 1
            @inbounds nextj_at += reds[j]
            @inbounds res = convert(T, cfs[j])
        else
            @inbounds i, l = rem[k]
            @inbounds res *= dv[i, l]
            k += 1
        end
    end
    @inbounds v[j] = res
    v
end

@inline function fillvalues!(
    values::AbstractMatrix{T},
    xs::AbstractVector,
    diffs) where T
    m, n = size(values)
    if n == 0
        return values
    end
    for i=1:m
        @inbounds l = diffs[i,1]
        xi = xs[i]
        v = l == 0 ? one(T) : pow(xi, l)
        @inbounds values[i, 1] = v
        for k=2:n
            @inbounds l = diffs[i,k]
            v = l == 0 ? v : v * pow(xi, l)
            @inbounds values[i,k] = v
        end
    end
    nothing
end

@inline function _evaluate(x::AbstractVector, cfg::PolyConfig{T}) where T
    v = cfg.reduced_values
    mds = cfg.monomials_delimiters
    ms = cfg.monomials

    N = length(ms)
    n = length(mds)
    k = 1
    j = 1
    @inbounds nextj_at = j + mds[j]
    @inbounds res = v[j]
    out = zero(T)
    while k ≤ N || j < n
        if k == nextj_at
            out += res
            j += 1
            @inbounds nextj_at += mds[j]
            @inbounds res = v[j]
        else
            res *= x[ms[k]]
            k += 1
        end
    end
    out += res

    out
end

@inline function _gradient!(u::AbstractVector, x::AbstractVector, cfg::PolyConfig{T}) where T
    v = cfg.reduced_values
    gmds = cfg.grad_monomials_delimiters
    gms = cfg.grad_monomials

    for i=1:length(gmds)
        @inbounds ui = _gradient(v, gmds[i], gms[i], x)
        u[i] = ui
    end
    u
end

"""
    gradient_row!(u::AbstractMatrix, x::AbstractVector, cfg::PolyConfig{T}, i)
Store the gradient in the i-th row of 'u'.
"""
@inline function gradient_row!(u::AbstractMatrix, x::AbstractVector, cfg::PolyConfig{T}, i) where T
    v = cfg.reduced_values
    gmds = cfg.grad_monomials_delimiters
    gms = cfg.grad_monomials

    for j=1:length(gmds)
        @inbounds uij = _gradient(v, gmds[j], gms[j], x)
        u[i, j] = uij
    end
    u
end

function _gradient(v::AbstractVector{T}, mds, ms, x) where T
    N = length(ms)
    n = length(mds)
    k = 1
    j = 1
    @inbounds delim, exponent = mds[j]
    nextj_at = j + delim
    @inbounds res = v[j] * exponent
    out = zero(T)
    while k ≤ N || j < n
        if k == nextj_at
            out += res
            j += 1
            @inbounds delim, exponent = mds[j]
            nextj_at += delim
            @inbounds res = exponent == 0 ? zero(T) : v[j] * exponent
        else
            res *= x[ms[k]]
            k += 1
        end
    end
    out += res
    out
end

"""
    GradientConfig(f::Polynomial{T}, [x::AbstractVector{S}])

A data structure with which the gradient of a `Polynomial` `f` can be
evaluated efficiently. Note that `x` is only used to determine the
output type of `f(x)`.

    GradientConfig(f::Polynomial{T}, [S])

Instead of a vector `x` a type can also be given directly.
"""
mutable struct GradientConfig{T}
    poly::PolyConfig{T}
    differences::Matrix{UInt8}
    differences_values::Matrix{T}
end

function GradientConfig(f::Polynomial{T}, ::AbstractArray{S}) where {T, S}
    GradientConfig(f, S)
end
GradientConfig(f::Polynomial{T}) where T = GradientConfig(f, T)
function GradientConfig(f::Polynomial{T}, ::Type{S}) where {T, S}
    diffs, diffs_values, big_lookup, reduced_exponents = differences(f, S)
    poly = PolyConfig(f, reduced_exponents, big_lookup, S)

    GradientConfig(poly, diffs, diffs_values)
end

"""
    config(F::Polynomial, x)

Construct a `GradientConfig` for the evaluation of `f` with values like `x`.
"""
function config(f::Polynomial{T}, x::AbstractVector{S}) where {S, T}
    GradientConfig(f, typeof(one(T) * one(S) + one(T) * one(S)))
end

function differences(f::Polynomial{T}, ::Type{S}) where {T, S}
    exponents = f.exponents
    reduced_exponents = convert.(UInt16, max.(exponents .- 1, 0))
    differences, big_lookup = computetables(reduced_exponents)
    differences_values = convert.(promote_type(T, S), differences)

    differences, differences_values, big_lookup, reduced_exponents
end

# function Base.deepcopy(cfg::GradientConfig)
#     GradientConfig(
#         deepcopy(cfg.poly),
#         deepcopy(cfg.differences),
#         deepcopy(cfg.differences_values))
# end

"""
    evaluate(g, x, cfg::GradientConfig [, precomputed=false])

Evaluate `g` at `x` using the precomputated values in `cfg`.
Note that this is usually signifcant faster than `evaluate(g, x)`.

### Example
```julia
cfg = GradientConfig(g)
evaluate(g, x, cfg)
```

With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
the result is only correct if you previouls called `evaluate`, or `gradient` with the same
`x`.
"""
function evaluate(g::Polynomial, x::AbstractVector, cfg::GradientConfig{T}, precomputed=false) where T
    if !precomputed
        fillvalues!(cfg.differences_values, x, cfg.differences)
        fillreduced_values!(cfg.poly, g, x, cfg.differences_values)
    end
    _evaluate(x, cfg.poly)
end

"""
    gradient(g, x, cfg::GradientConfig[, precomputed=false])

Compute the gradient of `g` at `x` using the precomputated values in `cfg`.
The return vector is constructed using `similar(x, T)`.

### Example
```julia
cfg = GradientConfig(g)
gradient(g, x, cfg)
```

With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
the result is only correct if you previouls called `evaluate`, or `gradient` with the same
`x`.
"""
function gradient(g::Polynomial, x::AbstractVector, cfg::GradientConfig{T}, precomputed=false) where T
    u = similar(x, T, nvariables(g))
    gradient!(u, g, x, cfg, precomputed)
    u
end

"""
    gradient!(u, g, x, cfg::GradientConfig [, precomputed=false])

Compute the gradient of `g` at `x` using the precomputated values in `cfg`
and store thre result in u.

### Example
```julia
cfg = GradientConfig(g)
gradient(u, g, x, cfg)
```

With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
the result is only correct if you previouls called `evaluate`, or `gradient` with the same
`x`.
"""
function gradient!(u::AbstractVector, g::Polynomial, x::AbstractVector, cfg::GradientConfig{T}, precomputed=false) where T
    if !precomputed
        fillvalues!(cfg.differences_values, x, cfg.differences)
        fillreduced_values!(cfg.poly, g, x, cfg.differences_values)
    end
    _gradient!(u, x, cfg.poly)
end

"""
    GradientDiffResult(cfg::GradientConfig)

During the computation of ``∇g(x)`` we compute nearly everything we need for the evaluation of
``g(x)``. GradientDiffResult allocates memory to hold both values.
This structure also signals `gradient!` to store ``g(x)`` and ``∇g(x)``.

### Example

```julia
cfg = GradientConfig(g, x)
r = GradientDiffResult(cfg)
gradient!(r, g, x, cfg)

value(r) == g(x)
gradient(r) == gradient(g, x, cfg)
```

    GradientDiffResult(grad::AbstractVector)

Allocate the memory to hold the gradient by yourself.
"""
mutable struct GradientDiffResult{T, AV<:AbstractVector{T}}
    value::T
    grad::AV
end

function GradientDiffResult(cfg::GradientConfig{T}) where T
    GradientDiffResult{T, Vector{T}}(zero(T), Vector{T}(undef, size(cfg.differences, 1)))
end
function GradientDiffResult(grad::AbstractVector{T}) where T
    GradientDiffResult{T, Vector{T}}(zero(T), grad)
end

"""
    value(r::GradientDiffResult)

Get the currently stored value in `r`.
"""
value(r::GradientDiffResult) = r.value

"""
    gradient(r::GradientDiffResult)

Get the currently stored gradient in `r`.
"""
gradient(r::GradientDiffResult) = r.grad

"""
    gradient!(r::GradientDiffResult, g, x, cfg::GradientConfig)

Compute ``g(x)`` and the gradient of `g` at `x` at once using the precomputated values in `cfg`
and store thre result in `r`. This is faster than calling both values separetely.

### Example
```julia
cfg = GradientConfig(g)
r = GradientDiffResult(r)
gradient!(r, g, x, cfg)

value(r) == g(x)
gradient(r) == gradient(g, x, cfg)
```
"""
function gradient!(diffresult::GradientDiffResult, g::Polynomial, x::AbstractVector, cfg::GradientConfig{T}) where T
    fillvalues!(cfg.differences_values, x, cfg.differences)
    fillreduced_values!(cfg.poly, g, x, cfg.differences_values)
    diffresult.value = _evaluate(x, cfg.poly)
    _gradient!(diffresult.grad, x, cfg.poly)
    diffresult
end


text:
#
#  y1 +------+
#     |      |
#     |      |
#     |      |
#  y0 +------+
#     x0     x1
#
immutable Rectangle
    x0::Float64
    x1::Float64
    y0::Float64
    y1::Float64
end

if VERSION < v"0.3-"
    Rectangle(x0::Real, x1::Real, y0::Real, y1::Real) =
        Rectangle(@compat Float64(x0), @compat Float64(x1), @compat Float64(y0), @compat Float64(y1))
end

Rectangle() = Rectangle(NaN, NaN, NaN, NaN)
Rectangle(lowerleft::Point, upperright::Point) =
    Rectangle(lowerleft.x, lowerleft.y, upperright.x, upperright.y)
function Rectangle(bb::BoundingBox, xflipped::Bool, yflipped::Bool)
    p = xflipped ? Point(bb.xmax, bb.xmin) : Point(bb.xmin, bb.xmax)
    q = yflipped ? Point(bb.ymax, bb.ymin) : Point(bb.ymin, bb.ymax)
    Rectangle(p, q)
end

BoundingBox(r::Rectangle) = BoundingBox(min(r.x0,r.x1), max(r.x0,r.x1),
                                        min(r.y0,r.y1), max(r.y0,r.y1))

isincomplete(bb::BoundingBox) = isnan(bb.xmin) || isnan(bb.xmax) ||
                                isnan(bb.ymin) || isnan(bb.ymax)

lowerleft(r::Rectangle) = Point(r.x0, r.y0)
upperleft(r::Rectangle) = Point(r.x0, r.y1)
lowerright(r::Rectangle) = Point(r.x1, r.y0)
upperright(r::Rectangle) = Point(r.x1, r.y1)

width(r::Rectangle) = abs(r.x1 - r.x0)
height(r::Rectangle) = abs(r.y1 - r.y0)
diagonal(r::Rectangle) = hypot(r.x1 - r.x0, r.y1 - r.y0)

function deform(r::Rectangle, dl::Real, dr::Real, dt::Real, db::Real)
    Rectangle(r.x0 + dl, r.x1 + dr, r.y0 + dt, r.y1 + db)
end

# shift center by (dx,dy), keeping width & height fixed
function shift(r::Rectangle, dx::Real, dy::Real)
    Rectangle(r.x0 + dx, r.x1 + dx, r.y0 + dy, r.y1 + dy)
end

# scale width & height, keeping center fixed
function (*)(r::Rectangle, s::Real)
    dw = 0.5*(s - 1.)*(r.x1 - r.x0)
    dh = 0.5*(s - 1.)*(r.y1 - r.y0)
    deform(r, -dw, dw, -dh, dh)
end
(*)(s::Real, r::Rectangle) = r*s

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

abstract AbstractProjection1
abstract AbstractProjection2

project(p::AbstractProjection2, pt::Point) = Point(project(p, pt.x, pt.y)...)

immutable LinearProjection <: AbstractProjection1
    a::Float64
    b::Float64
end

project(p::LinearProjection, u) = p.a .+ p.b .* u
deproject(p::LinearProjection, x) = (x .- p.a) ./ p.b

immutable LogProjection <: AbstractProjection1
    a::Float64
    b::Float64
end

project(p::LogProjection, u) = p.a .+ p.b .* log10(u)
deproject(p::LogProjection, x) = 10.0 .^ ((x .- p.a) ./ p.b)

immutable SeparableProjection2{P1<:AbstractProjection1,
                               P2<:AbstractProjection1} <: AbstractProjection2
    x::P1
    y::P2
end

project(p::SeparableProjection2, u, v) = (project(p.x,u), project(p.y,v))
deproject(p::SeparableProjection2, u, v) = (deproject(p.x,u), deproject(p.y,v))

immutable PolarProjection
    x0::Float64
    y0::Float64
    sx::Float64
    sy::Float64
end

function project(p::PolarProjection, r, θ)
    x = p.x0 .+ p.sx .* r .* cos(θ)
    y = p.y0 .+ p.sy .* r .* sin(θ)
    x, y
end

function deproject(p::PolarProjection, x, y)
    r = hypot(x, y)
    θ = atan2(y, x)
    hypot
end

function PlotGeometry(orig::Rectangle, dest::BoundingBox, xlog, ylog)
    x0 = orig.x0
    x1 = orig.x1
    y0 = orig.y0
    y1 = orig.y1

    px = LinearProjection
    py = LinearProjection
    if xlog
        x0 = log10(x0)
        x1 = log10(x1)
        px = LogProjection
    end
    if ylog
        y0 = log10(y0)
        y1 = log10(y1)
        py = LogProjection
    end

    sx = width(dest)/(x1 - x0)
    sy = height(dest)/(y1 - y0)
    tx = lowerleft(dest).x - sx * x0
    ty = lowerleft(dest).y - sy * y0

    SeparableProjection2(px(tx,sx), py(ty,sy))
end

PlotGeometry(orig::Rectangle, dest::BoundingBox) = PlotGeometry(orig, dest, false, false)


text:
# Testing of functions in norms.jl

using NumericExtensions
using Base.Test

x1 = randn(6)
y1 = randn(6)

x2 = randn(5, 6)
y2 = randn(5, 6)

x3 = randn(3, 4, 5)
y3 = randn(3, 4, 5)

xs = Any[x1, x2, x3]
ys = Any[y1, y2, y3]
ps = Any[1, 2, 3, Inf]

tdims = Any[
    Any[1],   # N = 1
    Any[1, 2, (1, 2)], # N = 2
    Any[1, 2, 3, (1, 2), (1, 3), (2, 3), (1, 2, 3)] # N = 3
]

safe_vnorm(x, p) = p == 1 ? sum(abs(x)) :
                   p == 2 ? sqrt(sum(abs2(x))) :
                   p == Inf ? maximum(abs(x)) :
                   sum(abs(x) .^ p) .^ inv(p)

safe_vnorm(x, p, reg) = p == 1 ? sum(abs(x), reg) :
                        p == 2 ? sqrt(sum(abs2(x), reg)) :
                        p == Inf ? maximum(abs(x), reg) :
                        sum(abs(x) .^ p, reg) .^ inv(p)

# vnorm

for x in xs, p in ps
    @test_approx_eq vnorm(x, p) safe_vnorm(x, p)
end

for x in xs, reg in tdims[ndims(x)], p in ps
    @test_approx_eq vnorm(x, p, reg) safe_vnorm(x, p, reg)
end

# vnormdiff

for (x, y) in zip(xs, ys), p in ps
    @test_approx_eq vnormdiff(x, y, p) safe_vnorm(x - y, p)
end

for (x, y) in zip(xs, ys), reg in tdims[ndims(x)], p in ps
    @test_approx_eq vnormdiff(x, y, p, reg) safe_vnorm(x - y, p, reg)
end

# normalize

for x in xs, p in ps
    r = x ./ vnorm(x, p)
    @test_approx_eq normalize(x, p) r

    xc = copy(x)
    normalize!(xc, p)
    @test_approx_eq xc r
end

for x in xs, reg in tdims[ndims(x)], p in ps
    r = x ./ vnorm(x, p, reg)
    @test_approx_eq normalize(x, p, reg) r

    xc = copy(x)
    normalize!(xc, p, reg)
    @test_approx_eq xc r
end