text:
using Test
using LinearAlgebra, SparseArrays
using SummationByPartsOperators

# Accuracy tests with Float32.
let T=Float32
    xmin = one(T)
    xmax = 2*one(T)
    N = 100
    x1 = compute_coefficients(identity, periodic_derivative_operator(1, 2, xmin, xmax, N))

    x0 = fill(one(eltype(x1)), length(x1))
    x2 = x1 .* x1
    x3 = x2 .* x1
    x4 = x2 .* x2
    x5 = x2 .* x3
    x6 = x3 .* x3
    x7 = x4 .* x3

    res = fill(zero(eltype(x0)), length(x0))

    # first derivative operators
    derivative_order = 1
    accuracy_order = 2
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    mul!(res, D, x0)
    @test all(i->abs(res[i]) < eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x1)
    @test all(i->res[i] ≈ x0[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x2)
    @test all(i->res[i] ≈ 2*x1[i], accuracy_order:length(res)-accuracy_order)

    accuracy_order = 4
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    mul!(res, D, x0)
    @test all(i->abs(res[i]) < 50*eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x1)
    @test all(i->res[i] ≈ x0[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x2)
    @test all(i->res[i] ≈ 2*x1[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x3)
    @test all(i->res[i] ≈ 3*x2[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x4)
    @test all(i->res[i] ≈ 4*x3[i], accuracy_order:length(res)-accuracy_order)

    accuracy_order = 6
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    mul!(res, D, x0)
    @test all(i->abs(res[i]) < 50*eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x1)
    @test all(i->res[i] ≈ x0[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x2)
    @test all(i->res[i] ≈ 2*x1[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x3)
    @test all(i->res[i] ≈ 3*x2[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x4)
    @test all(i->res[i] ≈ 4*x3[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x5)
    @test all(i->res[i] ≈ 5*x4[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x6)
    @test all(i->res[i] ≈ 6*x5[i], accuracy_order:length(res)-accuracy_order)

    # second derivative operators
    derivative_order = 2
    accuracy_order = 2
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == true
    mul!(res, D, x0)
    @test all(i->abs(res[i]) < eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x1)
    @test all(i->abs(res[i]) < 100N*eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x2)
    @test all(i->isapprox(res[i], 2*x0[i], atol=610N*eps(T)), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x3)
    @test all(i->isapprox(res[i], 6*x1[i], atol=2000N*eps(T)), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x4)
    @test any(i->!(res[i] ≈ 12*x2[i]), accuracy_order:length(res)-accuracy_order)

    accuracy_order = 4
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == true
    mul!(res, D, x0)
    @test all(i->abs(res[i]) < 50N*eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x1)
    @test all(i->abs(res[i]) < 1000N*eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x2)
    @test all(i->isapprox(res[i], 2*x0[i], atol=5000N*eps(T)), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x3)
    @test all(i->isapprox(res[i], 6*x1[i], atol=5000N*eps(T)), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x4)
    @test all(i->isapprox(res[i], 12*x2[i], atol=8500N*eps(T)), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x5)
    @test any(i->!(res[i] ≈ 30*x3[i]), accuracy_order:length(res)-accuracy_order)

    accuracy_order = 6
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == true
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 50N*eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 1000N*eps(T)
    mul!(res, D, x2); @test norm((res - 2 .* x0)[accuracy_order:end-accuracy_order], Inf) < 5000N*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x1)[accuracy_order:end-accuracy_order], Inf) < 5000N*eps(T)
    mul!(res, D, x4); @test norm((res - 12 .* x2)[accuracy_order:end-accuracy_order], Inf) < 8000N*eps(T)
    mul!(res, D, x5); @test norm((res - 20 .* x3)[accuracy_order:end-accuracy_order], Inf) < 50000N*eps(T)
    mul!(res, D, x6); @test norm((res - 30 .* x4)[accuracy_order:end-accuracy_order], Inf) < 52000N*eps(T)
    mul!(res, D, x7); @test norm((res - 42 .* x5)[accuracy_order:end-accuracy_order], Inf) < 92000N*eps(T)

    tmp = D * x7
    @test typeof(tmp) == typeof(res)
    @test norm(tmp - res, Inf) < eps(T)

    # third derivative operators
    derivative_order = 3
    accuracy_order = 2
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == true # because this operator is zero!
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 1400*eps(T)
    mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 7000*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x0)[accuracy_order:end-accuracy_order], Inf) > 7000*eps(T)

    accuracy_order = 4
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 400*eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 10000N*eps(T)
    mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 50000N*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x0)[accuracy_order:end-accuracy_order], Inf) < 500000N*eps(T)
    mul!(res, D, x4); @test norm((res - 24 .* x1)[accuracy_order:end-accuracy_order], Inf) < 600000N*eps(T)
    mul!(res, D, x5); @test norm((res - 60 .* x2)[accuracy_order:end-accuracy_order], Inf) > 600000*eps(T)

    accuracy_order = 6
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 300*eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 80000N*eps(T)
    mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 300000N*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x0)[accuracy_order:end-accuracy_order], Inf) < 700000N*eps(T)
    mul!(res, D, x4); @test norm((res - 24 .* x1)[accuracy_order:end-accuracy_order], Inf) < 1000000N*eps(T)
    mul!(res, D, x5); @test norm((res - 60 .* x2)[accuracy_order:end-accuracy_order], Inf) < 4000000N*eps(T)
    mul!(res, D, x6); @test norm((res - 120 .* x3)[accuracy_order:end-accuracy_order], Inf) < 8000000N*eps(T)
    mul!(res, D, x7); @test norm((res - 240 .* x4)[accuracy_order:end-accuracy_order], Inf) > 220000N*eps(T)

    tmp = D * x7
    @test typeof(tmp) == typeof(res)
    @test norm(tmp - res, Inf) < eps(T)

    # higher derivative operators
    accuracy_order = 2
    @test_throws ArgumentError periodic_central_derivative_operator(4, accuracy_order, xmin, xmax, N)
    @test_throws ArgumentError periodic_central_derivative_operator(5, accuracy_order, xmin, xmax, N)
    @test_throws ArgumentError periodic_central_derivative_operator(6, accuracy_order, xmin, xmax, N)
end

# Accuracy tests with Float64.
let T = Float64
    xmin = -one(T)
    xmax = 2*one(T)
    N = 100
    x1 = compute_coefficients(identity, periodic_derivative_operator(1, 2, xmin, xmax, N))

    x0 = fill(one(eltype(x1)), length(x1))
    x2 = x1 .* x1
    x3 = x2 .* x1
    x4 = x2 .* x2
    x5 = x2 .* x3
    x6 = x3 .* x3
    x7 = x4 .* x3
    x8 = x4 .* x4
    x9 = x4 .* x5

    res = fill(zero(eltype(x0)), length(x0))

    # first derivative operators
    derivative_order = 1
    accuracy_order = 2
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    mul!(res, D, x0)
    @test all(i->abs(res[i]) < eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x1)
    @test all(i->res[i] ≈ x0[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x2)
    @test all(i->res[i] ≈ 2*x1[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x3)
    @test any(i->!(res[i] ≈ 3*x2[i]), accuracy_order:length(res)-accuracy_order)

    accuracy_order = 4
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    accuracy_order = 4
    D = periodic_central_derivative_operator(1, accuracy_order, xmin, xmax, N)
    @test issymmetric(D) == false
    mul!(res, D, x0)
    @test all(i->abs(res[i]) < 10*eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x1)
    @test all(i->res[i] ≈ x0[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x2)
    @test all(i->res[i] ≈ 2*x1[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x3)
    @test all(i->res[i] ≈ 3*x2[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x4)
    @test all(i->res[i] ≈ 4*x3[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x5)
    @test any(i->!(res[i] ≈ 5*x4[i]), accuracy_order:length(res)-accuracy_order)

    accuracy_order = 6
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    mul!(res, D, x0)
    @test all(i->abs(res[i]) < 10*eps(T), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x1)
    @test all(i->res[i] ≈ x0[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x2)
    @test all(i->res[i] ≈ 2*x1[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x3)
    @test all(i->isapprox(res[i], 3*x2[i], atol=100N*eps(T)), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x4)
    @test all(i->res[i] ≈ 4*x3[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x5)
    @test all(i->isapprox(res[i], 5*x4[i], atol=100N*eps(T)), accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x6)
    @test all(i->res[i] ≈ 6*x5[i], accuracy_order:length(res)-accuracy_order)
    mul!(res, D, x7)
    @test any(i->!(res[i] ≈ 7*x6[i]), accuracy_order:length(res)-accuracy_order)

    tmp = D * x7
    @test typeof(tmp) == typeof(res)
    @test norm(tmp - res, Inf) < eps(T)

    # second derivative operators
    derivative_order = 2
    accuracy_order = 2
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == true
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 1400*eps(T)
    mul!(res, D, x2); @test norm((res - 2 .* x0)[accuracy_order:end-accuracy_order], Inf) < 7000*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x1)[accuracy_order:end-accuracy_order], Inf) > 7000*eps(T)

    accuracy_order = 4
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == true
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 400*eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 5000*eps(T)
    mul!(res, D, x2); @test norm((res - 2 .* x0)[accuracy_order:end-accuracy_order], Inf) < 11000*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x1)[accuracy_order:end-accuracy_order], Inf) < 23000*eps(T)
    mul!(res, D, x4); @test norm((res - 12 .* x2)[accuracy_order:end-accuracy_order], Inf) < 40000*eps(T)
    mul!(res, D, x5); @test norm((res - 20 .* x3)[accuracy_order:end-accuracy_order], Inf) > 50000*eps(T)

    accuracy_order = 6
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == true
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 300*eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 5000*eps(T)
    mul!(res, D, x2); @test norm((res - 2 .* x0)[accuracy_order:end-accuracy_order], Inf) < 10020*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x1)[accuracy_order:end-accuracy_order], Inf) < 30000*eps(T)
    mul!(res, D, x4); @test norm((res - 12 .* x2)[accuracy_order:end-accuracy_order], Inf) < 28000*eps(T)
    mul!(res, D, x5); @test norm((res - 20 .* x3)[accuracy_order:end-accuracy_order], Inf) < 100000*eps(T)
    mul!(res, D, x6); @test norm((res - 30 .* x4)[accuracy_order:end-accuracy_order], Inf) < 220000*eps(T)
    mul!(res, D, x7); @test norm((res - 42 .* x5)[accuracy_order:end-accuracy_order], Inf) > 220000*eps(T)

    tmp = D * x7
    @test typeof(tmp) == typeof(res)
    @test norm(tmp - res, Inf) < eps(T)

    # third derivative operators
    derivative_order = 3
    accuracy_order = 2
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == true # because this operator is zero!
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 1400*eps(T)
    mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 7000*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x0)[accuracy_order:end-accuracy_order], Inf) > 7000*eps(T)

    accuracy_order = 4
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 400*eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 20000*eps(T)
    mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 200000*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x0)[accuracy_order:end-accuracy_order], Inf) < 300000*eps(T)
    mul!(res, D, x4); @test norm((res - 24 .* x1)[accuracy_order:end-accuracy_order], Inf) < 600000*eps(T)
    mul!(res, D, x5); @test norm((res - 60 .* x2)[accuracy_order:end-accuracy_order], Inf) > 600000*eps(T)

    accuracy_order = 6
    D = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    println(devnull, D)
    @test SummationByPartsOperators.derivative_order(D) == derivative_order
    @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
    @test issymmetric(D) == false
    mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 300*eps(T)
    mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 70000*eps(T)
    mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 200000*eps(T)
    mul!(res, D, x3); @test norm((res - 6 .* x0)[accuracy_order:end-accuracy_order], Inf) < 600000*eps(T)
    mul!(res, D, x4); @test norm((res - 24 .* x1)[accuracy_order:end-accuracy_order], Inf) < 1000000*eps(T)
    mul!(res, D, x5); @test norm((res - 60 .* x2)[accuracy_order:end-accuracy_order], Inf) < 4000000*eps(T)
    mul!(res, D, x6); @test norm((res - 120 .* x3)[accuracy_order:end-accuracy_order], Inf) < 8000000*eps(T)
    mul!(res, D, x7); @test norm((res - 240 .* x4)[accuracy_order:end-accuracy_order], Inf) > 220000*eps(T)

    tmp = D * x7
    @test typeof(tmp) == typeof(res)
    @test norm(tmp - res, Inf) < eps(T)

    # higher derivative operators
    accuracy_order = 2
    @test_throws ArgumentError periodic_central_derivative_operator(4, accuracy_order, xmin, xmax, N)
    @test_throws ArgumentError periodic_central_derivative_operator(5, accuracy_order, xmin, xmax, N)
    @test_throws ArgumentError periodic_central_derivative_operator(6, accuracy_order, xmin, xmax, N)
end

# Compare Fornberg algorithm with exact representation of central derivative coefficients.
for T in (Float32, Float64), accuracy_order in 2:2:12, derivative_order in 1:3
    xmin = zero(T)
    xmax = one(T)
    N = 21
    D_central = periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    D_general = periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    @test norm(D_central.coefficients.lower_coef - D_general.coefficients.lower_coef) < 3*eps(T)
    @test abs(D_central.coefficients.central_coef - D_general.coefficients.central_coef) < 3*eps(T)
    @test norm(D_central.coefficients.upper_coef - D_general.coefficients.upper_coef) < 3*eps(T)
end

# Compare serial and threaded with full and sparse matrix vector products implementation.
for T in (Float32, Float64), accuracy_order in 1:10, derivative_order in 1:3
    xmin = zero(T)
    xmax = 5*one(T)
    N = 50
    D_serial = periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    D_threads= periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N, Val{:threads}())

    u = compute_coefficients(x->x^5, D_serial)
    dest_serial = fill(zero(eltype(u)), length(u))
    dest_threads= fill(zero(eltype(u)), length(u))
    dest_full   = fill(zero(eltype(u)), length(u))
    dest_sparse = fill(zero(eltype(u)), length(u))

    D_full = Matrix(D_serial)
    D_sparse = sparse(D_serial)
    mul!(dest_serial, D_serial, u)
    mul!(dest_threads, D_threads, u)
    mul!(dest_full, D_full, u)
    mul!(dest_sparse, D_sparse, u)
    @test all(i->dest_serial[i] ≈ dest_threads[i], eachindex(u))
    @test all(i->isapprox(dest_serial[i], dest_full[i], rtol=5*sqrt(eps(T))), eachindex(u))
    @test all(i->isapprox(dest_serial[i], dest_sparse[i], rtol=5*sqrt(eps(T))), eachindex(u))
end

# Compare mul! with β=0 and mul! without β.
for T in (Float32, Float64), accuracy_order in 1:10, derivative_order in 1:3
    xmin = zero(T)
    xmax = 5*one(T)
    N = 51
    D_serial = periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)
    D_threads= periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N, Val{:threads}())

    u = compute_coefficients(x->x^5, D_serial)
    dest1 = fill(zero(eltype(u)), length(u))
    dest2 = fill(zero(eltype(u)), length(u))

    mul!(dest1, D_serial, u, one(T), zero(T))
    mul!(dest2, D_serial, u, one(T))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
    mul!(dest1, D_threads, u, one(T), zero(T))
    mul!(dest2, D_threads, u, one(T))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
end

# Accuracy tests for periodic signals
let T = Float32
    xmin = T(-1)
    xmax = T(1)
    N = 50
    res = Array{T}(undef, N)

    # order 2
    accuracy_order = 2
    D = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N)

    factor=1; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 1.0e-2
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=2; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.8
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=5; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 1.2
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    # order 4
    accuracy_order = 4
    D = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N)

    factor=1; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 5.0e-5
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=2; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 1.0e-3
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=5; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.09
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=8; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.9
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    # order 6
    accuracy_order = 6
    D = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N)

    factor=1; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 3.0e-6
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=2; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 3.0e-5
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=5; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.008
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=8; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.3
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol
end
let T = Float64
    xmin = T(-1)
    xmax = T(1)
    N = 50
    res = Array{T}(undef, N)

    # order 2
    accuracy_order = 2
    D = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N)

    factor=1; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.01
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=2; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.09
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=5; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 1.2
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    # order 4
    accuracy_order = 4
    D = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N)

    factor=1; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 3.0e-5
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=2; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 1.0e-3
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=5; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.1
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=8; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.9
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    # order 6
    accuracy_order = 6
    D = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N)

    factor=1; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 1.0e-7
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=2; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 3.0e-5
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=5; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.01
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol

    factor=8; ufunc = x->sinpi(factor*x); dufunc = x->factor*T(π)*cospi(factor*x)
    tol = 0.3
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    res = D*u
    @test maximum(abs, res-du) < tol
    xplot, duplot = evaluate_coefficients(res, D)
    @test maximum(abs, duplot-dufunc.(xplot)) < tol
end

# Robust methods of Holoborodko
let T = Float32
    xmin = one(T)
    xmax = 2*one(T)
    N = 100
    x1 = compute_coefficients(identity, periodic_derivative_operator(1, 2, xmin, xmax, N))

    x0 = fill(one(eltype(x1)), length(x1))
    x2 = x1 .* x1
    x3 = x2 .* x1
    x4 = x2 .* x2
    x5 = x2 .* x3
    x6 = x3 .* x3
    x7 = x4 .* x3

    res = fill(zero(eltype(x0)), length(x0))

    # first derivative operators
    der_order = 1
    acc_order = 2
    for stencil_width in 5:2:11
        D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=stencil_width)
        println(devnull, D)
        @test derivative_order(D) == der_order
        @test accuracy_order(D) == acc_order
        @test issymmetric(D) == false
        mul!(res, D, x0)
        @test all(i->abs(res[i]) < eps(T), stencil_width:length(res)-stencil_width)
        mul!(res, D, x1)
        @test all(i->res[i] ≈ x0[i], stencil_width:length(res)-stencil_width)
        mul!(res, D, x2)
        @test all(i->res[i] ≈ 2*x1[i], stencil_width:length(res)-stencil_width)
    end
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=3)
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=13)

    acc_order = 4
    for stencil_width in 7:2:11
        D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=stencil_width)
        println(devnull, D)
        @test derivative_order(D) == der_order
        @test accuracy_order(D) == acc_order
        @test issymmetric(D) == false
        mul!(res, D, x0)
        @test all(i->abs(res[i]) < 50*eps(T),  stencil_width:length(res)-stencil_width)
        mul!(res, D, x1)
        @test all(i->res[i] ≈ x0[i],  stencil_width:length(res)-stencil_width)
        mul!(res, D, x2)
        @test all(i->res[i] ≈ 2*x1[i],  stencil_width:length(res)-stencil_width)
        mul!(res, D, x3)
        @test all(i->res[i] ≈ 3*x2[i],  stencil_width:length(res)-stencil_width)
        mul!(res, D, x4)
        @test all(i->res[i] ≈ 4*x3[i],  stencil_width:length(res)-stencil_width)
    end
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=5)
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=13)

    acc_order = 6
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=9)
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=13)

    # second derivative operators
    der_order = 2
    acc_order = 2
    for stencil_width in 5:2:9
        D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=stencil_width)
        println(devnull, D)
        @test derivative_order(D) == der_order
        @test accuracy_order(D) == acc_order
        @test issymmetric(D) == true
        mul!(res, D, x0)
        @test all(i->abs(res[i]) < eps(T),  stencil_width:length(res)-stencil_width)
        mul!(res, D, x1)
        @test all(i->abs(res[i]) < 100N*eps(T),  stencil_width:length(res)-stencil_width)
        mul!(res, D, x2)
        @test all(i->isapprox(res[i], 2*x0[i], atol=600N*eps(T)),  stencil_width:length(res)-stencil_width)
        mul!(res, D, x3)
        @test all(i->isapprox(res[i], 6*x1[i], atol=2000N*eps(T)),  stencil_width:length(res)-stencil_width)
    end
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=3)
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=11)

    acc_order = 4
    for stencil_width in 7:2:9
        D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=stencil_width)
        println(devnull, D)
        @test derivative_order(D) == der_order
        @test accuracy_order(D) == acc_order
        @test issymmetric(D) == true
        mul!(res, D, x0)
        @test all(i->abs(res[i]) < 50N*eps(T),  stencil_width:length(res)-stencil_width)
        mul!(res, D, x1)
        @test all(i->abs(res[i]) < 1000N*eps(T),  stencil_width:length(res)-stencil_width)
        mul!(res, D, x2)
        @test all(i->isapprox(res[i], 2*x0[i], atol=5000N*eps(T)),  stencil_width:length(res)-stencil_width)
        mul!(res, D, x3)
        @test all(i->isapprox(res[i], 6*x1[i], atol=5000N*eps(T)),  stencil_width:length(res)-stencil_width)
        mul!(res, D, x4)
        @test all(i->isapprox(res[i], 12*x2[i], atol=7000N*eps(T)),  stencil_width:length(res)-stencil_width)
        mul!(res, D, x5)
        @test any(i->!(res[i] ≈ 30*x3[i]),  stencil_width:length(res)-stencil_width)
    end
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=5)
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=11)

    acc_order = 6
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=9)
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=13)

    der_order = 3
    acc_order = 2
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=5)
    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
                                         xmin, xmax, N, stencil_width=13)
end


text:
# Dijkstra's algorithm for single-source shortest path

###################################################################
#
#   The type that capsulates the states of Dijkstra algorithm
#
###################################################################

mutable struct DijkstraStates{V,D<:Number,Heap,H}
    parents::Vector{V}
    parent_indices::Vector{Int}
    dists::Vector{D}
    colormap::Vector{Int}
    heap::Heap
    hmap::Vector{H}
end

struct DijkstraHEntry{V,D}
    vertex::V
    dist::D
end

<(e1::DijkstraHEntry, e2::DijkstraHEntry) = e1.dist < e2.dist

# create Dijkstra states

function create_dijkstra_states(g::AbstractGraph{V}, ::Type{D}) where {V,D<:Number}
    n = num_vertices(g)
    parents = Array{V}(undef, n)
    parent_indices = zeros(Int, n)
    dists = fill(typemax(D), n)
    colormap = zeros(Int, n)
    heap = mutable_binary_minheap(DijkstraHEntry{V,D})
    hmap = zeros(Int, n)

    DijkstraStates(parents, parent_indices, dists, colormap, heap, hmap)
end

###################################################################
#
#   visitors
#
###################################################################

abstract type AbstractDijkstraVisitor end

# invoked when a new vertex is first encountered
discover_vertex!(visitor::AbstractDijkstraVisitor, u, v, d) = nothing

# invoked when the distance of a vertex is determined
# (for each source vertex at the beginning, and when a vertex is popped from the heap)
# returns whether the algorithm should continue
include_vertex!(visitor::AbstractDijkstraVisitor, u, v, d) = true

# invoked when the distance to a vertex is updated (decreased)
update_vertex!(visitor::AbstractDijkstraVisitor, u, v, d) = nothing

# invoked when all neighbors of a vertex has been examined
close_vertex!(visitor::AbstractDijkstraVisitor, v) = nothing

# trivial visitor

mutable struct TrivialDijkstraVisitor <: AbstractDijkstraVisitor
end

# log visitor

mutable struct LogDijkstraVisitor <: AbstractDijkstraVisitor
    io::IO
end

function discover_vertex!(visitor::LogDijkstraVisitor, u, v, d)
    println(visitor.io, "discover vertex $v (parent = $u, dist = $d)")
end

function include_vertex!(visitor::LogDijkstraVisitor, u, v, d)
    println(visitor.io, "include vertex $v (parent = $u, dist = $d)")
    true
end

function update_vertex!(visitor::LogDijkstraVisitor, u, v, d)
    println(visitor.io, "update distance $v (parent = $u, dist = $d)")
end

function close_vertex!(visitor::LogDijkstraVisitor, v)
    println(visitor.io, "close vertex $v")
end

###################################################################
#
#   core algorithm implementation
#
###################################################################

function set_source!(state::DijkstraStates{V,D}, g::AbstractGraph{V}, s::V) where {V,D}
    i = vertex_index(s, g)
    state.parents[i] = s
    state.parent_indices[i] = i
    state.dists[i] = 0
    state.colormap[i] = 2
end

function process_neighbors!(
    state::DijkstraStates{V,D,Heap,H},
    graph::AbstractGraph{V},
    edge_dists::AbstractEdgePropertyInspector{D},
    u::V, du::D, visitor::AbstractDijkstraVisitor) where {V,D,Heap,H}

    dists::Vector{D} = state.dists
    parents::Vector{V} = state.parents
    parent_indices::Vector{Int} = state.parent_indices
    colormap::Vector{Int} = state.colormap
    heap::Heap = state.heap
    hmap::Vector{H} = state.hmap
    dv::D = zero(D)

    for e in out_edges(u, graph)
        v::V = target(e, graph)
        iv::Int = vertex_index(v, graph)
        v_color::Int = colormap[iv]

        if v_color == 0
            dists[iv] = dv = du + edge_property(edge_dists, e, graph)
            parents[iv] = u
            parent_indices[iv] = vertex_index(u, graph)
            colormap[iv] = 1
            discover_vertex!(visitor, u, v, dv)

            # push new vertex to the heap
            hmap[iv] = push!(heap, DijkstraHEntry(v, dv))

        elseif v_color == 1
            dv = du + edge_property(edge_dists, e, graph)
            if dv < dists[iv]
                dists[iv] = dv
                parents[iv] = u
                parent_indices[iv] = vertex_index(u, graph)

                # update the value on the heap
                update_vertex!(visitor, u, v, dv)
                update!(heap, hmap[iv], DijkstraHEntry(v, dv))
            end
        end
    end
end

function dijkstra_shortest_paths!(
    graph::AbstractGraph{V},                # the graph
    edge_dists::AbstractEdgePropertyInspector{D}, # distances associated with edges
    sources::AbstractVector{V},             # the sources
    visitor::AbstractDijkstraVisitor,       # visitor object
    state::DijkstraStates{V,D,Heap,H}) where {V, D, Heap, H}     # the states

    @graph_requires graph incidence_list vertex_map vertex_list
    edge_property_requirement(edge_dists, graph)

    # get state fields

    parents::Vector{V} = state.parents
    dists::Vector{D} = state.dists
    colormap::Vector{Int} = state.colormap
    heap::Heap = state.heap
    hmap::Vector{H} = state.hmap

    # initialize for sources

    d0 = zero(D)

    for s in sources
        set_source!(state, graph, s)
        if !include_vertex!(visitor, s, s, d0)
            return state
        end
    end

    # process direct neighbors of all sources

    for s in sources
        process_neighbors!(state, graph, edge_dists, s, d0, visitor)
        close_vertex!(visitor, s)
    end

    # main loop

    while !isempty(heap)

        # pick next vertex to include
        entry = pop!(heap)
        u::V = entry.vertex
        du::D = entry.dist

        ui = vertex_index(u, graph)
        colormap[ui] = 2
        if !include_vertex!(visitor, parents[ui], u, du)
            return state
        end

        # process u's neighbors

        process_neighbors!(state, graph, edge_dists, u, du, visitor)
        close_vertex!(visitor, u)
    end

    state
end

function dijkstra_shortest_paths(
    graph::AbstractGraph{V},                # the graph
    edge_len::AbstractEdgePropertyInspector{D}, # distances associated with edges
    sources::AbstractVector{V};
    visitor::AbstractDijkstraVisitor=TrivialDijkstraVisitor()) where {V, D}
    #
    state = create_dijkstra_states(graph, D)
    dijkstra_shortest_paths!(graph, edge_len, sources, visitor, state)
end

# Convenient functions

function dijkstra_shortest_paths(
    graph::AbstractGraph{V}, edge_dists::Vector{D}, s::V;
    visitor::AbstractDijkstraVisitor=TrivialDijkstraVisitor()) where {V,D}

    edge_len::AbstractEdgePropertyInspector{D} = VectorEdgePropertyInspector(edge_dists)
    state = create_dijkstra_states(graph, D)
    dijkstra_shortest_paths!(graph, edge_len, [s], visitor, state)
end

function dijkstra_shortest_paths(
    graph::AbstractGraph{V}, edge_dists::Vector{D}, sources::AbstractVector{V};
    visitor::AbstractDijkstraVisitor=TrivialDijkstraVisitor()) where {V,D}

    edge_len::AbstractEdgePropertyInspector{D} = VectorEdgePropertyInspector(edge_dists)
    state = create_dijkstra_states(graph, D)
    dijkstra_shortest_paths!(graph, edge_len, sources, visitor, state)
end

function dijkstra_shortest_paths_withlog(
    graph::AbstractGraph{V}, edge_dists::Vector{D}, s::V) where {V,D}
    #
    dijkstra_shortest_paths(graph, edge_dists, s, visitor=LogDijkstraVisitor(STDOUT))
end

function dijkstra_shortest_paths_withlog(
    graph::AbstractGraph{V}, edge_dists::Vector{D}, sources::AbstractVector{V}) where {V,D}
    #
    dijkstra_shortest_paths(graph, edge_dists, sources, visitor=LogDijkstraVisitor(STDOUT))
end

dijkstra_shortest_paths(
    graph::AbstractGraph{V}, s::V
) where {V} = dijkstra_shortest_paths(graph, ones(num_edges(graph)), s)

function enumerate_indices(parent_indices::Vector{Int}, dest_indices::Vector{Int})
    num_dest = length(dest_indices)
    all_paths = Array{Vector{Int}}(undef, num_dest)
    for i=1:num_dest
        all_paths[i] = Int[]
        index = dest_indices[i]
        if parent_indices[index] != 0
            while parent_indices[index] != index
                push!(all_paths[i], index)
                index = parent_indices[index]
            end
            push!(all_paths[i], index)
            reverse!(all_paths[i])
        end
    end
    all_paths
end

enumerate_indices(parent_indices::Vector{Int}, dest_index::Int) = enumerate_indices(parent_indices, Int[dest_index])[1]
enumerate_indices(parent_indices::Vector{Int}) = enumerate_indices(parent_indices, collect(1:length(parent_indices)))
enumerate_paths(vertices, parent_indices::Vector{Int}, dest_indices::Vector{Int}) = [vertices[i] for i in enumerate_indices(parent_indices, dest_indices)]
enumerate_paths(vertices, parent_indices::Vector{Int}, dest_index::Int) = enumerate_paths(vertices, parent_indices, [dest_index])[1]
enumerate_paths(vertices, parent_indices::Vector{Int}) = enumerate_paths(vertices, parent_indices, collect(1:length(parent_indices)))

ext: jl


text:
@testset "laguerre_horner" begin
    tol = 1e-10

    for α in 0:50, x in 0:0.1:10
        @test isapprox(laguerre_horner(0, α)(x), 1, atol=tol)
        @test isapprox(
            laguerre_horner(1, α)(x), -x + α + 1, atol=tol)
        @test isapprox(
            laguerre_horner(2, α)(x),
            x^2 / 2 -
            ((α + 2) * x) +
            ((α + 2) * (α + 1)) / 2,
            atol=tol
        )
        @test isapprox(
            laguerre_horner(3, α)(x),
            -x^3 / 6 +
            ((α + 3) * x^2) / 2 -
            ((α + 2) * (α + 3) * x) / 2 +
            ((α + 1) * (α + 2) * (α + 3)) / 6,
            atol=tol
        )
    end
end

@testset "laguerre" begin
    for n in 0:5, α in 0:5, x in 0:0.1:10
        # the precision of Horner's method is pretty bad when n is large
        thm_val = laguerre(n, α)(x)
        test_val = laguerre_horner(n, α)(x)
        @test abs(test_val - thm_val) < 1e-11
    end

    m_range = n_range = 1:35
    x_range = -1:0.1:1
    p_range = -0.6:0.1:0.6

    # n >= m
    for m in m_range, n in m:n_range.stop, x in x_range, p in p_range
        @test abs(laguerre(m, n, x, p) - laguerre(m, n-m, 2(x^2 + p^2))) < 1e-5
    end
    # n < m
    for n in n_range, m in n:m_range.stop, x in x_range, p in p_range
        @test abs(laguerre(m, n, x, p) - laguerre(n, m-n, 2(x^2 + p^2))) < 1e-5
    end

    ms = collect(m_range)
    ns = collect(n_range)
    xs = collect(x_range)
    ps = collect(p_range)
    @test laguerre(ms, ns)(xs, ps) == laguerre(ms, ns, xs, ps)
    @test laguerre(xs, ps)(ms, ns) == laguerre(ms, ns, xs, ps)
end


text:
"""
    Writers

Abstracts writing dimensioned data so that output can be to a NetCDF
file or another file format.
"""

module Writers

export AbstractWriter, NetCDFWriter, full_name, init_data, append_data

abstract type AbstractWriter end

"""
    full_name(
        writer,
        filename,
    )

Appends the appropriate (based on `writer`) extension to the specified
filename.
"""
function full_name end

"""
    init_data(
        writer,
        filename,
        no_overwrite,
        dims,
        vars,
    )

Creates the specified file, initializing it with the specified dimension
information. An unlimited `time` dimension is implicitly created. The
specified variables are also defined. This function must be called before
`append_data()`. Specialized by every `Writer` subtype.

# Arguments:
# - `writer`: instance of a subtype of `AbstractWriter`.
# - `filename`: into which to write data (without extension).
# - `no_overwrite`: if `true`, throw an error if `filename` exists and
#   would be overwritten.
# - `dims`: Dict of dimension name to 2-tuple of dimension values and Dict
#   of attributes.
# - `vars`: Dict of variable name to 3-tuple of a k-tuple of dimension
#   names, variable type, and Dict of attributes.
"""
function init_data end

"""
    append_data(
        writer,
        filename,
        varvals,
        simtime,
    )

Appends the specified variables to the specified file. The file must have
been previously created with `init_data()`. `simtime` is appended to the
`time` dimension variable. Specialized by every `Writer` subtype.

# Arguments:
# - `writer`: instance of a subtype of `AbstractWriter`.
# - `filename`: into which to write data (without extension).
# - `varvals`: Dict of variable name to k-dimensional array of values.
# - `simtime`: Current simulation time.
"""
function append_data end

include("netcdf_writer.jl")

end


text:
@testset "Check identifiability from given IO-equations" begin
    R, (x, y, z) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "z"])

    println(simplify_field_generators([[R(1), x + y, x^2 + y^2]]))

    println(simplify_field_generators([[R(1), (x + y)^2, x^3 + y^3, x^2 + y^2]]))
end


text:
using Distributions, SpecialFunctions, Random

# Functions for Gaussian task simulations (Section 2.3.1)

abstract type GLearner end

# Run one learner on one Gaussian environment
function run!(envtype::Gaussian, learnertype, learnerparam;
                    sigma = 0.1, changeprobability = 0.001,
                    nsteps = 10^3, iseedenvpolicy = 3, iseedlearner = 3)
    env, learner, callback = initialize(envtype, learnertype, learnerparam,
                                sigma = sigma,
                                changeprobability = changeprobability,
                                nsteps = nsteps,
                                iseedenvpolicy = iseedenvpolicy,
                                iseedlearner = iseedlearner)
    learn!(learner, env, callback)
    callback
end

# Initialize one Gaussian learner and one Gaussian environment
function initialize(envtype::Gaussian, learnertype, learnerparameterset;
                      sigma = 0.1, changeprobability = 0.001,
                      nsteps = 10^3, iseedenvpolicy = 3,
                      iseedlearner = 3)
    env = setupenv(envtype, sigma = sigma,
                      changeprobability = changeprobability,
                      nsteps = nsteps,
                      iseedenvpolicy = iseedenvpolicy)
    learner = setuplearner(envtype, learnertype, learnerparameterset,
                      iseedlearner = iseedlearner)
    callback = CallbackGaussian()
    env, learner, callback
end

# Set up a Gaussian environment
function setupenv(envtype::Gaussian; sigma = 0.1, changeprobability = 0.001,
                    nsteps = 10^3, iseedenvpolicy = 3)
    Random.seed!(ENV_RNG, iseedenvpolicy)
    env = EnvGaussian(sigma = sigma,
                    changeprobability = changeprobability,
                    nsteps = nsteps,
                    seed = iseedenvpolicy)
end

# Set up a Gaussian learner
function setuplearner(envtype::Gaussian, learnertype, learnerparameterset;
                        iseedlearner = 3)
    if learnertype == GParticleFilter
        learner = GParticleFilter(nparticles = learnerparameterset[1],
                                sigma = learnerparameterset[2],
                                changeprobability = learnerparameterset[3],
                                seed = iseedlearner,
                                Neffthrs = learnerparameterset[1]/2.) # 10
    elseif learnertype == GLeakyIntegrator
        learner = GLeakyIntegrator(etaleak = learnerparameterset[1])
    elseif learnertype == GSmileOriginal
        learner = GSmileOriginal(m = learnerparameterset[1],
                sigma = learnerparameterset[2])
    elseif learnertype == GVarSmile
        if length(learnerparameterset) == 2
            learner = GVarSmile(m = learnerparameterset[1],
                        sigma = learnerparameterset[2])
        elseif length(learnerparameterset) == 3 # Used for calculation of Shannon for Experimental Prediction
            learner = GVarSmile(m = learnerparameterset[1],
                                    sigma = learnerparameterset[2],
                                    changeprobability = learnerparameterset[3])
        end
    elseif learnertype == GBayesFilter
        learner = GBayesFilter(sigma = learnerparameterset[1],
                changeprobability = learnerparameterset[2],
                cutoff = 1e-16)
    elseif learnertype == GMPN
        learner = GMPN(N = learnerparameterset[1],
                    sigma = learnerparameterset[2],
                    changeprobability = learnerparameterset[3])
    elseif learnertype == GSOR
        learner = GSOR(N = learnerparameterset[1],
                    sigma = learnerparameterset[2],
                    changeprobability = learnerparameterset[3],
                    seed = iseedlearner)
    elseif learnertype == GNas10
        learner = GNas10(sigma = learnerparameterset[1],
                            changeprobability = learnerparameterset[2])
    elseif learnertype == GNas12
        learner = GNas12(sigma = learnerparameterset[1],
                            changeprobability = learnerparameterset[2])
    elseif learnertype == GNas10original
        learner = GNas10original(sigma = learnerparameterset[1],
                            changeprobability = learnerparameterset[2])
    elseif learnertype == GNas12original
        learner = GNas12original(sigma = learnerparameterset[1],
                            changeprobability = learnerparameterset[2])
    end
    learner
end


text:
using MPIDistributedArrays
using Test
using MPI

testdir = @__DIR__
istest(f) = endswith(f, ".jl") && startswith(f, "test_")
testfiles = sort(filter(istest, readdir(testdir)))

@testset "$f" for f in testfiles
#@testset "MPIDistributedArrays.jl" begin
    mpiexec() do cmd
        println("nproc = 1")
        run(`$cmd -n 1 $(Base.julia_cmd()) $(joinpath(testdir, f))`)
        @test true
        println("nproc = 2")
        run(`$cmd -n 2 $(Base.julia_cmd()) $(joinpath(testdir, f))`)
        @test true
    end
end


text:
"""
High-order singular value decomposition (HO-SVD).
"""
function hosvd{T,N}(tensor::StridedArray{T,N}, core_dims::NTuple{N, Int};
                    pad_zeros::Bool=false, compute_error::Bool=false)
    pad_zeros || _check_tensor(tensor, core_dims)

    factors = map(1:N) do i
        X = _col_unfold(tensor, i)
        f = eigs(X'X, nev=core_dims[i])[2]
        if pad_zeros && size(f, 2) < core_dims[i] # fill missing factors with zeros
            f = hcat(f, zeros(T, size(tensor, i), core_dims[i]-size(f, 2)))
        end
        mapslices(_check_sign, f, 1)
    end

    res = Tucker((factors...), tensorcontractmatrices(tensor, factors))
    if compute_error
        _set_rel_residue(res, tensor)
    end
    return res
end

hosvd{T,N}(tensor::StridedArray{T,N}, r::Int; compute_error::Bool=false) =
    hosvd(tensor, (fill(r, N)...); compute_error=compute_error)


text:
# basic structs, e.g. Ray, Material, HitRecord, Scatter, Sphere
struct Ray{T}
	origin::Vec3{T} # Point
	dir::Vec3{T} # Vec3 # direction (unit vector)
end

"An object that can be hit by Ray"
abstract type Hittable end

const HittableList = Vector{Hittable}

"""Materials tell us how rays interact with a surface"""
abstract type Material{T <: AbstractFloat} end

"Record a hit between a ray and an object's surface"
struct HitRecord{T <: AbstractFloat}
	t::T # distance from the ray's origin to the intersection with a surface.

	p::Vec3{T} # point of the intersection between an object's surface and a ray
	n⃗::Vec3{T} # surface's outward normal vector, points towards outside of object?

	# If true, our ray hit from outside to the front of the surface.
	# If false, the ray hit from within.
	front_face::Bool

	mat::Material{T}

	@inline HitRecord(t::T,p,n⃗,front_face,mat) where T = new{T}(t,p,n⃗,front_face,mat)
end

struct Sphere{T <: AbstractFloat} <: Hittable
	center::Vec3{T}
	radius::T
	mat::Material{T}
end

struct Scatter{T<: AbstractFloat}
	r::Ray{T}
	attenuation::Vec3{T}

	# claforte: TODO: rename to "absorbed?", i.e. not reflected/refracted?
	reflected::Bool # whether the scattered ray was reflected, or fully absorbed
	@inline Scatter(r::Ray{T},a::Vec3{T},reflected=true) where T = new{T}(r,a,reflected)
end


text:
sdt = sdt_indices(6, 7, 8, 9, adjusted=false)
@test sdt["dprime"] ≈ -0.023 atol=0.001
@test sdt["aprime"] ≈ 0.491 atol=0.001
@test sdt["beta"] ≈ 0.996 atol=0.001
@test sdt["c"] ≈ 0.169 atol=0.001
@test sdt["c_relative"] ≈ -7.43 atol=0.001
@test sdt["Xc"] ≈ 0.157 atol=0.001
@test sdt["bpp"] ≈ 0.263 atol=0.001

@test sdt["dprime_glm"] ≈ -0.0227 atol=0.001
@test sdt["Xc_glm"] ≈ 0.157 atol=0.001

sdt = sdt_indices(6, 7, 8, 9, adjusted=true)
@test sdt["dprime"] ≈ -0.023 atol=0.001
@test sdt["beta"] ≈ 0.995 atol=0.001

sdt = sdt_indices.([9, 1], [1, 9], [1, 9], [9, 1], adjusted=false)
@test length(sdt) == 2
@test sdt[1]["aprime"] ≈ 0.944 atol=0.001
@test sdt[1]["pr"] ≈ 0.800 atol=0.001
@test sdt[1]["br"] ≈ 0.500 atol=0.001
@test sdt[1]["hit_rate"] ≈ 0.900 atol=0.001
@test sdt[1]["fa_rate"] ≈ 0.100 atol=0.001


text:
@diffusion_process SineAux begin
    :parameters
    (a, b, c, σ, t, T, u, v) --> Float64

    :additional
    constdiff --> true
end

B(t, P::SinDiffusionAux) = @SMatrix [0 + (t/P.T)/5]
β(t, P::SinDiffusionAux) = ℝ{1}((P.v-P.u)/P.T)

σ(t, x, P::Sine) = ℝ{1}(P.σ)


text:
using Test
using GWAlpha

geno_sync_fname = "test.sync"
pheno_py_fname = "test.py"
pheno_csv_fname = "test.csv"

function runGWAlpha(filename_sync::String, filename_phen::String, MAF::Float64, DEPTH::Int64, MODEL::String)
    GWAlpha.PoolGPAS(filename_sync=filename_sync, filename_phen=filename_phen, maf=MAF, depth=DEPTH, model=MODEL, plot=true)
    return(0)
end

### GWAlpha for GWAS (non-parallel)
Test.@test runGWAlpha(geno_sync_fname,
                      pheno_py_fname,
                      0.01,
                      10,
                      "GWAlpha") == 0

### Linear mixed model (default: MLE)
Test.@test runGWAlpha(geno_sync_fname,
                      pheno_csv_fname,
                      0.01,
                      10,
                      "MIXED") == 0

### Linear mixed model (default: RR)
Test.@test runGWAlpha(geno_sync_fname,
                    pheno_csv_fname,
                    0.01,
                    10,
                    "GLMNET") == 0

### GWAlpha for GWAS (parallel)
using Distributed
Distributed.addprocs(length(Sys.cpu_info()))
@everywhere using GWAlpha
function runGWAlpha(filename_sync::String, filename_phen::String, MAF::Float64, DEPTH::Int64, MODEL::String)
    GWAlpha.PoolGPAS(filename_sync=filename_sync, filename_phen=filename_phen, maf=MAF, depth=DEPTH, model=MODEL, plot=true)
    return(0)
end
Test.@test runGWAlpha(geno_sync_fname,
                      pheno_py_fname,
                      0.01,
                      10,
                      "GWAlpha") == 0

max_stars_repo_name: jeffersonfparil/GWAlpha


text:
function identity!(args...); end;

"""
    generate_random_input(n [, nrep=5][; σ=1, seed=1]) -> u

Generate a vector ``u`` containing a random signal. The values are draw from a
zero-mean Gaussian distribution with standard deviation `σ` and held for `nrep`
samples. The random number generator is initialized with `seed`, during the
generation.
"""
function generate_random_input(n, nrep=5; σ=1, seed=1)
    srand(seed)
    u = σ*randn(Int(n//nrep));
    return repeat(u; inner=(nrep));
end

## Lag function
function lag!(result::Vector{Float64},
              x::Vector{Float64}, n::Int, N::Int)
    for i = 1:n
        @inbounds result[i] = x[1]
    end
    for i = (n+1):N
        @inbounds result[i] = x[i-n]
    end
end

function lag(x::Vector{Float64}, n::Int)
    N = length(x)
    result = copy(x)
    lag!(result, x, n, N)
    return result
end

function lag(x::Vector{Float64}, list_n::Vector{Int})
    hcat((lag(x, n) for n in list_n)...)
end

function lag(x::Matrix{Float64}, llist_n::Vector{Vector{Int}})
    hcat((lag(x[:, i], list_n) for (i, list_n) in enumerate(llist_n))...)
end

## matrix2subviews
function matrix2subviews{T<:Any}(X::Matrix{T})
    N, p = size(X)
    X_views = Vector{Vector{T}}(p)
    for i = 1:p
        X_views[i] = X[:, i]
    end
    return X_views, N, p
end

function subviews2matrix{T<:Any}(X_views::Vector{Vector{T}}, N, p)
    X = Matrix{T}(N, p)
    for i = 1:p
        for j = 1:N
            X[j, i] = X_views[i][j]
        end
    end
    return X
end

# Convert Compressed Predictor Matrix to matrix
if glmnet_loaded
    function Base.convert(::Type{Matrix{Float64}},
                          X::GLMNet.CompressedPredictorMatrix)
        mat = zeros(X.ni, length(X.nin))
        for b = 1:size(mat, 2), i = 1:X.nin[b]
            mat[X.ia[i], b] = X.ca[i, b]
        end
        return mat
    end
end


text:
#=
    Contents: Default HTML forms for access control.
=#

"""
Returns: default login form.
Ensures HTTP method is POST.
"""
function login_form()
    "<form onsubmit='login(); return false;' method='post'>
	 Username:

	 <input type='text' id='username' name='username'/>

	 Password:

	 <input type='password' id='password' name='password'/>

	 <input type='submit' value='Login'/>
     </form>
     "
end

function logout_link()
    "

Logout

"
end

function logout_pwdreset_links()
    "

Logout

Reset Password

"
end

"Returns: default pwdreset form for server-side session."
function pwdreset_form(session_id::AbstractString)
    form = pwdreset_form_no_nonce()
    add_nonce_to_form("pwdreset", form, session_id)
end

"Returns: default pwdreset form for client-side session."
function pwdreset_form(session::Dict)
    form = pwdreset_form_no_nonce()
    add_nonce_to_form("pwdreset", form, session)
end

function pwdreset_form_no_nonce()
    "
Password reset.

     <form id='pwdreset' onsubmit='pwdreset(); return false;' method='post'>
         Current password:

         <input type='password' id='current_pwd' name='current_pwd'/>

	 New password:

	 <input type='password' id='new_pwd' name='new_pwd'/>

	 Retype new password:

	 <input type='password' id='new_pwd2' name='new_pwd2'/>

         <input type='submit' value='Reset Password'/>
     </form>
     "
end

# EOF


text:
module Kafka

export
    KafkaClient,
    metadata,
    produce,
    fetch,
    list_offsets,
    _metadata,
    _produce,
    _fetch,
    _list_offsets,
    earliest_offset,
    latest_offset,
    api_versions

include("core.jl")

end # module


text:
@testset "650.2-keys-keyboard.jl" begin
    for (i, j) in
        zip([3, 480, 81, 758, 484, 656, 537, 844, 260, 361, 778, 638, 905, 619, 9],
            [3, 18, 12, 381, 26, 49, 182, 215, 22, 38, 391, 42, 186, 619, 6])
        @test min_steps650(i) == j
    end
end


text:
a(x) = (println("\t# Called a($x)"); return x)
b(x) = (println("\t# Called b($x)"); return x)

for i in [true,false], j in [true, false]
    println("\nCalculating: x = a($i) && b($j)"); x = a(i) && b(j)
    println("\tResult: x = $x")
    println("\nCalculating: y = a($i) || b($j)"); y = a(i) || b(j)
    println("\tResult: y = $y")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using EMSimJulia
using Documenter

DocMeta.setdocmeta!(EMSimJulia, :DocTestSetup, :(using EMSimJulia); recursive=true)

makedocs(;
    modules=[EMSimJulia],
    authors="Mohamed Kamal AbdElrahman",
    repo="https://github.com/MKAbdElrahman/EMSimJulia.jl/blob/{commit}{path}#{line}",
    sitename="EMSimJulia.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MKAbdElrahman.github.io/EMSimJulia.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/MKAbdElrahman/EMSimJulia.jl",
    devbranch="main",
)


text:
using Documenter, GeoData, HDF5, ArchGDAL, NCDatasets

makedocs(
    modules = [GeoData],
    sitename = "GeoData.jl",
)

deploydocs(
    repo = "github.com/rafaqz/GeoData.jl.git",
)


text:
__precompile__()
module HydroModels

# Standard library
using LinearAlgebra
using DelimitedFiles
using Dates
using Statistics
using Printf

# Packages
using JuMP
using StochasticPrograms
using GLPKMathProgInterface
using RecipesBase
using Parameters
using Reexport
using MacroTools
using MacroTools: postwalk, @q
using DelimitedFiles
using Flux
using BSON
using BSON: @save, @load

import Base: show
import Statistics: mean

export
    AbstractModelIndices,
    AbstractModelData,
    Horizon,
    Resolution,
    PlantCollection,
    plants,
    nplants,
    Skellefteälven,
    HydroPlantCollection,
    hydrodata,
    plants_in_river,
    plants_in_area,
    plants_in_areas_and_rivers,
    PolyhedralWaterValue,
    cuts,
    ncuts,
    cut_indices,
    nindices,
    cut_lb,
    PriceCurve,
    PriceData,
    Inflows,
    InflowSequence,
    local_inflows,
    local_inflow_sequence,
    Forecaster,
    PriceForecaster,
    FlowForecaster,
    forecast,
    mean_price,
    reload!,
    reinitialize!,
    plan!,
    status,
    production,
    plants,
    discharge,
    spillage,
    reservoir,
    power,
    revenue,
    totalrevenue,
    num_hours,
    num_days,
    num_weeks,
    num_periods,
    water_volume,
    marginal_production,
    mean_price,
    mean_flow,
    water_flow_time,
    historic_flow,
    overflow,
    DayAhead,
    WeekAhead

macro exportSPjl()
    Expr(:export, names(StochasticPrograms)...)
end
@exportSPjl

# Include files
include("data/data.jl")
include("models/model.jl")

# Models
#include("models/short-term/ShortTerm.jl")
include("models/day-ahead/DayAhead.jl")
include("models/week-ahead/WeekAhead.jl")
include("models/maintenance_scheduling/MaintenanceScheduling.jl")
include("models/capacity_expansion/CapacityExpansion.jl")

# Analysis
include("productionplan.jl")

end # module


text:
module TestBasic
using Test
using RData

@testset "Basic/Utility RData functions" begin

# check for Float64 NA
@testset "Detect R floating-point NAs" begin
    @test !RData.isna_float64(reinterpret(UInt64, 1.0))
    @test !RData.isna(1.0)
    @test !RData.isna(NaN)
    @test !RData.isna(Inf)
    @test !RData.isna(-Inf)
    @test RData.isna_float64(RData.R_NA_FLOAT64)
    # check that alternative NA is also recognized (#10)
    @test RData.isna_float64(reinterpret(UInt64, RData.R_NA_FLOAT64 | ((Base.significand_mask(Float64) + 1) >> 1)))
end

@testset "Unescape R strings" begin
    @test RData.unescape_rstring("\\\"") == "\""
    @test RData.unescape_rstring("\\'") == "'"
end

end

end # TestBasic


text:
export AbstractVar, Var, SubseqVar, Ignore
export @V_str
export AbstractBindings, EmptyBindings, Bindings
export lookupfirst, lookupall, lookupequiv, lookup
export ubind, toDict

abstract type AbstractVar end

"""
A Unification vatiable.
"""
struct Var <: AbstractVar
    name::Symbol

    Var(name::AbstractString) = Var(Symbol(name))
    Var(name::Symbol) = new(name)
end

# Do we need to include some notion of scope?
function same(var1::AbstractVar, var2::AbstractVar)::Bool
    typeof(var1) == typeof(var2) &&
        var1.name == var2.name
end

same(var1::AbstractVar, ::Any) = false
same(::Any, Var1::AbstractVar) = false

"""
A Unification variable for subsequences.
The associated value will be a SubArray.
The presence of more than one SubseqVar in
the same sequence allows for multiple
unifications, and combiniatorial explosion
in finding them.
"""
struct SubseqVar <: AbstractVar
    name::Symbol

    SubseqVar(name::AbstractString) = SubseqVar(Symbol(name))
    SubseqVar(name::Symbol) = new(name)
end

"""
Ignore match3es anything but captures nothing.
"""
struct Ignore end

macro V_str(name)
    @assert name isa AbstractString
    if length(name) == 0
        :(Ignore())
    elseif endswith(name, "...")
        :(SubseqVar($name))
    else
        :(Var($name))
    end
end

abstract type AbstractBindings end

"""
The endmost tail of a bindings chain.
"""
struct EmptyBindings <: AbstractBindings
end

Base.iterate(::EmptyBindings) = nothing
Base.iterate(::AbstractBindings, ::EmptyBindings) = nothing

Base.IteratorSize(::Type{<:AbstractBindings}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:AbstractBindings}) = Tuple{AbstractVar, Any}

struct Bindings <: AbstractBindings
    var::AbstractVar
    val::Any
    tail::AbstractBindings

    Bindings(var::AbstractVar, val::Any, tail::AbstractBindings) =
        new(var, val, tail)
end

Base.iterate(b::Bindings, state::Bindings=b) = (state.var, state.val), state.tail

"""
    lookupfirst(::AbstractBindings, ::AbstractVar)
Return the `val` of the first entry in the `Bindings` about the`AbstractVar`.
"""
function lookupfirst(bindings::AbstractBindings, var::AbstractVar)
    for (v, val) in bindings
        @debug "$v, $val"
        if same(v, var)
            return val, true
        end
        # Binding a var to another var is symetric:
        if same(val, var)
            return v, true
        end
    end
    return nothing, false
end

"""
    lookupall(::Bindings, ::AbstractVar)::Set{Any}
Return a `Vector` of all of the `val`s in the `Bindings` which
correspond with the specified `AbstractVar`.
"""
function lookupall(bindings::AbstractBindings, var::AbstractVar)::Set{Any}
    found = Set{Any}()
    for (v, val) in bindings
        if same(v, var)
            push!(found, val)
        elseif same(val, var)
            # Binding a variable to a variable is symetric:
            push!(found, v)
        end
    end
    return found
end

"""
    lookupequiv(::Bindings, ::AbstractVar)
Return a `Set` of all of the values associated with the `AbstractVar`
and another Set of all equivalent (transitively bound) `AbstractVar`s.
"""
function lookupequiv(bindings::AbstractBindings, var::AbstractVar)
    found = Set{Any}()
    queue = Vector{AbstractVar}()
    push!(queue, var)
    done = Set{AbstractVar}()
    while !isempty(queue)
        var = pop!(queue)
        if var in done    ##### Probably using == rather than same
            continue
        end
        push!(done, var)
        vals = lookupall(bindings, var)
        # With Set I don't think we can skip over the ones already returned.
        for val in vals
            if val isa AbstractVar
                push!(queue, val)
            else
                push!(found, val)
            end
        end
    end
    return found, done
end

"""
    lookup(bindings::AbstractBindings, var::AbstractVar)
Return the value assoiciated with `var` in `bindings`.
The second return value is `true` if a unique value is found
and `false` otherwise.
"""
function lookup(bindings::AbstractBindings, var::AbstractVar)
    values, _ = lookupequiv(bindings, var)
    if length(values) == 1
        return first(values), true
    end
    return nothing, false
end

"""
    ubind(continuation, var::AbstractVar, val::Any, [::AbstractBindings])
Call `continuation` with the binding of `var` to `val` added
to `bindings`.
"""
function ubind(continuation, var::AbstractVar, value::Any,
               # Maybe we shouldn't default bindings.  FOrgetting to
               # include bindings in a call to ubind seems to be a
               # common source of errors.
               bindings::AbstractBindings = EmptyBindings())
    continuation(Bindings(var, value, bindings))
end

"""
    ubind(continuation, pairs::Vector{Pair{<:AbstractVar, <:Any}}, ::AbstractBindings=EmptyBindings())
Call `continuation` on a chain of `Bindings` where the `Pair`s are added to
the head of the chain such that the first pair is at the
head of the chain.
"""
function ubind(continuation, pairs::Vector{<:Pair{<:AbstractVar, <:Any}},
               bindings::AbstractBindings=EmptyBindings())
    for i = lastindex(pairs) : -1 :firstindex(pairs)
        p = pairs[i]
        bindings = Bindings(p.first, p.second, bindings)
    end
    continuation(bindings)
end

#=
function toDict(bindings::AbstractBindings)::Dict{AbstractVar, Any}
    d = Dict{AbstractVar, Any}()
    b = bindings
    while true
        if haskey(d, b.var)
            if d[b.var] != b.val
                throw(ErrorException(
                    "More than one distinct value for $(b.var): $(d[b.var]), $(b.val)"))
            end
        else
            d[b.var] = b.val
        end
    end
    return d
end
=#


text:
aE(E::Float64,l::Int,r::Float64,v::Float64) = 2*(E - v - l*(l+1)/(2r*r))

aE(E::Float64,l::Int,r::AbstractVector{Float64},v::AbstractVector{Float64}) =
    Float64[aE(E,l,r[i],v[i]) for i in 1:length(r)]

aElog(E::Float64,l::Int,r::Float64,v::Float64) = 2r*r*(E-v) - (l+0.5)^2

aElog(E::Float64,l::Int,r::AbstractVector{Float64},v::AbstractVector{Float64}) =
    Float64[aElog(E,l,r[i],v[i]) for i in 1:length(r)]

doc"""
`radschrod1(r,v,n,l;reverse=false,Emin=-1.0,Emax=0.0,Elim=1e4)`

solves the radial Schrödinger equation

$\left\{ -{1\over 2}{d^2\over dr^2} + {l(l+1)\over 2r^2} + v(r) \right\} u(r) = \epsilon u(r)$

where $\varphi(r,\Omega) = {u(r)\over r} Y_l^m(\Omega)$.

for a single eigenvalue/eigenfunction, specified by `n,l`.

the eigenfunction is normalized such that $\int dr\, u(r)^2 = 1$

`r` - uniform grid (`lingrid`, `simpsgrid`) or logarithmic grid (`loggrid`)

`v` - the radial potential $v(r)$ evaluated on the grid `r`

`reverse` - if true, integrate in reverse.

`Emin,Emax` - initial `E` interval to search. this will be expanded automatically
 (up to `±Elim`) until the range includes all requested eigenvalues.

returns `ev,ef`, the requested eigenvalue $\epsilon$ and the eigenfunction $u(r)$

e.g. *hydrogen atom 1s orbital*
```julia
r = loggrid(-12,4,500)
v = -1./r
ev,ef = radschrod1(r,v,1,0,reverse=true)
```
"""
function radschrod1(r::UniformGrid{Float64}, v::AbstractVector{Float64},
        n::Int, l::Int; reverse=false, Emin=-1., Emax=0., Elim=1e4)
    r_ = reverse ? Base.reverse(r) : r
    v_ = reverse ? Base.reverse(v) : v

    evs, efs, q = eig(E -> aE(E,l,r_,v_), r.dx, n-l-1, n-l-1,
        Emin=Emin, Emax=Emax, Elim=Elim)

    efs = reverse ? Base.reverse(efs[:,1]) : efs[:,1]  # recover wavefunction
    efs ./= sqrt(integrate(r, efs.^2))  # normalization

    evs[], efs
end

function radschrod1(r::LogarithmicGrid{Float64}, v::AbstractVector{Float64},
        n::Int, l::Int; reverse=false, Emin=-1., Emax=0., Elim=1e4)
    r_ = reverse ? Base.reverse(r) : r
    v_ = reverse ? Base.reverse(v) : v

    evs, efs, q = eig(E -> aElog(E,l,r_,v_), r.dx, n-l-1, n-l-1, Emin=Emin, Emax=Emax, Elim=Elim)

    efs = sqrt(r) .* (reverse ? Base.reverse(efs[:,1]) : efs[:,1])  # recover wavefunction
    efs ./= sqrt(integrate(r, efs.^2))  # normalization

    evs[], efs
end

doc"""
`radschrod(r,v,nmax,lmax;reverse=false,Emin=-1.0,Emax=0.0,Elim=1e4)`

solves the radial Schrödinger equation

$\left\{ -{1\over 2}{d^2\over dr^2} + {l(l+1)\over 2r^2} + v(r) \right\} u(r) = \epsilon u(r)$

where $\varphi(r,\Omega) = {u(r)\over r} Y_l^m(\Omega)$.

for all eigenvalues/eigenfunctions up to `nmax,lmax` (`nmax > 0`, `lmax >= 0` and `lmax < nmax`).

the eigenfunctions are normalized such that $\int dr\, u(r)^2 = 1$

`r` - logarithmic grid (see `loggrid`)

`v` - the radial potential $v(r)$ evaluated on the grid `r`

`reverse` - if true, integrate in reverse.

`Emin,Emax` - initial `E` interval to search. this will be expanded automatically
 (up to `±Elim`) until the range includes all requested eigenvalues.

returns `evs,efs,q`, the requested eigenvalues, eigenfunctions and a list of
the corresponding quantum numbers `n,l`

e.g. *hydrogenic atom, Z=1*
```julia
r = loggrid(-12,5,1000)
v = -1./r
nmax,lmax = 3,2
evs,efs,q = radschrod(r,v,nmax,lmax,reverse=true)
```
"""
function radschrod(r::LogarithmicGrid{Float64}, v::AbstractVector{Float64},
        nmax::Int, lmax::Int; reverse=false, Emin=-1., Emax=0., Elim=1e4)
    r_ = reverse ? Base.reverse(r) : r
    v_ = reverse ? Base.reverse(v) : v

    evs = Float64[]  # eigenvalues
    efs = Vector{Float64}[]  # eigenfunctions
    q = Array{Int}[]  # (n,l) quantum numbers

    for l in 0:lmax
        q_ = eig!(evs, efs, E -> aElog(E,l,r_,v_), r.dx, 0, nmax-l-1,
            Emin=Emin, Emax=Emax, Elim=Elim)
        append!(q, [[k+l+1 l] for k in q_])
    end

    ind = sortperm(evs)
    evs = evs[ind]
    efs = efs[ind]
    q = vcat(q[ind]...)

    for i in 1:length(efs)
        l = q[i,2]
        efs[i] = sqrt(r) .* (reverse ? Base.reverse(efs[i]) : efs[i])  # recover wavefunction
        efs[i] ./= sqrt(integrate(r, efs[i].^2))  # normalization
    end

    efs = hcat(efs...)

    evs, efs, q
end

doc"""
`analytic_continuation!(r,f,g)`

same as `analytic_continuation` but mutates `f`.
"""
function analytic_continuation!(r::AbstractVector, f::AbstractVector, g::Function)
    q = f./g(r)
    i = findfirst(diff(abs(diff(q))) .> 0.0, true)
    c = q[i]
    f[1:i] = c*g(r[1:i])
    return
end

doc"""
`analytic_continuation(r,f,g)`

smoothly replaces $f(r)$ with the analytic form given by $g(r)$ as $r\to 0$

the transition point is determined as the point where
the derivatives of `f` and `g` match.

analytic continuation as $r\to\infty$ can be acheived by reversing
`r` and `f`.

e.g.
```julia
r = loggrid(-12,4.,1000)
n,l = 3,1
ev,ef = radschrod1(r,-1./r,n,l,reverse=true)
efc = analytic_continuation(r,ef,r -> r.^(l+1))
```
"""
function analytic_continuation(r,f,g)
    fc = copy(f)
    analytic_continuation!(r,fc,g)
    return fc
end

doc"""
`analytic_continuation_coulomb_nuclei!(r,f,l)`

analytically continues multiple functions $f(r)$  (columns of `f`) with
the analytic form $r^{l+1}$ as $r\to 0$.

`f` - a matrix whose columns are functions of `r`

`l` - corresponding `l`-quantum numbers
"""
function analytic_continuation_coulomb_nuclei!(r::AbstractVector{Float64},
        f::AbstractMatrix{Float64}, l::AbstractVector{Int})
    for i in 1:size(f,2)
        analytic_continuation!(r, slice(f,:,i), r -> r.^(l[i]+1))
    end
end

function analytic_continuation_coulomb_nuclei!(r::AbstractVector{Float64},
        f::AbstractVector{Float64}, l::Int)
    analytic_continuation!(r, f, r -> r.^(l+1))
end

# import GSL.sf_laguerre_n

# function R_ha(Z::Real,n::Integer,l::Integer,r::Real)
#     x = 2Z*r/n
#     N = factorial(2l+1+n-l-1)*sqrt((2Z/n)^3*factorial(n-l-1)/(2n)/factorial(n+l)^3)
#     N*x^l*exp(-x/2)*sf_laguerre_n(n-l-1,2l+1,x)
# end

# function R_ha(Z::Real,n::Integer,l::Integer,r::AbstractArray)
#     reshape([R_ha(Z,n,l,r_) for r_ in r], size(r))
# end

# u_ha(Z::Real,n::Integer,l::Integer,r::AbstractArray) = r.*R_ha(Z,n,l,r)

max_stars_repo_name: johncsnyder/Element


text:
"Entry point of command line script"
function main()
    parsed_args = parse_commandline()

    seed = parsed_args["seed"]
    if  seed != nothing
        Random.seed!(seed)
    end

    if parsed_args["log_size"] != nothing
        size = 10^parsed_args["log_size"]
    else
        size = parsed_args["size"]
    end

    tree_type = TREE_DICT[parsed_args["tree_type"]]
    tree_ins = missing
    if tree_type == DAryTree
        tree_ins = DAryTree(size, parsed_args["d"])
    elseif tree_type == FullDAryTree
        tree_ins = FullDAryTree(parsed_args["tall"], parsed_args["d"])
    else
        tree_ins = tree_type(size)
    end

    sim = missing
    if parsed_args["%COMMAND%"] == "kcut"
        sim = KcutSimulator(tree_ins, parsed_args["kcut"]["k"])
    elseif parsed_args["%COMMAND%"] == "log-prod"
        sim = LogProductSimulator(tree_ins, parsed_args["log-prod"]["pow"])
    elseif parsed_args["%COMMAND%"] == "size-power"
        sim = SizePowerSimulator(tree_ins, parsed_args["size-power"]["pow"])
    elseif parsed_args["%COMMAND%"] == "height"
        sim = HeightSimulator(tree_ins)
    elseif parsed_args["%COMMAND%"] == "total-path"
        sim = TotalPathSimulator(tree_ins)
    elseif parsed_args["%COMMAND%"] == "leaf"
        sim = LeafSimulator(tree_ins)
    end

    if !ismissing(sim)
        print_simulation(sim, parsed_args["num"])
    end
end


text:
module Utils

export Interval

type Interval{T<:Real}
    lower::T
    upper::T
end

Base.in{T<:Real}(x::T, interval::Interval{T}) = interval.lower <= x <= interval.upper
==(a::Interval, b::Interval) = a.lower == b.lower && a.upper == b.upper
⊆(a::Interval, b::Interval) = b.lower <= a.lower <= a.upper <= b.upper
⊈(a::Interval, b::Interval) = !(a ⊆ b)
⊊(a::Interval, b::Interval) = a ⊆ b && (b.lower < a.lower || a.upper < b.upper)

end  # module Utils


text:
# Taken from https://gist.github.com/JeffBezanson/6b7f1785bb7f2509cbd5d4ff1380556d
# By Jeff Bezanson

import Base: endof, sizeof, pointer, next
using Base: UTF_ERR_INVALID_INDEX, is_valid_continuation, utf8_trailing, utf8_offset

type Str <: AbstractString
    len::Int
    function Str(p::Union{Ptr{Int8},Ptr{UInt8}}, n::Integer)
        s = ccall(:jl_gc_allocobj, Any, (Csize_t,), n+sizeof(Int))
        ps = pointer_from_objref(s)
        unsafe_store!(convert(Ptr{Ptr{Void}}, ps), pointer_from_objref(Str), 0)
        unsafe_store!(convert(Ptr{Int}, ps), n)
        ccall(:memcpy, Void, (Ptr{Void}, Ptr{Void}, Csize_t),
              ps + sizeof(Int), p, n)
        return s
    end
end

Str(d::Union{Vector{Int8},Vector{UInt8}}) = Str(pointer(d), length(d))
Str(s::String, offset=0, endof=length(s.data)) = Str(pointer(s.data)+offset, endof)
Str(s::SubString) = Str(s.string, s.offset, s.endof) # TODO: handle substring of substrings

endof(s::Str) = s.len
sizeof(s::Str) = s.len

pointer(s::Str) = convert(Ptr{UInt8}, pointer_from_objref(s)+sizeof(Int))

@noinline function slow_utf8_next(p::Ptr{UInt8}, b::UInt8, i::Int, l::Int)
    if is_valid_continuation(b)
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, unsafe_load(p,i)))
    end
    trailing = utf8_trailing[b + 1]
    if l < i + trailing
        return '\ufffd', i+1
    end
    c::UInt32 = 0
    for j = 1:(trailing + 1)
        c <<= 6
        c += unsafe_load(p,i)
        i += 1
    end
    c -= utf8_offset[trailing + 1]
    return Char(c), i
end

@inline function next(s::Str, i::Int)
    if i < 1 || i > s.len
        throw(BoundsError(s,i))
    end
    p = pointer(s)
    b = unsafe_load(p, i)
    if b < 0x80
        return Char(b), i+1
    end
    return slow_utf8_next(p, b, i, s.len)
end


text:
function get_walker(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.5,
    contact_feet=true,
    contact_body=true,
    limits=true,
    spring=0.0,
    damper=0.0,
    parse_damper=true,
    joint_limits=[[  0.0,   0.0, -45.0,   0.0,   0.0, -45.0] * π / 180.0,
                  [150.0, 150.0,  45.0, 150.0, 150.0,  45.0] * π / 180.0],
    T=Float64)

    path = joinpath(@__DIR__, "../deps/walker.urdf")
    mech = Mechanism(path; floating=false, T,
        gravity,
        timestep,
        parse_damper)

    # Adding springs and dampers
    set_springs!(mech.joints, spring)
    set_dampers!(mech.joints, damper)

    # joint limits
    joints = deepcopy(mech.joints)

    if limits
        thigh = get_joint(mech, :thigh)
        joints[thigh.id] = add_limits(mech, thigh,
            rot_limits=[SVector{1}(joint_limits[1][1]), SVector{1}(joint_limits[2][1])])

        leg = get_joint(mech, :leg)
        joints[leg.id] = add_limits(mech, leg,
            rot_limits=[SVector{1}(joint_limits[1][2]), SVector{1}(joint_limits[2][2])])

        foot = get_joint(mech, :foot)
        joints[foot.id] = add_limits(mech, foot,
            rot_limits=[SVector{1}(joint_limits[1][3]), SVector{1}(joint_limits[2][3])])

        thigh_left = get_joint(mech, :thigh_left)
        joints[thigh_left.id] = add_limits(mech, thigh_left,
            rot_limits=[SVector{1}(joint_limits[1][4]), SVector{1}(joint_limits[2][4])])

        leg_left = get_joint(mech, :leg_left)
        joints[leg_left.id] = add_limits(mech, leg_left,
            rot_limits=[SVector{1}(joint_limits[1][5]), SVector{1}(joint_limits[2][5])])

        foot_left = get_joint(mech, :foot_left)
        joints[foot_left.id] = add_limits(mech, foot_left,
            rot_limits=[SVector{1}(joint_limits[1][6]), SVector{1}(joint_limits[2][6])])

        mech = Mechanism(Origin{T}(), [mech.bodies...], [joints...];
            gravity,
            timestep)
    end

    if contact_feet
        origin = Origin{T}()
        bodies = mech.bodies
        joints = mech.joints

        normal = [0.0; 0.0; 1.0]
        names = contact_body ? getfield.(mech.bodies, :name) : [:ffoot, :foot]
        models = []
        for name in names
            body = get_body(mech, name)
            if name in [:foot, :foot_left] # need special case for foot
                # torso
                pf = [0.0, 0.0, 0.5 * body.shape.rh[2]]
                pb = [0.0, 0.0, -0.5 * body.shape.rh[2]]
                o = body.shape.rh[1]
                push!(models, contact_constraint(body, normal;
                    friction_coefficient,
                    contact_origin=pf,
                    contact_radius=o))
                push!(models, contact_constraint(body, normal;
                    friction_coefficient,
                    contact_origin=pb,
                    contact_radius=o))
            else
                p = [0.0; 0.0; 0.5 * body.shape.rh[2]]
                o = body.shape.rh[1]
                push!(models, contact_constraint(body, normal;
                    friction_coefficient,
                    contact_origin=p,
                    contact_radius=o))
            end
        end

        set_minimal_coordinates!(mech, get_joint(mech, :floating_joint), [1.25, 0.0, 0.0])

        mech = Mechanism(origin, bodies, joints, [models...];
            gravity,
            timestep)
    end
    return mech
end

function initialize_walker!(mechanism::Mechanism;
    body_position=[0.0, 0.0],
    body_orientation=0.0) where T

    set_minimal_coordinates!(mechanism,
                 get_joint(mechanism, :floating_joint),
                 [body_position[2] + 1.25 , -body_position[1], -body_orientation])
    for joint in mechanism.joints
        (joint.name != :floating_joint) && set_minimal_coordinates!(mechanism, joint, zeros(input_dimension(joint)))
    end
    zero_velocity!(mechanism)
end


text:
module LIGHTPATH
using ..AGFFileReader: Glass, GlassID, AGF
export D_ZK3M, C0550, BD2, D_ZLAF52LA_M, ECO550, NICHIA_MELT1, BD1, ECO550_E, ECO550M, PBH71, D_ZLAF52LAM

const D_ZK3M = Glass(GlassID(AGF, 1639), 1, 2.481845, -0.01213254, 0.01179439, 0.000804573, -9.005814e-5, 4.898641e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 1.7034e-6, 1.3949e-8, -2.2949e-11, 5.5488e-7, 2.3501e-10, 0.12523, 20.0, -0.0066, 1.0, -1.0, 7.6, 2.5, 0, 4.2, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.06, 10.0), (0.31, 0.27, 10.0), (0.32, 0.53, 10.0), (0.334, 0.786, 10.0), (0.35, 0.922, 10.0), (0.365, 0.967, 10.0), (0.37, 0.975, 10.0), (0.38, 0.984, 10.0), (0.39, 0.989, 10.0), (0.4, 0.992, 10.0), (0.404, 0.992, 10.0), (0.42, 0.993, 10.0), (0.436, 0.993, 10.0), (0.46, 0.995, 10.0), (0.5, 0.998, 10.0), (0.546, 0.999, 10.0), (0.58, 0.998, 10.0), (0.62, 0.998, 10.0), (0.66, 0.998, 10.0), (0.7, 0.998, 10.0), (1.06, 0.997, 10.0), (1.529, 0.993, 10.0), (1.97, 0.98, 10.0), (2.325, 0.867, 10.0)], 1.58642, -1.0, -1.0, 0, 60.713872, 0, 2.83, 0)
const C0550 = Glass(GlassID(AGF, 1640), 1, 2.5275, -0.00824759, 0.0161384, 0.000710762, -7.35332e-5, 6.55454e-6, 0.0, 0.0, 0.0, 0.0, 0.32, 1.7, -2.2341e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.002901484, -1.0, -1.0, 15.0, -1.0, 0, -1.0, [(0.32, 0.0, 10.0), (0.33, 0.0, 10.0), (0.34, 0.164, 10.0), (0.35, 0.503, 10.0), (0.36, 0.748, 10.0), (0.37, 0.878, 10.0), (0.38, 0.939, 10.0), (0.39, 0.966, 10.0), (0.4, 0.975, 10.0), (0.42, 0.975, 10.0), (0.44, 0.974, 10.0), (0.45, 0.976, 10.0), (0.47, 0.98, 10.0), (0.5, 0.987, 10.0), (0.55, 0.989, 10.0), (0.65, 0.991, 10.0), (1.0, 0.996, 10.0), (1.4, 0.996, 10.0), (1.6, 0.965, 10.0)], 1.605004, -1.0, -1.0, 0, 50.419332, 0, 3.8, 0)
const BD2 = Glass(GlassID(AGF, 1641), 3, 2.61445159, 0.14913534, -0.287509165, -9.57306751e-5, -5.10891622e-7, 9.89417441e-10, 0.0, 0.0, 0.0, 0.0, 2.0, 15.0, 4.72e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 16.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.0, -1.0, -1.0, 0, 0.0, 0, 4.67, 0)
const D_ZLAF52LA_M = Glass(GlassID(AGF, 1642), 1, 3.1717081, -0.013708448, 0.029629026, 0.0011150095, -5.5576606e-5, 5.6941122e-6, 0.0, 0.0, 0.0, 0.0, 0.36, 1.55, 5.35e-6, 1.8526e-8, -7.872e-11, 7.1515e-7, 1.2648e-10, -0.23138, 20.0, -0.0078, -1.0, -1.0, 6.9, 1.0, 0, 3.0, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.0, 10.0), (0.31, 0.0, 10.0), (0.32, 0.0, 10.0), (0.33, 0.031, 10.0), (0.34, 0.076, 10.0), (0.35, 0.326, 10.0), (0.36, 0.584, 10.0), (0.37, 0.754, 10.0), (0.38, 0.854, 10.0), (0.39, 0.91, 10.0), (0.4, 0.943, 10.0), (0.42, 0.972, 10.0), (0.44, 0.981, 10.0), (0.46, 0.985, 10.0), (0.48, 0.99, 10.0), (0.5, 0.995, 10.0), (0.55, 0.996, 10.0), (0.6, 0.996, 10.0), (0.65, 0.996, 10.0), (0.7, 0.998, 10.0), (0.8, 0.998, 10.0), (1.06, 0.998, 10.0), (1.5, 0.998, 10.0), (2.0, 0.973, 10.0)], 1.805881, -1.0, -1.0, 0, 40.888056, 0, 4.56, 0)
const ECO550 = Glass(GlassID(AGF, 1643), 1, 2.52780621, -0.00782886154, 0.0171870201, 0.000300693162, 8.66557318e-6, 2.58488712e-7, 0.0, 0.0, 0.0, 0.0, 0.4, 1.55, 4.9178e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0027, -1.0, -1.0, 11.62, -1.0, 0, -1.0, [(0.29, 0.57, 10.0), (0.3, 0.82, 10.0), (0.31, 0.928, 10.0), (0.32, 0.97, 10.0), (0.33, 0.988, 10.0), (0.34, 0.996, 10.0), (0.35, 0.997, 10.0), (0.36, 0.981, 10.0), (0.37, 0.984, 10.0), (0.38, 0.983, 10.0), (0.39, 0.986, 10.0), (0.4, 0.988, 10.0), (0.42, 0.989, 10.0), (0.44, 0.99, 10.0), (0.46, 0.992, 10.0), (0.48, 0.993, 10.0), (0.5, 0.996, 10.0), (0.55, 0.998, 10.0), (0.6, 0.999, 10.0), (0.65, 0.999, 10.0), (0.7, 0.999, 10.0)], 1.605503, -1.0, -1.0, 0, 50.277746, 0, 3.31, 0)
const NICHIA_MELT1 = Glass(GlassID(AGF, 1644), 5, 1.47444837, 0.0103147698, 0.00026742387, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.334, 2.325, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.493724, -1.0, -1.0, 0, 64.358478, 0, 1.0, 0)
const BD1 = Glass(GlassID(AGF, 1645), 3, 2.51253618, 0.0686942636, 0.0211563734, -0.000155048679, -4.35296121e-8, 1.99858193e-10, 0.0, 0.0, 0.0, 0.0, 1.4, 15.0, 7.05e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 13.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.0, -1.0, -1.0, 0, 0.0, 0, 4.4, 0)
const ECO550_E = Glass(GlassID(AGF, 1646), 1, 2.50983687, 1.51635542e-6, 0.0148303868, 0.000981490801, 3.23688017e-6, 7.98861178e-8, 0.0, 0.0, 0.0, 0.0, 0.25, 2.05, 4.8817e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.04895, -1.0, -1.0, 11.62, -1.0, 0, -1.0, [], 1.600348, -1.0, -1.0, 0, 47.142136, 0, 3.32, 0)
const ECO550M = Glass(GlassID(AGF, 1647), 1, 2.52035773, -0.00825366386, 0.0165830522, 0.000476837067, -1.50343126e-5, 1.4015069e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 1.9701, 4.8817e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, -0.000308066, -1.0, 100.0, 11.62, 2.0, 0, 5.0, [(0.29, 0.57, 10.0), (0.3, 0.82, 10.0), (0.31, 0.928, 10.0), (0.32, 0.97, 10.0), (0.33, 0.988, 10.0), (0.34, 0.996, 10.0), (0.35, 0.997, 10.0), (0.36, 0.981, 10.0), (0.37, 0.984, 10.0), (0.38, 0.983, 10.0), (0.39, 0.986, 10.0), (0.4, 0.988, 10.0), (0.42, 0.989, 10.0), (0.44, 0.99, 10.0), (0.46, 0.992, 10.0), (0.48, 0.993, 10.0), (0.5, 0.996, 10.0), (0.55, 0.998, 10.0), (0.6, 0.999, 10.0), (0.65, 0.999, 10.0), (0.7, 0.999, 10.0)], 1.602897, -1.0, -1.0, 0, 50.021894, 0, 3.32, 9)
const PBH71 = Glass(GlassID(AGF, 1648), 2, 1.97039844, 0.0155780448, 0.520330498, 0.0632642123, 1.30696186, 142.045097, 0.0, 0.0, 0.0, 0.0, 0.3, 2.4, 7.8574e-6, 1.8258e-8, 1.11e-11, 1.9714e-6, 1.6131e-9, 0.30894, 25.0, 0.0204, 4.0, 2.0, 8.9, -1.0, 2, 53.4, [(0.3, 0.0, 10.0), (0.31, 0.0, 10.0), (0.32, 0.0, 10.0), (0.33, 0.0, 10.0), (0.34, 0.0, 10.0), (0.35, 0.0, 10.0), (0.36, 0.0, 10.0), (0.37, 0.02, 10.0), (0.38, 0.23, 10.0), (0.39, 0.56, 10.0), (0.4, 0.77, 10.0), (0.42, 0.925, 10.0), (0.44, 0.962, 10.0), (0.46, 0.976, 10.0), (0.48, 0.985, 10.0), (0.5, 0.991, 10.0), (0.55, 0.998, 10.0), (0.6, 0.998, 10.0), (0.65, 0.999, 10.0), (0.7, 0.999, 10.0), (0.8, 0.999, 10.0), (0.9, 0.998, 10.0), (1.0, 0.998, 10.0), (1.2, 0.997, 10.0), (1.4, 0.995, 10.0), (1.6, 0.993, 10.0), (1.8, 0.985, 10.0), (2.0, 0.975, 10.0), (2.2, 0.953, 10.0), (2.4, 0.918, 10.0)], 1.922864, -1.0, -1.0, 0, 21.288699, 0, 6.05, 5)
const D_ZLAF52LAM = Glass(GlassID(AGF, 1649), 1, 3.1607213, -0.0084677629, 0.037044808, -0.0010934271, 0.000243631811, -8.9247355e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 5.35e-6, 1.8526e-8, -7.872e-11, 7.1515e-7, 1.2648e-10, -0.23138, 20.0, -0.00587, -1.0, -1.0, 6.9, 1.0, 0, 3.0, [(2.4, 0.75, 10.0), (2.2, 0.924, 10.0), (2.0, 0.969, 10.0), (1.8, 0.992, 10.0), (1.6, 0.998, 10.0), (1.4, 0.996, 10.0), (1.2, 0.998, 10.0), (1.06, 0.998, 10.0), (1.0, 0.998, 10.0), (0.95, 0.998, 10.0), (0.9, 0.998, 10.0), (0.85, 0.998, 10.0), (0.8, 0.998, 10.0), (0.7, 0.997, 10.0), (0.65, 0.996, 10.0), (0.6, 0.995, 10.0), (0.55, 0.995, 10.0), (0.5, 0.99, 10.0), (0.48, 0.985, 10.0), (0.46, 0.98, 10.0), (0.44, 0.972, 10.0), (0.42, 0.961, 10.0), (0.4, 0.937, 10.0), (0.39, 0.913, 10.0), (0.38, 0.87, 10.0), (0.37, 0.79, 10.0), (0.36, 0.64, 10.0), (0.35, 0.38, 10.0), (0.34, 0.08, 10.0)], 1.805885, -1.0, -1.0, 0, 40.792034, 0, 4.56, 1)
end # module

max_stars_repo_name: rambunctiousapple/GlassCat

max_issues_repo_name: rambunctiousapple/GlassCat

max_forks_repo_name: rambunctiousapple/GlassCat