text:
#Abstract volume description
abstract type Volume end

"""Describing a volume with cubic voxels"""
type CubeVolume <: Volume
    cube::Array{Complex{Float64},3}
    cubesize::Int64
    rmax::Float64

    function CubeVolume(cube::Array{Complex{Float64},3}, cubesize::Int64, rmax::Float64)
        new(cube, cubesize, rmax)
    end
    function CubeVolume(cube::Array{Float64,3}, cubesize::Int64, rmax::Float64)
        new(complex(cube), cubesize, rmax)
    end
    function CubeVolume(cubesize::Int64, rmax::Float64)
        new(zeros(Complex{Float64}, cubesize, cubesize, cubesize), cubesize, rmax)
    end
end

"""Get voxel size of Cube Volume"""
function dr(cube::CubeVolume) return 2*cube.rmax/(cube.cubesize-1) end
function qmax(cube::CubeVolume) return pi/dr(cube) end

function +(a::CubeVolume, b::CubeVolume) return CubeVolume(a.cube + b.cube, a.cubesize, a.rmax) end
function -(a::CubeVolume, b::CubeVolume) return CubeVolume(a.cube - b.cube, a.cubesize, a.rmax) end
function *(a::CubeVolume, b::Number) return CubeVolume(a.cube * b, a.cubesize, a.rmax) end
function *(a::Number, b::CubeVolume) return CubeVolume(b.cube * a, b.cubesize, b.rmax) end
function /(a::CubeVolume, b::Number) return CubeVolume(a.cube / b, a.cubesize, a.rmax) end
function sum(abs,a::CubeVolume) return sum(abs, a.cube) end
function abs(a::CubeVolume) return CubeVolume(abs(a.cube), a.cubesize, a.rmax) end
function real(a::CubeVolume) return Cubevolume(real(a.cube), a.cubesize, a.rmax) end
function maximum(a::CubeVolume) return Base.maximum(real(a.cube)) end

"""Calculates the density cube from atomlist coordinates"""
function getDensityCube(atomlist, cubesize::Int64, rmax::Float64)
    r = linspace(-rmax, rmax, cubesize)
    cube = Complex{Float64}[getDensity_cartesian(atomlist, [x,y,z]) for x=r, y=r, z=r]
    return CubeVolume(cube,cubesize, rmax)
end

"""Average a list of cubes"""
function average(list::Array{CubeVolume})
    return reduce(+, list) / length(list)
end

"""Center of mass calculation for fitting in real space."""
function calculate_center_of_mass(volume::CubeVolume)
    r = linspace(-volume.rmax, volume.rmax, volume.cubesize)
    positions = Vector{Float64}[Float64[x,y,z] for x=r, y=r, z=r ]
    cube = real(volume.cube)
    convert(Vector{Float64},sum(positions .* cube)/sum(cube))
end

"""Translates the content of a cube into a certain direction"""
function translate_cube(volume::CubeVolume, direction::Vector{Float64})
    r = linspace(-volume.rmax, volume.rmax, volume.cubesize)
    return CubeVolume(Complex{Float64}[getVolumeInterpolated(volume, Float64[x,y,z]+direction) for x=r, y=r, z=r], volume.cubesize, volume.rmax)
end

"""Fits a cube into another cube"""
function center_cube(volume::CubeVolume)
    translated_vol = translate_cube(volume, calculate_center_of_mass(volume) )
    return translated_vol
end

"""Translates the volume to the center of mass of the reference"""
function displacement_fitting(volume::CubeVolume, reference::CubeVolume)
    translated_vol = translate_cube(volume, -(calculate_center_of_mass(reference) - calculate_center_of_mass(volume)) )
    return translated_vol,similarity(translated_vol, reference)
end

# """Calculates the Intensity cube from direct expression"""
# function getIntensityCube(cubesize)
#     r = linspace(-qmax, qmax, cubesize)
#     cube = Complex{Float64}[getIntensity_cartesian([x,y,z]) for x=r, y=r, z=r]
#     return cube
# end

"""cubic forward fourier transform"""
function forward(volume::CubeVolume)
    ftcube = fftshift(fft(ifftshift(volume.cube)))
    return CubeVolume(ftcube, volume.cubesize, qmax(volume))
end

"""Calculates the aboslute square of volume"""
function absoluteSquare(volume::CubeVolume)
    return CubeVolume(abs.(volume.cube).^2, volume.cubesize, volume.rmax)
end

"""Mirrors the cube along one axis"""
function mirrorCube(volume::CubeVolume)
    cubesize = volume.cubesize
    mx,my,mz = deepcopy(volume),deepcopy(volume),deepcopy(volume)

    # mx.cube = Complex{Float64}[volume.cube[cubesize-x+1,y,z] for x=1:cubesize,y=1:cubesize,z=1:cubesize]
    # my.cube = Complex{Float64}[volume.cube[x,cubesize-y+1,z] for x=1:cubesize,y=1:cubesize,z=1:cubesize]
    mz.cube = Complex{Float64}[volume.cube[x,y,cubesize-z+1] for x=1:cubesize,y=1:cubesize,z=1:cubesize]

    # return [mx,my,mz]
    return mz
end

function similarity(volume1::CubeVolume, volume2::CubeVolume, K::Int64)
    return similarity(volume1, volume2)
end

function similarity(volume1::CubeVolume, volume2::CubeVolume, K_range::UnitRange{Int64})
    return similarity(volume1, volume2)
end

"Calculates the similarity between two cubes"
function similarity(volume1::CubeVolume, volume2::CubeVolume)
    vec1 = real(reshape(volume1.cube, volume1.cubesize^3))
    vec2 = real(reshape(volume2.cube, volume2.cubesize^3))
    vec1 = vec1 / norm(vec1)
    vec2 = vec2 / norm(vec2)
    return cor(vec1, vec2)
end

"""Saves a cubic structure to a .mrc file"""
function saveCube(volume::CubeVolume, filename::String)
    cubesize = volume.cubesize
    rmax = volume.rmax
    cube = convert(Array{Float32}, real(volume.cube))
    cube = cube / sum(abs,cube)
    offset = floor(cubesize/2)

    println("Saving cube to $filename with cubesize=$cubesize")
    out = open(filename, "w+")
    write(out, Int32[cubesize, cubesize, cubesize]) #columns, rows, sections
    write(out, Int32[2]) #mode
    write(out, Int32[-offset,-offset,-offset]) #starting point
    write(out, Int32[cubesize,cubesize,cubesize]) #grid size
    write(out, Float32[2.0*rmax, 2.0*rmax, 2.0*rmax]) #cell size
    write(out, Float32[90.0, 90.0, 90.0]) #cell angles
    write(out, Int32[1,2,3]) #mapc mapr, maps
    write(out, Float32[minimum(cube), maximum(cube),mean(cube)])
    write(out, Int32[0, 0]) #ispg next
    write(out, zeros(Int8, 1024-96))

    write(out, cube)
    close(out)
end

"""Loads a cube from a .mrc file"""
function loadCube(filename, layers=0)
    out = open(filename, "r")
    cubesizes = convert(Array{Int64},read(out, Int32, 3)) #columns, rows, sections
    mode = read(out, Int32) #mode
    startingpoint = read(out, Int32, 3) #starting point
    gridsize = read(out, Int32, 3) #grid size
    cellsize = read(out, Float32, 3) #cell size
    cellangles = read(out, Float32, 3) #cell angles
    mapparameters = read(out, Int32, 3) #mapc mapr, maps
    maxmimean = read(out, Float32, 3)
    ispgnext = read(out, Int32, 2) #ispg next
    stuff = read(out, Int8, 1024-96)
    length = layers > 0 ? layers : cubesizes[3]
    data = read(out, Float32, cubesizes[1]*cubesizes[2]*length)
    close(out)
    return CubeVolume(reshape(convert(Array{Complex{Float64}}, data),cubesizes[1], cubesizes[2], length), cubesizes[1], convert(Float64,cellsize[1]/2.0) )
end

"""Gets the value within a cube via trilinear interpolation"""
function getVolumeInterpolated(volume::CubeVolume, vec::Vector{Float64})
    cubesize = volume.cubesize
    rmax = volume.rmax
    deltar = dr(volume)
    cube = volume.cube

    x = clamp((vec[1] + rmax)/deltar + 1.0, 1.0, Float64(cubesize))
    y = clamp((vec[2] + rmax)/deltar + 1.0, 1.0, Float64(cubesize))
    z = clamp((vec[3] + rmax)/deltar + 1.0, 1.0, Float64(cubesize))

    x0,x1 = floor(Int64, x), ceil(Int64, x)
    y0,y1 = floor(Int64, y), ceil(Int64, y)
    z0,z1 = floor(Int64, z), ceil(Int64, z)

    a = Float64(x-x0)
    b = Float64(y-y0)
    c = Float64(z-z0)

    V000 = (1.0-a)*(1.0-b)*(1.0-c)
    V100 = a *(1.0-b)*(1.0-c)
    V010 = (1.0-a)* b*(1.0-c)
    V001 = (1.0-a)*(1.0-b)*c
    V101 = a*(1.0-b)*c
    V011 = (1.0-a)*b*c
    V110 = a* b*(1.0-c)
    V111 = a* b* c

    @inbounds @fastmath return V000 * cube[x0,y0,z0] + V100 * cube[x1,y0,z0] + V010 * cube[x0,y1,z0] + V001 * cube[x0,y0,z1] + V101 * cube[x1,y0,z1] + V011 * cube[x0,y1,z1] + V110 * cube[x1,y1,z0] + V111 * cube[x1,y1,z1]
end

function rotateStructure(volume::CubeVolume, theta::Float64, phi::Float64, gamma::Float64, K::Int64, lrange::Range)
    return rotateStructure(volume, theta, phi, gamma)
end

function rotateStructure(volume::CubeVolume, theta::Float64, phi::Float64, gamma::Float64, K::Int64, L::Int64)
    return rotateStructure(volume, theta, phi, gamma)
end

"""Rotates a cubic volume"""
function rotateStructure(volume::CubeVolume, theta::Float64, phi::Float64, gamma::Float64)
    newvol = deepcopy(volume)
    r = linspace(-volume.rmax, volume.rmax, volume.cubesize)
    newvol.cube = Complex{Float64}[getVolumeInterpolated(volume, euler(theta, phi, gamma)*Float64[x,y,z]) for x=r, y=r, z=r]
    return newvol
end

"""Calculates the shell-wise correlation of two cubes"""
function shell_correlation(volume1::CubeVolume, volume2::CubeVolume, K::Int64=0, cor_method=cor_FSC; N::Int64=48)
    K = (K == 0 ? floor(Int64,volume1.cubesize/2) : K)
    fsc = zeros(Float64,K)

    #let's sample the sphere evenly
    for k = 1:K
        r = k*dr(volume1)
        surfA = Complex{Float64}[ getVolumeInterpolated(volume1, sphericalToCartesian(r, phi, theta)) for phi=linspace(0,2.0*pi,N), theta=acos.(linspace(-1.0,1.0,N)) ]
        surfB = Complex{Float64}[ getVolumeInterpolated(volume2, sphericalToCartesian(r, phi, theta)) for phi=linspace(0,2.0*pi,N), theta=acos.(linspace(-1.0,1.0,N)) ]

        surfA = reshape(surfA,N^2)
        surfB = reshape(surfB,N^2)
        fsc[k] = cor_method(surfA, surfB)
    end

    return fsc
end

function shell_correlation_ISC(volume1::CubeVolume, volume2::CubeVolume, K::Int64=0)
    shell_correlation(volume1, volume2, K, cor_ISC)
end

function shell_correlation_FSC(volume1::CubeVolume, volume2::CubeVolume, K::Int64=0)
    shell_correlation(volume1, volume2, K, cor_FSC)
end

"""Calculate the Intensity Shell Correlation between two cubes"""
function ISC(density::CubeVolume, intensity::CubeVolume, K::Int64=0)
    return shell_correlation_ISC(intensity, absoluteSquare(forward(density)), K)
end

"""Calculate the Fourier Shell Correlation between two cubes"""
function FSC(density::CubeVolume, fourier::CubeVolume, K::Int64=0)
    return shell_correlation_FSC(fourier, forward(density), K)
end


text:
module TestVariables

using TinyModia
using ModiaPlot
using Measurements
using StaticArrays
using Test

@testset "... Test variable declarations" begin

v1 = input | output | flow | Var(min=0, nominal=10, init=5)
@test v1 == (class = :Var, input = true, output = true, flow = true, min = 0, nominal = 10, init = 5)

v2 = parameter | 5 | Var(info="v2 info")
@test v2 == (class = :Var, parameter = true, value = 5, info = "v2 info")

v3 = "string value" | parameter | info"v3 info"
@test v3 == (class = :Var, parameter = true, value = "string value", info = "v3 info")

v4 = parameter | 5.5u"m/s" | Var(min=0)
@test v4 == (class = :Var, parameter = true, value = 5.5u"m*s^-1", min = 0)

v5 = parameter | Var(5, min=0)
@test v5 == (class = :Var, parameter = true, value = 5, min = 0)

Pin = Model( v = potential | Var(nominal=10), i = flow )
@test Pin == (class = :Model, v = (class = :Var, potential = true, nominal = 10), i = (class = :Var, flow = true,))

v6 = Var(2.0)
@test v6 == (class = :Var, value = 2.0,)

v7 = Var(2.0 ± 0.1)
@test v7 == (class = :Var, value = 2.0 ± 0.1,)

v8 = Var(SVector(1.0, 2.0, 3.0))
@test v8 == (class = :Var, value = [1.0, 2.0, 3.0],)

v9 = v4 | Var(parameter=false)
@test v9 == (class = :Var, parameter = false, value = 5.5u"m*s^-1", min = 0)

v10 = v4 | Var(parameter=nothing)
@test v10 == (class = :Var, value = 5.5u"m*s^-1", min = 0)

v11 = v4 | interval(1,10)
@test v11 == (class = :Var, parameter = true, value = 5.5u"m*s^-1", min = 1, max = 10)

TestVar1 = Model(
    p = parameter | -1 | info"A parameter",
    s = parameter | "string value" | info"A string parameter",
    x = Var(init=0.1),
    y = Var(start=0.1), # Not used
    equations = :[
        der(x) = p*x+1
    ]
)
#@showModel TestVar1
model = @instantiateModel(TestVar1, log=false, logCode=false)
simulate!(model, merge=Map(p=-2, x=0.2), requiredFinalStates = [0.4593994150057028]) # x.init is not changed
plot(model, "x")

end

end

ext: jl


text:
d = pwd()

function make_toc(s)
    lines = split(s, '\n')
    tocline = findfirst(startswith("# TOC"), lines)
    tocline === nothing && return s

    others = lines[tocline+1:end]
    toc_entries = map(others) do line
        m = match(r"# (#+) (.+)", line)
        m === nothing && return nothing
        n = length(m[1])
        title = strip(m[2])
        githubanchor = replace(title, r"\s+" => "-")
        githubanchor = replace(githubanchor, r"[^\-a-z0-9A-Z]+" => "")
        (title, n, githubanchor)
    end

    toc_entries = filter(!isnothing, toc_entries)
    s = map(toc_entries) do (title, n, githubanchor)
        # "# $(" "^n)- [$title](#$githubanchor)"
        "# - [$title](#$githubanchor)"
    end |> x -> join(x, "\n")
    lines[tocline] = s
    join(lines, "\n")
end

try
    cd(@__DIR__)
    Literate.markdown("README.jl",
        execute = true,
        documenter = false,
        credit = false,
        preprocess = make_toc)
    mv("README.md", "../README.md", force = true)
finally
    cd(d)
end


text:
@testset "PcapStreamWriter" begin
    io = IOBuffer()
    writer = PcapStreamWriter(io)
    write(writer, UnixTime(Dates.UTInstant(Nanosecond(1234))), [0x5, 0x6, 0x7, 0x8, 0x9])
    data = take!(io)
    @assert length(data) == sizeof(PcapHeader) + sizeof(RecordHeader) + 5
    GC.@preserve data begin
        p = pointer(data)
        h = unsafe_load(Ptr{PcapHeader}(p))
        @test h.magic == 0xa1b23c4d
        @test h.version_major == 2
        @test h.version_minor == 4
        @test h.thiszone == 0
        @test h.sigfigs == 0
        @test h.linktype == LINKTYPE_ETHERNET
        rh = unsafe_load(Ptr{RecordHeader}(p + sizeof(PcapHeader)))
        @test rh.ts_sec == 0
        @test rh.ts_usec == 1234
        @test rh.incl_len == 5
        @test rh.orig_len == 5
    end
    @test data[end-4:end] == [0x5, 0x6, 0x7, 0x8, 0x9]
end


text:
# ---
# title: 203. Remove Linked List Elements
# id: problem203
# author: zhwang
# date: 2022-02-16
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/remove-linked-list-elements/description/>
# hidden: true
# ---
#
# Remove all elements from a linked list of integers that have value **_val_**.
#
# **Example:**
#
#
#
#     Input:  1->2->6->3->4->5->6, _**val**_ = 6
#     Output: 1->2->3->4->5
#
#
#
## @lc code=start
using LeetCode

function remove_elements!(
    head::Union{ListNode,Nothing}, value::Int
)::Union{ListNode,Nothing}
    cur = fake_head = ListNode()
    fake_head.next = head
    while !isnothing(cur.next)
        if cur.next.val == value
            cur.next = cur.next.next
        else
            cur = cur.next
        end
    end
    return fake_head.next
end
## @lc code=end


text:
using MySample
using Test
import Random

@testset "product-form" begin
    N=1000
    A = reshape(collect(1:1.0:N^2), N, N)
    B = reshape(collect(1:1.0:N^2), N, N)
    C = zeros(N, N)
    @time naiveform!(C, A, B)
    # copy result
    refC = copy(C)

    C = zeros(N, N)
    @time ipform!(C, A, B)
    @test C ≈ refC

    C = zeros(N, N)
    @time opform!(C, A, B)
    @test C ≈ refC

    C = zeros(N, N)
    @time mpform!(C, A, B)
    @test C ≈ refC
end

@testset "unrolling" begin
    N = 1000
    A = reshape(convert.(Float64, collect(1:N^2)), N, N)
    B = reshape(convert.(Float64, collect(1:N^2)), N, N)
    C = zeros(N, N)
    @time refC = A * B
    C = zeros(N, N)
    @time unroll_naive!(C, A, B)
    C = zeros(N ,N)
    @time unroll_k!(C, A, B)
    @test C ≈ refC
    C = zeros(N, N)
    @time unroll_i!(C, A, B)
    @test C ≈ refC
    C = zeros(N, N)
    @time unroll_j!(C, A, B)
    @test C ≈ refC
    C = zeros(N, N)
    @time unroll_ij!(C, A, B)
    @test C ≈ refC
end


text:
using SymbolServer, Pkg
using Test

@testset "SymbolServer" begin
    server = SymbolServerProcess()
    @test server isa SymbolServerProcess
    @test server.context isa Pkg.Types.Context
    SymbolServer.get_context(server)
    @test server.context isa Pkg.Types.Context

    pkgid = SymbolServer.PackageID("SymbolServer", string(SymbolServer.context_deps(server.context)["SymbolServer"]))
    pkgdeps = SymbolServer.pkg_deps(pkgid, server.context)
    @test all(d in keys(pkgdeps) for d in ("Pkg", "SHA", "Serialization", "Test"))
    @test SymbolServer.pkg_ver(pkgid, server.context) isa String

    depot = SymbolServer.load_core()
    @test !isempty(depot["Base"].vals)
    @test !isempty(depot["Core"].vals)

    SymbolServer.load_core(server)

    SymbolServer.import_package_names(pkgid, depot, server.context)
    @test any(p->p[2].name == "SymbolServer", depot)
    @test any(p->p[2].name == "Pkg", depot)
    @test any(p->p[2].name == "SHA", depot)

    SymbolServer.load_package(server, pkgid)

    SymbolServer.load_core(server)
    loaded_pkgs = SymbolServer.load_all(server)

    @test any(p->p[2] == "SymbolServer", loaded_pkgs)
    @test any(p->p[2] == "Pkg", loaded_pkgs)
    @test any(p->p[2] == "SHA", loaded_pkgs)

    SymbolServer.getstore(server)
    # @info keys(server.depot)
    # @test haskey(server.depot, "Base")
    # @test haskey(server.depot, "Core")
    # @test haskey(server.depot, "SymbolServer")
    kill(server)
end


text:
using LinearAlgebra

eye(n) = Matrix{Float64}(I,n,n)
eye(n,m) = Matrix{Float64}(I,n,m)

function e(i, dims)
    E = Matrix{Float64}(I, dims, dims)
    return E[:,i]
end


text:
# Pinger ensures that pings are sent regularly. It is a thin wrapper around a timer and an
# AbstractClientProtocol object.
# Ponger ensures that pongs are received regularly. The logic tells Ponger whenever a pong response
# is received, and Ponger in return tells the logic when a pong has not been received in the
# expected timeframe.

import DandelionWebSockets: ClientPingRequest, handle, PongMissed

type FakeLogic <: AbstractClientProtocol
    ping_requests::Int
    pongs_missed::Int

    FakeLogic() = new(0, 0)
end

handle(logic::FakeLogic, ::ClientPingRequest) = logic.ping_requests += 1
handle(logic::FakeLogic, ::PongMissed) = logic.pongs_missed += 1

facts("Pinger/Ponger") do
    context("Sending pings") do
        logic = FakeLogic()
        interval = 0.1
        pinger = Pinger(interval)
        attach(pinger, logic)

        sleep(3 * interval + 0.02)
        stop(pinger)

        @fact logic.ping_requests --> 3

        # Check that the pinger stopped.
        sleep(3 * interval)
        @fact logic.ping_requests --> 3
    end

    context("Missed pong") do
        logic = FakeLogic()
        timeout = 0.1

        ponger = Ponger(timeout)
        attach(ponger, logic)

        for i in range(1, 3)
            ping_sent(ponger)
            sleep(timeout + 0.02)
        end

        @fact logic.pongs_missed --> 3

        sleep(3 * timeout)
        @fact logic.pongs_missed --> 3
    end

    context("Disconnect only after three pong misses") do
        logic = FakeLogic()
        timeout = 0.1

        ponger = Ponger(timeout; misses=3)
        attach(ponger, logic)

        ping_sent(ponger)
        sleep(timeout + 0.02)

        @fact logic.pongs_missed --> 0

        ping_sent(ponger)
        sleep(timeout + 0.02)
        @fact logic.pongs_missed --> 0

        ping_sent(ponger)
        sleep(timeout + 0.02)
        @fact logic.pongs_missed --> 1

        sleep(3 * timeout)
        @fact logic.pongs_missed --> 1

        for i in range(1, 6)
            ping_sent(ponger)
            sleep(timeout + 0.02)
        end

        @fact logic.pongs_missed --> 3
    end

    context("Only expect a pong after a ping") do
        logic = FakeLogic()
        timeout = 0.1

        ponger = Ponger(timeout)
        attach(ponger, logic)

        ping_sent(ponger)
        sleep(3 * timeout + 0.02)

        @fact logic.pongs_missed --> 1
    end

    context("Receiving pongs") do
        logic = FakeLogic()
        timeout = 0.1

        ponger = Ponger(timeout)
        attach(ponger, logic)

        for i in range(1, 5)
            ping_sent(ponger)
            sleep(timeout / 2.0)
            pong_received(ponger)
        end

        @fact logic.pongs_missed --> 0
    end
end


text:
module Dating

using Distributions, Roots

export DatingEstimate, coaltime, SimpleEstimate, SpeedDating, SDResult, upper,
       lower

"""
An abstract dating estimate type.

Types inheriting from DatingEstimate should have the following methods defined:

* lower: Returns the lower bound of the dating estimate.
* upper: Returns the upper bound of the dating estimate.

They may define their own show methods.
"""
abstract DatingEstimate

Base.print(io::IO, de::DatingEstimate) = println(io, "$(lower(de)) ... $(upper(de))")
Base.show(io::IO, de::DatingEstimate) = println(io, "Coalesence time estimate between two sequences:\nt lies between $(lower(de)) and $(upper(de))")
Base.in(val::Float64, de::DatingEstimate) = val <= upper(de) && val >= lower(de)

# Different coalescence time estimation algorithms.
# method coaltime is dispatched according to arguments or type of
immutable SimpleEstimate end

"""
    coaltime(len::Int, nmut::Int, mu::Float64, ::Type{SimpleEstimate})

Compute the coalescence time between two sequences by assuming a mutation rate,
and then assuming the divergence between two sequences is 2 * mu * t

# Examples
```julia
coaltime(50, 17, 10e-9, SimpleEstimate)
```
"""
function coaltime(len::Int, nmut::Int, mu::Float64, ::Type{SimpleEstimate})
    return (nmut / len) / (2 * mu)
end

immutable SpeedDating end

"""
SDResult is a simple datatype to store the results of the coaltime
function.

coaltime estimates date range that may be considered a 95% confidence range in
which the true coalsescence time between two sequences lies.

The type store the 5%, 50%, and 95% values for the range.
"""
immutable SDResult <: DatingEstimate
    lower::Float64
    middle::Float64
    upper::Float64
end

lower(x::Float64) = x
lower(x::SDResult) = x.lower
upper(c::Float64) = x
upper(x::SDResult) = x.upper

function Base.show(io::IO, de::SDResult)
    println(io, "Coalescence time estimate:\n5%: $(de.lower), 95%: $(de.upper)")
end

function binomzero(p0::Float64, N::Int, B::Int)
    f(p::Float64) = cdf(Binomial(N, p), B) - p0
    return fzero(f, 0, 1)
end

"""
    coaltime(len::Int, nmut::Int, mu::Float64, ::Type{SpeedDating})

Compute the coalescence time between two sequences by modelling the process of
mutation accumulation between two sequences as a bernoulli process.

This method was first described in the paper:
Ward, B. J., & van Oosterhout, C. (2016). Hybridcheck: Software for the rapid detection, visualization and dating of recombinant regions in genome sequence data. Molecular Ecology Resources, 16(2), 534–539.

len is the length of the two aligned sequences, nmut is the number of mutations
and mu is the assumed mutation rate.

# Examples
```julia
coaltime(50, 17, 10e-9, SpeedDating)
```
"""
function coaltime(len::Int, nmut::Int, mu::Float64, ::Type{SpeedDating})
    ninetyfive = ceil(binomzero(0.05, len, nmut) / (2 * mu))
    fifty = ceil(binomzero(0.5, len, nmut) / (2 * mu))
    five = ceil(binomzero(0.95, len, nmut) / (2 * mu))
    return SDResult(five, fifty, ninetyfive)
end

end


text:
export periapsis, apoapsis, semimajoraxis, eccentricity
export meananomaly, trueanomaly, eccentricanomaly
export orbitalperiod, orbitaldistance, orbit

"""
    periapsis(a, e)

Compute the periapsis (closest approach) distance using semi-major axis and eccentricity
"""
periapsis(a, e) = a*(1 - e)

"""
    apoapsis(a, e)

Compute the apoapsis (farthest distance) distance using semi-major axis and eccentricity
"""
apoapsis(a, e) = a*(1 + e)

"""
    semimajoraxis(T, m)

Compute the [semi-major axis](https://en.wikipedia.org/wiki/Semi-major_and_semi-minor_axes) of an orbit
"""
semimajoraxis(T, m) = (𝐆*m*T^2/(4π^2))^(1/3)

"""
    eccentricity(rₚ, rₐ)

Compute [eccentricity](https://en.wikipedia.org/wiki/Orbital_eccentricity)
"""
eccentricity(rₚ, rₐ) = (rₐ - rₚ)/(rₐ + rₚ)

"""
    meananomaly(E, e)

Compute the [mean anomaly](https://en.wikipedia.org/wiki/Mean_anomaly)
"""
meananomaly(E, e) = E - e*sin(E)

"""
    trueanomaly(E, e)

Compute the [true anomaly](https://en.wikipedia.org/wiki/True_anomaly)
"""
function trueanomaly(E, e)
    f = 2*atan(sqrt((1 + e)/(1 - e))*tan(E/2))
    #use [0,2π] instead of [-π,π]
    f < 0 ? f + 2π : f
end

"""
    trueanomaly(t, a, m, e)

Compute the [true anomaly](https://en.wikipedia.org/wiki/True_anomaly)
"""
trueanomaly(t, a, m, e) = trueanomaly(eccentricanomaly(t, a, m, e), e)

"""
    eccentricanomaly(t, a, m, e)

Numerically compute the [eccentric anomaly](https://en.wikipedia.org/wiki/Eccentric_anomaly) using [Kepler's equation](https://en.wikipedia.org/wiki/Kepler%27s_equation)
"""
function eccentricanomaly(t, a, m, e)
    @assert t >= 0 "time must be positive"
    #Kepler's Third Law
    T = orbitalperiod(a, m)
    #definition of mean anomaly
    M = 2π*rem(t, T)/T
    #eccentric anomaly must be found numerically
    E = regulafalsi((E,p)->meananomaly(E, e) - M, 0, 2π)
    return E
end

"""
    orbitalperiod(a, m)

[Kepler's Third Law](https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion#Third_law) describing the orbital period of an elliptical orbit
"""
orbitalperiod(a, m) = 2π*√(a^3/(𝐆*m))

"""
    orbitaldistance(a, f, e)

Compute the distance of a planet from its host
"""
orbitaldistance(a, f, e) = a*(1 - e^2)/(1 + e*cos(f))

"""
    orbitaldistance(t, a, m, e)

Compute the distance of a planet from its host, assuming the planet is at periapsis at t=0
"""
orbitaldistance(t, a, m, e) = orbitaldistance(a, trueanomaly(t, a, m, e), e)

"""
    orbit(a, m, e, N=1000)

Create a distance time-series of `N` points for an elliptical orbit, returning vectors for time, distance, and true anomaly
"""
function orbit(a, m, e, N::Int=1000)
    T = orbitalperiod(a, m)
    t = LinRange(0, T, N+1)[1:end-1]
    f = trueanomaly.(t, a, m, e)
    r = orbitaldistance.(a, f, e)
    return t, r, f
end


text:
using TimeSeries
using Dates
using PowerSystems

dates = collect(
    DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime(
        "1/1/2024  23:00:00",
        "d/m/y  H:M:S",
    ),
)

nodes4_tcrd() = [
    Bus(
        1,
        "Bus 1",
        "REF",
        0.0,
        1.0,
        (min = 0.9, max = 1.1),
        220,
        nothing,
        nothing),
    Bus(
        2,
        "Bus 2",
        "PV",
        0.0,
        1.0,
        (min = 0.9, max = 1.1),
        220,
        nothing,
        nothing,
    ),
    Bus(
        3,
        "Bus 3",
        "PV",
        0.0,
        1.0,
        (min = 0.9, max = 1.1),
        220,
        nothing,
        nothing,
    ),
    Bus(
        4,
        "Bus 4",
        "PV",
        0.0,
        1.0,
        (min = 0.9, max = 1.1),
        220,
        nothing,
        nothing,
    ),
]
branches4_tcrd(nodes4_tcrd) = [
    Line(
        "Line1",
        true,
        0.0,
        0.0,
        Arc(from = nodes4_tcrd[1], to = nodes4_tcrd[2]),
        0.0,
        0.1,
        (from = 0.0, to = 0.0),
        10.20,
        1.0,
    ),
    Line(
        "Line2",
        true,
        0.0,
        0.0,
        Arc(from = nodes4_tcrd[1], to = nodes4_tcrd[3]),
        0.0,
        0.1,
        (from = 0.0, to = 0.0),
        12.0,
        1.0,
    ),
    Line(
        "Line3",
        true,
        0.0,
        0.0,
        Arc(from = nodes4_tcrd[2], to = nodes4_tcrd[3]),
        0.0,
        0.1,
        (from = 0.0, to = 0.0),
        12.0,
        1.0,
    ),
    Line(
        "Line4",
        true,
        0.0,
        0.0,
        Arc(from = nodes4_tcrd[2], to = nodes4_tcrd[4]),
        0.0,
        0.1,
        (from = 0.0, to = 0.0),
        12.0,
        1.0,
    ),
    Line(
        "Line5",
        true,
        0.0,
        0.0,
        Arc(from = nodes4_tcrd[3], to = nodes4_tcrd[4]),
        0.0,
        0.1,
        (from = 0.0, to = 0.0),
        12.0,
        1.0,
    ),
]

thermal_generators4_tcrd(nodes4_tcrd) = [
    ThermalStandard(
        name = "GBus1",
        available = true,
        status = true,
        bus = nodes4_tcrd[1],
        active_power = 0.50,
        reactive_power = 0.10,
        rating = 0.500,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 1.0),
        reactive_power_limits = (min = -0.25, max = 0.25),
        time_limits = nothing,
        ramp_limits = nothing,
        operation_cost = ThreePartCost((0.175000, 10.0000), 0.0, 0.0, 0.0),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "GBus2",
        available = true,
        status = true,
        bus = nodes4_tcrd[2],
        active_power = 0.50,
        reactive_power = 0.10,
        rating = 0.500,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 0.50),
        reactive_power_limits = (min = -0.25, max = 0.25),
        time_limits = nothing,
        ramp_limits = nothing,
        operation_cost = ThreePartCost((0.497000, 10.0000), 0.0, 0.0, 0.0),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "GBus3",
        available = true,
        status = true,
        bus = nodes4_tcrd[3],
        active_power = 0.50,
        reactive_power = 0.10,
        rating = 2.00,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 2.00),
        reactive_power_limits = (min = -1.00, max = 1.00),
        time_limits = nothing,
        ramp_limits = nothing,
        operation_cost = ThreePartCost((0.260000, 20.000), 0.0, 0.0, 0.0),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "GBus4",
        available = true,
        status = true,
        bus = nodes4_tcrd[4],
        active_power = 0.50,
        reactive_power = 0.10,
        rating = 2.00,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 2.00),
        reactive_power_limits = (min = -1.00, max = 1.00),
        time_limits = nothing,
        ramp_limits = nothing,
        operation_cost = ThreePartCost((0.325000, 20.000), 0.0, 0.0, 0.0),
        base_power = 100.0,
    ),
]

#time series per zone
loadz1_ts = ones(Float64, 24)

loads4_tcrd(nodes4_tcrd) = [
    PowerLoad(
        "Bus3",
        true,
        nodes4_tcrd[3],
        LoadModels.ConstantPower,
        1.0,
        0.0,
        100.0,
        1.0,
        0.0,
    ),
    PowerLoad(
        "Bus4",
        true,
        nodes4_tcrd[4],
        LoadModels.ConstantPower,
        1.0,
        0.0,
        100.0,
        1.0,
        0.0,
    ),
]

timeseries_DA4_tcrd = [
    TimeArray(dates, loadz1_ts),
    TimeArray(dates, loadz1_ts),
]

function filter_kwargs(; kwargs...)
    system_kwargs = filter(x -> in(first(x), PSY.SYSTEM_KWARGS), kwargs)
    return system_kwargs
end

function build_c_sys4_tcrd(; kwargs...)
    sys_kwargs = filter_kwargs(; kwargs...)
    nodes = nodes4_tcrd()
    c_sys4_tcrd = PSY.System(
        100.0,
        nodes,
        thermal_generators4_tcrd(nodes),
        loads4_tcrd(nodes),
        branches4_tcrd(nodes);
        time_series_in_memory = get(sys_kwargs, :time_series_in_memory, true),
        sys_kwargs...,
    )

    if get(kwargs, :add_forecasts, true)
        forecast_data = SortedDict{Dates.DateTime, TimeSeries.TimeArray}()
        for (ix, l) in enumerate(PSY.get_components(PowerLoad, c_sys4_tcrd))
            ini_time = TimeSeries.timestamp(timeseries_DA4_tcrd[ix])[1]
           forecast_data[ini_time] = timeseries_DA4_tcrd[ix]
           PSY.add_time_series!(
               c_sys4_tcrd,
               l,
               PSY.Deterministic("max_active_power", forecast_data),
           )
       end
    end

    return c_sys4_tcrd
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
module NewMatteApp

using Matte

const title = "My New Matte App"

function ui()
    sidebar_layout(
        side_panel(
            h1("Hello World!")
        ),
        main_panel(
            h1("Output here")
        )
    )
end

module Server

function my_output(my_input)
    my_input
end

end

end


text:
@testset "Norms" begin
    x = [2.0im, 3.0-1im, 5.0+2.0im]
    y = [-2.0im, 3.0-1im, 5.0+2.0im]
    Norm = HC.InfNorm()
    @test Norm isa HC.AbstractNorm
    @test Norm(x) ≈ abs(5.0+2.0im)
    @test Norm(x, x) ≈ 0.0
    @test Norm(x, y) ≈ abs(4im)

    WeightedNorm = HC.WeightedNorm(Norm, x)
    @test HC.weights(WeightedNorm) == ones(3)
    WeightedNorm .= 2.0
    @test HC.weights(WeightedNorm) == [2.0, 2, 2]
    WeightedNorm[1] = 4.0
    @test WeightedNorm isa HC.AbstractNorm
    @test WeightedNorm(x) ≈ 0.5abs(5+2.0im)
    @test WeightedNorm(x, x) ≈ 0.0
    @test WeightedNorm(x, y) ≈ 0.25abs(4im)

    Norm = HC.EuclideanNorm()
    @test Norm isa HC.AbstractNorm

    @test Norm(x) ≈ LinearAlgebra.norm(x, 2)
    @test Norm(x, x) ≈ 0.0
    @test Norm(x, y) ≈ abs(4im)

    WeightedNorm = HC.WeightedNorm([4.0, 2.0, 2.0], Norm)
    @test WeightedNorm isa HC.AbstractNorm
    @test WeightedNorm(x) ≈ sqrt(abs2(x[1])/16+abs2(x[2])/4+abs2(x[3])/4)
    @test WeightedNorm(x, x) ≈ 0.0
    @test WeightedNorm(x, y) ≈ sqrt(abs2(4im)/16)
end


text:
# all of these must be added as a test dependency!!
using BenchmarkTools
using Random
using Mycelia

const SUITE = BenchmarkGroup()

SUITE["utf8"] = BenchmarkGroup(["string", "unicode"])
teststr = String(join(rand(MersenneTwister(1), 'a':'d', 10^4)))
SUITE["utf8"]["replace"] = @benchmarkable replace($teststr, "a" => "b")
SUITE["utf8"]["join"] = @benchmarkable join($teststr, $teststr)
SUITE["utf8"]["plots"] = BenchmarkGroup()

SUITE["trigonometry"] = BenchmarkGroup(["math", "triangles"])
SUITE["trigonometry"]["circular"] = BenchmarkGroup()
for f in (sin, cos, tan)
    for x in (0.0, pi)
        SUITE["trigonometry"]["circular"][string(f), x] = @benchmarkable ($f)($x)
    end
end

SUITE["trigonometry"]["hyperbolic"] = BenchmarkGroup()
for f in (sin, cos, tan)
    for x in (0.0, pi)
        SUITE["trigonometry"]["hyperbolic"][string(f), x] = @benchmarkable ($f)($x)
    end
end

# any time that I update the benchmarks or want to modify the tunings,
# just delete the tune.json file and then re run this script
# https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/doc/manual.md#caching-parameters
tuning_file = "$(@__DIR__)/tune.json"
if !isfile(tuning_file)

    # tune the suite to configure benchmark parameters
    tune!(SUITE);

    # save the suite's parameters using a thin wrapper
    # over JSON (this wrapper maintains compatibility
    # across BenchmarkTools versions)
    BenchmarkTools.save("tune.json", params(SUITE));
end

max_stars_repo_name: snystrom/Mycelia

max_issues_repo_name: snystrom/Mycelia

max_forks_repo_name: cjprybol/Eisenia


text:
#!/usr/bin/env julia

println("Loading dependencies...")

using Plots
using LibSerialPort

# Plot backend
gr(size=(900, 500))

address = length(ARGS) == 1 ? ARGS[1] : "/dev/cu.usbmodem1421"

mcu = open(address, 115200)

# Number of points to display
n = 500

# Values for horizontal (time) axis and data
x = Float64[]
y, s = Float64[], Float64[]

function readmcu(sp::SerialPort)
    # Get message as a String, split to array
    msg = readuntil(mcu, '\n', 100)
    data = zeros(Float64, 3)
    try
        data = [parse(Float64, x) for x in split(msg)]
    catch
        println("Skipping: $msg")
        return
    end
    length(data) == 3 || return
    time, y1, y2 = data

    # println("$time, $y1, $y2")

    # Manage arrays as FIFO queues
    push!(x, time/1000) # ms to seconds
    push!(y, y1/1023)
    push!(s, y2/1023)
    if length(x) > n
        for v in (x, y, s)
            shift!(v)
        end
    end
end

function rtplot(x, y, s)
    plot(x, [y,s],
         title="Data stream", # Not displaying with GR backend (?)
         xaxis="time [s]",
         yaxis=("example [units]", (-0.2, 1.2)), # title, limits
         labels=(["Measured","Filtered"]), # legend
         line=([1 2], [:steppre :path], ["gray" "crimson"], [1.0 0.5]), # lineweights, linetypes, colors, alphas
         fill=(0, [0.1 0], "dodgerblue"), # y-origin, alpha, colors
         )
    # Update the plot window
    gui()
end

println("Creating plot...")
readmcu(mcu)
rtplot(x,y,s)

println("Ready")

while true

    @async begin
        a = readline(STDIN)
        a == "q\n" && quit()
        write(mcu, "$a")
    end

    @async begin
        readmcu(mcu)
        rtplot(x,y,s)
    end

    # Give the queued tasks a chance to run
    sleep(0.001)
end

max_stars_repo_name: andrewadare/realtime

max_issues_repo_name: andrewadare/realtime

max_forks_repo_name: andrewadare/realtime


text:
# isdefined(Base, :__precompile__) && __precompile__()

module Complementarity

# package code goes here

using JuMP
using Base.Meta

import PATHSolver, MathProgBase, NLsolve

export  MCPModel, MCPData, ComplementarityType,
        complements, solveMCP,
        @complements

mpec_tol = 1e-8

include("mcp.jl")
include("mpec.jl")

end # module


text:
module GridapMHD

using Gridap
using Gridap.Geometry:get_grid_topology, get_node_coordinates, get_face_nodes

using GridapODEs
using GridapODEs.TransientFETools
using GridapODEs.ODETools
import GridapODEs.TransientFETools: ∂t

using GridapEmbedded
using GridapGmsh

using ForwardDiff
using TimerOutputs
using IncompleteLU
using LineSearches: BackTracking, Static

include("InductionlessMHD.jl")

include("LinearSolvers.jl")

using GridapMHD.LinearSolvers

# LinearSolvers exports
export GmresSolver
export symbolic_setup
export numerical_setup
export numerical_setup!
export solve!

## IO functions exports
export Verbosity
export startPVD
export write_timestep
export closePVD

# Handy functions exports
export compute_u_j_errors
export streching
export cuboid
export add_entity!

# Handy functions
_l2(v, dΩ) = sqrt(sum(∫(v⋅v)*dΩ ))
_h1(v, dΩ) = sqrt(sum(∫(v⋅v + inner(∇(v),∇(v)) )*dΩ ))
_hdiv(v, dΩ) = sqrt(sum(∫(v⋅v + inner((∇⋅v),(∇⋅v)) )*dΩ ))

function compute_u_j_errors(uh, jh, u, j, dΩ)
  eu = uh - u
  ej = jh - j

  eu_l2 = _l2(eu, dΩ)
  ej_l2 = _l2(ej, dΩ)

  return eu_l2, ej_l2
end

# IO functions
struct Verbosity
  items::Dict
  function Verbosity(level::Int = 1)
    items = Dict{String,Bool}()
    if level == 0 # Silent
      items["preprocessSteps"] = false
      items["linearSolver"] = false
      items["nonlinearSolver"] = false
      items["timeStep"] = false
    elseif level == 1 # Standard
      items["preprocessSteps"] = false
      items["linearSolver"] = false
      items["nonlinearSolver"] = true
      items["timeStep"] = true
    elseif level == 2 # Verbose
      items["preprocessSteps"] = false
      items["linearSolver"] = true
      items["nonlinearSolver"] = true
      items["timeStep"] = true
    elseif level == 3 # Debug
      items["preprocessSteps"] = true
      items["linearSolver"] = true
      items["nonlinearSolver"] = true
      items["timeStep"] = true
    else
      @error "Unavailable verbosity level. Choose level 0 to 3, from less to more verbosity"
    end
    new(items)
  end
  function Verbosity(items::Dict)
    new(items)
  end
end

function preprocessStepOutput(v::Verbosity,msg)
  if v.items["preprocessSteps"]
    println(msg)
  end
end

function linearSolverOutput(v::Verbosity)
  return v.items["linearSolver"]
end

function nonlinearSolverOutput(v::Verbosity)
  return v.items["nonlinearSolver"]
end

function timeStepOutput(v::Verbosity,dt,t,it)
  if v.items["timeStep"]
    println("######### TIME STEP $it #########")
    println("Current time: $t")
    println("Time step length: $dt")
    println("#################################")
  end
end

function startPVD(filename)
  rm(filename,force=true,recursive=true)
  mkdir(filename)
  pvdheader  = """<?xml version="1.0"?>
<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">
  <Collection>\n"""

  io = open(filename * ".pvd", "w")
  write(io,pvdheader)
  close(io)
end

function write_timestep(filename,t,trian;kwargs...)
  pvdcontent = """    <DataSet timestep=""" * '"'
  pvdcontent *= string(t) * '"' * """ group="" part="0" file=""" * '"'
  pvdcontent *= filename*"""/time_"""*string(t)*""".vtu"/>\n"""

  writevtk(trian, filename*"/time_"*string(t)*".vtu";kwargs...)

  io = open(filename * ".pvd", "a")
  write(io,pvdcontent)
  close(io)
end

function closePVD(filename)
  pvdfooter  = "  </Collection>\n</VTKFile>"
  io = open(filename * ".pvd", "a")
  write(io,pvdfooter)
  close(io)
end

function streching(coord;domain=(0.0,1.0,0.0,1.0,0.0,1.0),factors=(1.0,1.0,1.0),
  dirs=(1,2,3))
ncoord = collect(coord.data)
for (i,dir) in enumerate(dirs)
  ξ0 = domain[i*2-1]
  ξ1 = domain[i*2]
  l =  ξ1 - ξ0
  if ξ0 <= coord[dir] <= ξ1
    ξx = (coord[dir] - ξ0)/l
    if ξx < 1e-12
      ncoord[dir] = ξ0
    else
      ncoord[dir] = ((ξx/sqrt(abs(ξx)))^factors[i])*l + ξ0
    end
  end
end
VectorValue(ncoord)
end

# CSG geometry definitions
function cuboid(;dx=1,dy=1,dz=1,x0=Point(0,0,0),name="cuboid",
  faces=["face$i" for i in 1:6])
  e1 = VectorValue(1,0,0)
  e2 = VectorValue(0,1,0)
  e3 = VectorValue(0,0,1)

  plane1 = plane(x0=x0-0.5*dz*e3,v=-e3,name=faces[1])
  plane2 = plane(x0=x0+0.5*dz*e3,v=+e3,name=faces[2])
  plane3 = plane(x0=x0-0.5*dy*e2,v=-e2,name=faces[3])
  plane4 = plane(x0=x0+0.5*dy*e2,v=+e2,name=faces[4])
  plane5 = plane(x0=x0-0.5*dx*e1,v=-e1,name=faces[5])
  plane6 = plane(x0=x0+0.5*dx*e1,v=+e1,name=faces[6])

  geo12 = intersect(plane1,plane2)
  geo34 = intersect(plane3,plane4)
  geo56 = intersect(plane5,plane6)

  intersect(intersect(geo12,geo34),geo56,name=name)
end

# Function that given a filter set a name to those elements that satisfy the
# conditions of the filter
function add_entity!(model,in,name)
  labels = get_face_labeling(model)
  node_coordinates = get_node_coordinates(model)
  entity = num_entities(labels) + 1
  for d in 0:num_dims(model)-1
    facets = get_face_nodes(model,d)
    for (i,facet) in enumerate(facets)
      coord = sum(node_coordinates[facet])/length(facet)
      if in(coord)
        labels.d_to_dface_to_entity[d+1][i] = entity
      end
    end
  end
  add_tag!(labels,name,[entity])
end

# Benchmarks
include("Hunt.jl")

include("Shercliff.jl")

include("ConductiveThinWall.jl")

include("TransientDuctFlow.jl")

include("DriverInductionlessMHD.jl")

include("TransientDriverInductionlessMHD.jl")

end # module


text:
immutable TransposeMap{T}<:AbstractLinearMap{T}
    lmap::AbstractLinearMap{T}
end
immutable CTransposeMap{T}<:AbstractLinearMap{T}
    lmap::AbstractLinearMap{T}
end

# transposition behavior of AbstractLinearMap objects
Base.transpose(A::TransposeMap)=A.lmap
Base.ctranspose(A::CTransposeMap)=A.lmap

Base.transpose{T}(A::AbstractLinearMap{T})=issym(A) ? A : TransposeMap{T}(A)
Base.ctranspose{T<:Real}(A::AbstractLinearMap{T})=transpose(A)
Base.ctranspose{T}(A::AbstractLinearMap{T})=ishermitian(A) ? A : CTransposeMap{T}(A)

# properties
Base.size(A::Union{TransposeMap,CTransposeMap},n)=(n==1 ? size(A.lmap,2) : (n==2 ? size(A.lmap,1) : error("AbstractLinearMap objects have only 2 dimensions")))
Base.size(A::Union{TransposeMap,CTransposeMap})=(size(A.lmap,2),size(A.lmap,1))
Base.isreal(A::Union{TransposeMap,CTransposeMap})=isreal(A.lmap)
Base.issym(A::Union{TransposeMap,CTransposeMap})=issym(A.lmap)
Base.ishermitian(A::Union{TransposeMap,CTransposeMap})=ishermitian(A.lmap)
Base.isposdef(A::Union{TransposeMap,CTransposeMap})=isposdef(A.lmap)

# comparison of TransposeMap objects
==(A::TransposeMap,B::TransposeMap)=A.lmap==B.lmap
==(A::CTransposeMap,B::CTransposeMap)=A.lmap==B.lmap

# multiplication with vector: should be ok for all lmap values, also those for which Transpose/CTranpose should not have been created
Base.A_mul_B!(y::AbstractVector,A::TransposeMap,x::AbstractVector)=(issym(A.lmap) ? Base.A_mul_B!(y,A.lmap,x) : Base.At_mul_B!(y,A.lmap,x))
*(A::TransposeMap,x::AbstractVector)=(issym(A.lmap) ? *(A.lmap,x) : Base.At_mul_B(A.lmap,x))

Base.At_mul_B!(y::AbstractVector,A::TransposeMap,x::AbstractVector)=Base.A_mul_B!(y,A.lmap,x)
Base.At_mul_B(A::TransposeMap,x::AbstractVector)=*(A.lmap,x)

Base.Ac_mul_B!(y::AbstractVector,A::TransposeMap,x::AbstractVector)=isreal(A.lmap) ? Base.A_mul_B!(y,A.lmap,x) : (Base.A_mul_B!(y,A.lmap,conj(x));conj!(y))
Base.Ac_mul_B(A::TransposeMap,x::AbstractVector)=isreal(A.lmap) ? *(A.lmap,x) : conj!(*(A.lmap,conj(x)))

Base.A_mul_B!(y::AbstractVector,A::CTransposeMap,x::AbstractVector)=(ishermitian(A.lmap) ? Base.A_mul_B!(y,A.lmap,x) : Base.Ac_mul_B!(y,A.lmap,x))
*(A::CTransposeMap,x::AbstractVector)=(ishermitian(A.lmap) ? *(A.lmap,x) : Base.Ac_mul_B(A.lmap,x))

Base.At_mul_B!(y::AbstractVector,A::CTransposeMap,x::AbstractVector)=isreal(A.lmap) ? Base.A_mul_B!(y,A.lmap,x) : (Base.A_mul_B!(y,A.lmap,conj(x));conj!(y))
Base.At_mul_B(A::CTransposeMap,x::AbstractVector)=isreal(A.lmap) ? *(A.lmap,x) : conj!(*(A.lmap,conj(x)))

Base.Ac_mul_B!(y::AbstractVector,A::CTransposeMap,x::AbstractVector)=Base.A_mul_B!(y,A.lmap,x)
Base.Ac_mul_B(A::CTransposeMap,x::AbstractVector)=*(A.lmap,x)


text:
# High-level request functions which create and send a request,
# then return Seis.Events

const _Period = Union{Dates.Period, Dates.CompoundPeriod}

#
# get_stations
#

"""
    get_stations(; server=$(DEFAULT_SERVER), verbose=true, T=Float64, kwargs...) -> stations

Query `server` for stations matching the keyword arguments in `kwargs`,
which are passed to [`FDSNStation`](@ref).  `stations` is a
`Vector{Seis.GeogStation{T}}`, where various useful fields are added to the
`.meta` field of each station.  These keyword arguments are supported by
all `get_stations` methods.

For search options and the format of `kwargs`, see [`FDSNStation`](@ref).

The element type of the `Seis.Station`s returned is set by `T`.

# Example
Find all stations in the XM network between 2012 and 2014:
```
julia> using Dates

julia> get_stations(network="XM", starttime=DateTime(2012), endtime=DateTime(2015))
[ Info: Request status: Successful request, results follow
21-element Array{Seis.Station{Float64,Seis.Geographic{Float64}},1}:
 Station: XM.A01E.., lon: 38.78384, lat: 7.78823, elev: 1921.0, meta: 5 keys
 Station: XM.A02E.., lon: 38.799709, lat: 7.86046, elev: 1643.0, meta: 5 keys
 Station: XM.A03E.., lon: 38.76218, lat: 7.81199, elev: 1952.0, meta: 5 keys
 Station: XM.A04E.., lon: 38.688702, lat: 7.76335, elev: 1640.0, meta: 5 keys
 Station: XM.A05E.., lon: 38.726269, lat: 7.84464, elev: 1652.0, meta: 5 keys
...
```

# Populated `Station` fields

The amount of information contained within each `Station` returned
depends on the `level` of information requested and the `format` of
the response (see [`FDSNStation`](@ref)).

If `format` is `"text"`, then only a limited number of extra information
is added to the `meta` field of each station.  By default, the servers returns
information in StationXML format, which may be bare or comprehensive.
The StationXML information for the network, station or channel is added
to the `meta.stationxml` field.

If `level` is `"network"`, then stations are not real 'stations' and contain
only the network fields.  Likewise, for `level="station"`, the channel code
(`cha` field) of each station is empty and no information about channel
orientation is available.  Therefore, it is usually best to specify
`level="channel"` if you plan to associate the `stations` with time series
data.

If you use the `code` keyword argument to specify the stations for which to
search (e.g., `code="XM.A*..*Z"`; see [`FDSNStation`](@ref)), then `level`
is automatically set to `"channel"`.  Override this if desired by passing
a value to `level` explicitly.
"""
function get_stations(; server=DEFAULT_SERVER, verbose=true, T=Float64, kwargs...)
    request = FDSNStation(; kwargs...)
    response = get_request(request; server=server, verbose=verbose)
    stations = parse_station_response(T, request, response, server)
end

"""
    get_stations(event::Seis.GeogEvent, starttime, endtime; longitude=event.lon, latitude=event.lat,
                 starttime=event.time, endtime=event.time, kwargs...) -> stations

Find stations relative to an `event`, typically constraining the station locations
to be a certain distance from `event`, or requiring stations to be active at
the time of the event.

This method queries `server` for stations matching the options to [`FDSNStation`](@ref)
in `kwargs`, but with the following defaults set:

- `longitude` and `latitude` are set to the event's coordinates.  This means that
  automatically, setting `minradius` and `maxradius` will search around the event.
- `starttime` and `endtime` are set to the event origin time), meaning only
  stations or channels which were active at the time of the event will be included.

To override the default behaviour of only matching stations/channels which are
active at the time of the event, set `starttime` and `endtime` to the date
range of interest.

!!! note
    To limit stations to only those with time series information for the period
    specified, use `matchtimeseries=true`.  This is especially useful when
    setting `starttime` and `endtime` using other methods of `get_stations`.

# Example
Get vertical high-gain, broadband and high broadband channels in the II and IU
networks which were active during a deep-focus earthquake on 2018-09-21,
and for which data is present in the IRIS datacentre.  (In this case,
`event` might have been obtained elsewhere or by a call to [`get_events`](@ref),
but we create the event here for simplicity.)
```
julia> using Dates, Seis

julia> event = Event(lon=-179.9776, lat=-17.9071, dep=652.35, time=DateTime(2018, 9, 21, 3, 40, 40, 550))
Seis.Event{Float64,Seis.Geographic{Float64}}:
        lon: -179.9776
        lat: -17.9071
        dep: 652.35
       time: 2018-09-21T03:40:40.55
         id: missing
       meta:

julia> get_stations(event, network="II,IU", channel="BHZ,HHZ", matchtimeseries=true, server="IRIS", format="text")
[ Info: Request status: Successful request, results follow
331-element Array{Station{Float64,Seis.Geographic{Float64}},1}:
 Station: II.AAK.00.BHZ, lon: 74.4942, lat: 42.6375, dep: 0.03, elev: 1633.1, azi: 0.0, inc: 0.0, meta: 7 keys
 Station: II.AAK.10.BHZ, lon: 74.4942, lat: 42.6375, dep: 0.03, elev: 1633.1, azi: 0.0, inc: 0.0, meta: 7 keys
 ⋮
 Station: IU.YSS.00.BHZ, lon: 142.7604, lat: 46.9587, dep: 0.002, elev: 148.0, azi: 0.0, inc: 0.0, meta: 7 keys
 Station: IU.YSS.10.BHZ, lon: 142.7604, lat: 46.9587, dep: 0.002, elev: 148.0, azi: 0.0, inc: 0.0, meta: 7 keys
```

See also: [`FDSNStation`](@ref).
"""
function get_stations(event::GeogEvent, starttime::Union{AbstractString,DateTime}=event.time,
        endtime::Union{AbstractString,DateTime}=event.time; kwargs...)
    default_kwargs = if haskey(kwargs, :minradius) || haskey(kwargs, :maxradius)
        (longitude=event.lon, latitude=event.lat,
            starttime=starttime, endtime=endtime)
    else
        (starttime=starttime, endtime=endtime)
    end
    get_stations(; merge(default_kwargs, kwargs)...)
end

"""
    get_stations(event, start_offset, end_offset; kwargs...)

Specify the time range for data selection as a window around the event time.

`start_offset` and `end_offset` define a data selection time range as a
window around the event time.  They may be given in seconds, or using
`Dates.Period`s or `Dates.CompoundPeriod`s (e.g., `Hour(1)` or `Week(2) + Day(5)`).

Offsets are rounded down to the nearest millisecond.

!!! note
    Note that for `start_offset` and `end_offset`, all positive values indicate
    a time **later** than the event origin time.  Use negative times to mean a
    time before the event.

# Examples
Find broadand channels within 20° of the 2004-12-26 M=9 earthquake with data
between 20 s before and 1800 s after the origin:
```
julia> event = get_events(minmagnitude=8.5, starttime="2004-12-26", endtime="2004-12-27")
[ Info: Request status: Successful request, results follow
1-element Array{Seis.Event{Float64,Seis.Geographic{Float64}},1}:
 Event: lon: 95.9012, lat: 3.4125, dep: 26.1, time: 2004-12-26T00:58:52.05, id: smi:service.iris.edu/fdsnws/event/1/query?originid=3788623, meta: 7 keys

julia> stations = get_stations(first(event), -20, 1800, level="channel", channel="BH?", maxradius=20, matchtimeseries=true)
[ Info: Request status: Successful request, results follow
24-element Array{Seis.Station{Float64,Seis.Geographic{Float64}},1}:
 Station: GE.UGM..BHE, lon: 110.523102, lat: -7.9125, dep: 0.0, elev: 350.0, azi: 90.0, inc: 90.0, meta: 5 keys
 Station: GE.UGM..BHN, lon: 110.523102, lat: -7.9125, dep: 0.0, elev: 350.0, azi: 0.0, inc: 90.0, meta: 5 keys
 Station: GE.UGM..BHZ, lon: 110.523102, lat: -7.9125, dep: 0.0, elev: 350.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: II.COCO.00.BH1, lon: 96.8349, lat: -12.1901, dep: 0.07, elev: 1.0, azi: 141.3, inc: 90.0, meta: 5 keys
 Station: II.COCO.00.BH2, lon: 96.8349, lat: -12.1901, dep: 0.07, elev: 1.0, azi: 231.3, inc: 90.0, meta: 5 keys
 Station: II.COCO.00.BHZ, lon: 96.8349, lat: -12.1901, dep: 0.07, elev: 1.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: II.COCO.20.BHE, lon: 96.8349, lat: -12.1901, dep: 0.0013, elev: 1.0, azi: 90.0, inc: 90.0, meta: 5 keys
...
```

In this example, we used [`get_event`](@ref) to first find the event of interest,
but any existing `Seis.Event` object can be used.  The event and stations could
now be passed to [`get_data`](@ref) to download the data for these stations.

To perform the same query, but for data between one minute before and one hour
after the event:
```
julia> using Dates

julia> get_stations(first(event), -Minute(1), Hour(1), level="channel", channel="BH?", maxradius=20, matchtimeseries=true)
[ Info: Request status: Successful request, results follow
24-element Array{Seis.Station{Float64,Seis.Geographic{Float64}},1}:
 Station: GE.UGM..BHE, lon: 110.523102, lat: -7.9125, dep: 0.0, elev: 350.0, azi: 90.0, inc: 90.0, meta: 5 keys
 Station: GE.UGM..BHN, lon: 110.523102, lat: -7.9125, dep: 0.0, elev: 350.0, azi: 0.0, inc: 90.0, meta: 5 keys
 Station: GE.UGM..BHZ, lon: 110.523102, lat: -7.9125, dep: 0.0, elev: 350.0, azi: 0.0, inc: 0.0, meta: 5 keys
 Station: II.COCO.00.BH1, lon: 96.8349, lat: -12.1901, dep: 0.07, elev: 1.0, azi: 141.3, inc: 90.0, meta: 5 keys
 Station: II.COCO.00.BH2, lon: 96.8349, lat: -12.1901, dep: 0.07, elev: 1.0, azi: 231.3, inc: 90.0, meta: 5 keys
 Station: II.COCO.00.BHZ, lon: 96.8349, lat: -12.1901, dep: 0.07, elev: 1.0, azi: 0.0, inc: 0.0, meta: 5 keys
...
```
"""
get_stations(event::GeogEvent, start_seconds, end_seconds; kwargs...) =
    get_stations(event, event.time + seconds_milliseconds(start_seconds),
        event.time + seconds_milliseconds(end_seconds); kwargs...)

get_stations(event::GeogEvent, start_offset::_Period, end_offset::_Period; kwargs...) =
    get_stations(event, event.time + start_offset, event.time + end_offset; kwargs...)

"""
    get_stations!(stations; starttime=nothing, endtime=nothing, server=$(DEFAULT_SERVER), level="channel", T=Float64)

Fill in missing information about each station in `stations`, which should
be a set of `Seis.Station`s.

This function uses the `net`, `sta`, `loc` and `cha` fields of the stations
to search for metadata matching the channel code of the station by making
a [`FDSNStation`](@ref) request to `server`.  If a station contains
information about the recording period of a station in `meta.startdate`
and `meta.enddate`, then this is used to constrain the search for station
data.

If no operatng time information is available, then this can be passed with
the `starttime` and `endtime` keyword arguments.

If more than one station is a match, then an error is thrown.

If no data are found for one or more stations, a warning is printed but
the function does not throw an error.

# Example
```
julia> using Seis

julia> stas = [Station(net="II", sta="ALE", loc="", cha="BHZ"), Station(net="IU", sta="XMAS", loc="00", cha="BHZ")]

julia> get_stations!([sta])
```
"""
function get_stations!(stations::AbstractArray{<:Seis.GeogStation};
        level="channel", T=Float64,
        server=DEFAULT_SERVER, verbose=true,
        kwargs...
    )
    requests = [
        FDSNStation(;
            (network=s.net,
             station=s.sta,
             location=s.loc,
             channel=s.cha,
             # Use 1-01-01T00:00:00 and 9999-01-01T00:00:00 as impossibly small and
             # large dates, since typemax/typemin values are too large to be
             # parsed by most servers
             starttime=coalesce(s.meta.startdate, DateTime(1800)),
             endtime=coalesce(s.meta.enddate, DateTime(4000)),
             level=level,
             kwargs...)...
        ) for s in stations]
    response = post_request(requests; server=server, verbose=verbose)
    stations_response = parse_station_response(T, first(requests), response, server)
    stations = assign_stations!(stations, stations_response, requests; verbose=verbose)
    stations
end

"""
    assign_stations!(stations, stations_response, requests)

Fill in extra information in each station in `stations` with information in
`stations_response`, matching them up by channel code.  `requests` is a set of
[`FDSNStation`](@ref) requests matched to stations similarly.
"""
function assign_stations!(stations, stations_response, requests; verbose=true)
    for station in stations
        code = Seis.channel_code(station)
        inds = findall(x -> x.net == station.net && x.sta == station.sta &&
                            x.loc == station.loc && x.cha == station.cha,
            stations_response)
        # Skip stations with no matching stations_response
        if length(inds) == 0
            verbose && @warn("No matching station for $code")
            continue
        else
            # Always warn if there are more than two
            length(inds) > 1 &&
                @warn("more than one period found matching $code; taking latest")
            latest_ind = argmax(
                [coalesce(stations_response[i].meta.starttime, typemin(DateTime))
                                for i in inds])
            sta = stations_response[inds[latest_ind]]
            station.lon, station.lat, station.elev = sta.lon, sta.lat, sta.elev
            station.azi, station.inc = sta.azi, sta.inc
            merge!(station.meta, sta.meta)
            # Add corresponding request for this station, since
            # in the call to parse_station_response() we only passed
            # the first request and it will be wrong at this point
            request_ind = findfirst(x -> x.network == station.net &&
                                     x.station == station.sta &&
                                     x.location == station.loc &&
                                     x.channel == station.cha,
                requests)
            # This should be impossible since we create the requests using
            # `stations`, but just in case we are called incorrectly
            request_ind === nothing &&
                throw(ArgumentError("no request for $code"))
            station.meta.request = requests[request_ind]
        end
    end
    stations
end

"""
    parse_station_response(T, request, response, server) -> stations::Vector{Seis.Station}

Convert a `response` to an `FDSNStation` request obtained with
[`SeisRequests.post_request`](@ref) or [`SeisRequests.get_request`](@ref)
into a set of stations.

`T` is the element type of the `Vector{Seis.GeogStation{T}}` returned
and defaults to `Float64` if this function is called from
[`get_stations`](@ref).
"""
function parse_station_response(T, request::FDSNStation, response, server)
    str = String(response.body)
    content_type = HTTP.Messages.header(response, "Content-Type")
    station_type = GeogStation{T}
    # No results
    if response.status == request.nodata
        isempty(str) || @warn("unexpected data in response reporting no data")
        return station_type[]
    elseif isempty(str)
        error("empty response but server did not indicate no data")
    end
    # If asked for a text format, parse that internally
    if coalesce(request.format, "") == "text"
        _check_content_type(response, "text/plain")
        lines = (line for line in split(str, '\n') if !occursin(r"^(\s*#.*|\s*)$", line))
        if coalesce(request.level, "") == "network"
            fdsn_networks = [convert(FDSNNetworkTextResponse, line) for line in lines]
            stations = parse_stations(station_type, request, fdsn_networks, server)
        # Default is station level
        elseif request.level === missing || request.level == "station"
            fdsn_stations = [convert(FDSNStationTextResponse, line) for line in lines]
            stations = parse_stations(station_type, request, fdsn_stations, server)
        elseif coalesce(request.level, "") == "channel"
            fdsn_channels = [convert(FDSNChannelTextResponse, line) for line in lines]
            stations = parse_stations(station_type, request, fdsn_channels, server)
        end
    # Otherwise we get XML and parse with StationXML
    elseif request.format === missing || request.format == "xml"
        _check_content_type(response, "application/xml")
        stationxml = StationXML.readstring(str)
        stations = parse_stations(station_type, request, stationxml, server)
    end
    stations
end

"""
    parse_stations(T, request, fdsn_text_responses, server)

Convert `fdsn_text_responses`, a set of [`SeisRequests.FDSNTextResponse`](@ref)s,
into a set of `Seis.Station`s.
`T` is the type of `Seis.Station` returned.

This method is usually called by [`SeisRequests.parse_station_response`](@ref).
"""
function parse_stations(T, request, fdsn_networks::AbstractArray{FDSNNetworkTextResponse},
        server)
    stations = similar(Vector{T}, axes(fdsn_networks))
    for (i, n) in enumerate(fdsn_networks)
        stations[i] = T(
            net=n.network,
            meta=Dict{Symbol,Any}(
                :network_description=>n.description,
                :network_number_stations=>n.total_stations,
                :startdate=>n.starttime,
                :server=>server
                )
            )
        # Avoid a key with an actual missing value in station[i]'s meta SeisDict
        stations[i].meta.enddate = n.endtime
    end
    stations
end

function parse_stations(T, request, fdsn_stations::AbstractArray{FDSNStationTextResponse},
        server)
    stations = similar(Vector{T}, axes(fdsn_stations))
    for (i, s) in enumerate(fdsn_stations)
        stations[i] = T(;
            net=s.network,
            sta=s.station,
            lon=s.longitude,
            lat=s.latitude,
            elev=s.elevation,
            meta=Dict{Symbol,Any}(
                :sitename=>s.sitename,
                :startdate=>s.starttime,
                :server=>server
                )
             )
        # Avoid a key with an actual missing value in station[i]'s meta SeisDict
        stations[i].meta.enddate = s.endtime
    end
    stations
end

function parse_stations(T, request, fdsn_channels::AbstractArray{FDSNChannelTextResponse},
        server)
    stations = similar(Vector{T}, axes(fdsn_channels))
    for (i, c) in enumerate(fdsn_channels)
        stations[i] = T(;
            net=c.network,
            sta=c.station,
            loc=c.location,
            cha=c.channel,
            lon=c.longitude,
            lat=c.latitude,
            elev=c.elevation,
            dep=c.depth/1000,
            azi=c.azimuth,
            inc=c.dip+90,
            meta=Dict{Symbol,Any}(
                :sensor_description=>c.sensor_description,
                :scale=>c.scale,
                :scale_frequency=>c.scale_frequency,
                :scale_units=>c.scale_units,
                :sample_rate=>c.sample_rate,
                :startdate=>c.starttime,
                :server=>server
                )
            )
        # Avoid a key with an actual missing value in station[i]'s meta SeisDict
        stations[i].meta.enddate = c.endtime
    end
    stations
end

"""
    parse_stations(T, request, stationxml::FDSNStationXML, server)

Create stations from an FDSNStationXML object.

Several fields in each station's `.meta` field are filled, including
the `stationxml` field which contains a full `FDSNStationXML` object
for that network/channel/station (depending on the request `level`), but
only for the network/station/channel in question.

This method is usually called by [`SeisRequests.parse_station_response`](@ref).
"""
function parse_stations(T, request, stationxml::FDSNStationXML, server)
    stations = T[]
    for network in stationxml.network
        if request.level !== missing && request.level == "network"
            net = T(net=network.code)
            net.meta.startdate = network.start_date
            net.meta.enddate = network.end_date
            net.meta.stationxml = filter_stationxml(stationxml, network)
            net.meta.server = server
            net.meta.request = request
            push!(stations, net)
            continue
        end

        for station in network.station
            if request.level === missing || request.level == "station"
                sta = T(net=network.code, sta=station.code,
                    lon=station.longitude, lat=station.latitude,
                    elev=station.elevation)
                sta.meta.startdate = station.start_date
                sta.meta.enddate = station.end_date
                sta.meta.stationxml = filter_stationxml(stationxml, network, station)
                sta.meta.server = server
                sta.meta.request = request
                push!(stations, sta)
                continue
            end

            for channel in station.channel
                # Fill in information which may need adjusting first
                cha = T(net=network.code, sta=station.code,
                    loc=channel.location_code, cha=channel.code,
                    lon=channel.longitude, lat=channel.latitude,
                    elev=channel.elevation, dep=channel.depth,
                    azi=channel.azimuth, inc=channel.dip)
                # Then adjust information by converting to correct units, etc.
                cha.dep = cha.dep/1000
                cha.inc = cha.inc + 90
                cha.meta.startdate = channel.start_date
                cha.meta.enddate = channel.end_date
                cha.meta.stationxml = filter_stationxml(stationxml, network, station, channel)
                cha.meta.server = server
                cha.meta.request = request
                push!(stations, cha)
            end
        end
    end
    stations
end

"""
    filter_stationxml(sxml, network, station=nothing, channel=nothing) -> ::FDSNStationXML

Take a `FDSNStationXML` document and filter out all networks, stations
and channels apart from `network`, `station` and `channel`.
The returned copy only contains one network if the request level is `"network"`,
and maybe one station in that network if the request level is `"station"`,
and maybe one channel at that station if the request level is `"channel"`.
"""
function filter_stationxml(sxml, network, station=nothing, channel=nothing)
    filtered_sxml = deepcopy(sxml)
    empty!(filtered_sxml.network)
    push!(filtered_sxml.network, deepcopy(network))
    empty!(filtered_sxml.network[1].station)
    if station !== nothing
        push!(filtered_sxml.network[1].station, deepcopy(station))
        empty!(filtered_sxml.network[1].station[1].channel)
        if channel !== nothing
            push!(filtered_sxml.network[1].station[1].channel, channel)
        end
    end
    filtered_sxml
end

#
# get_data
#

"""
    get_data(; query_type=FDSNDataSelect, server=$(DEFAULT_SERVER), verbose=true, T=Float64, kwargs...)

Get time series data from `server`.

If `verbose` is `true` (the default), information is printed to the screen on
the success or otherwise of the request.

`T` specifies the float type of the seismic data returned and defaults to `Float64`.

In this form, a data request is created by passing `kwargs` to `query_type`,
by default creating a [`FDSNDataSelect`](@ref) request.

# Keyword arguments
- `query_type = FDSNDataSelect`: Type of request to send.
- `server = $(DEFAULT_SERVER)`: Server name or URI to which to send the request.
- `verbose = true`: If `true`, print information about a query's progress to
  the terminal.
- `T = Float64`: Specify the float type to use for the data returned.
- `kwargs...`: Remaining keyword arguments, other than the above, are passed
  to `query_type` (i.e., by default to [`FDSNDataSelect`](@ref)) to set the
  search parameters for data.

# Examples
Get 10 minutes of data for the high-gain vertical channels at station JSA
(Jersey, Channel Islands) from the beginning of the year 2020.
This sends the request to `"$DEFAULT_SERVER"` by default.
```
julia> get_data(network="GB", station="JSA", location="", channel="?HZ", starttime="2020-01-01", endtime="2020-01-01T00:10:00")
[ Info: Request status: Successful request, results follow
3-element Array{Trace{Float64,Array{Float64,1},Seis.Geographic{Float64}},1}:
 Seis.Trace(GB.JSA..BHZ: delta=0.02, b=0.0, nsamples=29737)
 Seis.Trace(GB.JSA..BHZ: delta=0.02, b=595.6801, nsamples=216)
 Seis.Trace(GB.JSA..HHZ: delta=0.01, b=0.0, nsamples=60000)
```

Get the same data but from the Orfeus server.
```
julia> using Dates

julia> start = DateTime(2020)
2020-01-01T00:00:00

julia> get_data(server="Orfeus", T=Float32, network="GB", station="JSA", location="", channel="?HZ", starttime=start, endtime=start+Minute(10))
1-element Array{Trace{Float32,Array{Float32,1},Seis.Geographic{Float32}},1}:
 Seis.Trace(GB.JSA..HHZ: delta=0.01, b=0.0, nsamples=60921)
```
"""
function get_data(; query_type=FDSNDataSelect, server=DEFAULT_SERVER, verbose=true, T=Float64, kwargs...)
    request = query_type(; kwargs...)
    post_data_request_and_parse_response(T, [request], server, verbose)
end

"""
    get_data(stations, starttime, endtime; kwargs...)

Query for data defined at a set of `Station`s `stations`.

In this form, the values for `network`, `station`, `location` and `channel`
are taken from each station in `stations`.  Keyword arguments are as for
the zero-argument method.
"""
function get_data(stations::AbstractArray{<:Seis.Station}, starttime, endtime;
        query_type=FDSNDataSelect, server=DEFAULT_SERVER, verbose=true,
        T=Float64, kwargs...)
    requests = [query_type(; (network=sta.net,
                              station=sta.sta,
                              location=sta.loc,
                              channel=sta.cha,
                              starttime=starttime,
                              endtime=endtime,
                              kwargs...)...)
                for sta in stations]
    post_data_request_and_parse_response(T, requests, server, verbose; stations=stations)
end

"""
    get_data(event, stations, start_offset, end_offset; kwargs...)

Get time series data at `stations` relative to `event`.

`event` is an `Event` whose origin time defines the time window to request.
The window is defined as `event.time + start_offset` to `event.time + end_offset`.
`start_offset` and `end_offset` can be given in real numbers of seconds, or
as a `Dates.Period`.

`stations` is a set of `Station`s whose channel codes determine which data
to request.

The returned `Trace`s have their `.sta` field set to the appropriate element
of `stations` and their `.evt` field is `event.

This method is particularly useful because Miniseed data (returned by default
from datacentres) does not contain any station metadata such as coordinates.

Keyword arguments are as for the no-argument method of `get_data`.

# Examples
Get a minute of data from stations HYB and SCZ in the G network, for the
high-gain broadand vertical channel only, at location `""`.
```
julia> stations = get_stations(code="G.HYB,SCZ..BHZ", verbose=false)
2-element Array{Station{Float64,Seis.Geographic{Float64}},1}:
 Station: G.HYB..BHZ, lon: 78.553, lat: 17.417, dep: 0.0, elev: 510.0, azi: 0.0, inc: 0.0, meta: 4 keys
 Station: G.SCZ..BHZ, lon: -121.403, lat: 36.598, dep: 0.0, elev: 261.0, azi: 0.0, inc: 0.0, meta: 5 keys

julia> traces = get_data(stations, DateTime(2001), DateTime(2001) + Minute(1))
[ Info: Request status: Successful request, results follow
2-element Array{Trace{Float64,Array{Float64,1},Seis.Geographic{Float64}},1}:
 Seis.Trace(G.HYB..BHZ: delta=0.05, b=0.0, nsamples=1200)
 Seis.Trace(G.SCZ..BHZ: delta=0.05, b=0.0, nsamples=1200)

julia> all(traces.sta .== stations)
true

julia> traces.sta.lon, traces.sta.lat
([78.553, -121.403], [17.417, 36.598])
```
"""
function get_data(event::Seis.Event, stations::AbstractArray{<:Seis.GeogStation},
        start_offset::_Period, end_offset::_Period;
        query_type=FDSNDataSelect, server=DEFAULT_SERVER, verbose=true,
        T=Float64, kwargs...)
    requests = [query_type(; (network=sta.net,
                              station=sta.sta,
                              location=sta.loc,
                              channel=sta.cha,
                              starttime=event.time + start_offset,
                              endtime=event.time + end_offset,
                              kwargs...)...)
                for sta in stations]
    post_data_request_and_parse_response(T, requests, server, verbose, event=event,
        stations=stations)
end

function get_data(event, stations, start_offset::Real, end_offset::Real; kwargs...)
    get_data(event, stations, seconds_milliseconds(start_offset),
        seconds_milliseconds(end_offset); kwargs...)
end

"""
    get_data(event, stations, start_phase, start_offset, end_phase, end_offset; kwargs...)
    get_data(event, stations, phase, start_offset, end_offset; kwargs...)

Define the data window requested by reference to a predicted seismic phase
arrival.

In the first form, specify `start_phase` and `start_offset`, and `end_phase`
and `end_offset`.  The start of the data window is the event origin time,
plus the predicted travel time for `start_phase`, plus `start_offset`.
The end time limit is defined similarly.

In the second form, both window start and end are relative to a single `phase`.

# Keyword arguments
In addition to the keyword arguments accepted for all `get_data` methods,
the following can be used with this method:

- `exact = false`: If `true`, only use phase names which are an exact match
  for `phase`, `start_phase` and `end_phase`.
- `model = "iasp91"`: Choose which model to use for travel time calculations.
  Available models are currently: `("1066a", "1066b", "ak135", "ak135f_no_mud",
  "herrin", "iasp91", "jb", "prem", "pwdk", "sp6")`

# Notes
!!! note
    To use this form of `get_data`, you must first have installed the
    [SeisTau.jl](https://github.com/anowacki/SeisTau.jl) package, and
    then done `using SeisTau` or `import SeisTau` in your session.

!!! note
    If calling this form of `get_data` from another module in which you
    have loaded `SeisTau`, pass the `SeisTau` module object to this
    method's `seistau` keyword argument, like so:
    ```
    module MyModule
    using SeisTau
    using SeisRequests

    function process(event, stations, phase, t0, t1)
        data = get_data(event, stations, phase, t0, t1; seistau=SeisTau)
        ...
    end
    end # module
    ```
"""
function get_data(event, stations, start_phase, start_offset::_Period, end_phase, end_offset::_Period;
        server=DEFAULT_SERVER, verbose=true, T=Float64,
        query_type=FDSNDataSelect,
        # SeisTau options
        seistau=nothing, exact=false, model="iasp91",
        # Query options
        kwargs...)
    # FIXME: Horrible hack to permit the use of SeisTau without actually
    #        depending on it.
    # seistau is the module SeisTau
    if seistau === nothing
        !isdefined(Main, :SeisTau) &&
            error("Using phase arrival times requires loading SeisTau.  " *
                  "Do `using SeisTau` or `import SeisTau` and try again.")
        seistau = getfield(Main, :SeisTau)
    end
    requests = query_type[]
    for sta in stations
        start_arrivals = seistau.travel_time(event, sta, start_phase; exact=exact, model=model)
        isempty(start_arrivals) &&
            throw(ArgumentError("no arrivals for start phase \"$start_phase\""))
        end_arrivals = start_phase == end_phase ? start_arrivals :
            seistau.travel_time(event, sta, end_phase; exact=exact, model=model)
        isempty(end_arrivals) &&
            throw(ArgumentError("no arrivals for end phase \"$end_phase\""))
        starttime = event.time + seconds_milliseconds(first(start_arrivals).time) +
            start_offset
        endtime = event.time + seconds_milliseconds(first(end_arrivals).time) +
            end_offset
        request = query_type(; (network=sta.net, station=sta.sta, location=sta.loc,
                    channel=sta.cha, starttime=starttime, endtime=endtime, kwargs...)...)
        push!(requests, request)
    end
    traces = post_data_request_and_parse_response(T, requests, server, verbose; event=event,
        stations=stations)
    traces
end

# Accept Real arguments and convert to DateTimes
get_data(event, stations, start_phase, start_offset::Real, end_phase, end_offset::Real; kwargs...) =
    get_data(event, stations, start_phase, seconds_milliseconds(start_offset),
        end_phase, seconds_milliseconds(end_offset); kwargs...)

# Use the same phase for start and end window time
get_data(event, stations, phase, start_offset, end_offset; kwargs...) =
    get_data(event, stations, phase, start_offset, phase, end_offset; kwargs...)

"""
    parse_data_response(T, requests, response, server; event=nothing, stations=nothing)

Convert the `response` from a data `server`, from a `request`, into
`Trace`s.  In this form, the response is from a single event recorded
at multiple stations.

If `event` is given, then assume that all traces should have this
event as their origin and shift their origin times accordingly.

If `stations` (an array of `Seis.Station`s) is given, then assume that these
stations were requested in `request`
"""
function parse_data_response(T, requests::AbstractArray{<:SeisRequest},
        response, server;
        event=nothing, stations=nothing)
    if response.status == first(requests).nodata
        isempty(response.body) || @warn("unexpected data in response reporting no data")
        return Seis.Trace{T, Vector{T}, Seis.Geographic{T}}[]
    elseif isempty(response.body)
        error("empty response but server did not indicate no data")
    end
    format = first(requests).format
    if format === missing || format == "miniseed"
        _check_content_type(response, "application/vnd.fdsn.mseed")
        traces = Seis.read_mseed(Seis.Trace{T, Vector{T}, Seis.Geographic{T}},
            response.body)
    else
        error("format \"$(format)\" is currently unsupported")
    end
    traces.meta.server = server
    if event !== nothing
        Seis.origin_time!.(traces, event.time)
        traces.evt = event
    end
    if stations !== nothing
        for station in stations
            code = Seis.channel_code(station)
            inds = findall(x -> Seis.channel_code(x) == code, traces)
            for i in inds
                new_meta = merge!(traces[i].sta.meta, station.meta)
                traces[i].sta = station
                traces[i].sta.meta = new_meta
            end
        end

    end
    traces
end

parse_data_response(T, request::SeisRequest, response, server; kwargs...) =
    parse_data_response(T, [request], response, server; kwargs...)

"""
    post_data_request_and_parse_response(T, requests, server, verbose; event=nothing, stations=nothing) -> traces

Send the vector of `SeisRequest`s `requests` to `server` via the HTTP POST method.
If `verbose` is `true`, then print information about the success of the request.

Optionally pass information about an `event` and a set of `stations` (e.g.,
from a previous call to [`get_event`](@ref) or [`get_stations`](@ref)) to fill
in the `traces`' information with this additional data.
"""
function post_data_request_and_parse_response(T, requests::AbstractArray{<:SeisRequest},
        server, verbose; event=nothing, stations=nothing)
    response = post_request(requests; server=server, verbose=verbose)
    traces = parse_data_response(T, requests, response, server, event=event,
        stations=stations)
    traces
end

#
# get_events
#

"""
    get_events(; server=$(DEFAULT_SERVER), verbose=true, T=Float64, kwargs...) -> events

Query `server` for events matching those specified in `kwargs`.
`events` is a `Vector{Seis.GeogEvent{T}}`, where various useful
fields are added to the `.meta` field of each event.

For search options, see [`FDSNEvent`](@ref).

The element type of the `Seis.Event`s returned is set by `T`.
"""
function get_events(; server=DEFAULT_SERVER, verbose=true, T=Float64, kwargs...)
    request = FDSNEvent(; kwargs...)
    response = get_request(request; server=server, verbose=verbose)
    events = parse_event_response(T, request, response, server)
end

"""
    get_events(station; longitude=station.lon, latitude=station.lat,
        starttime=station.meta.startdate, endtime=station.meta.enddate, kwargs...)

Search for events in relation to a `station`, typically constraining the
event locations to be a certain distance from `event`, or requiring events
to occur when the station was active.

Other keyword arguments `kwargs` are passed to the zero-argument method above.

# Example
Find all events between 90° and 120° from station TA.112A, of mangitude 6 and
above, greater than 100 km deep, which occurred whilst the station was active:
```
julia> sta = first(get_stations(network="TA", station="112A"));

julia> get_events(sta, minradius=90, maxradius=120, mindepth=100, minmagnitude=6)
[ Info: Request status: Successful request, results follow
5-element Array{Seis.Event{Float64,Seis.Geographic{Float64}},1}:
 Event: lon: 130.3051, lat: -5.9259, dep: 188.5, time: 2008-08-04T20:45:15.82, id: smi:service.iris.edu/fdsnws/event/1/query?originid=4672548, meta: 7 keys
 Event: lon: 127.8966, lat: -7.5592, dep: 121.7, time: 2008-06-06T13:42:49.06, id: smi:service.iris.edu/fdsnws/event/1/query?originid=4620129, meta: 7 keys
 Event: lon: -28.0966, lat: -56.0769, dep: 112.3, time: 2008-04-14T09:45:16.97, id: smi:service.iris.edu/fdsnws/event/1/query?originid=4577607, meta: 7 keys
 Event: lon: 127.5146, lat: -7.6192, dep: 181.0, time: 2007-12-15T08:03:16.44, id: smi:service.iris.edu/fdsnws/event/1/query?originid=4494689, meta: 7 keys
 Event: lon: 151.8587, lat: -4.5977, dep: 134.3, time: 2007-05-29T01:03:28.32, id: smi:service.iris.edu/fdsnws/event/1/query?originid=4359323, meta: 7 keys
```
"""
get_events(station::Seis.GeogStation; kwargs...) =
    get_events(; (longitude=station.lon, latitude=station.lat,
        starttime=station.meta.startdate, endtime=station.meta.enddate, kwargs...)...)

"""
    parse_event_response(T, request, response, server) -> events

Parse the output of a `get_request` or `post_request` call and return a
set of `events` containing the results.

`request` is a `FDSNEvent` request, `response` is a `HTTP.Response` with
the server's response to the request, and `server` specifies which server
the request was sent to.
"""
function parse_event_response(T::DataType, request::FDSNEvent, response, server)
    str = String(response.body)
    event_type = GeogEvent{T}
    if response.status == request.nodata
        isempty(str) || @warn("unexpected data in response reporting no data")
        return event_type[]
    elseif isempty(str)
        error("empty response but server did not indicate no data")
    end
    # If have requested text format, then parse that internally
    if coalesce(request.format, "") == "text"
        _check_content_type(response, "text/plain")
        fdsn_events = [convert(FDSNEventTextResponse, line)
                       for line in split(str, '\n') if !occursin(r"^(\s*#.*|\s*)$", line)]
        events = parse_events(event_type, request, fdsn_events, server)
    # ISF format from the ISC
    elseif coalesce(request.format, "") == "isf"
        error("parsing of ISF responses not yet implements")
    # Otherwise in QuakeML format (XML)
    else
        _check_content_type(response, "application/xml")
        quakeml = QuakeML.readstring(str)
        events = parse_events(event_type, request, quakeml, server)
    end
    events
end

parse_event_response(request, response, server) =
    parse_event_response(Float64, request, response, server)

"""
    parse_events(T, request, fdsn_events) -> events

Convert a set of `fdsn_events` which have been parsed from a text-format
`FDSNEvent` request into a `Vector` of `Seis.Event`s.

`T` is the specific type of `Seis.Event` (e.g.,
`Seis.Event{Float32, Seis.Geographic{Float32}}`)
"""
function parse_events(T, request, fdsn_events::AbstractArray{FDSNEventTextResponse}, server)
    n = length(fdsn_events)
    events = Vector{T}(undef, n)
    for (i, e) in enumerate(fdsn_events)
        events[i] = T(;
            id=e.event_id,
            time=e.time,
            lon=e.longitude,
            lat=e.latitude,
            dep=e.depth,
            meta=Dict{Symbol,Any}(
                :author=>e.author,
                :catalog=>e.catalog,
                :contributor=>e.contributor,
                :contributor_id=>e.contributor_id,
                :mag_type=>e.mag_type,
                :mag=>e.magnitude,
                :mag_author=>e.mag_author,
                :location_name=>e.event_location_name,
                :server=>server
                )
            )
            # Possibly missing fields
            meta = events[i].meta
            meta.type = e.event_type
    end
    events
end

"""
    parse_events(T, request, quakeml_events, server) -> events

Convert a `QuakeML.EventParameters` object into a set of `T`s, where
`T <: Seis.GeogEvent`.
"""
function parse_events(T, request, quakeml_events::QuakeML.EventParameters, server)
    n = length(quakeml_events.event)
    events = Vector{T}(undef, n)
    for (i, event) in enumerate(quakeml_events.event)
        if !QuakeML.has_origin(event)
            @warn("No origin associated with event $(event.public_id)")
            events[i] = T()
        else
            origin = QuakeML.preferred_origin(event)
            events[i] = T(;
                id=origin.public_id.value,
                time=origin.time.value,
                lon=origin.longitude.value,
                lat=origin.latitude.value,
                dep=_getifnotmissing(origin.depth, :value)/1000
                )
            meta = events[i].meta
            meta.author = _getifnotmissing(origin.creation_info, :author)
            meta.catalog = request.catalog
            meta.type = _getifnotmissing(event.type, :value)
            meta.origin_type = _getifnotmissing(origin.type, :value)
            meta.server = server
            meta.quakeml = event
            if QuakeML.has_focal_mechanism(event)
                focal_mechanism = QuakeML.preferred_focal_mechanism(event)
                meta.focal_mechanism = focal_mechanism
            end
            if QuakeML.has_magnitude(event)
                magnitude = QuakeML.preferred_magnitude(event)
                meta.mag = magnitude.mag.value
                meta.mag_type = magnitude.type
                if magnitude.creation_info !== missing
                    meta.mag_author = magnitude.creation_info.author
                end
            end
            if !isempty(event.description)
                desc = event.description[1]
                meta.description = desc.text *
                    (desc.type === missing ? "" : " ("*desc.type.value*")")
            end
        end
    end
    events
end

"Return `missing` if `val.field` is `missing`, and its value otherwise."
_getifnotmissing(val, field) = val === missing ? missing : getfield(val, field)

"""
    _check_content_type(response, expected; allowempty=true) -> nothing

Log a warning if the `"Content-Type"` header of `response` is not `expected`.

If `allowempty` is `true` (the default) then empty messages do
not throw.
"""
function _check_content_type(response, expected; allowempty=true)
    content_type = HTTP.Messages.header(response, "Content-Type")
    allowempty && isempty(content_type) && return
    content_type != expected &&
        @warn("content type of response is \"$content_type\", " *
              "not \"$expected\" as expected")
    return
end


text:
@testset "224.basic-calculator.jl" begin
    @test calculate("1 + 1") == 2
    @test calculate(" 2-1 + 2 ") == 3
    @test calculate(" 2-1 + 2 ") == 3
    @test calculate("(1+(4+5+2)-3)+(6+8)") == 23
end


text:
using Stheno: GPC, ZeroMean, ConstMean, CustomMean, ZeroKernel
using Stheno: EQ, Exp

@timedtestset "gp" begin

    # Ensure that basic functionality works as expected.
    @timedtestset "GP" begin
        rng, gpc, N, N′ = MersenneTwister(123456), GPC(), 5, 6
        m, k = CustomMean(sin), EQ()
        f = GP(m, k, gpc)
        x = collect(range(-1.0, 1.0; length=N))
        x′ = collect(range(-1.0, 1.0; length=N′))

        @test mean_vector(f, x) == ew(m, x)
        @test cov(f, x) == pw(k, x)
        @test cov_diag(f, x) == diag(cov(f, x))
        @test cov(f, x, x) == pw(k, x, x)
        @test cov(f, x, x′) == pw(k, x, x′)
        @test cov(f, x, x′) ≈ cov(f, x′, x)'
    end

    # Check that mean-function specialisations work as expected.
    @timedtestset "sugar" begin
        @test GP(5, EQ(), GPC()).m isa ConstMean
        @test GP(EQ(), GPC()).m isa ZeroMean
    end

    # Check that GP(kernel, mean, gpc) always works.
    @timedtestset "reversed construction" begin
        @test GP(eq(), CustomMean(sin), GPC()).m isa CustomMean{typeof(sin)}
        @test GP(eq(), CustomMean(sin), GPC()).k isa Stheno.EQ
        @test GP(eq(), 5.0, GPC()).m isa ConstMean
        @test GP(eq(), 5.0, GPC()).k isa Stheno.EQ
    end

    # Test the creation of indepenent GPs.
    @timedtestset "independent GPs" begin
        rng, N, N′ = MersenneTwister(123456), 5, 6
        x, x′ = randn(rng, N), randn(rng, N′)

        # Specification of two independent GPs
        gpc = GPC()
        m1, m2 = ZeroMean(), ConstMean(5)
        k1, k2 = EQ(), EQ()
        f1, f2 = GP(m1, k1, gpc), GP(m2, k2, gpc)

        @test mean_vector(f1, x) == ew(m1, x)
        @test mean_vector(f2, x) == ew(m2, x)

        @test cov(f1, f2, x, x′) == zeros(N, N′)
        @test cov_diag(f1, x) == ones(N)

        @test cov(f1, f1, x′, x) ≈ cov(f1, f1, x, x′)'
    end
end


text:
@testset "intrinsiccoregion" begin
    rng = MersenneTwister(123)

    dims = (in=3, out=2, obs=3)
    r = 1

    A = randn(dims.out, r)
    B = A * transpose(A) + Diagonal(rand(dims.out))

    # XIF = [(rand(dims.in), rand(1:(dims.out))) for i in 1:(dims.obs)]
    x = [rand(dims.in) for _ in 1:2]
    XIF = KernelFunctions.MOInputIsotopicByFeatures(x, dims.out)
    XIO = KernelFunctions.MOInputIsotopicByOutputs(x, dims.out)
    y = [rand(dims.in) for _ in 1:2]
    YIF = KernelFunctions.MOInputIsotopicByFeatures(y, dims.out)
    YIO = KernelFunctions.MOInputIsotopicByOutputs(y, dims.out)
    z = [rand(dims.in) for _ in 1:3]
    ZIF = KernelFunctions.MOInputIsotopicByFeatures(z, dims.out)
    ZIO = KernelFunctions.MOInputIsotopicByOutputs(z, dims.out)

    kernel = SqExponentialKernel()
    icoregionkernel = IntrinsicCoregionMOKernel(kernel, B)

    icoregionkernel2 = IntrinsicCoregionMOKernel(; kernel=kernel, B=B)
    @test icoregionkernel == icoregionkernel2

    @test icoregionkernel.B == B
    @test icoregionkernel.kernel == kernel
    @test icoregionkernel(XIF[1], XIF[1]) ≈
          B[XIF[1][2], XIF[1][2]] * kernel(XIF[1][1], XIF[1][1])
    @test icoregionkernel(XIF[1], XIF[end]) ≈
          B[XIF[1][2], XIF[end][2]] * kernel(XIF[1][1], XIF[end][1])

    # kernelmatrix
    KernelFunctions.TestUtils.test_interface(icoregionkernel, XIF, YIF, ZIF)

    KernelFunctions.TestUtils.test_interface(icoregionkernel, XIO, YIO, ZIO)

    KernelFunctions.TestUtils.test_interface(
        icoregionkernel, Vector{Tuple{Float64,Int}}; dim_out=dims.out
    )

    test_ADs(icoregionkernel; dims=dims)

    @test string(icoregionkernel) ==
          string("Intrinsic Coregion Kernel: ", kernel, " with ", dims.out, " outputs")
end


text:
using Pkg

Pkg.activate(@__DIR__)
Pkg.instantiate()

using JSON, Query

function main()
    db_x64 = Dict{String,Any}("AvailableVersions"=>Dict{String,Any}(), "AvailableChannels"=>Dict{String,Any}())
    db_x86 = Dict{String,Any}("AvailableVersions"=>Dict{String,Any}(), "AvailableChannels"=>Dict{String,Any}())

    old_db = JSON.parsefile(joinpath(@__DIR__, "..", "versions.json"), use_mmap=false)

    versions = map(old_db["OptionalJuliaPackages"]) do i
        v = VersionNumber(i["JuliaVersion"])
        return (version=v, major="$(v.major)", minor="$(v.major).$(v.minor)")
    end

    minor_channels = versions |> @groupby(_.minor) |> @map({channel=key(_), version=first(_ |> @orderby_descending(i->i.version)).version})
    major_channels = versions |> @groupby(_.major) |> @map({channel=key(_), version=first(_ |> @orderby_descending(i->i.version)).version})

    for v in versions
        db_x64["AvailableVersions"]["$(v.version)+0~x64"] = Dict("Url"=>"https://julialang-s3.julialang.org/bin/winnt/x64/$(v.minor)/julia-$(v.version)-win64.tar.gz")
        db_x64["AvailableVersions"]["$(v.version)+0~x86"] = Dict("Url"=>"https://julialang-s3.julialang.org/bin/winnt/x86/$(v.minor)/julia-$(v.version)-win32.tar.gz")
        db_x86["AvailableVersions"]["$(v.version)+0~x86"] = Dict("Url"=>"https://julialang-s3.julialang.org/bin/winnt/x86/$(v.minor)/julia-$(v.version)-win32.tar.gz")

        db_x64["AvailableChannels"]["$(v.version)~x64"] = Dict("Version" => "$(v.version)+0~x64")
        db_x64["AvailableChannels"]["$(v.version)~x86"] = Dict("Version" => "$(v.version)+0~x86")
        db_x86["AvailableChannels"]["$(v.version)~x86"] = Dict("Version" => "$(v.version)+0~x86")
        db_x64["AvailableChannels"]["$(v.version)"] = Dict("Version" => "$(v.version)+0~x64")
        db_x86["AvailableChannels"]["$(v.version)"] = Dict("Version" => "$(v.version)+0~x86")
    end

    db_x64["AvailableVersions"]["1.7.0-beta2+0~x64"] = Dict("Url"=>"https://julialang-s3.julialang.org/bin/winnt/x64/1.7/julia-1.7.0-beta2-win64.tar.gz")
    db_x64["AvailableVersions"]["1.7.0-beta2+0~x86"] = Dict("Url"=>"https://julialang-s3.julialang.org/bin/winnt/x86/1.7/julia-1.7.0-beta2-win32.tar.gz")
    db_x86["AvailableVersions"]["1.7.0-beta2+0~x86"] = Dict("Url"=>"https://julialang-s3.julialang.org/bin/winnt/x86/1.7/julia-1.7.0-beta2-win32.tar.gz")

    db_x64["AvailableVersions"]["1.7.0-beta3+0~x64"] = Dict("Url"=>"https://julialang-s3.julialang.org/bin/winnt/x64/1.7/julia-1.7.0-beta3-win64.tar.gz")
    db_x64["AvailableVersions"]["1.7.0-beta3+0~x86"] = Dict("Url"=>"https://julialang-s3.julialang.org/bin/winnt/x86/1.7/julia-1.7.0-beta3-win32.tar.gz")
    db_x86["AvailableVersions"]["1.7.0-beta3+0~x86"] = Dict("Url"=>"https://julialang-s3.julialang.org/bin/winnt/x86/1.7/julia-1.7.0-beta3-win32.tar.gz")

    for c in minor_channels
        db_x64["AvailableChannels"]["$(c.channel)~x64"] = Dict("Version" => "$(c.version)+0~x64")
        db_x64["AvailableChannels"]["$(c.channel)~x86"] = Dict("Version" => "$(c.version)+0~x86")
        db_x86["AvailableChannels"]["$(c.channel)~x86"] = Dict("Version" => "$(c.version)+0~x86")
        db_x64["AvailableChannels"]["$(c.channel)"] = Dict("Version" => "$(c.version)+0~x64")
        db_x86["AvailableChannels"]["$(c.channel)"] = Dict("Version" => "$(c.version)+0~x86")
    end

    for c in major_channels
        db_x64["AvailableChannels"]["$(c.channel)~x64"] = Dict("Version" => "$(c.version)+0~x64")
        db_x64["AvailableChannels"]["$(c.channel)~x86"] = Dict("Version" => "$(c.version)+0~x86")
        db_x86["AvailableChannels"]["$(c.channel)~x86"] = Dict("Version" => "$(c.version)+0~x86")
        db_x64["AvailableChannels"]["$(c.channel)"] = Dict("Version" => "$(c.version)+0~x64")
        db_x86["AvailableChannels"]["$(c.channel)"] = Dict("Version" => "$(c.version)+0~x86")
    end

    db_x64["AvailableChannels"]["release~x64"] = Dict("Version" => "1.6.2+0~x64")
    db_x64["AvailableChannels"]["release~x86"] = Dict("Version" => "1.6.2+0~x86")
    db_x86["AvailableChannels"]["release~x86"] = Dict("Version" => "1.6.2+0~x86")
    db_x64["AvailableChannels"]["release"] = Dict("Version" => "1.6.2+0~x64")
    db_x86["AvailableChannels"]["release"] = Dict("Version" => "1.6.2+0~x86")

    db_x64["AvailableChannels"]["lts~x64"] = Dict("Version" => "1.0.5+0~x64")
    db_x64["AvailableChannels"]["lts~x86"] = Dict("Version" => "1.0.5+0~x86")
    db_x86["AvailableChannels"]["lts~x86"] = Dict("Version" => "1.0.5+0~x86")
    db_x64["AvailableChannels"]["lts"] = Dict("Version" => "1.0.5+0~x64")
    db_x86["AvailableChannels"]["lts"] = Dict("Version" => "1.0.5+0~x86")

    db_x64["AvailableChannels"]["beta~x64"] = Dict("Version" => "1.7.0-beta3+0~x64")
    db_x64["AvailableChannels"]["beta~x86"] = Dict("Version" => "1.7.0-beta3+0~x86")
    db_x86["AvailableChannels"]["beta~x86"] = Dict("Version" => "1.7.0-beta3+0~x86")
    db_x64["AvailableChannels"]["beta"] = Dict("Version" => "1.7.0-beta3+0~x64")
    db_x86["AvailableChannels"]["beta"] = Dict("Version" => "1.7.0-beta3+0~x86")

    db_x64["AvailableChannels"]["rc~x64"] = Dict("Version" => "1.6.2+0~x64")
    db_x64["AvailableChannels"]["rc~x86"] = Dict("Version" => "1.6.2+0~x86")
    db_x86["AvailableChannels"]["rc~x86"] = Dict("Version" => "1.6.2+0~x86")
    db_x64["AvailableChannels"]["rc"] = Dict("Version" => "1.6.2+0~x64")
    db_x86["AvailableChannels"]["rc"] = Dict("Version" => "1.6.2+0~x86")

    mkpath(joinpath(@__DIR__, "..", "build", "versiondb"))

    open(joinpath(@__DIR__, "..", "build", "versiondb", "juliaup-versionsdb-winnt-x64.json"), "w") do f
        JSON.print(f, db_x64, 4)
    end

    open(joinpath(@__DIR__, "..", "build", "versiondb", "juliaup-versionsdb-winnt-arm64.json"), "w") do f
        JSON.print(f, db_x64, 4)
    end

    open(joinpath(@__DIR__, "..", "build", "versiondb", "juliaup-versionsdb-winnt-x86.json"), "w") do f
        JSON.print(f, db_x86, 4)
    end

    open(joinpath(@__DIR__, "..", "build", "versiondb", "juliaup-versionsdb-winnt-arm.json"), "w") do f
        JSON.print(f, db_x86, 4)
    end
end

main()

max_stars_repo_name: piever/juliaup

max_issues_repo_name: piever/juliaup

max_forks_repo_name: piever/juliaup