text:
__precompile__
struct TimeSolver
    alpha::Array{Float64,2};
    beta::Array{Float64,2};
    gamma::Array{Float64,1};
    rkStages::Int64;
    Nx::Int64;
    dt::Float64;

    KRK::Array{Float64,3};
    SRK::Array{Float64,3};
    LRK::Array{Float64,3};
    Krhs::Array{Float64,3};
    Srhs::Array{Float64,3};
    Lrhs::Array{Float64,3};

    # constructor
    function TimeSolver(settings::Settings)
        rkStages = settings.rkStages;
        Nx = settings.Nx;
        dt = settings.dt;

        # set up alpha and beta values for SSP RK method
        alpha, beta = SetUpRKTable(rkStages);

        # setup gamma according to https://gkyl.readthedocs.io/en/latest/dev/ssp-rk.html
        gamma = zeros(rkStages);
        if rkStages == 1
            gamma[1] = 0.0;
        elseif rkStages == 2
            gamma[1] = 0.0;
            gamma[2] = 1.0;
        elseif rkStages == 3
            gamma[1] = 0.0;
            gamma[2] = 1.0;
            gamma[3] = 0.5;
        end

        # allocate memory for matrices used in Update function
        KRK = zeros(rkStages+1,Nx,settings.r);
        SRK = zeros(rkStages+1,settings.r,settings.r);
        LRK = zeros(rkStages+1,settings.N^2,settings.r);
        Krhs = zeros(rkStages,Nx,settings.r);
        Srhs = zeros(rkStages,settings.r,settings.r);
        Lrhs = zeros(rkStages,settings.N^2,settings.r);

        new(alpha,beta,gamma,rkStages,Nx,dt,KRK,SRK,LRK,Krhs,Srhs,Lrhs)
    end
end

function SetUpRKTable(p::Int64)
    alpha = zeros(p,p);
    beta = zeros(p,p);
    if p == 1
        alpha[1,1] = 1.0;
        beta[1,1] = 1.0;
    elseif p == 2
        alpha[1,1] = 1.0;
        alpha[2,1] = 0.5;
        alpha[2,2] = 0.5;
        beta[1,1] = 1.0;
        beta[2,1] = 0.0;
        beta[2,2] = 0.5;
    elseif p == 3
        alpha[1,1] = 1.0;
        alpha[2,1] = 3.0/4.0;
        alpha[2,2] = 1.0/4.0;
        alpha[3,1] = 1.0/3.0;
        alpha[3,3] = 2.0/3.0;
        beta[1,1] = 1.0;
        beta[2,1] = 0.0;
        beta[2,2] = 1.0/4.0;
        beta[3,3] = 2.0/3.0;
    end
    return alpha, beta
end


text:
"""
    measure(ψ::MPS, measurement::Tuple{String,Int}, s::Vector{<:Index})
    measure(ψ::MPS, measurement::Tuple{String,Int})

Perform a measurement of a 1-local operator on an MPS ψ. The operator
is identifyed by a String (corresponding to a `gate`) and a site.
If an additional set of indices `s` is provided, the correct site is
extracted by comparing the MPS with the index order in `s`.
"""
function measure(ψ::MPS, measurement::Tuple{String,Int}, s::Vector{<:Index})
  site0 = measurement[2]
  site = findsite(ψ, s[site0])
  ϕ = orthogonalize!(copy(ψ), site)
  ϕs = ϕ[site]
  obs_op = gate(measurement[1], firstsiteind(ϕ, site))
  T = noprime(ϕs * obs_op)
  return real((dag(T) * ϕs)[])
end

measure(ψ::MPS, measurement::Tuple{String,Int}) = measure(ψ, measurement, siteinds(ψ))

function measure(ρ::MPO, args...)
  return error("Measurement of one-body operator on MPOs not yet implemented")
end
function measure(L::LPDO, args...)
  return error("Measurement of one-body operator on LPDOs not yet implemented")
end

"""
    measure(ψ::MPS, measurement::Tuple{String,Array{Int}}, s::Vector{<:Index})
    measure(ψ::MPS, measurement::Tuple{String,Array{Int}})
    measure(ψ::MPS, measurement::Tuple{String,AbstractRange}, s::Vector{<:Index})
    measure(ψ::MPS, measurement::Tuple{String,AbstractRange})
    measure(ψ::MPS, measurement::String, s::Vector{<:Index})
    measure(ψ::MPS, measurement::String)

Perform a measurement of a 1-local operator on an MPS ψ on a set of sites passed
as a vector. If an additional set of indices `s` is provided, the correct site is
extracted by comparing the MPS with the index order in `s`.
"""
function measure(ψ::MPS, measurement::Tuple{String,Array{Int}}, s::Vector{<:Index})
  result = []
  sites0 = measurement[2]
  ϕ = copy(ψ)
  for site0 in sites0
    site = findsite(ϕ, s[site0])
    orthogonalize!(ϕ, site)
    ϕs = ϕ[site]
    obs_op = gate(measurement[1], firstsiteind(ϕ, site))
    T = noprime(ϕs * obs_op)
    push!(result, real((dag(T) * ϕs)[]))
  end
  return result
end

function measure(ψ::MPS, measurement::Tuple{String,Array{Int}})
  return measure(ψ, measurement, siteinds(ψ))
end

# for a range of sites
function measure(ψ::MPS, measurement::Tuple{String,AbstractRange}, s::Vector{<:Index})
  return measure(ψ, (measurement[1], Array(measurement[2])), s)
end

function measure(ψ::MPS, measurement::Tuple{String,AbstractRange})
  return measure(ψ, (measurement[1], Array(measurement[2])), siteinds(ψ))
end

# for every sites
function measure(ψ::MPS, measurement::String, s::Vector{<:Index})
  return measure(ψ::MPS, (measurement, 1:length(ψ)), s)
end

## for every sites
function measure(ψ::MPS, measurement::String)
  return measure(ψ::MPS, (measurement, 1:length(ψ)), siteinds(ψ))
end

# at a given site
"""
    measure(ψ::MPS, measurement::Tuple{String,Int,String,Int}, s::Vector{<:Index})

Perform a measurement of a 2-body tensor-product operator on an MPS ψ. The two operators
are defined by Strings (for op name) and the sites. If an additional set of indices `s` is provided, the correct site is
extracted by comparing the MPS with the index order in `s`.
"""
function measure(ψ::MPS, measurement::Tuple{String,Int,String,Int}, s::Vector{<:Index})
  obsA = measurement[1]
  obsB = measurement[3]
  siteA0 = measurement[2]
  siteB0 = measurement[4]
  siteA = findsite(ψ, s[siteA0])
  siteB = findsite(ψ, s[siteB0])

  if siteA > siteB
    obsA, obsB = obsB, obsA
    siteA, siteB = siteB, siteA
  end
  ϕ = orthogonalize!(copy(ψ), siteA)
  ϕdag = prime(dag(ϕ); tags="Link")

  if siteA == siteB
    C = ϕ[siteA] * gate(obsA, firstsiteind(ϕ, siteA))
    C = noprime(C; tags="Site") * gate(obsA, firstsiteind(ϕ, siteA))
    C = noprime(C; tags="Site") * noprime(ϕdag[siteA])
    return real(C[])
  end
  if siteA == 1
    C = ϕ[siteA] * gate(obsA, firstsiteind(ϕ, siteA))
    C = noprime(C; tags="Site") * ϕdag[siteA]
  else
    C =
      prime(ϕ[siteA], commonind(ϕ[siteA], ϕ[siteA - 1])) *
      gate(obsA, firstsiteind(ϕ, siteA))
    C = noprime(C; tags="Site") * ϕdag[siteA]
  end
  for j in (siteA + 1):(siteB - 1)
    C = C * ϕ[j]
    C = C * ϕdag[j]
  end
  if siteB == length(ϕ)
    C = C * ϕ[siteB] * gate(obsB, firstsiteind(ϕ, siteB))
    C = noprime(C; tags="Site") * ϕdag[siteB]
  else
    C =
      C *
      prime(ϕ[siteB], commonind(ϕ[siteB], ϕ[siteB + 1])) *
      gate(obsB, firstsiteind(ϕ, siteB))
    C = noprime(C; tags="Site") * ϕdag[siteB]
  end
  return real(C[])
end

function measure(ψ::MPS, measurement::Tuple{String,Int,String,Int})
  return measure(ψ, measurement, siteinds(ψ))
end

function measure(ψ::MPS, measurement::Tuple{String,String}, s::Vector{<:Index})
  N = length(ψ)
  C = Matrix{Float64}(undef, N, N)
  for siteA in 1:N
    for siteB in 1:N
      m = (measurement[1], siteA, measurement[2], siteB)
      result = measure(ψ, m, s)
      C[siteA, siteB] = result
    end
  end
  return C
end

# for every sites
function measure(ψ::MPS, measurement::Tuple{String,String})
  return measure(ψ::MPS, measurement, siteinds(ψ))
end

"""
    entanglemententropy(ψ::MPS; bond = nothing)

Measure the entanglement entropy of an MPS `ψ` at `bond`.
"""
function entanglemententropy(ψ0::MPS; bond::Int=length(ψ0) ÷ 2)
  # make sure the state is normalized
  ψ = normalize!(copy(ψ0))

  # number of qubits
  N = length(ψ)
  @assert (bond < N)

  # gauge the MPS
  orthogonalize!(ψ, bond)

  # get singular values
  row_inds = (bond > 1 ? (linkind(ψ, bond - 1), siteind(ψ, bond)) : siteind(ψ, bond))
  u, s, v = svd(ψ[bond], row_inds)

  # Compute Von Neumann Entropy S = -Tr(ρ log(ρ))
  S = 0.0
  for n in 1:dim(s, 1)
    λ = s[n, n]^2
    S -= λ * log(λ + 1e-20)
  end
  return S
end

function entanglemententropy(ρ0::MPO; kwargs...)
  return error("Measurement of entanglement entropy for MPOs not yet implemented")
end
function entanglemententropy(ρ0::LPDO; kwargs...)
  return error("Measurement of entanglement entropy for LPDOs not yet implemented")
end

ext: jl


text:
struct Sobol <: GSAMethod
    order::Vector{Int}
    nboot::Int
    conf_level::Float64
end
Sobol(; order = [0, 1], nboot = 1, conf_level = 0.95) = Sobol(order, nboot, conf_level)

mutable struct SobolResult{T1, T2, T3, T4}
    S1::T1
    S1_Conf_Int::T2
    S2::T3
    S2_Conf_Int::T4
    ST::T1
    ST_Conf_Int::T2
end

function fuse_designs(A, B; second_order=false)
    d = size(A,1)
    Aᵦ = [copy(A) for i in 1:d]
    for i in 1:d
        Aᵦ[i][i,:] = @view(B[i,:])
    end
    if second_order
        Bₐ = [copy(B) for i in 1:d]
        for i in 1:d
            Bₐ[i][i,:] = @view(A[i,:])
        end
        return hcat(A,B,reduce(hcat,Aᵦ),reduce(hcat,Bₐ))
    end
    hcat(A,B,reduce(hcat,Aᵦ))
end

function gsa(f, method::Sobol, A::AbstractMatrix{TA}, B::AbstractMatrix;
             batch=false, Ei_estimator = :Jansen1999, distributed::Val{SHARED_ARRAY} = Val(false), kwargs...) where {TA, SHARED_ARRAY}
    d, n = size(A)
    nboot = method.nboot # load to help alias analysis
    n = n ÷ nboot
    multioutput = false
    Anb = Vector{Matrix{TA}}(undef, nboot)
    for i ∈ 1:nboot
        Anb[i] = A[:,n*(i-1)+1:n*(i)]
    end
    Bnb = Vector{Matrix{TA}}(undef, nboot)
    for i ∈ 1:nboot
        Bnb[i] = B[:,n*(i-1)+1:n*(i)]
    end
    _all_points = mapreduce((args...) -> fuse_designs(args...;second_order = 2 in method.order), hcat, Anb, Bnb)
    if SHARED_ARRAY && isbits(TA)
        all_points = SharedMatrix{TA}(size(_all_points))
        all_points .= _all_points
    else
        all_points = _all_points
    end

    if batch
        all_y = f(all_points)
        multioutput = all_y isa AbstractMatrix
        y_size = nothing
        gsa_sobol_all_y_analysis(method, all_y, d, n, Ei_estimator, y_size, Val(multioutput))
    else
        _y = [f(all_points[:, i]) for i in 1:size(all_points, 2)]
        multioutput = !(eltype(_y) <: Number)
        if eltype(_y) <: RecursiveArrayTools.AbstractVectorOfArray
            y_size = size(_y[1])
            _y = vec.(_y)
        else
            y_size = nothing
        end
        if multioutput
            gsa_sobol_all_y_analysis(method, reduce(hcat, _y), d, n, Ei_estimator, y_size, Val(true))
        else
            gsa_sobol_all_y_analysis(method, _y, d, n, Ei_estimator, y_size, Val(false))
        end
    end
end
function gsa_sobol_all_y_analysis(method, all_y::AbstractArray{T}, d, n, Ei_estimator, y_size, ::Val{multioutput}) where {T, multioutput}
    nboot = method.nboot
    Eys = multioutput ? Matrix{T}[] : T[]
    Varys = multioutput ? Matrix{T}[] : T[]
    Vᵢs = multioutput ? Matrix{T}[] : Vector{T}[]
    Vᵢⱼs = multioutput ?  Array{T, 3}[] : Matrix{T}[]
    Eᵢs = multioutput ? Matrix{T}[] : Vector{T}[]
    step = 2 in method.order ? 2*d+2 : d+2
    if !multioutput
        for i in 1:step:step*nboot
            push!(Eys,mean(all_y[(i-1)*n+1:(i+1)*n]))
            push!(Varys,var(all_y[(i-1)*n+1:(i+1)*n]))

            fA = all_y[(i-1)*n+1:i*n]
            fB = all_y[(i*n+1):(i+1)*n]
            fAⁱ= [all_y[(j*n+1):((j+1)*n)] for j in i+1:(i+d)]
            if 2 in method.order
                fBⁱ= [all_y[(j*n+1):((j+1)*n)] for j in i+d+1:(i+2*d)]
            end

            push!(Vᵢs, [sum(fB.*(fAⁱ[k].-fA)) for k in 1:d]./n)
            if 2 in method.order
                M = zeros(T, d, d)
                for k in 1:d
                    for j in k+1:d
                        M[k,j] = sum((fBⁱ[k] .* fAⁱ[j]) .- (fA .* fB))/n
                    end
                end
                push!(Vᵢⱼs, M)
            end
            if Ei_estimator === :Homma1996
                push!(Eᵢs,[Varys[i] .- sum(fA .* fAⁱ[k])./(n) + Eys[i].^2 for k in 1:d])
            elseif Ei_estimator === :Sobol2007
                push!(Eᵢs,[sum(abs2,fA-fAⁱ[k]) for k in 1:d]./(2n))
            elseif Ei_estimator === :Jansen1999
                push!(Eᵢs,[sum(fA.*(fA.-fAⁱ[k])) for k in 1:d]./(n))
            end
        end
    else
        for i in 1:step:step*nboot
            push!(Eys,mean(all_y[:, (i-1)*n+1:(i+1)*n],dims=2))
            push!(Varys,var(all_y[:, (i-1)*n+1:(i+1)*n],dims=2))

            fA = all_y[:, (i-1)*n+1:i*n]
            fB = all_y[:, (i*n+1):(i+1)*n]
            fAⁱ= [all_y[:, (j*n+1):((j+1)*n)] for j in i+1:(i+d)]
            if 2 in method.order
                fBⁱ= [all_y[:, (j*n+1):((j+1)*n)] for j in i+d+1:(i+2*d)]
            end

            push!(Vᵢs,reduce(hcat, [sum(fB.*(fAⁱ[k].-fA), dims=2) for k in 1:d]./n))

            if 2 in method.order
                M = zeros(T, d, d, length(Eys[1]))
                for k in 1:d
                    for j in k+1:d
                        Vₖⱼ = sum((fBⁱ[k] .* fAⁱ[j]) .- (fA .* fB), dims = 2)/n
                        for l in 1:length(Eys[1])
                            M[k,j,l] = Vₖⱼ[l]
                        end
                    end
                end
                push!(Vᵢⱼs,M)
            end
            if Ei_estimator === :Homma1996
                push!(Eᵢs,reduce(hcat, [Varys[i] .- sum(fA .* fAⁱ[k], dims=2)./(n) + Eys[i].^2 for k in 1:d]))
            elseif Ei_estimator === :Sobol2007
                push!(Eᵢs,reduce(hcat, [sum(abs2,fA-fAⁱ[k],dims=2) for k in 1:d]./(2n)))
            elseif Ei_estimator === :Jansen1999
                push!(Eᵢs,reduce(hcat, [sum(fA.*(fA.-fAⁱ[k]), dims=2) for k in 1:d]./(n)))
            end
        end
    end
    if 2 in method.order
        Sᵢⱼs = similar(Vᵢⱼs)
        for i ∈ 1:nboot
            if !multioutput
                M = zeros(T, d, d)
                for k in 1:d
                    for j in k+1:d
                        M[k,j] =Vᵢs[i][k] + Vᵢs[i][j]
                    end
                end
                Sᵢⱼs[i] = (Vᵢⱼs[i] - M) ./ Varys[i]
            else
                M = zeros(T, d, d, length(Eys[1]))
                for l in 1:length(Eys[1])
                    for k in 1:d
                        for j in k+1:d
                            M[k,j,l] = Vᵢs[i][l,k] + Vᵢs[i][l, j]
                        end
                    end
                end
                Sᵢⱼs[i] = cat([(Vᵢⱼs[i][:,:,l] - M[:,:,l]) ./ Varys[i][l] for l in 1:length(Eys[1])]...; dims = 3)
            end
        end
    end

    Sᵢs = [Vᵢs[i] ./Varys[i] for i in 1:nboot]
    Tᵢs = [Eᵢs[i] ./Varys[i] for i in 1:nboot]
    if nboot > 1
        size_ = size(Sᵢs[1])
        S1 = [[Sᵢ[i] for Sᵢ in Sᵢs] for i in 1:length(Sᵢs[1])]
        ST = [[Tᵢ[i] for Tᵢ in Tᵢs] for i in 1:length(Tᵢs[1])]

        function calc_ci(x,mean=nothing)
            alpha = (1 - method.conf_level)
            std(x,mean=mean)/sqrt(length(x))
        end
        S1_CI = map(calc_ci,S1)
        ST_CI = map(calc_ci,ST)

        if 2 in method.order
            size__= size(Sᵢⱼs[1])
            S2_CI = Array{T}(undef, size__)
            Sᵢⱼ = Array{T}(undef, size__)
            b = getindex.(Sᵢⱼs, 1)
            Sᵢⱼ[1] = b̄ = mean(b)
            S2_CI[1] = calc_ci(b, b̄)
            for i ∈ 2:length(Sᵢⱼs[1])
                b .= getindex.(Sᵢⱼs, i)
                Sᵢⱼ[i] = b̄ = mean(b)
                S2_CI[i] = calc_ci(b, b̄)
            end
        end
        Sᵢ = reshape(mean.(S1),size_...)
        Tᵢ = reshape(mean.(ST),size_...)
    else
        Sᵢ = Sᵢs[1]
        Tᵢ = Tᵢs[1]
        if 2 in method.order
            Sᵢⱼ = Sᵢⱼs[1]
        end
    end
    if isnothing(y_size)
        _Sᵢ = Sᵢ
        _Tᵢ = Tᵢ
    else
        f_shape = let y_size = y_size
            x -> [reshape(x[:,i],y_size) for i in 1:size(x,2)]
        end
        _Sᵢ = f_shape(Sᵢ)
        _Tᵢ = f_shape(Tᵢ)
    end
    return SobolResult(_Sᵢ,
                     nboot > 1 ? reshape(S1_CI,size_...) : nothing,
                     2 in method.order ? Sᵢⱼ : nothing,
                     nboot > 1 && 2 in method.order ? S2_CI : nothing,
                     _Tᵢ ,
                     nboot > 1 ? reshape(ST_CI,size_...) : nothing)
end

function gsa(f,method::Sobol,p_range::AbstractVector; N, kwargs...)
    AB = QuasiMonteCarlo.generate_design_matrices(N, [i[1] for i in p_range], [i[2] for i in p_range], QuasiMonteCarlo.SobolSample(),2*method.nboot)
    A = reduce(hcat, @view(AB[1:method.nboot]))
    B = reduce(hcat, @view(AB[method.nboot+1:end]))
    gsa(f, method, A, B; kwargs...)
end


text:
export platform_key, platform_key_abi, platform_dlext, valid_dl_path, arch, libc, call_abi, wordsize, triplet, choose_download,
       CompilerABI, Platform, UnknownPlatform, Linux, MacOS, Windows, FreeBSD
import Base: show

abstract type Platform end

struct UnknownPlatform <: Platform
end

# We need to track our compiler ABI compatibility.
struct CompilerABI
    # Major GCC version that we're locked into.
    # Can be [:gcc4, :gcc5, :gcc6, :gcc7, :gcc8, :gcc_any]
    gcc_version::Symbol

    # Whether we're using cxx11abi strings
    # Can be [:cxx03, :cxx11, :cxx_any]
    cxx_abi::Symbol

    function CompilerABI(gcc_version::Symbol = :gcc_any, cxx_abi::Symbol = :cxx_any)
        if !in(gcc_version, [:gcc4, :gcc5, :gcc6, :gcc7, :gcc8, :gcc_any])
            throw(ArgumentError("Unsupported GCC major version '$gcc_version'"))
        end

        if !in(cxx_abi, [:cxx03, :cxx11, :cxx_any])
            throw(ArgumentError("Unsupported string ABI '$cxx_abi'"))
        end

        return new(gcc_version, cxx_abi)
    end
end

function show(io::IO, cabi::CompilerABI)
    write(io, "CompilerABI(")
    if cabi.gcc_version != :gcc_any || cabi.cxx_abi != :cxx_any
        write(io, "$(repr(cabi.gcc_version))")
    end
    if cabi.cxx_abi != :cxx_any
        write(io, ", $(repr(cabi.cxx_abi))")
    end
    write(io, ")")
end

struct Linux <: Platform
    arch::Symbol
    libc::Symbol
    call_abi::Symbol
    compiler_abi::CompilerABI

    function Linux(arch::Symbol;
                   libc::Symbol=:glibc,
                   call_abi::Symbol=:default_abi,
                   compiler_abi::CompilerABI=CompilerABI())
        if !in(arch, [:i686, :x86_64, :aarch64, :powerpc64le, :armv7l])
            throw(ArgumentError("Unsupported architecture '$arch' for Linux"))
        end

        # The default libc on Linux is glibc
        if libc === :blank_libc
            libc = :glibc
        end

        if !in(libc, [:glibc, :musl])
            throw(ArgumentError("Unsupported libc '$libc' for Linux"))
        end

        # The default calling abi on Linux is blank (e.g. not "eabi" or "eabihf"
        # or anything like that, just "blank"), so map over to that by default
        # except on armv7l, where we map it over to :eabihf
        if call_abi === :default_abi
            if arch != :armv7l
                call_abi = :blank_abi
            else
                call_abi = :eabihf
            end
        end

        if !in(call_abi, [:eabihf, :blank_abi])
            throw(ArgumentError("Unsupported calling abi '$call_abi' for Linux"))
        end

        # If we're constructing for armv7l, we MUST have the eabihf abi
        if arch == :armv7l && call_abi != :eabihf
            throw(ArgumentError("armv7l Linux must use eabihf, not '$call_abi'"))
        end
        # ...and vice-versa
        if arch != :armv7l && call_abi == :eabihf
            throw(ArgumentError("eabihf Linux is only on armv7l, not '$arch'!"))
        end

        return new(arch, libc, call_abi, compiler_abi)
    end
end

struct MacOS <: Platform
    arch::Symbol
    libc::Symbol
    call_abi::Symbol
    compiler_abi::CompilerABI

    # Provide defaults for everything because there's really only one MacOS
    # target right now.  Maybe someday iOS.  :fingers_crossed:
    function MacOS(arch::Symbol=:x86_64;
                   libc::Symbol=:blank_libc,
                   call_abi::Symbol=:blank_abi,
                   compiler_abi::CompilerABI=CompilerABI())
        if arch !== :x86_64
            throw(ArgumentError("Unsupported architecture '$arch' for macOS"))
        end
        if libc !== :blank_libc
            throw(ArgumentError("Unsupported libc '$libc' for macOS"))
        end
        if call_abi !== :blank_abi
            throw(ArgumentError("Unsupported abi '$call_abi' for macOS"))
        end

        return new(arch, libc, call_abi, compiler_abi)
    end
end

struct Windows <: Platform
    arch::Symbol
    libc::Symbol
    call_abi::Symbol
    compiler_abi::CompilerABI

    function Windows(arch::Symbol;
                     libc::Symbol=:blank_libc,
                     call_abi::Symbol=:blank_abi,
                     compiler_abi::CompilerABI=CompilerABI())
        if !in(arch, [:i686, :x86_64])
            throw(ArgumentError("Unsupported architecture '$arch' for Windows"))
        end
        # We only support the one libc/abi on Windows, so no need to play
        # around with "default" values.
        if libc !== :blank_libc
            throw(ArgumentError("Unsupported libc '$libc' for Windows"))
        end
        if call_abi !== :blank_abi
            throw(ArgumentError("Unsupported abi '$call_abi' for Windows"))
        end

        return new(arch, libc, call_abi, compiler_abi)
    end
end

struct FreeBSD <: Platform
    arch::Symbol
    libc::Symbol
    call_abi::Symbol
    compiler_abi::CompilerABI

    function FreeBSD(arch::Symbol;
                     libc::Symbol=:blank_libc,
                     call_abi::Symbol=:default_abi,
                     compiler_abi::CompilerABI=CompilerABI())
        # `uname` on FreeBSD reports its architecture as amd64 and i386 instead of x86_64
        # and i686, respectively. In the off chance that Julia hasn't done the mapping for
        # us, we'll do it here just in case.
        if arch === :amd64
            arch = :x86_64
        elseif arch === :i386
            arch = :i686
        elseif !in(arch, [:i686, :x86_64, :aarch64, :powerpc64le, :armv7l])
            throw(ArgumentError("Unsupported architecture '$arch' for FreeBSD"))
        end

        # The only libc we support on FreeBSD is the blank libc, which corresponds to
        # FreeBSD's default libc
        if libc !== :blank_libc
            throw(ArgumentError("Unsupported libc '$libc' for FreeBSD"))
        end

        # The default abi on FreeBSD is blank, execpt on armv7l
        if call_abi === :default_abi
            if arch != :armv7l
                call_abi = :blank_abi
            else
                call_abi = :eabihf
            end
        end

        if !in(call_abi, [:eabihf, :blank_abi])
            throw(ArgumentError("Unsupported calling abi '$call_abi' for FreeBSD"))
        end

        # If we're constructing for armv7l, we MUST have the eabihf abi
        if arch == :armv7l && call_abi != :eabihf
            throw(ArgumentError("armv7l FreeBSD must use eabihf, no '$call_abi'"))
        end
        # ...and vice-versa
        if arch != :armv7l && call_abi == :eabihf
            throw(ArgumentError("eabihf FreeBSD is only on armv7l, not '$arch'!"))
        end

        return new(arch, libc, call_abi, compiler_abi)
    end
end

"""
    platform_name(p::Platform)

Get the "platform name" of the given platform.  E.g. returns "Linux" for a
`Linux` object, or "Windows" for a `Windows` object.
"""
platform_name(p::Linux) = "Linux"
platform_name(p::MacOS) = "MacOS"
platform_name(p::Windows) = "Windows"
platform_name(p::FreeBSD) = "FreeBSD"
platform_name(p::UnknownPlatform) = "UnknownPlatform"

"""
    arch(p::Platform)

Get the architecture for the given `Platform` object as a `Symbol`.

# Examples
```jldoctest
julia> arch(Linux(:aarch64))
:aarch64

julia> arch(MacOS())
:x86_64
```
"""
arch(p::Platform) = p.arch
arch(u::UnknownPlatform) = :unknown

"""
    libc(p::Platform)

Get the libc for the given `Platform` object as a `Symbol`.

# Examples
```jldoctest
julia> libc(Linux(:aarch64))
:glibc

julia> libc(FreeBSD(:x86_64))
:default_libc
```
"""
libc(p::Platform) = p.libc
libc(u::UnknownPlatform) = :unknown

"""
   call_abi(p::Platform)

Get the calling ABI for the given `Platform` object as a `Symbol`.

# Examples
```jldoctest
julia> call_abi(Linux(:x86_64))
:blank_abi

julia> call_abi(FreeBSD(:armv7l))
:eabihf
```
"""
call_abi(p::Platform) = p.call_abi
call_abi(u::UnknownPlatform) = :unknown

"""
    compiler_abi(p::Platform)

Get the compiler ABI object for the given `Platform`
# Examples
```jldoctest
julia> compiler_abi(Linux(:x86_64))
CompilerABI(:gcc_any, :cxx_any)

julia> compiler_abi(Linux(:x86_64; compiler_abi=CompilerABI(:gcc7)))
CompilerABI(:gcc7, :cxx_any)
```
"""
compiler_abi(p::Platform) = p.compiler_abi
compiler_abi(p::UnknownPlatform) = CompilerABI()

"""
    wordsize(platform)

Get the word size for the given `Platform` object.

# Examples
```jldoctest
julia> wordsize(Linux(:arm7vl))
32

julia> wordsize(MacOS())
64
```
"""
wordsize(p::Platform) = (arch(p) === :i686 || arch(p) === :armv7l) ? 32 : 64
wordsize(u::UnknownPlatform) = 0

"""
    triplet(platform)

Get the target triplet for the given `Platform` object as a `String`.

# Examples
```jldoctest
julia> triplet(MacOS())
"x86_64-apple-darwin14"

julia> triplet(Windows(:i686))
"i686-w64-mingw32"

julia> triplet(Linux(:armv7l, :default_libc, :default_abi, CompilerABI(:gcc4))
"arm-linux-gnueabihf-gcc4"
```
"""
triplet(p::Platform) = string(
    arch_str(p),
    vendor_str(p),
    libc_str(p),
    call_abi_str(p),
    compiler_abi_str(p),
)
vendor_str(p::Windows) = "-w64-mingw32"
vendor_str(p::MacOS) = "-apple-darwin14"
vendor_str(p::Linux) = "-linux"
vendor_str(p::FreeBSD) = "-unknown-freebsd11.1"

# Special-case UnknownPlatform
triplet(p::UnknownPlatform) = "unknown-unknown-unknown"

# Helper functions for Linux and FreeBSD libc/abi mishmashes
arch_str(p::Platform) = (arch(p) == :armv7l) ? "arm" : string(arch(p))
function libc_str(p::Platform)
    if libc(p) == :blank_libc
        return ""
    elseif libc(p) == :glibc
        return "-gnu"
    else
        return "-$(libc(p))"
    end
end
call_abi_str(p::Platform) = (call_abi(p) == :blank_abi) ? "" : string(call_abi(p))
function compiler_abi_str(cabi::CompilerABI)
    str = ""
    if cabi.gcc_version != :gcc_any
        str *= "-$(cabi.gcc_version)"
    end
    if cabi.cxx_abi != :cxx_any
        str *= "-$(cabi.cxx_abi)"
    end
    return str
end
compiler_abi_str(p::Platform) = compiler_abi_str(compiler_abi(p))

Sys.isapple(p::Platform) = p isa MacOS
Sys.islinux(p::Platform) = p isa Linux
Sys.iswindows(p::Platform) = p isa Windows
Sys.isbsd(p::Platform) = (p isa FreeBSD) || (p isa MacOS)

"""
    platform_key_abi(machine::AbstractString)

Returns the platform key for the current platform, or any other though the
the use of the `machine` parameter.
"""
function platform_key_abi(machine::AbstractString)
    # We're going to build a mondo regex here to parse everything:
    arch_mapping = Dict(
        :x86_64 => "(x86_|amd)64",
        :i686 => "i\\d86",
        :aarch64 => "aarch64",
        :armv7l => "arm(v7l)?",
        :powerpc64le => "p(ower)?pc64le",
    )
    platform_mapping = Dict(
        :darwin => "-apple-darwin[\\d\\.]*",
        :freebsd => "-(.*-)?freebsd[\\d\\.]*",
        :mingw32 => "-w64-mingw32",
        :linux => "-(.*-)?linux",
    )
    libc_mapping = Dict(
        :blank_libc => "",
        :glibc => "-gnu",
        :musl => "-musl",
    )
    call_abi_mapping = Dict(
        :blank_abi => "",
        :eabihf => "eabihf",
    )
    gcc_version_mapping = Dict(
        :gcc_any => "",
        :gcc4 => "-gcc4",
        :gcc5 => "-gcc5",
        :gcc6 => "-gcc6",
        :gcc7 => "-gcc7",
        :gcc8 => "-gcc8",
    )
    cxx_abi_mapping = Dict(
        :cxx_any => "",
        :cxx03 => "-cxx03",
        :cxx11 => "-cxx11",
    )

    # Helper function to collapse dictionary of mappings down into a regex of
    # named capture groups joined by "|" operators
    c(mapping) = string("(",join(["(?<$k>$v)" for (k, v) in mapping], "|"), ")")

    triplet_regex = Regex(string(
        "^",
        c(arch_mapping),
        c(platform_mapping),
        c(libc_mapping),
        c(call_abi_mapping),
        c(gcc_version_mapping),
        c(cxx_abi_mapping),
        "\$",
    ))

    m = match(triplet_regex, machine)
    if m != nothing
        # Helper function to find the single named field within the giant regex
        # that is not `nothing` for each mapping we give it.
        get_field(m, mapping) = begin
            for k in keys(mapping)
                if m[k] != nothing
                   return k
                end
            end
        end

        # Extract the information we're interested in:
        arch = get_field(m, arch_mapping)
        platform = get_field(m, platform_mapping)
        libc = get_field(m, libc_mapping)
        call_abi = get_field(m, call_abi_mapping)
        gcc_version = get_field(m, gcc_version_mapping)
        cxx_abi = get_field(m, cxx_abi_mapping)

        # First, figure out what platform we're dealing with, then sub that off
        # to the appropriate constructor.  If a constructor runs into trouble,
        # catch the error and return `UnknownPlatform()` here to be nicer to client code.
        ctors = Dict(:darwin => MacOS, :mingw32 => Windows, :freebsd => FreeBSD, :linux => Linux)
        try
            T = ctors[platform]
            compiler_abi = CompilerABI(gcc_version, cxx_abi)
            return T(arch, libc=libc, call_abi=call_abi, compiler_abi=compiler_abi)
        catch
        end
    end

    @warn("Platform `$(machine)` is not an officially supported platform")
    return UnknownPlatform()
end

# Define show() for these Platform objects for two reasons:
#  - I don't like the `BinaryProvider.` at the beginning of the types;
#    it's unnecessary as these are exported
#  - I don't like the :blank_*/:any arguments, they're unnecessary
function show(io::IO, p::Platform)
    write(io, "$(platform_name(p))($(repr(arch(p)))")

    if libc(p) != :blank_libc
        write(io, ", libc=$(repr(libc(p)))")
    end
    if call_abi(p) != :blank_abi
        write(io, ", call_abi=$(repr(call_abi(p)))")
    end
    cabi = compiler_abi(p)
    if cabi.gcc_version != :gcc_any || cabi.cxx_abi != :cxx_any
        write(io, ", compiler_abi=$(repr(cabi))")
    end
    write(io, ")")
end

"""
    platform_dlext(platform::Platform = platform_key_abi())

Return the dynamic library extension for the given platform, defaulting to the
currently running platform.  E.g. returns "so" for a Linux-based platform,
"dll" for a Windows-based platform, etc...
"""
platform_dlext(::Linux) = "so"
platform_dlext(::FreeBSD) = "so"
platform_dlext(::MacOS) = "dylib"
platform_dlext(::Windows) = "dll"
platform_dlext(::UnknownPlatform) = "unknown"
platform_dlext() = platform_dlext(platform_key_abi())

"""
    parse_dl_name_version(path::AbstractString, platform::Platform)

Given a path to a dynamic library, parse out what information we can
from the filename.  E.g. given something like "lib/libfoo.so.3.2",
this function returns `"libfoo", v"3.2"`.  If the path name is not a
valid dynamic library, this method throws an error.  If no soversion
can be extracted from the filename, as in "libbar.so" this method
returns `"libbar", nothing`.
"""
function parse_dl_name_version(path::AbstractString, platform::Platform)
    dlext_regexes = Dict(
        # On Linux, libraries look like `libnettle.so.6.3.0`
        "so" => r"^(.*?).so((?:\.[\d]+)*)$",
        # On OSX, libraries look like `libnettle.6.3.dylib`
        "dylib" => r"^(.*?)((?:\.[\d]+)*).dylib$",
        # On Windows, libraries look like `libnettle-6.dylib`
        "dll" => r"^(.*?)(?:-((?:[\.\d]+)*))?.dll$"
    )

    # Use the regex that matches this platform
    dlregex = dlext_regexes[platform_dlext(platform)]
    m = match(dlregex, basename(path))
    if m === nothing
        throw(ArgumentError("Invalid dynamic library path '$path'"))
    end

    # Extract name and version
    name = m.captures[1]
    version = m.captures[2]
    if version === nothing || isempty(version)
        version = nothing
    else
        version = VersionNumber(strip(version, '.'))
    end
    return name, version
end

"""
    valid_dl_path(path::AbstractString, platform::Platform)

Return `true` if the given `path` ends in a valid dynamic library filename.
E.g. returns `true` for a path like `"usr/lib/libfoo.so.3.5"`, but returns
`false` for a path like `"libbar.so.f.a"`.
"""
function valid_dl_path(path::AbstractString, platform::Platform)
    try
        parse_dl_name_version(path, platform)
        return true
    catch
        return false
    end
end

"""
    detect_libgfortran_abi(libgfortran_name::AbstractString)

Examines the given libgfortran SONAME to see what version of GCC corresponds
to the given libgfortran version.
"""
function detect_libgfortran_abi(libgfortran_name::AbstractString, platform::Platform = platform_key_abi(Sys.MACHINE))
    # Extract the version number from this libgfortran.  Ironically, parse_dl_name_version()
    # wants a Platform, but we may not have initialized the default platform key yet when we
    # run this method for the first time (since we need the output of this function to set
    # that default platform) so we manually pass in Sys.MACHINE.  :P
    name, version = parse_dl_name_version(libgfortran_name, platform)
    if version === nothing
        @warn("Unable to determine libgfortran version from '$(libgfortran_name)'; returning :gcc_any")
        return :gcc_any
    end
    libgfortran_to_gcc = Dict(
        3 => :gcc4,
        4 => :gcc7,
        5 => :gcc8,
    )
    if !in(version.major, keys(libgfortran_to_gcc))
        @warn("Unsupported libgfortran version '$version'; returning :gcc_any")
        return :gcc_any
    end
    return libgfortran_to_gcc[version.major]
end

"""
    detect_libgfortran_abi()

If no parameter is given, introspects the current Julia process to determine
the version of GCC this Julia was built with.
"""
function detect_libgfortran_abi()
    libgfortran_paths = filter(x -> occursin("libgfortran", x), dllist())
    if isempty(libgfortran_paths)
         # One day, I hope to not be linking against libgfortran in base Julia
        return :gcc_any
    end
    return detect_libgfortran_abi(first(libgfortran_paths))
end

"""
    detect_libstdcxx_abi()

Introspects the currently running Julia process to find out what version of libstdc++
it is linked to (if any), as a proxy for GCC version compatibility.  E.g. if we are
linked against libstdc++.so.19, binary dependencies built by GCC 8.1.0 will have linker
errors.  This method returns the maximum GCC abi that we can support.
"""
function detect_libstdcxx_abi()
    libstdcxx_paths = filter(x -> occursin("libstdc++", x), dllist())
    if isempty(libstdcxx_paths)
        # This can happen if we were built by clang, so we don't link against
        # libstdc++ at all.
        return :gcc_any
    end

    # Extract all pieces of `.gnu.version_d` from libstdc++.so, find the `GLIBCXX_*`
    # symbols, and use the maximum version of that to find the GLIBCXX ABI version number
    #version_symbols = readmeta(first(libstdcxx_paths)) do oh
    #    unique(vcat((x -> x.names).(ELFVersionData(oh))...))
    #end
    #version_symbols = filter(x -> startswith(x, "GLIBCXX_"), version_symbols)
    #max_version = maximum([VersionNumber(split(v, "_")[2]) for v in version_symbols])

    # ^^ Okay, that's really cool, but unfortunately it introduces a dependency on
    # ObjectFile which is unacceptable for us.  So instead we just brute-force it.
    max_version = v"3.4.0"
    hdl = dlopen(first(libstdcxx_paths))
    for minor_version in 1:26
        if dlsym_e(hdl, "GLIBCXX_3.4.$(minor_version)") != C_NULL
            max_version = VersionNumber("3.4.$(minor_version)")
        end
    end
    dlclose(hdl)

    # Full list available here: https://gcc.gnu.org/onlinedocs/libstdc++/manual/abi.html
    if max_version < v"3.4.18"
        @warn "Cannot make sense of autodetected libstdc++ ABI version ('$max_version')"
        return :gcc_any
    elseif max_version < v"3.4.23"
        # If we aren't up to 7.1.0, then we fall all the way back to 4.8.5
        return :gcc4
    elseif max_version < v"3.4.25"
        return :gcc7
    else
        return :gcc8
    end
end

"""
    detect_cxx11_string_abi()

Introspects the currently running Julia process to see what version of the C++11 string
ABI it was compiled with.  (In reality, it checks for symbols within LLVM, but that is
close enough for our purposes, as you can't mix linkages between Julia and LLVM if they
are not compiled in the same way).
"""
function detect_cxx11_string_abi()
    function open_libllvm()
        for lib_name in ("libLLVM", "LLVM", "libLLVMSupport")
            hdl = dlopen_e(lib_name)
            if hdl != C_NULL
                return hdl
            end
        end
        error("Unable to open libLLVM!")
    end

    hdl = open_libllvm()
    # Check for llvm::sys::getProcessTriple(), first without cxx11 tag:
    if dlsym_e(hdl, "_ZN4llvm3sys16getProcessTripleEv") != C_NULL
        return :cxx03
    elseif dlsym_e(hdl, "_ZN4llvm3sys16getProcessTripleB5cxx11Ev") != C_NULL
        return :cxx11
    else
        error("Unable to find llvm::sys::getProcessTriple() in libLLVM!")
    end
end

function detect_compiler_abi()
    gcc_version = detect_libgfortran_abi()
    cxx11_string_abi = detect_cxx11_string_abi()

    # If we have no constraint from libgfortran linkage (impossible within current
    # Julia, but let's be planning for the future here) then inspect libstdc++.
    if gcc_version == :gcc_any
        gcc_version = detect_libstdcxx_abi()
    end

    return CompilerABI(gcc_version, cxx11_string_abi)
end

# Cache the default platform_key_abi() since that's by far the most common way
# we call platform_key_abi(), and we don't want to parse the same thing over
# and over and over again.  Note that we manually slap on a compiler abi
# string onto the end, so as to encode that in Sys.MACHINE like we expect our
# triplets to be encoded.
default_platkey = platform_key_abi(string(
    Sys.MACHINE,
    compiler_abi_str(detect_compiler_abi()),
))
function platform_key_abi()
    global default_platkey
    return default_platkey
end

function platforms_match(a::Platform, b::Platform)
    # Check to see if a and b  satisfy the rigid constraints first, these are
    # things that are simple equality checks:
    function rigid_constraints(a, b)
        return (typeof(a) <: typeof(b) || typeof(b) <: typeof(a)) &&
               (arch(a) == arch(b)) && (libc(a) == libc(b)) &&
               (call_abi(a) == call_abi(b))
    end

    # The flexible constraints are ones that can do equals, but also have things
    # like "any" values, etc....
    function flexible_constraints(a, b)
        ac = compiler_abi(a)
        bc = compiler_abi(b)

        # Map from GCC version to libgfortran SO version
        gfmap = Dict(
            :gcc4 => 3,
            :gcc5 => 3,
            :gcc6 => 3,
            :gcc7 => 4,
            :gcc8 => 5,
        )

        # We consider two GCC versions to match if their libgfortran
        # versions match; e.g. :gcc4 and :gcc5 match.
        gcc_match = (ac.gcc_version == :gcc_any
                  || bc.gcc_version == :gcc_any
                  || gfmap[ac.gcc_version] == gfmap[bc.gcc_version])
        cxx_match = (ac.cxx_abi == :cxx_any
                  || bc.cxx_abi == :cxx_any
                  || ac.cxx_abi == bc.cxx_abi)
        return gcc_match && cxx_match
    end

    return rigid_constraints(a, b) && flexible_constraints(a, b)
end

"""
    choose_download(download_info::Dict, platform::Platform = platform_key_abi())

Given a `download_info` dictionary mapping platforms to some value, choose
the value whose key best matches `platform`, returning `nothing` if no matches
can be found.

Platform attributes such as architecture, libc, calling ABI, etc... must all
match exactly, however attributes such as compiler ABI can have wildcards
within them such as `:gcc_any` which matches any version of GCC.
"""
function choose_download(download_info::Dict, platform::Platform = platform_key_abi())
    ps = collect(filter(p -> platforms_match(p, platform), keys(download_info)))

    if isempty(ps)
        return nothing
    end

    # At this point, we may have multiple possibilities.  E.g. if, in the future,
    # Julia can be built without a direct dependency on libgfortran, we may match
    # multiple tarballs that vary only within their libgfortran ABI.  To narrow it
    # down, we just sort by triplet, then pick the last one.  This has the effect
    # of generally choosing the latest release (e.g. a `libgfortran5` tarball
    # rather than a `libgfortran3` tarball)
    p = last(sort(ps, by = p -> triplet(p)))
    return download_info[p]
end


text:
using ControlSystems, Plots
gr()

sys = feedback(1.0, ss(-1.0, 2, 1, 0) * (delay(2.0) + delay(3.0) + delay(2.5)))
sys = feedback(ss(-1.0, 1, 1, 0), delay(1.0))

t = 0:0.02:8
@time y, t, x = lsim(sys, t-> [t>=0 ? 1.0 : 0.0], t)
@time y, t, x = lsim(sys, [1.0], t)
@time y, t, x = lsim(sys, (out, t) -> (out .= (t>=0 ? 1.0 : 0.0)), t)
@time y, t, x = lsim(sys, (out, t) -> (out[1] = (t>=0 ? 1.0 : 0.0)), t)

function u0(out,t)
    if t > 0
        out[1] = 1
    else
        out[1] = 0
    end
    return
end

@time y, t, x = lsim(sys, u0, t)

plot(t, y')

s = tf("s")
P = delay(2.6)*ss((s+3.0)/(s^2+0.3*s+1))
C = 0.06 * ss(1.0 + 1/s);
P*C
T = feedback(P*C,1.0)

t = 0:0.1:70
y, t, x = lsim(T, t -> (t<0 ? 0 : 1 ), t)
plot(t, y, c = :blue)

w = 10 .^ (-2:0.01:2)
marginplot(P*C, w)
marginplot(P*C)

notch = ss(tf([1, 0.2, 1],[1, .8, 1]));
C = ss(0.05 * (1 + 1/s));
Tnotch = feedback(P*C*notch, 1.0)

stepplot(Tnotch)

y, t, x = step(C, method=:zoh)

y2, t2, x2 = step(Tnotch)
stepplot(Tnotch)

stepplot(Tnotch, 40, 0.1)

stepplot(T, 100)

G = delay(5)/(s+1)
T = feedback(G, 0.5)
w = 10 .^ (-2:0.01:3)
bodeplot(T, w, plotphase=false)

# Test conversion, promotion
delay(1,Int64) + 3.5

G = 1 + 0.5 * delay(3)
w = 10 .^(-2:0.001:2)
bodeplot(G, w, plotphase=false)

G = delay(1) * ((0.8*s^2+s+2)/(s^2+s))
T = feedback(G,1)
# Not possible with direct term
stepplot(T)

bodeplot(T)

G = 1/(s+1) + delay(4)
T = feedback(1,G)
# Not possible to lsim with direct term
stepplot(T)
bodeplot(T)

s = tf("s")


text:
using BinaryBuilder

# Collection of sources required to build Pixman
name = "Pixman"
version = v"0.34.0"
sources = [
    "https://www.cairographics.org/releases/pixman-0.34.0.tar.gz" =>
    "21b6b249b51c6800dc9553b65106e1e37d0e25df942c90531d4c3997aa20a88e",
    "patches",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd pixman-0.34.0/

# Apply patch for compilation with clang
patch < $WORKSPACE/srcdir/patches/clang.patch

# Apply patch for arm on musl
patch -p1 < $WORKSPACE/srcdir/patches/arm_musl.patch

aclocal
automake

./configure --prefix=$prefix --host=$target
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpixman", :libpixman)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using ConstructiveGeometry
using GLMakie
# using CairoMakie

hexagon = polygon([[cospi(t/3),sinpi(t/3)] for t in 0:5])
c1, c2, c3 = parse.(ConstructiveGeometry.Colorant, ( "#cb3c33", "#9558b2", "#389826"))

bolt = union(linear_extrude(5)*(8*hexagon),
	cylinder(15,4), rotate_extrude(7*360, slide=14)*translate([1,0])*square(4,1))

m = union(c1*bolt, [20,0,0]+c2*bolt, [10,17,0]+c3*bolt)

s = plot(m)
zoom!(s, Vec3f0(0,0,0), 5., false)
s.center = false
save("logo.png", s)


text:
@testset "Table Joins for Several Pairs" begin
    import TexTables: IndexedTable, TableCol, join_table
    # Tests for issue in PR #28 (makes sure that the recursive table join
    # implementation on pairs works properly)
    n = 5
    keys = ["x$i" for i in 1:5]
    cols = Vector{IndexedTable{1,1}}(undef, 4)
    for j in 1:4
        cols[j] = hcat(TableCol("mean", keys, rand(n)),
                       TableCol("std", keys, rand(n)))
    end

    # Test that the construction doesn't error out with 4 pairs
    for k in 1:4
        @test begin
            tab = join_table(collect("($j)" => cols[j] for j in 1:k)...)
            true
        end
    end
end


text:
using Test, MyPackage, Random

Random.seed!(0)

@test sprint(MyPackage.hello, "Mike") == "Hello Mike"
@test MyPackage.randstr() == "0IPrGg0J"


text:
using ExecutableSpecifications: @given, @when, @then, @expect

@given "some precondition which does nothing" begin

end

@when "we perform a no-op step" begin

end

@then "the scenario as a whole succeeds" begin

end

@when "a step has more than one matching step definition" begin

end

@when "a step has more than one matching step definition" begin

end

@when "foo is set to value 42" begin
    context[:foo] = 42
end

@when "foo is set to value -17" begin
    context[:foo] = -17
end

@then "foo is greater than zero" begin
    @expect context[:foo] > 0
end

@then "foo is less than zero" begin
    @expect context[:foo] < 0
end


text:
function little_2var()
    v = Var(2)
    pos = compile(PlainLogicCircuit, var2lit(v))
    neg = compile(PlainLogicCircuit, -var2lit(v))
    or1 = pos | neg
    or2 = pos | neg

    v = Var(1)
    pos = compile(PlainLogicCircuit, var2lit(v))
    neg = compile(PlainLogicCircuit, -var2lit(v))

    and1 = pos & or1
    and2 = neg & or2
    and1 | and2
end

function little_3var()
    or1 = little_2var()
    v = Var(3)

    pos = compile(PlainLogicCircuit,  var2lit(v))
    neg = compile(PlainLogicCircuit, -var2lit(v))

    or2 = disjoin(children(or1))

    and1 = pos & or1
    and2 = neg & or2
    and1 | and2
end

function little_3var_constants()
    or1 = little_2var()
    v = Var(3)

    t = compile(PlainLogicCircuit, true)
    f = compile(PlainLogicCircuit, false)

    pos = compile(PlainLogicCircuit,  var2lit(v)) & t
    neg = compile(PlainLogicCircuit, -var2lit(v)) & f

    or2 = disjoin(children(or1))

    and1 = pos & or1
    and2 = neg & or2
    and1 | and2
end

function little_4var()
    ors = map(1:4) do v
        v = Var(v)
        pos = compile(PlainLogicCircuit, var2lit(v))
        neg = compile(PlainLogicCircuit, - var2lit(v))
        or = pos | neg
    end
    and1 = ors[1] & ors[2]
    and2 = ors[3] & ors[4]
    or = and1 | and2
end

function little_5var()
    c_4var = little_4var()
    v = Var(5)
    pos = compile(PlainLogicCircuit, var2lit(v))
    neg = compile(PlainLogicCircuit, - var2lit(v))
    or = pos | neg
    and = c_4var & or
    Plain⋁Node([and])
end


text:
using StiffKinetic
using Test

max_stars_repo_name: tremelow/StiffKinetic

max_issues_repo_name: tremelow/StiffKinetic

max_forks_repo_name: tremelow/StiffKinetic


text:
module WeightedGraphAlgo

export
    WeightedUndirectedGraph,
    WeightedDirectedGraph,
    MST,
    sumweights,
    LazyPrimMST,
    EagerPrimMST,
    KruskalMST,
    SP,
    DijkstraSP,
    distanceto

include("./heapqueue.jl")
using .HeapQueueModule: MinHeapQueue, heapify!, IndexedMinHQ, decreasekey!,
    contains, insert!, isempty, popfirst!

include("./unionfind.jl")
using .UnionFindAlgo: UnionFind, isconnected, union!

abstract type Edge end
abstract type WeightedEdge <: Edge end

EdgeData = Tuple{Int, Int, Float64}

# Weighted Graph Data Structure
struct WeightedUndirectedEdge <: WeightedEdge
    v1::Int
    v2::Int
    weight::Float64
    WeightedUndirectedEdge(edgedata::EdgeData) = new(edgedata...)
end

struct WeightedDirectedEdge <: WeightedEdge
    v1::Int
    v2::Int
    weight::Float64
    WeightedDirectedEdge(edgedata::EdgeData) = new(edgedata...)
end

Base.:>(e1::T, e2::T) where T<:WeightedEdge = e1.weight > e2.weight
Base.:<(e1::T, e2::T) where T<:WeightedEdge = e1.weight < e2.weight
Base.:(==)(e1::T, e2::T) where T<:WeightedEdge = e1.weight == e2.weight
Base.show(io::IO, edge::WeightedEdge) = print(io, split(string(typeof(edge)), ".")[end],
                                              "($(edge.v1), $(edge.v2), $(edge.weight))")

mutable struct Graph{T<:WeightedEdge}
    numvertices::Int
    numedges::Int
    adjlist::Vector{Vector{T}}
    function Graph{T}(numvertices::Int) where T
        new(numvertices, 0, [T[] for _ = 1:numvertices])
    end
end
Base.show(io::IO, g::Graph{T}) where T = print(io,
                    "Graph{$(split(string(T), ".")[end])}",
                    "(V=$(g.numvertices), E=$(g.numedges))")

function WeightedUndirectedGraph(graphdata::NamedTuple)
    g = Graph{WeightedUndirectedEdge}(graphdata.nv)
    addedges!(g, graphdata.edges)
    g
end

function WeightedDirectedGraph(graphdata::NamedTuple)
    g = Graph{WeightedDirectedEdge}(graphdata.nv)
    addedges!(g, graphdata.edges)
    g
  end

function addedge!(g::Graph{T}, edge::T) where T<:WeightedUndirectedEdge
    push!(g.adjlist[edge.v1], edge)
    push!(g.adjlist[edge.v2], edge)
    g.numedges += 1
end

function addedge!(g::Graph{T}, edge::T) where T<:WeightedDirectedEdge
    push!(g.adjlist[edge.v1], edge)
    g.numedges += 1
end

function addedges!(g::Graph{T}, edgesdata::Vector{EdgeData}) where {T<:WeightedEdge}
    for edgedata in edgesdata
        edge = T(edgedata)
        addedge!(g, edge)
    end
end

"""
    getedges(g)

Get edges from the adjacent edge list.
"""
function getedges(g::Graph{T}) where T<:WeightedEdge
    edgelist = T[]
    for (v, adjedges) in enumerate(g.adjlist)
        for edge in adjedges
            pushedge!(edgelist, edge, v)
        end
    end
    edgelist
end

function pushedge!(edgelist::Vector{T}, edge::T, v) where T<:WeightedUndirectedEdge
    findadjv(edge, v) < v && push!(edgelist, edge)
end
pushedge!(edgelist::Vector{T}, edge::T, v) where T<:WeightedDirectedEdge =
    push!(edgelist, edge)

struct UnfoundVertex <: Exception
    v
    edge
end
Base.showerror(io::IO, e::UnfoundVertex) = print(io, "vertex $(e.v) not found on $(e.edge)")
function findadjv(edge::WeightedUndirectedEdge, v)
    v == edge.v1 && return edge.v2
    v == edge.v2 && return edge.v1
    throw(UnfoundVertex(v, edge))
end

"""
    findadje(g, v)

Find adjacenwt edges of vertex `v` in graph `g`.
"""
findadje(g::Graph, v) =  g.adjlist[v]

# Minimum Spanning Tree (MST) problem
abstract type MSTResult end

abstract type MSTAlgo end
abstract type PrimMSTAlgo <: MSTAlgo end
struct LazyPrimMSTAlgo <: PrimMSTAlgo end
struct EagerPrimMSTAlgo <: PrimMSTAlgo end
struct KruskalMSTAlgo <: MSTAlgo end

const LazyPrimMST = LazyPrimMSTAlgo()
const EagerPrimMST = EagerPrimMSTAlgo()
const KruskalMST = KruskalMSTAlgo()
struct MST
    result::MSTResult
    algo::MSTAlgo

    function MST(g, algo)
        result = createresult(g, algo)
        findmst!(g, result, algo)
        new(result, algo)
    end
end

sumweights(mst) = sum([e.weight for e in mst.result.mstedges if !isnothing(e)])

struct LazyPrimMSTResult <: MSTResult
    mstvertices::Vector{Bool}
    mstedges::Vector{WeightedUndirectedEdge}
    crossingedges::MinHeapQueue{WeightedUndirectedEdge}

    function LazyPrimMSTResult(g)
        mstvertices = [false for _ in 1:g.numvertices]
        mstedges = WeightedUndirectedEdge[]
        crossingedges = MinHeapQueue{WeightedUndirectedEdge}()
        new(mstvertices, mstedges, crossingedges)
    end
end

"""
Eager Prim algorithm result.

- `mstedges`: `mstedges[v]` is the edge connecting v to the MST.
- `distto`: `distto[v]` is the weight of `mstedges[v]`
- `vqueue`: an indexed minimum heap queue of vertex with weight as prority.
"""
struct EagerPrimMSTResult <: MSTResult
    mstvertices::Vector{Bool}
    mstedges::Vector{Union{WeightedUndirectedEdge, Nothing}}
    distto::Vector{Float64}
    vqueue::IndexedMinHQ{Float64}

    function EagerPrimMSTResult(g)
        mstvertices = [false for _ in 1:g.numvertices]
        mstedges = Vector{Union{WeightedUndirectedEdge, Nothing}}(nothing, g.numvertices)
        distto = [Inf for _ in 1:g.numvertices]
        vqueue = IndexedMinHQ{Float64}(g.numvertices)
        new(mstvertices, mstedges, distto, vqueue)
    end
end

struct KruskalMSTResult <: MSTResult
    mstvertices::UnionFind
    mstedges::Vector{WeightedUndirectedEdge}
    edgequeue::MinHeapQueue{WeightedUndirectedEdge}

    function KruskalMSTResult(g)
        mstvertices = UnionFind(g.numvertices)
        mstedges = WeightedUndirectedEdge[]
        # create a min heap queue from edge list
        edgequeue = heapify!(getedges(g))
        new(mstvertices, mstedges, edgequeue)
    end
end

createresult(g, ::LazyPrimMSTAlgo) = LazyPrimMSTResult(g)
createresult(g, ::EagerPrimMSTAlgo) = EagerPrimMSTResult(g)
createresult(g, ::KruskalMSTAlgo) = KruskalMSTResult(g)

# Lazy Prim algorithm
function findmst!(g, result, algo::LazyPrimMSTAlgo)
    visit!(g, 1, result)
    while !isempty(result.crossingedges)
        e = popfirst!(result.crossingedges)
        iseligible(result.mstvertices, e) || continue
        push!(result.mstedges, e)
        ismarked(result.mstvertices, e.v1) || visit!(g, e.v1, result)
        ismarked(result.mstvertices, e.v2) || visit!(g, e.v2, result)
    end
end

function visit!(g, v, result::LazyPrimMSTResult)
    result.mstvertices[v] = true
    for e in findadje(g, v)
        adjv = findadjv(e, v)
        # find the adjacent vertex on the same edge
        # if the adjacent vertex is not on the tree,
        # then this edge is crossing edge
        ismarked(result.mstvertices, adjv) || push!(result.crossingedges, e)
    end
end

iseligible(mstvertices, edge::Edge) = mstvertices[edge.v1] == false || mstvertices[edge.v2] == false
ismarked(mstvertices, v) = mstvertices[v] == true

# Eager Prim Algorithm
function findmst!(g, result, ::EagerPrimMSTAlgo)
    result.distto[1] = 0.0
    insert!(result.vqueue, 1, 0.0)
    while !isempty(result.vqueue)
        minv, _ = popfirst!(result.vqueue)
        visit!(g, minv, result)
    end
end

function visit!(g, v, result::EagerPrimMSTResult)
    result.mstvertices[v] = true
    for e in findadje(g, v)
        adjv = findadjv(e, v)
        # adjv is already in the MST, ineligible edge
        ismarked(result.mstvertices, adjv) && continue
        if e.weight < result.distto[adjv]
            result.mstedges[adjv] = e
            result.distto[adjv] = e.weight
            if contains(result.vqueue, adjv)
                # if adjv exists in pq, make sure the edge weight
                # is smaller than existing value
                decreasekey!(result.vqueue, adjv, result.distto[adjv])
            else
                insert!(result.vqueue, adjv, result.distto[adjv])
            end
        end
    end
end

# Kruskal's algorithm
function findmst!(g, result, algo::KruskalMSTAlgo)
    while !isempty(result.edgequeue) && length(result.mstedges) < g.numvertices - 1
        edge = popfirst!(result.edgequeue)
        # if v1 and v2 are not connected in the subtree of MST
        # that means the edge is a non-black crossing edge with minimum
        # weight (because it's from minimum priority queue)
        # according to greedy MST algorithm.
        # this edge must be on MST of the graph
        if !isconnected(result.mstvertices, edge.v1, edge.v2)
            push!(result.mstedges, edge)
            union!(result.mstvertices, edge.v1, edge.v2)
        end
    end
end

# Shortest Paths
abstract type ShortestPathsAlgo end
abstract type ShortestPathsResult end
struct DijkstraSPAlgo <: ShortestPathsAlgo end
const DijkstraSP = DijkstraSPAlgo()

struct SP
    result::ShortestPathsResult
    algo::ShortestPathsAlgo
    function SP(g, source, algo::ShortestPathsAlgo)
        validate(g, algo)
        result = createresult(g, source, algo)
        findsp!(g, source, result, algo)
        new(result, algo)
    end
end

function validate(g::Graph{WeightedDirectedEdge}, ::DijkstraSPAlgo)
    for (v, adjedges) in enumerate(g.adjlist)
        for edge in adjedges
            edge.weight < 0 && throw(ArgumentError(
                "$(edge.weight): Dijkastra shortst paths requires non-negative weight"))
        end
    end
end
createresult(g::Graph{WeightedDirectedEdge}, source, ::DijkstraSPAlgo) = DijkstraSPResult(g, source)

mutable struct DijkstraSPResult <: ShortestPathsResult
    marked::Vector{Bool}
    edgeto::Vector{Union{Int,Nothing}}
    distto::Vector{Union{Float64,Nothing}}
    vqueue::IndexedMinHQ{Float64}
    source::Int

    function DijkstraSPResult(g::Graph{T}, source) where T<:WeightedDirectedEdge
        marked = [false for _ in 1:g.numvertices]
        mstedges = Vector{Union{T, Nothing}}(nothing, g.numvertices)
        distto = [Inf for _ in 1:g.numvertices]
        vqueue = IndexedMinHQ{Float64}(g.numvertices)
        new(marked, mstedges, distto, vqueue, source)
    end
end

function findsp!(g::Graph{WeightedDirectedEdge}, source, result, ::DijkstraSPAlgo)
    disttosource = 0.0
    result.distto[source] = disttosource
    insert!(result.vqueue, source, disttosource)
    while !isempty(result.vqueue)
        minv, _ = popfirst!(result.vqueue)
        relax!(g, minv, result)
    end
end

function relax!(g::Graph{WeightedDirectedEdge}, v, result::DijkstraSPResult)
    result.marked[v] = true
    for e in findadje(g, v)
        tov = findtov(e, v)
        result.marked[tov] && continue
        newdist = result.distto[v] + e.weight
        if newdist < result.distto[tov]
            result.edgeto[tov] = v
            result.distto[tov] = newdist
            if !contains(result.vqueue, tov)
                insert!(result.vqueue, tov, newdist)
            else
                decreasekey!(result.vqueue, tov, newdist)
            end
        end
    end
end

function findtov(e::WeightedDirectedEdge, v)
    e.v1 == v || throw(InvalidArgumentError("$v is not the from vertice"))
    e.v2
end

distanceto(sp::SP, tov) = sp.result.distto[tov]

end # module

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
cheb_mesh(T::Type, N::Integer) = (-cos.((T(π) * (0:N)) ./ N) .+ 1) ./ 2

function cheb(T::Type, N::Integer)
    if N == 0
        return (ones(T, 1), zeros(T, 1, 1))
    else
        x = cheb_mesh(T, N)
        c = [2; ones(T, N - 1, 1); 2] .* (-1) .^ (0:N)
        dx = x .- x'
        D = (c * (1 ./ c)') ./ (dx + I)
        D[diagind(D)] .-= vec(sum(D; dims=2))
        return (x, D)
    end
end

struct CollocationSystem{T} <: AbstractDynamicalSystem
    # AbstractSystem requirements
    eqs::Vector{Equation}
    states::Vector
    ps::Vector
    name::Symbol
    defaults::Dict
    # Extra fields
    n_dim::Int
    n_coll::Int
    n_mesh::Int
    t_coll::Vector{T}
    D::Matrix{T}
    D⁻¹::Matrix{T}
    ode_ps::Vector
    ode_u₀::Vector
    ode_u₁::Vector
    ode_t::Vector
end

function CollocationSystem(
    T::Type, sys::ODESystem, n_mesh, n_coll; name=gensym(nameof(sys)), free_ps=false
)
    # Turn the ODESystem into a function we can evaluate
    iv = MTK.get_iv(sys)  # independent variable
    sts = states(sys)
    ps = parameters(sys)
    f = @RuntimeGeneratedFunction(
        build_function([x.rhs for x in equations(sys)], sts, ps, iv)[1]
    )
    n_dim = length(sts)

    # Namespace the core symbols to avoid name collisions
    @variables coll₊u[1:n_mesh, 1:(n_coll + 1), 1:n_dim]
    for i in 1:(n_mesh - 1)
        coll₊u[i, end, :] .= coll₊u[i + 1, 1, :]  # continuity
    end
    @parameters coll₊ka[1:n_mesh]
    @variables coll₊t[0:1]

    # Regenerate the parameters as variables
    if free_ps
    p = [Num(Variable(nameof(par))) for par in ps]
    else
        p = ps
    end

    # Generate the collocation equations for each mesh interval
    #   D*u ~ Tp*f(u, p, t)
    # where D is the differentiation matrix, and Tp is the time interval
    # the ka variables enable mesh interval scaling (assumes sum(ka) == 1)

    t_coll_full, D_full = cheb(T, n_coll)

    # Calculate the corresponding integration matrix (think about solving D*x = f(t) with appropriate boundary conditions)
    D_full⁻¹ = [zeros(1, n_coll + 1); zeros(n_coll, 1) inv(D_full[2:end, 2:end])]

    # End points aren't used in the collocation equations
    t_coll = t_coll_full[1:(end - 1)]
    D = D_full[1:(end - 1), :]

    # Interval length
    Tp = coll₊t[2] - coll₊t[1]

    # Collocation equations
    eqns = Vector{Equation}()
    for i in 1:n_mesh
        Dx = [D * coll₊u[i, :, k] for k in 1:n_dim]
        t_int =
            coll₊t[1] .+
            Tp .* (coll₊ka[i] .* t_coll .+ sum(coll₊ka[1:(i - 1)]; init=Num(0)))
        f_vals = [
            Tp * coll₊ka[i] * Base.invokelatest(f, coll₊u[i, j, :], p, t_int[j]) for
            j in eachindex(t_int)
        ]
        append!(
            eqns, [0 ~ Dx[k][j] - f_vals[j][k] for k in eachindex(Dx), j in eachindex(f_vals)]
        )
    end

    # Bindings to the inner variables
    sts_names = nameof.(MTK.operation.(sts))
    sts_begin = [Num(Variable(st, 0)) for st in sts_names]
    append!(eqns, 0 .~ coll₊u[1, 1, :] .- sts_begin)
    sts_end = [Num(Variable(st, 1)) for st in sts_names]
    append!(eqns, 0 .~ coll₊u[end, end, :] .- sts_end)

    # Defaults
    defaults = merge(
        MTK.get_defaults(sys), Dict(MTK.value(ka) => one(T) / n_mesh for ka in coll₊ka)
    )

    new_sts = [sts_begin; sts_end; unique(vec(coll₊u)); coll₊t]
    new_ps = copy(coll₊ka)
    if free_ps
        append!(new_sts, p)
    else
        append!(new_ps, p)
    end

    return CollocationSystem(
        eqns,  # eqs
        MTK.value.(new_sts),  # states
        MTK.value.(new_ps),  # ps
        name,  # name
        defaults,  # defaults
        n_dim,  # n_dim
        n_coll,  # n_coll
        n_mesh,  # n_mesh
        t_coll_full,  # t_coll
        D_full,  # D
        D_full⁻¹,  # D⁻¹
        p,  # ode_ps
        sts_begin,  # ode_u₀
        sts_end,  # ode_u₁
        coll₊t,  # ode_t
    )
end

function CollocationSystem(sys::ODESystem, args...; kwargs...)
    return CollocationSystem(Float64, sys, args...; kwargs...)
end

"""
    $SIGNATURES

Add equations that are a function of the boundary states and times to `coll`. The input `bc`
is assumed to be a function of the form `(u₀, t₀, u₁, t₁, p) -> equation` where `u₀` and
`u₁` are state variables at the start and end of the interval respectively, `t₀` and `t₁`
are the start and end times, and `p` are the underlying ODE system parameters.

The value returned from `bc` can be an equation or a vector of equations.

`u₀` and `u₁` have the states in the same order as the underlying ODE system.

## Example - periodic boundary conditions

```julia
add_boundarycondition!(coll, (u₀, t₀, u₁, t₁, p) -> u₀ .~ u₁)
```

Note that `.~` is used since `u₀` and `u₁` are vectors of variables.
"""
function add_boundarycondition!(coll::CollocationSystem, bc)
    eqn = bc(coll.ode_u₀, coll.ode_t[1], coll.ode_u₁, coll.ode_t[2], coll.ode_ps)
    if eqn isa Equation
        push!(coll.eqs, eqn)
    elseif eqn isa Vector{Equation}
        append!(coll.eqs, eqn)
    elseif eqn isa Num
        push!(coll.eqs, 0 ~ eqn)
    elseif eqn isa Vector{Num}
        append!(coll.eqs, 0 .~ eqn)
    else
        throw(ArgumentError("Input function should return a (vector of) Equation or Num"))
    end
    return coll
end

function generate_initialconditions(
    coll::CollocationSystem{T}, sol, tspan; mesh=nothing
) where {T}
    # Mesh intervals
    ka = mesh === nothing ? ones(T, coll.n_mesh) ./ coll.n_mesh : mesh
    @assert length(ka) == coll.n_mesh
    # Generate mesh with representation points
    Tp = tspan[end] - tspan[begin]
    t_coll = coll.t_coll[1:(end - 1)]
    t_int = similar(t_coll)
    # Determine the initial conditions
    ics = Vector{Pair{typeof(Variable(:x)), T}}()
    sizehint!(ics, (coll.n_mesh*coll.n_coll + 1)*coll.n_dim)
    basename = _strip_toplevel(MTK.renamespace(coll.name, ""))  # need to strip the top-level namespace since it isn't used in solve
    basename_u = Symbol(basename, :coll₊u)
    for i in eachindex(ka)
        t_int .= tspan[begin] .+ Tp .* (ka[i] .* t_coll .+ sum(ka[1:(i - 1)]))
        for j in eachindex(t_int)
            sts = sol(t_int[j])
            for k in eachindex(sts)
                push!(ics, Variable(basename_u, i, j, k) => sts[k])
            end
        end
    end
    sts = sol(tspan[end])
    for k in eachindex(sts)
        push!(ics, Variable(basename_u, coll.n_mesh, coll.n_coll + 1, k) => sts[k])
    end
    # User-defined state variables
    for (var, val) in zip(coll.ode_u₁, sts)  # assumes sts = sol(tspan[end])
        push!(ics, MTK.value(var) => val)
    end
    for (var, val) in zip(coll.ode_u₀, ics)  # FRAGILE: assumes first elements of ics correspond to sol(tspan[begin])
        push!(ics, MTK.value(var) => val[2])
    end
    # Start/end times
    push!(ics, Variable(Symbol(basename, :coll₊t), 0) => tspan[begin])
    push!(ics, Variable(Symbol(basename, :coll₊t), 1) => tspan[end])
    return ics
end

# CollocationSystem does not allow embedded systems (it's not meaningful)
ModelingToolkit.get_systems(::CollocationSystem) = NonlinearSystem[]

# Convert a CollocationSystem to a NonlinearSystem
function ModelingToolkit.NonlinearSystem(coll::CollocationSystem)
    return NonlinearSystem(
        MTK.get_eqs(coll),
        MTK.get_states(coll),
        MTK.get_ps(coll);
        name=nameof(coll),
        defaults=MTK.get_defaults(coll),
    )
end


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

"""
    @with_kw mutable struct Neighborhood

Determines the size of the neighborhood, \$k\$ is adjusted as a callback, and it is intended to affect previously inserted vertices.
The neighborhood is designed to consider two components \$k=in+out\$, i.e. _in_coming and _out_going edges for each vertex.
- The \$out\$ size is computed as \$minsize + \\log(logbase, n)\$ where \$n\$ is the current number of indexed elements; this is computed searching
for \$out\$  elements in the current index.
- The \$in\$ size is unbounded.
- reduce is intended to postprocess neighbors (after search process, i.e., once out edges are computed); do not change \$k\$ but always must return a copy of the reduced result set.

Note: Set \$logbase=Inf\$ to obtain a fixed number of \$in\$ nodes; and set \$minsize=0\$ to obtain a pure logarithmic growing neighborhood.

"""
@with_kw struct Neighborhood{Reduction<:NeighborhoodReduction}
    logbase::Float32 = 2
    minsize::Int32 = 2
    reduce::Reduction = SatNeighborhood()
end

Base.copy(N::Neighborhood; logbase=N.logbase, minsize=N.minsize, reduce=copy(N.reduce)) =
    Neighborhood(; logbase, minsize, reduce)

neighborhoodsize(N::Neighborhood, index::SearchGraph) = ceil(Int, N.minsize + log(N.logbase, length(index)))

"""
    find_neighborhood(index::SearchGraph{T}, item, neighborhood, pools; hints=index.hints)

Searches for `item` neighborhood in the index, i.e., if `item` were in the index whose items should be its neighbors (intenal function).
`res` is always reused since `reduce` creates a new KnnResult from it (a copy if `reduce` in its simpler terms)

# Arguments
- `index`: The search index.
- `item`: The item to be inserted.
- `neighborhood`: A [`Neighborhood`](@ref) object that describes how to compute item's neighborhood.
- `pools`: Cache pools to be used
- `hints`: Search hints
"""
function find_neighborhood(index::SearchGraph, item, neighborhood::Neighborhood, pools::SearchGraphPools; hints=index.hints)
    n = length(index)
    if n > 0
        ksearch = neighborhoodsize(neighborhood, index)
        res = getknnresult(ksearch, pools)
        search(index.search_algo, index, item, res, hints, pools)
        neighborhoodreduce(neighborhood.reduce, index, item, res, pools)
    else
        Int32[]
    end
end

"""
    push_neighborhood!(index::SearchGraph, item, neighbors, callbacks; push_item=true)

Inserts the object `item` into the index, i.e., creates an edge for each item in `neighbors` (internal function)

# Arguments

- `index`: The search index to be modified.
- `item`: The item that will be inserted.
- `neighbors`: An array of indices that will be connected to the new vertex.
- `callbacks`: A [`SearchGraphCallbacks`] object (callback list) that will be called after some insertions
- `push_item`: Specifies if the item must be inserted into the internal `db` (sometimes is already there like in [`index!`](@ref))
"""
function push_neighborhood!(index::SearchGraph, item, neighbors, callbacks; push_item=true)
    push_item && push!(index.db, item)
    push!(index.links, neighbors)
    push!(index.locks, Threads.SpinLock())
    n = length(index)
    n == 1 && return
    ## vstate = getvisitedvertices(index)
    @inbounds for id in neighbors
        push!(index.links[id], n)  # sat push?
    end

    callbacks !== nothing && execute_callbacks(callbacks, index)

    if index.verbose && length(index) % 100_000 == 0
        println(stderr, "added n=$(length(index)), neighborhood=$(length(neighbors)), $(string(index.search_algo)), $(Dates.now())")
    end
end

const GlobalSatKnnResult = [KnnResult(1)]
function __init__neighborhood()
    for _ in 2:Threads.nthreads()
        push!(GlobalSatKnnResult, KnnResult(1))
    end
end

"""
    SatNeighborhood()

New items are connected with a small set of items computed with a SAT like scheme (**cite**).
It starts with `k` near items that are reduced to a small neighborhood due to the SAT partitioning stage.
"""
struct SatNeighborhood <: NeighborhoodReduction end
Base.copy(::SatNeighborhood) = SatNeighborhood()

"""
    DistalSatNeighborhood()

New items are connected with a small set of items computed with a Distal SAT like scheme (**cite**).
It starts with `k` near items that are reduced to a small neighborhood due to the SAT partitioning stage but in reverse order of distance.
"""
struct DistalSatNeighborhood <: NeighborhoodReduction end
Base.copy(::DistalSatNeighborhood) = DistalSatNeighborhood()

"""
    struct IdentityNeighborhood

It does not modifies the given neighborhood
"""
struct IdentityNeighborhood <: NeighborhoodReduction end
Base.copy(::IdentityNeighborhood) = IdentityNeighborhood()

## functions

function sat_should_push(sat_neighborhood::T, index, item, id, dist, near::KnnResult) where T
    @inbounds obj = index[id]
    dist = dist < 0.0 ? evaluate(index.dist, item, obj) : dist
    push!(near, 0, dist)

    @inbounds for linkID in sat_neighborhood
        d = evaluate(index.dist, index[linkID], obj)
        push!(near, linkID, d)
    end

    argmin(near) == 0
end

function neighborhoodreduce(::IdentityNeighborhood, index::SearchGraph, item, res, pools::SearchGraphPools)
    copy(res.id)
end

"""
    reduce(sat::DistalSatNeighborhood, index::SearchGraph, item, res, pools)

Reduces `res` using the DistSAT strategy.
"""
@inline function neighborhoodreduce(::DistalSatNeighborhood, index::SearchGraph, item, res, pools::SearchGraphPools)
    N = Vector{Int32}(undef, 2)
    resize!(N, 0)
    @inbounds for i in length(res):-1:1  # DistSat => works a little better but also produces a bit larger neighborhoods
        id, dist = res[i]
        sat_should_push(N, index, item, id, dist, getsatknnresult(pools)) && push!(N, id)
    end

    N
end

@inline function neighborhoodreduce(::SatNeighborhood, index::SearchGraph, item, res, pools::SearchGraphPools)
    N = Vector{Int32}(undef, 2)
    resize!(N, 0)
    @inbounds for (id, dist) in res  # DistSat => works a little better but also produces a bit larger neighborhoods
        sat_should_push(N, index, item, id, dist, getsatknnresult(pools)) && push!(N, id)
    end

    N
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
export Point, print

"""
    Point(coords::SortedDict{Variable, Number})

Define a mathematical point, that is a pairing of variables and numbers.

### Attributes
- `coords` : a dictionary associating a `Variable` to a a number.

### Exemple
```julia
julia > a, b = Variable("a", Complex), Variable("b", Complex)
julia > pt = Point(SortedDict(a=>π, b=>e+7im))
```
"""
struct Point
    coords::SortedDict{Variable, Number}
    isdense::Bool

    function Point(dict::SortedDict; isdense=false)
        dict_pt = SortedDict{Variable, Number}()
        for (var, val) in dict
            if !isa(var, Variable) || !isa(val, Number)
                error(LOGGER, "Point(): Expected pair of (Variable, Number), got ($var, $val) of type ($typeof(var), $typeof(val)) instead.")
            end
            if isbool(var)
                booled = 0
                if val != 0
                    booled = Int((val/abs(val) + 1) / 2)
                end
                (booled ≠ val) && warn(LOGGER, "Point(): $var is $(var.kind), provided value is $val, $booled stored.")
                add_to_dict!(dict_pt, var, booled, isdense = isdense)
            elseif isreal(var)
                realed = real(val)
                (realed ≠ val) && warn(LOGGER, "Point(): $var is $(var.kind), provided value is $val, $realed stored.")
                add_to_dict!(dict_pt, var, realed, isdense = isdense)
            else
                add_to_dict!(dict_pt, var, val, isdense = isdense)
            end
        end
        return new(dict_pt, isdense)
    end
end

#############################
## Constructors
#############################
Point() = Point(SortedDict{Variable, Number}())

function Point(vars::Array{Variable}, vals::Array{<:Number})
    if length(vars) != length(vals)
        error(LOGGER, "Point(): input arrays must have same size.")
    end

    pt = Point()
    for i=1:length(vars)
        var, val = vars[i], vals[i]
        if isreal(var) val = real(val) end
        if isbool(var) && val != 0
            val = Int((val/abs(val)+1)/2)
        end
        add_coord!(pt, vars[i], val)
    end
    return pt
end

#############################
## Equality, hash
#############################
==(pt1::Point, pt2::Point) = (pt1.coords == pt2.coords) && (pt1.isdense == pt2.isdense)
!=(pt1::Point, pt2::Point) = !(pt1 == pt2)

function setindex!(pt::Point, val::Number, var::Variable)
    if val!=0 || pt.isdense
        pt.coords[var] = val
    end
    return
end

#############################
## Print
#############################
function Base.print(io::IO, pt::Point)
    maxvarlength = maximum(map(x->length(x.name), keys(pt)))
    for (var, val) in pt
        print(io, var.name); print(io, " "^(maxvarlength-length(var.name)))
        println(io, " ", val)
    end
end


text:
# In this file, you'll fill in some blanks, and the associated `@test` lines thereafter will check whether
# you've done it correctly. We'll cover testing systematically later in the course.
# For now, know that:
#   - `@test` checks whether the following expression is `true`
#   - `@testset` wraps a collection of tests and gathers statistics. Here, I'm using it just so you can
#     run all the tests even if you have some that fail. (Without `@testset`, execution stops at the first
#     failing test)
#
# Example: if you see
#
#     # Set `a` to be an array with two items, `1` and `7`
#     a =
#     @test length(a) == 2 && a[1] == 1 && a[2] == 7
#
# notice the blank after `a =`: you would edit this to
#
#     # Set `a` to be an array with two items, `1` and `7`
#     a = [1, 7]
#     @test length(a) == 2 && a[1] == 1 && a[2] == 7
#
# Ideally, running the test file passes all tests.
#
# Also, it may be helpful to note you can comment out a block of code using matching `#=`, `=#`, but
# in VS Code it's often just as easy to select the whole range and use Ctrl-/ (⌘/ on Mac) to toggle between
# commented and uncommented. A good working strategy is to start by commenting out all the tests and then
# uncomment them one-by-one as you address them. That way, you can run the test script for each new test
# and check whether all work so far.

using Test

# Try the "gimmee" now!

# Set `a` to be an array with two items, `1` and `7`
a =
@test length(a) == 2 && a[1] == 1 && a[2] == 7

# OK, now let's get serious!

@testset "Learning Julia 1" begin
    # Create a range of all integers from 2 to 5, inclusive
    r =
    @test isa(r, AbstractUnitRange) && first(r) == 2 && last(r) == 5

    # Create an array that stores `Float64`s but otherwise holds the same values as `r` above
    rf =
    @test eltype(rf) == Float64 && rf == r

    # Create a matrix equal to a 2x2 identity matrix (1 on the diagonal, zero on the corners)
    M =
    @test size(M) == (2,2) && M[1,1] == M[2,2] == 1 && M[1,2] == M[2,1] == 0

    # Create a vector of `String`s that has 3 undefined elements, then make the middle one equal to "Julia"
    vstr =
    # you will need one more line here
    @test length(vstr) == 3 && eltype(vstr) == String && !isassigned(vstr, 1) && !isassigned(vstr, 3) && vstr[2] == "Julia"

    # Create two ranges, `row` and `col`: row should go from 1 to 3 (integers) and `col` should contain the even numbers between 2 and 8.
    # Use broadcasting to form a matrix of their product.
    row =
    col =
    m =
    @test isa(row, AbstractUnitRange) && isa(col, AbstractRange) && first(row) == 1 && last(row) == 3 &&
          first(col) == 2 && last(col) == 8 && step(col) == 2
    @test m ==
    [ 2   4   6   8;
      4   8  12  16;
      6  12  18  24
    ]

    # Create a tuple with a `String`, an `Int`, and a `Float64` (of your choice) in that order
    t =
    @test isa(t, Tuple) && isa(t[1], String) && isa(t[2], Int) && isa(t[3], Float64)

    # Create a NamedTuple with fields `make` and `model` with values "Honda" and "Odyssey", respectively
    nt =
    @test nt.make == "Honda"
    @test nt.model == "Odyssey"
    @test isa(nt, NamedTuple)

    # Split this string into words
    str = "Advanced scientific computing"
    sstr =
    @test length(sstr) == 3 && sstr[1] == "Advanced" && sstr[2] == "scientific" && sstr[3] == "computing"

    # Now things will start getting a bit (just a bit) harder
    # If it's useful to write a `for` loop, feel free to look ahead to https://en.wikibooks.org/wiki/Introducing_Julia/Controlling_the_flow#For_loops_and_iteration
    # Some problems can be solved in a single line using a comprehension or generator, https://en.wikibooks.org/wiki/Introducing_Julia/Controlling_the_flow#Comprehensions

    # Alphabetize the words in this string, omitting any that contain a 'u'
    # Do not type the answer directly, generate it by exploiting Julia's string-manipulation functions
    animals = "dog cat opossum feline antelope chimp octopus salamander"
    aanimals =  # feel free to use multiple lines to solve this; `aanimals` should hold your final answer
    @test aanimals == ["antelope", "cat", "chimp", "dog", "feline", "salamander"]

    # Compute the sum of the ascii codes in `animals`, excluding spaces
    sascii =
    @test sascii == 5257

    # What is the most common letter, excluding spaces, in `animals` above?
    d =     # this variable should be a dictionary
    mostcommon =    # compute this from `d`
    @test keytype(d) == Char && valtype(d) == Int
    @test d['l'] == 3
    @test mostcommon == 'o'

    # Determine all the unique characters in `animals` using a Set (including the space ' ' this time)
    chars =
    @test isa(chars, Set) && eltype(chars) == Char && length(chars) == 18 && 't' ∈ chars && 'z' ∉ chars
end


text:
"""
    addNPDGeoRegions() -> nothing

Checks for the three GeoRegions (GPM, TRMM, TRMMLSM) required by NASAPrecipitation.jl, and adds them if they do not exist.
"""
function addNPDGeoRegions()
	@info "$(modulelog()) - Checking to see if GeoRegions required by NASAPrecipitation.jl have been added to the list of available GeoRegions"
    disable_logging(Logging.Warn)
	if !isGeoRegion("IMERG",throw=false) ||
	    !isGeoRegion("TRMM",throw=false) ||
	    !isGeoRegion("TRMMLSM",throw=false)
        disable_logging(Logging.Debug)
        @info "$(modulelog()) - At least one of the required three GeoRegions (IMERG, TRMM, TRMMLSM) has not been added, proceeding to add them again ..."
	    addGeoRegions(joinpath(@__DIR__,"NPDGeoRegions.txt"))
    else
        disable_logging(Logging.Debug)
        @info "$(modulelog()) - All of the required three GeoRegions (IMERG, TRMM, TRMMLSM) have been added"
	end
	return
end

## DateString Aliasing
yrmo2dir(date::TimeType) = Dates.format(date,dateformat"yyyy/mm")
yrmo2str(date::TimeType) = Dates.format(date,dateformat"yyyymm")
yr2str(date::TimeType)   = Dates.format(date,dateformat"yyyy")
ymd2str(date::TimeType)  = Dates.format(date,dateformat"yyyymmdd")

function imergrawfiles()

    fnc = Vector{AbstractString}(undef,48);

    for ii = 1 : 48

        hr = (ii-1)/2; mi = mod(hr,1); hr = hr - mi;
        hr = @sprintf("%02d",hr);
        id = @sprintf("%04d",(ii-1)*30);

        if iszero(mi)
              fnc[ii] = "S$(hr)0000-E$(hr)2959.$(id)"
        else; fnc[ii] = "S$(hr)3000-E$(hr)5959.$(id)"
        end

    end

    return fnc

end

function gpmlonlat()
    lon = convert(Array,-179.95:0.1:179.95)
    lat = convert(Array,-89.95:0.1:89.95)
    return lon,lat
end

function trmmlonlat(;full::Bool=false)

    lon = convert(Array,-179.875:0.25:179.875)
    if !full
		  lat = convert(Array,-49.875:0.25:49.875)
	else; lat = convert(Array,-89.875:0.25:89.875)
	end

    return lon,lat
end

function ncoffsetscale(data::AbstractArray{<:Real},init=0)

    dmax = init
    dmin = init
    for ii = 1 : length(data)
        dataii = data[ii]
        if !isnan(dataii)
            if dataii > dmax; dmax = dataii end
            if dataii < dmin; dmin = dataii end
        end
    end

    scale = (dmax-dmin) / 65533;
    offset = (dmax+dmin-scale) / 2;

    return scale,offset

end

function real2int16!(
    outarray :: AbstractArray{Int16},
    inarray  :: AbstractArray{<:Real},
    scale    :: Real,
    offset   :: Real
)

    for ii = 1 : length(inarray)

        idata = (inarray[ii] - offset) / scale
        if isnan(idata)
              outarray[ii] = -32767
        else; outarray[ii] = round(Int16,idata)
        end

    end

    return

end

function extractregionlsm!(
	outarray :: Array{<:Real,2},
	inarray  :: Array{<:Real,2},
	ginfo	 :: RectGrid
)

	iglon = ginfo.ilon; nglon = length(iglon)
	iglat = ginfo.ilat; nglat = length(iglat)
	for ilat = 1 : nglat, ilon = 1 : nglon
		outarray[ilon,ilat] = inarray[iglat[ilat],iglon[ilon]]
	end

end

function extractregionlsm!(
	outarray :: Array{<:Real,2},
	inarray  :: Array{<:Real,2},
	ginfo	 :: PolyGrid
)

	iglon = ginfo.ilon; nglon = length(iglon)
	iglat = ginfo.ilat; nglat = length(iglat)
	for ilat = 1 : nglat, ilon = 1 : nglon
		outarray[ilon,ilat] = inarray[iglat[ilat],iglon[ilon]] * ginfo.mask[ilon,ilat]
	end

end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

function add_boundary_condition!(case::Simulation, bc::NeumannBC)
    push!(case.neumann_boundary_conditions, bc)
end

"""
Add node to model and renumber for output
"""
function add_node!(model::Model, index::Union{Int64, String},
    coords::Vector{Float64})
    node = Node(index, coords)
    model.nodes[index] = node
    node_number = length(model.nodes)
    model.renum_nodes[index] = node_number
end

function add_boundary_condition!(case::Simulation, bc::DirichletBC)
    push!(case.dirichlet_boundary_conditions, bc)
end

function add_solver!(case::Simulation, solver)
    case.solver = solver
end

function add_material!(model::Model, set_name::String, material::Material)
    element_set = model.elsets[set_name]
    set_ids = element_set.elements
    for each in set_ids
        element = model.elements[each]
        element.material = material
    end
    element_set.material = material
end

function add_element!(model::Model, idx::Union{Int64, String},
    eltype::Symbol, node_ids::Vector{Int64})

    element = Element(idx, node_ids, eltype)
    model.elements[idx] = element
end

function add_element_set!(model::Model, elset::ElementSet)
    name = elset.name
    model.elsets[name] = elset
end

function add_element_set!(model::Model, name::String, ids::Vector{Int64})
    elset = ElementSet(name, ids)
    model.elsets[name] = elset
end

function add_node_set!(model::Model, name::String, ids::Vector{Int64})
    nset = NodeSet(name, ids)
    model.nsets[name] = nset
end

function add_element_set!(model::Model, name::String,
    elements::Vector{Element})
    elset = ElementSet(name, elements)
model.elsets[name] = elset
end

function add_element_set!(case::Simulation, name::String)
    push!(case.sets, name)
end

function add_simulation!(model::Model, name::String, case::Simulation)
    model.load_cases[name] = case
end
#function Base.convert{T<:AbstractFloat}(::Type{Node}, data::Vector{T})
#    Node(data)
#end:q
#
#
#function set_material!{S <: AbstractString}(model::Model, material::Material, set_name::S)
#    model.sets[set_name].material = material
#end
#
#
#function get_element_set(model::Model, set_name::String)
#    get_set = model.sets[set_name]
#    isa(get_set, ElementSet) ? get_set : err("Found set $(set_name)
#    but it is not a ElementSet. Check if you have used
#    dublicate set names.")
#end
#
#function add_boundary_condition!{B <: BoundaryCondition}(case::Simulation, bc::B)
#    push!(case.boundary_conditions, bc)
#end
#
#function add_boundary_condition!{B <: BoundaryCondition}(case::Simulation, bc::Vector{B})
#    map(x-> push!(case.boundary_conditions, x), bc)
#end
#
#function add_loadcase!(model::Model, case::Simulation)
#    push!(model.load_cases, case)
#end
#
#function add_loadcase!{T<:Simulation}(model::Model, case::Vector{T})
#    map(x-> push!(model.load_cases, x), case)
#end
#
#"""
#"""
#function build_core_elements()
#
#end
#
#"""
#"""
#function set_core_element_material()
#
#end
#
#"""
#element_has_type( ::Type{Val{:C3D4}}) = Tet4
#"""
#function create_problems()
#
#end
#
#"""
#"""
#function add_problems!()
#
#end
#
#"""
#Get set from model
#"""
#function get_set{S <: AbstractString}(model::Model, name::S)
#    try
#        return model.set[name]
#    catch
#        err("Given set: $(name) does not exist in Model")
#    end
#end
#
#"""
#"""
#function push!(model::Model, element::Element)
#    push!(model.elements, element...)
#end
#
#"""
#"""
#function push!(model::Model, elements::Vector{Element})
#    push!(model.element, elements...)
#end
#
#"""
#"""
#function push!(model::Model, node::Node)
#    push!(model.nodes, node)
#end
#
#function push!(model::Model, nodes::Vector{Node})
#    push!(mode, nodes...)
#end
#
#"""
#"""
#function add_set!(model::Model, set::Union{NodeSet, ElementSet})
#    model.sets[set.name] = set
#end
#
#function add_set!{S <: AbstractString}(model::Model, ::Type{Val{:NSET}},
#    name::S, ids::Vector{Integer})
#    new_set = NodeSet(name, ids)
#    add_set!(model, new_set)
#end
#
#function add_set!{S <: AbstractString}(model::Model, ::Type{Val{:ELSET}},
#    name::S, ids::Vector{Integer})
#    new_set = ElementSet(name, ids)
#    add_set!(model, new_set)
#end


text:
# Example for HINT layer (Kruse et al, 2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test

# X dimensions
nx1 = 32
nx2 = 32
nx_channel = 32
nx_hidden = 64
batchsize = 2

# Y dimensions
ny1 = 128
ny2 = 64
ny_channel = 1
ny_hidden = 32

# Linear operator
A = randn(Float32, ny1*ny2, nx1*nx2)

# Input image
X = glorot_uniform(nx1, nx2, nx_channel, batchsize)
Y = glorot_uniform(ny1, ny2, ny_channel, batchsize)

# Conditional i-SLIM layer
CI = ConditionalLayerSLIM(nx1, nx2, nx_channel, nx_hidden, ny1, ny2, ny_channel, ny_hidden, batchsize; type="affine")

# Forward/inverse
Zx, Zy, logdet = CI.forward(X, Y, A)
X_, Y_ = CI.inverse(Zx, Zy, A)

@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)
@test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-3)

# Forward/backward
Zx, Zy, logdet = CI.forward(X, Y, A)
X_, Y_ = CI.backward(0f0*Zx, 0f0*Zy, Zx, Zy, A)[3:4]

@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)
@test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-3)

# Forward/inverse Y-lane only
Zy = CI.forward_Y(Y)
Y_ = CI.inverse_Y(Zy)

@test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-6)


text:
@testset "Partition" begin
  @testset "Directional" begin
    # merge operation does not produce NaN
    dir = (0.286788, -0.496732, -0.819152)
    𝒟 = georef(CSV.File(joinpath(datadir,"nanlags.csv")), (:X,:Y,:Z))
    γ = DirectionalVariogram(dir, 𝒟, :Cu, dtol=45, maxlag=150, nlags=20)
    x, y, n = values(γ)
    @test !any(isnan.(x))
    @test !any(isnan.(y))
    @test !any(isnan.(n))

    # directional variogram and known anisotropy ratio
    img = readdlm(joinpath(datadir,"anisotropic.tsv"))
    sdata = georef((z=img,))
    γhor = DirectionalVariogram((1.,0.), sdata, :z, maxlag=50.)
    γver = DirectionalVariogram((0.,1.), sdata, :z, maxlag=50.)
    γₕ = fit(GaussianVariogram, γhor)
    γᵥ = fit(GaussianVariogram, γver)
    @test range(γₕ) / range(γᵥ) ≈ 3. atol=.1

    if visualtests
      p1 = plot(γhor, showbins=false, label="horizontal")
      plot!(γver, showbins=false, label="vertical")
      p2 = plot(γₕ, 0., 50., label="horizontal")
      plot!(γᵥ, 0., 50., label="vertical")
      plt = plot(p1, p2, layout=(2,1))
      @test_reference "data/directional.png" plt
    end
  end

  @testset "Planar" begin
    # directional equals planar rotated by 90 degrees in 2D
    img = readdlm(joinpath(datadir,"anisotropic.tsv"))
    sdata = georef((z=img,))
    γ₁ = PlanarVariogram((0.,1.), sdata, :z, maxlag=50.)
    γ₂ = DirectionalVariogram((1.,0.), sdata, :z, maxlag=50.)
    x₁, y₁, n₁ = values(γ₁)
    x₂, y₂, n₂ = values(γ₂)
    @test x₁ == x₂
    @test y₁ ≈ y₂
    @test n₁ == n₂
    γ₁ = PlanarVariogram((1.,0.), sdata, :z, maxlag=50.)
    γ₂ = DirectionalVariogram((0.,1.), sdata, :z, maxlag=50.)
    x₁, y₁, n₁ = values(γ₁)
    x₂, y₂, n₂ = values(γ₂)
    @test x₁ == x₂
    @test y₁ ≈ y₂
    @test n₁ == n₂

    # planar variogram and known anisotropy ratio
    γhor = PlanarVariogram((0.,1.), sdata, :z, maxlag=50.)
    γver = PlanarVariogram((1.,0.), sdata, :z, maxlag=50.)
    γₕ = fit(GaussianVariogram, γhor)
    γᵥ = fit(GaussianVariogram, γver)
    @test range(γₕ) / range(γᵥ) ≈ 3. atol=.1

    if visualtests
      p1 = plot(γhor, showbins=false, label="horizontal")
      plot!(γver, showbins=false, label="vertical")
      p2 = plot(γₕ, 0., 50., label="horizontal")
      plot!(γᵥ, 0., 50., label="vertical")
      plt = plot(p1, p2, layout=(2,1))
      @test_reference "data/planar.png" plt
    end
  end
end


text:
function pointidx(i::Integer, j::Integer, order)
    # I'm sorry
    # Quad: https://gitlab.kitware.com/nick.laurenson/vtk/-/blob/d8aa5b89c622cf04b3322112fda420afc9d2a16d/Common/DataModel/vtkLagrangeQuadrilateral.cxx#L558
    # Hex: https://gitlab.kitware.com/nick.laurenson/vtk/blob/d8aa5b89c622cf04b3322112fda420afc9d2a16d/Common/DataModel/vtkLagrangeHexahedron.cxx#L734
    # Implemented for quad here
    is_i_boundary = i == 0 || i == order[1]
    is_j_boundary = j == 0 || j == order[2]

    n_boundary = is_i_boundary + is_j_boundary

    i_pos = i > 0
    j_pos = j > 0

    if n_boundary == 2
        # Vertex DOF
        return (i_pos ? (j_pos ? 2 : 1) : (j_pos ? 3 : 0))
    end

    offset = 4
    if n_boundary == 1
        # Edge DOF
        if !is_i_boundary
            # On i axis
            return (i - 1) +
            (j_pos ? order[1] - 1 + order[2] - 1 : 0) +
            offset
        end
        if !is_j_boundary
            # On j axis
            return (j - 1) +
            (i_pos ? order[1] - 1 : 2 * (order[1] - 1) + order[2] - 1) +
            offset
        end
    end

      offset += 2 * (order[1] - 1 + order[2] - 1)
      # n_boundary == 0 -> Face DOF
      return offset +
        (i - 1) + (order[1] - 1) * (
          (j - 1))

end

function get_plot_points(order)
    dim = 2
    n_points_x = (order + 1)
    n_points = (order + 1) * (order + 1)
    n_points_x, n_points
    points = zeros(dim, n_points)
    xs = (gausslobatto(n_points_x)[1] .+ 1) ./ 2
    ys = xs
    for xidx=1:n_points_x
        for yidx=1:n_points_x
            idx = pointidx(xidx-1,yidx-1,(order, order)) + 1
            x = xs[xidx]
            y = ys[yidx]
            points[:,idx] = [x,y]
        end
    end
    return transpose(points), size(points, 2)
end

"""
    VTKPlotter(eq::Equation, scenario::Scenario, grid::Grid,
               filename::String)

Initialize a VTKPlotter for equation `eq` and scenario `scenario`,
defined on `grid`.
Output name is `filename`.
"""
mutable struct VTKPlotter
    eq::Equation
    scenario::Scenario
    grid::Grid
    filename::String
    collection::WriteVTK.CollectionFile
    plot_counter::Int64
    vtkcells::Array{MeshCell,1}
    vtkpoints::Array{Float64,2}
    cellpoints::Array{Float64,2}

    function VTKPlotter(eq::Equation, scenario::Scenario, grid::Grid,
                        filename::String)
        collection = paraview_collection(filename)
        plot_counter = 0

        order = 1
        cellpoints, num_points_per_cell = get_plot_points(order)

        vtkcells = Array{MeshCell,1}(undef, length(grid.cells))
        vtkpoints = Array{Float64, 2}(undef, (2, length(grid.cells)*num_points_per_cell))

        for (i,cell) in enumerate(grid.cells)
            offset = cell.center - [cell.size[1]/2, cell.size[2]/2]
            start = (i-1) * num_points_per_cell
            for j=1:num_points_per_cell
                vtkpoints[:, start+j] = offset .+ cellpoints[j,:] .* cell.size
            end
            vtkcells[i] = MeshCell(VTKCellTypes.VTK_LAGRANGE_QUADRILATERAL, Array(start+1:start+num_points_per_cell))
        end

        new(eq, scenario, grid, filename, collection, plot_counter, vtkcells, vtkpoints, cellpoints)
        end

end

function evaluate_dof_points(eq, grid, points, dofidx)
    num_points = size(points, 1)
    eval_data = Array{Float64, 2}(undef, (length(grid.cells)*num_points, 1))
    offset = 0
    for cell in grid.cells
        for i=1:num_points
            eval_data[offset+i,:] .= grid.dofs[dofidx, cell.dataidx]
        end
        offset += num_points
    end
    eval_data
end

"""
    plot(plotter::VTKPlotter)

Write output with `plotter` for timestep.
"""
function plot(plotter::VTKPlotter)
    grid = plotter.grid
    eq = plotter.eq

    currentfile = @sprintf("%s_%d", plotter.filename, plotter.plot_counter)
    plotter.plot_counter += 1
    vtkfile = vtk_grid(currentfile, plotter.vtkpoints, plotter.vtkcells)

    for var=1:get_ndofs(eq)
        evaluated_data = evaluate_dof_points(eq, grid, plotter.cellpoints, var)
        vtk_point_data(vtkfile, evaluated_data, string(get_variable_name(eq, var)))
    end
    collection_add_timestep(plotter.collection, vtkfile, grid.time)
end

"""
    save(plotter::VTKPlotter)

Save final output file for `plotter`.
"""
function save(plotter::VTKPlotter)
    vtk_save(plotter.collection)
end


text:
inv2(A::Array{Float64,2}) = 1/(A[1]*A[4]-A[2]*A[3])*[A[4] -A[2];-A[3] A[1]]

inv3(A::Array{Float64,2}) = 1/det(A)*[(A[2,2]*A[3,3]-A[3,2]*A[2,3]) (A[1,3]*A[3,2]-A[3,3]*A[1,2]) (A[1,2]*A[2,3]-A[2,2]*A[1,3])
                    (A[2,3]*A[3,1]-A[3,3]*A[2,1]) (A[1,1]*A[3,3]-A[3,1]*A[1,3]) (A[1,3]*A[2,1]-A[2,3]*A[1,1])
                    (A[2,1]*A[3,2]-A[3,1]*A[2,2]) (A[1,2]*A[3,1]-A[3,2]*A[1,1]) (A[1,1]*A[2,2]-A[2,1]*A[1,2])]

function inv4(m::Array{Float64,2})
    iA = zeros(Float64,4,4)
    iA[1] = (((m[6] * m[11] * m[16] - m[6] * m[12] * m[15]) - m[10] * m[7] * m[16]) + m[10] * m[8] * m[15] + m[14] * m[7] * m[12]) - m[14] * m[8] * m[11]
    iA[5] = (((-(m[5]) * m[11] * m[16] + m[5] * m[12] * m[15] + m[9] * m[7] * m[16]) - m[9] * m[8] * m[15]) - m[13] * m[7] * m[12]) + m[13] * m[8] * m[11]
    iA[9] = (((m[5] * m[10] * m[16] - m[5] * m[12] * m[14]) - m[9] * m[6] * m[16]) + m[9] * m[8] * m[14] + m[13] * m[6] * m[12]) - m[13] * m[8] * m[10]
    iA[13] = (((-(m[5]) * m[10] * m[15] + m[5] * m[11] * m[14] + m[9] * m[6] * m[15]) - m[9] * m[7] * m[14]) - m[13] * m[6] * m[11]) + m[13] * m[7] * m[10]
    iA[2] = (((-(m[2]) * m[11] * m[16] + m[2] * m[12] * m[15] + m[10] * m[3] * m[16]) - m[10] * m[4] * m[15]) - m[14] * m[3] * m[12]) + m[14] * m[4] * m[11]
    iA[6] = (((m[1] * m[11] * m[16] - m[1] * m[12] * m[15]) - m[9] * m[3] * m[16]) + m[9] * m[4] * m[15] + m[13] * m[3] * m[12]) - m[13] * m[4] * m[11]
    iA[10] = (((-(m[1]) * m[10] * m[16] + m[1] * m[12] * m[14] + m[9] * m[2] * m[16]) - m[9] * m[4] * m[14]) - m[13] * m[2] * m[12]) + m[13] * m[4] * m[10]
    iA[14] = (((m[1] * m[10] * m[15] - m[1] * m[11] * m[14]) - m[9] * m[2] * m[15]) + m[9] * m[3] * m[14] + m[13] * m[2] * m[11]) - m[13] * m[3] * m[10]
    iA[3] = (((m[2] * m[7] * m[16] - m[2] * m[8] * m[15]) - m[6] * m[3] * m[16]) + m[6] * m[4] * m[15] + m[14] * m[3] * m[8]) - m[14] * m[4] * m[7]
    iA[7] = (((-(m[1]) * m[7] * m[16] + m[1] * m[8] * m[15] + m[5] * m[3] * m[16]) - m[5] * m[4] * m[15]) - m[13] * m[3] * m[8]) + m[13] * m[4] * m[7]
    iA[11] = (((m[1] * m[6] * m[16] - m[1] * m[8] * m[14]) - m[5] * m[2] * m[16]) + m[5] * m[4] * m[14] + m[13] * m[2] * m[8]) - m[13] * m[4] * m[6]
    iA[15] = (((-(m[1]) * m[6] * m[15] + m[1] * m[7] * m[14] + m[5] * m[2] * m[15]) - m[5] * m[3] * m[14]) - m[13] * m[2] * m[7]) + m[13] * m[3] * m[6]
    iA[4] = (((-(m[2]) * m[7] * m[12] + m[2] * m[8] * m[11] + m[6] * m[3] * m[12]) - m[6] * m[4] * m[11]) - m[10] * m[3] * m[8]) + m[10] * m[4] * m[7]
    iA[8] = (((m[1] * m[7] * m[12] - m[1] * m[8] * m[11]) - m[5] * m[3] * m[12]) + m[5] * m[4] * m[11] + m[9] * m[3] * m[8]) - m[9] * m[4] * m[7]
    iA[12] = (((-(m[1]) * m[6] * m[12] + m[1] * m[8] * m[10] + m[5] * m[2] * m[12]) - m[5] * m[4] * m[10]) - m[9] * m[2] * m[8]) + m[9] * m[4] * m[6]
    iA[16] = (((m[1] * m[6] * m[11] - m[1] * m[7] * m[10]) - m[5] * m[2] * m[11]) + m[5] * m[3] * m[10] + m[9] * m[2] * m[7]) - m[9] * m[3] * m[6]

    d = m[1] * iA[1] + m[2] * iA[5] + m[3] * iA[9] + m[4] * iA[13]
    d==0 && error("Matrix is singular")

    d = 1.0 / d
    for i = 1:16
        iA[i] *= d
    end
    return iA
end


text:
using LichessBot
using Test

@time begin
    include("eval.jl")
    include("netcode.jl")
    # Write your tests here.
end


text:
module Diagnostics

export NaNChecker, CFL, AdvectiveCFL, DiffusiveCFL, WindowedSpatialAverage

using CUDA
using Oceananigans
using Oceananigans.Operators

using Oceananigans: AbstractDiagnostic
using Oceananigans.Utils: TimeInterval, IterationInterval, WallTimeInterval

import Oceananigans: run_diagnostic!

include("nan_checker.jl")
include("cfl.jl")
include("field_slicer.jl")
include("windowed_spatial_average.jl")

end