text:
using CodecZlib
using NIfTI
using Test
using TranscodingStreams

function extractto(gzname, out)
    open(out, "w") do io_out
        open(gzname, "r") do io_in
            gz = GzipDecompressorStream(io_in)
            write(io_out, read(gz))
        end
    end
end

const TEMP_DIR_NAME = mktempdir()

# empty file test
const EMPTY_NII = joinpath(dirname(@__FILE__), "data/empty.nii.gz")
@test_throws EOFError niread(EMPTY_NII)

# single file storage
const GZIPPED_NII = joinpath(dirname(@__FILE__), "data/example4d.nii.gz")
const NII = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
extractto(GZIPPED_NII, NII)

# dual file storage
const GZIPPED_HDR = joinpath(dirname(@__FILE__), "data/example4d.hdr.gz")
hdr_stem = joinpath(TEMP_DIR_NAME, tempname())
const HDR = "$hdr_stem.hdr"
const IMG = "$hdr_stem.img"
extractto(GZIPPED_HDR, HDR)
extractto(joinpath(dirname(@__FILE__), "data/example4d.img.gz"), IMG)

function image_tests(fname, mmap)
    img = niread(fname, mmap=mmap)

    # Header
    @test time_step(img.header) == 2000000 # Actually an error in the file AFAIK
    @test voxel_size(img.header) ≈ Float32[2.0, 2.0, 2.2]
    @test size(img) == (128, 96, 24, 2)

    # Content
    @test img.raw[65, 49, 13, :][:] == [265, 266]
    @test vox(img, 64, 48, 12, :)[:] == [265, 266]
    @test vox(img, 69, 56, 13, :)[:] == [502, 521]

    @assert maximum(img) == maximum(img.raw)

    @test getaffine(img) ≈ [
        -2.0                    6.714715653593746e-19  9.081024511081715e-18  117.8551025390625
        6.714715653593746e-19  1.9737114906311035    -0.35552823543548584   -35.72294235229492
        8.25548088896093e-18   0.3232076168060303     2.171081781387329     -7.248798370361328
        0.0                    0.0                    0.0                   1.0
    ]

    @test NIfTI.get_qform(img) ≈ [
     -2.0          7.75482f-26  -6.93824f-27  117.855
      7.75482f-26  1.97371      -0.355528     -35.7229
      6.30749f-27  0.323208      2.17108       -7.2488
      0.0          0.0           0.0            1.0
    ]
    @test NIfTI.orientation(img) == (:right, :posterior, :inferior)
end

image_tests(NII, false)
image_tests(NII, true)
image_tests(HDR, false)
image_tests(HDR, true)
image_tests(GZIPPED_NII, false)
image_tests(GZIPPED_HDR, false)

@testset "Header Field Accessors" begin
    img = niread(GZIPPED_NII)
    @test NIfTI.freqdim(img) == 1
    @test NIfTI.phasedim(img) == 2
    @test NIfTI.slicedim(img) == 3
    @test NIfTI.slice_start(img) == 1
    @test NIfTI.slice_end(img) == 24
    @test NIfTI.slice_duration(img) == 0
end

@test_throws ErrorException niread(GZIPPED_NII; mmap=true)
@test_throws ErrorException niread(GZIPPED_HDR; mmap=true)

# Test writing
const TEMP_FILE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
vol = NIVolume()
niwrite(TEMP_FILE, vol)
niread(TEMP_FILE)

const TEMP_GZIPPED_FILE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii.gz")
niwrite(TEMP_GZIPPED_FILE, vol)
niread(TEMP_GZIPPED_FILE)

# Write and read DT_BINARY
const BOOL_WRITE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
const BIT_WRITE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
mask = rand(Bool, 3, 5, 7) # Array{Bool}
mask_bitarray = BitArray(mask) # BitArray
niwrite(BOOL_WRITE, NIVolume(mask))
niwrite(BIT_WRITE, NIVolume(mask_bitarray))
@test niread(BOOL_WRITE).raw == mask
@test niread(BIT_WRITE).raw == mask_bitarray

# Open mmaped file for reading and writing
const WRITE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
const VERIFY_WRITE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
cp(NII, WRITE)
img = niread(WRITE; mmap=true, mode="r+")
img.raw[1,1,1,1] = 5
img.raw[:,2,1,1] = ones(size(img)[1])
cp(WRITE, VERIFY_WRITE)
@test niread(VERIFY_WRITE)[1,1,1,1] == 5
@test niread(VERIFY_WRITE)[:,2,1,1] == ones(size(img)[1])
# Site is currently down TODO: reintroduce this test when site is up
# Big endian
# const BE = "$(tempname()).nii"
# download("https://nifti.nimh.nih.gov/nifti-1/data/avg152T1_LR_nifti.nii.gz", BE)
img = niread(joinpath(dirname(@__FILE__), "data/avg152T1_LR_nifti.nii.gz"))
@test size(img) == (91,109,91)

GC.gc() # closes mmapped files

@test NIfTI._dir2ori(-1.0, 0.0, 0.0,
                      0.0, 1.0, 0.0,
                      0.0, 0.0, 1.0) == (:right, :posterior, :inferior)

@test NIfTI._dir2ori(1.0,  0.0, 0.0,
                     0.0, -1.0, 0.0,
                     0.0,  0.0, 1.0) == (:left, :anterior, :inferior)

@test NIfTI._dir2ori(1.0,  0.0, 0.0,
                     0.0, -1.0, 0.0,
                     0.0,  0.0, -1.0) == (:left, :anterior, :superior)


text:
############################################################################
#                        TEST FUNCIONES DE ACTIVACIÓN
# Comprobamos que satisfacen alguno de sus requisitos teóricos
# Así como que al evaluarlas no de ningún error
############################################################################
using Test

######################### TEST  #########################
@testset "Activations functions" begin

    id(x)=x
    # Propiedas asintóticas funciones clásicas
    @test @ThresholdFunction(id,0)(-1) ≈ -1
    @test @ThresholdFunction(id,0)(1) ≈ 1
    @test CosineSquasher(-100) ≈ 0
    @test CosineSquasher(10) ≈ 1
    @test @IndicatorFunction(0)(-0.1) == 0
    @test @IndicatorFunction(0)(0.1) == 1
    @test RampFunction(-1)== 0
    @test RampFunction(2) == 1
    @test ReLU(-1)== 0
    @test ReLU(2) >0
    @test Sigmoid(9999999) ≈ 1
    @test Sigmoid(-9999999) ≈ 0
    @test HardTanh(9999999) ≈ 1
    @test HardTanh(-9999999) ≈ -1
    @test @LReLU(0.01)(10) == 10
    @test @LReLU(-0.01)(-1) == 0.01

    # Medibles en la práctica equivale a que tome reales
    valores_compacto = [ (rand()-0.5)*1000 for _ in 1:10]

    funciones_activacion = [
        CosineSquasher,
        RampFunction,
        ReLU,
        Sigmoid,
        HardTanh
    ]
    # Comprobamos funciones de activación que no son macros
    for f in funciones_activacion
        for x in valores_compacto
            @test typeof(f(x)) <: Real
        end
    end
    # Las macro me sale error al meternas en la lista :(
    for x in valores_compacto
        @test typeof(@ThresholdFunction(id,0)(x)) <: Real
    end
    for x in valores_compacto
        @test typeof(@IndicatorFunction(0)(x)) <: Real
    end
    for x in valores_compacto
        @test typeof(@LReLU(0.01)(x)) <: Real
    end
end

ext: jl


text:
data = readline()

function increment(c::Char)
    if c == 'z'
        return 'a'
    end
    return c + 1
end

function increment(s::String)
    if s[end] == 'z'
        return increment(s[1:end-1]) * increment(s[end])
    end
    return s[1:end-1] * increment(s[end])
end

function check(s)
    rule_1 = false
    for i in 1:length(s)-2
        if s[i]+2 == s[i+1]+1 == s[i+2]
            rule_1 = true
        end
    end
    rule_2 = !occursin.("iol", Ref(s))
    rule_3 = length(eachmatch(r"(.)\1", s, overlap=false) |> collect) >= 2
    rule_1 && rule_2 && rule_3
end

function solve_1()
    s = data |> increment
    while !check(s)
        s = increment(s)
    end
    s
end

solve_1() |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
using Test
using ClimateMachine.Mesh.Topologies
using Combinatorics, MPI

@testset "cubedshellwarp tests" begin
    import ClimateMachine.Mesh.Topologies: cubedshellwarp

    @testset "check radius" begin
        @test hypot(cubedshellwarp(3.0, -2.2, 1.3)...) ≈ 3.0 rtol = eps()
        @test hypot(cubedshellwarp(-3.0, -2.2, 1.3)...) ≈ 3.0 rtol = eps()
        @test hypot(cubedshellwarp(1.1, -2.2, 3.0)...) ≈ 3.0 rtol = eps()
        @test hypot(cubedshellwarp(1.1, -2.2, -3.0)...) ≈ 3.0 rtol = eps()
        @test hypot(cubedshellwarp(1.1, 3.0, 0.0)...) ≈ 3.0 rtol = eps()
        @test hypot(cubedshellwarp(1.1, -3.0, 0.0)...) ≈ 3.0 rtol = eps()
    end

    @testset "check sign" begin
        @test sign.(cubedshellwarp(3.0, -2.2, 1.3)) == sign.((3.0, -2.2, 1.3))
        @test sign.(cubedshellwarp(-3.0, -2.2, 1.3)) == sign.((-3.0, -2.2, 1.3))
        @test sign.(cubedshellwarp(1.1, -2.2, 3.0)) == sign.((1.1, -2.2, 3.0))
        @test sign.(cubedshellwarp(1.1, -2.2, -3.0)) == sign.((1.1, -2.2, -3.0))
        @test sign.(cubedshellwarp(1.1, 3.0, 0.0)) == sign.((1.1, 3.0, 0.0))
        @test sign.(cubedshellwarp(1.1, -3.0, 0.0)) == sign.((1.1, -3.0, 0.0))
    end

    @testset "check continuity" begin
        for (u, v) in zip(
            permutations([3.0, 2.999999999, 1.3]),
            permutations([2.999999999, 3.0, 1.3]),
        )
            @test all(cubedshellwarp(u...) .≈ cubedshellwarp(v...))
        end
        for (u, v) in zip(
            permutations([3.0, -2.999999999, 1.3]),
            permutations([2.999999999, -3.0, 1.3]),
        )
            @test all(cubedshellwarp(u...) .≈ cubedshellwarp(v...))
        end
        for (u, v) in zip(
            permutations([-3.0, 2.999999999, 1.3]),
            permutations([-2.999999999, 3.0, 1.3]),
        )
            @test all(cubedshellwarp(u...) .≈ cubedshellwarp(v...))
        end
        for (u, v) in zip(
            permutations([-3.0, -2.999999999, 1.3]),
            permutations([-2.999999999, -3.0, 1.3]),
        )
            @test all(cubedshellwarp(u...) .≈ cubedshellwarp(v...))
        end
    end
end

@testset "cubedshellunwarp" begin
    import ClimateMachine.Mesh.Topologies: cubedshellwarp, cubedshellunwarp

    for u in permutations([3.0, 2.999999999, 1.3])
        @test all(cubedshellunwarp(cubedshellwarp(u...)...) .≈ u)
    end
    for u in permutations([3.0, -2.999999999, 1.3])
        @test all(cubedshellunwarp(cubedshellwarp(u...)...) .≈ u)
    end
    for u in permutations([-3.0, 2.999999999, 1.3])
        @test all(cubedshellunwarp(cubedshellwarp(u...)...) .≈ u)
    end
    for u in permutations([-3.0, -2.999999999, 1.3])
        @test all(cubedshellunwarp(cubedshellwarp(u...)...) .≈ u)
    end
end

@testset "grid1d" begin
    g = grid1d(0, 10, nelem = 10)
    @test eltype(g) == Float64
    @test length(g) == 11
    @test g[1] == 0
    @test g[end] == 10

    g = grid1d(10f0, 20f0, elemsize = 0.1)
    @test eltype(g) == Float32
    @test length(g) == 101
    @test g[1] == 10
    @test g[end] == 20

    g = grid1d(10f0, 20f0, InteriorStretching(0), elemsize = 0.1)
    @test eltype(g) == Float32
    @test length(g) == 101
    @test g[1] == 10
    @test g[end] == 20

    g = grid1d(10f0, 20f0, SingleExponentialStretching(2.5f0), elemsize = 0.1)
    @test eltype(g) == Float32
    @test length(g) == 101
    @test g[1] == 10
    @test g[end] == 20
end

@testset "BrickTopology tests" begin

    let
        comm = MPI.COMM_SELF

        elemrange = (0:10,)
        periodicity = (true,)

        topology = BrickTopology(comm, elemrange, periodicity = periodicity)

        nelem = length(elemrange[1]) - 1

        for e in 1:nelem
            @test topology.elemtocoord[:, :, e] == [e - 1 e]
        end

        @test topology.elemtoelem ==
              [nelem collect(1:(nelem - 1))'; collect(2:nelem)' 1]
        @test topology.elemtoface == repeat(2:-1:1, outer = (1, nelem))

        @test topology.elemtoordr == ones(Int, size(topology.elemtoordr))
        @test topology.elemtobndy == zeros(Int, size(topology.elemtoordr))

        @test topology.elems == 1:nelem
        @test topology.realelems == 1:nelem
        @test topology.ghostelems == nelem .+ (1:0)

        @test length(topology.sendelems) == 0
        @test length(topology.exteriorelems) == 0
        @test collect(topology.realelems) == topology.interiorelems

        @test topology.nabrtorank == Int[]
        @test topology.nabrtorecv == UnitRange{Int}[]
        @test topology.nabrtosend == UnitRange{Int}[]
    end

    let
        comm = MPI.COMM_SELF
        topology = BrickTopology(comm, (0:4, 5:9), periodicity = (false, true))

        nelem = 16

        @test topology.elemtocoord[:, :, 1] == [0 1 0 1; 5 5 6 6]
        @test topology.elemtocoord[:, :, 2] == [1 2 1 2; 5 5 6 6]
        @test topology.elemtocoord[:, :, 3] == [1 2 1 2; 6 6 7 7]
        @test topology.elemtocoord[:, :, 4] == [0 1 0 1; 6 6 7 7]
        @test topology.elemtocoord[:, :, 5] == [0 1 0 1; 7 7 8 8]
        @test topology.elemtocoord[:, :, 6] == [0 1 0 1; 8 8 9 9]
        @test topology.elemtocoord[:, :, 7] == [1 2 1 2; 8 8 9 9]
        @test topology.elemtocoord[:, :, 8] == [1 2 1 2; 7 7 8 8]
        @test topology.elemtocoord[:, :, 9] == [2 3 2 3; 7 7 8 8]
        @test topology.elemtocoord[:, :, 10] == [2 3 2 3; 8 8 9 9]
        @test topology.elemtocoord[:, :, 11] == [3 4 3 4; 8 8 9 9]
        @test topology.elemtocoord[:, :, 12] == [3 4 3 4; 7 7 8 8]
        @test topology.elemtocoord[:, :, 13] == [3 4 3 4; 6 6 7 7]
        @test topology.elemtocoord[:, :, 14] == [2 3 2 3; 6 6 7 7]
        @test topology.elemtocoord[:, :, 15] == [2 3 2 3; 5 5 6 6]
        @test topology.elemtocoord[:, :, 16] == [3 4 3 4; 5 5 6 6]

        @test topology.elemtoelem == [
            1 1 4 4 5 6 6 5 8 7 10 9 14 3 2 15
            2 15 14 3 8 7 10 9 12 11 11 12 13 13 16 16
            6 7 2 1 4 5 8 3 14 9 12 13 16 15 10 11
            4 3 8 5 6 1 2 7 10 15 16 11 12 9 14 13
        ]

        @test topology.elemtoface == [
            1 2 2 1 1 1 2 2 2 2 2 2 2 2 2 2
            1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 2
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
        ]

        @test topology.elemtoordr == ones(Int, size(topology.elemtoordr))

        @test topology.elems == 1:nelem
        @test topology.realelems == 1:nelem
        @test topology.ghostelems == nelem .+ (1:0)

        @test length(topology.sendelems) == 0
        @test length(topology.exteriorelems) == 0
        @test collect(topology.realelems) == topology.interiorelems

        @test topology.nabrtorank == Int[]
        @test topology.nabrtorecv == UnitRange{Int}[]
        @test topology.nabrtosend == UnitRange{Int}[]
    end

    let
        comm = MPI.COMM_SELF
        for px in (true, false)
            topology = BrickTopology(comm, (0:10,), periodicity = (px,))
            @test Topologies.hasboundary(topology) == !px
        end
        for py in (true, false), px in (true, false)
            topology = BrickTopology(comm, (0:10, 0:3), periodicity = (px, py))
            @test Topologies.hasboundary(topology) == !(px && py)
        end
        for pz in (true, false), py in (true, false), px in (true, false)
            topology = BrickTopology(
                comm,
                (0:10, 0:3, -1:3),
                periodicity = (px, py, pz),
            )
            @test Topologies.hasboundary(topology) == !(px && py && pz)
        end
    end
end

@testset "StackedBrickTopology tests" begin
    let
        comm = MPI.COMM_SELF
        topology = StackedBrickTopology(
            comm,
            (2:5, 4:6),
            periodicity = (false, true),
            boundary = ((1, 2), (3, 4)),
        )

        nelem = 6

        @test topology.elemtocoord[:, :, 1] == [2 3 2 3; 4 4 5 5]
        @test topology.elemtocoord[:, :, 2] == [2 3 2 3; 5 5 6 6]
        @test topology.elemtocoord[:, :, 3] == [3 4 3 4; 4 4 5 5]
        @test topology.elemtocoord[:, :, 4] == [3 4 3 4; 5 5 6 6]
        @test topology.elemtocoord[:, :, 5] == [4 5 4 5; 4 4 5 5]
        @test topology.elemtocoord[:, :, 6] == [4 5 4 5; 5 5 6 6]

        @test topology.elemtoelem == [
            1 2 1 2 3 4
            3 4 5 6 5 6
            2 1 4 3 6 5
            2 1 4 3 6 5
        ]

        @test topology.elemtoface == [
            1 1 2 2 2 2
            1 1 1 1 2 2
            4 4 4 4 4 4
            3 3 3 3 3 3
        ]

        @test topology.elemtoordr == ones(Int, size(topology.elemtoordr))

        @test topology.elemtobndy == [
            1 1 0 0 0 0
            0 0 0 0 2 2
            0 0 0 0 0 0
            0 0 0 0 0 0
        ]

        @test topology.elems == 1:nelem
        @test topology.realelems == 1:nelem
        @test topology.ghostelems == nelem .+ (1:0)

        @test length(topology.sendelems) == 0
        @test length(topology.exteriorelems) == 0
        @test collect(topology.realelems) == topology.interiorelems

        @test topology.nabrtorank == Int[]
        @test topology.nabrtorecv == UnitRange{Int}[]
        @test topology.nabrtosend == UnitRange{Int}[]
    end
    let
        comm = MPI.COMM_SELF
        for py in (true, false), px in (true, false)
            topology =
                StackedBrickTopology(comm, (0:10, 0:3), periodicity = (px, py))
            @test Topologies.hasboundary(topology) == !(px && py)
        end
        for pz in (true, false), py in (true, false), px in (true, false)
            topology = StackedBrickTopology(
                comm,
                (0:10, 0:3, -1:3),
                periodicity = (px, py, pz),
            )
            @test Topologies.hasboundary(topology) == !(px && py && pz)
        end
    end
end

@testset "StackedCubedSphereTopology tests" begin
    topology = StackedCubedSphereTopology(MPI.COMM_SELF, 3, 1.0:3.0)
    @test Topologies.hasboundary(topology)
end

@testset "CubedShellTopology tests" begin
    topology = CubedShellTopology(MPI.COMM_SELF, 3, Float64)
    @test !Topologies.hasboundary(topology)
end


text:
using CxxWrap
using MacroTools

macro debugcb(ex)
    seen = false
    ex = MacroTools.prewalk(ex) do x
        if !seen && MacroTools.@capture(MacroTools.shortdef(x), f_(xs__) = body_)
            seen = true
            def = MacroTools.splitdef(x)
            # debug_ex = :($Base.@info $f)
            # push!(debug_ex.args, def[:args]...)
            debug_ex = :(println($f, $(def[:args]...)))
            pushfirst!(def[:body].args, debug_ex)
            return MacroTools.combinedef(def)
        else
            return x
        end
    end
    return esc(ex)
end

# @debugcb function baz(x,y)
#     println("yo")
# end
println(prettify(@macroexpand @debugcb @cxxdereference function baz(x)
    println("yo")
end))
# @debugcb @cxxdereference function foo(x)
#     function swag(y)
#         println("waoh $y")
#     end
#     swag(10)
# end


text:
push!(LOAD_PATH, string(@__DIR__, "/.."))
using Documenter;
using NonlinearEigenproblems

makedocs(
    clean = true,
    doctest = false,
    sitename = "NEP-PACK",
    pages = Any[
        "Introduction" => "index.md",
        "Manual" => [
        "methods.md",
            "types.md",
        "compute_functions.md",
        "linsolvers.md",
        "innersolvers.md",
        "errmeasure.md",
        "logger.md",
        "transformations.md",
        "deflation.md",
            "gallery.md"],
        "Tutorials" => [
        "Tutorial 1 (ABC)" => "movebc_tutorial.md",
        "Tutorial 2 (Contour)" => "tutorial_contour.md",
        "Tutorial 3 (BEM)" => "bemtutorial.md",
        "Tutorial 4 (Deflation)" => "deflate_tutorial.md",
        "Tutorial 5 (Python 1)" => "tutorial_call_python.md",
        "Tutorial 6 (Python 2)" => "tutorial_python_call.md",
        "Tutorial 7 (MATLAB)" => "tutorial_matlab1.md",
        "Tutorial 8 (FORTRAN)" => "tutorial_fortran1.md",
        "Tutorial 9 (gmsh + nanophotonics)" => "tutorial_nano1.md",
        "Tutorial 10 (New solver)" => "tutorial_newmethod.md",
        "Tutorial 11 (Linear solvers)" => "tutorial_linsolve.md",
            "Tutorial 12 (Orr--Somerfeld)" => "hydrotutorial.md"
            ]
    ]
)

#deploydocs(
#    repo = "github.com/nep-pack/NonlinearEigenproblems.jl.git",
#    target = "build",
#)


text:
struct NoBackend <: AbstractBackend end

const _backendType = Dict{Symbol,DataType}(:none => NoBackend)
const _backendSymbol = Dict{DataType,Symbol}(NoBackend => :none)
const _backends = Symbol[]
const _initialized_backends = Set{Symbol}()
const _default_backends = (:none, :gr, :plotly)

const _backend_packages = Dict{Symbol,Symbol}()

"Returns a list of supported backends"
backends() = _backends

"Returns the name of the current backend"
backend_name() = CURRENT_BACKEND.sym

function _backend_instance(sym::Symbol)::AbstractBackend
    haskey(_backendType, sym) ? _backendType[sym]() : error("Unsupported backend $sym")
end

backend_package_name(sym::Symbol) = _backend_packages[sym]

macro init_backend(s)
    package_str = string(s)
    str = lowercase(package_str)
    sym = Symbol(str)
    T = Symbol(string(s) * "Backend")
    esc(quote
        struct $T <: AbstractBackend end
        export $sym
        $sym(; kw...) = (default(; reset = false, kw...); backend($T()))
        backend_name(::$T) = Symbol($str)
        backend_package_name(pkg::$T) = backend_package_name(Symbol($str))
        push!(_backends, Symbol($str))
        _backendType[Symbol($str)] = $T
        _backendSymbol[$T] = Symbol($str)
        _backend_packages[Symbol($str)] = Symbol($package_str)
    end)
end

# include("backends/web.jl")
# include("backends/supported.jl")

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

# don't do anything as a default
_create_backend_figure(plt::Plot) = nothing
_prepare_plot_object(plt::Plot) = nothing
_initialize_subplot(plt::Plot, sp::Subplot) = nothing

_series_added(plt::Plot, series::Series) = nothing
_series_updated(plt::Plot, series::Series) = nothing

_before_layout_calcs(plt::Plot) = nothing

title_padding(sp::Subplot) = sp[:title] == "" ? 0mm : sp[:titlefontsize] * pt
guide_padding(axis::Axis) = axis[:guide] == "" ? 0mm : axis[:guidefontsize] * pt

"Returns the (width,height) of a text label."
function text_size(lablen::Int, sz::Number, rot::Number = 0)
    # we need to compute the size of the ticks generically
    # this means computing the bounding box and then getting the width/height
    # note:
    ptsz = sz * pt
    width = 0.8lablen * ptsz

    # now compute the generalized "height" after rotation as the "opposite+adjacent" of 2 triangles
    height = abs(sind(rot)) * width + abs(cosd(rot)) * ptsz
    width = abs(sind(rot + 90)) * width + abs(cosd(rot + 90)) * ptsz
    width, height
end
text_size(lab::AbstractString, sz::Number, rot::Number = 0) =
    text_size(length(lab), sz, rot)
text_size(lab::PlotText, sz::Number, rot::Number = 0) = text_size(length(lab.str), sz, rot)

# account for the size/length/rotation of tick labels
function tick_padding(sp::Subplot, axis::Axis)
    ticks = get_ticks(sp, axis)
    if ticks === nothing
        0mm
    else
        vals, labs = ticks
        isempty(labs) && return 0mm
        # ptsz = axis[:tickfont].pointsize * pt
        longest_label = maximum(length(lab) for lab in labs)

        # generalize by "rotating" y labels
        rot = axis[:rotation] + (axis[:letter] == :y ? 90 : 0)

        # # we need to compute the size of the ticks generically
        # # this means computing the bounding box and then getting the width/height
        # labelwidth = 0.8longest_label * ptsz
        #
        #
        # # now compute the generalized "height" after rotation as the "opposite+adjacent" of 2 triangles
        # hgt = abs(sind(rot)) * labelwidth + abs(cosd(rot)) * ptsz + 1mm
        # hgt

        # get the height of the rotated label
        text_size(longest_label, axis[:tickfontsize], rot)[2]
    end
end

# Set the (left, top, right, bottom) minimum padding around the plot area
# to fit ticks, tick labels, guides, colorbars, etc.
function _update_min_padding!(sp::Subplot)
    # TODO: something different when `RecipesPipeline.is3d(sp) == true`
    leftpad   = tick_padding(sp, sp[:yaxis]) + sp[:left_margin] + guide_padding(sp[:yaxis])
    toppad    = sp[:top_margin] + title_padding(sp)
    rightpad  = sp[:right_margin]
    bottompad = tick_padding(sp, sp[:xaxis]) + sp[:bottom_margin] + guide_padding(sp[:xaxis])

    # switch them?
    if sp[:xaxis][:mirror]
        bottompad, toppad = toppad, bottompad
    end
    if sp[:yaxis][:mirror]
        leftpad, rightpad = rightpad, leftpad
    end

    # @show (leftpad, toppad, rightpad, bottompad)
    sp.minpad = (leftpad, toppad, rightpad, bottompad)
end

_update_plot_object(plt::Plot) = nothing

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

mutable struct CurrentBackend
    sym::Symbol
    pkg::AbstractBackend
end
CurrentBackend(sym::Symbol) = CurrentBackend(sym, _backend_instance(sym))

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

_fallback_default_backend() = backend(GRBackend())

function _pick_default_backend()
    env_default = get(ENV, "PLOTS_DEFAULT_BACKEND", "")
    if env_default != ""
        sym = Symbol(lowercase(env_default))
        if sym in _backends
            backend(sym)
        else
            @warn(
                "You have set PLOTS_DEFAULT_BACKEND=$env_default but it is not a valid backend package.  Choose from:\n\t" *
                join(sort(_backends), "\n\t")
            )
            _fallback_default_backend()
        end
    else
        _fallback_default_backend()
    end
end

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

"""
Returns the current plotting package name.  Initializes package on first call.
"""
function backend()
    if CURRENT_BACKEND.sym == :none
        _pick_default_backend()
    end

    CURRENT_BACKEND.pkg
end

"""
Set the plot backend.
"""
function backend(pkg::AbstractBackend)
    sym = backend_name(pkg)
    if !(sym in _initialized_backends)
        _initialize_backend(pkg)
        push!(_initialized_backends, sym)
    end
    CURRENT_BACKEND.sym = sym
    CURRENT_BACKEND.pkg = pkg
    pkg
end

function backend(sym::Symbol)
    if sym in _backends
        backend(_backend_instance(sym))
    else
        @warn("`:$sym` is not a supported backend.")
        backend()
    end
end

const _deprecated_backends =
    [:qwt, :winston, :bokeh, :gadfly, :immerse, :glvisualize, :pgfplots]

function warn_on_deprecated_backend(bsym::Symbol)
    if bsym in _deprecated_backends
        if bsym == :pgfplots
            @warn("Backend $bsym has been deprecated. Use pgfplotsx instead.")
        else
            @warn("Backend $bsym has been deprecated.")
        end
    end
end

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

# these are args which every backend supports because they're not used in the backend code
const _base_supported_args = [
    :color_palette,
    :background_color,
    :background_color_subplot,
    :foreground_color,
    :foreground_color_subplot,
    :group,
    :seriestype,
    :seriescolor,
    :seriesalpha,
    :smooth,
    :xerror,
    :yerror,
    :zerror,
    :subplot,
    :x,
    :y,
    :z,
    :show,
    :size,
    :margin,
    :left_margin,
    :right_margin,
    :top_margin,
    :bottom_margin,
    :html_output_format,
    :layout,
    :link,
    :primary,
    :series_annotations,
    :subplot_index,
    :discrete_values,
    :projection,
    :show_empty_bins,
    :z_order,
]

function merge_with_base_supported(v::AVec)
    v = vcat(v, _base_supported_args)
    for vi in v
        if haskey(_axis_defaults, vi)
            for letter in (:x, :y, :z)
                push!(v, get_attr_symbol(letter, vi))
            end
        end
    end
    Set(v)
end

@init_backend PyPlot
@init_backend UnicodePlots
@init_backend Plotly
@init_backend PlotlyJS
@init_backend GR
@init_backend PGFPlots
@init_backend PGFPlotsX
@init_backend InspectDR
@init_backend HDF5
@init_backend Gaston

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

# create the various `is_xxx_supported` and `supported_xxxs` methods
# by default they pass through to checking membership in `_gr_xxx`
for s in (:attr, :seriestype, :marker, :style, :scale)
    f = Symbol("is_", s, "_supported")
    f2 = Symbol("supported_", s, "s")
    @eval begin
        $f(::AbstractBackend, $s) = false
        $f(bend::AbstractBackend, $s::AbstractVector) = all(v -> $f(bend, v), $s)
        $f($s) = $f(backend(), $s)
        $f2() = $f2(backend())
    end

    for bend in backends()
        bend_type = typeof(_backend_instance(bend))
        v = Symbol("_", bend, "_", s)
        @eval begin
            $f(::$bend_type, $s::Symbol) = $s in $v
            $f2(::$bend_type) = sort(collect($v))
        end
    end
end

# is_subplot_supported(::AbstractBackend) = false
# is_subplot_supported() = is_subplot_supported(backend())

################################################################################
# initialize the backends

function _initialize_backend(pkg::AbstractBackend)
    sym = backend_package_name(pkg)
    @eval Main begin
        import $sym
        export $sym
        Plots._check_compat($sym)
    end
end

_initialize_backend(pkg::GRBackend) = nothing

function _initialize_backend(pkg::PlotlyBackend)
    try
        @eval Main begin
            import PlotlyBase
        end
        _check_compat(PlotlyBase)
    catch
        @info "For saving to png with the Plotly backend PlotlyBase has to be installed."
    end
end

# ------------------------------------------------------------------------------
# gr

const _gr_attr = merge_with_base_supported([
    :annotations,
    :legend_background_color,
    :background_color_inside,
    :background_color_outside,
    :legend_foreground_color,
    :foreground_color_grid,
    :foreground_color_axis,
    :foreground_color_text,
    :foreground_color_border,
    :label,
    :seriescolor,
    :seriesalpha,
    :linecolor,
    :linestyle,
    :linewidth,
    :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    :markeralpha,
    :markerstrokewidth,
    :markerstrokecolor,
    :markerstrokealpha,
    :fillrange,
    :fillcolor,
    :fillalpha,
    :bins,
    :layout,
    :title,
    :window_title,
    :guide,
    :lims,
    :ticks,
    :scale,
    :flip,
    :titlefontfamily,
    :titlefontsize,
    :titlefonthalign,
    :titlefontvalign,
    :titlefontrotation,
    :titlefontcolor,
    :legend_font_family,
    :legend_font_pointsize,
    :legend_font_halign,
    :legend_font_valign,
    :legend_font_rotation,
    :legend_font_color,
    :tickfontfamily,
    :tickfontsize,
    :tickfonthalign,
    :tickfontvalign,
    :tickfontrotation,
    :tickfontcolor,
    :guidefontfamily,
    :guidefontsize,
    :guidefonthalign,
    :guidefontvalign,
    :guidefontrotation,
    :guidefontcolor,
    :grid,
    :gridalpha,
    :gridstyle,
    :gridlinewidth,
    :legend_position,
    :legend_title,
    :colorbar,
    :colorbar_title,
    :colorbar_entry,
    :colorbar_scale,
    :clims,
    :fill_z,
    :line_z,
    :marker_z,
    :levels,
    :ribbon,
    :quiver,
    :orientation,
    :overwrite_figure,
    :polar,
    :aspect_ratio,
    :normalize,
    :weights,
    :inset_subplots,
    :bar_width,
    :arrow,
    :framestyle,
    :tick_direction,
    :camera,
    :contour_labels,
    :connections,
])
const _gr_seriestype = [
    :path,
    :scatter,
    :straightline,
    :heatmap,
    :image,
    :contour,
    :path3d,
    :scatter3d,
    :surface,
    :wireframe,
    :mesh3d,
    :volume,
    :shape,
]
const _gr_style = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
const _gr_marker = _allMarkers
const _gr_scale = [:identity, :log10]
is_marker_supported(::GRBackend, shape::Shape) = true

# ------------------------------------------------------------------------------
# plotly

const _plotly_attr = merge_with_base_supported([
    :annotations,
    :legend_background_color,
    :background_color_inside,
    :background_color_outside,
    :legend_foreground_color,
    :foreground_color_guide,
    :foreground_color_grid,
    :foreground_color_axis,
    :foreground_color_text,
    :foreground_color_border,
    :foreground_color_title,
    :label,
    :seriescolor,
    :seriesalpha,
    :linecolor,
    :linestyle,
    :linewidth,
    :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    :markeralpha,
    :markerstrokewidth,
    :markerstrokecolor,
    :markerstrokealpha,
    :markerstrokestyle,
    :fillrange,
    :fillcolor,
    :fillalpha,
    :bins,
    :title,
    :titlelocation,
    :titlefontfamily,
    :titlefontsize,
    :titlefonthalign,
    :titlefontvalign,
    :titlefontcolor,
    :legend_font_family,
    :legend_font_pointsize,
    :legend_font_color,
    :tickfontfamily,
    :tickfontsize,
    :tickfontcolor,
    :guidefontfamily,
    :guidefontsize,
    :guidefontcolor,
    :window_title,
    :guide,
    :lims,
    :ticks,
    :scale,
    :flip,
    :rotation,
    :tickfont,
    :guidefont,
    :legendfont,
    :grid,
    :gridalpha,
    :gridlinewidth,
    :legend,
    :colorbar,
    :colorbar_title,
    :colorbar_entry,
    :marker_z,
    :fill_z,
    :line_z,
    :levels,
    :ribbon,
    :quiver,
    :orientation,
    # :overwrite_figure,
    :polar,
    :normalize,
    :weights,
    # :contours,
    :aspect_ratio,
    :hover,
    :inset_subplots,
    :bar_width,
    :clims,
    :framestyle,
    :tick_direction,
    :camera,
    :contour_labels,
    :connections,
])

const _plotly_seriestype = [
    :path,
    :scatter,
    :heatmap,
    :contour,
    :surface,
    :wireframe,
    :path3d,
    :scatter3d,
    :shape,
    :scattergl,
    :straightline,
    :mesh3d,
]
const _plotly_style = [:auto, :solid, :dash, :dot, :dashdot]
const _plotly_marker = [
    :none,
    :auto,
    :circle,
    :rect,
    :diamond,
    :utriangle,
    :dtriangle,
    :cross,
    :xcross,
    :pentagon,
    :hexagon,
    :octagon,
    :vline,
    :hline,
]
const _plotly_scale = [:identity, :log10]

defaultOutputFormat(plt::Plot{Plots.PlotlyBackend}) = "html"

# ------------------------------------------------------------------------------
# pgfplots

const _pgfplots_attr = merge_with_base_supported([
    :annotations,
    :legend_background_color,
    :background_color_inside,
    # :background_color_outside,
    # :legend_foreground_color,
    :foreground_color_grid,
    :foreground_color_axis,
    :foreground_color_text,
    :foreground_color_border,
    :label,
    :seriescolor,
    :seriesalpha,
    :linecolor,
    :linestyle,
    :linewidth,
    :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    :markeralpha,
    :markerstrokewidth,
    :markerstrokecolor,
    :markerstrokealpha,
    :markerstrokestyle,
    :fillrange,
    :fillcolor,
    :fillalpha,
    :bins,
    # :bar_width, :bar_edges,
    :title,
    # :window_title,
    :guide,
    :guide_position,
    :lims,
    :ticks,
    :scale,
    :flip,
    :rotation,
    :tickfont,
    :guidefont,
    :legendfont,
    :grid,
    :legend,
    :colorbar,
    :colorbar_title,
    :fill_z,
    :line_z,
    :marker_z,
    :levels,
    # :ribbon, :quiver, :arrow,
    # :orientation,
    # :overwrite_figure,
    :polar,
    # :normalize, :weights, :contours,
    :aspect_ratio,
    :tick_direction,
    :framestyle,
    :camera,
    :contour_labels,
])
const _pgfplots_seriestype = [
    :path,
    :path3d,
    :scatter,
    :steppre,
    :stepmid,
    :steppost,
    :histogram2d,
    :ysticks,
    :xsticks,
    :contour,
    :shape,
    :straightline,
]
const _pgfplots_style = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
const _pgfplots_marker = [
    :none,
    :auto,
    :circle,
    :rect,
    :diamond,
    :utriangle,
    :dtriangle,
    :cross,
    :xcross,
    :star5,
    :pentagon,
    :hline,
    :vline,
] #vcat(_allMarkers, Shape)
const _pgfplots_scale = [:identity, :ln, :log2, :log10]

# ------------------------------------------------------------------------------
# plotlyjs

function _initialize_backend(pkg::PlotlyJSBackend)
    @eval Main begin
        import PlotlyJS
        export PlotlyJS
    end
end

const _plotlyjs_attr       = _plotly_attr
const _plotlyjs_seriestype = _plotly_seriestype
const _plotlyjs_style      = _plotly_style
const _plotlyjs_marker     = _plotly_marker
const _plotlyjs_scale      = _plotly_scale

# ------------------------------------------------------------------------------
# pyplot

function _initialize_backend(::PyPlotBackend)
    @eval Main begin
        import PyPlot

        export PyPlot

        # we don't want every command to update the figure
        PyPlot.ioff()
    end
end

const _pyplot_attr = merge_with_base_supported([
    :annotations,
    :legend_background_color,
    :background_color_inside,
    :background_color_outside,
    :foreground_color_grid,
    :legend_foreground_color,
    :foreground_color_title,
    :foreground_color_axis,
    :foreground_color_border,
    :foreground_color_guide,
    :foreground_color_text,
    :label,
    :linecolor,
    :linestyle,
    :linewidth,
    :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    :markeralpha,
    :markerstrokewidth,
    :markerstrokecolor,
    :markerstrokealpha,
    :fillrange,
    :fillcolor,
    :fillalpha,
    :bins,
    :bar_width,
    :bar_edges,
    :bar_position,
    :title,
    :titlelocation,
    :titlefont,
    :window_title,
    :guide,
    :guide_position,
    :lims,
    :ticks,
    :scale,
    :flip,
    :rotation,
    :titlefontfamily,
    :titlefontsize,
    :titlefontcolor,
    :legend_font_family,
    :legend_font_pointsize,
    :legend_font_color,
    :tickfontfamily,
    :tickfontsize,
    :tickfontcolor,
    :guidefontfamily,
    :guidefontsize,
    :guidefontcolor,
    :grid,
    :gridalpha,
    :gridstyle,
    :gridlinewidth,
    :legend_position,
    :legend_title,
    :colorbar,
    :colorbar_title,
    :colorbar_entry,
    :colorbar_ticks,
    :colorbar_tickfontfamily,
    :colorbar_tickfontsize,
    :colorbar_tickfonthalign,
    :colorbar_tickfontvalign,
    :colorbar_tickfontrotation,
    :colorbar_tickfontcolor,
    :colorbar_scale,
    :marker_z,
    :line_z,
    :fill_z,
    :levels,
    :ribbon,
    :quiver,
    :arrow,
    :orientation,
    :overwrite_figure,
    :polar,
    :normalize,
    :weights,
    :contours,
    :aspect_ratio,
    :clims,
    :inset_subplots,
    :dpi,
    :stride,
    :framestyle,
    :tick_direction,
    :camera,
    :contour_labels,
    :connections,
])
const _pyplot_seriestype = [
    :path,
    :steppre,
    :stepmid,
    :steppost,
    :shape,
    :straightline,
    :scatter,
    :hexbin,
    :heatmap,
    :image,
    :contour,
    :contour3d,
    :path3d,
    :scatter3d,
    :mesh3d,
    :surface,
    :wireframe,
]
const _pyplot_style = [:auto, :solid, :dash, :dot, :dashdot]
const _pyplot_marker = vcat(_allMarkers, :pixel)
const _pyplot_scale = [:identity, :ln, :log2, :log10]

# ------------------------------------------------------------------------------
# Gaston

function _initialize_backend(::GastonBackend)
    @eval Main begin
        import Gaston
        export Gaston
    end
end

const _gaston_attr = merge_with_base_supported([
    :annotations,
    # :background_color_legend,
    # :background_color_inside,
    # :background_color_outside,
    # :foreground_color_legend,
    # :foreground_color_grid, :foreground_color_axis,
    # :foreground_color_text, :foreground_color_border,
    :label,
    :seriescolor,
    :seriesalpha,
    :linecolor,
    :linestyle,
    :linewidth,
    :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    :markeralpha,
    # :markerstrokewidth, :markerstrokecolor, :markerstrokealpha, :markerstrokestyle,
    # :fillrange, :fillcolor, :fillalpha,
    # :bins,
    # :bar_width, :bar_edges,
    :title,
    :window_title,
    :guide,
    :guide_position,
    :lims,
    :ticks,
    :scale,
    :flip,
    :rotation,
    :tickfont,
    :guidefont,
    :legendfont,
    :grid,
    :legend,
    # :colorbar, :colorbar_title,
    # :fill_z, :line_z, :marker_z, :levels,
    # :ribbon,
    :quiver,
    :arrow,
    # :orientation, :overwrite_figure,
    :polar,
    # :normalize, :weights, :contours,
    :aspect_ratio,
    :tick_direction,
    # :framestyle,
    # :camera,
    # :contour_labels,
    :connections,
])

const _gaston_seriestype = [
    :path,
    :path3d,
    :scatter,
    :steppre,
    :stepmid,
    :steppost,
    :ysticks,
    :xsticks,
    :contour,
    :shape,
    :straightline,
    :scatter3d,
    :contour3d,
    :wireframe,
    :heatmap,
    :surface,
    :mesh3d,
    :image,
]

const _gaston_style = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]

const _gaston_marker = [
    :none,
    :auto,
    :pixel,
    :cross,
    :xcross,
    :+,
    :x,
    :star5,
    :rect,
    :circle,
    :utriangle,
    :dtriangle,
    :diamond,
    :pentagon,
    # :hline,
    # :vline,
]

const _gaston_scale = [:identity, :ln, :log2, :log10]

# ------------------------------------------------------------------------------
# unicodeplots

const _unicodeplots_attr = merge_with_base_supported([
    :annotations,
    :bins,
    :guide,
    :grid,
    :label,
    :layout,
    :legend,
    :lims,
    :linealpha,
    :linecolor,
    :linestyle,
    :markershape,
    :quiver,
    :arrow,
    :seriesalpha,
    :seriescolor,
    :scale,
    :flip,
    :title,
    # :marker_z,
    :line_z,
])
const _unicodeplots_seriestype = [
    :path,
    :path3d,
    :scatter,
    :scatter3d,
    :straightline,
    # :bar,
    :shape,
    :histogram2d,
    :heatmap,
    :contour,
    # :contour3d,
    :spy,
    :surface,
    :wireframe,
    :mesh3d,
]
const _unicodeplots_style = [:auto, :solid]
const _unicodeplots_marker = [
    :none,
    :auto,
    :pixel,
    # vvvvvvvvvv shapes
    :circle,
    :rect,
    :star5,
    :diamond,
    :hexagon,
    :cross,
    :xcross,
    :utriangle,
    :dtriangle,
    :rtriangle,
    :ltriangle,
    :pentagon,
    # :heptagon,
    # :octagon,
    :star4,
    :star6,
    # :star7,
    :star8,
    :vline,
    :hline,
    :+,
    :x,
]
const _unicodeplots_scale = [:identity, :ln, :log2, :log10]

# ------------------------------------------------------------------------------
# hdf5

const _hdf5_attr = merge_with_base_supported([
    :annotations,
    :legend_background_color,
    :background_color_inside,
    :background_color_outside,
    :foreground_color_grid,
    :legend_foreground_color,
    :foreground_color_title,
    :foreground_color_axis,
    :foreground_color_border,
    :foreground_color_guide,
    :foreground_color_text,
    :label,
    :linecolor,
    :linestyle,
    :linewidth,
    :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    :markeralpha,
    :markerstrokewidth,
    :markerstrokecolor,
    :markerstrokealpha,
    :fillrange,
    :fillcolor,
    :fillalpha,
    :bins,
    :bar_width,
    :bar_edges,
    :bar_position,
    :title,
    :titlelocation,
    :titlefont,
    :window_title,
    :guide,
    :lims,
    :ticks,
    :scale,
    :flip,
    :rotation,
    :tickfont,
    :guidefont,
    :legendfont,
    :grid,
    :legend,
    :colorbar,
    :marker_z,
    :line_z,
    :fill_z,
    :levels,
    :ribbon,
    :quiver,
    :arrow,
    :orientation,
    :overwrite_figure,
    :polar,
    :normalize,
    :weights,
    :contours,
    :aspect_ratio,
    :clims,
    :inset_subplots,
    :dpi,
    :colorbar_title,
])
const _hdf5_seriestype = [
    :path,
    :steppre,
    :stepmid,
    :steppost,
    :shape,
    :straightline,
    :scatter,
    :hexbin,
    :heatmap,
    :image,
    :contour,
    :contour3d,
    :path3d,
    :scatter3d,
    :surface,
    :wireframe,
]
const _hdf5_style = [:auto, :solid, :dash, :dot, :dashdot]
const _hdf5_marker = vcat(_allMarkers, :pixel)
const _hdf5_scale = [:identity, :ln, :log2, :log10]

# Additional constants
#Dict has problems using "Types" as keys.  Initialize in "_initialize_backend":
const HDF5PLOT_MAP_STR2TELEM = Dict{String,Type}()
const HDF5PLOT_MAP_TELEM2STR = Dict{Type,String}()

#Don't really like this global variable... Very hacky
mutable struct HDF5Plot_PlotRef
    ref::Union{Plot,Nothing}
end
const HDF5PLOT_PLOTREF = HDF5Plot_PlotRef(nothing)

# ------------------------------------------------------------------------------
# inspectdr

const _inspectdr_attr = merge_with_base_supported([
    :annotations,
    :legend_background_color,
    :background_color_inside,
    :background_color_outside,
    # :foreground_color_grid,
    :legend_foreground_color,
    :foreground_color_title,
    :foreground_color_axis,
    :foreground_color_border,
    :foreground_color_guide,
    :foreground_color_text,
    :label,
    :seriescolor,
    :seriesalpha,
    :linecolor,
    :linestyle,
    :linewidth,
    :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    :markeralpha,
    :markerstrokewidth,
    :markerstrokecolor,
    :markerstrokealpha,
    :markerstrokestyle, #Causes warning not to have it... what is this?
    :fillcolor,
    :fillalpha, #:fillrange,
    #    :bins, :bar_width, :bar_edges, :bar_position,
    :title,
    :titlelocation,
    :window_title,
    :guide,
    :lims,
    :scale, #:ticks, :flip, :rotation,
    :titlefontfamily,
    :titlefontsize,
    :titlefontcolor,
    :legend_font_family,
    :legend_font_pointsize,
    :legend_font_color,
    :tickfontfamily,
    :tickfontsize,
    :tickfontcolor,
    :guidefontfamily,
    :guidefontsize,
    :guidefontcolor,
    :grid,
    :legend_position, #:colorbar,
    #    :marker_z,
    #    :line_z,
    #    :levels,
    #   :ribbon, :quiver, :arrow,
    #    :orientation,
    :overwrite_figure,
    :polar,
    #    :normalize, :weights,
    #    :contours, :aspect_ratio,
    #    :clims,
    #    :inset_subplots,
    :dpi,
    #    :colorbar_title,
])
const _inspectdr_style = [:auto, :solid, :dash, :dot, :dashdot]
const _inspectdr_seriestype = [
    :path,
    :scatter,
    :shape,
    :straightline, #, :steppre, :stepmid, :steppost
]
#see: _allMarkers, _shape_keys
const _inspectdr_marker = Symbol[
    :none,
    :auto,
    :circle,
    :rect,
    :diamond,
    :cross,
    :xcross,
    :utriangle,
    :dtriangle,
    :rtriangle,
    :ltriangle,
    :pentagon,
    :hexagon,
    :heptagon,
    :octagon,
    :star4,
    :star5,
    :star6,
    :star7,
    :star8,
    :vline,
    :hline,
    :+,
    :x,
]

const _inspectdr_scale = [:identity, :ln, :log2, :log10]
# ------------------------------------------------------------------------------
# pgfplotsx

const _pgfplotsx_attr = merge_with_base_supported([
    :annotations,
    :legend_background_color,
    :background_color_inside,
    :background_color_outside,
    :legend_foreground_color,
    :foreground_color_grid,
    :foreground_color_axis,
    :foreground_color_text,
    :foreground_color_border,
    :label,
    :seriescolor,
    :seriesalpha,
    :linecolor,
    :linestyle,
    :linewidth,
    :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    :markeralpha,
    :markerstrokewidth,
    :markerstrokecolor,
    :markerstrokealpha,
    :fillrange,
    :fillcolor,
    :fillalpha,
    :bins,
    :layout,
    :title,
    :window_title,
    :guide,
    :lims,
    :ticks,
    :scale,
    :flip,
    :titlefontfamily,
    :titlefontsize,
    :titlefonthalign,
    :titlefontvalign,
    :titlefontrotation,
    :titlefontcolor,
    :legend_font_family,
    :legend_font_pointsize,
    :legend_font_halign,
    :legend_font_valign,
    :legend_font_rotation,
    :legend_font_color,
    :tickfontfamily,
    :tickfontsize,
    :tickfonthalign,
    :tickfontvalign,
    :tickfontrotation,
    :tickfontcolor,
    :guidefontfamily,
    :guidefontsize,
    :guidefonthalign,
    :guidefontvalign,
    :guidefontrotation,
    :guidefontcolor,
    :grid,
    :gridalpha,
    :gridstyle,
    :gridlinewidth,
    :legend_position,
    :legend_title,
    :colorbar,
    :colorbar_title,
    :colorbar_entry,
    :fill_z,
    :line_z,
    :marker_z,
    :levels,
    :ribbon,
    :quiver,
    :orientation,
    :overwrite_figure,
    :polar,
    :aspect_ratio,
    :normalize,
    :weights,
    :inset_subplots,
    :bar_width,
    :arrow,
    :framestyle,
    :tick_direction,
    :camera,
    :contour_labels,
    :connections,
])
const _pgfplotsx_seriestype = [
    :path,
    :scatter,
    :straightline,
    :path3d,
    :scatter3d,
    :surface,
    :wireframe,
    :heatmap,
    :mesh3d,
    :contour,
    :contour3d,
    :quiver,
    :shape,
    :steppre,
    :stepmid,
    :steppost,
    :ysticks,
    :xsticks,
]
const _pgfplotsx_style = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
const _pgfplotsx_marker = [
    :none,
    :auto,
    :circle,
    :rect,
    :diamond,
    :utriangle,
    :dtriangle,
    :ltriangle,
    :rtriangle,
    :cross,
    :xcross,
    :x,
    :star5,
    :pentagon,
    :hline,
    :vline,
]
const _pgfplotsx_scale = [:identity, :ln, :log2, :log10]
is_marker_supported(::PGFPlotsXBackend, shape::Shape) = true

# additional constants
const _pgfplotsx_series_ids = KW()


text:
"""
    centroids_radii(points::Array{Float64, 2}, indices_simplices::Array{Int, 1})
Compute the centroids and radii of the simplices in the triangulation in any dimension.

`points::Array{Float64, 2}` Points furnishing the triangulation (size = npoints x dim).
`indices_simplices::Array{Int, 1}` Indices of the vertices furnishing simplices of the
    triangulation. Each row refers to one simplex.
"""
function centroids_radii2(points, indices_simplices)
    dim = size(points, 2)
    nsimplices = size(indices_simplices, 1)
    centroids = zeros(nsimplices, dim)
    radii = zeros(nsimplices)

    for i = 1:nsimplices
        simplex = points[indices_simplices[i, :], :] # (dim + 1) x dim)
        centroid = sum(simplex, dims=1) / (dim + 1)
        centroid_matrix = repeat(centroid, dim + 1, 1)

        # Subtract centroid from each simplex
        subtracted = simplex - centroid_matrix
        radius = sqrt(maximum(sum(subtracted.^2, dims=2)))

        centroids[i, :] = centroid
        radii[i] = radius

    end

    return centroids, radii
end


text:
# This file is a part of NeighborhoodApproximationIndex.jl

using SimilaritySearch
using Test

function test_knr_vectors(perf, knr, lower)
    @time p = probe(perf, knr)
    @test p.macrorecall >= lower
    p
end

@testset "indexing vectors with Knr" begin
    # NOTE: The following algorithms are complex enough to say we are testing
    # it doesn't have syntax errors and achieves some expected quality
    # a more grained test functions are required
    k = 10
    n = 10000 # number of items in the dataset
    m = 100  # number of queries
    dim = 8  # vector's dimension
    db = [rand(Float32, dim) for i in 1:n]
    queries = [rand(Float32, dim) for i in 1:m]

    dist = L2Distance()
    seq = ExhaustiveSearch(dist, db)
    perf = Performance(seq, queries, k)

    knr = Knr(dist, db; numrefs=64, kbuild=2)
    p = test_knr_vectors(perf, knr, 0.7)
    @info "-- $(string(knr)) $(p)"
    optimize!(perf, knr; recall=0.95)
    p = test_knr_vectors(perf, knr, 0.95)
    @info "-- Optimized $(string(knr)) $(knr.opts) -- $(p)"
end


text:
include("balance_law_interface.jl")

using ..Courant
"""
    calculate_dt(dg, model, Q, Courant_number, direction, t)

For a given model, compute a time step satisying the nondiffusive Courant number
`Courant_number`
"""
function calculate_dt(dg, model, Q, Courant_number, t, direction)
    Δt = one(eltype(Q))
    CFL = courant(nondiffusive_courant, dg, model, Q, Δt, t, direction)
    return Courant_number / CFL
end

function create_conservative_state(balance_law::BalanceLaw, grid)
    topology = grid.topology
    # FIXME: Remove after updating CUDA
    h_vgeo = Array(grid.vgeo)
    FT = eltype(h_vgeo)
    Np = dofs_per_element(grid)
    DA = arraytype(grid)

    weights = view(h_vgeo, :, grid.Mid, :)
    weights = reshape(weights, size(weights, 1), 1, size(weights, 2))

    V = vars_state_conservative(balance_law, FT)
    state_conservative = MPIStateArray{FT, V}(
        topology.mpicomm,
        DA,
        Np,
        number_state_conservative(balance_law, FT),
        length(topology.elems),
        realelems = topology.realelems,
        ghostelems = topology.ghostelems,
        vmaprecv = grid.vmaprecv,
        vmapsend = grid.vmapsend,
        nabrtorank = topology.nabrtorank,
        nabrtovmaprecv = grid.nabrtovmaprecv,
        nabrtovmapsend = grid.nabrtovmapsend,
        weights = weights,
    )
    return state_conservative
end

function create_auxiliary_state(balance_law, grid)
    topology = grid.topology
    Np = dofs_per_element(grid)

    h_vgeo = Array(grid.vgeo)
    FT = eltype(h_vgeo)
    DA = arraytype(grid)

    weights = view(h_vgeo, :, grid.Mid, :)
    weights = reshape(weights, size(weights, 1), 1, size(weights, 2))

    V = vars_state_auxiliary(balance_law, FT)
    state_auxiliary = MPIStateArray{FT, V}(
        topology.mpicomm,
        DA,
        Np,
        number_state_auxiliary(balance_law, FT),
        length(topology.elems),
        realelems = topology.realelems,
        ghostelems = topology.ghostelems,
        vmaprecv = grid.vmaprecv,
        vmapsend = grid.vmapsend,
        nabrtorank = topology.nabrtorank,
        nabrtovmaprecv = grid.nabrtovmaprecv,
        nabrtovmapsend = grid.nabrtovmapsend,
        weights = weights,
    )

    dim = dimensionality(grid)
    polyorder = polynomialorder(grid)
    vgeo = grid.vgeo
    device = array_device(state_auxiliary)
    nrealelem = length(topology.realelems)
    event = Event(device)
    event = kernel_init_state_auxiliary!(device, min(Np, 1024), Np * nrealelem)(
        balance_law,
        Val(dim),
        Val(polyorder),
        state_auxiliary.data,
        vgeo,
        topology.realelems,
        dependencies = (event,),
    )
    event = MPIStateArrays.begin_ghost_exchange!(
        state_auxiliary;
        dependencies = event,
    )
    event = MPIStateArrays.end_ghost_exchange!(
        state_auxiliary;
        dependencies = event,
    )
    wait(device, event)

    return state_auxiliary
end

function create_gradient_state(balance_law, grid)
    topology = grid.topology
    Np = dofs_per_element(grid)

    h_vgeo = Array(grid.vgeo)
    FT = eltype(h_vgeo)
    DA = arraytype(grid)

    weights = view(h_vgeo, :, grid.Mid, :)
    weights = reshape(weights, size(weights, 1), 1, size(weights, 2))

    # TODO: Clean up this MPIStateArray interface...
    V = vars_state_gradient_flux(balance_law, FT)
    state_gradient_flux = MPIStateArray{FT, V}(
        topology.mpicomm,
        DA,
        Np,
        number_state_gradient_flux(balance_law, FT),
        length(topology.elems),
        realelems = topology.realelems,
        ghostelems = topology.ghostelems,
        vmaprecv = grid.vmaprecv,
        vmapsend = grid.vmapsend,
        nabrtorank = topology.nabrtorank,
        nabrtovmaprecv = grid.nabrtovmaprecv,
        nabrtovmapsend = grid.nabrtovmapsend,
        weights = weights,
    )

    return state_gradient_flux
end

function create_higher_order_states(balance_law, grid)
    topology = grid.topology
    Np = dofs_per_element(grid)

    h_vgeo = Array(grid.vgeo)
    FT = eltype(h_vgeo)
    DA = arraytype(grid)

    weights = view(h_vgeo, :, grid.Mid, :)
    weights = reshape(weights, size(weights, 1), 1, size(weights, 2))

    ngradlapstate = num_gradient_laplacian(balance_law, FT)
    # TODO: Clean up this MPIStateArray interface...
    V = vars_gradient_laplacian(balance_law, FT)
    Qhypervisc_grad = MPIStateArray{FT, V}(
        topology.mpicomm,
        DA,
        Np,
        3ngradlapstate,
        length(topology.elems),
        realelems = topology.realelems,
        ghostelems = topology.ghostelems,
        vmaprecv = grid.vmaprecv,
        vmapsend = grid.vmapsend,
        nabrtorank = topology.nabrtorank,
        nabrtovmaprecv = grid.nabrtovmaprecv,
        nabrtovmapsend = grid.nabrtovmapsend,
        weights = weights,
    )

    V = vars_hyperdiffusive(balance_law, FT)
    Qhypervisc_div = MPIStateArray{FT, V}(
        topology.mpicomm,
        DA,
        Np,
        ngradlapstate,
        length(topology.elems),
        realelems = topology.realelems,
        ghostelems = topology.ghostelems,
        vmaprecv = grid.vmaprecv,
        vmapsend = grid.vmapsend,
        nabrtorank = topology.nabrtorank,
        nabrtovmaprecv = grid.nabrtovmaprecv,
        nabrtovmapsend = grid.nabrtovmapsend,
        weights = weights,
    )
    return Qhypervisc_grad, Qhypervisc_div
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
"""
Abstract type to represent the structure and interconnectedness of the system
"""
abstract type Topology <: Component end

"""
Represents a geographical region of system components.

All subtypes must implement the method `get_aggregation_topology_accessor`.
"""
abstract type AggregationTopology <: Topology end

"""
Return the method to be called on a Bus to get its AggregationTopology value for this type.
"""
function get_aggregation_topology_accessor(::Type{T}) where {T <: AggregationTopology}
    error("get_aggregation_topology_accessor must be implemented for $T")
    return
end

function check_bus_params(
    number,
    name,
    bustype,
    angle,
    voltage,
    voltage_limits,
    base_voltage,
    area,
    load_zone,
    ext,
    internal,
)
    if !isnothing(bustype)
        if bustype == BusTypes.SLACK
            bustype = BusTypes.REF
            @debug "Changed bus type from SLACK to" _group = IS.LOG_GROUP_SYSTEM bustype
        elseif bustype == BusTypes.ISOLATED
            throw(DataFormatError("isolated buses are not supported; name=$name"))
        end
    end

    return number,
    name,
    bustype,
    angle,
    voltage,
    voltage_limits,
    base_voltage,
    area,
    load_zone,
    ext,
    internal
end


text:
# -*- coding: utf-8;-*-
using Distributions
using StatsBase
using LightGraphs
using DataStructures

import Base: rand

export GraphGenerator, SocialNE, GEOGD, SFGD, nodearrival

abstract type GraphGenerator end

"""
    rand([rng::AbstractRNG], gf::GraphGenerator, Ntot::Integer, Tmax::Number, N0::Integer=5)

# Arguments
- `Ntot` : total # of nodes
- `N0` : # of nodes at time 0
- `Tmax` : end of timespan

Generates a random network depending on the `GraphGenerator`
"""
rand(gf::GraphGenerator, Ntot::Integer, Tmax::Number, N0::Integer=5) =
    rand(Base.GLOBAL_RNG, gf, Ntot, Tmax, N0)

"""
Node arrival functions for the dynamic networks
Ntot -- total numnber of nodes over all time
N0 -- initial number of nodes at time 0
Tmax -- max time
Returns t(i) for i=1..Ntot, where t(i) is the arrival time for node i
t(i) is sorted from smallest to largest
    """
nodearrival(rng::AbstractRNG, af::Symbol, Ntot::Integer, N0::Integer, Tmax::Number) =
    nodearrival(Val{af}(), Ntot, N0, Tmax)
nodearrival(af::Symbol, Ntot::Integer, N0::Integer, Tmax::Number) =
    nodearrival(Base.GLOBAL_RNG, Val{af}(), Ntot, N0, Tmax)

"""Constant node arrival rate
"""
function nodearrival(rng::AbstractRNG, af::Val{:constant}, Ntot::Integer, N0::Integer, Tmax::Number)
    nodes = zeros(Float64, Ntot)
    nodes[N0+1:end] = rand(rng, Float64, Ntot-N0) .* Tmax
    nodes
end

function newton(x0::Number, f, fprime, args::Tuple=();
                tol::AbstractFloat=1e-8, maxiter::Integer=50)
    for iter in 1:maxiter
        yprime = fprime(x0, args...)
        if abs(yprime) < eps(typeof(yprime))
            warn("First derivative is zero")
            return x0
        end
        y = f(x0, args...)
        x1 = x0 - y/yprime
        if abs(x1-x0) < tol
            return x1
        end
        x0 = x1
    end
    error("Max iteration exceeded")
end

"""Linear node arrival rate
pdf(x) = B0 + 2(1 - B0) x
cdf(x) = B0 x + (1 - B0) x^2
where x is time
"""
function nodearrival(rng::AbstractRNG, af::Val{:linear}, Ntot::Integer, N0::Integer, Tmax::Number)
    nodes = zeros(Float64,Ntot)
    # using inversion sampling
    B0 = N0/Ntot
    f(x,u) = B0*x + (1 - B0) * x^2 - u
    fp(x,u) = B0 + 2*(1 - B0) * x
    for i = N0+1:Ntot
        u = rand(rng, Float64)
        x = newton(0.5, f, fp, (u,))
        nodes[i] = x * Tmax
    end
    # nodes[N0+1:Ntot] = sqrt.(rand(rng, Float64, Ntot-N0)) .* Tmax
    nodes
end

"""Quadratic node arrival rate
pdf(x) = B0 + 3(1 - B0) x^2
cdf(x) = B0 x + (1 - B0) x^3
"""
function nodearrival(rng::AbstractRNG, af::Val{:quad}, Ntot::Integer, N0::Integer, Tmax::Number)
    nodes = zeros(Float64,Ntot)
    # using inversion sampling
    B0 = N0/Ntot
    f(x,u) = B0*x + (1-B0)*x^3 - u
    fp(x,u) = B0 + 3*(1-B0)*x^2
    for i = N0+1:Ntot
        u = rand(rng, Float64)
        x = newton(0.5, f, fp, (u,))
        nodes[i] = x * Tmax
    end
    # nodes[N0+1:Ntot] = (rand(rng, Float64, Ntot-N0) .^ (1/3)) .* Tmax
    nodes
end

"""Exponential node arrival rate
"""
function nodearrival(rng::AbstractRNG, af::Val{:exp}, Ntot::Integer, N0::Integer, Tmax::Number)
    nodes = zeros(Float64,Ntot)
    nodes[N0+1:Ntot] = log.(1 .+ rand(rng, Float64, Ntot-N0)) ./ log(2) .* Tmax
    nodes
end

function init_nodearrival_pq(nodetimes)
    # node (name, isalive, birthtime, lifetime) -> node waketime
    Q = PriorityQueue{Tuple{Int,Bool,Float64,Float64},Float64}()
    u = 1
    for t in nodetimes
        enqueue!(Q, (u, false, t, 0.0), t)
        u += 1
    end
    Q
end

"Pick a direction in 3D uniformly at random"
function sample_direction(rng::AbstractRNG)
    θ = 2π * rand(rng)
    z = 2 * rand(rng) - 1
    s = sqrt(1-z^2)
    (s * cos(θ), s * sin(θ), z)
end

"Pick a 3D point in a sphere uniformly at random"
function sample_sphere_volume(rng::AbstractRNG, R)
    (x,y,z) = sample_direction(rng)
    r = R*(rand(rng)^(1/3))
    (r*x, r*y, r*z)
end

"Euclidean distance between two points in 3D"
dist3d(x,y) = sqrt((x[1]-y[1])^2 + (x[2]-y[2])^2 + (x[3]-y[3])^2)

"""
Geometric gene duplication with probability cutoff

- `p` : probability cutoff
- `ε` : distance (set this to 1)
- `arrival` : node arrival function (:quad, :linear, :exp, :constant)

Przulj, N., Kuchaiev, O., Stevanovic, A., and Hayes, W. (2010). Geometric evolutionary dynamics of protein interaction networks. In Proc. of the Pacific Symposium Biocomputing, pages 4–8.
"""
immutable GEOGD <: GraphGenerator
    p::Float64
    ε::Float64
    arrival::Symbol
end
function rand(rng::AbstractRNG, gf::GEOGD, Ntot::Integer, Tmax::Number, N0::Integer=5)
    nodetimes = nodearrival(gf.arrival, Ntot, N0, Tmax-1)
    # initialize nodes
    Q = init_nodearrival_pq(nodetimes)
    nnodes = length(Q)
    I = Int[]
    J = Int[]
    events = Events[]
    points = Tuple{Int,Tuple{Float64,Float64,Float64}}[] # should be a kdtree or sth but hey
    # main loop
    while !isempty(Q)
        ((u,isalive,birthtime,lifetime),curtime) = dequeue_pair!(Q)
        if curtime == 0
            locu = sample_sphere_volume(rng, gf.ε/2)
        else
            (v,(vx,vy,vz)) = rand(rng, points)
            (x,y,z) = sample_direction(rng)
            if rand(rng) < gf.p
                r = gf.ε * rand(rng)
            else
                r = 10gf.ε * rand(rng)
            end
            locu = (vx+r*x,vy+r*y,vz+r*z)
        end
        push!(points, (u, locu))
        for i = 1:length(points)-1
            (v,locv) = points[i]
            if dist3d(locu,locv) < gf.ε
                push!(I, u)
                push!(J, v)
                push!(events, Events([(curtime, Tmax)]))
            end
        end
    end
    sparse(vcat(I,J), vcat(J,I), vcat(events,events), nnodes, nnodes, fixevents)
end

function add_event!(E::Dict, u, v, t_s, t_e)
    u, v = min(u,v), max(u,v)
    if (u,v) in keys(E)
        push!(E[(u,v)], (t_s,t_e))
    else
        E[(u,v)] = [(t_s,t_e)]
    end
end

function stop_event!(E::Dict, u, v, t_enew)
    u, v = min(u,v), max(u,v)
    (t_s,t_e) = E[(u,v)][end]
    if t_s == t_enew
        pop!(E[(u,v)])
    elseif t_e >= t_enew
        E[(u,v)][end] = (t_s,t_enew)
    else
        error("algo error")
    end
end

"""
 Scale-free gene duplication

- `p`
- `q`
- `arrival` : node arrival function (:quad, :linear, :exp, :constant)

 Vazquez, Alexei and Flammini, Alessandro and Maritan, Amos and Vespignani, Alessandro 2003 Modeling of protein interaction networks Complexus 1 38–44
"""
immutable SFGD <: GraphGenerator
    p::Float64
    q::Float64
    arrival::Symbol
end

function rand(rng::AbstractRNG, gf::SFGD, Ntot::Integer, Tmax::Number, N0::Integer=5)
    nodetimes = nodearrival(gf.arrival, Ntot, N0, Tmax-1)
    # initialize nodes
    Q = init_nodearrival_pq(nodetimes)
    nnodes = length(Q)
    G = Graph(nnodes)
    E = Dict{Tuple{Int,Int},Vector{Tuple{Float64,Float64}}}() # value vector will be max size 1 always
    nodes = Int[]
    # main loop
    while !isempty(Q)
        ((u,isalive,birthtime,lifetime),curtime) = dequeue_pair!(Q)
        push!(nodes, u)
        if curtime == 0
            for v in nodes[1:end-1]
                add_edge!(G, u, v)
                add_event!(E, u, v, curtime, Tmax)
            end
        else
            v = rand(rng, nodes[1:end-1]) # select a node randomly and make u its duplicate
            W = copy(neighbors(G,v))
            for w in W
                add_edge!(G, u, w)
                add_event!(E, u, w, curtime, Tmax)
            end
            for w in W
                if rand(rng) < 0.5
                    if rand(rng) < gf.q
                        rem_edge!(G, u, w)
                        stop_event!(E, u, w, curtime)
                    end
                else
                    if rand(rng) < gf.q
                        rem_edge!(G, v, w)
                        stop_event!(E, v, w, curtime)
                    end
                end
            end
            if rand(rng) < gf.p
                add_edge!(G, u, v)
                add_event!(E, u, v, curtime, Tmax)
            end
        end
    end
    I = Int[]
    J = Int[]
    events = Events[]
    for ((u,v),evs) in E
        if !isempty(evs)
            push!(I, u)
            push!(J, v)
            push!(events, Events(evs))
        end
    end
    sparse(vcat(I,J), vcat(J,I), vcat(events,events), nnodes, nnodes, fixevents)
end

function random_random_neighbor(rng, G,u)
    v = rand(rng, neighbors(G,u))
    W = neighbors(G,v)
    w = rand(rng, W)
end

""" Social network evolution model
- `λ` : node active lifetime
- `α, β` : how active a node is at adding edges

Leskovec, J., Backstrom, L., Kumar, R., and Tomkins, A. (2008). Microscopic evolution of social networks. In Proc. of the 14th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD'08, pages 462–470.
"""
immutable SocialNE <: GraphGenerator
    λ::Float64
    α::Float64
    β::Float64
    arrival::Symbol # node arrival function, see nodearrival
end

function rand(rng::AbstractRNG, gf::SocialNE, Ntot::Integer, Tmax::Number, N0::Integer=5)
    nodetimes = nodearrival(gf.arrival, Ntot, N0, Tmax-1)
    # initialize nodes
    Q = init_nodearrival_pq(nodetimes)
    nnodes = length(Q)
    G = Graph(nnodes)
    I = Int[]
    J = Int[]
    events = Events[]
    explaw = Exponential(1/gf.λ)
    # main loop
    iter = 1
    while !isempty(Q)
        ((u,isalive,birthtime,lifetime),curtime) = dequeue_pair!(Q)
        # distinguish between birth and wake
        if !isalive # if the node is born now
            isalive = true
            lifetime = rand(rng, explaw) # node u's lifetime
            if iter != 1
                # add u's first edge to node v w/ prob. prop. to v's degree
                v = sample(rng, 1:nnodes, StatsBase.Weights(LightGraphs.degree(G)))
                add_edge!(G, u, v)
                push!(I, u)
                push!(J, v)
                push!(events, Events([(curtime,Tmax)]))
            end
        else # if the node is waking up
            # do random-random triangle closing
            w = random_random_neighbor(rng, G, u)
            if v != w
                add_edge!(G, v, w)
                push!(I, v)
                push!(J, w)
                push!(events, Events([(curtime,Tmax)]))
            end
        end
        d = LightGraphs.degree(G,u)
        # set u's sleep time
        powerlawexpcutoff = Gamma(1 - gf.α, 1/(gf.β * d))
        sleeptime = rand(powerlawexpcutoff)
        waketime = curtime + sleeptime
        # add next waketime to queue
        if (waketime <= birthtime + lifetime) && (waketime < Tmax)
            # add to queue only if node active before lifetime
            enqueue!(Q, (u,isalive,birthtime,lifetime), waketime)
        end
        iter += 1
    end
    sparse(vcat(I,J), vcat(J,I), vcat(events,events), nnodes, nnodes, fixevents)
end


text:
using Test
using CircStats

@testset "All tests" begin
    include("parameters.jl")
    include("hypothesistests.jl")
end


text:
using PyCall

# With conda-forge, the first install sometimes fails with a strange numpy
# problem.
# const asdf = pyimport_conda("asdf", "asdf", "conda-forge")
const asdf = pyimport_conda("asdf", "asdf", "astropy")

version = VersionNumber(asdf.__version__)
@info "Using Python asdf library version $version"
if version < v"2.0.3"
    @error "This version is unsupported"
    exit(1)
end


text:
using GLMakie, TestImages
coffee = testimage("coffee.png")
img = coffee[:, 80:480]
set_theme!(theme_dark())
fig = Figure(resolution = (1200,1200))
ax1 = Axis(fig[1,1])
ax2 = Axis3(fig[1,2]; aspect = (1,1,1), elevation = pi/6,
        perspectiveness = 0.5)
ax3 = Axis3(fig[2,1]; aspect = (1,1,1),
    azimuth = 4.726, elevation = 0.05, perspectiveness = 0.5)
ax4 = Axis3(fig[2,2]; aspect = (1,1,1), elevation = pi/6,
    perspectiveness = 0.5)
image!(ax1, rotr90(img))
image!(ax2, rotr90(img))
mesh!(ax3, Sphere(Point3f(0), 1), color = img)
image!(ax4, rotr90(img); transformation = (:xy, -1))
image!(ax4, rotr90(img); transformation = (:xz, 420))
image!(ax4, rotr90(img); transformation = (:yz, 420))
zlims!(ax4, -1,420)
axs = [ax1, ax2, ax3, ax4]
[hidedecorations!(ax; grid = false) for ax in axs]
colgap!(fig.layout, 0)
rowgap!(fig.layout, 0)
save("image.png", fig)
fig


text:
# Load modules
using Optim

##
function get_reflectance(ftype, λ, incident, emergent)
    # Raw measured spectrum stored in Utils
    Rexp = SpectraDB.sl2_exp_spectrum(beam.λ)
    # Reference measured spectrum stored in Utils
    Rref = SpectraDB.sl2_ref_spectrum(beam.λ)
    # Theoretical reflectance spectrum for the reference
    Rthe = theoretical_spectrum(ftype, beam, incident, emergent)
    # Calculate the absolute normalised measured spectra to fit
    Rexp_norm = normalize_reflectance(beam.λ, [beam.λ Rexp], [beam.λ Rthe], [beam.λ Rref])
end
##

# Wavelength range [nm]
λ = 250:900
# Angle of incidence [degrees]
θ = [5.]
# Polarisation (1.0 = p, 0.0 = s, between 0.0 and 1.0 = average)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Refractive indices of incident (0) and substrate (2)
incident = RIdb.air(beam.λ)
emergent = RIdb.silicon(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(incident), # 1
    ModelFit(:looyenga; N=(ninc=incident, nhost=emergent)), # 2
    LayerTMMO(emergent), # 3
]

# Get the spectrum to fit
Rexp = get_reflectance(Reflectance(), beam.λ, incident, emergent)

## Seed for the optimization algorithm

# Boundaries for the space solution and number of grid points
b = BoundariesFit(7000.0, 7300.0, 0.5, 0.65)

# Brute force search
sol = space_solution_ema(Reflectance(Rexp), b, beam, layers)

## Optimization using Optim

# Take the seed as the output from the 2D search
seed = vcat(sol.optThickness, sol.optParams)

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, store_trace=true, show_trace=true,
)

solOptim = fit_tmm_optics(
    Reflectance(Rexp), [seed], beam, layers;
    alg=SAMIN(), options=options, lb=[0.5.*seed], ub=[1.5.*seed],
)

solOptim2 = fit_tmm_optics(
    Reflectance(Rexp), [seed], beam, layers;
    alg=NelderMead(), options=options,
)


text:
"""
	grdcut(cmd0::String="", arg1=[], kwargs...)

Produce a new outgrid file which is a subregion of ingrid. The subregion is specified with
``limits`` (the -R); the specified range must not exceed the range of ingrid (but see ``extend``).

Full option list at [`grdcut`]($(GMTdoc)grdcut.html)

Parameters
----------

- **G** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdcut(....) form.
    ($(GMTdoc)grdcut.html#g)
- $(GMT.opt_J)
- **N** | **extend** :: [Type => Str or []]

    Allow grid to be extended if new region exceeds existing boundaries. Append nodata value
    to initialize nodes outside current region [Default is NaN].
    ($(GMTdoc)grdcut.html#n)
- $(GMT.opt_R)
- **S** | **circ_subregion** :: [Type => Str]    ``Arg = [n]lon/lat/radius[unit]``

    Specify an origin and radius; append a distance unit and we determine the corresponding
    rectangular region so that all grid nodes on or inside the circle are contained in the subset.
    ($(GMTdoc)grdcut.html#s)
- $(GMT.opt_V)
- **Z** | **z_subregion** :: [Type => Str]       ``Arg = [n|N |r][min/max]``

    Determine a new rectangular region so that all nodes outside this region are also outside
    the given z-range.
    ($(GMTdoc)grdcut.html#z)
- $(GMT.opt_f)
"""
function grdcut(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("grdcut", cmd0, arg1)

    d = KW(kwargs)
    symbs = (GMTver >= 6) ? [:R :J :V_params :f] : [:R :V_params :f]
	cmd = parse_common_opts(d, "", symbs)
	cmd = parse_these_opts(cmd, d, [[:G :outgrid], [:N :extend], [:S :circ_subregion], [:Z :z_subregion]])

	common_grd(d, cmd0, cmd, "grdcut ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdcut(arg1, cmd0::String=""; kw...) = grdcut(cmd0, arg1; kw...)


text:
facts("virtual_cat is like cat") do
    context("empty parameters") do
        @fact_throws virtual_cat()
        @fact_throws cat()
    end
    context("concatenating no arrays") do
        expected = cat(1)
        test = virtual_cat(1)

        @fact test --> expected
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact eachindex(test) --> eachindex(expected)
        @fact typeof(test).parameters[1] --> typeof(expected).parameters[1]
        @fact typeof(test).parameters[2] --> typeof(expected).parameters[2]
    end
    context("concatenating one 1 d array") do
        len = rand(1:10)

        a = rand(len)
        expected = cat(1,a)
        test = virtual_cat(1,a)

        @fact test.parents[1] --> a
        @fact test --> expected
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact eachindex(test) --> eachindex(expected)
        @fact typeof(test).parameters[1] --> typeof(expected).parameters[1]
        @fact typeof(test).parameters[2] --> typeof(expected).parameters[2]
    end
    context("concatenating two 1 d arrays") do
        len = rand(1:10)

        a = rand(len)
        b = rand(len)
        expected = cat(1,a,b)
        test = virtual_cat(1,a,b)

        @fact test.parents[1] --> a
        @fact test.parents[2] --> b
        @fact test --> expected
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact eachindex(test) --> eachindex(expected)
        @fact typeof(test).parameters[1] --> typeof(expected).parameters[1]
        @fact typeof(test).parameters[2] --> typeof(expected).parameters[2]
    end
end

facts("virtual_vcat is like vcat") do
    context("concatenating two 1 d arrays") do
        len = rand(1:10)

        a = rand(len)
        b = rand(len)
        expected = vcat(a,b)
        test = virtual_vcat(a,b)

        @fact test.parents[1] --> a
        @fact test.parents[2] --> b
        @fact test --> expected
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact eachindex(test) --> eachindex(expected)
        @fact typeof(test).parameters[1] --> typeof(expected).parameters[1]
        @fact typeof(test).parameters[2] --> typeof(expected).parameters[2]
    end
end

facts("virtual_hcat is like hcat") do
    context("concatenating two 1 d arrays") do
        len = rand(1:10)

        a = rand(len)
        b = rand(len)
        expected = hcat(a,b)
        test = virtual_hcat(a,b)

        @fact test.parents[1] --> a
        @fact test.parents[2] --> b
        @fact test --> expected
        @fact length(test) --> length(expected)
        @fact size(test) --> size(expected)
        @fact eachindex(test) --> eachindex(expected)
        @fact typeof(test).parameters[1] --> typeof(expected).parameters[1]
        @fact typeof(test).parameters[2] --> typeof(expected).parameters[2]
    end
end


text:
using OceanGrids, Test, Unitful

@testset "Small grid tests" begin
    nlat, nlon, ndepth = 2, 3, 4
    grid = OceanGrid(nlat, nlon, ndepth)
    @test grid isa OceanRectilinearGrid
    @test grid isa OceanGrid
    @test size(grid) == (nlat, nlon, ndepth)
    @test length(grid) == *(nlat, nlon, ndepth)
    @testset "OceanGridBox" for box in grid
        @test box isa OceanGridBox
        show(box)
        @test OceanGrids.area(box) isa Quantity
        @test unit(OceanGrids.area(box)) == u"km^2"
        @test unit(OceanGrids.volume(box)) == u"m^3"
    end
end

@testset "edges to grid" begin
    elat = -90:90:90
    elats = [-90:90:90,
             -90.0:90.0:90.0,
             [-90,0,90],
             [-90.0,0.0,90.0],
             (-90.0:90.0:90.0) * u"°",
             [-90.0,0.0,90.0] * u"°",
             (-90:90:90) * u"°",
             [-90,0,90] * u"°",]
    elons = [-180:90:180,
             [-180,-90,0,90,180],
             (-180:90:180) * u"°",
             [-180,-90,0,90,180] * u"°"]
    edepths = [[0, 100, 200, 500, 1000, 3000] * u"m",
               [0.0, 100.0, 200.0, 500.0, 1000.0, 3000.0] * u"m",
               [0.0, 100.0, 200.0, 500.0, 1000.0, 3000.0],
               [0, 100, 200, 500, 1000, 3000],
               [0, 0.1, 0.2, 0.5, 1, 3] * u"km"]
    base_grid = OceanGrid(elats[1], elons[1], edepths[1])
    @testset "types of edges" for elat in elats, elon in elons, edepth in edepths
        grid = OceanGrid(elat,elon,edepth)
        @test grid isa OceanRectilinearGrid
        @test grid isa OceanGrid
        @test grid == base_grid
        show(grid)
        @testset "OceanGridBox" for box in grid
            @test box isa OceanGridBox
        end
    end
end


text:
using Compat

myname = splitdir(@__FILE__)[end]

excluded = []

if is_windows() && Sys.WORD_SIZE == 32
  push!(excluded, "except.jl")
end

for fname in readdir()
  if fname != myname && endswith(fname, ".jl") && fname ∉ excluded
    println("running test ", fname, "...")
    include(fname)
  end
end


text: ## This READMe.jl is for part 2 of exercise 1a.


text:
using SExpressions
using EnglishText

function rem(s::String)
    Remarkable.Remark.tohtml(SExpressions.parseall(s))
end

@testset "Remark" begin

@testset "Just HTML" begin
@test sx"""
(html ([lang "en"])
  (head (title "Hello World!"))
  (body (p "This is my first Remark page")))
""" == SExpression((
        :html, ((:lang, "en"),),
        (:head, (:title, "Hello World!")),
        (:body, (:p, "This is my first Remark page"))))

@test rem("""
(html ([lang "en"])
  (head (title "Page"))
  (body
    (p "This is a poem" (br) "Line 2"
       (button ([disabled]) "disabled button"))))
""") == """
<!DOCTYPE html>
<title>Page</title>
This is a poem
Line 2<button disabled>disabled button</button>
"""
end

@test rem("""
(remark
  (define (foo-bar x y) (string (+ x y))))
(html ([lang "en"])
  (head (title "Page " (remark (foo-bar 1 1))))
  (body (p "This is page " (remark (foo-bar 1 1)) ".")))
""") == """
<!DOCTYPE html>
<title>Page 2</title>
This is page 2.
"""

@test rem("""
(remark
  (define (sqr x) (* x x))
  (define (n^4 x) (* (sqr x) (sqr x)))
  (define (test x) (string (n^4 x))))
(html ([lang "en"])
  (title (remark (test 10)))
  (p "test"))
""") == """
<!DOCTYPE html>
<title>10000</title>
test
"""

@testset "When" begin
@test rem("""
(remark (when (defined? x)
               `(p "yes 1")))
(remark (define (x y) y))
(remark (when (defined? x)
              `(p "yes 2")))
""") == """
<!DOCTYPE html>

yes 2
"""

@test rem("""
(p (remark (when (< 1 2) (define x 1) (+ x x))))
""") == """
<!DOCTYPE html>

2
"""

@test rem("""
(p (remark (when (> 1 2) "Hello")))
""") == """
<!DOCTYPE html>

"""
end

@testset "Remark" begin
@test rem("""
(remark
  (define x "Hello, World"))
(p (remark (string x "!")))
""") == """
<!DOCTYPE html>

Hello, World!
"""

@test rem("""
(remark
  (define x "Hello, ")
  (define y "World"))
(p (remark (string x y "!")))
""") == """
<!DOCTYPE html>

Hello, World!
"""

@test rem("""
(remark (define (foo x) 0))
(p (remark (foo 1)))
""") == """
<!DOCTYPE html>

0
"""

sprint() do io
    Remarkable.Remark.tohtml(io, SExpressions.parseall("""
    (p (remark (ItemList "x" "y")))
    """); modules=[EnglishText])
end == """
<!DOCTYPE html>

x and y
"""
end

@testset "Remarks" begin
@test rem("""
(remarks
 (define n 7000000000)
 `((p "Hello World")
   (p "All " ,n " are welcome!")))
""") == """
<!DOCTYPE html>

Hello World

All 7000000000 are welcome!
"""
end

@test Remarkable.Remark.tohtml("data/file1.rem") == """
<!DOCTYPE html>

File 2: 100

File 1: 20
"""

@test Remarkable.Remark.tohtml("data/test-dispatch.rem") == """
<!DOCTYPE html>

12
"""

@test Remarkable.Remark.tohtml("data/test-markdown.rem") == """
<!DOCTYPE html>

Some Markdown

Test.
"""

@testset "Each" begin
@test rem("""
(#:each x (List "x" "y" "z")
  `((p ,x)))
""") == """
<!DOCTYPE html>

x

y

z
"""
end

@testset "Files" begin
@test Remarkable.Remark.tohtml("data/literal-test.rem") == """
<!DOCTYPE html>
"""
end

@testset "Markdown Render" begin
@test Remarkable.Remark.tohtml("data/markdown-render.rem") == """
<!DOCTYPE html>

my link to
"""
end

@testset "Object Include" begin
@test Remarkable.Remark.tohtml("data/object-include.rem") == """
<!DOCTYPE html>

Hello, World!
"""
end

@testset "Attributes" begin
    @test rem("""
(html ([lang "en"])
  (body (textarea ([cols 10]))))
""") == """
<!DOCTYPE html>
<textarea cols="10"></textarea>"""
end

end


text:
import HTTP, JSON3
using Test
using DashBase
using Dash

@testset "Components" begin

    a_comp = html_a("test", id = "test-a")

    @test a_comp.id == "test-a"
    @test a_comp.children == "test"

    input_comp = dcc_input(id = "test-input", type="text")
    @test input_comp.id == "test-input"
    @test input_comp.type == "text"
    @test_throws ErrorException html_a(undefined_prop = "rrrr")

    component_with_children = html_div() do
        html_a("fffff"),
        html_h1("fffff")
    end

    @test !isnothing(component_with_children.children)
    @test component_with_children.children isa Tuple{Component, Component}
    @test length(component_with_children.children) == 2
    @test DashBase.get_type(component_with_children.children[1]) == "A"
    @test DashBase.get_type(component_with_children.children[2]) == "H1"

end

@testset "handler" begin
    app = dash(external_stylesheets=["test.css"])
    app.layout = html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            html_div(id = "my-div"),
            html_div(id = "my-div2")
        end
    callback!(app, Output("my-div","children"), Input("my-id","value")) do value
        return value
    end
    callback!(app, Output("my-div2","children"), Input("my-id","value")) do value
        return "v_$(value)"
    end

    handler = Dash.make_handler(app)

    request = HTTP.Request("GET", "/")
    response = HTTP.handle(handler, request)
    @test response.status == 200
    body_str = String(response.body)

    request = HTTP.Request("GET", "/_dash-layout")
    response = HTTP.handle(handler, request)
    @test response.status == 200
    body_str = String(response.body)
    @test body_str == JSON3.write(app.layout)

    request = HTTP.Request("GET", "/_dash-dependencies")
    response = HTTP.handle(handler, request)
    @test response.status == 200
    body_str = String(response.body)
    resp_json = JSON3.read(body_str)

    @test resp_json isa AbstractVector
    @test length(resp_json) == 2
    @test haskey(resp_json[1], :inputs)
    @test haskey(resp_json[1], :state)
    @test haskey(resp_json[1], :output)
    @test length(resp_json[1].inputs) == 1
    @test length(resp_json[1].state) == 0
    @test resp_json[1].inputs[1].id == "my-id"
    @test resp_json[1].inputs[1].property == "value"
    @test resp_json[1].output == "my-div.children"

    test_json = """{"output":"my-div.children","changedPropIds":["my-id.value"],"inputs":[{"id":"my-id","property":"value","value":"initial value3333"}]}"""
    request = HTTP.Request("POST", "/_dash-update-component", [], Vector{UInt8}(test_json))
    response = HTTP.handle(handler, request)
    @test response.status == 200
    body_str = String(response.body)

    @test body_str == """{"response":{"my-div":{"children":"initial value3333"}},"multi":true}"""
end

@testset "layout as function" begin
    app = dash()
    global_id = "my_div2"
    layout_func = () -> begin
        html_div() do
            dcc_input(id = "my-id", value="initial value", type = "text"),
            html_div(id = "my-div"),
            html_div(id = global_id)
        end
    end
    app.layout = layout_func
    handler = Dash.make_handler(app)
    request = HTTP.Request("GET", "/_dash-layout")

    response = HTTP.handle(handler, request)
    @test response.status == 200
    body_str = String(response.body)
    @test body_str == JSON3.write(layout_func())
    @test occursin("my_div2", body_str)

    global_id = "my_div3"
    response = HTTP.handle(handler, request)
    @test response.status == 200
    body_str = String(response.body)
    @test body_str == JSON3.write(layout_func())
    @test occursin("my_div3", body_str)

end

@testset "assets" begin
    app = dash(assets_folder = "assets")
    app.layout = html_div() do
            html_img(src = "assets/test.png")
        end
    @test Dash.get_assets_path(app) == joinpath(pwd(),"assets")
    handler = Dash.make_handler(app)
    request = HTTP.Request("GET", "/assets/test.png")
    response = HTTP.handle(handler, request)
    @test response.status == 200
    body_str = String(response.body)

    request = HTTP.Request("GET", "/assets/wrong.png")
    response = HTTP.handle(handler, request)
    @test response.status == 404
    body_str = String(response.body)

end

@testset "PreventUpdate and no_update" begin
    app = dash()

    app.layout = html_div() do
            html_div(10, id = "my-id"),
            html_div(id = "my-div")
        end
    callback!(app,  Output("my-div","children"), Input("my-id","children")) do value
        throw(PreventUpdate())
    end

    handler = Dash.make_handler(app)

    test_json = """{"output":"my-div.children","changedPropIds":["my-id.children"],"inputs":[{"id":"my-id","property":"children","value":10}]}"""

    request = HTTP.Request("POST", "/_dash-update-component", [], Vector{UInt8}(test_json))
    response = HTTP.handle(handler, request)
    @test response.status == 204
    @test length(response.body) == 0

    app = dash()
    app.layout = html_div() do
            html_div(10, id = "my-id"),
            html_div(id = "my-div"),
            html_div(id = "my-div2")
        end
    callback!(app, [Output("my-div","children"), Output("my-div2","children")], Input("my-id","children")) do value
        no_update(), "test"
    end

    handler = Dash.make_handler(app)
    test_json = """{"output":"..my-div.children...my-div2.children..","changedPropIds":["my-id.children"],"inputs":[{"id":"my-id","property":"children","value":10}]}"""

    request = HTTP.Request("POST", "/_dash-update-component", [], Vector{UInt8}(test_json))
    response = HTTP.handle(handler, request)

    @test response.status == 200
    result = JSON3.read(String(response.body))

    @test length(result[:response]) == 1
    @test haskey(result[:response], Symbol("my-div2"))
    @test !haskey(result[:response], Symbol("my-div"))
    @test result[:response][Symbol("my-div2")][:children] == "test"
end

@testset "layout validation" begin
    app = dash(assets_folder = "assets_compressed", compress = true)
    @test_throws ErrorException make_handler(app)
    app.layout = html_div(id="top") do
        html_div(id="second") do
            html_div("dsfsd", id = "inner1")
        end,
        html_div(id = "third") do
            html_div("dsfsd", id = "inner2")
        end
    end
    make_handler(app)

    app.layout = html_div(id="top") do
        html_div(id="second") do
            html_div("dsfsd", id = "inner1")
        end,
        html_div(id = "second") do
            html_div("dsfsd", id = "inner2")
        end
    end
    @test_throws ErrorException make_handler(app)
end


text:
function compute_boundary_t!(boundary_t::Boundary, bulkevol::BulkEvolved,
                             boundary::Boundary, gauge::Gauge, sys::System, ::EvolTest0)

    a4_t, fx2_t, fy2_t = unpack(boundary_t)
    # a4  , fx2  , fy2   = unpack(boundary)

    fill!(a4_t,  0)
    fill!(fx2_t, 0)
    fill!(fy2_t, 0)

    nothing
end

function compute_boundary_t!(boundary_t::Boundary, bulk::BulkEvolved,
                             boundary::Boundary, gauge::Gauge, sys::System{Inner},
                             evoleq::AffineNull)
    Du  = sys.Du
    Dx  = sys.Dx
    Dy  = sys.Dy

    phi0  = evoleq.phi0
    phi03 = phi0 * phi0 * phi0

    _, Nx, Ny = size(sys)

    a4_t, fx2_t, fy2_t = unpack(boundary_t)

    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            xi      = gauge.xi[1,i,j]
            xi3     = xi*xi*xi

            phi     = bulk.phi[1,i,j]

            xi_x    = Dx(gauge.xi, 1,i,j)
            xi_y    = Dy(gauge.xi, 1,i,j)

            phi_u   = Du(bulk.phi, 1,i,j)
            phi_x   = Dx(bulk.phi, 1,i,j)
            phi_y   = Dy(bulk.phi, 1,i,j)

            b14_x   = Dx(bulk.B1, 1,i,j)
            b14_y   = Dy(bulk.B1, 1,i,j)

            b24_x   = Dx(bulk.B2, 1,i,j)
            b24_y   = Dy(bulk.B2, 1,i,j)

            g4_x    = Dx(bulk.G, 1,i,j)
            g4_y    = Dy(bulk.G, 1,i,j)

            a4_x    = Dx(boundary.a4, 1,i,j)
            a4_y    = Dy(boundary.a4, 1,i,j)

            fx2_x   = Dx(boundary.fx2, 1,i,j)
            fy2_y   = Dy(boundary.fy2, 1,i,j)

            # phi2 = phi0^3 phi(u=0) - phi0 xi^2
            phi2    = phi03 * phi - phi0 * xi * xi
            phi2_x  = phi03 * phi_x - 2 * phi0 * xi * xi_x
            phi2_y  = phi03 * phi_y - 2 * phi0 * xi * xi_y

            # phi2_t = 3 xi phi2 + phi0 xi^3 + phi0^3 phi(u=0)_u
            phi2_t  = 3 * xi * phi2 + phi0 * xi3 + phi03 * phi_u

            a4_t[1,i,j]  = -4//3 * (fx2_x + fy2_y + phi0 * phi2_t)

            fx2_t[1,i,j] = g4_y - b14_x - b24_x - 0.25 * a4_x + 1//3 * phi0 * phi2_x
            fy2_t[1,i,j] = g4_x + b14_y - b24_y - 0.25 * a4_y + 1//3 * phi0 * phi2_y
        end
    end

    nothing
end


text:
#!/usr/bin/env julia
using CSV, DataFrames, PGFPlotsX
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{fontspec}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\setmainfont{Asana-Math}")
csvform = tempname()
run(pipeline(`./yyg_se_0804.sh`, stdout=csvform))
yyg_se_0804 = CSV.File(csvform) |> DataFrame
pds = [
    Dict(:c => :predicted_deaths_mean, :lab => "Medel",),
    Dict(:c => :predicted_deaths_lower, :lab => "2,5 percentil",),
    Dict(:c => :predicted_deaths_upper, :lab => "97,5 percentil",)]
p = @pgf Axis({title = "Döda Sverige 2020-08-04 YYG",
    date_coordinates_in = "x",
    x_tick_label_style = "{rotate=90}",
    y_tick_label_style = {"/pgf/number format/use comma"},
    xlabel = "Prediktionsdatum", ylabel = "Kumulativt antal döda",
    xmajorgrids, ymajorgrids})
for pd in pds
    @pgf push!(p, PlotInc({no_markers}, Table([yyg_se_0804[:date], yyg_se_0804[pd[:c]]])),
        LegendEntry(pd[:lab]))
end
pgfsave("../../images/yyg_se_0804.svg", p)

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
@testset verbose=true "Interpolation" begin
    @testset verbose=true "Unflip" begin
        Random.seed!(1234)
        function measure_flips(q)
            flips = 0
            for i in 2:size(q, 1)
                if q[i].w*q[i-1].w+q[i].x*q[i-1].x+q[i].y*q[i-1].y+q[i].z*q[i-1].z < 0
                    flips += 1
                end
            end
            flips
        end

        @testset verbose=true "$T" for T in FloatTypes
            q = randn(Quaternion{T}, 100_000);
            flips = measure_flips(q)
            @test flips > 0
            p = unflip(q)
            flips = measure_flips(q)
            @test flips > 0
            flips = measure_flips(p)
            @test flips == 0
            unflip!(q)
            flips = measure_flips(q)
            @test flips == 0

            q = randn(Quaternion{T}, 10_000, 3);
            flips = sum(measure_flips(q[:, i]) for i in 1:3)
            @test flips > 0
            p = unflip(q)
            flips = sum(measure_flips(q[:, i]) for i in 1:3)
            @test flips > 0
            flips = sum(measure_flips(p[:, i]) for i in 1:3)
            @test flips == 0
            unflip!(q)
            flips = sum(measure_flips(q[:, i]) for i in 1:3)
            @test flips == 0

            q = randn(Quaternion{T}, 3, 1_000, 3);
            flips = sum(measure_flips(q[i, :, j]) for i in 1:3 for j in 1:3)
            @test flips > 0
            p = unflip(q, dim=2)
            flips = sum(measure_flips(q[i, :, j]) for i in 1:3 for j in 1:3)
            @test flips > 0
            flips = sum(measure_flips(p[i, :, j]) for i in 1:3 for j in 1:3)
            @test flips == 0
            unflip!(q, dim=2)
            flips = sum(measure_flips(q[i, :, j]) for i in 1:3 for j in 1:3)
            @test flips == 0
        end
    end

    @testset verbose=true "Slerp" begin
        ϵ = 20eps()
        comp = [0, 1, -1]
        Rs = [
            [Rotor(w, x, y, z) for w in comp for x in comp for y in comp for z in comp][2:end]...;
            randn(RotorF64, 20)...
        ]
        directions = [Rotor(1), Rotor(imx), Rotor(imy), Rotor(imz), Rotor(-imx), Rotor(-imy), Rotor(-imz)]

        # Check extremes
        for Q1 in directions
            @test distance(slerp(Q1, Q1, 0.0), Q1) < ϵ
            @test distance(slerp(Q1, Q1, 1.0), Q1) < ϵ
            @test distance(slerp(Q1, -Q1, 0.0), Q1) < ϵ
            @test distance(slerp(Q1, -Q1, 1.0), Q1) < ϵ
            for Q2 in directions
                @test distance(slerp(Q1, Q2, 0.0), Q1) < ϵ
                @test distance(slerp(Q1, Q2, 1.0), Q2) < ϵ
                @test distance(slerp(Q1, -Q2, 0.0), Q1) < ϵ
                @test distance(slerp(Q1, -Q2, 1.0), -Q2) < ϵ
                @test distance(slerp(Q2, Q1, 0.0), Q2) < ϵ
                @test distance(slerp(Q2, Q1, 1.0), Q1) < ϵ
            end
        end

        # Test simple rotations about each axis
        for Q2 in directions[2:end]
            for t in LinRange(0.0, 1.0, 100)
                @test distance(
                    slerp(Rotor(1), Q2, t),
                    Rotor(cos(π * t / 2) + sin(π * t / 2) * Q2)
                ) < ϵ
            end
        end

        # Test that (slerp of rotated rotors) is (rotated slerp of rotors)
        for R in Rs
            for Q2 in directions[1:end]
                for t in LinRange(0.0, 1.0, 100)
                    @test distance(
                        R * slerp(Rotor(1), Q2, t),
                        slerp(R, R*Q2, t)
                    ) < ϵ
                end
            end
        end

        # Check that unflipping works
        for Q in directions
            for t in LinRange(0.0, 1.0, 10)
                @test distance(slerp(Q, -Q, t, unflip=true), Q) < ϵ
            end
        end

    end

    @testset verbose=true "Squad" begin
        ϵ = 20eps()
        comp = [0, 1, -1]
        Rs = [
            [Rotor(w, x, y, z) for w in comp for x in comp for y in comp for z in comp][2:end]...;
            randn(RotorF64, 20)...
        ]
        directions = [Rotor(1), Rotor(imx), Rotor(imy), Rotor(imz), Rotor(-imx), Rotor(-imy), Rotor(-imz)]

        ω = 0.1
        tin = collect(LinRange(-10, 10, 201))
        Rin = [exp(ω*ti*imz/2) for ti in tin]
        tout = (tin[1:end-1] + tin[2:end]) / 2
        Rout = squad(Rin, tin, tout, validate=true)
        Ravg = [exp(ω*ti*imz/2) for ti in tout]
        @test maximum(distance.(Rout, Ravg)) < 2eps()
        random_signs = [1; rand([-1, 1], length(Rin)-1)...]  # First one must be 1
        Rin_flipped = [Rotor(r*R) for (r, R) in zip(random_signs, Rin)]
        @test Rout == squad(Rin_flipped, tin, tout, validate=true, unflip=true)
        Rout = squad(Rin, tin, tin, validate=true)
        @test maximum(distance.(Rout, Rin)) < 2eps()
        @test distance(Rotor(1), squad(Rin, tin, 0, validate=true)) == 0
        @test distance(Rin[1], squad(Rin, tin, tin[1], validate=true)) == 0
        @test distance(Rin[end], squad(Rin, tin, tin[end], validate=true)) == 0
        @test squad(Rin, tin, Vector{eltype(tin)}()) == Vector{eltype(Rin)}()

        t_in = LinRange(0.0, 1.0, 13)
        t_out = LinRange(0.0, 1.0, 37)
        t_out2 = sort(rand(59))
        # squad interpolated onto the inputs should be the identity
        for R1 in Rs
            for R2 in Rs
                R_in = [slerp(R1, R2, t) for t in t_in]
                @test squad(R_in, t_in, t_in) ≈ R_in atol=ϵ
            end
        end

        # squad should be the same as slerp for linear interpolation
        for R in directions
            R_in = [slerp(Rotor(1), R, t) for t in t_in]
            R_out_squad = squad(R_in, t_in, t_out)
            R_out_slerp = [slerp(Rotor(1), R, t) for t in t_out]
            @test R_out_squad ≈ R_out_slerp atol=ϵ
            R_out_squad = squad(R_in, t_in, t_out2)
            R_out_slerp = [slerp(Rotor(1), R, t) for t in t_out2]
            @test R_out_squad ≈ R_out_slerp atol=ϵ
        end
    end

    @testset verbose=true "Squad derivatives" begin
        R, ω⃗, Ṙ = precessing_nutating_example()
        tin = collect(range(0, 1, length=1000))
        for tout in [tin, tin[498]]
            Rsquad, ω⃗squad, Ṙsquad = squad(
                Rotor.(R.(tin)), tin, tout,
                compute_angular_velocity=true, compute_derivative=true
            )
            ω⃗eval = ω⃗.(tout)
            Ṙeval = Ṙ.(tout)
            @test ω⃗squad ≈ ω⃗eval rtol=1e-6
            @test Ṙsquad ≈ Ṙeval rtol=1e-6
            Rsquad2, ω⃗squad2 = squad(
                Rotor.(R.(tin)), tin, tout,
                compute_angular_velocity=true
            )
            @test Rsquad == Rsquad2
            @test ω⃗squad == ω⃗squad2
            Rsquad3, Ṙsquad3 = squad(
                Rotor.(R.(tin)), tin, tout,
                compute_derivative=true
            )
            @test Rsquad == Rsquad3
            @test Ṙsquad == Ṙsquad3
        end
    end
end