text:
using Images, Colors
import ImageView

z = ones(10,50);
y = 8; x = 2;
z[y,x] = 0
zimg = convert(Image, z)
imgc, img2 = ImageView.view(zimg,pixelspacing=[1,1]);
Tk.set_size(ImageView.toplevel(imgc), 200, 200)
idx = ImageView.annotate!(imgc, img2, ImageView.AnnotationText(x, y, "x", color=RGB(0,0,1), fontsize=3))
idx2 = ImageView.annotate!(imgc, img2, ImageView.AnnotationPoint(x+10, y, shape='.', size=4, color=RGB(1,0,0)))
idx3 = ImageView.annotate!(imgc, img2, ImageView.AnnotationPoint(x+20, y-6, shape='.', size=1, color=RGB(1,0,0), linecolor=RGB(0,0,0), scale=true))
idx4 = ImageView.annotate!(imgc, img2, ImageView.AnnotationLine(x+10, y, x+20, y-6, linewidth=2, color=RGB(0,1,0)))
idx5 = ImageView.annotate!(imgc, img2, ImageView.AnnotationBox(x+10, y, x+20, y-6, linewidth=2, color=RGB(0,0,1)))
idx6 = ImageView.annotate!(imgc, img2, ImageView.AnnotationPoints([(x+25, y-3), (x+24, y-3), (x+23, y-3)],
                                                                  shape='o', size=3, color=RGB(1,0,1),
                                                                  linecolor=RGB(0,0,0), scale=true))

# Also test the following:
# ImageView.delete!(imgc, idx)


text:
"""
    Electromagnetic

Represents the physical properties for a homogenous isotropic electromagnetic medium. Produces
a three dimensional vector field.
"""
struct Electromagnetic{Dim,T} <: PhysicalMedium{T,Dim,3}
    μ::Complex{T} # Permeability
    ε::Complex{T} # Permittivity
    σ::Complex{T} # Conductivity
end

name(e::Electromagnetic{T,Dim}) where {T,Dim} = "$(Dim)D Electromagnetic"

ext: jl


text:
export VaxFloatG

struct VaxFloatG <: VaxFloat
    x::UInt64

    VaxFloatG(x::UInt64) = new(ltoh(x))
end

function VaxFloatG(x::T) where T <: Real
    y = reinterpret(UInt64, convert(Float64, x))

    part1 = y & bmask32
    part2 = (y >>> 32) & bmask32

    if ENDIAN_BOM === 0x04030201
        vaxpart2 = part1
        ieeepart1 = part2
    else
        vaxpart2 = part2
        ieeepart1 = part1
    end

    e = reinterpret(Int64, ieeepart1 & IEEE_T_EXPONENT_MASK)

    if ieeepart1 & ~SIGN_BIT_64 === zero(UInt64)
        # ±0.0 becomes 0.0
        return zero(VaxFloatG)
    elseif e === IEEE_T_EXPONENT_MASK
        # Vax types don't support ±Inf or NaN
        throw(InexactError(:VaxFloatG, VaxFloatG, x))
    else
        e >>>= VAX_G_MANTISSA_SIZE
        m = ieeepart1 & VAX_G_MANTISSA_MASK

        if e === zero(Int64)
            m = (m << 1) | (vaxpart2 >>> 31)
            vaxpart2 <<= 1
            while m & VAX_G_HIDDEN_BIT === zero(UInt64)
                m = (m << 1) | (vaxpart2 >>> 31)
                vaxpart2 <<= 1
                e -= one(Int64)
            end
            m &= VAX_G_MANTISSA_MASK
        end

        e += one(Int64) + VAX_G_EXPONENT_BIAS - IEEE_T_EXPONENT_BIAS
        if e <= zero(Int64)
            # Silent underflow
            return zero(VaxFloatG)
        elseif e > (2*VAX_G_EXPONENT_BIAS - 1)
            # Overflow
            throw(InexactError(:VaxFloatG, VaxFloatG, x))
        else
            vaxpart = (ieeepart1 & SIGN_BIT_64) | (e << VAX_G_MANTISSA_SIZE) | m
        end
    end

    vaxpart_1 = vaxpart & bmask16
    vaxpart_2 = (vaxpart >>> 16) & bmask16

    vaxpart_3 = vaxpart2 & bmask16
    vaxpart_4 = (vaxpart2 >>> 16) & bmask16

    res = htol((vaxpart_3 << 48) |
          (vaxpart_4 << 32) |
          (vaxpart_1 << 16) |
          vaxpart_2)

    return VaxFloatG(res)
end

Base.typemax(::Type{VaxFloatG}) = VaxFloatG(0xffffffffffff7fff)
Base.floatmax(::Type{VaxFloatG}) = VaxFloatG(0xffffffffffff7fff)
Base.floatmin(::Type{VaxFloatG}) = VaxFloatG(0x0000000000000010)
Base.typemin(::Type{VaxFloatG}) = VaxFloatG(0x0000000000000010)

Base.zero(::Type{VaxFloatG}) = VaxFloatG(0x0000000000000000)
Base.one(::Type{VaxFloatG}) = VaxFloatG(0x0000000000004010)

function Base.convert(::Type{Float64}, x::VaxFloatG)
    y = ltoh(x.x)

    vaxpart_1 = y & bmask16
    vaxpart_2 = (y >>> 16) & bmask16
    vaxpart1 = (vaxpart_1 << 16) | vaxpart_2

    vaxpart_3 = (y >>> 32) & bmask16
    vaxpart_4 = (y >>> 48) & bmask16
    vaxpart2 = (vaxpart_3 << 16) | vaxpart_4

    e = reinterpret(Int64, vaxpart1 & VAX_G_EXPONENT_MASK)
    if e === zero(Int64)
        if vaxpart1 & SIGN_BIT_64 === SIGN_BIT_64
            # Reserved floating-point reserved operand
            throw(InexactError(:convert, Float64, x))
        end

        # Dirty zero
        return zero(Float64)
    else
        e >>>= VAX_G_MANTISSA_SIZE

        e -= one(Int64) + VAX_G_EXPONENT_BIAS - IEEE_T_EXPONENT_BIAS
        if e > zero(Int64)
            ieeepart1 = vaxpart1  - ((UNO64 + VAX_G_EXPONENT_BIAS - IEEE_T_EXPONENT_BIAS) << IEEE_T_MANTISSA_SIZE)
            ieeepart2 = vaxpart2
        else
            # Subnormal result
            vaxpart1 = (vaxpart1 & (SIGN_BIT_64 | VAX_G_MANTISSA_MASK)) | VAX_G_HIDDEN_BIT
            ieeepart1 = (vaxpart1 & SIGN_BIT_64) | ((vaxpart1 & (VAX_G_HIDDEN_BIT | VAX_G_MANTISSA_MASK)) >>> (1 - e))
            ieeepart2 = (vaxpart1 << (31 + e)) | (vaxpart2 >>> (1 - e))
        end

        if ENDIAN_BOM === 0x04030201
            out1 = ieeepart2
            out2 = ieeepart1
        else
            out1 = ieeepart1
            out2 = ieeepart2
        end
    end

    res = (out2 << 32) | out1

    return reinterpret(Float64, res)
end

function Base.convert(::Type{T}, x::VaxFloatG) where T <: Union{Float16,Float32,BigFloat,Integer}
    return convert(T, convert(Float64, x))
end


text:
@symbol_func function P_W(cur_reactor::AbstractReactor)
  cur_P = K_W(cur_reactor)

  cur_P *= safe_symbol(cur_reactor, :n_bar) ^ 2

  cur_P *= safe_symbol(cur_reactor, :R_0)

  cur_P *= cur_reactor.sigma_v

  cur_P
end


text:
function tomography(probabilities::Dict{Tuple,<:Dict}, sites::Vector{<:Index};
                    method::String="linear_inversion",
                    fillzeros::Bool=true,
                    process::Bool = false,
                    trρ::Number = 1.0,
                    max_iters::Int=10000,
                    kwargs...)

  # Generate the projector matrix corresponding to the probabilities.
  A, p = design_matrix(probabilities; return_probs = true, process = process)

  if (method == "LI"  || method == "linear_inversion")
    # Invert the Born rule and reshape
    ρ_vec = pinv(A) * p
    d = Int(sqrt(size(ρ_vec,1)))
    ρ̂ = reshape(ρ_vec,(d,d))

    ρ̂ .= ρ̂ * (trρ / tr(ρ̂))
    # Make PSD
    ρ̂ = make_PSD(ρ̂)
  else
    d = Int(sqrt(size(A,2)))
    N = Int(sqrt(d))
    n = N ÷ 2

    # Variational density matrix
    ρ = Convex.ComplexVariable(d,d)

    if (method == "LS"  || method == "least_squares")
      # Minimize the cost function C = ||A ρ⃗ - p̂||²
      cost_function = Convex.norm(A * vec(ρ) - p)
    elseif (method == "MLE" || method == "maximum_likelihood")
      # Minimize the negative log likelihood:
      cost_function = - p' * Convex.log(real(A * vec(ρ)) + 1e-10)
    else
      error("Tomography method not recognized
       Currently available methods: - LI  : linear inversion
                                      LS  : least squares
                                      MLS : maximum likelihood")
    end

    # Contrained the trace and enforce positivity and hermitianity
    if process
      function tracepreserving(ρ)
        for j in 1:n
          subsystem_dims = [2 for _ in 1:(N+1-j)]
          ρ = Convex.partialtrace(ρ,j+1,subsystem_dims)
        end
        return ρ
      end

      constraints = [Convex.tr(ρ) == (1<<n) * trρ
                    Convex.isposdef(ρ)
                    tracepreserving(ρ) == Matrix{Float64}(I,1<<n,1<<n)
                    ρ == ρ']
    else
      constraints = [Convex.tr(ρ) == trρ
                     Convex.isposdef(ρ)
                     ρ == ρ']
    end
    # Use Convex.jl to solve the optimization
    problem = Convex.minimize(cost_function,constraints)
    Convex.solve!(problem, () -> SCS.Optimizer(verbose=false,max_iters=max_iters),verbose=false)
    ρ̂ = ρ.value
  end
  return PastaQ.itensor(ρ̂, sites)
end

tomography(data::Matrix{Pair{String, Int}}; method::String="linear_inversion", fillzeros::Bool=true, kwargs...) =
  tomography(empirical_probabilities(data; fillzeros=fillzeros), siteinds("Qubit", size(data,2)); method = method, kwargs...)

tomography(data::Matrix{Pair{String, Int}}, sites::Vector{<:Index}; method::String="linear_inversion", fillzeros::Bool=true, kwargs...) =
  tomography(empirical_probabilities(data; fillzeros=fillzeros), sites; method = method, kwargs...)

tomography(data::Matrix{Pair{String,Pair{String, Int}}}; kwargs...) =
  tomography(data, siteinds("Qubit", size(data,2)); kwargs...)

function tomography(data::Matrix{Pair{String,Pair{String, Int}}}, sites::Vector{<:Index}; method::String="linear_inversion", fillzeros::Bool=true, kwargs...)
  sites_in  = addtags.(sites, "Input")
  sites_out = addtags.(sites, "Output")
  process_sites = Index[]
  for j in 1:size(data,2)
    push!(process_sites, sites_in[j])
    push!(process_sites, sites_out[j])
  end
  tomography(empirical_probabilities(data; fillzeros=fillzeros), process_sites; method = method, process = true, kwargs...)
end

"""
    tomography(train_data::Matrix{Pair{String,Pair{String, Int}}}, L::LPDO;
               optimizer::Optimizer,
               observer! = nothing,
               batchsize::Int64 = 100,
               epochs::Int64 = 1000,
               kwargs...)

Run quantum process tomography using a variational model `L` to fit `train_data`.
The model can be either a unitary circuit (MPO) or a Choi matrix (LPDO).

# Arguments:
- `train_data`: pairs of preparation/ (basis/outcome): `("X+"=>"X"=>0, "Z-"=>"Y"=>1, "Y+"=>"Z"=>0, …)`.
 - `L`: variational model (MPO/LPDO).
 - `optimizer`: algorithm used to update the model parameters.
 - `observer!`: if provided, keep track of training metrics.
 - `batch_size`: number of samples used for one gradient update.
 - `epochs`: number of training iterations.
 - `target`: target quantum process (if provided, compute fidelities).
 - `test_data`: data for computing cross-validation.
 - `outputpath`: if provided, save metrics on file.
"""
function tomography(
  train_data::AbstractMatrix, L::LPDO; (observer!)=nothing, kwargs...
)
  # Read arguments
  opt = get(kwargs, :optimizer, Optimisers.Descent(0.01))
  batchsize::Int64 = get(kwargs, :batchsize, 100)
  epochs::Int64 = get(kwargs, :epochs, 1000)
  trace_preserving_regularizer = get(kwargs, :trace_preserving_regularizer, 0.0)
  observe_step::Int64 = get(kwargs, :observe_step, 1)
  test_data = get(kwargs, :test_data, nothing)
  outputpath = get(kwargs, :fout, nothing)
  print_metrics = get(kwargs, :print_metrics, [])
  outputpath = get(kwargs, :outputpath, nothing)
  outputlevel = get(kwargs, :outputlevel, 1)
  savestate = get(kwargs, :savestate, false)
  earlystop = get(kwargs, :earlystop, false)

  model = copy(L)
  isqpt = train_data isa Matrix{Pair{String,Pair{String,Int}}}
  localnorm = isqpt ? 2.0 : 1.0

  # observer is not passed but earlystop is called
  observer! = (isnothing(observer!) && earlystop) ? Observer() : observer!

  # observer is defined
  if !isnothing(observer!)
    observer!["train_loss"] = nothing
    if !isnothing(test_data)
      observer!["test_loss"] = nothing
    end
    # add the standard early stop function to the observer
    if earlystop
      stop_if(; loss::Vector) = stopif_loss(; loss = loss, ϵ = 1e-3, min_iter = 10)
      observer!["earlystop"] = stopif
    end
  end

  # initialize optimizer
  st = PastaQ.state(opt, model)
  optimizer = (opt, st)

  @assert size(train_data, 2) == length(model)
  !isnothing(test_data) && @assert size(test_data)[2] == length(model)

  batchsize = min(size(train_data)[1], batchsize)
  num_batches = Int(floor(size(train_data)[1] / batchsize))

  tot_time = 0.0
  best_model = nothing
  best_testloss = 1000.0
  test_loss = nothing

  # Training iterations
  for ep in 1:epochs
    ep_time = @elapsed begin
      train_data = train_data[shuffle(1:end), :]
      train_loss = 0.0

      # Sweep over the data set
      for b in 1:num_batches
        batch = train_data[((b - 1) * batchsize + 1):(b * batchsize), :]

        normalized_model = copy(model)
        sqrt_localnorms = []
        normalize!(normalized_model;
                   (sqrt_localnorms!)=sqrt_localnorms,
                   localnorm=localnorm)

        grads, loss = gradients(
          normalized_model,
          batch;
          sqrt_localnorms=sqrt_localnorms,
          trace_preserving_regularizer=trace_preserving_regularizer,
        )

        nupdate = ep * num_batches + b
        train_loss += loss / Float64(num_batches)
        update!(model, grads, optimizer)
      end
    end # end @elapsed
    !isnothing(observer!) && push!(last(observer!["train_loss"]), train_loss)
    tot_time += ep_time

    # measurement stage
    if ep % observe_step == 0
      normalized_model = copy(model)
      sqrt_localnorms = []
      normalize!(normalized_model;
                 (sqrt_localnorms!)=sqrt_localnorms,
                 localnorm=localnorm)

      if !isnothing(test_data)
        test_loss = nll(normalized_model, test_data)
        !isnothing(observer!) && push!(last(observer!["test_loss"]), test_loss)
        if test_loss ≤ best_testloss
          best_testloss = test_loss
          best_model = copy(normalized_model)
        end
      else
        best_model = copy(normalized_model)
      end

      # update observer
      if !isnothing(observer!)
        loss = (!isnothing(test_data) ? results(observer!, "test_loss") :
                                        results(observer!, "train_loss"))

        model_to_observe = (isqpt && (normalized_model isa LPDO{MPS}) ? choi_mps_to_unitary_mpo(normalized_model) :
                                  !isqpt && (normalized_model isa LPDO{MPS}) ? normalized_model.X :
                                                                               normalized_model)
        update!(observer!, model_to_observe; train_loss = train_loss,
                                             test_loss  = test_loss)
      end

      # printing
      if outputlevel ≥ 1
        @printf("%-4d  ", ep)
        @printf("⟨logP⟩ = %-4.4f  ", train_loss)
        if !isnothing(test_data)
          @printf("(%.4f)  ", test_loss)
        end
        # TODO: add the trace preserving cost function here for QPT
        !isnothing(observer!) && printobserver(observer!, print_metrics)
        @printf("elapsed = %-4.3fs", ep_time)
        println()
      end
      # saving
      if !isnothing(outputpath)
        observerpath = outputpath * "_observer.jld2"
        save(observerpath, observer!)
        if savestate
          if isqpt
            model_to_be_saved = model isa LPDO{MPS} ? choi_mps_to_unitary_mpo(best_model) : best_model
          else
            model_to_be_saved = model isa LPDO{MPS} ? best_model.X : best_model
          end
          statepath = outputpath * "_state.h5"
          h5rewrite(statepath) do fout
            write(fout, "state", model_to_be_saved)
          end
        end
      end
    end
    !isnothing(observer!) && haskey(observer!.data,"earlystop") && results(observer!, "earlystop")[end] && break
  end
  return best_model
end

# QST
tomography(data::Matrix{Pair{String,Int}}, ψ::MPS; kwargs...) =
  tomography(data, LPDO(ψ); kwargs...).X

tomography(train_data::Vector{<:Vector{Pair{String,Int}}}, args...; kwargs...) =
  tomography(permutedims(hcat(train_data...)), args...; kwargs...)

# QPT
tomography(data::Matrix{Pair{String,Pair{String,Int}}}, U::MPO; kwargs...) =
  choi_mps_to_unitary_mpo(tomography(data, LPDO(unitary_mpo_to_choi_mps(U)); kwargs...))

tomography(train_data::Vector{<:Vector{Pair{String,Pair{String,Int}}}}, args...; kwargs...) =
  tomography(permutedims(hcat(train_data...)), args...; kwargs...)

"""
EARLY STOPPING FUNCTIONS
"""

#stopif_fidelity(M1, M2; ϵ::Number, kwargs...)
#  fidelity(M1,M2) ≤ ϵ

function stopif_loss(; loss::Vector, ϵ::Number, min_iter::Number)
  length(loss) < min_iter+1 && return false
  avgloss = StatsBase.mean(loss[end-size:end])
  Δ = StatsBase.sem(historyloss[end-size:end])
  return Δ/avgloss < ϵ
end


text:
"""
# PSH

Module PSH converts seismograms in the ZNE frame into the PSH frame using the free-surface
transform.
"""
module PSH

using SAC

export psh, psh!

"""
    psh_matrix!(α, β, p, M) -> M

Fill in the PSH transfer matrix `M`, given the P-wave velocity `α`, S-wave velocity `β`
and the slowness `p`.  Velocities are in km/s, slowness in s/km.
"""
function psh_matrix!{T}(a, b, p, M::AbstractArray{T,2})
    size(M) == (3,3) || error("Transfer matrix must have size (3,3)")
    isevanescent(a, b, p) &&
        warn("psh_matrix!: Evanescent wave (vp=$a, vs=$b, p=$p)")
    qa = sqrt(1/a^2 - p^2)
    qb = sqrt(1/b^2 - p^2)
    # Kennett, GJI, 1991
    Vpz = -(1 - 2b^2*p^2)/(2a*qa)
    Vpr = p*b^2/a
    Vsz = p*b
    Vsr = (1 - 2b^2*p^2)/(2b*qb)
    Vht = 1/2
    M[:,:] .= [Vpz Vsz  0
               Vpr Vsr  0
                0   0  Vht]
end

"""
    psh_matrix(α, β, p) -> M

Return the PSH transfer matrix M, given the P-wave velocity `α`, S-wave velocity `β`
and the slowness `p`.  Velocities are in km/s, slowness in s/km.
"""
psh_matrix(a, b, p) = psh_matrix!(a, b, p, Array{Complex{Float64}}(3,3))

"""
    isevanescent(α, β, p) -> ::Bool

Return true if the combination of Vp `α`, Vs `β` and slowness `p` lead to an evanescent
wave at the surface.
"""
function isevanescent(a, b, p)
    p = Float64(p)
    a, b = Float64(a), Float64(b)
    1/a^2 - p^2 <= 0.0 || 1/b^2 - p^2 <= 0.0
end

"""
    psh!(Z, R, T, α, β, p) -> P, S, H

Converts three arrays, `SACtr`s or arrays of `SACtr`s in ZRT orientation into
PSH orientation, using the supplied Vp `α`, Vs `β` (km/s) and slowness `p` (s/km).

Traces are updated to be in the order P (like vertical), S (like radial), H (like transverse).
"""
function psh!{F}(Z::Array{F,1}, R::Array{F,1}, T::Array{F,1}, a, b, p)
	length(Z) == length(R) == length(T) || error("psh!: Traces must be same length")
	M = psh_matrix(a, b, p)
    # Kennett 1991
    Z, R, T = M[1,1]*Z + M[1,2]*R + M[1,3]*T,
              M[2,1]*Z + M[2,2]*R + M[2,3]*T,
              M[3,1]*Z + M[3,2]*R + M[3,3]*T
    Z, R, T
end
function psh!(Z::SACtr, R::SACtr, T::SACtr, a, b, p)
    P, S, H = psh(Z.t, R.t, T.t, a, b, p)
    Z.t .= real.(P)
    R.t .= real.(S)
    T.t .= real.(H)
    for t in (Z, R, T) SAC.update_headers!(t) end
    Z.kcmpnm, R.kcmpnm, T.kcmpnm = "P", "S", "H"
    Z, R, T
end
function psh!(Z::Array{SACtr}, R::Array{SACtr}, T::Array{SACtr}, a, b, p)
    for (z, r, t) in zip(Z, R, T) psh!(z, r, t, a, b, p) end
    Z, R, T
end

"""
    psh(Z, R, T, α, β, p) -> P, S, H

Convert three arrays, `SACtr`s or arrays of `SACtr`s in ZRT orientation into
PSH orientation, using the supplied Vp `α`, Vs `β` (km/s) and slowness `p` (s/km),
where `P` is parallel to P-wave motion positive away from the source, `S` is parallel
to SV and `H` is the transverse direction.
"""
psh(Z, R, T, a, b, p) = psh!(deepcopy(Z), deepcopy(R), deepcopy(T), a, b, p)

end # module

max_stars_repo_name: anowacki/assorted-julia-modules

max_issues_repo_name: anowacki/assorted-julia-modules

max_forks_repo_name: anowacki/assorted-julia-modules


text:
using Flux
using Flux: @functor
using Parameters: @with_kw
using FastAI
import FastAI.Models.tabular_embedding_backbone
include("../layers/categorical_embeddings.jl")

@with_kw mutable struct MLPArgs
    embedding_dims::Vector{Tuple{Int64, Int64}}
    cont_input_dim::Int64
    output_dim::Int64 = 2
    lr::Float64 = 1e-2		# learning rate
    epochs::Int64 = 10        # number of epochs
    use_cuda::Bool = true   # use gpu (if cuda available)
    dropout::Float64 = 0.10 # dropout from dense layers
    hidden_dims::Tuple = (128, 64) # Size of dense hidden layers
    embedding_dropout::Float64 = 0.10 # dropout for categorical embedding
    batchnorm::Bool = true # batchnorm on dense layers
    linear_first::Bool = true # linear layer before or after batch norm
    activation = Flux.relu
    output_activation = Flux.sigmoid
end

"""
Create a series of dense layers following the
embedding and continious data vcat

https://github.com/FluxML/FastAI.jl/blob/master/src/models/tabularmodel.jl
"""
function dense_layers(args)
    cat_input_dim = sum([i[1] for i in args.embedding_dims])

    layers = []
    first_layer = linbndrop(cat_input_dim + args.cont_input_dim, first(args.hidden_dims);
                    use_bn=args.batchnorm, p=args.dropout, lin_first=args.linear_first,
                    act=args.activation)
    push!(layers, first_layer)

    for (isize, osize) in zip(args.hidden_dims[1:(end-1)], args.hidden_dims[2:end])
        layer = linbndrop(isize, osize; use_bn=args.batchnorm, p=args.dropout,
                        lin_first=args.linear_first, act=args.activation)
        push!(layers, layer)
    end

    return Chain(layers...)

end

struct MLP
    cat_backbone
    cont_backbone
    dense
    output
end
@functor MLP

MLP(args::MLPArgs) = MLP(
    tabular_embedding_backbone(args.embedding_dims, args.embedding_dropout),
    BatchNorm(args.cont_input_dim),
    dense_layers(args),
    Dense(args.hidden_dims[end], args.output_dim, args.output_activation)
)

function (mlp::MLP)(X_cat, X_cont)
    h_cat = mlp.cat_backbone(X_cat...)
    h_cont = mlp.cont_backbone(X_cont)
    h2 = mlp.dense(vcat(h_cat, h_cont))
    mlp.output(h2)
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    wclogptype,
    wclogphydrogentype,
    wclogpcontrib

const WCLOGP_TABLE = YAML.load(open(
    joinpath(dirname(@__FILE__), "..", "assets", "const", "wclogp.yaml")
))["logP"]

const ALIPH_HETERO = Set([:H, :C, :N, :O, :P, :S, :F, :Cl, :Br, :I])
const AROM_HETERO = Dict(
    :H  => :C18, :C => :C21, :N => :C22, :O => :C23, :S => :C24,
    :F => :C14, :Cl => :C15, :Br => :C16, :I => :C17
)
# TODO:
const P_BLOCK = (:Al, :B, :Si, :Ga, :Ge, :As, :Se, :Sn, :Te, :Pb,
                 :Ne, :Ar, :Kr, :Xe, :Rn)
const D_BLOCK = (:Fe, :Co, :Cu, :Zn, :Tc, :Cd, :Pt, :Au, :Hg, :Gd)

"""
    wclogptype(mol::GraphMol)

Return Wildman-Crippen LogP atom types.
"""
@cache function wclogptype(mol::GraphMol)
    atomtypes = fill(:undef, nodecount(mol))
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    bondorder_ = bondorder(mol)
    nodedegree_ = nodedegree(mol)
    hcount_ = hcount(mol)
    pielectron_ = pielectron(mol)
    isaromatic_ = isaromatic(mol)
    isaromaticbond_ = isaromaticbond(mol)
    # Carbons
    for i in findall(atomsymbol_ .== :C)
        if pielectron_[i] == 0
            # Aliphatic (C1-4,8-12,27)
            nbrs = collect(adjacencies(mol, i))
            if !isempty(setdiff(atomsymbol_[nbrs], ALIPH_HETERO))
                atomtypes[i] = :C27 # Adjacent to inorganic atoms
            elseif all(isaromatic_[nbrs] .== false)
                if all(atomsymbol_[nbrs] .== :C)
                    # Aliphatic carbon (C1,2)
                    if nodedegree_[i] <= 2
                        atomtypes[i] = :C1
                    else
                        atomtypes[i] = :C2
                    end
                else
                    # Adjacent to heteroatoms (C3,4)
                    if nodedegree_[i] <= 2
                        atomtypes[i] = :C3
                    else
                        atomtypes[i] = :C4
                    end
                end
            else
                # Adjacent to aromatic atoms (C8-12)
                if nodedegree_[i] == 1
                    nbr = nbrs[1]
                    if atomsymbol_[nbr] == :C
                        atomtypes[i] = :C8
                    else
                        atomtypes[i] = :C9
                    end
                elseif nodedegree_[i] == 2
                    atomtypes[i] = :C10
                elseif nodedegree_[i] == 3
                    atomtypes[i] = :C11
                elseif nodedegree_[i] == 4
                    atomtypes[i] = :C12
                end
            end
        elseif isaromatic_[i]
            # Aromatic (C13-25)
            if degree(mol, i) == 2
                atomtypes[i] = :C18 # Aromatic non-substituted
                continue
            end
            subst = nothing
            substbond = nothing
            aromcnt = 0
            for (inc, adj) in neighbors(mol, i)
                if isaromaticbond_[inc]
                    aromcnt += 1
                else
                    subst = adj
                    substbond = inc
                end
            end
            if aromcnt == 3
                atomtypes[i] = :C19 # Bridgehead
            elseif !haskey(AROM_HETERO, atomsymbol_[subst])
                atomtypes[i] = :C13 # Inorganic substituent
            elseif isaromatic_[subst]
                atomtypes[i] = :C20 # Aromatic substituent
            elseif bondorder_[substbond] == 2
                atomtypes[i] = :C25 # Double bond substituent
            else
                # Typical substituent (C14-17,21-24)
                atomtypes[i] = AROM_HETERO[atomsymbol_[subst]]
            end
        else
            # Aliphatic multiple bond (C5-7,26)
            nbrs = collect(adjacencies(mol, i))
            bonds = collect(incidences(mol, i))
            if any(bondorder_[bonds] .== 3)
                atomtypes[i] = :C7 # Alkyne, Nitrile
            elseif any((pielectron_[nbrs] .> 0) .* (atomsymbol_[nbrs] .!= :C))
                atomtypes[i] = :C5 # Double bond to non-C
            elseif any(isaromatic_[nbrs] .== true)
                atomtypes[i] = :C26 # Adjacent to aromatic
            else
                atomtypes[i] = :C6 # Double bond to C (including allene)
            end
        end
        if atomtypes[i] == :undef
            atomtypes[i] = :CS # Not found
        end
    end

    # Nitrogens
    for i in findall(atomsymbol_ .== :N)
        if charge_[i] > 0
            nbrs = collect(adjacencies(mol, i))
            if isaromatic_[i]
                atomtypes[i] = :N12 # Charged aromatic nitrogen
            elseif hcount_[i] == 0
                if nodedegree_[i] == 2 && all(atomsymbol_[nbrs] .== :N)
                    atomtypes[i] = :N14 # Azide is exceptionally N14
                else
                    atomtypes[i] = :N13 # Quart-ammonium
                end
            elseif pielectron_[i] == 0
                atomtypes[i] = :N10 # Protonated amine
            else
                atomtypes[i] = :N14 # Other charged
            end
        elseif charge_[i] < 0
            atomtypes[i] = :N14 # Neg charged
        elseif isaromatic_[i]
            atomtypes[i] = :N11 # Neutral aromatic
        elseif pielectron_[i] == 2
            atomtypes[i] = :N9 # Nitrile
        elseif pielectron_[i] == 1
            # Imine (N5,6)
            if nodedegree_[i] == 1
                atomtypes[i] = :N5
            elseif nodedegree_[i] == 2
                atomtypes[i] = :N6
            end
        else
            nbrs = collect(adjacencies(mol, i))
            if all(isaromatic_[nbrs] .== false)
                # Aliphatic amine (N1,2,7)
                if nodedegree_[i] == 1
                    atomtypes[i] = :N1
                elseif nodedegree_[i] == 2
                    atomtypes[i] = :N2
                elseif nodedegree_[i] == 3
                    atomtypes[i] = :N7
                end
            else
                # Aromatic amine (N3,4,8)
                if nodedegree_[i] == 1
                    atomtypes[i] = :N3
                elseif nodedegree_[i] == 2
                    atomtypes[i] = :N4
                elseif nodedegree_[i] == 3
                    atomtypes[i] = :N8
                end
            end
        end
        if atomtypes[i] == :undef
            # Not found
            atomtypes[i] = :NS
        end
    end

    # Oxygens
    for i in findall(atomsymbol_ .== :O)
        if isaromatic_[i]
            # Aromatic oxygen (O1)
            atomtypes[i] = :O1
            continue
        end
        nbrs = collect(adjacencies(mol, i))
        bonds = collect(incidences(mol, i))
        if nodedegree_[i] == 1
            # Hydroxyl (O2,5-12)
            nbr = nbrs[1]
            bond = bonds[1]
            if atomtypes[nbr] == :C5 && charge_[i] < 0
                # Acid (O12)
                atomtypes[i] = :O12
            elseif bondorder_[bond] == 2 || charge_[i] < 0
                # Carbonyl or oxide (O5-11)
                if atomtypes[nbr] == :C25
                    atomtypes[i] = :O8 # Aromatic carbonyl
                elseif atomtypes[nbr] == :C5
                    # Aliphatic carbonyl
                    cnbrs = adjacencies(mol, nbr)
                    pop!(cnbrs, i)
                    cnbrs = collect(cnbrs)
                    if all(atomsymbol_[cnbrs] .!= :C)
                        atomtypes[i] = :O11 # Carbonyl heteroatom
                    elseif any(isaromatic_[cnbrs] .== true)
                        atomtypes[i] = :O10 # Carbonyl aromatic
                    else
                        atomtypes[i] = :O9 # Carbonyl aliphatic
                    end
                elseif atomsymbol_[nbr] in (:O, :N)
                    atomtypes[i] = :O5 # O2? or N-oxide
                elseif atomsymbol_[nbr] == :S
                    atomtypes[i] = :O6 # S-oxide
                else
                    atomtypes[i] = :O7 # Other oxide
                end
            else
                # Alcohol (O2)
                atomtypes[i] = :O2
            end
        elseif nodedegree_[i] == 2
            # Ether (O3,4)
            if all(isaromatic_[nbrs] .== false)
                atomtypes[i] = :O3 # Aliphatic ether
            else
                atomtypes[i] = :O4 # Aromatic ether
            end
        end
    end

    # Others
    for i in findall(atomtypes .== :undef)
        if atomsymbol_[i] in (:F, :Cl, :Br, :I)
            if charge_[i] == 0
                atomtypes[i] = atomsymbol_[i]
            else
                atomtypes[i] = :Hal
            end
        elseif atomsymbol_[i] == :P
            atomtypes[i] = :P
        elseif atomsymbol_[i] == :S
            if isaromatic_[i]
                atomtypes[i] = :S3
            elseif charge_[i] == 0
                atomtypes[i] = :S1
            else
                atomtypes[i] = :S2
            end
        elseif atomsymbol_[i] == :H
            nbr = iterate(adjacencies(mol, i))[1]
            atomtypes[i] = wclogphydrogentype(mol, nbr)
        elseif atomsymbol_[i] in P_BLOCK
            atomtypes[i] = :Me1
        elseif atomsymbol_[i] in D_BLOCK
            atomtypes[i] = :Me2
        end
    end
    return atomtypes
end

function wclogphydrogentype(mol::GraphMol, i)
    atomsymbol_ = atomsymbol(mol)
    pielectron_ = pielectron(mol)
    isaromatic_ = isaromatic(mol)
    if atomsymbol_[i] == :C
        return :H1 # Hydrocarbon
    elseif atomsymbol_[i] == :N
        return :H3 # Amine
    elseif atomsymbol_[i] == :O
        nbr = iterate(adjacencies(mol, i))[1]
        if atomsymbol_[nbr] == :N
            return :H3 # Hydroxyamine
        elseif atomsymbol_[nbr] in (:O, :S)
            return :H4 # Peroxide, sulfoxide
        elseif atomsymbol_[nbr] == :C
            if pielectron_[nbr] == 1 && !isaromatic_[nbr]
                return :H4 # Acid
            else
                return :H2 # Alcohol
            end
        else
            return :H2 # Other hydroxyl
        end
    else
        return :H2 # Other hydrate
    end
    return :HS # HS may be impossible
end

function wclogpcontrib(mol::GraphMol)
    contrib = zeros(nodecount(mol))
    atomsymbol_ = atomsymbol(mol)
    hcount_ = hcount(mol)
    atomtypes = wclogptype(mol)
    for i in 1:nodecount(mol)
        hcnt = hcount_[i]
        cont = get(WCLOGP_TABLE, string(atomtypes[i]), 0)
        if atomsymbol_[i] == :H
            contrib[i] = 0  # avoid double count
        elseif hcnt == 0
            contrib[i] = cont
        else
            htype = wclogphydrogentype(mol, i)
            hcont = WCLOGP_TABLE[string(htype)] * hcnt
            contrib[i] = cont + hcont
        end
    end
    return contrib
end


text:
timeout(sim::Simulation, delay::Number=0; priority::Bool=false, value::Any=nothing) = schedule(sim, Event(), delay, priority=priority, value=value)
timeout(sim::Simulation, delay::Period; priority::Bool=false, value::Any=nothing) = schedule(sim, Event(), delay, priority=priority, value=value)


text:
using DataFrames
using DataFramesMeta
using CSV

blastfile = ARGS[1]
dadatax = ARGS[2]
outstem = ARGS[3]

blast_fields = "qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qlen"
dada_fields = "Contig Kingdom Phylum Class Order Family Genus"
df = CSV.File(blastfile,header = Array{String,1}(split(blast_fields)))
df2 = CSV.File(dadatax,header = Array{String,1}(split(dada_fields)), skipto = 2)
df = DataFrame(df)
df2 = DataFrame(df2)
#ensure that contig names are strings
df.qseqid = string.(df.qseqid)
df2.Contig = string.(df2.Contig)
tt = getindex.(split.(df2.Contig," "),1)

df2.Contig = tt
contig_genus_map = Dict{String,String}()
contigs = collect(df2.Contig)
genus = df2.Genus
for i in 1:length(contigs)
    contig_genus_map[df2.Contig[i]] = df2.Genus[i]
end
sizes = fill(0,size(df)[1])
species = fill("",size(df)[1])
read_name = fill("",size(df)[1])
genus = fill("",size(df)[1])
genus_species = fill("",size(df)[1])
for i in 1:size(df)[1]
    if occursin("=",df.qseqid[i]) && occursin(";",df.qseqid[i])
        sizes[i] = parse(Int64,split(split(df.qseqid[i],"=")[2],";")[1]) #swarm attaches ; to the end of name
    else
        sizes[i] = 1
    end
    read_name[i] = split(df.qseqid[i],";")[1]

    species[i] = String(split(df.sseqid[i],":")[1])

    genus[i] = contig_genus_map[df.qseqid[i]]
    genus_species[i] = replace(genus[i]*"||"*species[i],"/"=>"-")
end
df.sizes = sizes
df.species = species
df.read_name = read_name
df.genus = genus
df.genus_species = genus_species

dfout = DataFrame(Genus_Species = Array{String,1}(),
    Maximum_length = Array{Int64,1}(),
    Number_of_sequences = Array{Int64,1}(),
    Maximum_abundance = Array{Int64,1}(), #parses size argument in name, not the real number od sequences
    Total_abundance = Array{Int64,1}(), #parses size argument in name, not the real number od sequences
    Samples = Array{String,1}())
for gs in unique(df.genus_species)
    dfs = @linq df |>
               where(:genus_species .== gs) |>
               orderby(-:sizes)
    dfs_ids = @linq dfs |>
        select(:read_name)
    out = outstem*"_ids_of_"*replace(gs,"/"=>"|")
    CSV.write(out,dfs_ids,header=[""])
    maxl = maximum(dfs.qlen)
    maxnb = maximum(dfs.sizes)
    totalnb = sum(dfs.sizes)
    smp = split(outstem,'/')[end]
    push!(dfout,[gs,maxl,size(dfs)[1],maxnb, totalnb, smp])
end
out = outstem*"_blast_summary.tsv"
CSV.write(out,dfout,delim='\t')
#group by genus
dfout = DataFrame(Genus = Array{String,1}(),
    Maximum_length = Array{Int64,1}(),
    Number_of_sequences = Array{Int64,1}(),
    Maximum_abundance = Array{Int64,1}(), #parses size argument in name, not the real number od sequences
    Total_abundance = Array{Int64,1}(), #parses size argument in name, not the real number od sequences
    Samples = Array{String,1}())
for g in unique(df.genus)
    dfs = @linq df |>
               where(:genus .== g) |>
               orderby(-:sizes)
    dfs_ids = @linq dfs |>
        select(:read_name)
    out = outstem*"_GENUS_ids_of_"*replace(g,"/"=>"|")
    CSV.write(out,dfs_ids,header=[""])
    maxl = maximum(dfs.qlen)
    maxnb = maximum(dfs.sizes)
    totalnb = sum(dfs.sizes)
    smp = split(outstem,'/')[end]
    push!(dfout,[g, maxl, size(dfs)[1], maxnb, totalnb, smp])
end
out = outstem*"_blast_summary_genus.tsv"
CSV.write(out,dfout,delim='\t')


text:
using Flux
using Transformers.Basic

struct Encoder
    layers::Array{Transformer}
end

Flux.@functor Encoder
Flux.trainable(encoder::Encoder) = (encoder.layers,)

function (encoder::Encoder)(embedding)
    for transformer ∈ encoder.layers
        embedding = transformer(embedding)
    end
    return embedding
end

function Encoder(
    size::Integer,
    head::Integer,
    hs::Integer,
    ps::Integer,
    layer::Integer;
    act=relu,
    pdrop=0.1
)::Encoder
    return Encoder([
        Transformer(size, head, hs, ps, act=act, pdrop=pdrop)
        for i ∈ 1:layer
    ])
end


text:
using Gadfly

function testMultipleStarts(tsize::Int=300,npoints::Int=20)
    SpeedImage = ones(tsize,tsize)
    SourcePoint = rand(2,npoints).*tsize.+1
    FastMarching.msfm(SpeedImage,SourcePoint,true,true)
end  # function testMultipleStarts

t1 = testMultipleStarts()
draw(PNG("multiplestart.png"),plot(z=t1,Geom.contour))


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    CrossValidation(k, shuffle)

Compare estimation solvers using k-fold cross validation.

The result of the comparison stores the errors for each
variable of the problem.

## Parameters

* k       - number of folds for cross-validation
* shuffle - whether or not to shuffle the data

## Examples

Compare `solver₁` and `solver₂` on a `problem` with variable
`:var` using 10 folds. Plot error distribution:

```julia
julia> result = compare([solver₁, solver₂], problem, CrossValidation(10))
julia> plot(result, bins=50)
```
"""
struct CrossValidation <: AbstractEstimSolverComparison
  k::Int
  shuffle::Bool

  function CrossValidation(k::Int, shuffle::Bool)
    @assert k > 1 "number of folds must be greater than 1"
    new(k, shuffle)
  end
end

CrossValidation(k) = CrossValidation(k, true)
CrossValidation() = CrossValidation(10)

"""
    CrossValidationResult(errors, solvernames)

Stores the cross-validation errors for each variable of the
problem and the solver names used in the comparison.
"""
struct CrossValidationResult
  errors4var::Dict{Symbol,Vector}
  solvernames::Vector{String}
end

function compare(solvers::AbstractVector{S}, problem::EstimationProblem,
                 cmp::CrossValidation) where {S<:AbstractEstimationSolver}
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)
  pmapper = mapper(problem)

  nfolds = cmp.k

  # save results in a dictionary
  errors4var = Dict{Symbol,Vector}()

  for (var,V) in variables(problem)
    # mappings from domain to data locations
    varmap = datamap(problem, var)
    domlocs = collect(keys(varmap))
    datlocs = collect(values(varmap))

    # number of points for variable
    npts = length(varmap)

    # points per fold
    m = npts ÷ nfolds

    @assert nfolds ≤ npts "number of folds must be smaller or equal to number of points"

    # shuffle points if necessary
    perm = cmp.shuffle ? shuffle(1:npts) : sortperm(datlocs)
    domlocs = domlocs[perm]
    datlocs = datlocs[perm]

    # validation errors for each solver
    errors = [V[] for s in 1:length(solvers)]

    # k-fold validation loop
    for i in 0:nfolds-1
      iₛ, iₑ = i*m + 1, (i+1)*m

      # holdout set
      domhold = domlocs[iₛ:iₑ]
      dathold = datlocs[iₛ:iₑ]

      # training set
      train = vcat(datlocs[1:iₛ-1], datlocs[iₑ+1:end])

      subproblem = EstimationProblem(view(pdata, train), pdomain, var, mapper=pmapper)

      if nworkers() > 1
        # run solvers in parallel
        λ = solver -> solve(subproblem, solver)
        solutions = pmap(λ, solvers)
      else
        # fallback to serial execution
        solutions = [solve(subproblem, solver) for solver in solvers]
      end

      for (s, solution) in enumerate(solutions)
        # get solver estimates at holdout locations
        estimates = [solution.mean[var][j] for j in domhold]

        # get true holdout values
        observations = [value(pdata, j, var) for j in dathold]

        # save error and continue
        append!(errors[s], estimates - observations)
      end
    end

    push!(errors4var, var => errors)
  end

  CrossValidationResult(errors4var, [string(s) for s in solvers])
end

@recipe function f(result::CrossValidationResult; labels=nothing)
  # retrieve solver names
  solvers = result.solvernames
  nsolvers = length(solvers)

  # provide default names if labels are not specified
  labels == nothing && (labels = solvers)

  for (var, errs) in result.errors4var
    layout := (nsolvers, 1)
    link := :x
    for (s, err) in enumerate(errs)
      @series begin
        subplot := s
        seriestype := :histogram
        label := labels[s]
        seriescolor --> :black
        fillalpha --> 0.7
        linecolor --> :white
        err
      end
      @series begin
        subplot := s
        seriestype := :vline
        label := "mean error"
        linecolor := :green
        linestyle := :dash
        linewidth := 2
        [sum(err)/length(err)]
      end
      @series begin
        subplot := s
        seriestype := :vline
        label := "interquartile"
        linecolor := :orange
        linestyle := :dot
        linewidth := 2
        quantile(err, [.25,.75])
      end
      xlabel --> "errors = estimates - observations"
      ylabel --> "histogram"
    end
  end
end


text:
# Modules represent the top-level structure in an LLVM program.

export dispose,
       name, name!,
       triple, triple!,
       datalayout, datalayout!,
       context, inline_asm!,
       set_used!, set_compiler_used!

# forward definition of Module in src/core/value/constant.jl
reftype(::Type{Module}) = API.LLVMModuleRef

# forward declarations
@checked struct DataLayout
    ref::API.LLVMTargetDataRef
end

Module(name::String) = Module(API.LLVMModuleCreateWithName(name))
Module(name::String, ctx::Context) =
    Module(API.LLVMModuleCreateWithNameInContext(name, ref(ctx)))
Module(mod::Module) = Module(API.LLVMCloneModule(ref(mod)))

dispose(mod::Module) = API.LLVMDisposeModule(ref(mod))

function Module(f::Core.Function, args...)
    mod = Module(args...)
    try
        f(mod)
    finally
        dispose(mod)
    end
end

function Base.show(io::IO, mod::Module)
    output = unsafe_message(API.LLVMPrintModuleToString(ref(mod)))
    print(io, output)
end

function name(mod::Module)
    out_len = Ref{Csize_t}()
    ptr = convert(Ptr{UInt8}, API.LLVMGetModuleIdentifier(ref(mod), out_len))
    return unsafe_string(ptr, out_len[])
end
name!(mod::Module, str::String) =
    API.LLVMSetModuleIdentifier(ref(mod), str, Csize_t(length(str)))

triple(mod::Module) = unsafe_string(API.LLVMGetTarget(ref(mod)))
triple!(mod::Module, triple) = API.LLVMSetTarget(ref(mod), triple)

datalayout(mod::Module) = DataLayout(API.LLVMGetModuleDataLayout(ref(mod)))
datalayout!(mod::Module, layout::String) = API.LLVMSetDataLayout(ref(mod), layout)
datalayout!(mod::Module, layout::DataLayout) =
    API.LLVMSetModuleDataLayout(ref(mod), ref(layout))

inline_asm!(mod::Module, asm::String) =
    API.LLVMSetModuleInlineAsm(ref(mod), asm)

context(mod::Module) = Context(API.LLVMGetModuleContext(ref(mod)))

if VERSION >= v"1.5" && !(v"1.6-" <= VERSION < v"1.6.0-DEV.90")

set_used!(mod::Module, values::GlobalVariable...) =
    API.LLVMExtraAppendToUsed(ref(mod), collect(ref.(values)), length(values))

set_compiler_used!(mod::Module, values::GlobalVariable...) =
    API.LLVMExtraAppendToCompilerUsed(ref(mod), collect(ref.(values)), length(values))

else
set_used!(mod::Module, values::GlobalVariable...) = nothing
set_compiler_used!(mod::Module, values::GlobalVariable...) = nothing
end

## type iteration

export types

struct ModuleTypeDict <: AbstractDict{String,LLVMType}
    mod::Module
end

types(mod::Module) = ModuleTypeDict(mod)

function Base.haskey(iter::ModuleTypeDict, name::String)
    return API.LLVMGetTypeByName(ref(iter.mod), name) != C_NULL
end

function Base.getindex(iter::ModuleTypeDict, name::String)
    objref = API.LLVMGetTypeByName(ref(iter.mod), name)
    objref == C_NULL && throw(KeyError(name))
    return LLVMType(objref)
end

## metadata iteration

export metadata

struct ModuleMetadataDict <: AbstractDict{String,Vector{MetadataAsValue}}
    mod::Module
end

metadata(mod::Module) = ModuleMetadataDict(mod)

function Base.haskey(iter::ModuleMetadataDict, name::String)
    return API.LLVMGetNamedMetadataNumOperands(ref(iter.mod), name) != 0
end

function Base.getindex(iter::ModuleMetadataDict, name::String)
    nops = API.LLVMGetNamedMetadataNumOperands(ref(iter.mod), name)
    nops == 0 && throw(KeyError(name))
    ops = Vector{API.LLVMValueRef}(undef, nops)
    API.LLVMGetNamedMetadataOperands(ref(iter.mod), name, ops)
    return MetadataAsValue[MetadataAsValue(op) for op in ops]
end

Base.push!(iter::ModuleMetadataDict, name::String, val::MetadataAsValue) =
    API.LLVMAddNamedMetadataOperand(ref(iter.mod), name, ref(val))

## global variable iteration

export globals

struct ModuleGlobalSet
    mod::Module
end

globals(mod::Module) = ModuleGlobalSet(mod)

Base.eltype(::ModuleGlobalSet) = GlobalVariable

function Base.iterate(iter::ModuleGlobalSet, state=API.LLVMGetFirstGlobal(ref(iter.mod)))
    state == C_NULL ? nothing : (GlobalVariable(state), API.LLVMGetNextGlobal(state))
end

Base.last(iter::ModuleGlobalSet) =
    GlobalVariable(API.LLVMGetLastGlobal(ref(iter.mod)))

Base.isempty(iter::ModuleGlobalSet) =
    API.LLVMGetLastGlobal(ref(iter.mod)) == C_NULL

Base.IteratorSize(::ModuleGlobalSet) = Base.SizeUnknown()

# partial associative interface

function Base.haskey(iter::ModuleGlobalSet, name::String)
    return API.LLVMGetNamedGlobal(ref(iter.mod), name) != C_NULL
end

function Base.getindex(iter::ModuleGlobalSet, name::String)
    objref = API.LLVMGetNamedGlobal(ref(iter.mod), name)
    objref == C_NULL && throw(KeyError(name))
    return GlobalVariable(objref)
end

## function iteration

export functions

struct ModuleFunctionSet
    mod::Module
end

functions(mod::Module) = ModuleFunctionSet(mod)

Base.eltype(::ModuleFunctionSet) = Function

function Base.iterate(iter::ModuleFunctionSet, state=API.LLVMGetFirstFunction(ref(iter.mod)))
    state == C_NULL ? nothing : (Function(state), API.LLVMGetNextFunction(state))
end

Base.last(iter::ModuleFunctionSet) =
    Function(API.LLVMGetLastFunction(ref(iter.mod)))

Base.isempty(iter::ModuleFunctionSet) =
    API.LLVMGetLastFunction(ref(iter.mod)) == C_NULL

Base.IteratorSize(::ModuleFunctionSet) = Base.SizeUnknown()

# partial associative interface

function Base.haskey(iter::ModuleFunctionSet, name::String)
    return API.LLVMGetNamedFunction(ref(iter.mod), name) != C_NULL
end

function Base.getindex(iter::ModuleFunctionSet, name::String)
    objref = API.LLVMGetNamedFunction(ref(iter.mod), name)
    objref == C_NULL && throw(KeyError(name))
    return Function(objref)
end

## module flag iteration
# TODO: doesn't actually iterate, since we can't list the available keys

if libllvm_version >= v"8.0"

export flags

struct ModuleFlagDict <: AbstractDict{String,Metadata}
    mod::Module
end

flags(mod::Module) = ModuleFlagDict(mod)

Base.haskey(iter::ModuleFlagDict, name::String) =
    API.LLVMGetModuleFlag(ref(iter.mod), name, length(name)) != C_NULL

function Base.getindex(iter::ModuleFlagDict, name::String)
    objref = API.LLVMGetModuleFlag(ref(iter.mod), name, length(name))
    objref == C_NULL && throw(KeyError(name))
    return Metadata(objref)
end

function Base.setindex!(iter::ModuleFlagDict, val::Metadata,
                        (name, behavior)::Tuple{String, LLVM.API.LLVMModuleFlagBehavior})
    API.LLVMAddModuleFlag(ref(iter.mod), behavior, name, length(name), ref(val))
end

end


text:
const MSB = 0x80
const MASK7 = 0x7f
const MASK8 = 0xff

const _wfbuf = (Vector{UInt8}(undef, 1), Vector{UInt8}(undef, 2), Vector{UInt8}(undef, 4), Vector{UInt8}(undef, 8), Vector{UInt8}(undef, 16))

const TIO = Union{IO, TTransport}

function _write_fixed(io::TIO, ux::T, bigendian::Bool) where T <: Unsigned
    N = sizeof(ux)
    _write_fixed(io, ux, _wfbuf[Int(log2(N))+1], bigendian ? (N:-1:1) : (1:N))
end

function _write_fixed(io::TIO, ux::T, a::Vector{UInt8}, r::R) where {T <: Unsigned, R <: AbstractRange}
    for n in r
        a[n] = UInt8(ux & MASK8)
        ux >>>= 8
    end
    write(io, a)
end

function _read_fixed(io::TIO, ret::T, N::Int, bigendian::Bool) where T <: Unsigned
    r = bigendian ? ((N-1):-1:0) : (0:1:(N-1))
    for n in r
        byte = convert(T, read(io, UInt8))
        ret |= (byte << (8*n))
    end
    ret
end

function _write_uleb(io::TIO, x::T) where T <: Integer
    nw = 0
    cont = true
    while cont
        byte = x & MASK7
        if (x >>>= 7) != 0
            byte |= MSB
        else
            cont = false
        end
        nw += write(io, UInt8(byte))
    end
    nw
end

function _read_uleb(io::TIO, typ::Type{T}) where T <: Integer
    res = convert(typ, 0)
    n = 0
    byte = UInt8(MSB)
    while (byte & MSB) != 0
        byte = read(io, UInt8)
        res |= (convert(typ, byte & MASK7) << (7*n))
        n += 1
    end
    res
end

function _write_zigzag(io::TIO, x::T) where T <: Integer
    nbits = 8*sizeof(x)
    zx = xor((x << 1), (x >> (nbits-1)))
    _write_uleb(io, zx)
end

function _read_zigzag(io::TIO, typ::Type{T}) where T <: Integer
    zx = _read_uleb(io, UInt64)
    # result is positive if zx is even
    convert(typ, iseven(zx) ? (zx >>> 1) : -((zx+1) >>> 1))
end


text:
#abstract type AbstractSimpleValueEdge{T, U} <: AbstractEdge{T} end

struct SimpleValueEdge{V<:Integer, E_VAL} <: AbstractEdge{V}
    src::V
    dst::V
    value::E_VAL
end

SimpleValueEdge(src, dst) = SimpleValueEdge(src, dst, nothing)

src(e::SimpleValueEdge) = e.src
dst(e::SimpleValueEdge) = e.dst
val(e::SimpleValueEdge) = e.value

reverse(e::SimpleValueEdge) = SimpleValueEdge(dst(e), src(e), val(e))

show(io::IO, e::SimpleValueEdge) = print(io, "Edge $(src(e)) => $(dst(e)) with value $(repr(val(e)))")

Tuple(e::SimpleValueEdge) = (src(e), dst(e), val(e))


text:
using Documenter, OmniSci

makedocs(
    # options
    modules = [OmniSci],
    clean = false,
    format = Documenter.HTML(assets = ["assets/favicon.ico"]),
    sitename = "OmniSci.jl",
    pages = [
        "Getting Started" => "index.md",
        "Functions" => "functions.md",
        "API Reference" => "apireference.md",
        "Contributing / Developer Docs" => "devdocs.md"

    ]
)

deploydocs(
    repo = "github.com/omnisci/OmniSci.jl.git",
    #julia = "1.0",
    deps=nothing,
    make=nothing,
    target = "build"
)


text:
basetype(x) = basetype(typeof(x))
basetype(t::Type) = t.name.wrapper
basetype(t::UnionAll) = t

f <| x = f(x)

# This shouldn't be hard coded, but it makes plots tolerable for now
shorten(x::AbstractFloat) = round(x, sigdigits=3)
shorten(x) = x

# Nothing doesn't string
getstring(::Nothing) = ""
getstring(x) = string(x)


text:
#=================================== SUITS ===================================#
"""
The root abstract type for suit types.
This is also used to represent "pseudo-suit" types, like the winds, dragons, flowers, etc.
"""
abstract type Suit end

abstract type Character <: Suit end
abstract type Bamboo <: Suit end
abstract type Circle <: Suit end
abstract type Wind <: Suit end
abstract type Dragon <: Suit end
abstract type Flower <: Suit end
abstract type Season <: Suit end

@doc """
    is_sequential(suit_type)

Check whether this suit type (or pseudo-suit) is a sequential one (for the purposes of hands).
For example, typically the characters suit is sequential, but the winds suit is not.
"""
is_sequential(::Type{<:Suit}) = false
is_sequential(s::Suit) = is_sequential(typeof(s))

is_sequential(::Type{Character}) = true
is_sequential(::Type{Bamboo}) = true
is_sequential(::Type{Circle}) = true

@doc """
    max_number(suit_type)

The maximum ("index") number for this suit type (or pseudo-suit).
For example, the traditional suited tiles (characters, bamboo, circle)
have a max of 9 (starting at 1),
whereas the winds pseudo-suit has a max of 4.
"""
max_number(::Type{<:Suit}) = 1
max_number(s::Suit) = max_number(typeof(s))

max_number(::Type{Character}) = 9
max_number(::Type{Bamboo}) = 9
max_number(::Type{Circle}) = 9
max_number(::Type{Wind}) = 4
max_number(::Type{Dragon}) = 3
max_number(::Type{Flower}) = 4
max_number(::Type{Season}) = 4

@doc """
    first_char(suit_type)

The character for the first tile in this suit (in Unicode).

## Examples
```julia-repl
julia> first_char(Character)
🀇

```
"""
first_char(::Type{<:Suit}) = throw("not implemented")
first_char(::Type{Character}) = '🀇'
first_char(::Type{Bamboo}) = '🀐'
first_char(::Type{Circle}) = '🀙'
first_char(::Type{Wind}) = '🀀'
first_char(::Type{Dragon}) = '🀄'
first_char(::Type{Flower}) = '🀢'
first_char(::Type{Season}) = '🀦'

@doc """
    typical_duplication(suit_type)

The number of duplicate tiles usually used for this suit type
(this is usually 1 or 4).
"""
typical_duplication(::Type{<:Suit}) = 1
typical_duplication(::Type{Character}) = 4
typical_duplication(::Type{Bamboo}) = 4
typical_duplication(::Type{Circle}) = 4
typical_duplication(::Type{Wind}) = 4
typical_duplication(::Type{Dragon}) = 4

Base.isless(x::Type{<:Suit}, y::Type{<:Suit}) = first_char(x) < first_char(y)

#=================================== TILES ===================================#

"""
The root tile type for Mahjong tiles.
However, the contract for what the Tile interface entails is not well-defined,
so most functions dispatch on the `SuitedTile` and not on `Tile`.
"""
abstract type Tile end

"""
A tile with a suit (or pseudo-suit).
"""
struct SuitedTile{SuitType <: Suit} <: Tile
    number::Int8

    """
        SuitedTile(suit_type, number=1)

    Create a tile of a particular suit `suit_type` and `number`
    (use the default `1` for a "suit" that has only 1 tile).
    """
    function SuitedTile(SuitType::Type{<:Suit}, number::Integer = 1)
        number < 1 && throw("number must be positive")
        number > max_number(SuitType) && throw("number is greater than the max allowed for this suit")

        return new{SuitType}(number)
    end
end

suit(::SuitedTile{T}) where T = T
number(tile::SuitedTile) = tile.number
is_suit(::SuitedTile{T}, suit_type::Type{<:Suit}) where T = T == suit_type
is_number(tile::SuitedTile, num::Integer) = tile.number == num

function Base.:(==)(x::SuitedTile{T}, y::SuitedTile{S}) where {T, S}
    return T == S && x.number == y.number
end
Base.hash(tile::SuitedTile{SuitType}, h::UInt) where SuitType = hash(SuitType, hash(tile.number, h))

function Base.isless(left::SuitedTile{LeftSuit}, right::SuitedTile{RightSuit}) where {LeftSuit, RightSuit}
    if LeftSuit != RightSuit
        return isless(LeftSuit, RightSuit)
    else
        return left.number < right.number
    end
end

function Base.show(io::IO, ::MIME"text/plain", tile::SuitedTile{T}) where T
    print(io, first_char(T) + tile.number - 1)
end
function Base.show(io::IO, tile::SuitedTile{T}) where T
    show(io::IO, MIME"text/plain"(), tile)
end

#============================== TILE CONSTANTS ===============================#

"""
List of the "normal" suits: character, bamboo, circles.

```@repl
using MahjongTiles

MahjongTiles.suits
```
"""
const suits = [Character, Bamboo, Circle]

"""
List of all suits, including "pseudo-suits" like the winds and flowers.
"""
const all_suits = [suits..., Wind, Dragon, Flower, Season]

for suit in [Character, Bamboo, Circle, Wind, Dragon, Flower, Season], i in 1:max_number(suit)
    tile_name = Symbol(first_char(suit) + i - 1)
    @eval MahjongTiles, const $tile_name = SuitedTile($suit, $i)
    @eval MahjongTiles, export $tile_name
end

@doc """
    character(index)

Creates a tile of the character suit with the particular number.

## Example

```jldoctest
julia> character(5)
🀋
```

"""
character(index::Integer) = SuitedTile(Character, index)

@doc """
    bamboo(index)

Creates a tile of the bamboo suit with the particular number.

## Example

```jldoctest
julia> bamboo(2)
🀑
```
"""
bamboo(index::Integer) = SuitedTile(Bamboo, index)

@doc """
    circle(index)

Creates a tile of the circle suit with the particular number.

## Example

```jldoctest
julia> circle(9)
🀡
```
"""
circle(index::Integer) = SuitedTile(Circle, index)

@doc """
    wind(index)
    wind(direction)

Creates a wind tile.
If an integer is used,
the standard Chinese ordering of the cardinal directions is used
(i.e. `1` -> east, `2` -> south, `3` -> west, `4` -> north).
Otherwise, symbols can also be used for the names of the directions
(e.g. `:east`, `:north`).

## Example

```jldoctest
julia> wind(:north)
🀃

julia> wind(2)
🀁
```
"""
wind(index::Integer) = SuitedTile(Wind, index)
function wind(name::Symbol)
    if name === :east
        return SuitedTile(Wind, 1)
    elseif name === :south
        return SuitedTile(Wind, 2)
    elseif name === :west
        return SuitedTile(Wind, 3)
    elseif name === :north
        return SuitedTile(Wind, 4)
    else
        throw("unrecognized cardinal direction")
    end
end

@doc """
    dragon(index)
    dragon(name)

Creates a wind tile.
If an integer is used,
the standard ordering is used
(i.e. `1` -> red dragon, `2` -> green dragon, `3` -> white dragon).
Otherwise, symbols can also be used for their names/colors
(e.g. `:red`, `:green`).

## Example

```jldoctest
julia> dragon(:red)
🀄

julia> dragon(2)
🀅
```
"""
dragon(index::Integer) = SuitedTile(Dragon, index)
function dragon(name::Symbol)
    if name === :red
        return SuitedTile(Dragon, 1)
    elseif name === :green
        return SuitedTile(Dragon, 2)
    elseif name === :white
        return SuitedTile(Dragon, 3)
    else
        throw("unrecognized dragon color")
    end
end

@doc """
    flower(index)

Creates a flower tile with the particular number.

!!! note

    While this function (and the `season` function) don't allow using the flower (season) names,
    they may have this functionality in the future.

## Example

```jldoctest
julia> flower(2)
🀣
```
"""
flower(index::Integer) = SuitedTile(Flower, index)

@doc """
    season(index)

Creates a season (flower) tile with the particular number.

!!! note

    While this function (and the `flower` function) don't allow using the season (flower) names,
    they may have this functionality in the future.

## Example

```jldoctest
julia> season(4)
🀩
```
"""
season(index::Integer) = SuitedTile(Season, index)


text:
using SatelliteToolbox
using Plots
gr()

Lat = -π/2:10/180*π:π/2
x = zeros(18,3)
println(x)
anim = @animate for i=1:18
    println(i)
    x[i,:] = GeodetictoECEF(Lat[i],0,4000000)
    println(x[i,:])
    plot(x[1:i,1],x[1:i,3])
end
gif(anim, "test6_anim.gif", fps = 5)

#y = ECEFtoGeodetic([6367000, 0, -4000000])
#println(y)

#=
Lat = 30/180*π
Lon = 120/180*π
h = 4000000

sin_lat, cos_lat = sincos(Lat)
sin_lon, cos_lon = sincos(Lon)

a_wgs84  = 6378137.0
f_wgs84  = 1/298.257223563
b_wgs84  = a_wgs84*(1-f_wgs84)

N =a_wgs84/sqrt(1 - e_wgs84^2*sin_lat^2 )

x = [(N+h)*cos_lat*cos_lon, (N+h)*cos_lat*sin_lon, ((b_wgs84/a_wgs84)^2*N+h)*sin_lat]
println(x)
=#

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
import SAMPLER
SAMPLER.main()

max_stars_repo_name: SeisSol/SAMPLER

max_issues_repo_name: SeisSol/SAMPLER

max_forks_repo_name: SeisSol/SAMPLER


text:
module Quine
    export generate_quine
    function escape(string)
        string = replace(string, "\\" => "\\\\")
        return replace(string, "\"" => "\\\"")
    end

    function generate_quine(n)
        variable_names = collect('a':'z')
        lines = []
        if n > 26
            println("sorry only supports 26 variables max because of variable name")
            println("...and because the author's lazy :p")
        elseif n == 1
            println("a = :(println(\"a = :(\" * string(a) * \"); eval(a)\"); eval(a)")
            return
        end
        push!(lines, "a = :(eval(b))")
        for i in 2:n - 1
            previous_string = escape(lines[i - 1])
            push!(lines, "$(variable_names[i]) = :(eval($(variable_names[i + 1])); println(\"$(previous_string)\"))")
        end
        last_string = escape(lines[n - 1])
        push!(lines, "$(variable_names[n]) = :(println(\"$(variable_names[n]) = :(\" * string($(variable_names[n])) * \")\\n$(last_string)\"))")
        for i in n:-1:1
            println(lines[i])
        end
        println("eval(a)")
    end
end


text:
using TestPkgTemplates
using Test

@testset "TestPkgTemplates.jl" begin
    # Write your tests here.
    @test testadd(1,2) == 3
end


text:
@testset "readTopology Hybrid Node isMajor Tests" begin
	# n1 previously failed (issue #44): both hybrid parent edges were minor
	n1 = readTopology("(A,((B,#H1),(C,(D)#H1:::0.2)));")
	@test writeTopology(n1) == "(A,((B,(D)#H1:::0.8),(C,#H1:::0.2)));"
	hybridParents = [x for x in n1.hybrid[1].edge if x.hybrid]
	# Exclusive or to check exactly one parent edge is major
	@test hybridParents[1].isMajor ⊻ hybridParents[2].isMajor
	# Gammas must sum to 1 and be non-negative for hybrids
	@test hybridParents[1].gamma + hybridParents[2].gamma == 1.0
	@test hybridParents[1].gamma >= 0.0
	@test hybridParents[2].gamma >= 0.0
	# n2 previously failed: both hybrid parent edges were major
	n2 = readTopology("(A,((B,#H1:::0.5),(C,(D)#H1:::0.5)));")
	@test writeTopology(n2) == "(A,((B,#H1:::0.5),(C,(D)#H1:::0.5)));"
	hybridParents = [x for x in n2.hybrid[1].edge if x.hybrid]
	# Exclusive or to check exactly one parent edge is major
	@test hybridParents[1].isMajor ⊻ hybridParents[2].isMajor
	# Gammas must sum to 1 and be non-negative for hybrids
	@test hybridParents[1].gamma + hybridParents[2].gamma == 1.0
	@test hybridParents[1].gamma >= 0.0
	@test hybridParents[2].gamma >= 0.0
	n3 = readTopology("(A,((C,(D)#H1:::0.5),(B,#H1:::0.5)));")
	@test writeTopology(n3) =="(A,((C,(D)#H1:::0.5),(B,#H1:::0.5)));"
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LLVM_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LLVM")
JLLWrappers.@generate_main_file("LLVM", UUID("86de99a1-58d6-5da7-8064-bd56ce2e322c"))
end  # module LLVM_jll


text:
function insertion_sort!(A)
    for j in 2:length(A)
        key = A[j]

        i = j - 1
        while i > 0 && A[i] > key
            A[i+1] = A[i]
            i -= 1
        end
        A[i+1] = key
    end
end

max_stars_repo_name: blallen/CodePractice

max_issues_repo_name: blallen/CodePractice

max_forks_repo_name: blallen/CodePractice


text:
function single_task_components(xm, vm, task::Task{<:BaseTaskMap}, CM::Chart{I,M},
        CN::Chart{J,N}) where {M,N,I,J}
    m, n = dim(M), dim(N)

    σ, σdot = xm, vm
    xn = task_map_chart(σ, task, CM, CN)
    if active_weight_position_chart(xn, task, CN)
        Jf = task_jacobian_chart(σ, task, CM, CN)
        vn = Jf*vm
        W = weight_metric_chart(xn, vn, task, CN)
        if W != @SMatrix zeros(n,n)
            Jfdot = task_jacobian_chart_dot(σ, σdot, task, CM, CN)

            Γ = christoffel_symbols(xn, task, CN)
            if any(isinf.(Γ)) || any(isnan.(Γ))
                Γ = eltype(xm).(christoffel_symbols(BigFloat.(xn), task, CN))
            end
            g = metric_chart(xn, task, CN)
            ginv = inv(g)

            ℱ_pot = potential_force_chart(xn, task, CN)
            ℱ_dis = dissipative_forces_chart(xn, vn, task, CN)
            ℱ = ℱ_pot + ℱ_dis

            m_inds, n_inds = static(1):static(m), static(1):static(n)
            Ξ = SMatrix{n,m,eltype(xm)}([sum(Tuple(Jf[l,j]*Γ[l,h,k]*Jf[h,r]*σdot[r]
                for l=n_inds, h=n_inds, r=m_inds)) for k=n_inds, j=m_inds])

            JftWJf = Jf'*W*Jf
            JftWA = Jf'*W*(ginv*ℱ - (Jfdot + Ξ)*σdot)
        else
            JftWJf = @SMatrix zeros(m,m)
            JftWA = @SVector zeros(m)
        end
    else
        JftWJf = @SMatrix zeros(m,m)
        JftWA = @SVector zeros(m)
    end

    JftWJf, JftWA
end

function single_task_acceleration(xm, vm, task::Task{<:BaseTaskMap}, CM::Chart{I,M},
        CN::Chart{J,N}, robot_coord_rep=ChartRep()) where {M,N,I,J}
    m, n = dim(M), dim(N)
    (xm, vm) = robot_coord_rep == EmbRep() ? emb_to_chart_differential(xm, vm, CM) : (xm, vm)
    CN = isglobal(CN) ? CN : choose_chart_chart(xm, task, CM, CN)
    JftWJf, JftWA = single_task_components(xm, vm, task, CM, CN)
    σxddot = SMatrix{m,m,eltype(xm)}(pinv(Matrix(JftWJf)))*JftWA
    if robot_coord_rep == EmbRep()
        xme, vme, σxddot = chart_to_emb_differential(xm, vm, σxddot, CM)
    end
    σxddot, CN
end

function multiple_task_acceleration(xm, vm, tasks::TaskList, CM::Chart{I,M}, CNs::ChartList,
        robot_coord_rep=ChartRep(); log_task_chart=false) where {M,I}
    m = dim(M)
    JftWJf_sum = zeros(m,m)
    JftWA_sum = zeros(m)
    CNs_out = ChartList()

    (xm, vm) = robot_coord_rep == EmbRep() ? emb_to_chart_differential(xm, vm, CM) : (xm, vm)
    for i in 1:length(tasks)
        CN = isglobal(CNs[i]) ? CNs[i] : choose_chart_chart(xm, tasks[i], CM, CNs[i])
        JftWJf, JftWA = single_task_components(xm, vm, tasks[i], CM, CN)
        JftWJf_sum += JftWJf
        JftWA_sum += JftWA
        log_task_chart && push!(CNs_out, CN)
    end

    σxddot = SMatrix{m,m,eltype(xm)}(pinv(Matrix(JftWJf_sum)))*JftWA_sum
    if robot_coord_rep == EmbRep()
        xme, vme, σxddot = chart_to_emb_differential(xm, vm, σxddot, CM)
    end
    σxddot, CNs_out
end

function task_map_pullback!(JftWJfm_sum, Am_sum, Bm_sum, JftWginvℱm_sum, Γm_sum, xm, vm,
        node::TreeNode{<:BaseTaskMap}, CM::Chart{I,M}, CN::Chart{J,N};
        log_tasks=false) where {M,N,I,J}
    m, n = dim(M), dim(N)
    JftWJfn_sum, An_sum, Bn_sum = zeros(n,n), zeros(n,n), zeros(n,n)
    JftWginvℱn_sum = zeros(n)
    Γn_sum = zeros(n,n,n)

    task_map = node.data

    Jf = task_jacobian_chart(xm, task_map, CM, CN)
    xn = task_map_chart(xm, task_map, CM, CN)
    vn = Jf*vm

    for child in children(node)
        if !isglobal(child.chart)
            CN_child = choose_chart_chart(xn, child.data, CN, child.chart)
        else
            CN_child = child.chart
        end
        task_map_pullback!(JftWJfn_sum, An_sum, Bn_sum, JftWginvℱn_sum, Γn_sum,
            xn, vn, child, CN, CN_child; log_tasks)
    end

    if Bn_sum != @SMatrix zeros(n,n)
        Jfdot = task_jacobian_chart_dot(xm, vm, task_map, CM, CN)
        Am = Jf'*(An_sum*Jf + Bn_sum*Jfdot)
        Bm = Jf'*Bn_sum*Jf
        Bm_sum .+= Bm
        log_tasks && push!(node.traj_log.B, Bm)
    else
        Am = Jf'*An_sum*Jf
    end
    Am_sum .+= Am

    if JftWJfn_sum != @SMatrix zeros(n,n)
        JftWJfm = Jf'*JftWJfn_sum*Jf
        JftWJfm_sum .+= JftWJfm
        log_tasks && push!(node.traj_log.JftWJf, JftWJfm)
        if JftWginvℱn_sum != @SVector zeros(n)
            JftWginvℱm = Jf'*JftWginvℱn_sum
            JftWginvℱm_sum .+= JftWginvℱm
            log_tasks && push!(node.traj_log.JftWginvℱ, JftWginvℱm)
        end
    end

    if Γn_sum != @SArray zeros(n,n,n)
        @tullio Γmnn[i,h,k] := Jf[l,i]*Γn_sum[l,h,k]
        @tullio Γmmn[l,s,k] := Jf[h,s]*Γmnn[l,h,k]
        @tullio Γm[l,h,q] := Jf[k,q]*Γmmn[l,h,k]
        Γm_sum .+= Γm
        log_tasks && push!(node.traj_log.Γ, Γm)
    end

    if log_tasks
        if node.coord_rep == ChartRep()
            push!(node.traj_log.x, xn)
            push!(node.traj_log.v, vn)
        else
            xne, vne = chart_to_emb_differential(xn, vn, CN)
            push!(node.traj_log.x, xne)
            push!(node.traj_log.v, vne)
        end
        push!(node.traj_log.chart, CN)
        push!(node.traj_log.A, Am)
    end

    nothing
end

function task_components!(JftWJfm_sum, Am_sum, Bm_sum, JftWginvℱm_sum, Γm_sum, xm, vm,
        node::TreeNode{<:Task{<:BaseTaskMap}}, CM::Chart{I,M}, CN::Chart{J,N};
        log_tasks=false) where {M,N,I,J}
    m, n = dim(M), dim(N)
    task = node.data

    xn = task_map_chart(xm, task, CM, CN)
    if active_weight_position_chart(xn, task, CN) || log_tasks
        Jf = task_jacobian_chart(xm, task, CM, CN)
        vn = Jf*vm
        W = weight_metric_chart(xn, vn, task, CN)

        if log_tasks
            if node.coord_rep == ChartRep()
                push!(node.traj_log.x, xn)
                push!(node.traj_log.v, vn)
            else
                xne, vne = chart_to_emb_differential(xn, vn, CN)
                push!(node.traj_log.x, xne)
                push!(node.traj_log.v, vne)
            end
            push!(node.traj_log.chart, CN)
        end

        if W != @SMatrix zeros(n,n)
            JftW = Jf'*W
            Jfdot = task_jacobian_chart_dot(xm, vm, task, CM, CN)
            Γn = christoffel_symbols(xn, task, CN)
            if any(isinf.(Γn)) || any(isnan.(Γn))
                Γn = eltype(xm).(christoffel_symbols(BigFloat.(xn), task, CN))
            end
            g = metric_chart(xn, task, CN)

            ℱ_pot = potential_force_chart(xn, task, CN)
            ℱ_dis = dissipative_forces_chart(xn, vn, task, CN)
            ℱ = ℱ_pot + ℱ_dis
            JftWJf = Jf'*W*Jf
            JftWJfm_sum .+= JftWJf
            Bm = JftWJf
            Bm_sum .+= Bm

            if ℱ != @SVector zeros(n)
                JftWginvℱm = Jf'*W*inv(g)*ℱ
                JftWginvℱm_sum .+= JftWginvℱm
                log_tasks && push!(node.traj_log.JftWginvℱ, JftWginvℱm)
            end
            if Γn != @SArray zeros(n,n,n)
                @tullio Γmnn[i,h,k] := Jf[l,i]*Γn[l,h,k]
                @tullio Γmmn[l,s,k] := Jf[h,s]*Γmnn[l,h,k]
                @tullio Γm[l,h,q] := JftW[q,k]*Γmmn[l,h,k]
                Γm_sum .+= Γm
                log_tasks && push!(node.traj_log.Γ, Γm)
            end

            Am = Jf'*W*Jfdot
            Am_sum .+= Am

            if log_tasks
                push!(node.traj_log.g, g)
                push!(node.traj_log.ginv, inv(g))
                push!(node.traj_log.JftWJf, JftWJf)
                push!(node.traj_log.A, Am)
                push!(node.traj_log.B, Bm)
            end

        elseif log_tasks
            push!(node.traj_log.g, 0)
            push!(node.traj_log.ginv, 0)
            push!(node.traj_log.Γ, 0)
            push!(node.traj_log.JftWJf, 0)
            push!(node.traj_log.JftWginvℱ, 0)
            push!(node.traj_log.A, 0)
            push!(node.traj_log.B, 0)
        end
    end

    nothing
end

function task_map_pullback!(JftWJfm_sum, Am_sum, Bm_sum, JftWginvℱm_sum, Γm_sum, xm, vm,
        node::TreeNode{<:TaskX}, CM, CN; log_tasks=false)
    task_components!(JftWJfm_sum, Am_sum, Bm_sum, JftWginvℱm_sum, Γm_sum, xm, vm,
        node, CM, CN; log_tasks)
end

function task_acceleration(xm, vm, node::TreeNode{M}, CM::Chart{I,M}, robot_coord_rep=ChartRep();
        log_tasks=false) where {M<:Manifold,I}
    m = dim(M)
    SVm, SMmm = SVector{m,eltype(xm)}, SMatrix{m,m,eltype(xm)}
    JftWJf_sum, A_sum, B_sum = zeros(m,m), zeros(m,m), zeros(m,m)
    JftWginvℱ_sum = zeros(m)
    Γ_sum = zeros(m,m,m)

    (xm, vm) = robot_coord_rep == EmbRep() ? emb_to_chart_differential(xm, vm, CM) : (xm, vm)
    for child in children(node)
        if !isglobal(child.chart)
            CN_child = choose_chart_chart(xm, child.data, CM, child.chart)
        else
            CN_child = child.chart
        end
        task_map_pullback!(JftWJf_sum, A_sum, B_sum, JftWginvℱ_sum, Γ_sum,
            xm, vm, child, CM, CN_child; log_tasks)
    end

    m_inds = static(1):static(m)
    vmΓm = SMmm([sum(Tuple(vm[l]*Γ_sum[l,h,k] for l=m_inds)) for k=m_inds, h=m_inds])
    vmΓmvm = vmΓm*vm

    JftWginvℱ_sum_st = SVm(JftWginvℱ_sum)
    A_sum_st = SMmm(A_sum)

    if log_tasks
        push!(node.traj_log.x, xm)
        push!(node.traj_log.v, vm)
        push!(node.traj_log.chart, CM)
        push!(node.traj_log.JftWJf, JftWJf_sum)
        push!(node.traj_log.JftWginvℱ, JftWginvℱ_sum_st)
        push!(node.traj_log.A, A_sum_st)
        push!(node.traj_log.vmΓmvm, vmΓmvm)
    end

    σxddot = SMmm(pinv(JftWJf_sum))*(JftWginvℱ_sum_st - A_sum_st*vm - vmΓmvm)
    if robot_coord_rep == EmbRep()
        xme, vme, σxddot = chart_to_emb_differential(xm, vm, σxddot, CM)
    end
    σxddot
end

# Version where mechanism state needs to be updated
function task_acceleration(xm, vm, node::TreeNode{M}, CM::Chart{I,M}, state::MechanismState,
    robot_coord_rep=ChartRep(); log_tasks=false) where {M<:Manifold,I}
    set_configuration!(state, xm)
    task_acceleration(xm, vm, node, CM, robot_coord_rep; log_tasks)
end

# Version where cache is also present
function task_acceleration(xm, vm, node::TreeNode{M}, CM::Chart{I,M}, state::MechanismState,
        cache::ControllerCache, robot_coord_rep=ChartRep(); log_tasks=false) where {M<:Manifold,I}
    set_configuration!(state, xm)
    setdirty!(cache)
    task_acceleration(xm, vm, node, CM, robot_coord_rep; log_tasks)
end


text:
module ChiralWaveguide

using Reexport

@reexport using QuantumOptics

include("WavePackets/WavePackets.jl")
export WavePacket, ContinuousWave, NoPacket

include("WaveguideProblem/WaveguideProblem.jl")
export WaveguideProblem

import SciMLBase: solve
include("WaveguideProblem/solve.jl")
export solve

# ----------------------------------------
# Predefined Systems
include("Systems/TwoLevelChain.jl")
export TwoLevelChain

include("Systems/DissipativeLambdaChain.jl")
export DissipativeLambdaChain

end


text:
using StateSpaceRoutines, Test, HDF5, JLD2, FileIO, LinearAlgebra, PDMats, Distributions, Random, DSGE

my_tests = [
            "kalman_filter"
            "chand_recursion"
            "tempered_particle_filter"
            "smoothers"
            "dyn_measure_tpf"
            "poolmodel_tpf"
            ]

for test in my_tests
    @testset "$test.jl" begin
        include(string("$test.jl"))
    end
end


text:
module ACC

using NeuralNetworkAnalysis, LaTeXStrings
using NeuralNetworkAnalysis: FunctionPreprocessing
import Plots, DifferentialEquations
using Plots: plot, plot!, lens!, bbox, savefig, font, Measures.mm, annotations

# problem

u = 0.0001
a_lead = -2.0

@taylorize function ACC!(dx, x, p, t)
    v_lead = x[2]
    γ_lead = x[3]
    v_ego = x[5]
    γ_ego = x[6]
    a_ego = x[7]

    dx[1] = v_lead
    dx[2] = γ_lead
    dx[3] = 2 * (a_lead - γ_lead) - u * v_lead^2
    dx[4] = v_ego
    dx[5] = γ_ego
    dx[6] = 2 * (a_ego - γ_ego) - u * v_ego^2
    dx[7] = zero(a_ego)
    return dx
end

controller = read_nnet_sherlock(@modelpath "" "controller_ACC_sherlock")

X₀ = Hyperrectangle(low= [ 90, 32,   0, 10, 30,   0],
                    high=[110, 32.2, 0, 11, 30.2, 0])
U₀ = ZeroSet(1)
vars_idx = Dict(:state_vars=>1:6, :control_vars=>7)
ivp = @ivp(x' = ACC!(x), dim: 7, x(0) ∈ X₀ × U₀)

v_set = 30.0
T_gap = 1.4
M = zeros(3, 6)
M[1, 5] = 1.0
M[2, 1] = 1.0
M[2, 4] = -1.0
M[3, 2] = 1.0
M[3, 5] = -1.0
function preprocess(X::LazySet)
    Y1 = Singleton([v_set, T_gap])
    Y2 = linear_map(M, X)
    return cartesian_product(Y1, Y2)
end
function preprocess(X::AbstractVector)
    Y1 = [v_set, T_gap]
    Y2 = M * X
    return vcat(Y1, Y2)
end
control_preprocessing = FunctionPreprocessing(preprocess)

period = 0.1
k = 50
T = k * period
T_warmup = 2 * period

prob = ControlledPlant(ivp, controller, vars_idx, period;
                       preprocessing=control_preprocessing)

# simulation

println("simulation")
res = @timed simulate(prob, T=T, trajectories=10, include_vertices=true)
sim = res.value
print_timed(res)

# reachability analysis

alg = TMJets(abstol=1e-6, orderT=6, orderQ=1)
alg_nn = Ai2()

function benchmark(; T=T, silent::Bool=false)
    silent || println("flowpipe construction")
    res = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg)
    sol = res.value
    silent || print_timed(res)

    return sol
end

benchmark(T=T_warmup, silent=true)
res = @timed benchmark()
sol = res.value

# plotting

vars = (1, 4)
fig = plot(xlab=L"x_{lead}", ylab=L"x_{ego}",
           xlims=(50, 300), ylims=(0, 200),
           leg=:topleft,
           legendfontsize=25,
           tickfont=font(25, "Times"),
           guidefontsize=25,
           xguidefont=font(25, "Times"),
           yguidefont=font(25, "Times"),
           bottom_margin=0mm, left_margin=0mm, right_margin=5mm, top_margin=2mm, size=(500, 450))
plot!(fig, sol, vars=vars, color=:orange, lab="")
plot!(project(X₀, vars), c=:blue, lab=L"X_0", alpha=0.6)
plot_simulation!(fig, sim; vars=vars, color=:black, lab="")
lens!([89, 111], [9.7, 11.5], inset=(1, bbox(0.65, 0.6, 0.3, 0.25)), xtick=90:10:110, ytick=10:1:11, tickfontsize=15, subplot=2)
savefig("acc_juliareach.png")

end
nothing


text:
using BinDeps

@BinDeps.setup

superlu = library_dependency("libsuperlu")

provides(Sources,URI("http://crd-legacy.lbl.gov/~xiaoye/SuperLU/superlu_4.3.tar.gz"),superlu)

@linux_only begin
    provides(AptGet, "libsuperlu4", superlu)  # created for Debian, not yet on Ubuntu
    provides(Yum, "SuperLU", superlu)
end

@osx_only begin
    if Pkg.installed("Homebrew") === nothing
            error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    provides( Homebrew.HB, "superlu", superlu, os = :Darwin ) # doesn't exist at present
end

@BinDeps.install


text:
function convolve_beam(NN::Int, pix_size::Float64, beam_waist::Float64,
                       Map::Array{Float64, 2})

    gaussian = make_2d_gaus(NN, pix_size, beam_waist)

    FT_beam = fft(fftshift(gaussian))
    FT_Map  = fft(fftshift(Map))
    conv = FT_beam .* FT_Map

    conv_map = fftshift(real.(ifft(conv)))

    return conv_map

end

function make_2d_gaus(N::Int, pix_size::Float64, beam_waist::Float64)

    inds = ((0:(N-1)) .+ 0.5 .- (N/2)) * pix_size
    X = reshape(repeat(inds, N), N ,N)
    Y = X'
    R = (X .^2 .+ Y .^2) .^0.5

    sigma = beam_waist / sqrt(8. * log(2))
    gauss = exp.(-0.5 .* ((R ./ sigma).^2.0))
    Norm_Fact = sum(gauss)
    gauss_norm = gauss ./ Norm_Fact

    return gauss_norm

end


text:
function MPIFiles.saveasMDF(store::DatasetStore, scanner::MPIScanner, sequence::Sequence, data::Array{Float32,4},
                            params::Dict; bgdata::Union{Array{Float32,4},Nothing}=nothing)

  mdf = MDFv2InMemory()
  # / subgroup
  mdf.root = defaultMDFv2Root()

  name = params["studyName"]
  date = params["studyDate"]
  study = Study(store, name; date=date)

  fillMDFStudy(mdf, study)
  fillMDFExperiment(mdf, study, params)
  fillMDFScanner(mdf, scanner, params)
  fillMDFTracer(mdf, params)

  fillMDFMeasurement(mdf, sequence, data, bgdata)
  fillMDFAcquisition(mdf, scanner, sequence)

  filename = getNewExperimentPath(study)

  return saveasMDF(filename, mdf)
end

function MPIFiles.saveasMDF(store::DatasetStore, scanner::MPIScanner, sequence::Sequence, data::Array{Float32,4},
  positions::Positions, isBackgroundFrame::Vector{Bool}, params::Dict)

  mdf = MDFv2InMemory()
  # / subgroup
  mdf.root = defaultMDFv2Root()

  if params["storeAsSystemMatrix"]
    study = MPIFiles.getCalibStudy(store)
  else
    name = params["studyName"]
    date = params["studyDate"]
    study = Study(store, name; date=date)
  end

  fillMDFStudy(mdf, study)
  fillMDFExperiment(mdf, study, params)
  fillMDFScanner(mdf, scanner, params)
  fillMDFTracer(mdf, params)

  fillMDFMeasurement(mdf, sequence, data, isBackgroundFrame)
  fillMDFAcquisition(mdf, scanner, sequence)
  fillMDFCalibration(mdf, positions, params)

  filename = getNewExperimentPath(study)

  return saveasMDF(filename, mdf)
end

function fillMDFCalibration(mdf::MDFv2InMemory, positions::Positions, params::Dict)

  # /calibration/ subgroup

  deltaSampleSize = haskey(params, "calibDeltaSampleSize") ?
       Float64.(ustrip.(uconvert.(Unitful.m, params["calibDeltaSampleSize"]))) : nothing

  subgrid = isa(positions,BreakpointGridPositions) ? positions.grid : positions

  # TODO: THIS NEEDS TO BE DEFINED IN THE MDF! we otherwise cannot store these grids
  # params["calibIsMeanderingGrid"] = isa(subgrid,MeanderingGridPositions)

  fov = Float64.(ustrip.(uconvert.(Unitful.m, fieldOfView(subgrid))))
  fovCenter = Float64.(ustrip.(uconvert.(Unitful.m, fieldOfViewCenter(subgrid))))
  size = shape(subgrid)

  method = "robot"
  offsetFields = nothing
  order = "xyz"
  positions = nothing
  snr = nothing

  mdf.calibration = MDFv2Calibration(;
      deltaSampleSize = deltaSampleSize,
      fieldOfView = fov,
      fieldOfViewCenter = fovCenter,
      method = method,
      offsetFields = offsetFields,
      order = order,
      positions = positions,
      size = size,
      snr = snr)

  return
end

function fillMDFScanner(mdf::MDFv2InMemory, scanner::MPIScanner, params::Dict)
  # /scanner/ subgroup
  mdf.scanner = MDFv2Scanner(
    boreSize = ustrip(u"m", scannerBoreSize(scanner)),
    facility = scannerFacility(scanner),
    manufacturer = scannerManufacturer(scanner),
    name = scannerName(scanner),
    operator = get(params, "scannerOperator", "default"),
    topology = scannerTopology(scanner)
  )
  return
end

function fillMDFStudy(mdf::MDFv2InMemory, study::Study)
  # /study/ subgroup
  mdf.study = MDFv2Study(;
    description = "n.a.",
    name = study.name,
    number = 0, # FIXME: This is never set!!!!!!
    time = study.date,
    uuid = study.uuid
  )
  return
end

function fillMDFExperiment(mdf::MDFv2InMemory, study::Study, params::Dict)
  # /experiment/ subgroup

  expNum = getNewExperimentNum(study)

  mdf.experiment = MDFv2Experiment(;
      description = get(params,"experimentDescription","n.a."),
      isSimulation = false,
      name = get(params,"experimentName","default"),
      number = expNum,
      subject = get(params,"experimentSubject","n.a."),
      uuid = uuid4()
    )
  return
end

function fillMDFTracer(mdf::MDFv2InMemory, params::Dict)
  # /tracer/ subgroup
  mdf.tracer = MDFv2Tracer(;
      batch = get(params,"tracerBatch",["n.a"]),
      concentration = get(params,"tracerConcentration",[0.0]),
      injectionTime = [t for t in get(params,"tracerInjectionTime", [Dates.unix2datetime(time())]) ],
      name = get(params,"tracerName",["n.a"]),
      solute = get(params,"tracerSolute",["Fe"]),
      vendor = get(params,"tracerVendor",["n.a"]),
      volume = get(params,"tracerVolume",[0.0]))
  return
end

function fillMDFMeasurement(mdf::MDFv2InMemory, sequence::Sequence, data::Array{Float32,4},
                            bgdata::Union{Array{Float32,4}, Nothing})

  # /measurement/ subgroup
  numFrames = acqNumFrames(sequence)

  if bgdata == nothing
    isBackgroundFrame = zeros(Bool, numFrames)
    data_ = data
  else
    numBGFrames = size(bgdata,4)
    data_ = cat(bgdata, data, dims=4)
    isBackgroundFrame = cat(ones(Bool,numBGFrames), zeros(Bool,numFrames), dims=1)
    numFrames = numFrames + numBGFrames
  end

  return fillMDFMeasurement(mdf, sequence, data_, isBackgroundFrame)
end

function fillMDFMeasurement(mdf::MDFv2InMemory, sequence::Sequence, data::Array{Float32,4},
  isBackgroundFrame::Vector{Bool})

  # /measurement/ subgroup
  numFrames = size(data,4)

  measData(mdf, data)
  measIsBackgroundCorrected(mdf, false)
  measIsBackgroundFrame(mdf, isBackgroundFrame)
  measIsFastFrameAxis(mdf, false)
  measIsFourierTransformed(mdf, false)
  measIsFramePermutation(mdf, false)
  measIsFrequencySelection(mdf, false)
  measIsSparsityTransformed(mdf, false)
  measIsSpectralLeakageCorrected(mdf, false)
  measIsTransferFunctionCorrected(mdf, false)
  return
end

function fillMDFAcquisition(mdf::MDFv2InMemory, scanner::MPIScanner, sequence::Sequence)
  # Needs to be filled after(!) measurement group
  numPeriodsPerFrame_ = acqNumPeriodsPerFrame(sequence)
  numRxChannels_ = rxNumChannels(sequence)
  numSamplingPoints_ = rxNumSamplingPoints(sequence)

  # /acquisition/ subgroup
  acqGradient(mdf, acqGradient(sequence))
  acqNumAverages(mdf, acqNumAverages(sequence))
  acqNumFrames(mdf, length(measIsBackgroundFrame(mdf))) # important since we might have added BG frames
  acqNumPeriodsPerFrame(mdf, acqNumPeriodsPerFrame(sequence))
  acqOffsetField(mdf, acqOffsetField(sequence))
  acqStartTime(mdf, Dates.unix2datetime(time())) #seqCont.startTime)

  # /acquisition/drivefield/ subgroup
  dfBaseFrequency(mdf, ustrip(u"Hz", dfBaseFrequency(sequence)))
  dfCycle(mdf, ustrip(u"s", dfCycle(sequence)))
  dfDivider(mdf, dfDivider(sequence))
  dfNumChannels(mdf, dfNumChannels(sequence))
  dfPhase(mdf, ustrip.(u"rad", dfPhase(sequence)))
  dfStrength(mdf, ustrip.(u"T", dfStrength(sequence)))
  dfWaveform(mdf, fromWaveform.(dfWaveform(sequence)))

  # /acquisition/receiver/ subgroup
  rxBandwidth(mdf, ustrip(u"Hz", rxBandwidth(sequence)))
  convFactor = zeros(2,numRxChannels_)
  convFactor[1,:] .= 1.0
  rxDataConversionFactor(mdf, convFactor)
  rxNumChannels(mdf, numRxChannels_)
  rxNumSamplingPoints(mdf, numSamplingPoints_)
  rxUnit(mdf, "V")

  # transferFunction
  if hasTransferFunction(scanner)
    numFreq = div(numSamplingPoints_,2)+1
    freq = collect(0:(numFreq-1))./(numFreq-1).*ustrip(u"Hz", rxBandwidth(sequence))
    tf_ =  TransferFunction(scanner)
    tf = tf_[freq,1:numRxChannels_]
    rxTransferFunction(mdf, tf)
    rxInductionFactor(mdf, tf_.inductionFactor)
  end

end