text:
using PencilFFTs
using .Transforms: binv, is_inplace

import FFTW
using MPI

using LinearAlgebra
using Random
using Test
using TimerOutputs

const DATA_DIMS = (16, 12, 6)

const FAST_TESTS = !("--all" in ARGS)

# Test all possible r2r transforms.
const TEST_KINDS_R2R = Transforms.R2R_SUPPORTED_KINDS

function test_transform_types(size_in)
    @testset "r2c transforms" begin
        transforms = (Transforms.RFFT(), Transforms.FFT(), Transforms.FFT())
        fft_params = PencilFFTs.GlobalFFTParams(size_in, transforms)

        @test fft_params isa PencilFFTs.GlobalFFTParams{Float64, 3, false,
                                                        typeof(transforms)}
        @test binv(Transforms.RFFT(), 10) === Transforms.BRFFT(10)
        @test binv(Transforms.RFFT(), 11) === Transforms.BRFFT(11)

        transforms_binv = binv.(transforms, size_in)
        size_out = Transforms.length_output.(transforms, size_in)

        @test transforms_binv ===
            (Transforms.BRFFT(size_in[1]), Transforms.BFFT(), Transforms.BFFT())
        @test size_out === (size_in[1] ÷ 2 + 1, size_in[2:end]...)
        @test Transforms.length_output.(transforms_binv, size_out) === size_in

        @test PencilFFTs.input_data_type(fft_params) === Float64
    end

    @testset "c2c transforms" begin
        n = 4
        @test binv(Transforms.FFT(), n) === Transforms.BFFT()
        @test binv(Transforms.BFFT(), n) === Transforms.FFT()
        @test binv(Transforms.FFT!(), n) === Transforms.BFFT!()
        @test binv(Transforms.BFFT!(), n) === Transforms.FFT!()
    end

    @testset "NoTransform" begin
        transform = Transforms.NoTransform()
        transform! = Transforms.NoTransform!()

        @test binv(transform, 42) === transform
        @test binv(transform!, 42) === transform!

        @test !is_inplace(transform)
        @test is_inplace(transform!)

        x = rand(4)
        p = Transforms.plan(transform, x)
        p! = Transforms.plan(transform!, x)

        @test p * x !== x  # creates copy
        @test p * x == x
        @test p \ x !== x  # creates copy
        @test p \ x == x

        @test p! * x === x
        @test p! \ x === x

        y = similar(x)
        @test mul!(y, p, x) === y == x
        @test mul!(x, p, x) === x  # this is also allowed

        rand!(x)
        @test ldiv!(x, p, y) === x == y

        # in-place IdentityPlan applied to out-of-place data
        @test_throws ArgumentError mul!(y, p!, x)
        @test mul!(x, p!, x) === x
        @test ldiv!(x, p!, x) === x
        @test mul!(x, p, x) === x
        @test ldiv!(x, p, x) === x
    end

    # Test type stability of generated plan_r2r (which, as defined in FFTW.jl,
    # is type unstable!). See comments of `plan` in src/Transforms/r2r.jl.
    @testset "r2r transforms" begin
        kind = FFTW.REDFT01
        transform = Transforms.R2R(kind)
        transform! = Transforms.R2R!(kind)

        @inferred (() -> Transforms.R2R(FFTW.REDFT10))()
        @inferred (() -> Transforms.R2R!(FFTW.REDFT10))()

        let kind_inv = FFTW.REDFT10
            @test binv(transform, 42) === Transforms.R2R(kind_inv)
            @test binv(transform!, 42) === Transforms.R2R!(kind_inv)
        end

        A = zeros(4, 6, 8)
        for tr in (transform, transform!)
            @inferred Transforms.plan(tr, A, 2)
            @inferred Transforms.plan(tr, A, (1, 3))
            @inferred Transforms.plan(tr, A)

            # This will fail because length(2:3) is not known by the compiler.
            @test_throws ErrorException @inferred Transforms.plan(tr, A, 2:3)
        end

        for kind in (FFTW.R2HC, FFTW.HC2R), T in (Transforms.R2R, Transforms.R2R!)
            # Unsupported r2r kinds.
            @test_throws ArgumentError T(kind)
        end
    end

    @testset "In-place transforms 1D" begin
        FFT = Transforms.FFT()
        FFT! = Transforms.FFT!()

        @inferred Transforms.is_inplace(FFT, FFT, FFT!)
        @inferred Transforms.is_inplace(FFT!, FFT, FFT)
        @inferred Transforms.is_inplace(FFT, FFT, FFT)
        @inferred Transforms.is_inplace(FFT!, FFT!, FFT!)

        @test Transforms.is_inplace(FFT, FFT, FFT!) === nothing
        @test Transforms.is_inplace(FFT, FFT!, FFT) === nothing
        @test Transforms.is_inplace(FFT, FFT!, FFT!) === nothing
        @test Transforms.is_inplace(FFT!, FFT, FFT!) === nothing
        @test Transforms.is_inplace(FFT!, FFT!, FFT!) === true
        @test Transforms.is_inplace(FFT, FFT, FFT) === false

        @inferred PencilFFTs.GlobalFFTParams(size_in, (FFT!, FFT!, FFT!))

        # Cannot combine in-place and out-of-place transforms.
        @test_throws ArgumentError PencilFFTs.GlobalFFTParams(size_in, (FFT, FFT!, FFT!))
    end

    @testset "Transforms internals" begin
        FFT = Transforms.FFT()
        x = zeros(ComplexF32, 3, 4)
        @test Transforms.scale_factor(FFT, x) == length(x)
    end

    nothing
end

function test_inplace_fft(
        ::Type{T}, comm, proc_dims, size_in;
        extra_dims=(),
    ) where {T}
    transforms = Transforms.FFT!()  # in-place c2c FFT
    plan = PencilFFTPlan(size_in, transforms, proc_dims, comm, T;
                         extra_dims=extra_dims)

    # Out-of-place plan, just for verifying that we throw errors.
    plan_oop = PencilFFTPlan(size_in, Transforms.FFT(), proc_dims, comm, T;
                             extra_dims=extra_dims)

    dims_fft = 1:length(size_in)

    @testset "In-place transforms 3D" begin
        test_transform(plan, x -> FFTW.plan_fft!(x, dims_fft), plan_oop)
    end

    nothing
end

function test_transform(plan::PencilFFTPlan, args...; kwargs...)
    println("\n", "-"^60, "\n\n", plan, "\n")

    @inferred allocate_input(plan)
    @inferred allocate_input(plan, 2, 3)
    @inferred allocate_input(plan, Val(3))
    @inferred allocate_output(plan)

    test_transform(Val(is_inplace(plan)), plan, args...; kwargs...)
end

function test_transform(inplace::Val{true}, plan::PencilFFTPlan,
                        serial_planner::Function, plan_oop::PencilFFTPlan;
                        root=0)
    @assert !is_inplace(plan_oop)

    vi = allocate_input(plan)
    @test vi isa PencilArrays.ManyPencilArray

    let vo = allocate_output(plan)
        @test typeof(vi) === typeof(vo)
    end

    u = first(vi)  # input PencilArray
    v = last(vi)   # output PencilArray

    randn!(u)
    u_initial = copy(u)
    ug = gather(u, root)  # for comparison with serial FFT

    # Input array type ManyPencilArray{...} incompatible with out-of-place
    # plans.
    @test_throws ArgumentError plan_oop * vi

    # Out-of-place plan applied to in-place data.
    @test_throws ArgumentError mul!(v, plan_oop, u)

    let vi_other = allocate_input(plan)
        # Input and output arrays for in-place plan must be the same.
        @test_throws ArgumentError mul!(vi_other, plan, vi)
    end

    # Input array type incompatible with in-place plans.
    @test_throws ArgumentError plan * u
    @test_throws ArgumentError plan \ v

    @assert PencilFFTs.is_inplace(plan)
    plan * vi               # apply in-place forward transform
    @test isempty(u) || !(u ≈ u_initial)  # `u` was modified!

    # Now `v` contains the transformed data.
    vg = gather(v, root)
    if ug !== nothing && vg !== nothing
        p = serial_planner(ug)
        p * ug  # apply serial in-place FFT
        @test ug ≈ vg
    end

    plan \ vi  # apply in-place backward transform

    # Now `u` contains the initial data (approximately).
    @test u ≈ u_initial

    ug_again = gather(u, root)
    if ug !== nothing && ug_again !== nothing
        p \ ug  # apply serial in-place FFT
        @test ug ≈ ug_again
    end

    let components = ((Val(3), ), (3, 2))
        @testset "Components: $comp" for comp in components
            vi = allocate_input(plan, comp...)
            u = first.(vi)
            v = last.(vi)
            randn!.(u)
            u_initial = copy.(u)

            # In some cases, generally when data is split among too many
            # processes, the local process may have no data.
            empty = isempty(first(u))

            plan * vi
            @test empty || !all(u_initial .≈ u)

            plan \ vi
            @test all(u_initial .≈ u)
        end
    end

    nothing
end

function test_transform(
        inplace::Val{false}, plan::PencilFFTPlan, serial_planner;
        root = 0, is_c2r = false,
    )
    u = allocate_input(plan)
    v = allocate_output(plan)

    @test u isa PencilArray

    if is_c2r
        # Generate input via inverse r2c transform, to avoid issues with input
        # not respecting required symmetry.
        randn!(v)
        ldiv!(u, plan, v)
        fill!(v, 0)
    else
        randn!(u)
    end

    mul!(v, plan, u)
    uprime = similar(u)
    ldiv!(uprime, plan, v)

    @test u ≈ uprime

    # Compare result with serial FFT.
    ug = gather(u, root)
    vg = gather(v, root)

    if ug !== nothing && vg !== nothing && serial_planner !== nothing
        let p = serial_planner(ug)
            vg_serial = p * ug
            if is_c2r
                @test size(vg) != size(vg_serial)
            else
                @test vg ≈ vg_serial
            end
        end
    end

    nothing
end

function test_transforms(::Type{T}, comm, proc_dims, size_in;
                         extra_dims=()) where {T}
    plan_kw = (:extra_dims => extra_dims, )
    N = length(size_in)

    make_plan(planner, args...; dims=1:N) = x -> planner(x, args..., dims)

    pair_r2r(tr::Transforms.R2R) =
        tr => make_plan(FFTW.plan_r2r, Transforms.kind(tr))
    pairs_r2r = (pair_r2r(Transforms.R2R(k)) for k in TEST_KINDS_R2R)

    pairs = if FAST_TESTS &&
            (T === Float32 || !isempty(extra_dims) || length(proc_dims) == 1)
        # Only test one transform with Float32 / extra_dims / 1D decomposition.
        (Transforms.RFFT() => make_plan(FFTW.plan_rfft), )
    else
        (
         Transforms.FFT() => make_plan(FFTW.plan_fft),
         Transforms.RFFT() => make_plan(FFTW.plan_rfft),
         Transforms.BFFT() => make_plan(FFTW.plan_bfft),
         Transforms.NoTransform() => (x -> Transforms.IdentityPlan()),
         pairs_r2r...,
         (Transforms.NoTransform(), Transforms.RFFT(), Transforms.FFT())
             => make_plan(FFTW.plan_rfft, dims=2:3),
         (Transforms.FFT(), Transforms.NoTransform(), Transforms.FFT())
             => make_plan(FFTW.plan_fft, dims=(1, 3)),
         (Transforms.FFT(), Transforms.NoTransform(), Transforms.NoTransform())
             => make_plan(FFTW.plan_fft, dims=1),

         # TODO compare BRFFT with serial equivalent?
         # The special case of BRFFT is a bit complicated, because
         # multidimensional FFTW plans returned by `plan_brfft` perform the
         # actual c2r transform along the first dimension. In PencilFFTs we do
         # the opposite: the c2r transform is applied along the last dimension.
         Transforms.BRFFT(size_in) => nothing,
        )
    end

    @testset "$(p.first) -- $T" for p in pairs
        transform, fftw_planner = p
        is_c2r = transform isa Transforms.BRFFT
        plan = @inferred PencilFFTPlan(
            size_in, transform, proc_dims, comm, T; plan_kw...)
        test_transform(plan, fftw_planner; is_c2r)
    end

    nothing
end

function test_pencil_plans(size_in::Tuple, pdims::Tuple, comm)
    N = length(size_in)
    @assert N >= 3

    @inferred PencilFFTPlan(size_in, Transforms.RFFT(), pdims, comm, Float64)

    let to = TimerOutput()
        plan = PencilFFTPlan(size_in, Transforms.RFFT(), pdims, comm, Float64,
                             timer=to)
        @test timer(plan) === to
    end

    @testset "Transform types" begin
        let transforms = (Transforms.RFFT(), Transforms.FFT(), Transforms.FFT())
            @inferred PencilFFTPlan(size_in, transforms, pdims, comm)
            @inferred PencilFFTs._input_data_type(Float64, transforms...)
        end

        let transforms = (Transforms.NoTransform(), Transforms.FFT())
            @test PencilFFTs._input_data_type(Float32, transforms...) ===
                ComplexF32
            @inferred PencilFFTs._input_data_type(Float32, transforms...)
        end

        let transforms = (Transforms.NoTransform(), Transforms.NoTransform())
            @test PencilFFTs._input_data_type(Float32, transforms...) ===
                Float32
            @inferred PencilFFTs._input_data_type(Float32, transforms...)
        end
    end

    if FAST_TESTS && length(pdims) == 1
        # Only test one case for 1D decomposition.
        types = (Float64, )
        extra_dims = ((), )
    else
        types = (Float64, Float32)
        extra_dims = ((), (3, ))
    end

    for T in types, edims in extra_dims
        test_inplace_fft(T, comm, pdims, size_in, extra_dims=edims)
        test_transforms(T, comm, pdims, size_in, extra_dims=edims)
    end

    @testset "FFT! + NoTransform!" begin
        transforms = (
            ntuple(_ -> Transforms.FFT!(), N - 1)...,
            Transforms.NoTransform!(),
        )
        transforms_oop = (
            ntuple(_ -> Transforms.FFT(), N - 1)...,
            Transforms.NoTransform(),
        )
        plan = PencilFFTPlan(size_in, transforms, pdims, comm)
        plan_oop = PencilFFTPlan(size_in, transforms_oop, pdims, comm)
        fftw_planner(x) = FFTW.plan_fft!(x, 1:(N - 1))
        test_transform(plan, fftw_planner, plan_oop)
    end

    nothing
end

# Test N-dimensional transforms decomposing along M dimensions.
function test_dimensionality(dims::Dims{N}, ::Val{M}, comm;
                             plan_kw...) where {N, M}
    @assert M < N
    pdims = make_pdims(Val(M), MPI.Comm_size(comm))

    @testset "Decompose $M/$N dims" begin
        # Out-of-place transform.
        let transform = Transforms.RFFT()
            plan = PencilFFTPlan(dims, transform, pdims, comm; plan_kw...)
            test_transform(plan, FFTW.plan_rfft)
        end

        # In-place transform.
        let transform = Transforms.FFT!()
            plan = PencilFFTPlan(dims, transform, pdims, comm; plan_kw...)
            plan_oop = PencilFFTPlan(dims, Transforms.FFT(), pdims, comm;
                                     plan_kw...)
            test_transform(plan, FFTW.plan_fft!, plan_oop)
        end

    end

    nothing
end

function test_dimensionality(comm)
    # 1D decomposition of 2D problem.
    test_dimensionality((11, 15), Val(1), comm)

    let dims = (11, 7, 13)
        test_dimensionality(dims, Val(1), comm)  # slab decomposition
        test_dimensionality(dims, Val(2), comm)  # pencil decomposition
    end

    let dims = (9, 7, 5, 12)
        test_dimensionality(dims, Val(1), comm)
        test_dimensionality(dims, Val(2), comm)
        test_dimensionality(dims, Val(3), comm)  # 3D decomposition of 4D problem

        # Same with some non-default options for the plans.
        test_dimensionality(
            dims, Val(3), comm,
            permute_dims=Val(false),
            transpose_method=Transpositions.Alltoallv(),
        )
    end

    nothing
end

# Test incompatibilities between plans and inputs.
function test_incompatibility(comm)
    pdims = (MPI.Comm_size(comm), )
    dims = (10, 8)
    dims_other = (6, 8)

    @testset "Incompatibility" begin
        plan = PencilFFTPlan(dims, Transforms.FFT(), pdims, comm)
        plan! = PencilFFTPlan(dims, Transforms.FFT!(), pdims, comm)

        u = allocate_input(plan)  :: PencilArray
        v = allocate_output(plan) :: PencilArray

        # "input array type PencilArray{...} incompatible with in-place plans"
        @test_throws ArgumentError plan! * u

        # "input array type ManyPencilArray{...} incompatible with out-of-place plans"
        M! = allocate_input(plan!) :: ManyPencilArray
        @test_throws ArgumentError plan * M!

        # "array types (...) incompatible with in-place plans"
        @test_throws ArgumentError mul!(v, plan!, u)

        # "array types (...) incompatible with out-of-place plan"
        @test_throws ArgumentError mul!(M!, plan, M!)

        # "out-of-place plan applied to aliased data"
        @test_throws ArgumentError mul!(last(M!), plan, first(M!))

        # "collections have different lengths: 3 ≠ 2"
        u3 = allocate_input(plan, Val(3))
        v2 = allocate_output(plan, Val(2))
        @test_throws ArgumentError mul!(v2, plan, u3)

        let plan_other = PencilFFTPlan(dims_other, Transforms.FFT(), pdims, comm)
            # "unexpected dimensions of input data"
            @test_throws ArgumentError plan_other * u

            # "unexpected dimensions of output data"
            v_other = allocate_output(plan_other)
            @test_throws ArgumentError mul!(v_other, plan, u)
        end
    end

    nothing
end

function make_pdims(::Val{M}, Nproc) where {M}
    # Let MPI.Dims_create! choose the decomposition.
    pdims = zeros(Int, M)
    MPI.Dims_create!(Nproc, pdims)
    ntuple(d -> pdims[d], Val(M))
end

MPI.Init()

size_in = DATA_DIMS
comm = MPI.COMM_WORLD
Nproc = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
rank == 0 || redirect_stdout(devnull)

test_transform_types(size_in)
test_incompatibility(comm)
test_dimensionality(comm)

pdims_1d = (Nproc, )  # 1D ("slab") decomposition
pdims_2d = make_pdims(Val(2), Nproc)

for p in (pdims_1d, pdims_2d)
    test_pencil_plans(size_in, p, comm)
end


text:
Random.seed!(0)

n = 4
diags = Subspace([ Array{ComplexF64, 2}(basis_vec((n,n), (i,i))) for i in 1:n ])
S = S0Graph([1 n], diags)

w = random_bounded(n)
λ = dsw(S, w).λ

X = HermitianSemidefinite(n)
problem = maximize(real(tr(X * w)), [ X[i,i] == 1 for i in 1:n ])
solve!(problem, () -> SCS.Optimizer(verbose=0, eps=solver_eps))
@test λ ≈ problem.optval  atol=tol

ext: jl

max_stars_repo_name: dstahlke/NoncommutativeGraphs

max_issues_repo_name: dstahlke/NoncommutativeGraphs

max_forks_repo_name: dstahlke/NoncommutativeGraphs


text:
# glo_functions.jl
using Dates
using Unitful
using UnitfulAstro
import GPLinearODEMaker

"""
A structure that holds all of the relevant information for constructing the
model used in the Jones et al. 2017+ paper (https://arxiv.org/pdf/1711.01318.pdf).
"""
struct GLO_RV
	GLO::GPLinearODEMaker.GLO  # kernel function
	time::Vector{T} where T<:Unitful.Time
	rv_factor::Unitful.Velocity

	function GLO_RV(
		glo::GPLinearODEMaker.GLO,
		time_unit::Unitful.Time,
		rv_factor::Unitful.Velocity)
		return new(glo, glo.x_obs .* time_unit, rv_factor)
	end
end

get_rv(glo::GLO_RV) = glo.GLO.y_obs[1:glo.GLO.n_out:end] .* glo.rv_factor
get_rv_noise(glo::GLO_RV) = glo.GLO.noise[1:glo.GLO.n_out:end] .* glo.rv_factor


text:
"""
`SpectralBVPProblem` - Initalize the data structure
"""
function SpectralBVPProblem(a, b, n, f, bcl::TBCL, bcr::TBCR) where {TBCL<:PeriodicBC, TBCR<:PeriodicBC}
    Δx = (b-a)/(n+1);
    x = LinRange(a, b, n+2)[1:end-1];
    N = n+1;

    # allocate the matrix
    λ = zeros(n+1);
    rhs = zeros(n+1);

    return SpectralBVPProblem(a, b, Δx, n, N, x, f, λ, rhs, bcl, bcr)
end

function SpectralBVPProblem(a, b, n, f, bcl::TBCL, bcr::TBCR) where {TBCL<:DirichletBC, TBCR<:DirichletBC}
    Δx = (b-a)/(n+1);
    x = LinRange(a, b, n+2)[2:end-1];
    N = 2*(n+1);

    # allocate the matrix
    λ = zeros(N);
    rhs = zeros(N);

    return SpectralBVPProblem(a, b, Δx, n, N, x, f, λ, rhs, bcl, bcr)
end

function SpectralBVPProblem(a, b, n, f, bcl::TBCL, bcr::TBCR) where {TBCL<:NeumannBC, TBCR<:NeumannBC}
    Δx = (b-a)/(n+1);
    x = LinRange(a, b, n+2);
    N = 2*(n+1);

    # allocate the matrix
    λ = zeros(N);
    rhs = zeros(N);

    return SpectralBVPProblem(a, b, Δx, n, N, x, f, λ, rhs, bcl, bcr)
end


text:
type TcpListener
    ptr::Ptr{Void}

    function TcpListener(ptr::Ptr{Void})
        t = new(ptr)
        finalizer(t, destroy)
        t
    end
end

function TcpListener()
    TcpListener(ccall((:sfTcpListener_create, libcsfml_network), Ptr{Void}, ()))
end

function destroy(listener::TcpListener)
    ccall((:sfTcpListener_destroy, libcsfml_network), Void, (Ptr{Void},), listener.ptr)
end

function set_blocking(listener::TcpListener, blocking::Bool)
    ccall((:sfTcpListener_setBlocking, libcsfml_network), Void, (Ptr{Void}, Int32,), listener.ptr, blocking)
end

function is_blocking(listener::TcpListener)
    return ccall((:sfTcpListener_isBlocking, libcsfml_network), Int32, (Ptr{Void},), listener.ptr) == 1
end

function get_localport(listener::TcpListener)
    return ccall((:sfTcpListener_getLocalPort, libcsfml_network), UInt16, (Ptr{Void},), listener.ptr)
end

function listen(listener::TcpListener, port::Integer)
    SocketStatus(ccall((:sfTcpListener_listen, libcsfml_network), Int32, (Ptr{Void}, UInt16,), listener.ptr, port))
end

function accept(listener::TcpListener, socket::SocketTCP)
    nstruct = ccall((:sjTcpListener_accept, "libjuliasfml"), NetworkStruct, (Ptr{Void}, Ptr{Void},), listener.ptr, socket.ptr)
    socket.ptr = nstruct.ptr
    return SocketStatus(nstruct.status)
end


text:
# BAM Auxiliary Data
# ==================

struct AuxData <: AbstractDict{String,Any}
    data::Vector{UInt8}
end

function Base.getindex(aux::AuxData, tag::AbstractString)
    checkauxtag(tag)
    return getauxvalue(aux.data, 1, length(aux.data), UInt8(tag[1]), UInt8(tag[2]))
end

function Base.length(aux::AuxData)
    data = aux.data
    p = 1
    len = 0
    while p ≤ length(data)
        len += 1
        p = next_tag_position(data, p)
    end
    return len
end

function Base.iterate(aux::AuxData, pos=1)
    if pos > length(aux.data)
        return nothing
    end

    data = aux.data
    @label doit
    t1 = data[pos]
    t2 = data[pos+1]
    pos, typ = loadauxtype(data, pos + 2)
    pos, value = loadauxvalue(data, pos, typ)
    if t1 == t2 == 0xff
        @goto doit
    end
    return Pair{String,Any}(String([t1, t2]), value), pos
end

# Internals
# ---------

function checkauxtag(tag::AbstractString)
    if sizeof(tag) != 2
        throw(ArgumentError("tag length must be 2"))
    end
end

function getauxvalue(data::Vector{UInt8}, start::Int, stop::Int, t1::UInt8, t2::UInt8)
    pos = findauxtag(data, start, stop, t1, t2)
    if pos == 0
        throw(KeyError(String([t1, t2])))
    end
    pos, T = loadauxtype(data, pos + 2)
    _, val = loadauxvalue(data, pos, T)
    return val
end

function loadauxtype(data::Vector{UInt8}, p::Int)
    function auxtype(b)
        return (
            b == UInt8('A') ? Char  :
            b == UInt8('c') ? Int8  :
            b == UInt8('C') ? UInt8 :
            b == UInt8('s') ? Int16 :
            b == UInt8('S') ? UInt16 :
            b == UInt8('i') ? Int32 :
            b == UInt8('I') ? UInt32 :
            b == UInt8('f') ? Float32 :
            b == UInt8('Z') ? String :
            error("invalid type tag: '$(Char(b))'"))
    end
    t = data[p]
    if t == UInt8('B')
        return p + 2, Vector{auxtype(data[p+1])}
    else
        return p + 1, auxtype(t)
    end
end

function loadauxvalue(data::Vector{UInt8}, p::Int, ::Type{T}) where T
    return p + sizeof(T), unsafe_load(Ptr{T}(pointer(data, p)))
end

function loadauxvalue(data::Vector{UInt8}, p::Int, ::Type{Char})
    return p + 1, Char(unsafe_load(pointer(data, p)))
end

function loadauxvalue(data::Vector{UInt8}, p::Int, ::Type{Vector{T}}) where T
    n = unsafe_load(Ptr{Int32}(pointer(data, p)))
    p += 4
    xs = Vector{T}(undef, n)
    unsafe_copyto!(pointer(xs), Ptr{T}(pointer(data, p)), n)
    return p + n * sizeof(T), xs
end

function loadauxvalue(data::Vector{UInt8}, p::Int, ::Type{String})
    dataptr = pointer(data, p)
    endptr = ccall(:memchr, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t),
                   dataptr, '\0', length(data) - p + 1)
    q::Int = p + (endptr - dataptr) - 1
    return q + 2, String(data[p:q])
end

function findauxtag(data::Vector{UInt8}, start::Int, stop::Int, t1::UInt8, t2::UInt8)
    pos = start
    while pos ≤ stop && !(data[pos] == t1 && data[pos+1] == t2)
        pos = next_tag_position(data, pos)
    end
    if pos > stop
        return 0
    else
        return pos
    end
end

# Find the starting position of a next tag in `data` after `p`.
# `(data[p], data[p+1])` is supposed to be a current tag.
function next_tag_position(data::Vector{UInt8}, p::Int)
    typ = Char(data[p+2])
    p += 3
    if typ == 'A'
        p += 1
    elseif typ == 'c' || typ == 'C'
        p += 1
    elseif typ == 's' || typ == 'S'
        p += 2
    elseif typ == 'i' || typ == 'I'
        p += 4
    elseif typ == 'f'
        p += 4
    elseif typ == 'd'
        p += 8
    elseif typ == 'Z' || typ == 'H'
        while data[p] != 0x00  # NULL-terminalted string
            p += 1
        end
        p += 1
    elseif typ == 'B'
        eltyp = Char(data[p])
        elsize = eltyp == 'c' || eltyp == 'C'                 ? 1 :
                 eltyp == 's' || eltyp == 'S'                 ? 2 :
                 eltyp == 'i' || eltye == 'I' || eltyp == 'f' ? 4 :
                 error("invalid type tag: '$(Char(eltyp))'")
        p += 1
        n = unsafe_load(Ptr{Int32}(pointer(data, p)))
        p += 4 + elsize * n
    else
        error("invalid type tag: '$(Char(typ))'")
    end
    return p
end


text:
struct Benchmark{D, T, K1, K2}
    filename::String
    name::String
    complex::UnionAll
    data::D
    threshold::T
    kwargs::K1
    extra_kwargs::K2
end

Base.show(io::IO, b::Benchmark) = print(io, "Benchmark($(b.name))")

function Benchmark(filename; threshold=nothing, extra::K=(;), kwargs...) where K
    name, ext = splitext(filename)
    filename = joinpath(@__DIR__, "..", "datasets", filename)
    data = load_data(filename)
    complex = if ext == ".dipha"
        Cubical
    elseif ext == ".alpha"
        Alpha
    else
        Rips
    end
    T = typeof(threshold)
    return Benchmark{typeof(data), T, typeof(kwargs), K}(
        filename, name, complex, data, threshold, kwargs, extra
    )
end

Ripserer.nv(b::Benchmark) = Ripserer.nv(b.complex(b.data))
get_meta(arg) = function(b)
    return string(if arg == :threshold
        !isnothing(b.threshold) ? b.threshold : ""
    elseif arg == :size
        Ripserer.nv(b)
    elseif arg == :complex
        nameof(b.complex)
    elseif haskey(b.kwargs, arg)
        b.kwargs[arg]
    elseif haskey(b.extra_kwargs, arg)
        b.extra_kwargs[arg]
    else
        ""
    end)
end


text:
using Convex1d
using Test

@testset "minimize" begin
    res = @inferred minimize(x -> x^2, (-1.1, 1), atol=1e-3)
    @test res.minimizer ≈ 0 atol=1e-3
    @test res.minimum === res.minimizer^2

    x = [1,2]
    y = [2,4]
    f = c -> sum(abs, c*x - y)
    sol = minimize(f, (-10, 10), atol=1e-10)
    @test sol.minimizer ≈ 2 atol=1e-10

    for _ in 1:100
        x_opt = randn()
        f_opt = randn()
        s = 0.1 + 10rand()
        p = 1 + rand()
        f = x -> s*abs(x - x_opt)^p + f_opt
        @assert f(x_opt) ≈ f_opt
        atol = eps(x_opt)^(1/3)
        res = minimize(f, (x_opt - 10rand(), x_opt + 10rand()), atol=atol)
        @test res.nevals < 100 # quick convergence
        @test res.xbounds[1] <= res.minimizer <= res.xbounds[2]
        @test res.xbounds[2] - res.xbounds[1] <= atol
        # @test res.ybounds[1] <= res.minimum <= res.ybounds[2]
        # @test res.ybounds[1] <= f_opt <= res.ybounds[2]

        x_res, f_res = res
        @test x_opt ≈ x_res atol=atol
        @test f(x_res) == f_res

        res = minimize(f, (x_opt, x_opt + 10rand()), atol=atol)
        x_res, f_res = res
        @test x_opt ≈ x_res atol=atol
        @test f(x_res) == f_res

        I = (x_opt-10rand(), x_opt)
        res = minimize(f, I, atol=atol)
        x_res, f_res = res
        if !(isapprox(x_opt, x_res, atol=atol))
            @show x_opt
            @show f(x_opt)
            @show f_opt
            @show s
            @show p
            @show res
            @show I
            @show atol
            error()
        end
        @test x_opt ≈ x_res atol=atol
        @test f(x_res) == f_res
    end
end

@testset "minimize no domain" begin
    for scale in [10^x for x in -10.0:10.0]
        x_opt = scale*randn()
        f = x -> abs(x - x_opt)
        x_lo, x_hi = Convex1d.find_initial_domain(f)
        @test x_lo <= x_opt <= x_hi
        @test x_hi - x_lo < 8*(scale+1)
        sol = minimize(f, atol=sqrt(eps(scale)))
        @test sol.minimizer ≈ x_opt atol=sqrt(eps(scale))
    end
end


text:
using ClimateMachine.Grids

using KernelAbstractions
using KernelAbstractions.Extras: @unroll
using StaticArrays

# Split off into separate parameter struct
const _ξ1x1, _ξ2x1, _ξ3x1 = Grids._ξ1x1, Grids._ξ2x1, Grids._ξ3x1
const _ξ1x2, _ξ2x2, _ξ3x2 = Grids._ξ1x2, Grids._ξ2x2, Grids._ξ3x2
const _ξ1x3, _ξ2x3, _ξ3x3 = Grids._ξ1x3, Grids._ξ2x3, Grids._ξ3x3
const _M, _MI = Grids._M, Grids._MI
const _x1, _x2, _x3 = Grids._x1, Grids._x2, Grids._x3
const _JcV = Grids._JcV

const _n1, _n2, _n3 = Grids._n1, Grids._n2, Grids._n3
const _sM, _vMI = Grids._sM, Grids._vMI

function launch_volume_divergence!(flux_divergence, flux, grid::G; dependencies = nothing) where {G <: AbstractGrid}

    nrealelem = length(grid.interiorelems)
    device = array_device(flux)
    N = polynomialorders(grid)
    dim = length(grid.D)

    FT = eltype(flux)
    workgroup = (N[1]+1, N[2]+1) # same size as shared memory tiles
    ndrange = ( (N[1]+1)*nrealelem, (N[2]+1))
    comp_stream = dependencies
    flux_divergence.realdata[:] .= 0 # necesary for now
    # If the model direction is EveryDirection, we need to perform
    # both horizontal AND vertical kernel calls; otherwise, we only
    # call the kernel corresponding to the model direction `dg.diffusion_direction`
    comp_stream = volume_divergence_kernel!(device, workgroup)(
        flux_divergence.data,
        flux.data,
        Val(dim),
        Val(N[1]),
        grid.vgeo,
        grid.D[1],
        ndrange = ndrange,
        dependencies = comp_stream,
    )

    return comp_stream
end

@kernel function volume_divergence_kernel!(
    flux_divergence,
    flux,
    ::Val{dim},
    ::Val{polyorder},
    vgeo,
    D,
) where {dim, polyorder}
    # This allows for the variables to remain after @synchronize
    # for CPU computation
    @uniform begin
        N = polyorder
        FT = eltype(flux)
        Nq = N + 1
        Nqk = (dim == 2) ? 1 : Nq
        local_flux = MArray{Tuple{3, 1}, FT}(undef)
        local_flux_3 = MArray{Tuple{1, 1}, FT}(undef)
    end

    # Arrays for F, and the differentiation matrix D
    shared_flux = @localmem FT (2, Nq, Nq) # shared memory on the gpu
    s_D = @localmem FT (Nq, Nq)            # shared memory on the gpu (can remove)

    # Storage for tendency and mass inverse M⁻¹, **perhaps in @uniform block**
    local_tendency = @private FT (Nqk,)   # thread private memory
    local_MI = @private FT (Nqk,)         # thread private memory

    # Grab the index associated with the current element `e` and the
    # horizontal quadrature indices `i` (in the ξ1-direction),
    # `j` (in the ξ2-direction) [directions on the reference element].
    # Parallelize over elements, then over columns
    e = @index(Group, Linear)    # Group is the index of the work group
    i, j = @index(Local, NTuple) # i, j is the index within a work group

    @inbounds begin
        # load differentiation matrix into local memory
        s_D[i, j] = D[i, j]
        @unroll for k in 1:Nqk # if 2D Nqk = 1, if 3D Nq
            ijk = i + Nq * ((j - 1) + Nq * (k - 1)) # convert to linear indices
            # initialize local tendency
            local_tendency[k] = zero(FT)
            # read in mass matrix inverse for element `e` in each plane k
            local_MI[k] = vgeo[ijk, _MI, e]
        end
        # end of boiler plate loading

        @unroll for k in 1:Nqk # search
            @synchronize # perhaps moved up outside the loop
            ijk = i + Nq * ((j - 1) + Nq * (k - 1))
            M = vgeo[ijk, _M, e]
            # Extract Jacobian terms
            # ∂/∂x¹ = ∂ξ¹/ ∂x¹ ∂/∂ξ¹ + ∂ξ²/ ∂x¹ ∂/∂ξ² + ∂ξ³/ ∂x¹ ∂/∂ξ³
            # ∂/∂x² = ∂ξ¹/ ∂x² ∂/∂ξ¹ + ∂ξ²/ ∂x² ∂/∂ξ² + ∂ξ³/ ∂x² ∂/∂ξ³
            # ∂/∂x³ = ∂ξ¹/ ∂x³ ∂/∂ξ¹ + ∂ξ²/ ∂x³ ∂/∂ξ² + ∂ξ³/ ∂x³ ∂/∂ξ³
            # The integration by parts is then done with respect to the
            # (ξ¹, ξ², ξ³) coordinate system
            # (wich is why terms associated with the columns appear)
            ξ1x1 = vgeo[ijk, _ξ1x1, e]
            ξ1x2 = vgeo[ijk, _ξ1x2, e]
            ξ1x3 = vgeo[ijk, _ξ1x3, e]

            ξ2x1 = vgeo[ijk, _ξ2x1, e]
            ξ2x2 = vgeo[ijk, _ξ2x2, e]
            ξ2x3 = vgeo[ijk, _ξ2x3, e]

            ξ3x1 = vgeo[ijk, _ξ3x1, e]
            ξ3x2 = vgeo[ijk, _ξ3x2, e]
            ξ3x3 = vgeo[ijk, _ξ3x3, e]
            # ∂ᵗθ = ∇⋅(u⃗θ -κ∇θ), F1 = uθ , F2 = vθ, F3 = wθ, where u⃗ = (u,v,w)
            # Computes the local inviscid fluxes Fⁱⁿᵛ
            # Need to load in user passed in flux aka make flux into local flux
            fill!(local_flux, -zero(eltype(local_flux))) # can remove this

            # probably becomes flux[ijk, e, 1]
            shared_flux[1, i, j] = flux[ijk, e, 1] # local_flux[1]

            # probably becomes flux[ijk, e, 2]
            shared_flux[2, i, j] = flux[ijk, e, 2] # local_flux[2]

            # not taking derivative in third direction
            local_flux_3         = flux[ijk, e, 3] # local_flux[3]

            # Build "inside metrics" flux
            F1 = shared_flux[1, i, j]
            F2 = shared_flux[2, i, j]
            F3 = local_flux_3

            shared_flux[1, i, j] = M * (ξ1x1 * F1 + ξ1x2 * F2 + ξ1x3 * F3)
            shared_flux[2, i, j] = M * (ξ2x1 * F1 + ξ2x2 * F2 + ξ2x3 * F3)
            local_flux_3         = M * (ξ3x1 * F1 + ξ3x2 * F2 + ξ3x3 * F3)
            ## Everything above is point-wise multiplication
            # Fix Me

            @unroll for n in 1:Nqk
                MI = local_MI[n]
                local_tendency[n] += MI * s_D[k, n] * local_flux_3
            end

            @synchronize

            # Weak "inside metrics" derivative.
            # Computes the rest of the volume term: M⁻¹SᵀF (note the transpose)
            # note: (1) The indices are reversed from the usual matrix mult
            #       (2) The stiffness matrix comes in since shared flux has
            #           the mass matrix already multiplied through
            MI = local_MI[k]
            @unroll for n in 1:Nq
                # ξ1-grid lines
                local_tendency[k] += MI * s_D[n, i] * shared_flux[1, n, j]
                # ξ2-grid lines
                local_tendency[k] += MI * s_D[n, j] * shared_flux[2, i, n]
            end

            #FIXME: ijk does not persist after @synchronize and must be redefined here
            #       Maybe try @private and change the variable to a tuple
            #ijk = i + Nq * ((j - 1) + Nq * (k - 1))
            #flux_divergence[ijk, e, 1] += local_tendency[k]
        end

        @unroll for k in 1:Nqk # search
            ijk = i + Nq * ((j - 1) + Nq * (k - 1))
            flux_divergence[ijk, e, 1] += local_tendency[k]
        end
    end
end

"""
    launch_interface_gradients!(dg, state_prognostic, t; surface::Symbol, dependencies)

Launches horizontal and vertical kernels for computing the interface gradients.
The argument `surface` is either `:interior` or `:exterior`, which denotes whether
we are computing interface gradients on boundaries which are interior (exterior resp.)
to the _parallel_ boundary.
"""
function launch_interface_divergence!(
    flux_divergence,
    flux,
    grid::G;
    dependencies = nothing,
) where {G <: AbstractGrid}
    # MPI
    # @assert surface === :interior || surface === :exterior

    FT = eltype(flux)

    device = array_device(flux)
    dim = length(grid.D)
    N = polynomialorders(grid)

    # Assumes poly order is the same in every direction
    # TODO: copy from new kernel with variable polynomial order
    if dim == 1
        Np = (N[1] + 1)
        Nfp = 1
        nface = 2
    elseif dim == 2
        Np = (N[1] + 1) * (N[2] + 1)
        Nfp = (N[1] + 1)
        nface = 4
    elseif dim == 3
        Np = (N[1] + 1) * (N[2] + 1) * (N[3] + 1)
        Nfp = (N[1] + 1) * (N[2] + 1)
        nface = 6
    end
    flux_divergence.realdata[:] .= 0
    workgroup = Nfp
    elems = grid.interiorelems
    ndrange = Nfp * length(grid.interiorelems)

    comp_stream = dependencies

    #TODO: fix with appropriate function calls
    # If the model direction is EveryDirection, we need to perform
    # both horizontal AND vertical kernel calls; otherwise, we only
    # call the kernel corresponding to the model direction `dg.diffusion_direction`
    comp_stream = interface_divergence_kernel!(device, workgroup)(
        flux_divergence.data,
        flux.data,
        Val(dim),
        Val(N[1]),
        grid.vgeo,
        grid.sgeo,
        grid.vmap⁻,
        grid.vmap⁺,
        grid.elemtobndy,
        elems,
        nface,
        Np,
        ndrange = ndrange,
        dependencies = comp_stream,
    )
    return comp_stream
end

@doc """
    function interface_divergence_kernel!(
        flux_divergence,
        flux,
        ::Val{dim},
        ::Val{polyorder},
        direction, # we will probably want to get rid of this
        numerical_flux_first_order,  # this is metadata, we need to specify what numerical flux we are using
        tendency,
        vgeo,
        sgeo,
        vmap⁻,
        vmap⁺,
        elemtobndy,
        elems,
    )

Compute kernel for evaluating the interface tendencies for the
DG form:

∫ₑ ψ⋅ ∂q/∂t dV - ∫ₑ ∇ψ⋅F⃗ dV + ∮ₑ ψF⃗⋅n̂ dS,

or equivalently in matrix form:

dQ/dt = M⁻¹(MS + DᵀM F⃗ + ∑ᶠ LᵀMᶠF⃗ ).

This kernel computes the surface terms: M⁻¹ ∑ᶠ LᵀMᶠ(F)),
where M is the mass matrix, Mf is the face mass matrix, L is an interpolator
from volume to face, and Fⁱⁿᵛ⋆, Fᵛⁱˢᶜ⋆
are the numerical fluxes for the inviscid and viscous
fluxes, respectively.  interface_tendency!
"""
@kernel function interface_divergence_kernel!(
    flux_divergence,
    flux,
    ::Val{dim},
    ::Val{polyorder},
    vgeo,
    sgeo,
    vmap⁻,
    vmap⁺,
    elemtobndy,
    elems,
    nface,
    Np,
) where {dim, polyorder}
    @uniform begin
        N = polyorder
        FT = eltype(flux)

        # faces: (1, 2, 3, 4, 5, 6)  = (West, East, South, North, Bottom, Top)
        # We think...
        faces = 1:nface

        Nq = N + 1
        Nqk = dim == 2 ? 1 : Nq
    end
    eI = @index(Group, Linear)
    n = @index(Local, Linear)

    e = @private Int (1,)

    e[1] = elems[eI]
    # goes over too many points right now
    for f in faces
        # Inter-element index
        e⁻ = e[1]
        normal_vector = SVector(
            sgeo[_n1, n, f, e⁻],
            sgeo[_n2, n, f, e⁻],
            sgeo[_n3, n, f, e⁻],
        )
        # Get surface mass, volume mass inverse

        # sM = surface mass matrix
        # vMI = volume mass inverse matrix
        sM, vMI = sgeo[_sM, n, f, e⁻], sgeo[_vMI, n, f, e⁻]

        # Ids corresponding to the workitem on the face
        # vmap⁻ and vmap⁺ are the linear indices
        id⁻, id⁺ = vmap⁻[n, f, e⁻], vmap⁺[n, f, e⁻]

        # neighbor element index
        e⁺ = ((id⁺ - 1) ÷ Np) + 1 # this is why 1 indexing is failure

        # vid's are cartesian indices of the GL points
        # needs to be modified by the faces?
        vid⁻, vid⁺ = ((id⁻ - 1) % Np) + 1, ((id⁺ - 1) % Np) + 1

        # No Flux BC and Central Fluxes for the interfaces allows the kernel
        # to be extensible
        # TODO: make it no flux on the boundary for now and add the flux later
        # bctype = elemtobndy[f, e⁻]

        # flux_divergence[vmap⁻] = (flux[vmap⁻] + flux[vmap⁺])/2
        # ∫dA F⃗⋅n̂
        # normal_vector[1] * (flux[vmap⁻,1] + flux[vmap⁺,1])/2
        # normal_vector[2] * (flux[vmap⁻,2] + flux[vmap⁺,2])/2
        # normal_vector[3] * (flux[vmap⁻,3] + flux[vmap⁺,3])/2
        # # ∂ᵗ∫dVρ = - ∫dS Φ⃗⋅n̂ = - flux on boundary

        local_flux  = normal_vector[1] * (flux[vid⁻, e⁻, 1] + flux[vid⁺, e⁺, 1])/2
        local_flux += normal_vector[2] * (flux[vid⁻, e⁻, 2] + flux[vid⁺, e⁺, 2])/2
        local_flux += normal_vector[3] * (flux[vid⁻, e⁻, 3] + flux[vid⁺, e⁺, 3])/2

        # Update RHS (in outer face loop): M⁻¹ Mfᵀ(Fⁱⁿᵛ⋆ ))
        flux_divergence[vid⁻, e⁻, 1] += vMI * sM * local_flux # Warning, need to think about this
        # @show vid⁻ e⁻ id⁻ local_flux flux_divergence[vid⁻, e⁻, 1] vMI sM f
        # Need to wait after even faces to avoid race conditions
        # Note: only need synchronization for F%2 == 0, but this crashed...
        @synchronize
    end
end

function launch_volume_gradient!(gradient, state, grid::G; dependencies = nothing) where {G <: AbstractGrid}

    nrealelem = length(grid.interiorelems)
    device = array_device(gradient)
    N = polynomialorders(grid)
    dim = length(grid.D)

    FT = eltype(state)
    workgroup = (N[1]+1, N[2]+1) # same size as shared memory tiles
    ndrange = ( (N[1]+1)*nrealelem, (N[2]+1))
    comp_stream = dependencies
    gradient.realdata[:] .= 0 # necessary for now
    # If the model direction is EveryDirection, we need to perform
    # both horizontal AND vertical kernel calls; otherwise, we only
    # call the kernel corresponding to the model direction `dg.diffusion_direction`
    comp_stream = volume_gradient_kernel!(device, workgroup)(
        gradient.data,
        state.data,
        Val(dim),
        Val(N),
        grid.vgeo,
        grid.D[1],
        ndrange = ndrange,
        dependencies = comp_stream,
    )

    return comp_stream
end

# TODO: should be similar to vectorgradients kernel
@kernel function volume_gradient_kernel!(
    gradient,
    state,
    ::Val{dim},
    ::Val{polyorder},
    vgeo,
    D,
) where {dim, polyorder}
    # This allows for the variables to remain after @synchronize
    # for CPU computation
    @uniform begin
        N = polyorder
        FT = eltype(state)
        @inbounds Nq1 = N[1] + 1
        @inbounds Nq2 = N[2] + 1
        Nq3 = (dim == 2) ? 1 : N[3] + 1
        local_flux = MArray{Tuple{3, 1}, FT}(undef)
        local_flux_3 = MArray{Tuple{1, 1}, FT}(undef)
    end

    # Arrays for F, and the differentiation matrix D
    shared_flux = @localmem FT (2, Nq1, Nq2) # shared memory on the gpu

    # Storage for tendency and mass inverse M⁻¹, **perhaps in @uniform block**
    local_tendency = @private FT (Nq3, dim)   # thread private memory
    local_MI = @private FT (Nq3,)         # thread private memory

    # Grab the index associated with the current element `e` and the
    # horizontal quadrature indices `i` (in the ξ1-direction),
    # `j` (in the ξ2-direction) [directions on the reference element].
    # Parallelize over elements, then over columns
    e = @index(Group, Linear)    # Group is the index of the work group
    i, j = @index(Local, NTuple) # i, j is the index within a work group

    @inbounds begin
        # load differentiation matrix into local memory
        @unroll for k in 1:Nq3 # if 2D Nq3 = 1, if 3D Nq
            ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1)) # convert to linear indices
            # initialize local tendency
            local_tendency[k, 1] = zero(FT)
            local_tendency[k, 2] = zero(FT)
            local_tendency[k, 3] = zero(FT)
            # read in mass matrix inverse for element `e` in each plane k
            local_MI[k] = vgeo[ijk, _MI, e]
        end
        # end of boiler plate loading

        @unroll for k in 1:Nq3 # search
            @synchronize # perhaps moved up outside the loop
            ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1))
            # Extract Jacobian terms
            # ∂/∂x¹ = ∂ξ¹/ ∂x¹ ∂/∂ξ¹ + ∂ξ²/ ∂x¹ ∂/∂ξ² + ∂ξ³/ ∂x¹ ∂/∂ξ³
            # ∂/∂x² = ∂ξ¹/ ∂x² ∂/∂ξ¹ + ∂ξ²/ ∂x² ∂/∂ξ² + ∂ξ³/ ∂x² ∂/∂ξ³
            # ∂/∂x³ = ∂ξ¹/ ∂x³ ∂/∂ξ¹ + ∂ξ²/ ∂x³ ∂/∂ξ² + ∂ξ³/ ∂x³ ∂/∂ξ³
            # The integration by parts is then done with respect to the
            # (ξ¹, ξ², ξ³) coordinate system
            # (which is why terms associated with the columns appear)
            # maybe load gradient into shared memory
            shared_flux[1, i, j] = state[ijk, e, 1] # local_flux[1]

            ξ3x1 = vgeo[ijk, _ξ3x1, e]
            ξ3x2 = vgeo[ijk, _ξ3x2, e]
            ξ3x3 = vgeo[ijk, _ξ3x3, e]

            @unroll for n in 1:Nq3
                # ξ3-grid lines
                ∂Q∂ξ³ = D[n, k] * shared_flux[1, i, j]
                local_tendency[n, 1] += ξ3x1 * ∂Q∂ξ³
                local_tendency[n, 2] += ξ3x2 * ∂Q∂ξ³
                local_tendency[n, 3] += ξ3x3 * ∂Q∂ξ³
            end

            @synchronize
            ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1))
            ξ1x1 = vgeo[ijk, _ξ1x1, e]
            ξ1x2 = vgeo[ijk, _ξ1x2, e]
            ξ1x3 = vgeo[ijk, _ξ1x3, e]

            ξ2x1 = vgeo[ijk, _ξ2x1, e]
            ξ2x2 = vgeo[ijk, _ξ2x2, e]
            ξ2x3 = vgeo[ijk, _ξ2x3, e]

            @unroll for n in 1:Nq1
                # ξ1-grid lines
                ∂Q∂ξ¹ = D[i, n] * shared_flux[1, n, j]
                gradient[ijk, e, 1] += ξ1x1 * ∂Q∂ξ¹
                gradient[ijk, e, 2] += ξ1x2 * ∂Q∂ξ¹
                gradient[ijk, e, 3] += ξ1x3 * ∂Q∂ξ¹
            end

            @unroll for n in 1:Nq2
                ∂Q∂ξ² = D[j, n] * shared_flux[1, i, n]
                # ξ2-grid lines
                gradient[ijk, e, 1] += ξ2x1 * ∂Q∂ξ²
                gradient[ijk, e, 2] += ξ2x2 * ∂Q∂ξ²
                gradient[ijk, e, 3] += ξ2x3 * ∂Q∂ξ²
            end
        end

        @unroll for k in 1:Nq3
            # ξ3-grid lines
            ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1))
            gradient[ijk, e, 1] += local_tendency[k, 1]
            gradient[ijk, e, 2] += local_tendency[k, 2]
            gradient[ijk, e, 3] += local_tendency[k, 3]
        end
    end
end

"""
    launch_interface_gradients!(dg, state_prognostic, t; surface::Symbol, dependencies)

Launches horizontal and vertical kernels for computing the interface gradients.
The argument `surface` is either `:interior` or `:exterior`, which denotes whether
we are computing interface gradients on boundaries which are interior (exterior resp.)
to the _parallel_ boundary.
"""
function launch_interface_gradient!(
    gradient,
    state,
    grid::G;
    dependencies = nothing,
) where {G <: AbstractGrid}
    # MPI
    # @assert surface === :interior || surface === :exterior

    FT = eltype(state)
    device = array_device(state)
    dim = length(grid.D)
    N = polynomialorders(grid)

    # Assumes poly order is the same in every direction
    # TODO: copy from new kernel with variable polynomial order
    if dim == 1
        Np = (N[1] + 1)
        Nfp = 1
        nface = 2
    elseif dim == 2
        Np = (N[1] + 1) * (N[2] + 1)
        Nfp = (N[1] + 1)
        nface = 4
    elseif dim == 3
        Np = (N[1] + 1) * (N[2] + 1) * (N[3] + 1)
        Nfp = (N[1] + 1) * (N[2] + 1)
        nface = 6
    end
    gradient.realdata[:] .= 0
    workgroup = Nfp
    elems = grid.interiorelems
    ndrange = Nfp * length(grid.interiorelems)

    comp_stream = dependencies

    #TODO: fix with appropriate function calls
    # If the model direction is EveryDirection, we need to perform
    # both horizontal AND vertical kernel calls; otherwise, we only
    # call the kernel corresponding to the model direction `dg.diffusion_direction`
    comp_stream = interface_gradient_kernel!(device, workgroup)(
        gradient.data,
        state.data,
        Val(dim),
        Val(N[1]),
        grid.vgeo,
        grid.sgeo,
        grid.vmap⁻,
        grid.vmap⁺,
        grid.elemtobndy,
        elems,
        nface,
        Np,
        ndrange = ndrange,
        dependencies = comp_stream,
    )
    return comp_stream
end

@doc """
    function interface_gradient_kernel!(
        gradient,
        state,
        ::Val{dim},
        ::Val{polyorder},
        direction, # we will probably want to get rid of this
        numerical_state_first_order,  # this is metadata, we need to specify what numerical state we are using
        tendency,
        vgeo,
        sgeo,
        vmap⁻,
        vmap⁺,
        elemtobndy,
        elems,
    )

Compute kernel for evaluating the interface tendencies for the
DG form:

∫ₑ ψ⋅ ∂q/∂t dV - ∫ₑ ∇ψ⋅F⃗ dV + ∮ₑ ψF⃗⋅n̂ dS,

or equivalently in matrix form:

dQ/dt = M⁻¹(MS + DᵀM F⃗ + ∑ᶠ LᵀMᶠF⃗ ).

This kernel computes the surface terms: M⁻¹ ∑ᶠ LᵀMᶠ(F)),
where M is the mass matrix, Mf is the face mass matrix, L is an interpolator
from volume to face, and Fⁱⁿᵛ⋆, Fᵛⁱˢᶜ⋆
are the numerical statees for the inviscid and viscous
statees, respectively.  interface_tendency!
"""
@kernel function interface_gradient_kernel!(
    gradient,
    state,
    ::Val{dim},
    ::Val{polyorder},
    vgeo,
    sgeo,
    vmap⁻,
    vmap⁺,
    elemtobndy,
    elems,
    nface,
    Np,
) where {dim, polyorder}
    @uniform begin
        N = polyorder
        FT = eltype(state)

        # faces: (1, 2, 3, 4, 5, 6)  = (West, East, South, North, Bottom, Top)
        # We think...
        faces = 1:nface

        Nq = N + 1
        Nqk = dim == 2 ? 1 : Nq
    end
    eI = @index(Group, Linear)
    n = @index(Local, Linear)

    e = @private Int (1,)

    e[1] = elems[eI]
    # goes over too many points right now
    for f in faces
        # Inter-element index
        e⁻ = e[1]
        normal_vector = SVector(
            sgeo[_n1, n, f, e⁻],
            sgeo[_n2, n, f, e⁻],
            sgeo[_n3, n, f, e⁻],
        )
        # Get surface mass, volume mass inverse

        # sM = surface mass matrix
        # vMI = volume mass inverse matrix
        sM, vMI = sgeo[_sM, n, f, e⁻], sgeo[_vMI, n, f, e⁻]

        # Ids corresponding to the workitem on the face
        # vmap⁻ and vmap⁺ are the linear indices
        id⁻, id⁺ = vmap⁻[n, f, e⁻], vmap⁺[n, f, e⁻]

        # neighbor element index
        e⁺ = ((id⁺ - 1) ÷ Np) + 1 # this is why 1 indexing is failure

        # vid's are cartesian indices of the GL points
        # needs to be modified by the faces?
        vid⁻, vid⁺ = ((id⁻ - 1) % Np) + 1, ((id⁺ - 1) % Np) + 1

        # No Flux BC and Central Fluxes for the interfaces allows the kernel
        # to be extensible
        # TODO: make it no state on the boundary for now and add the state later
        # bctype = elemtobndy[f, e⁻]

        # gradient[vmap⁻] = (state[vmap⁻] + state[vmap⁺])/2
        # ∫dA F⃗⋅n̂
        # normal_vector[1] * (state[vmap⁻,1] + state[vmap⁺,1])/2
        # normal_vector[2] * (state[vmap⁻,2] + state[vmap⁺,2])/2
        # normal_vector[3] * (state[vmap⁻,3] + state[vmap⁺,3])/2
        # # ∂ᵗ∫dVρ = - ∫dS Φ⃗⋅n̂ = - state on boundary

        local_state_1 = normal_vector[1] * (state[vid⁺, e⁺] - state[vid⁻, e⁻])/2
        local_state_2 = normal_vector[2] * (state[vid⁺, e⁺] - state[vid⁻, e⁻])/2
        local_state_3 = normal_vector[3] * (state[vid⁺, e⁺] - state[vid⁻, e⁻])/2

        # Update RHS (in outer face loop): M⁻¹ Mfᵀ(Fⁱⁿᵛ⋆ ))
        gradient[vid⁻, e⁻, 1] += vMI * sM * local_state_1
        gradient[vid⁻, e⁻, 2] += vMI * sM * local_state_2
        gradient[vid⁻, e⁻, 3] += vMI * sM * local_state_3

        # @show vid⁻ e⁻ id⁻ local_state gradient[vid⁻, e⁻, 1] vMI sM f
        # Need to wait after even faces to avoid race conditions
        # Note: only need synchronization for F%2 == 0, but this crashed...
        @synchronize
    end
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using DelimitedFiles
function load_data(f :: String)
  a = open(readdlm, f)
  x2 = a[a[:,1] .== 2, 2:end]
  x3 = a[a[:,1] .== 3, 2:end]
  y2 = ones(size(x2,1))
  y3 = (-1)*ones(size(x3,1))
  x = vcat(x2, x3)
  x = hcat(ones(size(x,1)), x)
  y = vcat(y2, y3)
  (x, y)
end

function lr_train(x_train :: Array{Float64,2}, y_train :: Array{Float64,1})
  inv(transpose(x_train) * x_train) * (transpose(x_train) * y_train)
end

function lr_test(x_test :: Array{Float64,2}, y_test :: Array{Float64,1}, beta :: Array{Float64,1})
  y_result = x_test * beta
  y_compare = y_result .* y_test
  length(y_compare[y_compare .> 0]) / length(y_result)
end


text:
macro inferred_return_type(expr)
  types = Base.gen_call_with_extracted_types(Base.return_types, expr)
  quote
    @assert length($types)==1
    $types[1]
  end
end

#macro inferred_isa(expr)

"""
Impertinent copy of `Base.inferred` with type equality replaced with `isa`
"""
macro inferred_isa(ex)
    if Meta.isexpr(ex, :ref)
        ex = Expr(:call, :getindex, ex.args...)
    end
    Meta.isexpr(ex, :call)|| error("@inferred_isa requires a call expression")

    Base.remove_linenums!(quote
        $(if any(a->(Meta.isexpr(a, :kw) || Meta.isexpr(a, :parameters)), ex.args)
            # Has keywords
            args = gensym()
            kwargs = gensym()
            quote
                $(esc(args)), $(esc(kwargs)), result = $(esc(Expr(:call, _args_and_call, ex.args[2:end]..., ex.args[1])))
                inftypes = $(Base.gen_call_with_extracted_types(Base.return_types, :($(ex.args[1])($(args)...; $(kwargs)...))))
            end
        else
            # No keywords
            quote
                args = ($([esc(ex.args[i]) for i = 2:length(ex.args)]...),)
                result = $(esc(ex.args[1]))(args...)
                inftypes = Base.return_types($(esc(ex.args[1])), Base.typesof(args...))
            end
        end)
        @assert length(inftypes) == 1
        rettype = isa(result, Type) ? Type{result} : typeof(result)
        rettype <: inftypes[1] || error("return type $rettype is not a subtype of the inferred return type $(inftypes[1])")
        result
    end)
end


text:
using SimLynx
using Documenter

"""
Loops through the files in the examples folder and adds them (with any header comments) to the examples.md
markdown file.
"""
function generateExamples()
    f = open(joinpath(@__DIR__, "src/examples.md"), "w")
    write(
        f,
        "```@meta
    # NOTE this file is autogenerated, do not edit examples.md directly. To make an example, upload the .jl file to the examples folder. Header comments may be included at the top of the file using \"\"\" syntax
``` ",
    )
    write(f, "\n")
    write(f, "# Examples")
    write(
        f,
        "\nBelow are some useful examples to give you an idea of how this package can be leveraged. The code for these examples can also be found on Github
 in the `docs/examples` folder.",
    )
    write(f, "\n")
    for (root, dirs, files) in walkdir(joinpath(@__DIR__, "examples"))

        for file in files
            println(file)
            #extract title from example
            write(f, "\n")
            title = file
            title = replace(title, "_" => " ")
            title = replace(title, ".jl" => "")
            title = titlecase(title)
            title = "## " * title * "\n"
            write(f, title)
            #open each file and read contents
            opened = open(joinpath(@__DIR__, "examples/") * file)
            lines = readlines(opened, keep = true)
            index = 1
            #find doc string intro if exists
            if "\"\"\"\n" in lines
                index = findall(isequal("\"\"\"\n"), lines)[2]
                print(index)
                for i = 2:index-1
                    write(f, lines[i])
                end
                lines = lines[index+1:end]
            end

            write(f, "```julia")
            write(f, "\n")
            for line in lines
                write(f, line)
            end
            write(f, "\n")
            write(f, "```")
            close(opened)
        end
    end
    close(f)
end

generateExamples()

makedocs(;
    modules=[SimLynx],
    authors="Trystan Kaes <trystanblkaes@gmail.com>, Anthony Dupont <anthony.dupont@ucdenver.edu>, Doug Williams <milton.williams@ucdenver.edu>",
    repo="https://github.com/LynxUCD/SimLynx.jl/blob/{commit}{path}#L{line}",
    sitename="SimLynx.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://LynxUCD.github.io/SimLynx.jl",
        assets=String[],
    ),
    # devbranch = "master",
    # devurl = "dev",
    versions = ["stable" => "v^", "v#.#", "dev" => "dev"],
    pages=[
        "Home" => "index.md",
        "Discrete Simulations" => "discrete.md",
        "Continuous Simulations" => "continuous.md",
        "Examples" => "examples.md",
    ],
    push_preview=true,
)

deploydocs(;
    repo="github.com/LynxUCD/SimLynx.jl",
)


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 1, Zadanie 4.
#

fl = Float64

# define the step
δ = fl(2^(-52))

x = fl(1 + δ)

# go through all the consecutive numbers
for i in range(1, length=2^52-1)
    if fl(x * fl(1/x)) != 1
        # if the condition is fulfilled print out the number in question and break the loop
        println(
            bitstring(x), "\n",
            bitstring(fl(1/x)), "\n",
            fl(x * fl(1/x)) != 1, "\n",
            x, "\n",
            fl(1/x), "\n",
            )
        break
    end
    # otherwise continue
    global x = fl(x + fl(δ * i))
end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
"""
Utilities for storing and manipulating electric (vector) and scalar fields
"""
module VectorFields

export VectorField
export ScalarField
export ElectricField
export update!
export component

dimensional_files = (
    "1D/VectorFields1D.jl",
    "2D/VectorFields2D.jl",
    # "3D/VectorFields3D.jl",
    )

using ..Points
using Interpolations
using RecipesBase
using StaticArrays

"""
    VectorFields{N,M} <: AbstractMatrix{ComplexF64}

`N`is the dimension, `M` is the number of field components
"""
struct VectorField{N,M,T,TC} <: AbstractMatrix{T} # N = Dimension, M = number of field components
    positions::Vector{Point{N,Cartesian}}
    values::Matrix{T}
    start::Point{N,Cartesian}
    stop::Point{N,Cartesian}
    start_inds::Vector{Int}
    stop_inds::Vector{Int}
    components::SVector{M,TC}

    # function VectorField{N,M}(pos::Vector{TP}, val::AbstractMatrix), start::Point{N}, stop::Point{N}, start_inds::Vector, stop_inds::Vector) where {N,M,TP<:Point{N}}
    function VectorField{N,M}(pos::Vector{TP}, val::AbstractMatrix{T}) where {N,M,TP<:Point{N},T}
        start = Point(ntuple(i->minimum(map(p->p[i],pos)), N))::Point{N}
        stop  = Point(ntuple(i->maximum(map(p->p[i],pos)), N))::Point{N}
        start_inds = map(i->findfirst(isequal(start[i]),map(p->p[i],pos))::Int, 1:N)
        stop_inds = map(i->findfirst(isequal(stop[i]),map(p->p[i],pos))::Int, 1:N)
        n = size(val,1)
        components = map(i -> view(val,(i-1)n÷M+1:(i)n÷M,:), 1:M)
        n==M*length(pos) || throw("provided matrix has size(matrix,1)=$(size(val,1))≠$M*length(pos)=$(M*length(pos))")

        return new{N,M,T,eltype(components)}(pos,val,start,stop,start_inds,stop_inds,components)
    end
end

(f::VectorField)(i::Integer) = VectorField(f,view(f.values,:,i))

"""
    VectorField{M}(pos,val,start,stop,start_inds,stop_inds) -> field

`field` is an `M`-component vector field in `N` dimensions, where `pos::Vector{Point{N}}`.
`val` is a vector or matrix of field values.
`start` and `stop` are `Point{N}`'s that define the smallest and largest corners of the hyperrectangle in `N`-dimsensions that all the elements of `pos` reside in
`start_inds` and `stop_inds` give the indices
"""
VectorField{M}(pos::Vector{TP}, args...) where {M,TP<:Point{N}} where N = VectorField{N,M}(pos, args...)
VectorField{N,M}(pos::Vector{TP},val::AbstractVector) where {N,M,TP<:Point{N}} = VectorField{N,M}(pos,reshape(val,length(val),1))
VectorField{N,M}(pos::Vector{TP}, m::Integer, ::Type{T}=ComplexF64) where {N,M,TP<:Point{N},T} = VectorField{N,M}(pos,zeros(T,M*length(pos),m))

"""
    VectorField(field,vals) -> newfield

`newfield` has field values `vals` on the same sites as `field`
"""
VectorField(f::VectorField{N,M},vals::AbstractVecOrMat) where {N,M} = VectorField{M}(f,vals)
VectorField{M}(f::VectorField{N,M},vals::AbstractVecOrMat) where {N,M} = VectorField{N,M}(f.positions,vals)

"""
    update!(field,vals)

update `field::VectorField` with values `vals` in place.
"""
update!(f::VectorField,vals::AbstractVecOrMat) = copyto!(f.values,vals)

"""
    component(i,field)

view of `i`-th component of `field::VectorField`
"""
function component(i::Integer,f::VectorField{N,M}; verbose::Bool=true) where {N,M}
    1 ≤ i ≤ M || (if throw("attempt to access i=$i component of $M-component field") end)
    return f.components[i]
end

Base.conj(f::VectorField) = VectorField(f,conj(f.values))
Base.conj!(f::VectorField) = begin conj!(f.values); f end

fnames = (:+,:-)
for fname ∈ fnames
    @eval begin
        function Base.$fname(f1::VectorField{N,M},f2::VectorField{N,M}) where {N,M}
            f1.positions == f2.positions || throw("cannot binary operate `VectorField`s that do not share field points")
            return VectorField(f1,$fname(f1.values,f2.values))
        end
    end
end
Base.:*(f::VectorField,a::Number) = a*f
Base.:*(a::Number,f::VectorField) = VectorField(f,a*f.values)
Base.:/(f::VectorField,a::Number) = a\f
Base.:\(a::Number,f::VectorField) = VectorField(f,a\f.values)

Base.getindex(f::VectorField,inds...) = f.values[inds...]
Base.setindex!(f::VectorField,v,inds...) = setindex!(f.values,v,inds...)
fnames = (:firstindex, :lastindex, :length, :size, :eltype, :iterate)
for fname ∈ fnames @eval Base.$fname(f::VectorField) = $fname(f.values) end
Base.iterate(f::VectorField,state) = iterate(f.values,state)
Base.IndexStyle(::Type{<:VectorField}) = IndexLinear()

"""
    ScalarField{N} <: AbstractMatrix{ComplexF64}

Access a given field (mode) `i` with `S(i)`.
"""
ScalarField{N} = VectorField{N,1}
ScalarField(args...) = VectorField{1}(args...)

"""
    ElectricField{N} <: AbstractMatrix{ComplexF64}

`E.x`, `E.y`, `E.z` create views into components of electric fields.

Access a given field (mode) `i` with `E(i)`.
"""
ElectricField{N} = VectorField{N,3}
ElectricField(args...) = VectorField{3}(args...)

function Base.getproperty(s::ScalarField{N}, sym::Symbol) where N
    p = _getdimensionalproperty(s,sym)
    if isnothing(p)
        n,m = size(getfield(s,:values))
        if sym == :pos
            return getfield(s,:positions)
        elseif sym == :val
            return getfield(s,:values)
        else
            return getfield(s,sym)
        end
    else
        return p
    end
end

function Base.getproperty(e::ElectricField{N}, sym::Symbol) where N
    p = _getdimensionalproperty(e,sym)
    if isnothing(p)
        n,m = size(getfield(e,:values))
        if sym == :pos
            return getfield(e,:positions)
        elseif sym == :val
            return getfield(e,:values)
        elseif Base.sym_in(sym,(:x,:X))
            return view(getfield(e,:values),0n÷3+1:1n÷3,:)
        elseif Base.sym_in(sym,(:y,:Y))
            return view(getfield(e,:values),1n÷3+1:2n÷3,:)
        elseif Base.sym_in(sym,(:z,:Z))
            return view(getfield(e,:values),2n÷3+1:3n÷3,:)
        else
            return getfield(e,sym)
        end
    else
        return p
    end
end

_getdimensionalproperty(v::VectorField{N,M}, sym::Symbol) where {N,M} = nothing

foreach(include,dimensional_files)

################################################################################
# extending INTERPOLATIONS
#
# const AbInt = AbstractInterpolation
#
# """
# Interpolated electric field.
#
# Access interpolations through `E.x`, `E.y`, `E.z`. Evaluate with `E.x(::Real)`, etc.
# """
# struct InterpolatedVectorField{N,M,TIX,TIY,TIZ}
#     positions::Vector{Point{N}}
#     x_interpolation::TIX
#     y_interpolation::TIY
#     z_interpolation::TIZ
#     InterpolatedVectorField{M}(pos::Vector{Point{N}},xi::TIX,yi::TIY,zi::TIZ) where {N,TIX<:AbInt,TIY<:AbInt,TIZ<:AbInt} =
#             new{N,M,TIX,TIY,TIZ}(pos,xi,yi,zi)
# end
#
# InterpolatedScalarField{N} = InterpolatedVectorField{N,1}
# InterpolatedVectorField{N} = InterpolatedVectorField{N,1}
#
# function Base.getproperty(ie::InterpolatedElectricField,sym::Symbol)
#     if sym==:x
#         return getfield(ie,:x_interpolation)
#     elseif sym==:y
#         return getfield(ie,:y_interpolation)
#     elseif sym==:z
#         return getfield(ie,:z_interpolation)
#     elseif sym==:pos
#         return getfield(ie,:positions)
#     else
#         return getfield(ie,sym)
#     end
# end
#
# function Base.propertynames(::InterpolatedElectricField{N},private=false) where N
#     if private
#         return fieldnames(InterpolatedElectricField{N})
#     else
#         return propertynames(ElectricField{N})
#     end
# end

################################################################################
# Pretty Printing

import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_DARK
import ..PRINTED_COLOR_INSTRUCTION
import ..PRINTED_COLOR_VARIABLE

function Base.show(io::IO,mime::MIME"text/plain",s::VectorField{N,M}) where {N,M}
    ds = displaysize(io)
    ioc = IOContext(io, :displaysize => (ds[1],ds[2]), :typeinfo => ComplexF64)
    show(io,s)
    println(io)
    show(ioc,mime,s.values)
end

function Base.show(io::IO,mime::MIME"text/plain",s::ScalarField)
    ds = displaysize(io)
    ioc = IOContext(io, :displaysize => (ds[1],ds[2]))
    show(io,s)
    println(io)
    show(ioc,mime,s.values)
end

function Base.show(io::IO,mime::MIME"text/plain",e::ElectricField)
    ds = displaysize(io)
    ioc = IOContext(io, :displaysize => (ds[1]÷3,ds[2]))
    show(io,e)
    println(io)
    printstyled(io,"\nX-component\n",color=PRINTED_COLOR_VARIABLE)
    show(ioc,mime,e.x)
    printstyled(io,"\nY-component\n",color=PRINTED_COLOR_VARIABLE)
    show(ioc,mime,e.y)
    printstyled(io,"\nZ-component\n",color=PRINTED_COLOR_VARIABLE)
    show(ioc,mime,e.z)
end

function Base.show(io::IO,f::VectorField{N,M}) where {N,M}
    n,m = size(f)
    if M==1
        printstyled(io,"ScalarField",color=PRINTED_COLOR_DARK)
    elseif M==3
        printstyled(io,"ElectricField",color=PRINTED_COLOR_DARK)
    else
        printstyled(io,"$M-component VectorField",color=PRINTED_COLOR_DARK)
    end
    print(io," (")
    if m > 1
        printstyled(io,m,color=PRINTED_COLOR_NUMBER)
        print(io," field")
        m > 1 ? print(io,"s, ") : nothing
    end
    printstyled(io,n÷M,color=PRINTED_COLOR_NUMBER)
    print(io," points)")
end

# function Base.show(io::IO,ie::InterpolatedScalarField{N}) where N
#     print(io,"$(N)D ")
#     printstyled(io,"Interpolated ScalarField ",color=PRINTED_COLOR_DARK)
#     printstyled(io,"(access fields ",color=PRINTED_COLOR_INSTRUCTION)
#     printstyled(io,".values",color=PRINTED_COLOR_VARIABLE)
#     printstyled(io,")",color=PRINTED_COLOR_INSTRUCTION)
# end
#
# function Base.show(io::IO,ie::InterpolatedElectricField{N}) where N
#     print(io,"$(N)D ")
#     printstyled(io,"Interpolated VectorField ",color=PRINTED_COLOR_DARK)
#     printstyled(io,"(access fields ",color=PRINTED_COLOR_INSTRUCTION)
#     printstyled(io,".x",color=PRINTED_COLOR_VARIABLE)
#     printstyled(io,", ",color=PRINTED_COLOR_INSTRUCTION)
#     printstyled(io,".y",color=PRINTED_COLOR_VARIABLE)
#     printstyled(io,", ",color=PRINTED_COLOR_INSTRUCTION)
#     printstyled(io,".z",color=PRINTED_COLOR_VARIABLE)
#     printstyled(io,")",color=PRINTED_COLOR_INSTRUCTION)
# end

@recipe f(e::VectorField; by=abs2) = e, by
@recipe f(by::Function, e::VectorField) = e,by

end #module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties <: SwaggerModel
    groupId::Any # spec type: Union{ Nothing, String } # spec name: groupId
    requiredMemberName::Any # spec type: Union{ Nothing, String } # spec name: requiredMemberName
    fqdns::Any # spec type: Union{ Nothing, Vector{String} } # spec name: fqdns

    function NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties(;groupId=nothing, requiredMemberName=nothing, fqdns=nothing)
        o = new()
        validate_property(NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties, Symbol("groupId"), groupId)
        setfield!(o, Symbol("groupId"), groupId)
        validate_property(NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties, Symbol("requiredMemberName"), requiredMemberName)
        setfield!(o, Symbol("requiredMemberName"), requiredMemberName)
        validate_property(NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties, Symbol("fqdns"), fqdns)
        setfield!(o, Symbol("fqdns"), fqdns)
        o
    end
end # type NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties

const _property_map_NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties = Dict{Symbol,Symbol}(Symbol("groupId")=>Symbol("groupId"), Symbol("requiredMemberName")=>Symbol("requiredMemberName"), Symbol("fqdns")=>Symbol("fqdns"))
const _property_types_NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties = Dict{Symbol,String}(Symbol("groupId")=>"String", Symbol("requiredMemberName")=>"String", Symbol("fqdns")=>"Vector{String}")
Base.propertynames(::Type{ NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties }) = collect(keys(_property_map_NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties))
Swagger.property_type(::Type{ NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties[name]))}
Swagger.field_name(::Type{ NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties }, property_name::Symbol) =  _property_map_NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties[property_name]

function check_required(o::NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties)
    true
end

function validate_property(::Type{ NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties }, name::Symbol, val)
end


text:
coord1 = EuclideanCoordinate(1, 2)
coord2 = EuclideanCoordinate(1, 3)

station1 = EuclideanStation(1, coord1, "test_2", false)
station2 = EuclideanStation(2, coord2, "test", false)
label_node3 = EuclideanStation(3, coord2, "test", true)
label_node4 = EuclideanStation(4, coord2, "test", true)

line1 = Line("S1")
line2 = Line("DummyEdge2")
line3 = Line("DummyEdge3")
edge = ProcessedEdge(station1, station2, line1, 0, 1, false)

edge2 = ProcessedEdge(station1, label_node3, line2, 0, 7, true)
edge3 = ProcessedEdge(station2, label_node4, line2, 0, 7, true)

transit_map = InputGraph([station1, station2, label_node3, label_node4],
    [edge, edge2, edge3],
    [line1, line2, line3])

result = optimize(GLPKSolverMIP(msg_lev = 3), transit_map, false)

@test result.stations[1].coordinate.y ≈ result.stations[2].coordinate.y

# label one should not be in sector 2
@test !(result.stations[3].coordinate.x ≈ result.stations[2].coordinate.x &&
    result.stations[3].coordinate.y > result.stations[2].coordinate.y)

# label one should not be in sector 6
@test !(result.stations[3].coordinate.x ≈ result.stations[2].coordinate.x &&
    result.stations[3].coordinate.y < result.stations[2].coordinate.y)

max_stars_repo_name: brainfo/MetroMap

max_issues_repo_name: brainfo/MetroMap

max_forks_repo_name: brainfo/MetroMap


text:
function ρ(X, X_sim; metric = "absolute")
    """
    Calculate the sum of the absolute difference at each time point between observed data
    X and simulated data X_sim noting that each row of these corresponds to a realisation.
    that
    """

    ϵ = zeros(size(X, 2))

    if metric == "absolute"
        ϵ .= sum(abs.(X - X_sim), dims = 1)[:]
    elseif metric == "MSE"
        ϵ .= sum((X - X_sim) .^ 2, dims = 1)[:]
    elseif metric == "RMSE"
        N = size(X, 1)
        ϵ .= sqrt.(1 / N * sum((X - X_sim) .^ 2, dims = 1)[:])
    end

    return ϵ

end

function final_check(D_good, local_cases; p = 0.1, metric = "absolute")
    """
    Determines the top p×100% of simulations based on some metric ρ(X, X_sim) and returns
    those indices.
    """
    # get only the observed cases over the period we have cases for
    D_observed = D_good[begin:length(local_cases), 1, :] +
        D_good[begin:length(local_cases), 2, :]
    # tile local_cases to easily calculate the distance ρ(X, X_sim)
    local_cases_mat = repeat(local_cases, 1, size(D_observed, 2))

    # calculate distance between data and the simulations
    ϵ = ρ(local_cases_mat, D_observed, metric = metric)
    # determine the sim indices of the top p×100%
    good_inds = ϵ .< quantile(ϵ, p)

    return good_inds

end

function count_cases2!(
    case_counts,
    forecast::Forecast,
    sim,
)

    D = forecast.sim_realisations.D

    day = 1
    tmp = 0

    for i in 1:length(case_counts)
        case_counts[day] = D[i, 1, sim] + D[i, 2, sim]
        day += 1
    end

    return nothing

end

function calculate_bounds2(local_cases, τ, state)

    # tolerance for when we consider cases to be stable
    ϵ = 0.0

    # multipliers on the n-day average
    (ℓ, u) = (0.5, 2.5)

    # observation period
    T = length(local_cases)
    # the slope over the n-day period
    m = 0.0
    Mₜ = zeros(Float64, T)

    Lₜ = zeros(Float64, T)
    Uₜ = zeros(Float64, T)

    # consider τ = 3 and t = (0, 1, 2), clearly n = 2 - 0 = 2
    n = τ - 1

    n2 = 7

    for t in 1:T

        # if state == "NSW"
        #     if t < 30
        #         (ℓ, u) = (0, 10)
        #     elseif t < T - 7
        #         (ℓ, u) = (0.5, 2.0)
        #     else
        #         (ℓ, u) = (0.5, 2.0)
        #     end
        # elseif state == "NT"
        #     if t < 30
        #         (ℓ, u) = (0, 4.0)
        #     elseif t < T - 7
        #         (ℓ, u) = (0.25, 4.0)
        #     else
        #         (ℓ, u) = (0.5, 3.0)
        #     end
        # end

        if t < T - 7
            (ℓ, u) = (0.25, 2.0)
        else
            (ℓ, u) = (0.5, 2.0)
        end

        # approximate the slope naively
        if t < T - n
            m = 1 / n * (
                local_cases[t + n] -
                local_cases[t]
            )
        else
            m = 1 / n * (
                local_cases[t] -
                local_cases[t - n]
            )
        end

        Mₜ[t] = m

        # depending on the sign of the slope, take average of future or past cases
        if m < -ϵ
            Lₜ[t] = ℓ * mean(
                local_cases[t:min(T, t + n)]
            )
            Uₜ[t] = u * mean(
                local_cases[max(1, t - n):t]
            )
            # Lₜ[t] = mean(local_cases[t:min(T, t + n)])
            # Uₜ[t] = mean(local_cases[max(1, t - n):t])
        elseif m > ϵ
            Lₜ[t] = ℓ * mean(
                local_cases[max(1, t - n):t]
            )
            Uₜ[t] = u * mean(
                local_cases[t:min(T, t + n)]
            )
            # Lₜ[t] = mean(local_cases[max(1, t - 2):t])
            # Uₜ[t] = mean(local_cases[t:min(T, t + 2)])
        else
            n2 = n ÷ 2
            Lₜ[t] = ℓ * mean(
                local_cases[
                    max(1, t - n2):min(T, t + n2)
                ]
            )
            Uₜ[t] = u * mean(
                local_cases[
                    max(1, t - n2):min(T, t + n2)
                ]
            )
            # Lₜ[t] = mean(local_cases[max(1, t - n2):min(T, t + n2)])
            # Uₜ[t] = mean(local_cases[max(1, t - n2):min(T, t + n2)])
        end

        # adjust the bounds for periods with low cases
        if Lₜ[t] < 50
            Lₜ[t] = 0
        end

        if Uₜ[t] < 50
            Uₜ[t] = 50
        end

        if isnan(Lₜ[t])
            Lₜ[t] = Lₜ[t-1]
        end

        if isnan(Uₜ[t])
            Uₜ[t] = Uₜ[t-1]
        end
    end

    return (Lₜ, Uₜ)

end


text:
using Test
using SCS

tests = ["direct.jl",
         "indirect.jl",
         "options.jl"]

include("test_problems.jl")

for curtest in tests
    @testset "$curtest" begin
        include(curtest)
    end
end

include("MPB_wrapper.jl")
include("MOI_wrapper.jl")


text:
using PolynomialAmoebas
using Test

@testset "PolynomialAmoebas" begin
    include("fibers.jl")
    include("membership_test.jl")
    include("newton_test.jl")

    include("newton_polygon_test.jl")
    include("tropical_test.jl")
    include("tropicalcurve_test.jl")

    include("simple_grid_test.jl")
    include("amoeba_archtrop_test.jl")
    include("contour_test.jl")
    include("greedy_grid_test.jl")
    include("spine_test.jl")
    include("polygonal_amoeba_test.jl")
end


text:
using ExecutableSpecifications

exitcode = runspec() ? 0 : - 1
exit(exitcode)


text:
# approximation algorithm from ``Algorithm 917: Complex Double-Precision
# Evaluation of the Wright ω Function'' by Piers W. Laurence, Robert M. Corless
# and David J. Jeffrey
# Series approximations can be found in the book chapter ``The Wright ω Function''
# by Corless, R. M. and Jeffrey, D. J.

function wrightomega(z::T) where {T <: Real}
    z = float(z)
    if z < -2
        t = exp(z)
        w = t * (1 + t * (-1 + t * (T(3) / 2 + t * (-T(8) / 3 + T(125) / 24 * t))))
    elseif z < T(1) + π
        z1 = z - 1
        w = 1 + z1 / 2 * (1 + z1 / 8 * (1 + z1 / 12 * (-1 + z1 / 16 *
            (-1 + z1 * T(13) / 20))))
    else
        lz = log(z)
        w = z + lz * (-1 + (1 + (lz / 2 - 1 + (lz * (lz / 3 - T(3) / 2) + 1) /
            z) / z) / z)
    end
    # refinement steps
    for _ in 1:5
        r = z - w - log(w)
        w1 = w + 1
        t = w1 * (w1 + T(2) / 3 * r)
        w *= 1 + r / w1 * (t - r / 2) / (t - r)
        fscn = abs(r^4 * (2 * w * (w - 4) - 1))
        if fscn < eps(float(T)) * 72 * w1^6
            break
        end
    end
    return w::float(T)
end

function rootnewton(
    f::Function,
    g::Function;
    lower::T = -Inf,
    upper::T = Inf,
    init::T = (lower + upper) / 2,
    increasing::Bool = true,
    ) where {T <: Real}
    curr = init
    f_new = f(curr)
    iter = 0
    while abs(f_new) > 1000eps()
        candidate = curr - f_new / g(BigFloat(curr))
        if (candidate < lower) || (candidate > upper)
            curr = (lower + upper) / 2
        else
            curr = candidate
        end
        f_new = f(curr)
        if (f_new < 0 && increasing) || (f_new >= 0 && !increasing)
            lower = curr
        else
            upper = curr
        end
        iter += 1
        if iter > 200
            @warn "too many iters in dual grad"
            break
        end
    end
    # @show iter
    return (T(curr), iter)
end

max_stars_repo_name: lkapelevich/conjgradexperiments

max_issues_repo_name: lkapelevich/conjgradexperiments

max_forks_repo_name: lkapelevich/conjgradexperiments


text:
import PredictMD
import Test

Test.@test( isa(Test.@inferred(PredictMD.version_codename()), String) )

Test.@test( length(Test.@inferred(PredictMD.version_codename())) > 0 )

Test.@test(
    Test.@inferred(PredictMD.version_codename()) ==
        Test.@inferred(PredictMD.version_codename(PredictMD))
    )

Test.@test(
    Test.@inferred(PredictMD.version_codename()) ==
        Test.@inferred(
            PredictMD.version_codename(first(methods(PredictMD.eval)))
            )
    )

Test.@test(
    Test.@inferred(PredictMD.version_codename()) ==
        Test.@inferred(PredictMD.version_codename(PredictMD.eval))
    )

Test.@test(
    Test.@inferred(PredictMD.version_codename()) ==
        Test.@inferred(
            PredictMD.version_codename(PredictMD.eval, (Any,),)
            )
    )

Test.@test(
    Test.@inferred(PredictMD.version_codename(TestModuleA)) ==
        "Code name for TestModuleA"
    )

Test.@test(
    Test.@inferred(PredictMD.version_codename(TestModuleB)) ==
        "TestModuleB code name"
    )

Test.@test_throws(
    ErrorException,
    PredictMD.TomlFile(joinpath(PredictMD.maketempdir(),"1","2","3","4")),
    )


text:
# test pd matrix types
using LinearAlgebra, PDMats, SparseArrays, SuiteSparse
using Test

call_test_pdmat(p::AbstractPDMat,m::Matrix) = test_pdmat(p,m,cmat_eq=true,verbose=1)

for T in [Float64, Float32]
    #test that all external constructors are accessible
    m = Matrix{T}(I, 2, 2)
    @test PDMat(m, cholesky(m)).mat == PDMat(Symmetric(m)).mat == PDMat(m).mat == PDMat(cholesky(m)).mat
    d = ones(T,2)
    @test PDiagMat(d,d).inv_diag == PDiagMat(d).inv_diag
    x = one(T)
    @test ScalMat(2,x,x).inv_value == ScalMat(2,x).inv_value
    s = SparseMatrixCSC{T}(I, 2, 2)
    @test PDSparseMat(s, cholesky(s)).mat == PDSparseMat(s).mat == PDSparseMat(cholesky(s)).mat

    #test the functionality
    M = convert(Array{T,2}, [4. -2. -1.; -2. 5. -1.; -1. -1. 6.])
    V = convert(Array{T,1}, [1.5, 2.5, 2.0])
    X = convert(T,2.0)

    call_test_pdmat(PDMat(M), M) #tests of PDMat
    call_test_pdmat(PDiagMat(V), Matrix(Diagonal(V))) #tests of PDiagMat
    call_test_pdmat(ScalMat(3,x), x*Matrix{T}(I, 3, 3)) #tests of ScalMat
    call_test_pdmat(PDSparseMat(sparse(M)), M)
end

m = Matrix{Float32}(I, 2, 2)
@test convert(PDMat{Float64}, PDMat(m)).mat == PDMat(convert(Array{Float64}, m)).mat
@test convert(AbstractArray{Float64}, PDMat(m)).mat == PDMat(convert(Array{Float64}, m)).mat
m = ones(Float32,2)
@test convert(PDiagMat{Float64}, PDiagMat(m)).diag == PDiagMat(convert(Array{Float64}, m)).diag
@test convert(AbstractArray{Float64}, PDiagMat(m)).diag == PDiagMat(convert(Array{Float64}, m)).diag
x = one(Float32); d = 4
@test convert(ScalMat{Float64}, ScalMat(d, x)).value == ScalMat(d, convert(Float64, x)).value
@test convert(AbstractArray{Float64}, ScalMat(d, x)).value == ScalMat(d, convert(Float64, x)).value
s = SparseMatrixCSC{Float32}(I, 2, 2)
@test convert(PDSparseMat{Float64}, PDSparseMat(s)).mat == PDSparseMat(convert(SparseMatrixCSC{Float64}, s)).mat

# no-op conversion with correct eltype (#101)
X = PDMat((Y->Y'Y)(randn(Float32, 4, 4)))
@test convert(AbstractArray{Float32}, X) === X
@test convert(AbstractArray{Float64}, X) !== X

# type stability of whiten! and unwhiten!
a = PDMat([1 0.5; 0.5 1])
@inferred whiten!(ones(2), a, ones(2))
@inferred unwhiten!(ones(2), a, ones(2))
@inferred whiten(a, ones(2))
@inferred unwhiten(a, ones(2))


text:
using Base.Test
using Klara

println("    Testing logistic function...")

@test logistic(0.7, 3, 4, 2.1, 1.4) == 1.4110527196983078


text:
using Seismic
using Base.Test

d,e = SeisLinearEvents()
ddec = SeisDecimate(d,perc=50)

dpocs = SeisPOCS(ddec,dt=0.004,fmax=100,Niter=100,p=1.5)

# test that quality factor is greater than 10 Decibels
quality_factor = 10*log10(norm(d[:],2)/norm(dpocs[:]-d[:],2))
println("Quality factor = ",quality_factor)
@test quality_factor > 5

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
@testset "Estimation ........................... " begin

  function simonsDesign(r1, n1, r, n)
    nvec = [[n1 for x1 in 0:r1]; [n for x1 in (r1 + 1):n1]]
    cvec = [[Inf for x1 in 0:r1]; [r for x1 in (r1 + 1):n1]]
    return Design(nvec, cvec)
  end

  # Simon's designs for beta = .2, alpha = .05, p1 = p0 +0.2
  # cf. Simon, R "Optimal Two-Stage Designs for Phase II Clinical Trials",
  # Controlled Clinical Trials 10:1-10 (1989). (p. 4)
  sd = simonsDesign( 7, 16, 23, 46) # p0 = 0.4

  mle    = MLEstimator(sd)
  ube    = RBEstimator(sd)
  supp   = support(sd)
  @test all(
    estimate.(mle, supp[:, 1], supp[:, 2]) .==
      (supp[:, 1] + supp[:, 2])./samplesize.(sd, supp[:, 1])
  )
  @test maximum(abs.(bias.(ube, linspace(0, 1)))) ≈ 0.0 atol=1e-10

  ss = SampleSpace(5:20, 50)

  params = MESS(
    ss, .2, .4, alpha = .05, beta = .2, pess = .4
  )

  design = optimaldesign(10, params, solver)

  cpe    = OCEstimator(design, qsolver)

  @test maximum(rmse.(cpe, linspace(0, 1, 10))) <= 0.12

  @test length(incompatibleoutcomes(cpe)) == 0

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
function fibsq08()
#
#  test polynomial suggested by Goedecker
#
    n = 8;
    p = fib(n)
    p = conv(p,p)
    z = [-.8762862300182460;
         -.6416053894725764-.6063952206057722*im; -.6416053894725764+.6063952206057722*im;
         -.04694032461065018-.9030234661229178*im; -.04694032461065018+.9030234661229178*im;
         .6286732392246423-.7084725692273577*im; .6286732392246423+.7084725692273577*im;
         1.996031179735415];
    z = [z 2*ones(n)]
    p, PolyZeros(z)
end


text:
const metre = 1e2
const meter = metre
const cm = metre / 1e2
const mm = metre / 1e3
const um = metre / 1e6
const nm = metre / 1e9

const second = 1e3
const ms = second / 1e3
const Hz = 1 / second

const voltage = 1e3
const mV = voltage / 1e3

const ampere = 1e6
const mA = ampere / 1e3
const uA = ampere / 1e6
const nA = ampere / 1e9

const farad=1e6
const uF = farad/1e6
const ufarad = uF

const siemens = 1e3
const mS = siemens / 1e3
const msiemens = mS
const nS = siemens / 1e9
const nsiemens = nS

macro load_units()
    exs = map((:metre, :meter, :cm, :mm, :um, :nm, :second, :ms, :Hz, :voltage, :mV, :ampere,
        :mA, :uA, :nA, :farad, :uF, :ufarad, :siemens, :mS, :msiemens, :nS, :nsiemens)) do s
        :($s = getfield($@__MODULE__, $(QuoteNode(s))))
    end
    ex = Expr(:block, exs...)
    esc(ex)
end


text:
abstract type AbstractTasklet{T} end

struct Tasklet{T} <: AbstractTasklet{T}
    thunk::OpaqueClosure{Tuple{},T}
    promise::Promise{T}
end

Tasklet(thunk::OpaqueClosure{Tuple{},T}) where {T} = Tasklet{T}(thunk, Promise{T}())

#=
struct TypedTasklet{T,Thunk} <: AbstractTasklet{T}
    thunk::Thunk
    promise::Promise{T}
end

TypedTasklet{T}(thunk::Thunk) where {T,Thunk} = TypedTasklet{T,Thunk}(thunk, Promise{T}())
=#

macro tasklet(thunk)
    thunk = Expr(:block, __source__, thunk)
    ex = :($Tasklet($Base.Experimental.@opaque () -> $thunk))
    return esc(ex)
end

(tasklet::AbstractTasklet)() = race_put_with!(tasklet.thunk, tasklet.promise)
Base.fetch(tasklet::AbstractTasklet) = fetch(tasklet.promise)
Base.wait(tasklet::AbstractTasklet) = wait(tasklet.promise)
ConcurrentUtils.try_race_fetch(tasklet::AbstractTasklet) = try_race_fetch(tasklet.promise)


text:
# Helper functions
# last version of Julia does not have support for the
using SpecialFunctions

# reference Helmholtz solution for a centered Gaussian
@inline function solRefHelmholtz(x::Array{Float64,1},y::Array{Float64,1},z::Array{Float64,1}, sigma::Float64, k::Float64)
  r  = sqrt( x.^2 + y.^2 + z.^2)
  u  = (exp.(-sigma^2*k^2/2)./(4*pi.*r)).*(real( exp.(-1im.*k.*r).*erf( (2*sigma^2*1im*k - 2.*r)./(2*sqrt(2*sigma^2) ) )) - 1im*sin(k.*r))
  return u
end

@inline function Gtruncated2D{T<:Number,N}(L::Float64,k::Float64,s::Array{T,N})
    return (1 + (1im*pi/2*L*hankelh1.(0,L*k))*(s.*besselj.(1,L*s)) - (1im*pi/2*L*k*hankelh1.(1,L*k))*besselj.(0,L*s)  )./(s.^2 - k^2)
end

@inline function Gtruncated3D{T<:Number,N}(L::Float64,k::Float64,s::Array{T,N})
    return (-1 + exp.(1im*L*k).*( cos(L*s) - (1im*k*L*sinc(L*s/pi)) ) )./( (k^2  - s.^2))
end

@inline function Gtruncated3D(L::Float64,k::Float64,s::Float64)
    return (-1 + exp.(1im*L*k)*( cos(L*s) - (1im*k*L*sinc(L*s/pi)) ))/( (k^2 - s^2))
end

# This is Not used in the release version
type changeVariable
    coef::Array{Float64,1}
    order::Int64
    delta::Float64

    function changeVariable(deltta::Float64)
        A = [   delta^5    delta^4   delta^3;
              5*delta^4  4*delta^3 3*delta^2;
             20*delta^3 12*delta^2 6*delta^1];

        val = A\[ delta^2 0 0].'
        coeff= vcat( val, [0 0 0].')
        return new(coeff, 5, delta)
    end
end

@inline function changeVar(changeVariable, x)
    # function for the change of variable to integrate the smoth cut-off
    # the derivative at 0 is zero and at delta is equal to 1.
    xInt = x.*(x.<changeVariable.delta)
    powInt = copy(xInt)
    ans = zeros(xInt)
    ans[fin(x.<changeVariable.delta)] = changeVariable.coef[end]

    for ii = 1:length(changeVariable.coef)-1
        ans += changeVariable.coef[end-ii]*powInt;
        powInt = powInt.*xInt;
    end

    ans[fin(x.>=changeVariable.delta)] =  x[x.<changeVariable.delta] - (changeVariable.delta)^2;

    return ans
end

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
module MaxHelpingHandFlagTouchSwitch

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/FlagTouchSwitch" FlagTouchSwitch(x::Integer, y::Integer,
    flag::String="flag_touch_switch", icon::String="vanilla", persistent::Bool=false,
    inactiveColor::String="5FCDE4", activeColor::String="FFFFFF", finishColor::String="F141DF", smoke::Bool=true,
    inverted::Bool=false, allowDisable::Bool=false, hitSound::String="event:/game/general/touchswitch_any",
    completeSoundFromSwitch::String="event:/game/general/touchswitch_last_cutoff", completeSoundFromScene::String="event:/game/general/touchswitch_last_oneshot")

const bundledIcons = String["vanilla", "tall", "triangle", "circle", "diamond", "double", "heart", "square", "wide", "winged"]

const placements = Ahorn.PlacementDict(
    "Flag Touch Switch (max480's Helping Hand)" => Ahorn.EntityPlacement(
        FlagTouchSwitch
    )
)

Ahorn.editingOrder(entity::FlagTouchSwitch) = String["x", "y", "width", "height", "inactiveColor", "activeColor", "finishColor", "hitSound", "completeSoundFromSwitch", "completeSoundFromScene"]

Ahorn.editingOptions(entity::FlagTouchSwitch) = Dict{String,Any}(
    "icon" => bundledIcons
)

function Ahorn.selection(entity::FlagTouchSwitch)
    x, y = Ahorn.position(entity)

    return  Ahorn.Rectangle(x - 7, y - 7, 14, 14)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagTouchSwitch, room::Maple.Room)
    Ahorn.drawSprite(ctx, "objects/touchswitch/container.png", 0, 0)

    icon = get(entity.data, "icon", "vanilla")

    iconPath = "objects/touchswitch/icon00.png"
    if icon != "vanilla"
        iconPath = "objects/MaxHelpingHand/flagTouchSwitch/$(icon)/icon00.png"
    end

    Ahorn.drawSprite(ctx, iconPath, 0, 0)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
using JuMP, EAGO, Cbc

function make_model()
    sb = SubSolvers(;relaxed_optimizer= Cbc.Optimizer())
    m = Model(() -> EAGO.Optimizer(subsolver_block = sb))
    set_optimizer_attribute(m, "output_iterations", 1)
    set_optimizer_attribute(m, "iteration_limit", 5)
    set_optimizer_attribute(m, "verbosity", 0)

    # Define bounded variables
    xL = [10.0; 0.0; 0.0; 0.0; 0.0; 85.0; 90.0; 3.0; 1.2; 145.0]
    xU = [2000.0; 16000.0; 120.0; 5000.0; 2000.0; 93.0; 95.0; 12.0; 4.0; 162.0]
    @variable(m, xL[1] <= x[1] <= xU[1])
    @NLobjective(m, Max, -2.0 * x[1])
    return m
end

function f(m)
    JuMP.optimize!(m)
end

m = make_model()
m1 = make_model()

println("EAGO first call ")
@time f(m)
@time f(m1)

#=
m1 = Model(EAGO.Optimizer)
set_optimizer_attribute(m1, "output_iterations", 1)
set_optimizer_attribute(m1, "iteration_limit", 5)
set_optimizer_attribute(m1, "verbosity", 0)

# Define bounded variables
xL = [10.0; 0.0; 0.0; 0.0; 0.0; 85.0; 90.0; 3.0; 1.2; 145.0]
xU = [2000.0; 16000.0; 120.0; 5000.0; 2000.0; 93.0; 95.0; 12.0; 4.0; 162.0]
@variable(m1, xL[1] <= x[1] <= xU[1])

@constraint(m, e2, -x[1]+1.22*x[4]-x[5] == 0.0)
@constraint(m, e6, x[9]+0.222*x[10] == 35.82)
@constraint(m, e7, -3*x[7]+x[10] == -133.0)

# Define nonlinear constraints
@NLconstraint(m, e1, -x[1]*(1.12+0.13167*x[8]-0.00667* (x[8])^2)+x[4] == 0.0)
@NLconstraint(m, e3, -0.001*x[4]*x[9]*x[6]/(98-x[6])+x[3] == 0.0)

###Original constraints
@NLconstraint(m, e4, -(1.098*x[8]-0.038* (x[8])^2)-0.325*x[6]+x[7] == 57.425)
@NLconstraint(m, e5, -(x[2]+x[5])/x[1]+x[8] == 0.0)

@NLobjective(m1, Max, -2.0 * x[1])
JuMP.optimize!(m1)
=#
#=
=#
# Define nonlinear constraints
#@NLconstraint(m, e1, -x[1]*(1.12+0.13167*x[8]-0.00667* (x[8])^2)+x[4] == 0.0)
#@NLconstraint(m, e3, -0.001*x[4]*x[9]*x[6]/(98-x[6])+x[3] == 0.0)

###Original constraints
#@NLconstraint(m, e4, -(1.098*x[8]-0.038* (x[8])^2)-0.325*x[6]+x[7] == 57.425)
#@NLconstraint(m, e5, -(x[2]+x[5])/x[1]+x[8] == 0.0)

###Rewritten constraints
#=
@NLexpressions(m, begin
 ex1, -(x[2]+x[5])/x[1]+x[8]
 ex2, -(1.098*x[8]-0.038* (x[8])^2)-0.325*x[6]+x[7]
end)
@NLconstraint(m, e4, ex2 == 57.425)
@NLconstraint(m, e5, ex1 == 0.0)
=#

# Define linear constraints
#@constraint(m, e2, -x[1]+1.22*x[4]-x[5] == 0.0)
#@constraint(m, e6, x[9]+0.222*x[10] == 35.82)
#@constraint(m, e7, -3*x[7]+x[10] == -133.0)

# Define nonlinear objective
#@NLobjective(m, Max, 0.063*x[4]*x[7] - 5.04*x[1] - 0.035*x[2] - 10*x[3] - 3.36*x[5])
#@NLobjective(m, Min, -(x[4]*x[7]))
#@NLobjective(m, Min, (0.063*x[4]*x[7]^2 - 5.04*x[1]^2 - 0.035*x[2] - 10*x[3] - 3.36*x[5]))

#=
@NLobjective(m, Max, -2.0*x[1])

JuMP.optimize!(m)
# Solve the optimization problem
#using ProfileView
#@profview f(m)
#@profview f(m)

@show objective_value(m)
@show objective_bound(m)
@show value(x[1])
=#

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
module PluginWrapper

using JavaCall

const Class = @jimport java.lang.Class
const DefaultDataContextPlugin = @jimport org.orekit.compiler.plugin.DefaultDataContextPlugin
const JString = @jimport java.lang.String
const JavacTask = @jimport com.sun.source.util.JavacTask
const Object = @jimport java.lang.Object
const Plugin = @jimport com.sun.source.util.Plugin
const TaskEvent = @jimport com.sun.source.util.TaskEvent

export DefaultDataContextPlugin
export auto_start
export equals
export finished
export get_class
export get_name
export hash_code
export init
export notify
export notify_all
export started
export to_string
export wait

include("default_data_context_plugin.jl")

end


text:
struct InpContent{dim, TF, N, TI}
    node_coords::Vector{NTuple{dim,TF}}
    celltype::String
    cells::Vector{NTuple{N,TI}}
    nodesets::Dict{String,Vector{TI}}
    cellsets::Dict{String,Vector{TI}}
    E::TF
    ν::TF
    density::TF
    nodedbcs::Dict{String, Vector{Tuple{TI,TF}}}
    cloads::Dict{Int, Vector{TF}}
    facesets::Dict{String, Vector{Tuple{TI,TI}}}
    dloads::Dict{String, TF}
end

const stopping_pattern = r"^\*[^\*]"

function extract_inp(filepath_with_ext)
    file = open(filepath_with_ext, "r")

    local node_coords
    local celltype, cells, offset
    nodesets = Dict{String,Vector{Int}}()
    cellsets = Dict{String,Vector{Int}}()
    local E, mu
    nodedbcs = Dict{String, Vector{Tuple{Int,Float64}}}()
    cloads = Dict{Int, Vector{Float64}}()
    facesets = Dict{String, Vector{Tuple{Int,Int}}}()
    dloads = Dict{String, Float64}()
    density = 0. # Should extract from the file

    node_heading_pattern = r"\*Node\s*,\s*NSET\s*=\s*([^,]*)"
    cell_heading_pattern = r"\*Element\s*,\s*TYPE\s*=\s*([^,]*)\s*,\s*ELSET\s*=\s*([^,]*)"
    nodeset_heading_pattern = r"\*NSET\s*,\s*NSET\s*=\s*([^,]*)"
    cellset_heading_pattern = r"\*ELSET\s*,\s*ELSET\s*=\s*([^,]*)"
    material_heading_pattern = r"\*MATERIAL\s*,\s*NAME\s*=\s*([^\s]*)"
    boundary_heading_pattern = r"\*BOUNDARY"
    cload_heading_pattern = r"\*CLOAD"
    dload_heading_pattern = r"\*DLOAD"

    line = readline(file)
    while !eof(file)
        m = match(node_heading_pattern, line)
        if m != nothing && m[1] == "Nall"
            node_coords, line = extract_nodes(file)
            dim = length(node_coords[1])
            continue
        end
        m = match(cell_heading_pattern, line)
        if m != nothing
            celltype = String(m[1])
            cellsetname = String(m[2])
            cells, offset, line = extract_cells(file)
            cellsets[cellsetname] = collect(1:length(cells))
            continue
        end
        m = match(nodeset_heading_pattern, line)
        if m != nothing
            nodesetname = String(m[1])
            line = extract_set!(nodesets, nodesetname, file)
            continue
        end
        m = match(cellset_heading_pattern, line)
        if m != nothing
            cellsetname = String(m[1])
            line = extract_set!(cellsets, cellsetname, file, offset)
            continue
        end
        m = match(material_heading_pattern, line)
        if m != nothing
            material_name = String(m[1])
            E, mu, line = extract_material(file)
            continue
        end
        m = match(boundary_heading_pattern, line)
        if m != nothing
            line = extract_nodedbcs!(nodedbcs, file)
            continue
        end
        m = match(cload_heading_pattern, line)
        if m != nothing
            line = extract_cload!(cloads, file, Val{dim})
            continue
        end
        m = match(dload_heading_pattern, line)
        if m != nothing
            line = extract_dload!(dloads, facesets, file, Val{dim}, offset)
            continue
        end
        line = readline(file)
    end

    close(file)

    return InpContent(node_coords, celltype, cells, nodesets, cellsets, E, mu, density, nodedbcs, cloads, facesets, dloads)
end


text:
@testset "counts" begin
    for _ in 1:100
        a1, a2 = rand(1:10), rand(1:10)
        b1, b2 = a1 + rand(0:10), a2 + rand(0:10)
        c1 = rand(1:20, b1 - a1 + 1)
        c2 = rand(1:20, b2 - a2 + 1)
        cap = rand((a1+a2):(b1+b2))
        D = Dict{Int,Int}()
        for (i1, c1) in zip(a1:b1, c1)
            for (i2, c2) in zip(a2:b2, c2)
                k = i1 + i2
                if k ≤ cap
                    D[k] = get(D, k, 0) + c1 * c2
                end
            end
        end
        C = count_combinations(cap, Counts(a1:b1, c1), Counts(a2:b2, c2))
        for (i, c) in pairs(C)
            @test get(D, i, 0) == c
        end
        i1, i2 = extrema(keys(D))
        @test firstindex(C) ≤ i1
        @test i2 ≤ lastindex(C) ≤ cap
    end
end

"Simple function to test iterator."
function naive_capped_indices(cap, I)
    [Tuple(ι) for ι in CartesianIndices(map(i -> 1:i, I)) if sum(Tuple(ι)) ≤ cap]
end

@testset "capped cartesian indices" begin
    for _ in 1:100
        I = Tuple(rand(1:5, rand(1:5)))
        cap = length(I) + rand(0:5)
        iter = CappedCartesianIndices(cap, I)
        naive = naive_capped_indices(cap, I)
        @test length(naive) == @inferred length(iter)
        @test collect(iter) == naive
        @test @inferred(largest_indices(iter)) ==
            Tuple(maximum(mapreduce(x -> [x...], hcat, naive); dims = 2))
    end
    @test_throws ArgumentError CappedCartesianIndices(1, (3, 3,)) # sum too low
    @test_throws ArgumentError CappedCartesianIndices(5, (0, 3,)) # negative
    T = typeof(CappedCartesianIndices(5, (3, 3,)))
    @test Base.IteratorSize(T) ≡ Base.HasLength()
    @test Base.IteratorEltype(T) ≡ Base.HasEltype()
    @test Base.eltype(T) ≡ Tuple{Int,Int}
end