text:
# Linear operator implementation for Petsc matrix-explicit

"""
  Petsc matrix-explicit linear operator.

  Note that if an interface common to both Petsc and regular matrices is used
  when accessing the `Ap` field, it is possible to write functions that
  operate on both regular and Petsc linear operators.

  **Public Fields**

   * A: a PetscMat

"""
mutable struct PetscMatLO <: AbstractPetscMatLO
  A::PetscMat
  xtmp::PetscVec  # these are shared with the PC if possible
  btmp::PetscVec
  is_setup::Array{Bool, 1}  # is matrix assembled (shared with PC.is_assembled)
                            # if A is shared
  is_shared::Bool
  is_finalized::Bool
  nassemblies::Array{Int, 1}
  nsolves::Int
  ntsolves::Int

  # MPI stuff
  comm::MPI.Comm
  myrank::Int
  commsize::Int
end

"""
  Outer constructor for PetscMatLO

  **Inputs**

   * pc: a Petsc preconditioner of some kind
   * mesh
   * sbp
   * eqn
   * opts
"""
function PetscMatLO(pc::AbstractPetscPC, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  pc2 = getBasePC(pc)

  # share matrix if possible
  if !opts["use_jac_precond"] && !(typeof(pc2) <: PetscMatFreePC)
    A = pc2.A
    is_setup = pc2.is_assembled  # share the indicator array
    is_shared = true
    nassemblies = pc2.nassemblies
  else
    A = createPetscMat(mesh, sbp, eqn, opts)
    is_setup = Bool[false]
    is_shared = false
    nassemblies = Int[0]
  end

  # update the PC so it knows if the LO used its matrix or not
  getBasePC(pc).is_shared = is_shared

  if typeof(pc2) <: PetscMatPC  # if pc2 has Petsc vectors inside it
    xtmp = pc2.xtmp
    btmp = pc2.btmp
  else
    xtmp = createPetscVec(mesh, sbp, eqn, opts)
    btmp = createPetscVec(mesh, sbp, eqn, opts)
  end

  is_finalized = false
  nsolves = 0
  ntsolves = 0
  comm = eqn.comm
  myrank = eqn.myrank
  commsize = eqn.commsize

  return PetscMatLO(A, xtmp, btmp, is_setup, is_shared, is_finalized,
                    nassemblies, nsolves, ntsolves, comm, myrank, commsize)
end

function free(lo::PetscMatLO)

  if !lo.is_finalized
    if lo.A.pobj != C_NULL
      PetscDestroy(lo.A)
      lo.A.pobj = C_NULL
    end

    if lo.xtmp.pobj != C_NULL
      PetscDestroy(lo.xtmp)
      lo.xtmp.pobj = C_NULL
    end

    if lo.btmp.pobj != C_NULL
      PetscDestroy(lo.btmp)
      lo.btmp.pobj = C_NULL
    end
  end

  lo.is_finalized = true

  return nothing
end

function calcLinearOperator(lo::PetscMatLO, mesh::AbstractMesh,
                            sbp::AbstractOperator, eqn::AbstractSolutionData,
                            opts::Dict, ctx_residual, t)

#  physicsJac(mesh, sbp, eqn, opts, pc.A, ctx_residual, t)
  setIsSetup(lo, false)

  return nothing
end

function applyLinearOperator(lo::AbstractPetscMatLO, mesh::AbstractMesh,
                             sbp::AbstractOperator, eqn::AbstractSolutionData,
                             opts::Dict, ctx_residual, t, x::AbstractVector,
                             b::AbstractVector)

  lo2 = getBaseLO(lo)

  # assemble Ap (if needed)
  assemblePetscData(lo2, x, lo2.xtmp)

  MatMult(lo2.A, lo2.xtmp, lo2.btmp)

  btmp = VecGetArrayRead(lo2.btmp)
  copy!(b, btmp)
  VecRestoreArrayRead(lo2.btmp, btmp)

  return nothing
end

function applyLinearOperatorTranspose(lo::AbstractPetscMatLO,
                             mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::AbstractSolutionData, opts::Dict,
                             ctx_residual, t, x::AbstractVector,
                             b::AbstractVector)
  lo2 = getBaseLO(lo)

  # assemble Ap (if needed)
  assemblePetscData(lo2, x, lo.xtmp)

  MatMultTranspose(lo2.A, lo2.xtmp, lo2.btmp)

  btmp = VecGetArrayRead(lo2.btmp)
  copy!(b, btmp)
  VecRestoreArrayRead(lo2.btmp, btmp)

  return nothing
end

function assemblePetscData(lo::PetscMatLO, b::AbstractVector,
                           b_petsc::PetscVec)

  myrank = lo.myrank

  # assemble things
  if !getIsSetup(lo)
    MatAssemblyBegin(lo.A, MAT_FINAL_ASSEMBLY)
  end

  # copy values into the vector
  btmp = VecGetArray(b_petsc)
  copy!(btmp, b)
  VecRestoreArray(b_petsc, btmp)

  if !getIsSetup(lo)
    MatAssemblyEnd(lo.A, MAT_FINAL_ASSEMBLY)
    setIsSetup(lo, true)
    lo.nassemblies[1] += 1
    matinfo = MatGetInfo(lo.A, PETSc2.MAT_LOCAL)
    if matinfo.mallocs > 0.5  # if any mallocs
      println(BSTDERR, "Warning: non-zero number of mallocs for A on process $myrank: $(matinfo.mallocs) mallocs")
    end
  end

  return nothing
end

# this function works for all linear operators
function getIsSetup(lo::AbstractLO)

  lo2 = getBaseLO(lo)
  return lo2.is_setup[1]
end

function setIsSetup(lo::AbstractLO, val::Bool)

  lo2 = getBaseLO(lo)
  lo2.is_setup[1] = val

  return nothing
end

function getBaseLO(lo::PetscMatLO)

  # this is the bottom of the recursion tree
  return lo
end


text:
@cache struct PDIRK44Cache{uType,rateType,N,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::Array{rateType}
  k2::Array{rateType}
  nlsolver::N
  tab::TabType
end

struct PDIRK44ConstantCache{N,TabType} <: OrdinaryDiffEqConstantCache
  nlsolver::N
  tab::TabType
end

struct PDIRK44Tableau{T,T2}
  γs::SVector{2,T2}
  cs::SVector{4,T2}
  α1::SVector{2,T}
  α2::SVector{2,T}
  b1::T
  b2::T
  b3::T
  b4::T
end

function PDIRK44Tableau(T, T2)
  γ1 = convert(T2,  1//2)
  γ2 = convert(T2,  2//3)
  γs = SVector(γ1, γ2)
  c1 = convert(T2,  1//2)
  c2 = convert(T2,  2//3)
  c3 = convert(T2,  1//2)
  c4 = convert(T2,  1//3)
  cs = SVector(c1,c2,c3,c4)
  α11 = convert(T, -5//2)
  α12 = convert(T, -5//3)
  α1 = SVector(α11, α12)
  α21 = convert(T,  5//2)
  α22 = convert(T,  4//3)
  α2 = SVector(α21, α22)
  b1 = convert(T, -1//1)
  b2 = convert(T, -1//1)
  b3 = convert(T,  3//2)
  b4 = convert(T,  3//2)
  PDIRK44Tableau{T,T2}(γs,cs,α1,α2,b1,b2,b3,b4)
end

function alg_cache(alg::PDIRK44,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  γ, c = 1.0, 1.0
  if alg.threading
    nlsolver1 = build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(true))
    nlsolver2 = build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(true))
    nlsolver = [nlsolver1, nlsolver2]
  else
    _nlsolver = build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(true))
    nlsolver = [_nlsolver]
  end
  tab = PDIRK44Tableau(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  k1 = [zero(rate_prototype) for i in 1:2 ]
  k2 = [zero(rate_prototype) for i in 1:2 ]
  PDIRK44Cache(u,uprev,k1,k2,nlsolver,tab)
end

function alg_cache(alg::PDIRK44,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  γ, c = 1.0, 1.0
  if alg.threading
    nlsolver1 = build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(false))
    nlsolver2 = build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(false))
    nlsolver = [nlsolver1, nlsolver2]
  else
    _nlsolver = build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(false))
    nlsolver = [_nlsolver]
  end
  tab = PDIRK44Tableau(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  PDIRK44ConstantCache(nlsolver,tab)
end

ext: jl


text:
using Ipopt
using ForwardDiff

model = TrajectoryOptimization.Dynamics.pendulum[1]
model! = TrajectoryOptimization.Dynamics.pendulum![1]
obj = TrajectoryOptimization.Dynamics.pendulum[2]
obj = TrajectoryOptimization.ConstrainedObjective(obj, u_min=-2, u_max=2)
R = obj.R; Q = obj.Q; Qf = obj.Qf; xf = obj.xf
f = model.f
n, m = 2,1
N = 50
dt = 0.1

f_aug! = f_augmented!(model!.f, model.n, model.m)
zdot = zeros(2)
F(z) = ForwardDiff.jacobian(f_aug!,zdot,z)
F([1,2,1])
vec(F([1,2,1]))

weights = ones(N)
weights[[1,end]] = 0.5

function packZ(X,U)
    n, N = size(X)
    m = size(U,1)
    Z = zeros(n+m,N)
    Z[1:n,:] .= X
    Z[n+1:end,1:end] .= U
    Z = reshape(Z,1,(n+m)N)
end

function unpackZ(Z)
    Z = reshape(Z,n+m,N)
    X = Z[1:n,:]
    U = Z[n+1:end,:]
    return X,U
end

function cost(X,U)
    J = 0
    for k = 1:N-1
      J += dt/2*(cost_k(X[:,k],U[:,k]) + cost_k(X[:,k+1],U[:,k+1]))
    end
    J += 0.5*(X[:,N] - xf)'*Qf*(X[:,N] - xf)
end

function cost_k(x,u)
    0.5*(x-xf)'Q*(x-xf) + 0.5*u'R*u
end

function eval_f(Z)
    X,U = unpackZ(Z)
    cost(X,U)
end

function eval_g(Z, g)
    X,U = unpackZ(Z)
    g = reshape(g,n,N-1)
    @show size(g)
    for k = 1:N-1
        g[:,k] = ( f(X[:,k+1],U[:,k+1]) + f(X[:,k],U[:,k]) )/2 - X[:,k+1] + X[:,k]
    end
    # reshape(g,n*(N-1),1)
    return nothing
end

function eval_grad_f(Z, grad_f)
    X, Z = unpackZ(Z)
    grad_f = reshape(grad_f, n+m, N)
    for k = 1:N
        grad_f[1:n,k] = weights[k]*Q*X[:,k]
        grad_f[n+1:end,k] = weights[k]*R*U[:,k]
    end
    grad_f[1:n,N] += Qf*X[:,N]
    return nothing
end

function eval_jac_g(x, mode, rows, cols, vals)
    z = reshape(x,n+m,N)
    if mode == :Structure
        ind = 1
        for k = 1:N-1 # loop over blocks
            for j = 1:2(n+m) # loop over columns of each block
                for i = 1:n # loop over rows of each block
                    row[ind] = i + (k-1)*n
                    col[ind] = j + (k-1)*(n+m)
                    ind += 1
                end
            end
        end
    else
        n_blk = n*2(n+m)
        n_blk2 = n*(n+m)
        fz = F(z[:,1])
        vals[1:n_blk2] .= vec(fz+eye(n,n+m))
        for k = 2:N-1
            fz = F(z[:,k])
            vals[(1:n_blk2) + (k-1)*n_blk - n_blk2] .= vec(fz+eye(n,n+m))
            vals[(1:n_blk2) + (k-1)*n_blk] .= vec(fz-eye(n,n+m))
        end
        fz = F(z[:,N])
        vals[end-n_blk2+1:end] .= vec(fz-eye(n,n+m))
    end
end

function eval_h(x, mode, rows, cols, obj_factor, lambda, values)
    if mode == :Structure
        ind = 1
        for k = 1:N
            n_blk = (n^2 + m^2)
            for j = 1:n
                for i = 1:n
                    row[ind] = i + (k-1)*n_blk
                    col[ind] = j + (k-1)*n_blk
                    ind += 1
                end
            end
            for j = 1:m
                for i = 1:m
                    row[ind] = i + (k-1)*n_blk + n^2
                    col[ind] = j + (k-1)*n_blk + n^2
                    ind += 1
                end
            end
        end
    end
end

X = rand(n,N)
U = rand(m,N)
Z = packZ(X,U)

X2,U2 = unpackZ(Z)
g = zeros(n*(N-1))
eval_g(Z,g)
g
grad_f = zeros((n+m)*N)
eval_grad_f(Z,grad_f)

nnz_jac_g = (N-1)*n*2(n+m)
row = zeros(Int,nnz_jac_g)
col = zeros(Int,nnz_jac_g)
val = zeros(nnz_jac_g)
eval_jac_g(Z,:Structure, row, col, val)

jac_g0 = zeros(Int,(N-1)*n,N*(n+m))
inds = sub2ind(jac_g0,row,col)
jac_g0[inds] .= 1:nnz_jac_g
jac_g0

eval_jac_g(Z,:None,row,col,val)
val

nnz_h = (n^2 + m^2)*N
h = zeros(N*(n+m), N*(n+m))
row = zeros(Int,nnz_h)
col = zeros(Int,nnz_h)
vals = zeros(nnz_h)
lambda = zeros(n, N-1)
eval_h(Z, :Structure, row, col, 1, lambda, vals)

row

col
inds = CartesianIndex.(row[1:10],col[1:10])
h[inds] .= 1:10
h


text:
using GD
using GD.Storage: Compressor, Store, GDFile, compress!, extract, patch, unpatch,
    setup_api_endpoint, validate_remote!

"""
    GDPayload(data, minval, maxval)

Payload generated by the GD compression. Since the GD compression scheme is
built on top of quantization, the same metadata are presnts, while the structure
of the data is under the form of a small `GDFile`.
"""
mutable struct GDPayload
    data::GDFile
    minval::Float32
    maxval::Float32
end

"""
    GDPayloadSerde(chunksize, fingerprint, msbsize, store_host, store_port, is_client)

Apply both quantization and GD to the weights before serialization in order to
compress the payload. Quantization is a lossy process, thus a loss of
information is to be expected after deserialization. Depending on the chosen
fingerprint, collisions can occur, also leading to a loss of information.
"""
mutable struct GDPayloadSerde{T <: Unsigned} <: PayloadSerde
    store::Store
    is_client::Bool
    gdfile::GDFile

    GDPayloadSerde{T}(
        chunksize::Int,
        fingerprint::Function,
        msbsize::T,
        store_host::String,
        store_port::Int,
        is_client::Bool
    ) where T <: Unsigned = begin
        quantizer = GD.Transform.Quantizer{T}(chunksize, msbsize)
        compressor = Compressor(chunksize, quantizer, fingerprint)
        store = Store(compressor, Dict(), 0, 0)
        @async setup_api_endpoint(store, store_host, store_port)

        return new(store, is_client, GDFile(Vector(undef, 0), Vector(undef, 0), 0))
    end
end

"""
    serialize_payload(::GDPayloadSerde, weights)

Serialize `weights` with the `GDPayloadSerde` where quantization and
generalized deduplication are applied before serialization.
"""
function serialize_payload(
    p::GDPayloadSerde{T},
    weights::Vector{Float32}
)::Vector{UInt8} where T <: Unsigned
    # quantize weights
    q = Quantizer{T}(weights)
    qweights = [quantize(q, w) for w in weights]

    # gd compression
    gdfile = compress!(p.store, qweights)

    # patching
    if p.is_client
        gdfile = GD.patch(gdfile, p.gdfile)
    else
        p.gdfile = gdfile
    end

    STATS.base.req_data = gdfile

    payload = GDPayload(gdfile, q.minval, q.maxval)

    return pack(payload)
end

"""
    deserialize_payload(GDPayloadSerde, data, from)

Deserialize `data` with the `GDPayloadSerde` where generalized deduplication
and dequantization are applied before deserialization.
"""
function deserialize_payload(
    p::GDPayloadSerde{T},
    data::Vector{UInt8},
    from::String
)::Vector{Float32} where T <: Unsigned
    payload = unpack(data)

    # patching
    if p.is_client
        p.gdfile = payload.data
    else
        num_identical_hash = sum([1 for el in payload.data.hashes if el == [0x00]])
        num_identical_dev = sum([1 for el in payload.data.deviations if el == [0x00]])
        STATS.network.num_identical_hashes += num_identical_hash
        STATS.network.num_identical_devs += num_identical_dev

        payload.data = GD.unpatch(payload.data, p.gdfile)
    end

    push!(STATS.base.res_data, payload.data)

    # gd decompression
    STATS.network.num_unknown_bases += validate_remote!(p.store, payload.data, from)
    # since deserializing a client payload means the client is done with his work, we
    # can also update the number of requested bases, to make sure the one requested by
    # this client are taken into account.
    STATS.network.num_requested_bases = p.store.num_requested_bases

    qweights = extract(p.store, payload.data)

    # dequantize weights
    q = Quantizer{T}(payload.minval, payload.maxval)
    weights = [dequantize(q, w) for w in qweights]

    return weights
end


text:
#### Land sources

function land_source!(
    f::Function,
    land::LandModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    f(land, source, state, diffusive, aux, t, direction)
end

function land_source!(
    ::Nothing,
    land::LandModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
) end

# sources are applied additively

@generated function land_source!(
    stuple::Tuple,
    land::LandModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    N = fieldcount(stuple)
    return quote
        Base.Cartesian.@nexprs $N i -> land_source!(
            stuple[i],
            land,
            source,
            state,
            diffusive,
            aux,
            t,
            direction,
        )
        return nothing
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module Server
export serve

using HTTP, Sockets, JSON
using ..Types, ..Execute
import ..Types.GraphQLServer

function parsePayload(payload)
  return ExecutionArgs("TODO")
end

function GraphQLServer(config=ServerConfiguration())
  server = Sockets.listen(Sockets.InetAddr(Sockets.localhost, config.port))
  @async HTTP.serve(server=server) do request::HTTP.Request
    try
      payload = String(HTTP.payload(request))
      args = parsePayload(payload)
      result = execute(args)
      return HTTP.Response(json(result))
    catch e
      return HTTP.Response(500, "Something went wrong. Please open an issue to resolve.")
    end
 end
 return GraphQLServer(server,"$(Sockets.localhost):$(config.port)",config)
end

end


text:
"""
activation(layer)

Return activation function of the layer.
In case the layer is unknown or no activation function is found, `nothing` is returned.
"""
activation(l::Dense) = l.σ
activation(l::Conv) = l.σ
activation(l::BatchNorm) = l.λ
activation(layer) = nothing # default for all other layer types

"""
    flatten_model(c)

Flatten a Flux chain containing Flux chains.
"""
function flatten_model(chain::Chain)
    if any(isa.(chain.layers, Chain))
        flatchain = Chain(vcat(_flatten_model.(chain.layers)...)...)
        return flatten_model(flatchain)
    end
    return chain
end
@deprecate flatten_chain(c) flatten_model(c)

_flatten_model(x) = x
_flatten_model(c::Chain) = [c.layers...]

is_softmax(x) = x isa SoftmaxActivation
has_output_softmax(x) = is_softmax(x) || is_softmax(activation(x))
has_output_softmax(model::Chain) = has_output_softmax(model[end])

"""
    check_output_softmax(model)

Check whether model has softmax activation on output.
Return the model if it doesn't, throw error otherwise.
"""
function check_output_softmax(model::Chain)
    if has_output_softmax(model)
        throw(ArgumentError("""Model contains softmax activation on output.
                            Call `strip_softmax` on your model first."""))
    end
    return model
end

"""
    strip_softmax(model)

Remove softmax activation on model output if it exists.
"""
function strip_softmax(model::Chain)
    if has_output_softmax(model)
        model = flatten_model(model)
        if is_softmax(model[end])
            return Chain(model.layers[1:(end - 1)]...)
        end
        return Chain(model.layers[1:(end - 1)]..., strip_softmax(model[end]))
    end
    return model
end
strip_softmax(l::Union{Dense,Conv}) = set_params(l, l.weight, l.bias, identity)

# helper function to work around `bias=false` (Flux v0.13) and `bias=Flux.Zeros` (v0.12)
function get_params(layer)
    W = layer.weight
    b = layer.bias
    if b == false || typeof(b) <: Flux.Zeros
        b = zeros(eltype(W), size(W, 1))
    end
    return W, b
end

"""
    set_params(layer, W, b)

Duplicate layer using weights W, b.
"""
set_params(l::Conv, W, b, σ=l.σ) = Conv(σ, W, b, l.stride, l.pad, l.dilation, l.groups)
set_params(l::Dense, W, b, σ=l.σ) = Dense(W, b, σ)


text:
import Dance

include("./utils/main.jl")

Dance.start_project("demo")
project_settings()
routes("json")

## Test all routes with pending slash, to ensure is removed ##
@testset "HTTP.listen" begin
    @async include(joinpath(abspath(pwd()), "dance.jl"))
    sleep(1)

    make_and_test_request_get("/dict/", 200, Dict("Content-Type" => "application/json"), 9, true, Dict("a" => 123))

    make_and_test_request_get(
        "/dataframe/",
        200,
        Dict("Content-Type" => "application/json"),
        43,
        true,
        [
            ["A", "B"],
            [1, "M"],
            [2, "F"],
            [3, "F"],
            [4, "M"]
        ]
    )

    # Test int url param
    # Test using OPTIONS method before POST (see https://github.com/axios/axios/issues/475)
    make_and_test_request_options("/post/dict/12")
    make_and_test_request_post("/post/dict/12", Dict("b" => "abc"), 200, Dict("Content-Type" => "application/json"), 8, true, Dict("b" => 12))

    # Test setting Header value in backend
    make_and_test_request_post(
        "/post/dataframe/",
        [
            ["A", "B"],
            [1, "M"],
            [2, "F"],
            [3, "G"],
            [4, "Z"]
        ],
        200,
        Dict("Content-Type" => "application/json", "foo" => "bar"),
        43,
        true,
        [
            ["A", "B"],
            [1, "M"],
            [2, "F"],
            [3, "G"],
            [4, "Z"]
        ]
    )
end

delete_project()


text:
abstract type DTree{K} end
struct Leaf{K} <: DTree{K}
	term::Any
end
struct Fail{K} <: DTree{K} end

struct Switch{K, T} <: DTree{K}
	occ::Vector{Int}
	cases::Vector{Pair{K, DTree{K}}}
	default::Union{DTree{K}, Nothing}
end

const MatchMatrix{K} = Vector{Pair{V, T}} where {P<:Pat, V<:Vector{P}, T}

preprocess_variables(occ::Vector{Int}, p::StarPat, occMap::Dict{Vector{Int}, Symbol}) = (p, Dict{Symbol, Symbol}())
preprocess_variables(occ::Vector{Int}, v::VarPat, occMap::Dict{Vector{Int}, Symbol}) = (StarPat(), Dict{Symbol, Symbol}(v.var=>get!(()->gensym(), occMap, occ)))
function preprocess_variables(occ::Vector{Int}, v::CstrPat, occMap::Dict{Vector{Int}, Symbol})
	out = Dict{Symbol, Symbol}()
	out_pats = Pat[]
	for (i,arg) in enumerate(v.args)
		pat, vars = preprocess_variables([occ; i], arg, occMap)
		merge!(out, vars)
		push!(out_pats, pat)
	end
	return (CstrPat(v.cstr, out_pats), out)
end

function preprocess_variables(P::Vector{Vector{Pat}}, A::Vector{T}) where {T}
	occMap = Dict{Vector{Int}, Symbol}()
	out_pats = Vector{Pat}[]
	out_exprs = Expr[]
	for (pattern, expr) in zip(P, A)
		varMap = Dict{Symbol, Symbol}()
		pats = Pat[]
		for (i, pat_elem) in enumerate(pattern)
			pat,vars = preprocess_variables([i], pat_elem, occMap)
			merge!(varMap, vars)
			push!(pats, pat)
		end
		final_expr = Expr(:block, (Expr.(:(=), keys(varMap), values(varMap)))..., expr)
		push!(out_exprs, final_expr)
		push!(out_pats, pats)
	end
	return out_pats, out_exprs, occMap
end
function preprocess_variables(P::Vector{Pair{Vector{Pat}, T}}) where {T}
	pats, exprs, occMap = preprocess_variables(first.(P), last.(P))
	return Pair.(pats, exprs), occMap
end

function S(constructors, t, cstr_head, i, P::Vector{Vector{Pat}}, A::Vector{T}) where {T}
	rowty = Pair{Vector{Pat}, T}
	patlen = length(first(P))
	@assert all(l->l==patlen, length.(P))
	@assert length(P) == length(A)
	arguments = constructors(t)[cstr_head]

	adj_prefix = getindex.(P, (1:(i-1),))
	patterns = getindex.(P, i)
	adj_postfix = getindex.(P, (i+1:patlen, ))
	makerows(prefix, c::CstrPat, postfix, a) =
		if c.cstr == cstr_head
			@assert length(c.args) == length(arguments)
			[[prefix; c.args; postfix] => a]
		else
			rowty[]
		end
	makerows(prefix, ::StarPat, postfix, a) = rowty[[prefix; repeat([StarPat()], length(arguments)); postfix] => a]
	makerows(prefix, o::OrPat, postfix, a) = [S(constructors, t, cstr_head, i, [[prefix; o.left; postfix]], [a]); S(constructors, t, cstr_head, i, [[prefix; o.right; postfix]], [a])]
	out = vcat((makerows.(adj_prefix, patterns, adj_postfix, A))...)
	return out
end
function D(i, P::Vector{Vector{Pat}}, A::Vector{T}) where {T}
	rowty = Pair{Vector{Pat{K}}, T}
	makerows(prefix, ::CstrPat, postfix, a) = rowty[]
	makerows(prefix, ::StarPat, postfix, a) = rowty[[prefix; postfix] => a]
	makerows(prefix, o::OrPat, postfix, a) = [D([[prefix; o.left; postfix]], a); D([[prefix; o.right; postfix]], a)]

	patlen = length(first(P))
	patterns = getindex.(P, i)
	adj_prefix = getindex.(P, (1:(i-1),))
	patterns = getindex.(P, i)
	adj_postfix = getindex.(P, (i+1:patlen, ))
	return vcat((makerows.(adj_prefix, patterns, adj_postfix, A))...)
end

collect_binders(pats::Vector{P} where P<:Pat) = collect_binders.(pats)
collect_binders(::StarPat) = nothing
function cc(o::Vector{Vector{Int}}, P::MatchMatrix{K}, types::Vector{T} where T<:Type) where K
	if length(P) == 0 return Fail{K}() end
	firstpat = first(P)
	if all(iswildcard, first(firstpat)) return Leaf{K}(last(firstpat)) end

	i=first(filter((!)∘isnothing, findfirst.((!)∘iswildcard, first.(P))))
	elem_type = types[i]
	cstr_heads = reduce(union, heads.(getindex.(first.(P), i)); init=Set{K}())
	cstrs = constructors(elem_type)
	possible_heads = keys(cstrs)
	if issubset(possible_heads, cstr_heads) && !issubset(cstr_heads,possible_heads)
		throw(Exception("Invalid head!"))
	end

	Pats = first.(P)
	A = last.(P)
	switch = Pair{K, DTree{K}}[head => cc([o[1:i-1]; vcat.((o[i],), 1:length(cstrs[head])); o[i+1:end]], S(elem_type, head, i, Pats, A), vcat(types[1:i-1], cstrs[head]..., types[i+1:end])) for head in cstr_heads]
	default = nothing
	println(cstr_heads)
	println(possible_heads)
	if !issetequal(cstr_heads, possible_heads)
		default = cc(vcat(o[1:i-1], o[i+1:end]), D(i, Pats, A), vcat(types[1:i-1], types[i+1:end]))
	end
	return Switch{K, elem_type}(o[i], switch, default)
end


text:
export ⋅, absdot, magnitude², magnitude, normalize,
×, coordinate_system_from, permute, face_forward

# General functions for VectorLikes.

@inline ⋅(u::VectorLike, v::VectorLike) = sum(u.coordinates .* v.coordinates)

@inline absdot(u::VectorLike, v::VectorLike) = abs(u ⋅ v)

@inline magnitude²(v::VectorLike) = v ⋅ v

@inline magnitude(v::VectorLike) = √(magnitude²(v))

@inline normalize(v::VectorLike) = v / magnitude(v)

@inline Base.min(a::Vect{N,T}, b::Vect{N,U}) where {N,T,U} = promote_type(typeof(a), typeof(b))(min.(a.coordinates, b.coordinates)...)

@inline Base.max(a::Vect{N,T}, b::Vect{N,U}) where {N,T,U} = promote_type(typeof(a), typeof(b))(max.(a.coordinates, b.coordinates)...)

# Vector3-specific functions.

@inline ×(a::VectorLike{3,T}, b::VectorLike{3,U}) where {T,U} =
        Vector3{promote_type(T,U)}(
            a.y * b.z - a.z * b.y,
            a.z * b.x - a.x * b.z,
            a.x * b.y - a.y * b.x)

@inline face_forward(u::VectorLike{3,T}, v::VectorLike{3,U}) where {T,U} =
        u ⋅ v < zero(promote_type(T,U)) ? -u : u

@inline function coordinate_system_from(v₁::Vector3{T}) where T
    if abs(v₁.x) > abs(v₁.y)
        v₂ = typeof(v₁)(-v₁.z, 0.0, v₁.x) / √(v₁.x * v₁.x + v₁.z * v₁.z)
    else
        v₂ = typeof(v₂)(0.0, v₁.z, -v₁.y) / √(v₁.y * v₁.y + v₁.z * v₁.z)
    end
    v₁, v₂, v₁ × v₂
end

@inline permute(a::Vector3{T}, (x, y, z)) where T = typeof(a)(a[x], a[y], a[z])


text:
# Autogenerated wrapper script for Box2D_jll for aarch64-apple-darwin
export libbox2d

JLLWrappers.@generate_wrapper_header("Box2D")
JLLWrappers.@declare_library_product(libbox2d, "@rpath/libbox2d.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libbox2d,
        "lib/libbox2d.2.4.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using SymbolicUtils, SymbolicUtils.Code, Test
@testset "CSE" begin
    @syms x
    t = cse(hypot(hypot(cos(x), sin(x)), atan(cos(x), sin(x))))

    @test t isa Let
    @test length(t.pairs) == 2
    @test occursin(t.pairs[1].lhs, t.body)
    @test occursin(t.pairs[2].lhs, t.body)
end


text:
export Invite

"""
An invite to a [`Guild`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/invite#invite-object).
"""
struct Invite
    code::String
    guild::Union{Guild, Missing}
    channel::DiscordChannel
    approximate_presence_cound::Union{Int, Missing}
    approximate_member_count::Union{Int, Missing}
end
@boilerplate Invite :dict :docs :lower :merge


text:
import Base: deleteat!, insert!, findfirst
function deleteat!(a::CircularArray{T,1,A}, r::UnitRange{U}) where {T, A, U <: Integer}
    ind = mod1.(r, length(a))
    deleteat!(a.data, sort(ind))
    a
end

function insert!(a::CircularArray{T,1,A}, i::Integer, item) where {T,A}
    ind = mod1(i, length(a))
    insert!(a.data, ind, item)
end
function insert!(a::CircularArray{T,1,A}, i::Integer, items::AbstractArray{T,1}) where {T,A}
    ind = mod1(i, length(a))
    insert!.( (a.data,), ind, reverse(items))
    a
end
function findfirst(predicate::Function, a::CircularArray)
    findfirst(predicate, a.data)
end


text:
# libraries

const libeng = Ref{Ptr{Cvoid}}()
const libmx  = Ref{Ptr{Cvoid}}()
const libmat = Ref{Ptr{Cvoid}}()

# matlab engine functions

const eng_open          = Ref{Ptr{Cvoid}}()
const eng_close         = Ref{Ptr{Cvoid}}()
const eng_set_visible   = Ref{Ptr{Cvoid}}()
const eng_get_visible   = Ref{Ptr{Cvoid}}()
const eng_output_buffer = Ref{Ptr{Cvoid}}()
const eng_eval_string   = Ref{Ptr{Cvoid}}()
const eng_put_variable  = Ref{Ptr{Cvoid}}()
const eng_get_variable  = Ref{Ptr{Cvoid}}()

# mxarray functions

const mx_destroy_array   = Ref{Ptr{Cvoid}}()
const mx_duplicate_array = Ref{Ptr{Cvoid}}()

# functions to access mxarray

const mx_free         = Ref{Ptr{Cvoid}}()

const mx_get_classid  = Ref{Ptr{Cvoid}}()
const mx_get_m        = Ref{Ptr{Cvoid}}()
const mx_get_n        = Ref{Ptr{Cvoid}}()
const mx_get_nelems   = Ref{Ptr{Cvoid}}()
const mx_get_ndims    = Ref{Ptr{Cvoid}}()
const mx_get_elemsize = Ref{Ptr{Cvoid}}()
const mx_get_data     = Ref{Ptr{Cvoid}}()
const mx_get_dims     = Ref{Ptr{Cvoid}}()
const mx_get_nfields  = Ref{Ptr{Cvoid}}()
const mx_get_pr       = Ref{Ptr{Cvoid}}()
const mx_get_pi       = Ref{Ptr{Cvoid}}()
const mx_get_ir       = Ref{Ptr{Cvoid}}()
const mx_get_jc       = Ref{Ptr{Cvoid}}()

const mx_is_double    = Ref{Ptr{Cvoid}}()
const mx_is_single    = Ref{Ptr{Cvoid}}()
const mx_is_int64     = Ref{Ptr{Cvoid}}()
const mx_is_uint64    = Ref{Ptr{Cvoid}}()
const mx_is_int32     = Ref{Ptr{Cvoid}}()
const mx_is_uint32    = Ref{Ptr{Cvoid}}()
const mx_is_int16     = Ref{Ptr{Cvoid}}()
const mx_is_uint16    = Ref{Ptr{Cvoid}}()
const mx_is_int8      = Ref{Ptr{Cvoid}}()
const mx_is_uint8     = Ref{Ptr{Cvoid}}()
const mx_is_char      = Ref{Ptr{Cvoid}}()

const mx_is_numeric   = Ref{Ptr{Cvoid}}()
const mx_is_logical   = Ref{Ptr{Cvoid}}()
const mx_is_complex   = Ref{Ptr{Cvoid}}()
const mx_is_sparse    = Ref{Ptr{Cvoid}}()
const mx_is_empty     = Ref{Ptr{Cvoid}}()
const mx_is_struct    = Ref{Ptr{Cvoid}}()
const mx_is_cell      = Ref{Ptr{Cvoid}}()

# functions to create & delete MATLAB arrays

const mx_create_numeric_matrix = Ref{Ptr{Cvoid}}()
const mx_create_numeric_array  = Ref{Ptr{Cvoid}}()

const mx_create_double_scalar  = Ref{Ptr{Cvoid}}()
const mx_create_logical_scalar = Ref{Ptr{Cvoid}}()

const mx_create_sparse         = Ref{Ptr{Cvoid}}()
const mx_create_sparse_logical = Ref{Ptr{Cvoid}}()

const mx_create_string         = Ref{Ptr{Cvoid}}()
const mx_create_char_array     = Ref{Ptr{Cvoid}}()

const mx_create_cell_array     = Ref{Ptr{Cvoid}}()

const mx_create_struct_matrix  = Ref{Ptr{Cvoid}}()
const mx_create_struct_array   = Ref{Ptr{Cvoid}}()

const mx_get_cell              = Ref{Ptr{Cvoid}}()
const mx_set_cell              = Ref{Ptr{Cvoid}}()

const mx_get_field             = Ref{Ptr{Cvoid}}()
const mx_set_field             = Ref{Ptr{Cvoid}}()
const mx_get_field_bynum       = Ref{Ptr{Cvoid}}()
const mx_get_fieldname         = Ref{Ptr{Cvoid}}()

const mx_get_string = Ref{Ptr{Cvoid}}()

# load I/O mat functions

const mat_open         = Ref{Ptr{Cvoid}}()
const mat_close        = Ref{Ptr{Cvoid}}()
const mat_get_variable = Ref{Ptr{Cvoid}}()
const mat_put_variable = Ref{Ptr{Cvoid}}()
const mat_get_dir      = Ref{Ptr{Cvoid}}()


text:
"""
mpcc_solve: Relaxation method for the MPCC.

`mpcc_solve(:: MPCCStopping; verbose :: Bool = true, kwargs...)`

Note: kwargs are passed to the ParamMPCC structure

See also: *ParamMPCC*, *solveIpopt*
"""
function mpcc_solve(stp :: MPCCStopping; verbose :: Bool = true, kwargs...)

 mod = stp.pb
 parammpcc = ParamMPCC(2*(mod.meta.nvar + mod.meta.ncon + mod.meta.ncc); kwargs...)

 fill_in!(stp, stp.current_state.x)
 stp_rlx = _rlx_init(stp, parammpcc)

 OK = start!(stp)

 while !OK

  stp_rlx = rlx_solve(stp_rlx, stp.current_state.x)

  #stopping
  _reforward!(stp_rlx, stp)
  OK = stop!(stp)

  #additional stop if the parameters are too small
  stp.meta.suboptimal = max(stp_rlx.pb.r, stp_rlx.pb.s, stp_rlx.pb.t) < parammpcc.paramin
  #update the sub-problem
  _rlx_update!(stp_rlx, parammpcc)

  #output
  verbose && printstyled("$(stp.meta.nb_of_stop) Iterate: $(stp.current_state.x), r = $(stp_rlx.pb.r), s = $(stp_rlx.pb.s), t = $(stp_rlx.pb.t) \n", color = :red)
  #verbose && print_with_color(:red, "End - Min: l=$(ma.ractif.iter) |x|=$(norm(xjkl,Inf)) |c(x)|=$(ma.ractif.norm_feas) |L'|=$(ma.sts.optimality)  ρ=$(norm(ma.pen.ρ,Inf)) prec=$(ma.sts.atol) \n")

 end

 return stp
end

"""
Refill the mpccstp.current_state using stp.current_state

typeof(stp_rlx.pb) has to be a RlxMPCC
"""
function _reforward!(stp_rlx :: NLPStopping, stp_mpcc :: MPCCStopping)
 #fill_in!(stp_mpcc, stp_nlp.current_state.x)
 ncon, ncc = stp_mpcc.pb.meta.ncon, stp_mpcc.pb.meta.ncc
 r, s, t, tb = stp_rlx.pb.r, stp_rlx.pb.s, stp_rlx.pb.t, stp_rlx.pb.tb

 state = stp_rlx.current_state
 lG, lH = - state.lambda[ncon+1:ncon+ncc], - state.lambda[ncon+ncc+1:ncon+2*ncc]
 lphi   = state.lambda[ncon+2*ncc+1:ncon+3*ncc]
 dp = dphi(state.cx[ncon+1:ncon+ncc], state.cx[ncon+ncc+1:ncon+2*ncc], r, s, t)
 update!(stp_mpcc.current_state, x   = state.x,
                                 fx  = state.fx,
                                 gx  = state.gx,
                                 Hx  = state.Hx,
                                 mu  = state.mu,
                                 cx  = state.cx[1:ncon],
                                 Jx  = state.Jx[1:ncon,:],
                                 lambda = state.lambda[1:ncon],
                                 cGx = state.cx[ncon+1:ncon+ncc] .- tb,
                                 JGx = state.Jx[ncon+1:ncon+ncc,:],
                                 lambdaG = lG - lphi .* dp[1:ncc],
                                 cHx = state.cx[ncon+ncc+1:ncon+2*ncc] .- tb,
                                 JHx = state.Jx[ncon+ncc+1:ncon+2*ncc,:],
                                 lambdaH = lH - lphi .* dp[ncc+1:2*ncc])
 return stp_mpcc
end

"""
rlx_solve: Solve the regularized problem

`rlx_solve(stp :: NLPStopping, x0 :: AbstractVector)`

"""
function rlx_solve(stp :: NLPStopping, x0 :: AbstractVector)
    reinit!(stp, rstate = true, x = x0)
    solveIpopt(stp)
 return stp
end

"""
Initialize the solve sub-problem structure
"""
function _rlx_init(stp :: MPCCStopping, parammpcc :: ParamMPCC)

 #Initialize the sub-problem
 (r,s,t) = parammpcc.initrst()
 tb      = parammpcc.tb(r, s, t)
 #ρ       = parammpcc.rho_init
 prec    = parammpcc.prec_oracle(r, s, t, parammpcc.precmpcc)

 rlx = RlxMPCC(stp.pb, r, s, t, tb)

 ncon, ncc = stp.pb.meta.ncon, stp.pb.meta.ncc
 state = stp.current_state
 dp = dphi(state.cGx, state.cHx, r, s, t)
 rlx_state = NLPAtX(rlx.meta.x0, vcat(state.lambda, state.lambdaG, state.lambdaH, zeros(ncc)),
                                 fx = state.fx,
                                 gx = state.gx,
                                 Hx = state.Hx,
                                 mu = state.mu,
                                 cx = vcat(state.cx,
                                           state.cGx .+ tb,
                                           state.cHx .+ tb,
                                           phi(state.cGx, state.cHx, r, s, t)),
                                 Jx = vcat(state.Jx,
                                           state.JGx,
                                           state.JHx,
                                           diagm(0 => dp[1:ncc]) * state.JGx
                                         + diagm(0 => dp[ncc+1:2*ncc]) * state.JHx))

 stp_rlx = NLPStopping(rlx, rlx_state, optimality_check = KKT, atol = prec, max_cntrs = stp.meta.max_cntrs)

 return stp_rlx
end

"""
Update of the solve sub-problem structure

update atol in the Stopping.
"""
function _rlx_update!(stp :: NLPStopping, parammpcc :: ParamMPCC)

  rlx = stp.pb

  r,s,t = parammpcc.updaterst(rlx.r, rlx.s, rlx.t)
  tb    = parammpcc.tb(r, s, t)
  prec  = parammpcc.prec_oracle(r, s, t, parammpcc.precmpcc)

 return update_rlx!(rlx, r, s, t, tb)
end


text:
struct HeadersFrame
    is_end_stream::Bool
    is_end_headers::Bool
    is_priority::Bool
    stream_identifier::UInt32
    exclusive::Union{Nothing,Bool}
    dependent_stream_identifier::Union{Nothing,UInt32}
    weight::Union{Nothing,UInt8}
    fragment::Array{UInt8, 1}
end

==(a::HeadersFrame, b::HeadersFrame) =
    a.is_end_stream == b.is_end_stream &&
    a.is_end_headers == b.is_end_headers &&
    a.is_priority == b.is_priority &&
    a.stream_identifier == b.stream_identifier &&
    ((a.exclusive === nothing) || a.exclusive == b.exclusive) &&
    ((a.exclusive === nothing) || a.dependent_stream_identifier == b.dependent_stream_identifier) &&
    ((a.exclusive === nothing) || a.weight == b.weight) &&
    a.fragment == b.fragment

function decode_headers(header, payload)
    is_end_stream = header.flags & 0x1 == 0x1
    is_end_headers = header.flags & 0x4 == 0x4
    is_priority = header.flags & 0x20 == 0x20

    payload = remove_padding(header, payload)

    if is_priority
        dependent_stream_identifier = (UInt32(payload[1]) << 24 + UInt32(payload[2]) << 16 +
                                       UInt32(payload[3]) << 8 + UInt32(payload[4])) & 0x7fffffff
        exclusive = payload[1] & 0x80 == 0x80
        weight = payload[5]

        return HeadersFrame(is_end_stream, is_end_headers, is_priority, header.stream_identifier,
                           exclusive, dependent_stream_identifier,
                           weight, getindex(payload, 6:length(payload)))
    else
        return HeadersFrame(is_end_stream, is_end_headers, is_priority, header.stream_identifier,
                           nothing, nothing, nothing, payload)
    end
end

function encode_headers(frame)
    typ = HEADERS
    flags = 0x0 | (frame.is_end_stream ? 0x1 : 0x0) |
        (frame.is_end_headers ? 0x4 : 0x0) |
        (frame.is_priority ? 0x20 : 0x0)

    if frame.is_priority
        payload::Array{UInt8, 1} = [ UInt8(frame.dependent_stream_identifier >> 24) & 0x7f;
                                     UInt8(frame.dependent_stream_identifier >> 16 & 0x000000ff);
                                     UInt8(frame.dependent_stream_identifier >> 8 & 0x000000ff);
                                     UInt8(frame.dependent_stream_identifier & 0x000000ff) ]
        payload[1] = frame.exclusive ? (payload[1] | 0x80 ) : payload[1]
        push!(payload, frame.weight)
        append!(payload, frame.fragment)
    else
        payload = frame.fragment
    end

    return wrap_payload(payload, typ, flags, frame.stream_identifier)
end


text:
using Test
using AbstractPlotting
using MakieRipserer
using Ripserer

import GeometryBasics
using Ripserer: chain_element_type

using MakieRipserer:
    ChainPlot, FilteredVertices, FilteredEdges, FilteredTriangles, FilteredChain

@testset "FilteredChain" begin
    @testset "From filtration" begin
        data = [(sin(t), cos(t)) for t in range(0, 2π; length=13)[1:12]]
        filtration = Rips(data)

        chain = FilteredChain(filtration, data)

        vs = chain.vertices
        es = chain.edges
        ts = chain.triangles

        @test length(vs[0.0]) == 12
        @test length(vs[-0.1]) == 0
        @test length(es[0]) == 0
        for _ in 1:20
            @test iseven(length(es[rand()]))
        end
        @test isempty(ts[0])
        @test ts[0] isa GeometryBasics.Mesh
        @test ts[1] isa GeometryBasics.Mesh
    end
end

@testset "plottype" begin
    data = [(sin(t), cos(t)) for t in range(0, 2π; length=13)[1:12]]
    sx = Simplex{2}(1, 1)
    ch = chain_element_type(typeof(sx), Mod{2})
    @test AbstractPlotting.plottype(sx, data) == ChainPlot
    @test AbstractPlotting.plottype([sx], data) == ChainPlot
    @test AbstractPlotting.plottype(ch(sx), data) == ChainPlot
    @test AbstractPlotting.plottype([ch(sx)], data) == ChainPlot
    @test AbstractPlotting.plottype(Rips(data), data) == ChainPlot
end

@testset "convert_arguments" begin
    data = [(sin(t), cos(t)) for t in range(0, 2π; length=13)[1:12]]
    sx = Simplex{2}(1, 1)
    ch = chain_element_type(typeof(sx), Mod{2})
    @test convert_arguments(ChainPlot, sx, data)[1] isa FilteredChain
    @test convert_arguments(ChainPlot, [sx], data)[1] isa FilteredChain
    @test convert_arguments(ChainPlot, ch(sx), data)[1] isa FilteredChain
    @test convert_arguments(ChainPlot, [ch(sx)], data)[1] isa FilteredChain
    @test convert_arguments(ChainPlot, Rips(data), data)[1] isa FilteredChain
end

@testset "No data error" begin
    data = [(sin(t), cos(t)) for t in range(0, 2π; length=13)[1:12]]
    sx = Simplex{2}(1, 1)
    ch = chain_element_type(typeof(sx), Mod{2})
    @test_throws ArgumentError AbstractPlotting.plot(sx)
    @test_throws ArgumentError AbstractPlotting.plot([sx])
    @test_throws ArgumentError AbstractPlotting.plot(ch(sx))
    @test_throws ArgumentError AbstractPlotting.plot([ch(sx)])
    @test_throws ArgumentError AbstractPlotting.plot(Rips(data))
end


text:
####
# Comparing Bounds With Various Shape Constraints
# ARGS[1] i either MAD or CV
# ARGS[2] is outPathStub
####

using DelimitedFiles
include("src.jl")  #load up vopp

#the default values were used to generate data for experiment
function runMADExp(file_out_name;
	S=2, M=1, mu=1, mode=1,
	dev_grid = range(0., stop=.25, length=100))

	#Create file_out and header
	outPath = "$(file_out_name)__$(S)_$(M)_$(mu)_$(mode)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["Dev" "Shape" "isLB" "Bound" "Time"])

	for dev in dev_grid
		#Dumb bound for plotting ease.
		writedlm(f, [dev "None" true 1.0 0.])

		t = @elapsed bound = vopp.vopp_ub_MAD(S, M, dev)[1]
		writedlm(f, [dev "None" false bound t])

		t = @elapsed bound = vopp.vopp_ub_unimodal_MAD(mu, S, M, dev, mode)
		writedlm(f, [dev "Unimodal" false bound t])

		# t = @elapsed bound = vopp.vopp_ub_symmetric_MAD(mu, S, M, dev, mode)
		# writedlm(f, [dev "Symmetric" false bound t])

		t = @elapsed bound = vopp.vopp_lb_unimodal_MAD(mu, S, M, dev, mode, method=:Formula)
		writedlm(f, [dev "UniAnalytic" true bound t])

		t = @elapsed bound = vopp.vopp_lb_unimodal_MAD(mu, S, M, dev, mode, safe_fail=false)[1]
		writedlm(f, [dev "Unimodal" true bound t])

		# t = @elapsed bound = vopp.vopp_lb_symmetric_MAD(mu, S, M, dev, safe_fail=false)[1]
		# writedlm(f, [dev "Symmetric" true bound t])

		flush(f)
	end
	close(f)
end

function runCVExp(file_out_name;
	S=2, M=1, mu=1, mode=1,
	cv_grid = range(0., stop=.5773502691896258, length=100))

	#Create file_out and header
	outPath = "$(file_out_name)__$(S)_$(M)_$(mu)_$(mode)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["CV" "Shape" "isLB" "Bound" "Time"])

	for C in cv_grid
		#add a dumb bound for ease in plotting
		writedlm(f, [C "None" true 1.0 0.])

		t = @elapsed bound = vopp.vopp_ub_CV(S, M, mu, C)
		writedlm(f, [C "None" false bound t])

		t = @elapsed bound = vopp.vopp_ub_unimodal_CV(mu, S, M, C, mode)
		writedlm(f, [C "Unimodal" false bound t])

		# t = @elapsed bound = vopp.vopp_ub_symmetric_CV(mu, S, M, C, mode)
		# writedlm(f, [C "Symmetric" false bound t])

		t = @elapsed bound = vopp.vopp_lb_unimodal_CV(mu, S, M, C, mode)[1]
		writedlm(f, [C "Unimodal" true bound t])

		# t = @elapsed bound = vopp.vopp_lb_symmetric_CV(mu, S, M, C, safe_fail=true)[1]
		# writedlm(f, [C "Symmetric" true bound t])

		flush(f)
	end
	close(f)
end

function runGMExp(file_out_name;
	S=2, M=.9, mu=1, mode=1,
	B_grid = range(.81168, stop=1., length=100))

	#Create file_out and header
	outPath = "$(file_out_name)__$(S)_$(M)_$(mu)_$(mode)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["GM" "Shape" "isLB" "Bound" "Time"])

	for B in B_grid
		#add a dumb bound for ease in plotting
		writedlm(f, [B "None" true 1.0 0.])

		t = @elapsed bound = vopp.vopp_ub_GM(S, M, mu, B, useSep=true)
		writedlm(f, [B "None" false bound t])

		t = @elapsed bound = vopp.vopp_ub_unimodal_GM(mu, S, M, B, mode)
		writedlm(f, [B "Unimodal" false bound t])

		t = @elapsed bound = vopp.vopp_lb_unimodal_GM(mu, S, M, B, mode)[1]
		writedlm(f, [B "Unimodal" true bound t])

		flush(f)
	end
	close(f)
end

function runICExp(file_out_name;
	S=2, M=.9, mu=1., mode=1., phat = .8)
	qmin, qmax = vopp.min_max_q_IC_uni(S, M, mu, mode, phat, 500)
	q_grid = range(qmin, stop=qmax, length=100)

	#Create file_out and header
	outPath = "$(file_out_name)__$(S)_$(M)_$(mu)_$(mode)_$(phat)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["IC" "Shape" "isLB" "Bound" "Time"])

	for q in q_grid
		#add a dumb bound for ease in plotting
		writedlm(f, [q "None" true 1.0 0.])

		t = @elapsed bound = vopp.vopp_ub_IC(S, M, mu, phat, q)
		writedlm(f, [q "None" false bound t])

		t = @elapsed bound = vopp.vopp_ub_unimodal_IC(mu, S, M, phat, q, mode)
		writedlm(f, [q "Unimodal" false bound t])

		t = @elapsed bound = vopp.vopp_lb_unimodal_IC(mu, S, M, phat, q, mode)[1]
		writedlm(f, [q "Unimodal" true bound t])

		flush(f)
	end
	close(f)
end

if ARGS[1] == "MAD"
	runMADExp(ARGS[2])
elseif ARGS[1] == "MAD2"
	runMADExp(ARGS[2], M=.9, dev_grid = range(0., stop=.2368245, length=100))
elseif ARGS[1]== "CV"
	runCVExp(ARGS[2])
elseif ARGS[1] == "CV2"
	runCVExp(ARGS[2], M=.9, cv_grid = range(0., stop=.5477125, length=100))
elseif ARGS[1] == "GM"
	runGMExp(ARGS[2])
elseif ARGS[1] == "IC"
	runICExp(ARGS[2])
end


text:
module InstrumentConfig

using Downloads

export
    initialize,
    terminate

using YAML

mutable struct Config
   name::String
   example::String
   file_locations::Array
   config::Union{Dict, Nothing}
   loaded_file::String

   function Config(name; example = "")
       dot = new()
       dot.name = name
       dot.example = example
       dot.file_locations = [joinpath(i, name) for i in [pwd(), homedir()]]
       dot.config = nothing
       dot.loaded_file = ""
       return dot
   end

   function Config(name, mod; example = name)
       module_path = pathof(mod)
       pkg_dir = isnothing(module_path) ? pwd() : dirname(dirname(module_path))
       example_path = joinpath(pkg_dir, example)
       return Config(name; example=example_path)
   end
end

function create_config(config; dir=homedir())
    file_path = joinpath(dir, config.name)
    # Local
    if !occursin("http", config.example)
        if Sys.iswindows()
            run(`powershell cp $(config.example) $(file_path)`)
        else
            run(`cp $(config.example) $(file_path)`)
        end
        load_config(config)
        return
    end
    if isempty(config.example)
        @info """
        No example config specified creating empty config file at:
            $file_path
        """
        Base.run(`touch $file_path`)
    else
        Downloads.download(config.example, file_path)
    end
    load_config(config)
end

function edit_config(config)
    config.config == nothing && error("""
    No config loaded.
    Cannot edit a config file that doesn't exists.

    Create a new empty config from your shell by running:
        touch $(config.name)

    Or

    To load the latest config use the create_config function
    """)
    editor, file = ENV["EDITOR"], config.loaded_file
    if Sys.iswindows()
        Base.run(`powershell "$editor $file"`)
    else
        Base.run(`$editor $file`)
    end
    load_config(config)
end

function load_config(config)
    for f in config.file_locations
        if isfile(f)
            config.loaded_file = f
            break
        end
    end
    if isempty(config.loaded_file)
        config.config = Dict()
        @info "No configuration file found:\n$(config.file_locations)"
        return
    end
    try
        config.config = YAML.load_file(config.loaded_file)
        @info "$(config.loaded_file) ~ loaded"
    catch
        @info "$(config.loaded_file) is empty"
        config.config = Dict()
    end
end

function get_config(config)
    if config.config == nothing
        load_config(config)
    end
    return config.config
end

initialize(instr) =
    error("Initialize not implemented for $(typeof(instr))")

terminate(instr) =
    error("Terminate not implemented for $(typeof(instr))")

end


text:
using UltraDark
using PencilFFTs
using NPZ

function add_soliton(grids, mass, position, velocity, phase, t0)

    profile = npzread(joinpath(@__DIR__, "initial_f.npy"))
    delta_x = 0.00001
    alpha = (mass / 3.883)^2
    beta = 2.454

    if typeof(grids) <: PencilGrids
        ψ_glob = global_view(grids.ψx)
    elseif typeof(grids) <: Grids
        ψ_glob = grids.ψx
    else
        throw("Unrecognised grids type")
    end

    for I in CartesianIndices(ψ_glob)

        i, j, k = Tuple(I)  # unpack indices

        x = grids.x[i]
        y = grids.y[j]
        z = grids.z[k]

        vx = velocity[1]
        vy = velocity[2]
        vz = velocity[3]

        distfromcentre = ((x - position[1])^2 + (y - position[2])^2 + (z - position[3])^2)^0.5

        if alpha^0.5 * distfromcentre <= 5.6
            f = alpha * profile[trunc(Int, alpha^0.5 * (distfromcentre / delta_x + 1))]
            f *= exp(im * (alpha * beta * t0 + (x * vx + y * vy + z * vz) - 0.5 * (vx^2 + vy^2 + vz^2) * t0 ))
            f *= exp(im * phase)
            ψ_glob[I] += f
        else
            ψ_glob[I] += 0
        end
    end
end


text:
function Base.show(io::IO, results::Results{Fl}) where Fl
    if !isempty(results)
        println(io, "                         Results                        ")
        println(io, "========================================================")
        println(io, "Number of observations:       ", results.num_observations)
        println(io, "Number of unknown parameters: ", results.num_hyperparameters)
        println(io, "Log-likelihood:               ", @sprintf("%.3f", results.llk))
        println(io, "AIC:                          ", @sprintf("%.3f", results.aic))
        println(io, "BIC:                          ", @sprintf("%.3f", results.bic))
        print_coef_table(io, results.coef_table)
    else
        # Up to discussion of what should be the behaviour
        error("You must call fit! to get results.")
    end
    return nothing
end

function Base.show(io::IO, model::StateSpaceModel)
    print(io, typeof(model), " model")
    return nothing
end

function print_coef_table(io::IO, coef_table::CoefficientTable{Fl}) where Fl
    println(io, "--------------------------------------------------------")
    println(io, "Parameter      Estimate   Std.Error     z stat   p-value")
    offset = 1
    for i in 1:length(coef_table)
        p_c, p_std, p_z_stat, p_p_val = print_coef(coef_table, offset)
        println(io, build_print(p_c, p_std, p_z_stat, p_p_val, coef_table.names[offset]))
        offset += 1
    end
    return nothing
end

function print_coef(coef_table::CoefficientTable{Fl}, offset::Int) where Fl
    p_c = @sprintf("%.4f", coef_table.coef[offset])
    p_std = if isnan(coef_table.std_err[offset])
        " - "
    else
        @sprintf("%.3f", coef_table.std_err[offset])
    end
    p_z_stat = isnan(coef_table.z[offset]) ? " - " : @sprintf("%.3f", coef_table.z[offset])
    p_p_val = if isnan(coef_table.p_value[offset])
        " - "
    else
        @sprintf("%.3f", coef_table.p_value[offset])
    end
    return p_c, p_std, p_z_stat, p_p_val
end

function build_print(p_c, p_std, p_z_stat, p_p_val, param)
    p = ""
    p *= param
    p *= " "^max(0, 23 - length(p_c) - length(param))
    p *= p_c
    p *= " "^max(0, 12 - length(p_std))
    p *= p_std
    p *= " "^max(0, 11 - length(p_z_stat))
    p *= p_z_stat
    p *= " "^max(0, 10 - length(p_p_val))
    p *= p_p_val
    return p
end


text:
using Bukdu # ApplicationController Conn routes get
using HTTP.Messages: setheader

struct RESTController <: ApplicationController
    conn::Conn
end

function index(c::RESTController)
    setheader(c.conn.request.response, "Access-Control-Allow-Origin" => "*")
    render(JSON, 42)
end

routes() do
    get("/", RESTController, index)
end

Bukdu.start(8080)


text:
using Jl2Py
using Documenter

DocMeta.setdocmeta!(Jl2Py, :DocTestSetup, :(using Jl2Py); recursive=true)

makedocs(;
    modules=[Jl2Py],
    authors="Gabriel Wu <wuzihua@pku.edu.cn> and contributors",
    repo="https://github.com/lucifer1004/Jl2Py.jl/blob/{commit}{path}#{line}",
    sitename="Jl2Py.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://lucifer1004.github.io/Jl2Py.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/lucifer1004/Jl2Py.jl",
    devbranch="main",
)


text:
#eigfreesum.jl
#Algorithm 7.1 of Random Eigenvalues by Alan Edelman

#Experiment:    n by n diagonal matrices with random entries +1 or -1
#Plot:          classical and free sums
#Theory:        arcsine distribution

## Parameters
t = 1           # number of trials
n = 1000        # size of matrix
dx = .1

## Experiment
function eigfreesum_experiment(t, n, dx)
    v = Float64[]
    w = Float64[]
    for i = 1:t
        a = sign(randn(n))
        b = sign(randn(n))
        Q = full(QR(randn(n,n))[:Q])                    # random orthogonal matrix
        append!(v, a+b)                                 # classical sum
        append!(w, eigvals(diagm(a) + Q'*diagm(b)*Q))   # free sum
    end
    x = -2-dx/2:dx:2+dx/2
    _,countw = hist(w, x)
    _,countv = hist(v, x)
    x, countw, countv
end
grid, countw, countv = eigfreesum_experiment(t, n, dx)

## Theory
x = -2+dx/20:dx/10:2-dx/20
y = (1/pi) ./ sqrt(4 - x.^2)

## Plot
using Winston
p = FramedPlot()
hw = Histogram(countw*pi/(n*t*dx), step(grid), "color", "blue")
hw.x0 = first(grid)
add(p, hw)
hv = Histogram(countv*pi/(n*t), step(grid), "color", "green")
hv.x0 = first(grid)
add(p, hv)
add(p, Curve(x, y, "color", "red", "linewidth", 2))
if isinteractive()
    Winston.display(p)
else
    file(p, "eigfreesum.png")
end


text:
# This file contains function to generate an example
# of a file containing the generator parameters

"""
    example(file::AbstractString)

Generate an example of a file containing the generator parameters.

# Usage

```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.Gen.example(file)

# output

```
"""
function example(file::AbstractString)

    # Strip the string
    file = strip(file)

    # Check is it's a directory
    if isdir(file)
        throw(ScatsGenIsADir(file))
    end

    open(file, "w") do io
        println(io,
        """
        Sample size
        230

        Sample step
        1.0

        Significance level
        0.01

        Parameter α of a linear trend
        0.1

        Parameter β of a linear trend
        0.05

        Number of harmonics
        1

        Amplitudes array
        1.0

        Frequencies array
        0.1

        Phase shifts array
        0.0

        «Signal-to-noise» ratio
        0.50"""
        )
    end

end


text:
module TurbulenceStatistics

using Oceananigans
using Oceananigans.AbstractOperations
using Oceananigans.Architectures
using Oceananigans.Diagnostics
using Oceananigans.Grids
using Oceananigans.Fields
using Oceananigans.Operators
using Oceanostics

ViscousDissipation(model; data=nothing) =
    Oceanostics.IsotropicViscousDissipationRate(model, model.velocities..., model.diffusivities.νₑ)

include("first_through_third_order.jl")
include("subfilter_fluxes.jl")
include("turbulent_kinetic_energy_budget.jl")

end # module


text:
using PyPlot, PyCall, IJulia

"""
    t = delaunay(p)

Delaunay triangulation `t` of N x 2 node array `p`.
"""
function delaunay(p)
    tri = pyimport("matplotlib.tri")
    t = tri[:Triangulation](p[:,1], p[:,2])
    return Int64.(t[:triangles] .+ 1)
end

"""
    edges, boundary_indices = all_edges(t)

Find all unique edges in the triangulation `t` (ne x 2 array)
Second output is indices to the boundary edges.
"""
function all_edges(t)
    etag = vcat(t[:,[1,2]], t[:,[2,3]], t[:,[3,1]])
    etag = hcat(sort(etag, dims=2), 1:3*size(t,1))
    etag = sortslices(etag, dims=1)
    dup = all(etag[2:end,1:2] - etag[1:end-1,1:2] .== 0, dims=2)[:]
    keep = .![false;dup]
    edges = etag[keep,1:2]
    emap = cumsum(keep)
    invpermute!(emap, etag[:,3])
    emap = reshape(emap,:,3)
    dup = [dup;false]
    dup = dup[keep]
    bndix = findall(.!dup)
    return edges, bndix, emap
end

"""
    e = boundary_nodes(t)

Find all boundary nodes in the triangulation `t`.
"""
function boundary_nodes(t)
    edges, boundary_indices = all_edges(t)
    return unique(edges[boundary_indices,:][:])
end

"""
    tplot(p, t, u=nothing)

If `u` == nothing: Plot triangular mesh with nodes `p` and triangles `t`.
If `u` == solution vector: Plot filled contour color plot of solution `u`.
"""
function tplot(p, t, u=nothing)
    clf()
    axis("equal")
    if u == nothing
        tripcolor(p[:,1], p[:,2], t .- 1, 0*t[:,1],
                  cmap="Set3", edgecolors="k", linewidth=1)
    else
        tricontourf(p[:,1], p[:,2], t .- 1, u, 20)
    end
    draw()
end

"""
    inside = inpolygon(p, pv)

Determine if each point in the N x 2 node array `p` is inside the polygon
described by the NE x 2 node array `pv`.
"""
function inpolygon(p::Array{Float64,2}, pv::Array{Float64,2})
    path = pyimport("matplotlib.path")
    poly = path[:Path](pv)
    inside = [poly[:contains_point](p[ip,:]) for ip = 1:size(p,1)]
end

function remove_outside_tris(p, t, pv)
    pmid = dropdims(sum(p[t,:], dims=2), dims=2) / 3
    is_inside = inpolygon(pmid, pv)
    t = t[is_inside,:]
end

function triarea(p, t)
    d12 = @. p[t[:,2],:] - p[t[:,1],:]
    d13 = @. p[t[:,3],:] - p[t[:,1],:]
    @. abs(d12[:,1] * d13[:,2] - d12[:,2] * d13[:,1]) / 2
end

function remove_tiny_tris(p, t)
    t = t[triarea(p,t) .> 1e-14,:]
end

function circumcenter(p)
    dp1 = @. p[2,:] - p[1,:]
    dp2 = @. p[3,:] - p[1,:]

    mid1 = @. (p[1,:] + p[2,:]) / 2
    mid2 = @. (p[1,:] + p[3,:]) / 2

    s = [ -dp1[2] dp2[2]; dp1[1] -dp2[1]] \ (-mid1 .+ mid2)
    pc = @. mid1' + s[1] * [-dp1[2] dp1[1]]
end

function edge_midpoints(p, t)
    pmid = reshape(p[t,:] + p[t[:,[2,3,1]],:], :, 2) / 2
    pmid = unique(pmid, dims=1)
end

function pmesh(pv, hmax, nref)
    p = zeros(Float64, 0, 2)
    for i = 1:size(pv,1) - 1
        pp = pv[i:i+1,:]
        L = sqrt(sum(diff(pp, dims=1).^2, dims=2))[1]
        if L > hmax
            n = ceil(Int, L / hmax)
            ss = (0:n) / n
            pp = [1 .- ss ss] * pp
        end
        p = [p; pp[1:end-1,:]]
    end

    t = zeros(Int64, 0, 3)
    while true
        t = delaunay(p)
        t = remove_tiny_tris(p, t)
        t = remove_outside_tris(p, t, pv)
        # tplot(p,t), pause(1e-3)

        area = triarea(p, t)
        maxarea, ix = findmax(area)
        if maxarea < hmax^2 / 2
            break
        end
        pc = circumcenter(p[t[ix,:],:])
        p = [p; pc]
    end

    for iref = 1:nref
        p = [p; edge_midpoints(p, t)]
        t = delaunay(p)
        t = remove_tiny_tris(p, t)
        t = remove_outside_tris(p, t, pv)
        # tplot(p, t), pause(1e-3)
    end

    e = boundary_nodes(t)
    p, t, e
end

function pmesh_test()
    pv = [0 0; 1 0; .5 .5; 1 1; 0 1; 0 0]
    p,t,e = pmesh(pv, 0.2, 1);
    tplot(p,t)
end

using SparseArrays
using LinearAlgebra
function coefficients(RHS, p, t, k)
    LHS = [1 p[t[k,1],1] p[t[k,1],2]; 1 p[t[k,2],1] p[t[k,2],2]; 1 p[t[k,3],1] p[t[k,3],2]]
    coeff = LHS\RHS'
    return coeff'
end

function element(coeffmatrix, a, b)
    return coeffmatrix[a,2]*coeffmatrix[b,2] + coeffmatrix[a,3]*coeffmatrix[b,3]
end

function fempoi(p, t, e)
    areas = triarea(p,t)
    n = size(p)[1]
    A = spzeros(n, n); b = zeros(n);
    for k = 1:size(t)[1]
        cm = [coefficients([1 0 0] , p, t, k) ;
              coefficients([0 1 0] , p, t, k) ;
              coefficients([0 0 1] , p, t, k) ]
        insert = [element(cm, 1, 1) element(cm, 1, 2) element(cm, 1, 3);
                  element(cm, 2, 1) element(cm, 2, 2) element(cm, 2, 3);
                  element(cm, 3, 1) element(cm, 3, 2) element(cm, 3, 3)]
        insert *=  areas[k]

        A[t[k,:],t[k,:]] = A[t[k,:],t[k,:]]+insert
        #aa = areas[k]/3.0
        b[t[k,:]] = b[t[k,:]].+ areas[k]/3.0
    end
    for k = 1:size(e)[1]
        i=e[k]
        A[i,:].=0.0
        A[i,i]=1.0
        b[i]=0.0
    end
    dropzeros!(A)
    println(b)
    println(A)

    return A \ b
end

function poiconv(pv, hmax, nrefmax)
   pmax, tmax, emax = pmesh(pv, hmax, nrefmax)
   umax = fempoi(pmax, tmax, emax)
   errors = zeros(nrefmax)
   for ncur = 0:(nrefmax-1)
       pprox, tprox, eprox = pmesh(pv, hmax, ncur)
       uprox = fempoi(pprox, tprox, eprox)
       errorarray = maximum(abs.(umax[1:size(uprox)[1]] - uprox))
       #errors = [errors ; errorarray]
       errors[ncur+1]=errorarray
   end

   return errors
end

function p2mesh(p,t)
    "npoints = edge_midpoints(p,t)"
    "p2 = [p; npoints]"
    "e2 = zeroes()"
    m = size(p)[1]

    elist, bedges, emap = all_edges(t)
    mid_point=Dict()
    for k = 1:size(elist)[1]
        x=(p[elist[k,1],1]+p[elist[k,2],1])/2.0
        y=(p[elist[k,1],2]+p[elist[k,2],2])/2.0
        p=[p; x y]
        mid_point[elist[k,2],elist[k,1]]=size(p)[1]
        mid_point[elist[k,1],elist[k,2]]=size(p)[1]
    end
    t2=zeros(Int64,0,6)
    for i = 1:size(t)[1]
       t12=mid_point[t[i,1],t[i,2]]
       t23=mid_point[t[i,2],t[i,3]]
       t31=mid_point[t[i,3],t[i,1]]
       t2=[t2; t[i,:]' t12 t23 t31]
    end

    e2 = boundary_nodes(t)
    for i=1:size(bedges)[1]
       e=elist[bedges[i],:]
       md=mid_point[e[1],e[2]]
       e2=vcat(e2,md)
    end
    return p,t2,e2
end

using SymPy

function gquad(k, i, j, cmx, p2, t2, ak)
    f(x,y) = 4x^2*(cmx[5,i]*cmx[5,j]) + 2x*y*(cmx[5,i]*cmx[4,j] + cmx[5,j]*cmx[4,i]) + y^2*(cmx[4,i]*cmx[4,j]) + 2x*(cmx[5,i]*cmx[2,j] + cmx[5,j]*cmx[2,i]) + y*(cmx[2,i]*cmx[4,j] + cmx[2,j]*cmx[4,i]) + cmx[2,i]*cmx[2,j]
    g(x,y) = 4y^2*(cmx[6,i]*cmx[6,j]) + 2x*y*(cmx[6,i]*cmx[4,j] + cmx[6,j]*cmx[4,i]) + x^2*(cmx[4,i]*cmx[4,j]) + 2y*(cmx[6,i]*cmx[3,j] + cmx[6,j]*cmx[3,i]) + x*(cmx[3,i]*cmx[4,j] + cmx[3,j]*cmx[4,i]) + cmx[3,i]*cmx[3,j]

    result = 0
    for n = 1:3
	base = (p2[t2[k,1],:] + p2[t2[k,2],:] + p2[t2[k,3],:])/6
        num = base + p2[t2[k,n],:]/2
        result = result + f(num[1],num[2]) + g(num[1],num[2])
    end
    result = (result*ak)/3
    return result
end

function gquad2(k, i, cmx, p2, t2, ak)
    h(x,y) = cmx[1,i] + x*cmx[2,i] + y*cmx[3,i] + x*y*cmx[4,i] + x^2*cmx[5,i] + y^2*cmx[6,i]
    result = 0
    for n = 1:3
	base = (p2[t2[k,1],:] + p2[t2[k,2],:] + p2[t2[k,3],:])/6
        num = base + p2[t2[k,n],:]/2
        result = result + h(num[1],num[2])
    end
    result = (result*ak)/3
    return result
end

function coeff2(RHS, p2, t2, k)
    nm = p2[t2[k,1:6],:]
    LHS = zeros(6,6)
    for n = 1:6
       LHS[n,:] = [1 nm[n,1] nm[n,2] nm[n,1]*nm[n,2] nm[n,1]^2 nm[n,2]^2]
    end
    return (LHS\RHS), LHS
end

function fempoi2(p2, t2, e2)
    areas = triarea(p2,t2)
    n = size(p2)[1]
    A = spzeros(n,n); b = zeros(n);

    for k = 1:size(t2)[1]
        insert = zeros(6,6)
        insert2 = zeros(6)
        cmx, lhss = coeff2(Matrix{Float64}(I, 6, 6), p2, t2, k)
        for i = 1:6
            for j = 1:6
                insert[i,j] = gquad(k,i,j,cmx,p2,t2,areas[k])
            end
            insert2[i] = gquad2(k,i,cmx,p2,t2,areas[k])
        end
        A[t2[k,:],t2[k,:]] = A[t2[k,:],t2[k,:]]+insert
        b[t2[k,:]] = b[t2[k,:]]+ insert2
    end
    for k = 1:size(e2)[1]
        i=e2[k]
        A[i,:].=0.0
        A[i,i]=1.0
        b[i]=0.0
    end

    return A\b

end

function convtest(pv, hmax, nrefmax)
    p, t, e = pmesh(pv,hmax,nrefmax)
    p2, t2, e2 = p2mesh(p,t)
    umax = fempoi2(p2,t2,e2)
    errors = zeros(nrefmax)
    for ncur = 0:(nrefmax-1)
       pprox, tprox, eprox = pmesh(pv, hmax, ncur)
       pprox2, tprox2, eprox2 = p2mesh(pprox, tprox)
       uprox = fempoi2(pprox2, tprox2, eprox2)
       errorarray = maximum(abs.(umax[1:size(uprox)[1]] - uprox))
       #errors = [errors ; errorarray]
       errors[ncur+1]=errorarray
    end

    return errors
end

"errors = convtest(pv,hmax,nrefmax)

clf()
loglog(hmax ./ [1,2,4,8], errors)
rates = @. log2(errors[end-1,:]) - log2(errors[end,:])

@. log2(errors[end-1,:])"


text:
module TestFindminmax
test() = include(joinpath(@__DIR__, "examples/findminmax.jl"))
end  # module


text:
"""
    w = ds_hann(n)

A Hann window of length `n`. Does not smear tones located exactly in a bin.
"""
function ds_hann(n)
    return 0.5 * (1 .- cos.(2*pi*(0:n-1)/n))
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
using Test
using DotNET
import DotNET:unbox

@testset "Type Loader" begin
    @test T"System.Int64".Name == "Int64"
    @test T"System.String".Name == "String"
    @test T"System.Console, mscorlib".Name == "Console"
    @test_throws CLRException T"BadType.KUUPXXCL"
end # testset

@testset "Marshaller" begin
    uturn(x) = @test (unbox(convert(CLRObject, x)) === x);true
    uturn(Int8(42))
    uturn(UInt8(42))
    uturn(Int16(42))
    uturn(UInt16(42))
    uturn(Int32(42))
    uturn(UInt32(42))
    uturn(Int64(42))
    uturn(UInt64(42))
    uturn(Float32(42))
    uturn(Float64(42))
    uturn('A')
    uturn("Hello, World!")
    array = convert(CLRObject, [1,2,3])
    @test collect(array) == [1,2,3]
    for (i, x) in enumerate(array)
        @test x == i
    end
end

@testset "Member Invocation" begin
    ArrayList = T"System.Collections.ArrayList, mscorlib"
    li = ArrayList.new()
    li.Add(42)
    @test li.Count == 1
    ListT = T"System.Collections.Generic.List`1, mscorlib"
    li = ListT.new[T"System.Int64"]()
    li.Add(42)
    @test li.Count == 1
    @test isclrtype(T"System.Array".Empty[T"System.Int64"](), T"System.Int64[]")
    Console = T"System.Console, mscorlib"
    oldpos = Console.CursorTop
    Console.CursorTop = Int32(0)
    @test Console.CursorTop == Int32(0)
    Console.CursorTop = oldpos
    WeakReference = T"System.WeakReference, mscorlib"
    ref = WeakReference.new(CLRObject(0))
    ref.Target = WeakReference
    @test T"System.Object, mscorlib".Equals(ref.Target, WeakReference)
end

@testset "Type System" begin
    li = convert(CLRObject, Int64[1])
    @test eltype(collect(li)) == Int64
    List = T"System.Collections.Generic.List`1"
    li = List.new[T"System.Int64"]()
    li.Add(Int64(1))
    @test eltype(collect(li)) == Int64
end


text:
function test12()
  filename = "test12.asc"
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)

  v1steps = [1.0, 1.1487, 1.31951, 1.51572, 1.7411,
            2.0, 2.2974, 2.63902, 3.03143, 3.4822,
            4.0, 4.59479, 5.27803, 6.06287, 6.9644,
            8.0, 9.18959, 10.5561, 12.1257, 13.9288,
            16, 18.3792, 20]
  bsteps = [1.0, 3.0, 5.0, 7.0, 9.0, 10.0]
  csteps = [4.0, 5.0, 6.0]

  @test stepvalues(sim) == (v1steps, bsteps, csteps)
  @test measurementnames(sim) == ()
  @test stepnames(sim) == ("V1","b","c")
  #@test measurementvalues(sim) == [23.0,6.0,3.0,0.0]

  @test length(collect(perlineiterator(sim))) == length(v1steps)*length(bsteps)*length(csteps)

  pli = perlineiterator(sim,header=true)
  (header,state) = iterate(pli)
  @test header == ("V1","b","c","a")
  show(IOBuffer(),sim)
end
test12()


text:
using LearnedImageBinarization
using Test

@testset "LearnedImageBinarization.jl" begin
    # Write your own tests here.
end


text:
"""
Create a named tuple observation.
"""
observation_from_full_traj(times::Vector{Float64}, states::Vector{Int64}) = (y = sum(times),a = last(states))

"""
QQQQ
"""
function sufficient_stat_from_trajectory(d::MAPHDist, sojourn_times::Array{Float64}, states::Array{Int})::MAPHSufficientStats
    p,q = model_size(d)
    transient_states = (q+1):(q+p)

    ss =  MAPHSufficientStats(d)

    for s = 1:p
        if states[1] == transient_states[s]
            ss.B[s] +=1
        end
    end

    for i = 1:(length(states)-1)
        ss.Z[states[i]-q] += sojourn_times[i]
        ss.N[states[i]-q, states[i+1]] += 1
    end

    return ss
end

"""
QQQQ
"""
#Maybe rename this function to "maximization step..." or similar
function maximum_likelihood_estimate(p::Int, q::Int, ss::MAPHSufficientStats)
    α_next = max.(ss.B,0)
    t_next = max.(ss.N[:,1:q] ./ ss.Z,0)

    t_next[isnan.(t_next)] .= 0

    T_next = zeros(p,p)
    T0_next = zeros(p,q)

    for i = 1:p
        T_next[i,:] = max.(ss.N[i,(q+1):(q+p)]./ss.Z[i],0)
        T_next[i,isnan.(T_next[i,:])].=0
        T_next[i,i] = -(sum(t_next[i,:]) + sum(T_next[i,:])) #-(t_next[i]+sum(T_next[i,:]))
        T0_next[i,:] = max.(ss.N[i,1:q]./ss.Z[i],0)
        T0_next[i,isnan.(T0_next[i,:])].=0
    end
    #################

    T = ss.N[:,(q+1):(q+p)]./ss.Z
    T0 = ss.N[:,1:q]./ss.Z
    for i = 1:p
        T[i,i] += -sum(T[i,:])-sum(T0[i,:])
    end

    display(T0)
    display(T0_next)

    return α_next, T_next, T0_next
end

"""
Fits ... QQQQ
"""
function fit_maph(times::Vector{Float64}, absorbing_states::Vector{Int}, p::Int; max_iter::Int = 100)
    @assert length(times) == length(absorbing_states) "Vector of times and absorbing states mismatch in length"
    unique_absorbing_states = unique(absorbing_states)
    q = length(unique_absorbing_states)
    @assert minimum(unique_absorbing_states) == 1 && maximum(unique_absorbing_states) == q  "Mismatch with absorbing states"
    n = length(times)
    print("Fitting MAPH with p=$p hidden states, q=$q absorbing states, and n=$n observations")

    ds = compute_descriptive_stats(times, absorbing_states)

    @show ds
    dist = MAPHDist(p, ds...) #Moment based initilization...

    @show absorption_probs(dist)

    iter = 1
    while iter < max_iter

        #QQQQ - refactor next couple of lines...

        ## Expectation steps
        computed_ss = MAPHSufficientStats[]
        for i = 1:n
            obs = (y = times[i], a = absorbing_states[i])
            ss_i = sufficient_stats(obs,dist)
            push!(computed_ss, ss_i)
        end

        ## Maximisation steps
        mean_ss = mean(computed_ss)

        α_next, T_next, T0_next = maximum_likelihood_estimate(p, q, mean_ss)
        #remove numerical instabilities

        #@show T_next, T0_next

        α_next = max.(α_next, 0) #QQQQ check why - maybe goes slightly negative
        α_next /= sum(α_next)
        dist = MAPHDist(α_next', T_next, T0_next)

        iter +=1
    end

    return dist
end

"""
Computes some descriptive statistics of a data sequence
"""
function compute_descriptive_stats(times::Vector{Float64}, abosrbing_states::Vector{Int})
    #QQQQ - good to have a test function for this function (in tests)
    @assert length(times) == length(abosrbing_states) "Vector of times and absorbing states mismatch in length"
    unique_absorbing_states = unique(abosrbing_states)
    q = length(unique_absorbing_states)
    @assert minimum(unique_absorbing_states) == 1 && maximum(unique_absorbing_states) == q  "Mismatch with absorbing states"
    n = length(times)

    probs = [mean(abosrbing_states .== i) for i in 1:q]
    means = [mean(times[abosrbing_states .== i]) for i in 1:q]
    scvs = [var(times[abosrbing_states .== i]) for i in 1:q] ./ (means .^ 2)

    @show probs

    return probs, means, scvs
end


text:
using CompScienceMeshes
using SauterSchwabQuadrature

pI = point(1,5,3)
pII = point(2,5,3)
pIII = point(7,1,0)
pIV = point(5,1,-3)
pV = point(0,0,0)

Sourcechart = simplex(pI,pIII,pII)
Testchart = simplex(pV,pIV,pI)

Accuracy = 12

function integrand(x,y)
			return(((x-pI)'*(y-pV))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
 end

result = sauterschwabintegral(Sourcechart, Testchart, integrand, Accuracy, Accuracy)
println(result)

#=For those who want to test the sauterschwab_nonparameterized() function,
may uncomment the following five lines=#

#sourcechart = simplex(pI,pIII,pII)
#testchart = simplex(pI,pIV,pV)
#cv = CommonVertex(Accuracy)
#result2 = sauterschwab_nonparameterized(sourcechart, testchart, integrand, cv)
#println(result2)

#=The common vertex is the first input argument of both simplex() functions,
hence the required condition is fulfilled. The user may also compare the two
results and see that both are equal=#


text:
include("..\\libs\\FunctionsRobot.jl")

function putmarkerInCorner!(r::Robot)
    for side in instances(HorizonSide)
        putmarker!(r)
        while !isborder(r, side)
            move!(r, side)
        end
    end
end

function ugly!(r::Robot)
    sides = moveAndReturnDirections!(r)
    putmarkerInCorner!(r)
    moveToStartplace!(r)
    moveToBeginplace!(r, sides)

end


text:
# arf_get_d and arf_get_si clash
get_d(x::ArfLike; rnd::Union{arb_rnd,RoundingMode} = RoundNearest) = get_d(x, rnd)
function get_d(x::ArfLike, rnd::Union{arb_rnd,RoundingMode})
    ccall(@libarb(arf_get_d), Float64, (Ref{arf_struct}, arb_rnd), x, rnd)
end
get_si(x::ArfLike; rnd::Union{arb_rnd,RoundingMode} = RoundNearest) = get_si(x, rnd)
function get_si(x::ArfLike, rnd::Union{arb_rnd,RoundingMode})
    ccall(@libarb(arf_get_si), Int, (Ref{arf_struct}, arb_rnd), x, rnd)
end