text:
#
# """
# Take one step of iLQR algorithm (non-allocating)
# """
# function step!(solver::iLQRSolver, J)
#     Z = solver.Z
#     state_diff_jacobian!(solver.G, solver.model, Z)
#     discrete_jacobian!(solver.∇F, solver.model, Z)
#     cost_expansion!(solver.Q, solver.G, solver.obj, solver.model, solver.Z)
#     ΔV = backwardpass!(solver)
#     forwardpass!(solver, ΔV, J)
# end
#
#
# """
# $(SIGNATURES)
# Calculates the optimal feedback gains K,d as well as the 2nd Order approximation of the
# Cost-to-Go, using a backward Riccati-style recursion. (non-allocating)
# """
# function backwardpass!(solver::StaticiLQRSolver{T,QUAD}) where {T,QUAD<:QuadratureRule}
#     n,m,N = size(solver)
#
#     # Objective
#     obj = solver.obj
#     model = solver.model
#
#     # Extract variables
#     Z = solver.Z; K = solver.K; d = solver.d;
#     G = solver.G
#     S = solver.S
#     Q = solver.Q
#
#     # Terminal cost-to-go
#     S.xx[N] = Q.xx[N]
#     S.x[N] = Q.x[N]
#
#     # Initialize expecte change in cost-to-go
#     ΔV = @SVector zeros(2)
#
#
#     k = N-1
#     while k > 0
#         ix = Z[k]._x
#         iu = Z[k]._u
#
#         # fdx = G[k+1]'solver.∇F[k][ix,ix]*G[k]
#         # fdu = G[k+1]'solver.∇F[k][ix,iu]
#         fdx,fdu = dynamics_expansion(solver.∇F[k], G[k], G[k+1], model, Z[k])
#         # fdx, fdu = dynamics_expansion(QUAD, model, Z[k])
#
#         Qx =  Q.x[k] + fdx'S.x[k+1]
#         Qu =  Q.u[k] + fdu'S.x[k+1]
#         Qxx = Q.xx[k] + fdx'S.xx[k+1]*fdx
#         Quu = Q.uu[k] + fdu'S.xx[k+1]*fdu
#         Qux = Q.ux[k] + fdu'S.xx[k+1]*fdx
#
#         if solver.opts.bp_reg_type == :state
#             Quu_reg = Quu + solver.ρ[1]*fdu'fdu
#             Qux_reg = Qux + solver.ρ[1]*fdu'fdx
#         elseif solver.opts.bp_reg_type == :control
#             Quu_reg = Quu + solver.ρ[1]*I
#             Qux_reg = Qux
#         end
#
#         # Regularization
#         if solver.opts.bp_reg
#             vals = eigvals(Hermitian(Quu_reg))
#             if minimum(vals) <= 0
#                 @warn "Backward pass regularized"
#                 regularization_update!(solver, :increase)
#                 k = N-1
#                 ΔV = @SVector zeros(2)
#                 continue
#             end
#         end
#
#         # Compute gains
#         K[k] = -(Quu_reg\Qux_reg)
#         d[k] = -(Quu_reg\Qu)
#
#         # Calculate cost-to-go (using unregularized Quu and Qux)
#         S.x[k]  =  Qx + K[k]'*Quu*d[k] + K[k]'* Qu + Qux'd[k]
#         S.xx[k] = Qxx + K[k]'*Quu*K[k] + K[k]'*Qux + Qux'K[k]
#         S.xx[k] = 0.5*(S.xx[k] + S.xx[k]')
#
#         # calculated change is cost-to-go over entire trajectory
#         ΔV += @SVector [d[k]'*Qu, 0.5*d[k]'*Quu*d[k]]
#
#         k -= 1
#     end
#
#     regularization_update!(solver, :decrease)
#
#     return ΔV
#
# end
#

#
# """
# $(SIGNATURES)
# Simulate forward the system with the optimal feedback gains from the iLQR backward pass.
# (non-allocating)
# """
# function rollout!(solver::StaticiLQRSolver{T,Q}, α) where {T,Q}
#     Z = solver.Z; Z̄ = solver.Z̄
#     K = solver.K; d = solver.d;
#
#     Z̄[1].z = [solver.x0; control(Z[1])]
#
#     temp = 0.0
#
#
#     for k = 1:solver.N-1
#         δx = state_diff(solver.model, state(Z̄[k]), state(Z[k]))
#         ū = control(Z[k]) + K[k]*δx + α*d[k]
#         set_control!(Z̄[k], ū)
#
#         # Z̄[k].z = [state(Z̄[k]); control(Z[k]) + δu]
#         Z̄[k+1].z = [discrete_dynamics(Q, solver.model, Z̄[k]);
#             control(Z[k+1])]
#
#         temp = norm(Z̄[k+1].z)
#         if temp > solver.opts.max_state_value
#             return false
#         end
#     end
#     return true
# end


text:
# io functions for neid solar spectra

#make_manifest copied from EchelleInstruments/src/files.jl and modified to walk through neid file structure
function get_filenames_NEID(data_path::String, r::Regex; verbose::Bool = true)
   if !isdir(data_path)
      @error "Can't access data directory ", data_path, "."
   end
   if verbose  println("# Creating manifest of files to process.")    end

   #dir_filelist = readdir(data_path,join=true)
   dir_filelist = Vector{String}()
   for (root, dirs, files) in walkdir(data_path, follow_symlinks=true)
      for file in files
          append!(dir_filelist,[joinpath(root, file)]) # path to files
      end
   end

   idx_spectra = map(fn->occursin(r, last(split(fn,'/')) ),dir_filelist)
   #spectra_filelist = dir_filelist[idx_spectra]

   df_files = DataFrame(:Filename => dir_filelist[idx_spectra])

   if size(df_files,1) < 1
      @error("Did not find any files in ", data_path ,".")
   end

   if match(r"(\d{8})[T_](\d{6})", last(split(df_files.Filename[1],'/')) ) != nothing
      results = map(fn->match(r"(\d{8})[T_](\d{6})", last(split(fn,'/')) ), df_files.Filename)
      @assert all(map(r->length(r.captures),results).==2)
      df_files.date_str = map(r->String(r.captures[1]),results)
      df_files.time_str = map(r->String(r.captures[2]),results)
   end

   return df_files
end

"""Create Dataframe containing filenames and key data for all files neid*.fits in directory"""
function make_manifest_NEID(data_path::String ; max_spectra_to_use::Int = 1000 )
    df_filenames = get_filenames_NEID(data_path,r"^neidL[1,2]_\d+[T\.]\d+\.fits$")
    df_files = DataFrame(NEID.read_metadata(df_filenames.Filename[1]))
    keys = propertynames(df_files)
    allowmissing!(df_files, keys[map(k->k∉[:Filename, :bjd, :target, :airmass],keys)] )
    #println("df_files = ", df_files)

    #if length(spectra_filelist) >= 2
    if length(df_filenames.Filename) >= 2
        max_idx = min(max_spectra_to_use,length(df_filenames.Filename))
        map(fn->NEID.add_metadata_from_fits!(df_files,fn),df_filenames.Filename[2:max_idx])
    end
    df_files

end

ext: jl

max_stars_repo_name: alexander-wise/RvLineList

max_issues_repo_name: alexander-wise/RvLineList

max_forks_repo_name: alexander-wise/RvLineList


text:
module MakieThemes

include("GGThemr/GGThemr.jl")

using .GGThemr
export ggthemr, show_ggthemr, ggthemr_colorthemes

end # module


text:
module BenchGEMM

using ArgCheck
using BenchmarkTools
using LinearAlgebra
using Referenceables: referenceable
using Transducers
using Transducers: @simd_if, @next, complete, maybe_usesimd, BottomRF, SideEffect

_size(x) = size.(axes(x), 1)
_size(x, i) = size(axes(x)[i], 1)

@inline function ij_ik_kj_foldable(C, A, B)
    @argcheck _size(C) === (_size(A, 1), _size(B, 2))
    @argcheck _size(A, 2) === _size(B, 1)
    @argcheck !any(Base.has_offset_axes, (A, B, C))
    return AdHocFoldable() do rf, acc, _
        Base.@_inline_meta
        for j in 1:size(C, 2), k in 1:size(A, 2)
            b = @inbounds B[k, j]
            @simd_if rf for i in 1:size(A, 1)
                c = @inbounds C[i, j]
                a = @inbounds A[i, k]
                acc = @next(rf, acc, (c, a, b))
            end
        end
        return complete(rf, acc)
    end
end

function xfmul!(C, A, B, simd=Val(false))
    fill!(C, 0)

    CAB = ij_ik_kj_foldable(referenceable(C), A, B)

    # Passing `simd::Val` as a keyword argument introduces an
    # `invoke`, even if I add `@inline` annotation to `foreach`...
    #=
    foreach(CAB; simd=simd) do (c, a, b)
        c[] = muladd(a, b, c[])
    end
    =#

    rf = SideEffect() do (c, a, b)
        c[] = muladd(a, b, c[])
        return  # for type stability
    end
    transduce(maybe_usesimd(BottomRF(rf), simd), nothing, CAB)

    return C
end

macro _simd_if(simd, loop)
    ex = quote
        if $simd === Val(:ivdep)
            $Base.@simd ivdep $loop
        elseif $simd === Val(true)
            $Base.@simd $loop
        else
            $loop
        end
    end
    return esc(ex)
end

function manmul!(C, A, B, simd=Val(true))
    @argcheck size.((C, C, A), (1, 2, 2)) === size.((A, B, B), (1, 2, 1))
    @argcheck !any(Base.has_offset_axes, (A, B, C))
    fill!(C, 0)
    for j in 1:size(C, 2), k in 1:size(A, 2)
        b = @inbounds B[k, j]
        aofs = LinearIndices(A)[1, k] - 1
        cofs = LinearIndices(C)[1, j] - 1
        @_simd_if simd for i in 1:size(A, 1)
            c = @inbounds C[i + cofs]
            a = @inbounds A[i + aofs]
            @inbounds C[i + cofs] = muladd(a, b, c)
        end
    end
    C
end

let
    A = randn(3, 3)
    B = randn(3, 3)
    desired = A * B
    C = zero(desired)

    xfmul!(C, A, B)
    @assert maximum(abs.(C .- desired)) < 1e-8
end

suite = BenchmarkGroup()

let sub = suite["mul"] = BenchmarkGroup()
    for impl in ["xf", "man"]
        sub[impl] = BenchmarkGroup()
        sub[impl]["false"] = BenchmarkGroup()
        sub[impl]["true"] = BenchmarkGroup()
        sub[impl]["ivdep"] = BenchmarkGroup()
    end

    sub["linalg"] = BenchmarkGroup()
end

for n in [8, 32, 256]
    A = randn(n, n)
    B = randn(n, n)
    desired = A * B
    C = zero(desired)

    sub = suite["mul"]
    for (impl, f) in [("xf", xfmul!), ("man", manmul!)]
        for simd in [false, true, :ivdep]
            sub[impl][string(simd)][n] = @benchmarkable(
                $f($C, $A, $B, $(Val(simd)))
            )
        end
    end
    sub["linalg"][n] = @benchmarkable mul!($C, $A, $B)
end

function fusedxfmul!(C, A, B1, B2)
    fill!(C, 0)

    B = Broadcast.instantiate(Broadcast.broadcasted(tuple, B1, B2))
    CAB = ij_ik_kj_foldable(referenceable(C), A, B)

    foreach(CAB; simd=Val(:ivdep)) do (c, a, (b1, b2))
        c0 = c[]
        c1 = muladd(a, b1, c0)
        c2 = a^2 * b2
        c[] = c1 + c2
        return  # for type stability
    end

    return C
end

function nonfusedmul!(C, A1, A2, B1, B2)
    Ⅰ = one(eltype(C))
    ∅ = zero(eltype(C))
    A2 .= A1 .^ 2
    BLAS.gemm!('N', 'N', Ⅰ, A1, B1, ∅, C)
    BLAS.gemm!('N', 'N', Ⅰ, A2, B2, Ⅰ, C)
    return C
end

let sub = suite["fusedmul"] = BenchmarkGroup(),
    n = 1024

    sub["xf"] = BenchmarkGroup()
    sub["blas"] = BenchmarkGroup()

    for m in [2, 8, 16, 32]
        A1 = randn(n, n)
        A2 = A1 .^ 2
        B1 = randn(n, m)
        B2 = randn(n, m)
        C = zero(A1 * B2)

        @assert isapprox(
            copy(fusedxfmul!(C, A1, B1, B2)),
            copy(nonfusedmul!(C, A1, A2, B1, B2)),
        )

        sub["xf"][m] = @benchmarkable fusedxfmul!($C, $A1, $B1, $B2)
        sub["blas"][m] = @benchmarkable nonfusedmul!($C, $A1, $A2, $B1, $B2)
    end
end

end  # module
BenchGEMM.suite


text:
############################################################
# judiPDE ####################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiPDE, judiPDEexception, judiPDEadjoint, judiPDEadjointException, subsample

############################################################

# Type for linear operator representing  Pr*A(m)^-1*Ps,
# i.e. it includes source and receiver projections

# Forward
struct judiPDE{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    model::Model
    geometry::Geometry
    options::Options
    fop::Function              # forward
    fop_T::Union{Function, Nothing}  # transpose
end

mutable struct judiPDEexception <: Exception
    msg :: String
end

# Adjoint
struct judiPDEadjoint{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    model::Model
    geometry::Geometry
    options::Options
    fop::Function              # forward
    fop_T::Union{Function, Nothing}  # transpose
end

mutable struct judiPDEadjointException <: Exception
    msg :: String
end

############################################################
## outer constructors

function judiPDE(name::String,info::Info,model::Model, geometry::Geometry; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
    (DDT == Float32 && RDT == Float32) || throw(judiPDEexception("Domain and range types not supported"))
    if typeof(geometry) == GeometryOOC
        m = sum(geometry.nsamples)
    else
        m = 0
        for j=1:info.nsrc m += length(geometry.xloc[j])*geometry.nt[j] end
    end
    n = info.n * sum(info.nt)
    if info.nsrc > 1
        srcnum = 1:info.nsrc
    else
        srcnum = 1
    end
    return F = judiPDE{Float32,Float32}(name, m, n, info, model, geometry, options,
                              args -> time_modeling(model, args[1], args[2], args[3], args[4], args[5], srcnum, 'F', 1, options),
                              args_T -> time_modeling(model, args_T[1], args_T[2], args_T[3], args_T[4], args_T[5], srcnum, 'F', -1, options)
                              )
end

function judiPDEadjoint(name::String,info::Info,model::Model, geometry::Geometry; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
    (DDT == Float32 && RDT == Float32) || throw(judiPDEadjointException("Domain and range types not supported"))
    if typeof(geometry) == GeometryOOC
        m = sum(geometry.nsamples)
    else
        m = 0
        for j=1:info.nsrc m += length(geometry.xloc[j])*geometry.nt[j] end
    end
    n = info.n * sum(info.nt)
    if info.nsrc > 1
        srcnum = 1:info.nsrc
    else
        srcnum = 1
    end
    return F = judiPDEadjoint{Float32,Float32}(name, m, n, info, model, geometry, options,
                                args_T -> time_modeling(model, args_T[1], args_T[2], args_T[3], args_T[4], args_T[5], srcnum, 'F', -1, options),
                                args -> time_modeling(model, args[1], args[2], args[3], args[4], args[5], srcnum, 'F', 1, options),
                                )
end

############################################################
## overloaded Base functions

# conj(judiPDE)
conj(A::judiPDE{DDT,RDT}) where {DDT,RDT} =
    judiPDE{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        A.fop,
        A.fop_T
        )

# transpose(jo)
transpose(A::judiPDE{DDT,RDT}) where {DDT,RDT} =
    judiPDEadjoint{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
        A.fop_T,
        A.fop
        )

# adjoint(jo)
adjoint(A::judiPDE{DDT,RDT}) where {DDT,RDT} =
	judiPDEadjoint{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
		A.fop_T,
		A.fop
		)

# conj(jo)
conj(A::judiPDEadjoint{DDT,RDT}) where {DDT,RDT}=
    judiPDEadjoint{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        A.fop,
        A.fop_T
        )

# transpose(jo)
transpose(A::judiPDEadjoint{DDT,RDT}) where {DDT,RDT} =
    judiPDE{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
        A.fop_T,
        A.fop
        )
# adjoint(jo)
adjoint(A::judiPDEadjoint{DDT,RDT}) where {DDT,RDT} =
    judiPDE{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
        A.fop_T,
        A.fop
        )

############################################################
## overloaded Base *(...judi...)

# *(judiPDE,judiRHS)
function *(A::judiPDE{ADDT,ARDT},v::judiRHS{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(v,1) || throw(judiPDEexception("shape mismatch"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiRHS):",A.name,typeof(A),vDT]," / "))
	args = (v.geometry,v.data,A.geometry,nothing,nothing)
    V = A.fop(args)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

function *(A::judiPDEadjoint{ADDT,ARDT},v::judiRHS{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(v,1) || throw(judiPDEadjointException("shape mismatch"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiRHS):",A.name,typeof(A),vDT]," / "))
	args = (A.geometry,nothing,v.geometry,v.data,nothing)
    V = A.fop(args)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(judiPDE,judiWavefield)
function *(A::judiPDE{ADDT,ARDT},v::judiWavefield{vDT}) where {ADDT,ARDT,vDT}
	A.n == size(v,1) || throw(judiPDEexception("shape mismatch"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiWavefield):",A.name,typeof(A),vDT]," / "))
	args = (nothing,v.data,A.geometry,nothing,nothing)
	V = A.fop(args)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

function *(A::judiPDEadjoint{ADDT,ARDT},v::judiWavefield{vDT}) where {ADDT,ARDT,vDT}
	A.n == size(v,1) || throw(judiPDEadjointException("shape mismatch"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiWavefield):",A.name,typeof(A),vDT]," / "))
	args = (A.geometry,nothing,nothing,v.data,nothing)
	V = A.fop(args)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# *(judiPDE,judiProjection)
function *(A::judiPDE{CDT,ARDT},B::judiProjection{BDDT,CDT}) where {ARDT,BDDT,CDT}
    A.n == size(B,1) || throw(judiPDEexception("shape mismatch"))
    return judiModeling(A.info,A.model,B.geometry,A.geometry;options=A.options,DDT=CDT,RDT=ARDT)
end

function *(A::judiPDEadjoint{CDT,ARDT},B::judiProjection{BDDT,CDT}) where {ARDT,BDDT,CDT}
    A.n == size(B,1) || throw(judiPDEadjointException("shape mismatch"))
    return adjoint(judiModeling(A.info,A.model,A.geometry,B.geometry,options=A.options,DDT=CDT,RDT=ARDT))
end

# *(num,judiPDE)
function *(a::Number,A::judiPDE{ADDT,ARDT}) where {ADDT,ARDT}
    return judiPDE{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1 -> jo_convert(ARDT,a*A.fop(v1),false),
        v2 -> jo_convert(ADDT,a*A.fop_T(v2),false)
        )
end

function *(a::Number,A::judiPDEadjoint{ADDT,ARDT}) where {ADDT,ARDT}
    return judiPDEadjoint{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1 -> jo_convert(ARDT,a*A.fop(v1),false),
        v2 -> jo_convert(ADDT,a*A.fop_T(v2),false)
        )
end

############################################################
## overloaded Basees +(...judiPDE...), -(...judiPDE...)

# +(judiPDE,num)
function +(A::judiPDE{ADDT,ARDT},b::Number) where {ADDT,ARDT}
    return judiPDE{ADDT,ARDT}("("*A.name*"+N)",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1 -> A.fop(v1)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
        v2 -> A.fop_T(v2)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
        )
end

function +(A::judiPDEadjoint{ADDT,ARDT},b::Number) where {ADDT,ARDT}
    return judiPDEadjoint{ADDT,ARDT}("("*A.name*"+N)",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1->A.fop(v1)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
        v2->A.fop_T(v2)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
        )
end

# -(judiPDE,num)
function -(A::judiPDE{ADDT,ARDT},b::Number) where {ADDT,ARDT}
    return judiPDE{ADDT,ARDT}("("*A.name*"-N)",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1 -> A.fop(v1)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
        v2 -> A.fop_T(v2)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
        )
end

function -(A::judiPDEadjoint{ADDT,ARDT},b::Number) where {ADDT,ARDT}
    return judiPDEadjoint{ADDT,ARDT}("("*A.name*"-N)",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1->A.fop(v1)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
        v2->A.fop_T(v2)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
        )
end

# -(judiPDE)
-(A::judiPDE{DDT,RDT}) where {DDT,RDT} =
    judiPDE{DDT,RDT}("(-"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1->-A.fop(v1),
        v2->-A.fop_T(v2)
        )

-(A::judiPDEadjoint{DDT,RDT}) where {DDT,RDT} =
    judiPDEadjoint{DDT,RDT}("(-"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1->-A.fop(v1),
        v2->-A.fop_T(v2)
        )

############################################################
## Additional overloaded functions

# Subsample Modeling operator
function subsample(F::judiPDE{ADDT,ARDT}, srcnum) where {ADDT,ARDT}
    geometry = subsample(F.geometry,srcnum)     # Geometry of subsampled data container
    info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
    return judiPDE(F.name, info, F.model, geometry; options=F.options, DDT=ADDT, RDT=ARDT)
end

function subsample(F::judiPDEadjoint{ADDT,ARDT}, srcnum) where {ADDT,ARDT}
    geometry = subsample(F.geometry,srcnum)     # Geometry of subsampled data container
    info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
    return judiPDEadjoint(F.name, info, F.model, geometry; options=F.options, DDT=ADDT, RDT=ARDT)
end

getindex(F::judiPDE,a) = subsample(F,a)
getindex(F::judiPDEadjoint,a) = subsample(F,a)


text:
module ModelingToolkit

using DiffEqBase
using StaticArrays, LinearAlgebra

using MacroTools
import MacroTools: splitdef, combinedef

abstract type Expression <: Number end
abstract type AbstractOperation <: Expression end
abstract type AbstractComponent <: Expression end
abstract type AbstractSystem end

include("variables.jl")

Base.promote_rule(::Type{T},::Type{T2}) where {T<:Number,T2<:Expression} = Expression
Base.zero(::Type{<:Expression}) = Constant(0)
Base.one(::Type{<:Expression}) = Constant(1)
Base.convert(::Type{Variable},x::Int64) = Constant(x)

function caclulate_jacobian end

@enum FunctionVersions ArrayFunction=1 SArrayFunction=2

include("operations.jl")
include("differentials.jl")
include("equations.jl")
include("systems/diffeqs/diffeqsystem.jl")
include("systems/diffeqs/first_order_transform.jl")
include("systems/nonlinear/nonlinear_system.jl")
include("function_registration.jl")
include("simplify.jl")
include("utils.jl")

export Operation, Expression, AbstractComponent, AbstractDomain
export @register
end # module


text:
export fast_pair!, full_pair!

"only find mates for primary reads (flag & 0x900 == 0)"
function fast_pair!(bam::Bam)
    namedict = Dict{String, Read}()
    for r in bam.reads
        r.flag & 0x900 == 0 || continue
        if r.qname in keys(namedict)
            r.mate = namedict[r.qname]
            namedict[r.qname].mate = r
            delete!(namedict, r.qname)
        else
            namedict[r.qname] = r
        end
    end
    bam
end

"mates are primary reads"
function full_pair!(bam::Bam)
    error("TODO")
end


text:
using Conda, PyCall, Logging
try
	Conda.add("numpy")
	Conda.add("scipy")
	Conda.add("cython")
	Conda.add("pip")
	run(`$(PyCall.pyprogramname) -m pip install --user ttpy`)
	@info "Python TT Toolbox installed"
catch;
	@warn "Python TT Toolbox installation failed; the TTPy submodule will not be usable"
end


text:
export alloc_special, alloc_local

"Allocates on-device memory statically from the specified address space."
@generated function alloc_special(::Val{id}, ::Type{T}, ::Val{as}, ::Val{len}) where {id,T,as,len}
    JuliaContext() do ctx
        eltyp = convert(LLVMType, T, ctx)

        # old versions of GPUArrays invoke _shmem with an integer id; make sure those are unique
        if !isa(id, String) || !isa(id, Symbol)
            id = "alloc_special_$id"
        end

        T_ptr = convert(LLVMType, DevicePtr{T,as}, ctx)

        # create a function
        llvm_f, _ = create_function(T_ptr)

        # create the global variable
        mod = LLVM.parent(llvm_f)
        gv_typ = LLVM.ArrayType(eltyp, len)
        gv = GlobalVariable(mod, gv_typ, string(id), convert(Int, as))
        if len > 0
            linkage!(gv, LLVM.API.LLVMInternalLinkage)
            #initializer!(gv, null(gv_typ))
        end

        # by requesting a larger-than-datatype alignment, we might be able to vectorize.
        # TODO: Make the alignment configurable
        alignment!(gv, Base.max(32, Base.datatype_alignment(T)))

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry", ctx)
            position!(builder, entry)

            ptr_with_as = gep!(builder, gv, [ConstantInt(0, ctx),
                                             ConstantInt(0, ctx)])

            val = ptrtoint!(builder, ptr_with_as, T_ptr)
            ret!(builder, val)
        end

        call_function(llvm_f, DevicePtr{T,as})
    end
end

@inline alloc_local(id, T, len) = alloc_special(Val(id), Val(T), Val(AS.Local), Val(len))

@inline @generated function alloc_string(::Val{str}) where str
    JuliaContext() do ctx
        T_pint8_generic = LLVM.PointerType(LLVM.Int8Type(ctx), convert(Int, AS.Generic))
        llvm_f, _ = create_function(LLVM.Int64Type(ctx))
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry", ctx)
            position!(builder, entry)
            str_ptr = globalstring_ptr!(builder, String(str))
            str_ptr_i64 = ptrtoint!(builder, str_ptr, LLVM.Int64Type(ctx))
            ret!(builder, str_ptr_i64)
        end
        call_function(llvm_f, DevicePtr{UInt8,AS.Generic})
    end
end

# TODO: Support various types of len
@inline @generated function memcpy!(dest_ptr::DevicePtr{UInt8,DestAS}, src_ptr::DevicePtr{UInt8,SrcAS}, len::LT) where {DestAS,SrcAS,LT<:Union{Int64,UInt64}}
    JuliaContext() do ctx
        T_nothing = LLVM.VoidType(ctx)
        dest_as = convert(Int, DestAS)
        src_as = convert(Int, SrcAS)
        T_int8 = LLVM.Int8Type(ctx)
        T_int64 = LLVM.Int64Type(ctx)
        T_pint8_dest = LLVM.PointerType(T_int8, dest_as)
        T_pint64_dest = LLVM.PointerType(T_int64, dest_as)
        T_pint8_src = LLVM.PointerType(T_int8, src_as)
        T_pint64_src = LLVM.PointerType(T_int64, src_as)
        T_int1 = LLVM.Int1Type(ctx)

        llvm_f, _ = create_function(T_nothing, [T_int64, T_int64, T_int64])
        mod = LLVM.parent(llvm_f)
        T_intr = LLVM.FunctionType(T_nothing, [T_pint8_dest, T_pint8_src, T_int64, T_int1])
        intr = LLVM.Function(mod, "llvm.memcpy.p$(dest_as)i8.p$(src_as)i8.i64", T_intr)
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry", ctx)
            position!(builder, entry)

            dest_ptr_i64 = inttoptr!(builder, parameters(llvm_f)[1], T_pint64_dest)
            dest_ptr_i8 = bitcast!(builder, dest_ptr_i64, T_pint8_dest)

            src_ptr_i64 = inttoptr!(builder, parameters(llvm_f)[2], T_pint64_src)
            src_ptr_i8 = bitcast!(builder, src_ptr_i64, T_pint8_src)

            call!(builder, intr, [dest_ptr_i8, src_ptr_i8, parameters(llvm_f)[3], ConstantInt(T_int1, 0)])
            ret!(builder)
        end
        call_function(llvm_f, Nothing, Tuple{DevicePtr{UInt8,DestAS},DevicePtr{UInt8,SrcAS},LT}, :((dest_ptr, src_ptr, len)))
    end
end
memcpy!(dest_ptr::DevicePtr{T,DestAS}, src_ptr::DevicePtr{T,SrcAS}, len::Integer) where {T,DestAS,SrcAS} =
    memcpy!(convert(DevicePtr{UInt8,DestAS}, dest_ptr), convert(DevicePtr{UInt8,SrcAS}, src_ptr), UInt64(len))
@inline @generated function memset!(dest_ptr::DevicePtr{UInt8,DestAS}, value::UInt8, len::LT) where {DestAS,LT<:Union{Int64,UInt64}}
    JuliaContext() do ctx
        T_nothing = LLVM.VoidType(ctx)
        dest_as = convert(Int, DestAS)
        T_int8 = LLVM.Int8Type(ctx)
        T_int64 = LLVM.Int64Type(ctx)
        T_pint8_dest = LLVM.PointerType(T_int8, dest_as)
        T_pint64_dest = LLVM.PointerType(T_int64, dest_as)
        T_int1 = LLVM.Int1Type(ctx)

        llvm_f, _ = create_function(T_nothing, [T_int64, T_int8, T_int64])
        mod = LLVM.parent(llvm_f)
        T_intr = LLVM.FunctionType(T_nothing, [T_pint8_dest, T_int8, T_int64, T_int1])
        intr = LLVM.Function(mod, "llvm.memset.p$(dest_as)i8.i64", T_intr)
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry", ctx)
            position!(builder, entry)

            dest_ptr_i64 = inttoptr!(builder, parameters(llvm_f)[1], T_pint64_dest)
            dest_ptr_i8 = bitcast!(builder, dest_ptr_i64, T_pint8_dest)

            call!(builder, intr, [dest_ptr_i8, parameters(llvm_f)[2], parameters(llvm_f)[3], ConstantInt(T_int1, 0)])
            ret!(builder)
        end
        call_function(llvm_f, Nothing, Tuple{DevicePtr{UInt8,DestAS},UInt8,LT}, :((dest_ptr, value, len)))
    end
end
memset!(dest_ptr::DevicePtr{T,DestAS}, value::UInt8, len::Integer) where {T,DestAS} =
    memset!(convert(DevicePtr{UInt8,DestAS}, dest_ptr), value, UInt64(len))


text:
using DelimitedFiles
using LinearAlgebra

@testset "Celerite Tests" begin
    @testset "Celerite Kalman and Raw Log-Likelihoods Agree" begin
        # Test that the raw and Kalman log-likelihood agree
        # The data are a snippet from KIC 011615890
        data = readdlm("test-data.dat")
        test_ts = data[:,1]
        test_fs = data[:,2]
        test_dfs = data[:,3]

        post = Celerite.CeleriteKalmanPosterior(test_ts, test_fs, test_dfs,
                                                1, 1,
                                                210900, 211000,
                                                10.0, 20.0,
                                                0.1, 10.0,
                                                0.5, 1000.0,
                                                1.0, 1000.0,
                                                1.0, 1000.0)

        drw_rms = [sqrt(1000.0)]
        drw_rates = [1.0]
        osc_rms = [sqrt(1432.0)]
        osc_freqs = [18.0]
        osc_Qs = [100.0]

        mu = 210933.0
        nu = 1.0

        p = Celerite.CeleriteKalmanParams(mu, nu, drw_rms, drw_rates, osc_rms, osc_freqs, osc_Qs)
        x = Celerite.to_array(post, p)

        N = size(test_ts,1)

        filt = Celerite.CeleriteKalmanFilter(mu, drw_rms, drw_rates, osc_rms, osc_freqs, osc_Qs)
        ll_filt = Celerite.log_likelihood(filt, test_ts, test_fs, test_dfs)
        ll_params = Celerite.log_likelihood(post, x)
        cov = Celerite.raw_covariance(test_ts, test_dfs, drw_rms, drw_rates, osc_rms, osc_freqs, osc_Qs)

        F = cholesky(cov)
        L = F.L
        logdet = 0.0
        for i in 1:size(cov,1)
            logdet += log(L[i,i])
        end
        ll_raw = -0.5*N*log(2.0*pi) - logdet - 0.5*dot((test_fs .- mu), F \ (test_fs .- mu))

        @test abs(ll_filt - ll_raw[1]) < 1e-10
        @test abs(ll_filt - ll_params) < 1e-10
    end

    @testset "Check that PSD can be run" begin
        data = readdlm("test-data.dat")
        test_ts = data[:,1]
        test_fs = data[:,2]
        test_dfs = data[:,3]

        drw_rms = [sqrt(1000.0)]
        drw_rates = [1.0]
        osc_rms = [sqrt(1432.0)]
        osc_freqs = [18.0]
        osc_Qs = [100.0]

        mu = 210933.0

        N = size(test_ts,1)

        filt = Celerite.CeleriteKalmanFilter(mu, drw_rms, drw_rates, osc_rms, osc_freqs, osc_Qs)

        T = test_ts[end]-test_ts[1]
        dtmin = minimum(diff(test_ts))

        psdfs = exp.(range(-log(T), stop=-log(2.0*dtmin), length=1000))

        psd = Celerite.psd(filt, psdfs)

        @test psd[1] > 0.0
    end
end


text:
"""
A wrapper for the interpolation methods of DataInterpolations.jl.

$(SIGNATURES)

Wraps the methods in such a way that they are callable as `f(u,t)` to
create and return an interpolation of `u` over `t`.
The first argument of the constructor always defines the interpolation method,
all following arguments will be used in the interpolation.

The additional keyword `crop = false` indicates to discard the first and last element of the time series.

# Example

```julia
# Create the wrapper struct
itp_method = InterpolationMethod(QuadraticSpline)
# Create a callable interpolation
itp = itp_method(u,t)
# Return u[2]
itp(t[2])
```
"""
struct InterpolationMethod{T} <: AbstractInterpolationMethod
  itp::T
  args

  function InterpolationMethod(itp, args...)
    return new{typeof(itp)}(itp, args)
  end

end

(x::InterpolationMethod)(u, t) = x.itp(u,t,x.args...)

# TODO Wrap all types
# Wrap the common itps
InterpolationMethod() = InterpolationMethod(LinearInterpolation)

# Taken from DiffEqFlux
# https://github.com/SciML/DiffEqFlux.jl/blob/master/src/collocation.jl
# On 3-11-2021

struct EpanechnikovKernel <: CollocationKernel end
struct UniformKernel <: CollocationKernel end
struct TriangularKernel <: CollocationKernel end
struct QuarticKernel <: CollocationKernel end
struct TriweightKernel <: CollocationKernel end
struct TricubeKernel <: CollocationKernel end
struct GaussianKernel <: CollocationKernel end
struct CosineKernel <: CollocationKernel end
struct LogisticKernel <: CollocationKernel end
struct SigmoidKernel <: CollocationKernel end
struct SilvermanKernel <: CollocationKernel end

function calckernel(::EpanechnikovKernel,t)
    if abs(t) > 1
        return 0
    else
        return 0.75*(1-t^2)
    end
end

function calckernel(::UniformKernel,t)
    if abs(t) > 1
        return 0
    else
        return 0.5
    end
end

function calckernel(::TriangularKernel,t)
    if abs(t) > 1
        return 0
    else
        return (1-abs(t))
    end
end

function calckernel(::QuarticKernel,t)
  if abs(t)>0
    return 0
  else
    return (15*(1-t^2)^2)/16
  end
end

function calckernel(::TriweightKernel,t)
  if abs(t)>0
    return 0
  else
    return (35*(1-t^2)^3)/32
  end
end

function calckernel(::TricubeKernel,t)
  if abs(t)>0
    return 0
  else
    return (70*(1-abs(t)^3)^3)/80
  end
end

function calckernel(::GaussianKernel,t)
  exp(-0.5*t^2)/(sqrt(2*π))
end

function calckernel(::CosineKernel,t)
  if abs(t)>0
    return 0
  else
    return (π*cos(π*t/2))/4
  end
end

function calckernel(::LogisticKernel,t)
  1/(exp(t)+2+exp(-t))
end

function calckernel(::SigmoidKernel,t)
  2/(π*(exp(t)+exp(-t)))
end

function calckernel(::SilvermanKernel,t)
  sin(abs(t)/2+π/4)*0.5*exp(-abs(t)/sqrt(2))
end

function construct_t1(t,tpoints)
    hcat(ones(eltype(tpoints),length(tpoints)),tpoints.-t)
end

function construct_t2(t,tpoints)
  hcat(ones(eltype(tpoints),length(tpoints)),tpoints.-t,(tpoints.-t).^2)
end

function construct_w(t,tpoints,h,kernel)
    W = @. calckernel((kernel,),(tpoints-t)/h)/h
    Diagonal(W)
end

"""
$(SIGNATURES)

Unified interface for collocation techniques. The input can either be
a `CollocationKernel` (see list below) or a wrapped `InterpolationMethod` from
[DataInterpolations.jl](https://github.com/PumasAI/DataInterpolations.jl).

Computes a non-parametrically smoothed estimate of `u'` and `u`
given the `data`, where each column is a snapshot of the timeseries at
`tpoints[i]`.

# Examples
```julia
u′,u,t = collocate_data(data,tpoints,kernel=SigmoidKernel())
u′,u,t = collocate_data(data,tpoints,tpoints_sample,interp,args...)
u′,u,t = collocate_data(data,tpoints,interp)
```

# Collocation Kernels
See [this paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2631937/) for more information.
+ EpanechnikovKernel
+ UniformKernel
+ TriangularKernel
+ QuarticKernel
+ TriweightKernel
+ TricubeKernel
+ GaussianKernel
+ CosineKernel
+ LogisticKernel
+ SigmoidKernel
+ SilvermanKernel

# Interpolation Methods
See [DataInterpolations.jl](https://github.com/PumasAI/DataInterpolations.jl) for more information.
+ ConstantInterpolation
+ LinearInterpolation
+ QuadraticInterpolation
+ LagrangeInterpolation
+ QuadraticSpline
+ CubicSpline
+ BSplineInterpolation
+ BSplineApprox
+ Curvefit
"""
function collocate_data(data,tpoints,kernel=TriangularKernel(); crop = false, kwargs...)
  _one = oneunit(first(data))
  _zero = zero(first(data))
  e1 = [_one;_zero]
  e2 = [_zero;_one;_zero]
  n = length(tpoints)
  h = (n^(-1/5))*(n^(-3/35))*((log(n))^(-1/16))

  Wd = similar(data, n, size(data,1))
  WT1 = similar(data, n, 2)
  WT2 = similar(data, n, 3)
  x = map(tpoints) do _t
      T1 = construct_t1(_t,tpoints)
      T2 = construct_t2(_t,tpoints)
      W = construct_w(_t,tpoints,h,kernel)
      mul!(Wd,W,data')
      mul!(WT1,W,T1)
      mul!(WT2,W,T2)
      (e2'*((T2'*WT2)\T2'))*Wd,(e1'*((T1'*WT1)\T1'))*Wd
  end
  estimated_derivative = reduce(hcat,transpose.(first.(x)))
  estimated_solution = reduce(hcat,transpose.(last.(x)))
  crop && return estimated_derivative[:,2:end-1], estimated_derivative[:, 2:end-1], tpoints[2:end-1]
  estimated_derivative,estimated_solution, tpoints
end

# Adapted to dispatch on InterpolationMethod
collocate_data(data, tpoints, interp::InterpolationMethod; kwargs...) = collocate_data(data, tpoints, tpoints, interp; kwargs...)

function collocate_data(data::AbstractVector,tpoints::AbstractVector,tpoints_sample::AbstractVector,
                        interp::InterpolationMethod; kwargs...)
  u, du, tpoints = collocate_data(reshape(data, 1, :),tpoints,tpoints_sample,interp; kwargs...)
  return du[1, :], u[1, :] , tpoints
end

# Adapted to dispatch on InterpolationMethod
function collocate_data(data::AbstractMatrix{T},tpoints::AbstractVector{T},
                        tpoints_sample::AbstractVector{T},interp::InterpolationMethod; crop = false, kwargs...) where T

  u = zeros(T,size(data, 1),length(tpoints_sample))
  du = zeros(T,size(data, 1),length(tpoints_sample))
  for d1 in 1:size(data,1)
    interpolation = interp(data[d1,:],tpoints)
    u[d1,:] .= interpolation.(tpoints_sample)
    du[d1,:] .= DataInterpolations.derivative.((interpolation,), tpoints_sample)
  end
  crop && return du[:, 2:end-1], u[:, 2:end-1], tpoints[2:end-1]
  return du, u, tpoints
end


text:
"""
    connectID, maxconnectID = linkclusters!(connectID::Int64,
        maxconnectID::Int64, updateind::Vector{Int64},
        assignment::Vector{Int64})

Force localizations linked by `assignment` to share same the same `connectID`
"""
function linkclusters!(connectID::Vector{Int64}, maxconnectID::Int64,
    updateind::Union{Vector{Int64},Int64}, assignment::Vector{Int64})

    # Initialize each localization being updated as a new cluster.
    nlocalizations = div(length(assignment), 2)
    connectID[updateind] = maxconnectID .+ collect(1:nlocalizations)

    # Loop through the first `nlocalizations` entries of assignment and make
    # the appropriate assignments in `connectID`.
    # NOTE: The cost matrix organization means that the last `nlocalizations`
    #       entries of assignment don't require changes in `connectID`.
    connectID_current = connectID[updateind]
    connectID_copy = deepcopy(connectID_current)
    for nn = 1:nlocalizations
        # If this is a cluster "birth", no action is needed.
        if assignment[nn] > nlocalizations
            continue
        end

        # Determine which ID should be attributed with this assignment.
        newID = minimum([connectID_current[nn]
            connectID_current[assignment[nn]]
            connectID_copy[nn]
            connectID_copy[assignment[nn]]])

        # Update the reassigned values of `connectID`.
        connectID_current[[nn assignment[nn]]] .= newID
    end

    # Store the updated `connectID`s ensuring that they are all greater than
    # the `maxconnectID`.
    connectID[updateind] = connectID_current .+ maxconnectID
    maxconnectID += nlocalizations

    return connectID, maxconnectID
end

"""
    connectID, maxconnectID = linkclusters(
        connectID::Int64, maxconnectID::Int64,
        updateind::Vector{Int64}, assignment::Vector{Int64})

Force localizations linked by `assignment` to share same the same `connectID`
"""
function linkclusters(connectID::Vector{Int64}, maxconnectID::Int64,
    updateind::Union{Vector{Int64},Int64}, assignment::Vector{Int64})

    return linkclusters!(deepcopy(connectID), deepcopy(maxconnectID), updateind, assignment)
end


text:
@testset "Allele frequenies" begin
    @testset "Gene frequencies" begin

        sequences = ["ATCGATCG", "AGGGGG", "CCCCCCCCCCCCCCC", "TTTTTCCCC",
                     "ATCGATCG", "AGGGGG", "ATCGATCG", "ATCGATCG"]

        answer = Dict{String, Float64}("TTTTTCCCC" => 0.125,
                                       "CCCCCCCCCCCCCCC" => 0.125,
                                       "ATCGATCG" => 0.5,
                                       "AGGGGG" => 0.25)

        function test_gene_frequencies(genes, answer, seqtype)
            test_genes = convert(Vector{seqtype}, sequences)
            test_answer = Dict{seqtype, Float64}(convert(seqtype, key) => val for (key, val) in answer)
            @test gene_frequencies(test_genes) == test_answer
            @test GeneticVariation._gene_frequencies(test_genes, seqtype, Base.SizeUnknown()) == test_answer
            @test GeneticVariation._gene_frequencies(test_genes, seqtype, Base.HasShape{1}()) == test_answer
        end

        for n in (2, 4)
            test_gene_frequencies(sequences, answer, BioSequence{DNAAlphabet{n}})
        end
    end
end


text:
@testset "Testset 3 Read-in" begin
    data = joinpath("data", "testset3")
    result = AnymodResult(data)

    @test typeof(result) == AnymodResult
    @test result.scenarios == ["integriert"]
    @test size(result.summarytable) == (10670, 11)
    @test size(result.costtable) == (3058, 11)
    @test size(result.exchangetable) == (3162, 10)
end


text:
include("dubinsmaneuver2d.jl")

function Vertical(qi, qf, rhomin, pitchmax)
    maneuver = DubinsManeuver2D(qi, qf, rhomin, DubinsStruct(0.0, 0.0, 0.0, Inf, ""))

    dx = maneuver.qf[1] - maneuver.qi[1]
    dy = maneuver.qf[2] - maneuver.qi[2]
    D = sqrt(dx^2 + dy^2)

    # Distance normalization
    d = D/maneuver.rhomin

    # Normalize the problem using rotation
    rotationAngle = mod2pi(atan(dy, dx))
    a = mod2pi(maneuver.qi[3] - rotationAngle)
    b = mod2pi(maneuver.qf[3] - rotationAngle)

    # CSC
    pathLSL = _LSL(maneuver)
    pathRSR = _RSR(maneuver)
    pathLSR = _LSR(maneuver, pitchmax)
    pathRSL = _RSL(maneuver, pitchmax)
    _paths = [pathLSR ,pathLSL, pathRSR, pathRSL]

    a(x) = x.length
    sort!(_paths, by=a)

    for p in _paths
        # chech if the turns are too long (do not meet pitch constraint)
        if abs(p.t) < pi && abs(p.q) < pi
            # check the inclination based on pitch constraint
            center_angle = maneuver.qi[3] + ((p.case[1] == 'L') ? p.t : - p.t)
            if ((center_angle < pitchmax[1]) || (center_angle > pitchmax[2]))
                continue
            end
            maneuver.maneuver = p
            break
        end
    end

    if maneuver.maneuver === nothing
        maneuver.maneuver = DubinsStruct(Inf, Inf, Inf, Inf, "XXX")
    end

    return maneuver
end

########## LSL ##########
function _LSL(self)
    theta1 = self.qi[3]
    theta2 = self.qf[3]

    if theta1 <= theta2
        # start/end points
        p1 = self.qi[1:2]
        p2 = self.qf[1:2]

        radius = self.rhomin

        c1, s1 = radius * cos(theta1), radius * sin(theta1)
        c2, s2 = radius * cos(theta2), radius * sin(theta2)

        # origins of the turns
        o1 = p1 + [ -s1, c1]
        o2 = p2 + [ -s2, c2]

        diff = o2 - o1
        center_distance = sqrt(diff[1]^2 + diff[2]^2)
        centerAngle = atan(diff[2], diff[1])

        t = mod2pi(-theta1 + centerAngle)
        p = center_distance / radius
        q = mod2pi(theta2 - centerAngle)

        if t > pi
            t = 0.0
            q = theta2 - theta1
            turn_end_y = o2[2] - radius * cos(theta1)
            diff_y = turn_end_y - p1[2]
            if abs(theta1) > 1e-5 && (diff_y < 0 == theta1 < 0)
                p = diff_y / sin(theta1) / radius
            else
                t = p = q = Inf
            end
        end
        if q > pi
            t = theta2 - theta1
            q = 0.0
            turn_end_y = o1[2] - radius * cos(theta2)
            diff_y = p2[2] - turn_end_y
            if abs(theta2) > 1e-5 && (diff_y < 0 == theta2 < 0)
                p = diff_y / sin(theta2) / radius
            else
                t = p = q = Inf
            end
        end
    else
        t = p = q = Inf
    end

    length = (t+p+q) * self.rhomin
    case = "LSL"

    return DubinsStruct(t, p, q, length, case)
end

########## RSR ##########
function _RSR(self)
    theta1 = self.qi[3]
    theta2 = self.qf[3]

    if theta2 <= theta1
        # start/end points
        p1 = self.qi[1:2]
        p2 = self.qf[1:2]

        radius = self.rhomin

        c1, s1 = radius * cos(theta1), radius * sin(theta1)
        c2, s2 = radius * cos(theta2), radius * sin(theta2)

        # origins of the turns
        o1 = p1 + [ s1, -c1]
        o2 = p2 + [ s2, -c2]

        diff = o2 - o1
        center_distance = sqrt(diff[1]^2 + diff[2]^2)
        centerAngle = atan(diff[2], diff[1])

        t = mod2pi(theta1 - centerAngle)
        p = center_distance / radius
        q = mod2pi(-theta2 + centerAngle)

        if t > pi
            t = 0.0
            q = -theta2 + theta1
            turn_end_y = o2[2] + radius * cos(theta1)
            diff_y = turn_end_y - p1[2]
            if abs(theta1) > 1e-5 && (diff_y < 0 == theta1 < 0)
                p = diff_y / sin(theta1) / radius
            else
                t = p = q = Inf
            end
        end
        if q > pi
            t = -theta2 + theta1
            q = 0.0
            turn_end_y = o1[2] + radius * cos(theta2)
            diff_y = p2[2] - turn_end_y
            if abs(theta2) > 1e-5 && (diff_y < 0 == theta2 < 0)
                p = diff_y / sin(theta2) / radius
            else
                t = p = q = Inf
            end
        end
    else
        t = p = q = Inf
    end

    length = (t+p+q) * self.rhomin
    case = "RSR"

    return DubinsStruct(t, p, q, length, case)
end

########## LSR ##########
function _LSR(self, pitchmax)
    theta1 = self.qi[3]
    theta2 = self.qf[3]

    # start/end points
    p1 = self.qi[1:2]
    p2 = self.qf[1:2]

    radius = self.rhomin

    c1, s1 = radius * cos(theta1), radius * sin(theta1)
    c2, s2 = radius * cos(theta2), radius * sin(theta2)

    # origins of the turns
    o1 = p1 + [-s1,  c1]
    o2 = p2 + [ s2, -c2]

    diff = o2 - o1
    center_distance = sqrt(diff[1]^2 + diff[2]^2)

    # not constructible
    if center_distance < 2 * radius
        diff[1] = sqrt(4.0 * radius * radius - diff[2] * diff[2])
        alpha = pi/2.0
    else
        alpha = asin(2.0 * radius / center_distance)
    end

    centerAngle = atan(diff[2], diff[1]) + alpha

    if centerAngle < pitchmax[2]
        t = mod2pi(-theta1 + centerAngle)
        p = sqrt(max(0.0, center_distance * center_distance - 4.0 * radius * radius)) / radius
        q = mod2pi(-theta2 + centerAngle)
    else
        centerAngle = pitchmax[2]
        t = mod2pi(-theta1 + centerAngle)
        q = mod2pi(-theta2 + centerAngle)

        # points on boundary between C and S segments
        c1, s1 = radius * cos(centerAngle), radius * sin(centerAngle)
        c2, s2 = radius * cos(centerAngle), radius * sin(centerAngle)
        w1 = o1 - [-s1,  c1]
        w2 = o2 - [ s2, -c2]

        p = (w2[2] - w1[2]) / sin(centerAngle) / radius
    end

    length = (t+p+q) * self.rhomin
    case = "LSR"

    return DubinsStruct(t, p, q, length, case)
end

########## RSL ##########
function _RSL(self, pitchmax)
    theta1 = self.qi[3]
    theta2 = self.qf[3]

    # start/end points
    p1 = self.qi[1:2]
    p2 = self.qf[1:2]

    radius = self.rhomin

    c1, s1 = radius * cos(theta1), radius * sin(theta1)
    c2, s2 = radius * cos(theta2), radius * sin(theta2)

    # origins of the turns
    o1 = p1 + [ s1, -c1]
    o2 = p2 + [-s2,  c2]

    diff = o2 - o1
    center_distance = sqrt(diff[1]^2 + diff[2]^2)

    # not constructible
    if center_distance < 2 * radius
        diff[1] = sqrt(4.0 * radius * radius - diff[2] * diff[2])
        alpha = pi/2.0
    else
        alpha = asin(2.0 * radius / center_distance)
    end

    centerAngle = atan(diff[2], diff[1]) - alpha

    if centerAngle > pitchmax[1]
        t = mod2pi(theta1 - centerAngle)
        p = sqrt(max(0.0, center_distance * center_distance - 4.0 * radius * radius)) / radius
        q = mod2pi(theta2 - centerAngle)
    else
        centerAngle = pitchmax[1]
        t = mod2pi(theta1 - centerAngle)
        q = mod2pi(theta2 - centerAngle)

        # points on boundary between C and S segments
        c1, s1 = radius * cos(centerAngle), radius * sin(centerAngle)
        c2, s2 = radius * cos(centerAngle), radius * sin(centerAngle)
        w1 = o1 - [ s1, -c1]
        w2 = o2 - [-s2,  c2]

        p = (w2[2] - w1[2]) / sin(centerAngle) / radius
    end

    length = (t+p+q) * self.rhomin
    case = "RSL"

    return DubinsStruct(t, p, q, length, case)
end


text:
mutable struct foot_flag_t <: LCMType
    timestamp::Int64
    right::Bool
    left::Bool
end

@lcmtypesetup(foot_flag_t)


text:
include("attributes.jl")
include("subproblem.jl")
include("AbstractProgressiveHedging.jl")


text:
macro shared_fields_stanmodels()
  return esc(:(
    name::AbstractString;              # Name of the Stan program
    model::AbstractString;             # Stan language model program
    n_chains::Vector{Int64};           # Number of chains
    seed::StanBase.RandomSeed;         # Seed section of cmd to run cmdstan
    init::StanBase.Init;               # Init section of cmd to run cmdstan
    output::StanBase.Output;           # Output section of cmd to run cmdstan
    tmpdir::AbstractString;            # Holds all below copied/created files
    output_base::AbstractString;       # Used for naming failes to be created
    exec_path::AbstractString;         # Path to the cmdstan excutable
    data_file::Vector{String};         # Array of data files input to cmdstan
    init_file::Vector{String};         # Array of init files input to cmdstan
    cmds::Vector{Cmd};                 # Array of cmds to be spawned/pipelined
    sample_file::Vector{String};       # Sample file array created by cmdstan
    log_file::Vector{String};          # Log file array created by cmdstan
    diagnostic_file::Vector{String};   # Diagnostic file array created by cmdstan
    summary::Bool;                     # Store cmdstan's summary as a .csv file
    printsummary::Bool;                # Print the summary
    cmdstan_home::AbstractString;      # Directory where cmdstan can be found
  ))
end


text:
N = 2^12
M = 2^9

A = rand(Complex128, N)
B = rand(Complex128, N)

function f1(M::Int, A::Vector{Complex128}, B::Vector{Complex128})
    counter1 = 1
    counter2 = 1
    for i=1:M
        counter1 += 2
        counter2 += 2
        B[counter2] = A[counter1]
    end
end

function f2(M::Int, A::Vector{Complex128}, B::Vector{Complex128})
    ptr1 = pointer(A)
    ptr2 = pointer(B)
    for i=1:M
        ptr1 += 2*16
        ptr2 += 2*16
        unsafe_copy!(ptr2, ptr1, 1)
    end
end

@time f1(M, A, B)
@time f1(M, A, B)
@time f2(M, A, B)
@time f2(M, A, B)


text:
module AutoDrivers

include(Pkg.dir("AutoDrivers", "src", "GaussianMixtureRegressionDrivers", "GaussianMixtureRegressionDrivers.jl"))
# include(Pkg.dir("AutoDrivers", "src", "ConditionalLinearGaussianDrivers", "ConditionalLinearGaussianDrivers.jl"))

end # module


text:
#Regression Constants
export PRAS,INTRSPEEDS,THETAS,PSIS,OWNSPEEDS,ACTIONS,VERTICAL_TAUS,RANGES, nstates
### STATE CUTPOINTS ###
const RANGES = [0.0,25.0,50.0,75.0,100.0,150.0,200.0,300.0,400.0,500.0,510.0,750.0,1000.0,1500.0,2000.0,3000.0,4000.0,5000.0,7000.0,9000.0,11000.0,13000.0,15000.0,17000.0,19000.0,21000.0,25000.0,30000.0,35000.0,40000.0,48000.0,56000.0]
const THETAS = Array(LinRange(-pi,pi,41))
const PSIS   = Array(LinRange(-pi,pi,41))
const OWNSPEEDS = [200.0] #[100.0,200.0,300.0,400.0]  #ft/s
const INTRSPEEDS = [200.0] #[100.0,200.0,300.0,400.0]  #ft/s
const VERTICAL_TAUS   = Array(0:40)
const PRAS    = [0, 1, 2, 3, 4]

const RANGEMAX = maximum(RANGES)/1000.0 #kft
const NSTATES = length(RANGES)*length(THETAS)*length(PSIS)*length(OWNSPEEDS)*length(INTRSPEEDS)
const ACTIONS = deg2rad.([0.0, -1.5, 1.5, -3.0, 3.0])

max_stars_repo_name: cherrywoods/HorizontalCAS

max_issues_repo_name: ahmed-irfan/HorizontalCAS

max_forks_repo_name: ahmed-irfan/HorizontalCAS


text:
# Matern 1/2 isotropic covariance Function

@doc """
# Description
Constructor for the isotropic Matern 1/2 kernel (covariance)

k(x,x') = σ²exp(-d/ℓ), where d=|x-x'|
# Arguments:
* `ll::Float64`: Log of the length scale ℓ
* `lσ::Float64`: Log of the signal standard deviation σ
""" ->
type Mat12Iso <: Isotropic
    ℓ::Float64     # Length scale
    σ2::Float64    # Signal std
    Mat12Iso(ll::Float64, lσ::Float64) = new(exp(ll),exp(2*lσ))
end

function set_params!(mat::Mat12Iso, hyp::Vector{Float64})
    length(hyp) == 2 || throw(ArgumentError("Matern 1/2 covariance function only has two parameters"))
    mat.ℓ, mat.σ2 = exp(hyp[1]), exp(2.0*hyp[2])
end

get_params(mat::Mat12Iso) = Float64[log(mat.ℓ), log(mat.σ2)/2.0]
get_param_names(mat::Mat12Iso) = [:ll, :lσ]
num_params(mat::Mat12Iso) = 2

metric(mat::Mat12Iso) = Euclidean()
cov(mat::Mat12Iso, r::Float64) = mat.σ2*exp(-r/mat.ℓ)

function grad_kern(mat::Mat12Iso, x::Vector{Float64}, y::Vector{Float64})
    r = distance(mat,x,y)
    exp_r = exp(-r/mat.ℓ)

    g1 = mat.σ2*r/mat.ℓ*exp_r  #dK_d(log ℓ)
    g2 = 2.0*mat.σ2*exp_r      #dK_d(log σ)
    return [g1,g2]
end

function grad_stack!(stack::AbstractArray, mat::Mat12Iso, X::Matrix{Float64}, data::IsotropicData)
    nobsv = size(X,2)
    R = distance(mat, X, data)
    exp_R = exp(-R/mat.ℓ)

    for i in 1:nobsv, j in 1:i
        @inbounds stack[i,j,1] = mat.σ2*R[i,j]/mat.ℓ*exp_R[i,j] # dK_dℓ
        @inbounds stack[j,i,1] = stack[i,j,1]
        @inbounds stack[i,j,2] = 2.0 * mat.σ2 * exp_R[i,j]    # dK_dσ
        @inbounds stack[j,i,2] = stack[i,j,2]
    end

    return stack
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later
using CpuId

if cpufeature(:AVX2)
    avxlib = "avx2"
    println("AVX2 support detected, vml_avx2 selected")
else
    avxlib = "avx"
    println("AVX2 support missing, vml_avx selected")
end

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

baseproducts = [("mkl_core", :libmkl_core),
                ("mkl_rt", :libmkl_rt),
                ("mkl_vml_" * avxlib, :libmkl_vml_avx)]
# if Sys.iswindows()
#     push!(baseproducts, ("mkl_intel_thread", :libmkl_intel_thread))
# end

products = [LibraryProduct(prefix, ["lib"*bp[1], bp[1]], bp[2]) for bp in baseproducts]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaPackaging/Yggdrasil/releases/download/MKL-v2019.0.117"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:x86_64)   => ("$bin_prefix/MKL.v2019.0.117.x86_64-linux-gnu.tar.gz", "9a496908c05eccdb331218f4cf32229b024790b30bf7bc0ca0f5587e030d34e6"),
    Linux(:i686)     => ("$bin_prefix/MKL.v2019.0.117.i686-linux-gnu.tar.gz", "ccdce675bf48738f28878bc831231498c7c9560a94f756da8c114f664790ffee"),
    MacOS(:x86_64)   => ("$bin_prefix/MKL.v2019.0.117.x86_64-apple-darwin14.tar.gz", "605da525b16f61837bc35834f3a6ff90609b1d0a5f1606faa25bf45180ced6e9"),
    Windows(:x86_64) => ("$bin_prefix/MKL.v2019.0.117.x86_64-w64-mingw32.tar.gz", "0a9aaac254421fde0f26a95856b595deb8cce85d039050d8167293d049fd716d"),
    Windows(:i686)   => ("$bin_prefix/MKL.v2019.0.117.i686-w64-mingw32.tar.gz", "53cac3a29bbb2acdc383297c256e3b68b3d4712e1b83417226da4971582e80cd"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)

# cleanup for lean install
libfiles = readdir("usr/lib")
idx = .!mapreduce(p -> occursin.(p[1], libfiles), (x,y) -> x.|y, baseproducts)
rm.("usr/lib/" .* libfiles[idx])

rm("usr/downloads/", recursive = true)


text:
module StaticPolynomials

using TypedPolynomials: mergesorted
import TypedPolynomials: AbstractVariable,
                         AbstractMonomial,
                         AbstractTerm,
                         AbstractPolynomial,
                         name,
                         exponents,
                         exponent,
                         variables,
                         coefficient,
                         monomial,
                         subs,
                         terms

import Base: literal_pow, +, -, *, /, ==, isless,
copy, promote_rule, convert

export @polyvar,
       Variable,
       Monomial,
       Term,
       Polynomial,
       name,
       exponents,
       variables,
       coefficient,
       monomial,
       terms,
       degree,
       subs

include("types.jl")
include("conversion.jl")
include("promotion.jl")
include("operators.jl")
include("substitution.jl")

end


text:
using OffsetArrays

#
#
#=========================================================================#
function bc1(meqn::Int, mbc::Int, mx::Int, xlower::Float64, dx::Float64,
             q::OffsetArray{Float64,2}, maux::Int, aux::OffsetArray{Float64,2},
             t::Float64, dt::Float64, mthbc::Array{Int,1})
#=========================================================================#
#
#     # Standard boundary condition choices for claw2
#
#     # At each boundary  k = 1 (left),  2 (right):
#     #   mthbc(k) =  0  for user-supplied BC's (must be inserted!)
#     #            =  1  for zero-order extrapolation
#     #            =  2  for periodic boundary coniditions
#     #            =  3  for solid walls, assuming this can be implemented
#     #                  by reflecting the data about the boundary and then
#     #                  negating the 2'nd component of q.
#     ------------------------------------------------
#
#     # Extend the data from the computational region
#     #      i = 1, 2, ..., mx2
#     # to the virtual cells outside the region, with
#     #      i = 1-ibc  and   i = mx+ibc   for ibc=1,...,mbc
#
#
#
#-------------------------------------------------------
#     # left boundary:
#-------------------------------------------------------
#
    range::UnitRange = 1:meqn
    if mthbc[1] == 0

        # user-specified boundary conditions go here in place of error output
        println("*** ERROR *** mthbc(1)=0 and no BCs specified in bc1")
        exit(-1)

    elseif mthbc[1] == 1
        # zero-order extrapolation:
        for ibc=1:mbc
            q[range, 1-ibc] = q[range, 1]
        end
    elseif mthbc[1] == 2
        # periodic:
        for ibc=1:mbc
            q[range, 1-ibc] = q[range, mx+1-ibc]
        end
    elseif mthbc[1] == 3
        # solid wall (assumes 2'nd component is velocity or momentum in x):
        for ibc=1:mbc
            q[range, 1-ibc] = q[range, ibc]
            q[2,1-ibc] = -q[2,ibc]
        end
        # negate the normal velocity:
    end # if mthbc[1] == ...

#
#-------------------------------------------------------
#     # right boundary:
#-------------------------------------------------------

    if mthbc[2] == 0
        # user-specified boundary conditions go here in place of error output
        println("*** ERROR *** mthbc(2)=0 and no BCs specified in bc2")
        exit(-1)

    elseif mthbc[2] == 1
        # zero-order extrapolation:
        for ibc=1:mbc
            q[range, mx+ibc] = q[range, mx]
        end
    elseif mthbc[2] == 2
        # periodic:
        for ibc=1:mbc
            q[range, mx+ibc] = q[range, ibc]
        end
    elseif mthbc[2] == 3
        # solid wall (assumes 2'nd component is velocity or momentum in x):
        for ibc=1:mbc
            q[range, mx+ibc] = q[range, mx+1-ibc]
            q[2,mx+ibc] = -q[2,mx+1-ibc]
        end
    end # if mthbc[2] == ...

end


text:
using GLFW

mutable struct GlfwOutput <: IFeature
    mVkInstance::vk.VkInstance
    mPhysicalDevice::vk.VkPhysicalDevice
    mDevice::Device

    mTempWindow::GLFW.Window
    mTempSurface::vk.VkSurfaceKHR
    mChainFormat::vk.VkSurfaceFormatKHR

    mPresentIndex::UInt32

    function GlfwOutput()
        this = new()
        GLFW.Init()
        return this
    end
end

@class GlfwOutput

function LavaCore.:destroy(this::GlfwOutput)
    if isdefined(this, :mTempSurface)
        vk.vkDestroySurfaceKHR(getLogicalDevice(this.mDevice), this.mTempSurface, C_NULL)
    end
    if isdefined(this, :mTempWindow)
        GLFW.DestroyWindow(this.mTempWindow)
    end
end

function create(::Type{GlfwOutput})
    return GlfwOutput()
end

function openWindow(this::GlfwOutput, width::UInt32 = UInt32(800), height::UInt32 = UInt32(600), resizable::Bool = false, title::String = "Lava Window")::GlfwWindow
    return GlfwWindow(this.mDevice, this.mChainFormat, width, height, resizable, title)
end

function format(this::GlfwOutput)::vk.VkFormat
    return this.mChainFormat.format
end

function bestFormat(dev::vk.VkPhysicalDevice, surface::vk.VkSurfaceKHR)::vk.VkSurfaceFormatKHR
    formats = VkExt.getSurfaceFormatsKHR(dev, surface)

    if (length(formats) == 1) && (formats[1].format == vk.VK_FORMAT_UNDEFINED)
        return vk.VkSurfaceFormatKHR(
                                        vk.VK_FORMAT_R8G8B8A8_SRGB, #format::VkFormat
                                        vk.VK_COLOR_SPACE_SRGB_NONLINEAR_KHR#colorSpace::VkColorSpaceKHR
                                    )
    end

    for f in formats
        if (f.format == vk.VK_FORMAT_B8G8R8A8_SRGB ||
            f.format == vk.VK_FORMAT_R8G8B8A8_SRGB ||
            f.format == vk.VK_FORMAT_B8G8R8_SRGB ||
            f.format == vk.VK_FORMAT_R8G8B8_SRGB)
            return f
        end
    end

    for f in formats
        if (f.format == vk.VK_FORMAT_B8G8R8A8_UNORM ||
            f.format == vk.VK_FORMAT_R8G8B8A8_UNORM ||
            f.format == vk.VK_FORMAT_B8G8R8_UNORM ||
            f.format == vk.VK_FORMAT_R8G8B8_UNORM)
            return f
        end
    end

    error("No suitable format found!")
end

########## override begin ##########
function LavaCore.:layers(this::GlfwOutput, available::Vector{String})::Vector{String}
    return []
end

function LavaCore.:instanceExtensions(this::GlfwOutput, available::Vector{String})::Vector{String}
    ret = Vector{String}()
    glfwReqExts = GLFW.GetRequiredInstanceExtensions()
    extCount = length(glfwReqExts)

    # TODO: check availability

    append!(ret, glfwReqExts)
    return ret
end

function LavaCore.:deviceExtensions(this::GlfwOutput)::Vector{String}
    return [vk.VK_KHR_SWAPCHAIN_EXTENSION_NAME]
end

function LavaCore.:onInstanceCreated(this::GlfwOutput, vkInstance::vk.VkInstance)
    this.mVkInstance = vkInstance

    GLFW.WindowHint(GLFW.CLIENT_API, GLFW.NO_API)
    GLFW.WindowHint(GLFW.VISIBLE, false)
    #GLFW.WindowHint(GLFW.RESIZABLE, false)

    this.mTempWindow = GLFW.CreateWindow(100, 100, "Vulkan")
    this.mTempSurface = GLFW.CreateWindowSurface(this.mVkInstance, this.mTempWindow)
end

function LavaCore.:onLogicalDeviceCreated(this::GlfwOutput, device::Device)
    this.mDevice = device
end

function LavaCore.:onPhysicalDeviceSelected(this::GlfwOutput, phy::vk.VkPhysicalDevice)
    this.mPhysicalDevice = phy
    this.mChainFormat = bestFormat(phy, this.mTempSurface)
end

function LavaCore.:supportsDevice(this::GlfwOutput, dev::vk.VkPhysicalDevice)::Bool
    families = VkExt.getQueueFamilyProperties(dev)
    for i::UInt32 = 0 : length(families) - 1 #queueFamilyIndex should start from 0
        if VkExt.getSurfaceSupportKHR(dev, i, this.mTempSurface) == vk.VK_TRUE
            return true
        end
    end
    return false
end

function LavaCore.:queueRequests(this::features.GlfwOutput, families::Vector{vk.VkQueueFamilyProperties})
    result = Vector{QueueRequest}()
    for i::UInt32 = 0 : length(families) - 1 #queueFamilyIndex should start from 0
        if VkExt.getSurfaceSupportKHR(this.mPhysicalDevice, i, this.mTempSurface) == vk.VK_TRUE
            this.mPresentIndex = i
            push!(result, createByFamily(QueueRequest, "present", i, 1.0f0))
            break
        end
    end
    if length(result) == 0
        error("Device can't present to this surface.")
    end

    return result
end

# function addPhysicalDeviceFeatures(this::GlfwOutput, outDeviceFeatures::VkExt.VkPhysicalDeviceFeatures)
# end

########## override end ##########

function pollEvent()
    GLFW.PollEvents()
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("findmode_simple.jl")
include("findmode_optim.jl")
include("marginalmode.jl")


text:
"""
  Get the sub-triangulation for an element of the specified order element
  the first 3 values indicate the verticies, later values refer to the nodes
  triangulation must be a 3 x n array of Int32s, so when it gets transposed
  by passing it to C, it becomes a n x 3 array of ints

  This function should only be called for shape types that support
  subtriangulation, otherwise an error is thrown

  Inputs:

    order: the degree of the operator
    shape_type: the integer shape type

  Outputs:

    triangulation: an Int32 array, 2 x numNodesPerElement for CG meshes
                   that have nodes on the vertices, or
                   2 x numNodesPerElement + 3 for DG meshes that do not have
                   nodes on vertices are supported.

  Currently only shape_type == 1 (CG) and shape_type == 2 (DG SBP Omega) are
  supported.
"""
function getTriangulation(order::Integer, shape_type::Integer)
  if shape_type == 1
    if order == 1  # no need to subtriangulate this case
      triangulation = zeros(Int32, 0, 0)
    elseif order == 2
      triangulation = Int32[1 1 4 2 5 6; 4 7 2 5 3 7; 7 6 7 7 7 3]
    elseif order == 3
      triangulation = Int32[1 1 4 5 5 2 6 7 10 12 12 10 9; 4 10 5 11 2 6 7 12 11 7 3 12 10; 10 9 10 10 11 11 11 11 12 3 8 8 8]
    elseif order == 4
      triangulation = Int32[1 1 4 5 5 6 6 2 7 7 15 14 14 8 8 9 17 11 18 13 13 18; 4 13 5 14 6 15 2 7 8 16 16 16 18 17 9 3 3 17 17 18 11 16; 13 12 13 13 14 14 15 15 16 15 14 18 13 16 17 17 10 10 11 11 12 17]
    else
      throw(ErrorException("unsupported triangulation requested"))
    end

  elseif shape_type == 2
    if order == 1
      triangulation = Int32[1 2 2 6 1 5 5; 2 6 3 3 5 4 6; 5 5 6 4 3 3 4]
    elseif order == 2
      triangulation = Int32[1 8 4 1 1 9 2 5 6 1 8 6 1; 4 4 9 9 2 2 3 3 7 6 5 5 8; 8 5 5 4 9 5 5 7 3 3 6 7 6]
    elseif order == 3
      triangulation = Int32[7 1 8 8 8 7 13 7 5 12 8 5 13 13 12 11 12 10 3 11 6;12 5 1 3 4 3 8 4 10 6 13 1 12 9 13 13 11 1 12 1 11;3 9 9 1 3 4 9 8 9 2 7 10 7 10 11 10 6 11 2 2 2]
    #  triangulation = Int32[1 1 1 1 2 6 10 10 5 5 1 2 2 2 3 3 3 9 9 13; 10 11 6 2 12 12 11 12 10 13 5 7 4 3 8 9 1 7 7 13 12; 5 10 11 6 6 11 12 13 13 9 9 12 7 4 4 8 9 4 8 7 7]
    elseif order == 4
      triangulation = Int32[6 11 18 3 12 5 15 12 1 1 4 1 13 3 9 3 16 11 14 2 14 10 18 7 12 11 15 12 7 8 6;
      2 9 7 14 7 15 5 17 16 5 14 9 14 13 10 9 15 10 13 8 10 11 6 11 8 7 1 18 1 2 18;
      17 16 2 8 18 16 1 8 9 16 3 3 4 4 13 13 11 9 10 17 8 8 17 15 11 12 7 17 2 3 2]
    else
      throw(ErrorException("unsupported triangulation requested"))
    end
  else
      throw(ErrorException("unsupported triangulation requested"))
  end  # end if shape_type

  return triangulation
end

function getNodeMaps(order::Integer, shape_type::Integer, numNodesPerElement, dim, isDG::Bool)
# get the mappings between the SBP and Pumi node orderings
# having to do the mapping at all is inelegent to say the least
# store mappings in both directions in case they are needed
# use UInt8s to save cache space during loops

  if !isDG  # CG
    if dim == 2
      if shape_type == 1
        if order == 1
          sbpToPumi = UInt8[1,2,3]
          pumiToSbp = UInt8[1,2,3]
        elseif order == 2
          sbpToPumi = UInt8[1,2,3,4,5,6,7]
          pumiToSbp = UInt8[1,2,3,4,5,6,7]
        elseif order == 3
          # TODO: is it possible to change the apf::FieldShape object to
          # make the sbp and Pumi node orderings the same?
          sbpToPumi = UInt8[1,2,3,4,5,6,7,9,8,12,10,11]
          pumiToSbp= UInt8[1,2,3,4,5,6,7,9,8,11,12,10]
        elseif order == 4
          sbpToPumi = UInt8[1,2,3,4,5,6,7,8,9,12,11,10,17,13,15,14,16,18]
          pumiToSbp = UInt8[1,2,3,4,5,6,7,8,9,12,11,10,14,16,15,17,13,18]
        else
          throw(ErrorException("Unsupported element order $order requestion in getNodeMaps"))
        end  # end if order
      else
        throw(ErrorException("CG only supports shape type 1, requested $(shape_type)"))
      end  # end if shape_type
    else  # 3D
      throw(ErrorException("3 dimensional continuous Galerkin not supported"))
    end  # end if dim

  else # DG
    # the sbp and Pumi nodemaps are the same in this case
    sbpToPumi = collect(UInt8, 1:numNodesPerElement)
    pumiToSbp = collect(UInt8, 1:numNodesPerElement)
  end

  return sbpToPumi, pumiToSbp
end  # end getNodeMaps

function createSubtriangulatedMesh(mesh::AbstractMesh, opts)
# this function is used to figure out how to do subtriangulation/visualization
# for 2D meshes only
# 3D meshes don't subtriangulate

  @assert mesh.dim == 2

  order = mesh.order
  shape_type = mesh.shape_type
  dofpernode = mesh.numDofPerNode

  if !mesh.isDG  #mesh.shape_type == 1
    if order >= 3
      mesh.mnew_ptr, mesh.fnew_ptr, mesh.fnewshape_ptr = createCGSubMesh(mesh)

    else
      mesh.triangulation = zeros(Int32, 0, 0)
      # maybe makes these equal f_ptr and m_ptr for consistency?
      mesh.mnew_ptr = C_NULL
      mesh.fnew_ptr = C_NULL
    end
  else
    # I don't think subtriangulation will work in this case, so create a
    # field with the same degree as a coordinate field on the existing mesh,
    # and interpolate the solution onto it
    # it appears the pumi vtu writing bumps the field up the same as the
    # coordinate field, so we might as well do  it ourselves
    fshape_new = apf.getFieldShape(0, mesh.coord_order, mesh.dim)

    mesh.mnew_ptr = mesh.m_ptr
    # create a new field to store the solution if it does not already exist
    # (reloading mesh after mesh adaptation)
    mesh.fnew_ptr = apf.createPackedField(mesh, "solution_field_interp", dofpernode)
    mesh.fnewshape_ptr = fshape_new

    if (mesh.shape_type == 2) && opts["exact_visualization"]
      mesh.mexact_ptr, mesh.fexact_ptr, mesh.fexactshape_ptr = createDGSubMesh(mesh)
    else
      mesh.mexact_ptr = C_NULL
      mesh.fexact_ptr = C_NULL
      mesh.fexactshape_ptr = C_NULL
    end

    if mesh.shape_type != 2 && opts["exact_visualization"]
      throw(ErrorException("exact visualization only supported for SBP Omega"))
    end
  end  # end if mesh.isDG

  if mesh.myrank == 0
    println("finished createSubtriangulatedMesh")
  end

  return nothing
end

function createCGSubMesh(mesh::PumiMeshCG)
  order = mesh.order
  shape_type = mesh.shape_type
  dofpernode = mesh.numDofPerNode

  mesh.triangulation = getTriangulation(order, shape_type)
  mnew_ptr = apf.createSubMesh(mesh.m_ptr, mesh.triangulation, mesh.elementNodeOffsets, mesh.typeOffsetsPerElement_, mesh.entity_Nptrs)

  fnew_ptr = apf.createPackedField(mnew_ptr, "solution_field", dofpernode)

  fnewshape_ptr = apf.getMeshShapePtr(mnew_ptr)

  return mnew_ptr, fnew_ptr, fnewshape_ptr
end

function createDGSubMesh(mesh::PumiMeshDG)
  order = mesh.order
  shape_type = mesh.shape_type
  dofpernode = mesh.numDofPerNode

  mesh.triangulation = getTriangulation(order, shape_type)
  mnew_ptr = apf.createSubMeshDG(mesh.m_ptr, mesh.mshape_ptr, mesh.triangulation, mesh.elementNodeOffsets, mesh.typeOffsetsPerElement_, mesh.nodemapPumiToSbp, mesh.entity_Nptrs, real(mesh.coords))

  fnew_ptr = apf.createPackedField(mnew_ptr, "solution_field", dofpernode)
  fnewshape_ptr = apf.getMeshShapePtr(mnew_ptr)

  return mnew_ptr, fnew_ptr, fnewshape_ptr
end

function transferFieldToSubmesh(mesh::AbstractMesh, u, mnew_ptr=mesh.mnew_ptr,
                                fnew_ptr=mesh.fnew_ptr)

  if mesh.isInterpolated
    if mesh.shape_type != 3
      apf.transferFieldDG(mesh.m_ptr, mnew_ptr, mesh.triangulation,
                      mesh.elementNodeOffsets, mesh.typeOffsetsPerElement_,
                      mesh.entity_Nptrs, mesh.f_ptr, mesh.interp_op.', fnew_ptr)
    else
      throw(ErrorException("Submesh not supported for SBPGamma"))
    end

      # this is a bit wasteful, because the un-interpolated solution was
      # already saved to mesh.f_ptr
      # interpolateToMesh(mesh, u)
    #end
  else  # CG mesh
    if mesh.order >= 3
      apf.transferField(mesh.m_ptr, mesh.mnew_ptr, mesh.triangulation, mesh.elementNodeOffsets, mesh.typeOffsetsPerElement_, mesh.entity_Nptrs, mesh.f_ptr, mesh.fnew_ptr)
    end
  end

  return nothing
end

"""
  This function gets the xi coordinates of the nodes of a lagrangian entity
  in order verts, edges, faces, regions. 1D, 2D, and 3D entities are supported,
  which is useful because it allows getting the xi coordinates of both element
  coordinates and its faces.

  Note: this matches the ordering of the entities returned by
  [`apf.getNodeEntities`](@ref)

  Inputs:
    order: order of the element
    dim: dimensionality of the entity

  Outputs:
    node_xi: a dim x numnodes array containing the xi coordinates of each node

  Note that this follows the Pumi convention that the reference element is
  a right triangle with the right angle in the lower left corner.  Barycentric
  coordinates are used and span from 0 to 1.  The right angle corner has
  coordinate xi3, the next corner counter-clowckwise is xi1, and the next is
  xi2.  The array contains xi1 and xi2 along the first dimension.  However,
  the ordering of the vertices is counter-clockwise from the bottom left,
  ie. the 1st vertex has coordinate xi3.  This is rather confusing and the
  source of many problems.

  The easiest way to see the Pumi convention is to look at the shape functions
  for the lagrangian elements in apfShape.cc
"""
function getXiCoords(order::Integer, dim::Integer)

  if dim == 1
    if order == 1
      xi = [0. 1;]
    elseif order == 2
      xicoords = [0. 1 0.5;]
    else
      throw(ErrorException("unsupported order $order for dimension $dim"))
    end
  elseif dim == 2
    if order == 1
      xicoords = [0. 1 0;
                  0 0 1]
    elseif order == 2
      xicoords = [0. 1 0 0.5 0.5 0.0;
                  0 0 1 0.0 0.5 0.5]
    else
      throw(ErrorException("unsupported order $order for dimension $dim"))
    end
  elseif dim == 3
    if order == 1
      xicoords = [0. 1 0 0;
                  0 0 1 0;
                  0  0 0 1]
    elseif order == 2
      xicoords = [0. 1 0 0 0.5 0.5 0.0 0.0 0.5 0.0;
                  0 0 1 0 0.0 0.5 0.5 0.0 0.0 0.5;
                  0 0 0 1 0.0 0.0 0.0 0.5 0.5 0.5]
    else
      throw(ErrorException("unsupported order $order for dimension $dim"))
    end
  else
    throw(ErrorException("unsupported dimension $dim"))
  end

end


text:
# Types and functions
# -------------------
abstract type AbstractELM end

struct ELM{T,I, S<:AbstractArray{T}, F} <: AbstractELM

    n_hidden_neurons::I
    weight_matrix::S
    bias_vector::S
    activation::F
    output_weights::S
    β::T
end

"""
    ELM(n_hidden_neurons::Int,input_data::AbstractArray{T},output_data::AbstractArray{T}; activation::Function = sigmoid, regularization::T = zero(T)) where T<:AbstractFloat
Construct an ELM passing a number of neurons, the inputs and outputs.
As a keyword argument, you can also pass a different activation function (Default = sigmoid).
Inputs should be in the format below.
| Feature1| Feature2|
|---------|---------|
|  entry1 |  entry1 |
|  entry2 |  entry2 |
"""
function ELM(n_hidden_neurons::Int,
    input_data::AbstractArray{T},
    output_data::AbstractArray{T};
    activation::Function = sigmoid,
    regularization::T = zero(T)) where T<:AbstractFloat

in_size = size(input_data, 2)
out_size = size(output_data, 2)
weight_matrix = rand(T, in_size, n_hidden_neurons) .* 2 .- 1
output_weights = zeros(T, n_hidden_neurons,  out_size)
bias_vector = rand(T, 1, n_hidden_neurons)

return ELM{T,typeof(n_hidden_neurons),
        typeof(weight_matrix),
        typeof(activation)
        }(n_hidden_neurons,
        weight_matrix,
        bias_vector,
        activation,
        output_weights,
        regularization)
end

"""
    fit!(elm::ELM, x::AbstractArray{T}, y::AbstractArray{T}) where T<:AbstractFloat
Train the ELM output weights.
"""

function fit!(elm::ELM,
    x::AbstractArray{T},
    y::AbstractArray{T}) where T<:AbstractFloat

hidden_matrix = elm.activation.(x*elm.weight_matrix .+ elm.bias_vector)
elm.output_weights .= hidden_matrix\y
nothing
end
"""
    predict(elm::ELM,x::AbstractArray{T}) where T<:AbstractFloat
Predict new values
"""
function predict(elm::ELM,
    x::AbstractArray{T}) where T<:AbstractFloat

hidden_matrix = elm.activation.(x*elm.weight_matrix .+ elm.bias_vector)
prediction = hidden_matrix*elm.output_weights
end

struct ELMrow{T,I, S<:AbstractArray{T}, F} <: AbstractELM

    n_hidden_neurons::I
    weight_matrix::S
    bias_vector::S
    activation::F
    output_weights::S
    β::T
end

function ELMrow(n_hidden_neurons::Int,
        input_data::AbstractArray{T},
        output_data::AbstractArray{T};
        activation::Function = sigmoid,
        regularization::T = zero(T)) where T<:AbstractFloat

    in_size = size(input_data, 1)
    out_size = size(output_data, 1)
    weight_matrix = rand(T,n_hidden_neurons, in_size) .* 2 .- 1
    output_weights = zeros(T,out_size,n_hidden_neurons)
    bias_vector = rand(T,n_hidden_neurons,1)
    return ELM{T,typeof(n_hidden_neurons),
            typeof(weight_matrix),
            typeof(activation)
            }(n_hidden_neurons,
            weight_matrix,
            bias_vector,
            activation,
            output_weights,
            regularization)

end

function fit!(elm::ELMrow,
            x::AbstractArray{T},
            y::AbstractArray{T}) where T<:AbstractFloat

    hidden_matrix = elm.activation.(elm.weight_matrix*x .+ elm.bias_vector)
    # elm.output_weights .= y*pinv(hidden_matrix)
    elm.output_weights .= hidden_matrix\y
    nothing
end

function predict(elm::ELMrow,
                x::AbstractArray{T}) where T<:AbstractFloat

    hidden_matrix = elm.activation.(elm.weight_matrix*x .+ elm.bias_vector)
    prediction = elm.output_weights*hidden_matrix
end


text:
# Autogenerated wrapper script for SoXResampler_jll for x86_64-apple-darwin
export libsoxr, libsoxr_lsr

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SoXResampler")
JLLWrappers.@declare_library_product(libsoxr, "@rpath/libsoxr.0.dylib")
JLLWrappers.@declare_library_product(libsoxr_lsr, "@rpath/libsoxr-lsr.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsoxr,
        "lib/libsoxr.0.1.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsoxr_lsr,
        "lib/libsoxr-lsr.0.1.9.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function gather(input::AbstractArray{T,N}, index::AbstractArray{<:Integer,N}, dims::Integer;
                out::AbstractArray{T,N}=similar(index, T)) where {T,N}
    @assert dims <= N "Specified dimensions must lower or equal to the rank of input matrix."

    @inbounds for x = CartesianIndices(out)
        tup = collect(Tuple(x))
        tup[dims] = index[x]
        out[x] = input[tup...]
    end
    return out
end

function gather(input::Matrix{T}, index::Array{Int}) where T
    out = Array{T}(undef, size(input,1), size(index)...)
    @inbounds for ind = CartesianIndices(index)
        out[:, ind] = input[:, index[ind]]
    end
    return out
end

identity(; kwargs...) = kwargs.data

struct GraphInfo{A,T<:Integer}
    adj::AbstractVector{A}
    edge_idx::A
    V::T
    E::T

    function GraphInfo(adj::AbstractVector{<:AbstractVector{<:Integer}})
        V = length(adj)
        edge_idx = edge_index_table(adj, V)
        E = edge_idx[end]
        new{typeof(edge_idx),typeof(V)}(adj, edge_idx, V, E)
    end
end

function edge_index_table(adj::AbstractVector{<:AbstractVector{<:Integer}},
                          N::Integer=size(adj,1))
    y = similar(adj[1], N+1)
    y .= 0, cumsum(map(length, adj))...
    y
end