text:
function onedBilinear( kernel, ϕ, test, wg )
#function onedBilinear( kernel:: Array(Float64,1),
#                       ϕ     :: Array(Float64,2),
#                       test  :: Array(Float64,2),
#                       wg    :: Array(Float64,1) )
#= """
  onedBilinear - routine to compute (kernel*ϕ,test) = \int{ kernel*ϕ*test }
                 returns a matrix,
                 rows correspond to different test functions,
                 columns correspond to different local unknowns (in ϕ)

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    M = onedBilinear(kernel, ϕ, test, wg)
  ```

  Arguments:
  - `kernel`: Kernel function in the integral evaluated at the Gauss points
  - `ϕ`: Matrix of element test functions evaluated at the Gauss points
           (dim: n_gauss, n_dof)
  - `test`: Matrix of test functions evaluated at the Gauss points
            (dim: n_gauss, n_test)
  - `wg`: Column vector of Gauss weights

  Output arguments:
  - `M`: the current element's contribution to the system matrix
""" =#

  wk = wg.*kernel

  M = test'*(LinearAlgebra.Diagonal(wk)*ϕ)

  return M

end


text:
using Libdl
function feature_string()
    llvmlib_path = VERSION ≥ v"1.6.0-DEV.1429" ? Base.libllvm_path() : only(filter(lib->occursin(r"LLVM\b", basename(lib)), Libdl.dllist()))
    libllvm = Libdl.dlopen(llvmlib_path)
    gethostcpufeatures = Libdl.dlsym(libllvm, :LLVMGetHostCPUFeatures)
    features_cstring = ccall(gethostcpufeatures, Cstring, ())
    features = filter(ext -> (m = match(r"\d", ext); isnothing(m) ? true : m.offset != 2 ), split(unsafe_string(features_cstring), ','))
    features, features_cstring
end

# const FEATURE_DICT = Dict{String,Bool}()
# has_feature(str) = get(FEATURE_DICT, str, false)

archstr() = Sys.ARCH === :i686 ? "x86_64_" : string(Sys.ARCH) * '_'

# feature_name(ext) = archstr() * replace(ext[2:end], r"\." => "_")
feature_name(ext) = Symbol(archstr() * ext[2:end])
process_feature(ext) = (feature_name(ext), first(ext) == '+')

has_feature(_) = False()
function set_featue(feature::Symbol, has::Bool)
    featqn = QuoteNode(feature)
    if has
        @eval has_feature(::Val{$featqn}) = True()
    else
        @eval has_feature(::Val{$featqn}) = False()
    end
end
function set_features!()
    features, features_cstring = feature_string()
    for ext ∈ features
        feature, has = process_feature(ext)
        set_featue(feature, has)
    end
    Libc.free(features_cstring)
end
set_features!()

function reset_features!()
    features, features_cstring = feature_string()
    for ext ∈ features
        feature, has = process_feature(ext)
        if Bool(has_feature(Val(feature)))::Bool !== has
            @info "Defining $(has ? "presence" : "absense") of feature $feature."
            set_featue(feature, has)
        end
    end
    Libc.free(features_cstring)
end

register_size() = ifelse(
    has_feature(Val(:x86_64_avx512f)),
    StaticInt{64}(),
    ifelse(
        has_feature(Val(:x86_64_avx)),
        StaticInt{32}(),
        StaticInt{16}()
    )
)
simd_integer_register_size() = ifelse(
    has_feature(Val(:x86_64_avx2)),
    register_size(),
    ifelse(
        has_feature(Val(:x86_64_sse2)),
        StaticInt{16}(),
        StaticInt{8}()
    )
)
fma_fast() = has_feature(Val(:x86_64_fma)) | has_feature(Val(:x86_64_fma4))
if Sys.ARCH === :i686
    register_count() = StaticInt{8}()
else
    register_count() = ifelse(has_feature(Val(:x86_64_avx512f)), StaticInt{32}(), StaticInt{16}())
end
has_opmask_registers() = has_feature(Val(:x86_64_avx512f))

register_size(::Type{T}) where {T} = register_size()
register_size(::Type{T}) where {T<:Union{Signed,Unsigned}} = simd_integer_register_size()

fast_int64_to_double() = has_feature(Val(:x86_64_avx512dq))

function define_cpu_name()
    cpu = QuoteNode(Symbol(Sys.CPU_NAME::String))
    @eval cpu_name() = Val{$cpu}()
end
define_cpu_name()

ext: jl


text:
@enum Verbosity VerbosityTOP VerbosityCLOSEST VerbosityALL

"""
    transfer_casing_for_similar_text(text_w_casing, text_wo_casing)
"""
function transfer_casing_for_similar_text(text_w_casing, text_wo_casing)
    isempty(text_wo_casing) && return text_wo_casing
    isempty(text_w_casing) && throw(ArgumentError("Empty 'text_w_casing', we need to know what casing to transfer"))

    # we will collect the case_text:
    res = ""

    # get the operation codes describing the differences between the
    # two strings and handle them based on the per operation code rules
    transfer(c1, c2) = isuppercase(c1) ? uppercase(c2) : lowercase(c2)
    for (tag, i1, i2, j1, j2) in get_opcodes(lowercase(text_w_casing), text_wo_casing)
        if tag == "insert"
            # if this is the first character and so there is no
            # character on the left of this or the left of it a space
            # then take the casing from the following character
            # otherwise just take the casing from the prior
            # character
            idx = i1 == 1 || text_w_casing[i1 - 1] == ' ' ? i1 : i1 - 1
            res *= transfer(text_w_casing[idx], text_wo_casing[j1:j2])
        elseif tag == "delete"
            # for deleted characters we don't need to do anything
        elseif tag == "equal"
            # for 'equal' we just transfer the text from the
            # text_w_casing, as anyhow they are equal (without the
            # casing)
            res *= text_w_casing[i1:i2]
        elseif tag == "replace"
            w_casing = text_w_casing[i1:i2]
            wo_casing = text_wo_casing[j1:j2]

            # if they are the same length, the transfer is easy
            if length(w_casing) == length(wo_casing)
                res *= transfer_casing_for_matching_text(w_casing, wo_casing)
            else
                # if the replaced has a different length, then we
                # transfer the casing character-by-character and using
                # the last casing to continue if we run out of the
                # sequence
                last_case = 'a'
                for (w, wo) in zip(w_casing, wo_casing)
                    last_case = w
                    res *= transfer(w, wo)
                end
                # once we ran out of 'w', we will carry over
                # the last casing to any additional 'wo'
                # characters
                res *= transfer(last_case, wo_casing[length(w_casing) + 1:end])
            end
        end
    end

    return res
end

"""
    transfer_casing_for_matching_text(text_w_casing, text_wo_casing)

Transferring the casing from one text to another - assuming that
they are 'matching' texts, alias they have the same length.
Parameters
----------
text_w_casing : str
    Text with varied casing
text_wo_casing : str
    Text that is in lowercase only
Returns
-------
str
    Text with the content of `text_wo_casing` and the casing of
    `text_w_casing`
Raises
------
ArgumentError
    If the input texts have different lengths
"""
function transfer_casing_for_matching_text(text_w_casing, text_wo_casing)
    if length(text_w_casing) != length(text_w_casing)
        throw(ArgumentError("The 'text_w_casing' and 'text_wo_casing' don't have the same length,
        so you can't use them with this method, you should be using the more general
        transfer_casing_similar_text() method."))
    end

    return String(collect(isuppercase(x) ? uppercase(y) : lowercase(y) for (x, y) in zip(text_w_casing, text_wo_casing)))
end

# Based on the functionality of the python difflib get_opcodes()
function get_opcodes(s1::S1, s2::S2) where {S1 <: AbstractString, S2 <: AbstractString}
    res = Tuple{String, Int, Int, Int, Int}[]
    blocks = sort!(collect(matching_blocks(s1, s2)), by = x -> x[1])
    i1, i2 = 1, 1
    for block in blocks
        if i1 < block[1] && i2 == block[2]
            push!(res, ("delete", i1, block[1] - 1, i2, block[2] - 1))
        elseif i1 == block[1] && i2 < block[2]
            push!(res, ("insert", i1, block[1] - 1, i2, block[2] - 1))
        elseif i1 < block[1] && i2 < block[2]
            push!(res, ("replace", i1, block[1] - 1, i2, block[2] - 1))
        end
        i1 = block[1] + block[3]
        i2 = block[2] + block[3]
        push!(res, ("equal", block[1], i1 - 1, block[2], i2 - 1))
    end
    if i1 <= length(s1) && i2 > length(s2)
        push!(res, ("delete", i1, length(s1), i2, i2 - 1))
    elseif i1 > length(s1) && i2 <= length(s2)
        push!(res, ("insert", i1, i1 - 1, i2, length(s2)))
    elseif i1 <= length(s1) && i2 <= length(s2)
        push!(res, ("replace", i1, length(s1), i2, length(s2)))
    end

    return res
end


text:
@testset "INT algorithm" begin
    prob, tspan = exponential_1d()
    sol = solve(prob, tspan=tspan, INT(δ=0.01))
    @test isa(sol.alg, INT)
    @test dim(sol) == 1
    @test length(sol) == 100
    @test setrep(sol) <: Interval
    @test setrep(sol) == Interval{Float64,IntervalArithmetic.Interval{Float64}}

    prob, tspan = exponential_1d(invariant=HalfSpace([-1.0], -0.3)) # x >= 0.3
    sol_inv = solve(prob, tspan=tspan, INT(δ=0.01))
    @test length(sol_inv) == 52
    @test [0.3] ∈ sol_inv[end]
    # check that the following reach-set escapes the invariant
    @test [0.3] ∈ sol[52] && [0.3] ∉ sol[53]

    # check NSTEPS option
    sol = solve(prob, NSTEPS=10, INT(δ=0.01))
    @test length(sol) == 10

    # doesn't work for higher dimensional systems
    prob, tspan = linear5D_homog()
    @test_throws ArgumentError solve(prob, tspan=tspan, INT(δ=0.01))
end

# Issue #378
@testset "Scalar affine ODE" begin
    prob = @ivp(x' = x + [1.0], x(0) ∈ Singleton([1.0]))
    sol = solve(prob, T=4.0)
end


text:
for f in ["ufixed.jl", "fixed.jl"]
    println("Testing $f")
    include(f)
end


text:
struct ODEAdjointSensitivityFunction{rateType,uType,uType2,UF,PF,G,JC,GC,A,DG,TJ,PJT,PJC,CP,SType,INT,CV} <: SensitivityFunction
  uf::UF
  pf::PF
  g::G
  J::TJ
  pJ::PJT
  dg_val::uType
  jac_config::JC
  g_grad_config::GC
  paramjac_config::PJC
  alg::A
  f_cache::rateType
  discrete::Bool
  y::uType2
  sol::SType
  dg::DG
  checkpoints::CP
  integrator::INT
  colorvec::CV
end

@noinline function ODEAdjointSensitivityFunction(g,u0,p,alg,discrete,sol,dg,checkpoints,tspan,colorvec)
  numparams = length(p)
  numindvar = length(u0)
  # if there is an analytical Jacobian provided, we are not going to do automatic `jac*vec`
  f = sol.prob.f
  isautojacvec = DiffEqBase.has_jac(f) ? false : get_jacvec(alg)
  J = isautojacvec ? nothing : similar(u0, numindvar, numindvar)

  if !discrete
    if dg != nothing || isautojacvec
      pg = nothing
      pg_config = nothing
    else
      pg = UGradientWrapper(g,tspan[2],p)
      pg_config = build_grad_config(alg,pg,u0,p)
    end
  else
    pg = nothing
    pg_config = nothing
  end

  if DiffEqBase.has_jac(f) || isautojacvec
    jac_config = nothing
    uf = nothing
  else
    uf = DiffEqDiffTools.UJacobianWrapper(f,tspan[2],p)
    jac_config = build_jac_config(alg,uf,u0)
  end

  y = copy(sol(tspan[1])) # TODO: Has to start at interpolation value!
  paramjac_config = nothing
  pf = nothing
  if !isquad(alg)
    if DiffEqBase.has_paramjac(f) || isautojacvec
      paramjac_config = nothing
    else
      pf = DiffEqDiffTools.ParamJacobianWrapper(f,tspan[1],y)
      paramjac_config = build_param_jac_config(alg,pf,y,p)
    end
  end

  pJ = if !isquad(alg)
    isautojacvec ? nothing : similar(sol.prob.u0, numindvar, numparams)
  else
    nothing
  end
  integrator = if ischeckpointing(alg)
    integ = init(sol.prob, sol.alg, save_on=false)
    integ.u = y
    integ
  else
    nothing
  end
  dg_val = similar(u0, numindvar) # number of funcs size
  f_cache = deepcopy(u0)

  return ODEAdjointSensitivityFunction(uf,pf,pg,J,pJ,dg_val,
                               jac_config,pg_config,paramjac_config,
                               alg,f_cache,
                               discrete,y,sol,dg,checkpoints,integrator,colorvec)
end

# u = λ'
function (S::ODEAdjointSensitivityFunction)(du,u,p,t)
  @unpack y, sol, J, uf, alg, f_cache, jac_config, discrete, dg, dg_val, g, g_grad_config = S
  idx = length(y)
  f = sol.prob.f
  isautojacvec = DiffEqBase.has_jac(f) ? false : get_jacvec(alg)

  if isbcksol(alg)
    λ     = @view u[1:idx]
    dλ    = @view du[1:idx]
    grad  = @view u[idx+1:end-idx]
    dgrad = @view du[idx+1:end-idx]
    _y    = @view u[end-idx+1:end]
    dy    = @view du[end-idx+1:end]
    copyto!(vec(y), _y)
    isautojacvec || f(dy, _y, p, t)
  else
    if ischeckpointing(alg)
      @unpack integrator, checkpoints = S
      # assuming that in the forward direction `t0` < `t1`, and the
      # `checkpoints` vector is sorted with respect to the forward direction
      tidx = findlast(x->x <= t, checkpoints)
      t0 = checkpoints[tidx]
      dt = t-t0
      if abs(dt) > integrator.opts.dtmin
        sol(integrator.u, t0)
        copyto!(integrator.uprev, integrator.u)
        integrator.t = t0
        # set `iter` to some arbitrary integer so that there won't be max maxiters error
        integrator.iter=100
        u_modified!(integrator, true)
        step!(integrator, dt, true)
        # `integrator.u` is aliased to `y`
      else
        sol(y,t)
      end
    else
      sol(y,t)
    end
    if isquad(alg)
      λ     = u
      dλ    = du
    else
      λ     = @view u[1:idx]
      dλ    = @view du[1:idx]
      grad  = @view u[idx+1:end]
      dgrad = @view du[idx+1:end]
    end
  end

  if !isautojacvec
    if DiffEqBase.has_jac(f)
      f.jac(J,y,p,t) # Calculate the Jacobian into J
    else
      uf.t = t
      jacobian!(J, uf, y, f_cache, alg, jac_config)
    end
    mul!(dλ',λ',J)
  elseif isquad(alg)
    _dy, back = Tracker.forward(y) do u
      if DiffEqBase.isinplace(sol.prob)
        out_ = map(zero, u)
        f(out_, u, p, t)
        Tracker.collect(out_)
      else
        vec(f(u, p, t))
      end
    end
    dλ[:] = Tracker.data(back(λ)[1])
    isbcksol(alg) && (dy[:] = vec(Tracker.data(_dy)))
  else
    _dy, back = Tracker.forward(y, sol.prob.p) do u, p
      if DiffEqBase.isinplace(sol.prob)
        out_ = map(zero, u)
        f(out_, u, p, t)
        Tracker.collect(out_)
      else
        vec(f(u, p, t))
      end
    end
    dλ[:], dgrad[:] = Tracker.data.(back(λ))
    isbcksol(alg) && (dy[:] = vec(Tracker.data(_dy)))
  end

  dλ .*= -one(eltype(λ))

  if !discrete
    if dg != nothing
      dg(dg_val,y,p,t)
    else
      g.t = t
      gradient!(dg_val, g, y, alg, g_grad_config)
    end
    dλ .+= dg_val
  end

  if !isquad(alg) && !isautojacvec
    @unpack pJ, pf, paramjac_config = S
    if DiffEqBase.has_paramjac(f)
      f.paramjac(pJ,y,sol.prob.p,t) # Calculate the parameter Jacobian into pJ
    else
      jacobian!(pJ, pf, sol.prob.p, f_cache, alg, paramjac_config)
    end
    mul!(dgrad',λ',pJ)
  end
  nothing
end

# g is either g(t,u,p) or discrete g(t,u,i)
@noinline function ODEAdjointProblem(sol,g,t=nothing,dg=nothing,
                           alg=SensitivityAlg();
                           checkpoints=sol.t,
                           callback=CallbackSet())
  f = sol.prob.f
  tspan = (sol.prob.tspan[2],sol.prob.tspan[1])
  discrete = t != nothing

  p = sol.prob.p
  p === DiffEqBase.NullParameters() && error("Your model does not have parameters, and thus it is impossible to calculate the derivative of the solution with respect to the parameters. Your model must have parameters to use parameter sensitivity calculations!")

  u0 = zero(sol.prob.u0)

  len = isquad(alg) ? length(u0) : length(u0)+length(p)
  λ = similar(u0, len)
  sense = ODEAdjointSensitivityFunction(g,u0,
                                        p,alg,discrete,
                                        sol,dg,checkpoints,tspan,f.colorvec)

  init_cb = t !== nothing && sol.prob.tspan[2] == t[end]
  cb = generate_callbacks(sense, g, λ, t, callback, init_cb)
  z0 = isbcksol(alg) ? [vec(zero(λ)); vec(sense.y)] : vec(zero(λ))
  ODEProblem(sense,z0,tspan,p,callback=cb)
end

function generate_callbacks(sensefun, g, λ, t, callback, init_cb)
  if sensefun.discrete
    @unpack alg, y, sol = sensefun
    prob = sol.prob
    cur_time = Ref(length(t))
    function time_choice(integrator)
      cur_time[] > 0 ? t[cur_time[]] : nothing
    end
    affect! = let isq = isquad(alg), λ=λ, t=t, y=y, cur_time=cur_time, idx=length(prob.u0)
      function (integrator)
        p, u = integrator.p, integrator.u
        λ  = isq ? λ : @view(λ[1:idx])
        g(λ,y,p,t[cur_time[]],cur_time[])
        if isq
          u .+= λ
        else
          u = @view u[1:idx]
          u .= λ .+ @view integrator.u[1:idx]
        end
        u_modified!(integrator,true)
        cur_time[] -= 1
      end
    end
    cb = IterativeCallback(time_choice,affect!,eltype(prob.tspan);initial_affect=init_cb)

    _cb = CallbackSet(cb,callback)
  else
    _cb = callback
  end
  return _cb
end

struct AdjointSensitivityIntegrand{pType,uType,rateType,S,AS,PF,PJC,A,PJT}
  sol::S
  adj_sol::AS
  p::pType
  y::uType
  λ::uType
  pf::PF
  f_cache::rateType
  pJ::PJT
  paramjac_config::PJC
  alg::A
end

function AdjointSensitivityIntegrand(sol,adj_sol,alg=SensitivityAlg())
  prob = sol.prob
  @unpack f, p, tspan, u0 = prob
  y = similar(sol.prob.u0)
  λ = similar(adj_sol.prob.u0)
  # we need to alias `y`
  pf = DiffEqDiffTools.ParamJacobianWrapper(f,tspan[1],y)
  f_cache = similar(y)
  isautojacvec = DiffEqBase.has_paramjac(f) ? false : get_jacvec(alg)
  pJ = isautojacvec ? nothing : similar(u0,length(u0),length(p))

  if DiffEqBase.has_paramjac(f) || isautojacvec
    paramjac_config = nothing
  else
    paramjac_config = build_param_jac_config(alg,pf,y,p)
  end
  AdjointSensitivityIntegrand(sol,adj_sol,p,y,λ,pf,f_cache,pJ,paramjac_config,alg)
end

function (S::AdjointSensitivityIntegrand)(out,t)
  @unpack y, λ, pJ, pf, p, f_cache, paramjac_config, alg, sol, adj_sol = S
  f = sol.prob.f
  sol(y,t)
  adj_sol(λ,t)
  λ .*= -one(eltype(λ))
  isautojacvec = DiffEqBase.has_paramjac(f) ? false : get_jacvec(alg)
  # y is aliased
  pf.t = t

  if !isautojacvec
    if DiffEqBase.has_paramjac(f)
      f.paramjac(pJ,y,p,t) # Calculate the parameter Jacobian into pJ
    else
      jacobian!(pJ, pf, p, f_cache, alg, paramjac_config)
    end
    mul!(out',λ',pJ)
  else
    _, back = Tracker.forward(y, p) do u, p
      if DiffEqBase.isinplace(sol.prob)
        out_ = map(zero, u)
        f(out_, u, p, t)
        Tracker.collect(out_)
      else
        vec(f(u, p, t))
      end
    end
    out[:] = vec(Tracker.data(back(λ)[2]))
  end
  out'
end

function (S::AdjointSensitivityIntegrand)(t)
  out = similar(S.p)
  S(out,t)
end

function adjoint_sensitivities_u0(sol,alg,g,t=nothing,dg=nothing;
                                  abstol=1e-6,reltol=1e-3,
                                  iabstol=abstol, ireltol=reltol,sensealg=SensitivityAlg(checkpointing=!sol.dense,quad=false),
                                  checkpoints=sol.t,
                                  kwargs...)
  isquad(sensealg) && error("Can't get sensitivities of u0 with quadrature.")
  adj_prob = ODEAdjointProblem(sol,g,t,dg,sensealg,checkpoints=checkpoints)
  adj_sol = solve(adj_prob,alg;abstol=abstol,reltol=reltol,kwargs...,save_everystep=false,save_start=false,saveat=eltype(sol[1])[])

  -adj_sol[end][1:length(sol.prob.u0)],
    adj_sol[end][(1:length(sol.prob.p)) .+ length(sol.prob.u0)]'
end

function adjoint_sensitivities(sol,alg,g,t=nothing,dg=nothing;
                               abstol=1e-6,reltol=1e-3,
                               iabstol=abstol, ireltol=reltol,sensealg=SensitivityAlg(checkpointing=!sol.dense),
                               checkpoints=sol.t,
                               kwargs...)
  adj_prob = ODEAdjointProblem(sol,g,t,dg,sensealg,checkpoints=checkpoints)
  isq = isquad(sensealg)
  adj_sol = solve(adj_prob,alg;abstol=abstol,reltol=reltol,
                               save_everystep=isq,save_start=isq,kwargs...)
  !isq && return adj_sol[end][(1:length(sol.prob.p)) .+ length(sol.prob.u0)]'
  integrand = AdjointSensitivityIntegrand(sol,adj_sol,sensealg)

  if t === nothing
    res,err = quadgk(integrand,sol.prob.tspan[1],sol.prob.tspan[2],
                   atol=iabstol,rtol=ireltol)
  else
    res = zero(integrand.p)'
    for i in 1:length(t)-1
      res .+= quadgk(integrand,t[i],t[i+1],
                     atol=iabstol,rtol=ireltol)[1]
    end
    if t[1] != sol.prob.tspan[1]
      res .+= quadgk(integrand,sol.prob.tspan[1],t[1],
                     atol=iabstol,rtol=ireltol)[1]
    end
  end
  res
end


text:
""""""
function robustintervalscale(Fₕ, Fₗ=zeros(size(Fₕ)))
    σₗ = mapslices(iqr, Fₗ, dims=2) |> Array
    σₕ = mapslices(iqr, Fₕ, dims=2) |> Array
    function scale(F)
        σ = mapslices(iqr, F, dims=2) |> Array
        σᵢ = [rampOn(0, 1, σₗ[i], σₕ[i])(σ[i]) for i ∈ 1:length(σ)] # No need to worry about means?
        idxs = vec(isnan.(σᵢ) .| isnan.(σ))
        σᵢ[idxs] .= 0
        σ[idxs] .= Inf
        return F.*σᵢ./σ
    end
    return scale
end
export robustintervalscale


text:
using CSV
using DataFrames
import DimensionalData.formatdims

export RadiationPattern
export gain
export read_HFSS_pattern

@dim ϕ "Azimuth"
@dim θ "Elevation"
@dim f "Frequency"

# Single frequency constructor
function RadiationPattern(data::AbstractArray,azimuth,elevation ;name=NoName(), refdims=(), metadata=NoMetadata())
    RadiationPattern(data, formatdims(data, (ϕ(azimuth),θ(elevation))), refdims, name, metadata)
end

# Multi frequency constructor
function RadiationPattern(data::AbstractArray,azimuth,elevation,frequencies;name=NoName(), refdims=(), metadata=NoMetadata())
    RadiationPattern(data, formatdims(data, (ϕ(azimuth),θ(elevation),f(frequencies))), refdims, name, metadata)
end

# Arbitray dims
function RadiationPattern(data::AbstractArray,dims;name=NoName(), refdims=(), metadata=NoMetadata())
    if length(filter(x->x isa Union{ϕ,θ},dims)) != 2
        error("`dims` must contain at least azimuth and elevation")
    end
    RadiationPattern(data, formatdims(data, dims), refdims, name, metadata)
end

# From kwargs
function RadiationPattern(; data, azimuth,elevation, frequencies,refdims=(), name=NoName(), metadata=NoMetadata())
    RadiationPattern(data, formatdims(data, (ϕ(azimuth),θ(elevation),f(frequencies))), refdims, name, metadata)
end

# Construct from another AbstractDimArray
function RadiationPattern(A::AbstractDimArray; data=data(A), dims=dims(A), refdims=refdims(A), name=name(A), metadata=metadata(A))
    if length(filter(x->x isa Union{ϕ,θ},dims)) != 2
        error("`dims` must contain at least azimuth and elevation")
    end
    RadiationPattern(data, formatdims(data, dims), refdims, name, metadata)
end

"""
    RadiationPattern(f::Function, dim::Dimension [, name])
Apply function `f` across the values of the dimension `dim`
(using `broadcast`), and return the result as a RadiationPattern with
the given dimension. Optionally provide a name for the result.
"""
function RadiationPattern(f::Function, dim::Union{ϕ,θ}, name=Symbol(nameof(f), "(", name(dim), ")"))
    RadiationPattern(f.(val(dim)), (dim,), name)
end

@inline function DimensionalData.rebuild(
    A::RadiationPattern, data::AbstractArray, dims::Tuple, refdims::Tuple, name, metadata
)
    RadiationPattern(data, dims, refdims, name, metadata)
end

dB = Gain{Unitful.LogInfo{:Decibel, 10, 10}, :?, Float64}

function gain(pattern::AbstractArray{T}) where {T <: Gain}
    pattern
end

function gain(pattern::AbstractArray{<:Real})
    @. 20log10(pattern) * u"dB"
end

function gain(pattern::AbstractArray{<:Complex})
    @. 20log10(abs(pattern)) * u"dB"
end

function gain(rp::RadiationPattern)
    gain(rp.data)
end

"""
    read_HFSS_pattern("myAntenna.csv")
Utility function to read the exported fields from HFSS into a `RadiationPattern`.
"""
function read_HFSS_pattern(file)
    # Read Pattern
    pattern = CSV.File(file,header=["ϕ","θ","gain"],skipto=2) |> DataFrame
    # HFSS requires the data to be sorted, thank goodness
    az = pattern[!,:ϕ] |> unique |> sort!
    el = pattern[!,:θ] |> unique |> sort!
    az_range = az[1]:az[2]-az[1]:az[end]
    el_range = el[1]:el[2]-el[1]:el[end]
    gain = reshape(pattern[!,:gain],(length(az_range),length(el_range))) .* u"dB"
    # Create pattern
    RadiationPattern(gain,az_range,el_range;name="Gain (dBi)")
end


text:
using DiffEqBase: CallbackSet

using NQCDynamics: InitialConditions, DynamicsOutputs

abstract type AbstractSelection end

"""
Select the initial conditions from the distribution in order.
"""
struct OrderedSelection{D,I} <: AbstractSelection
    "Distribution that is sampled."
    distribution::D
    indices::I
end

"""
Obtain initial conditions by randomly sampling the distribution.
"""
struct RandomSelection{D} <: AbstractSelection
    "Distribution that is sampled."
    distribution::D
end

function Selection(distribution, selection::AbstractVector)
    @info "Sampling the provided distribution in the range $selection."
    OrderedSelection(distribution, selection)
end

function Selection(distribution, ::Any)
    @info "Sampling randomly from provided distribution."
    RandomSelection(distribution)
end

function (select::OrderedSelection)(prob, i, repeat)
    j = select.indices[i]
    u0 = sample_distribution(prob.p, select.distribution, j)
    DynamicsMethods.create_problem(u0, prob.tspan, prob.p)
end

function (select::RandomSelection)(prob, i, repeat)
    j = rand(1:lastindex(select.distribution))
    u0 = sample_distribution(prob.p, select.distribution, j)
    DynamicsMethods.create_problem(u0, prob.tspan, prob.p)
end

struct SelectWithCallbacks{S<:AbstractSelection,C1,C2,V} <: AbstractSelection
    selection::S
    standard_callbacks::C1
    changing_callbacks::C2
    values::V
    function SelectWithCallbacks(selection, standard_callbacks, output, trajectories; saveat=[])

        callbacks = []
        values = []
        for i=1:trajectories
            cb, vals = DynamicsOutputs.create_saving_callback(output; saveat=saveat)
            push!(callbacks, cb)
            push!(values, vals)
        end

        new{typeof(selection), typeof(standard_callbacks), typeof(callbacks), typeof(values)}(
            selection, standard_callbacks, callbacks, values)
    end
end

function (select::SelectWithCallbacks)(prob, i, repeat)
    prob = select.selection(prob, i, repeat)
    remake(prob, callback=CallbackSet(select.standard_callbacks, select.changing_callbacks[i]))
end


text:
# # Tutorial: missing value handling
#
# This tutorial illustrates the usage of Transducers.jl by stepping
# through various handling of missing values.

using Transducers

# ## Dot product
#
# Here is a simple way to compute a dot product using
# [`mapfoldl`](@ref) and [`MapSplat`](@ref):

mapfoldl(MapSplat(*), +, zip(1:3, 10:2:14))

# Let's see what it does step by step.  First we create a "printer"
# transducer using the following function (see [`Map`](@ref)):

xf_printer(label) = Map() do x
    println(label, ": ", x)
    return x  # just return it as-is
end
nothing  # hide

# This transducer just pass-through the input while printing its value
# (prefixed by a `label`).  Let's sandwich the previous `MapSplat(*)`
# with it:

mapfoldl(
    xf_printer(" input") |> MapSplat(*) |> xf_printer("output"),
    +, zip(1:3, 10:2:14))
flush(stdout)  # hide

# You can see that the input tuple `(1, 10)` is splatted into function
# `*` by `MapSplat` which then outputs `10`.  This is repeated for all
# inputs.

# Perhaps unfortunately, this way of computing a dot product
# propagates any missing values contained in the input arrays to the
# result (which may actually be desired in certain cases).

xs = [1, missing, 3, 2]
ys = [10, 14, missing, 12]
mapfoldl(MapSplat(*), +, zip(xs, ys))

# However, it is very simple to ignore any missing values using
# [`OfType`](@ref):

xf_mdot = OfType(Tuple{Vararg{Number}}) |> MapSplat(*)
mapfoldl(xf_mdot, +, zip(xs, ys))

# Here, `Tuple{Vararg{Number}}` is a type that matches with a tuple of
# any length with numbers.  It does not match with a tuple if it has a
# `missing`.

@assert (1, 0.5) isa Tuple{Vararg{Number}}
@assert (1, 0.5, 2im) isa Tuple{Vararg{Number}}
@assert !((1, missing) isa Tuple{Vararg{Number}})

# ## Covariance
#
# Transducer `xf_mdot` above can also be used to compute the
# covariance.  First, we need the number of pairs of elements in `xs`
# and `ys` that _both_ of them are not `missing`:

nonmissings = mapfoldl(OfType(Tuple{Vararg{Number}}) |> Count(),
                       right,
                       zip(xs, ys);
                       init = 0)
nonmissings  # hide
#-
@assert nonmissings == 2  # hide

# We do this by using [`Count`](@ref) and [`right`](@ref).  `Count`
# ignores input and count the number of times the input is provided.
# Since `OfType(Tuple{Vararg{Number}})` provides the inputs to the
# downstream transducer only if there is no `missing` values, this
# correctly counts the number of non-missing pairs.  Function `right`
# is simply defined as `right(l, r) = r` (and `right(r) = r`).  Thus,
# the whole `mapfoldl` returns the last output of `Count`.  In case
# `Count` never gets called (i.e., there are no non-missing pairs), we
# pass `init=0`.

ans =  # hide
mapfoldl(OfType(Tuple{Vararg{Number}}) |> Count(),
         right,
         zip(Int[], Int[]);
         init = 0)
ans  # hide
#-
@assert ans == 0  # hide

# Finally, we have to pre-process the input to `xf_mdot` by
# subtracting the average.  It's easy to do with `Map`:

using Statistics: mean

function xf_demean(xs, ys)
    xmean = mean(skipmissing(xs))
    ymean = mean(skipmissing(ys))
    return Map(((x, y),) -> (x - xmean, y - ymean))
end

mapfoldl(xf_demean(xs, ys) |> xf_mdot, +, zip(xs, ys)) / nonmissings

# ## Addition
#
# How do we use transducers for vector-to-vector transformation?  Here
# is a function to calculate ``y = x + y`` while ignoring missing
# values in ``x``.  First, mandatory input shape check:

function add_skipmissing!(ys, xs)
    length(ys) == length(xs) || error("length(ys) != length(xs)")
    firstindex(ys) == 1 || error("firstindex(ys) != 1")
#+
# For filtering out missing values from `xs` while tracking indices,
# we use [`Enumerate`](@ref) and [`Filter`](@ref).  To iterate over
# the output of the transducer, [`foreach`](@ref) is used instead of
# [`mapfoldl`](@ref) since mutating an array is better expressed as a
# side-effect than a fold.

    foreach(Enumerate() |> Filter(!(ismissing ∘ last)), xs) do (i, xi)
        @inbounds ys[i] += xi
    end
#+
# We then return the mutated value to behave like the rest of Julia
# functions (`push!`, `mul!`, etc.):

    return ys
end
nothing  # hide

# Example:

ans =  # hide
add_skipmissing!([100, 110, 120], [1, missing, 2])
#-
@assert ans == [101, 110, 122]  # hide

# ## Vectorized reduction
#
# `foldl`, `mapfoldl`, etc. in `Base` support `dims` argument.
# Transducers.jl does not support this argument as of writing.
# However, this can easily be emulated using `eachcol`, `eachrow`, or
# `eachslice` iterators in `Base`.

xs = [
    0       missing 1       2
    3       4       5       missing
    missing 6       7       missing
]

function xf_sum_columns(prototype)
    T = Base.nonmissingtype(eltype(prototype)) # subtract Missing from type
    dims = size(prototype)
    return Scan(add_skipmissing!, CopyInit(zeros(T, dims)))
end
nothing  # hide

# We use [`CopyInit`](@ref) here to allocate the "output array"
# into which the columns are added by `add_skipmissing!`.

if VERSION >= v"1.1-"  # eachcol not in Julia 1.0  #src
ans =  # hide
mapfoldl(xf_sum_columns(xs[:, 1]), right, eachcol(xs))
#-
@assert ans == [3, 12, 13]  # hide

# Above computation returns the sum over each row without taking into
# account the relationship within a column.  Another possibly useful
# reduction is the sum of the columns with no missing values.  This
# can easily be done by prepending a filter:

ans =  # hide
mapfoldl(Filter(x -> !any(ismissing, x)) |> xf_sum_columns(xs[:, 1]),
         right, eachcol(xs))
#md ans  # hide
#src `#md ans` is a hack to avoid Literate.jl to put `continued = true`.
#-
@assert ans == [1, 5, 7]  # hide

# Note that above combination of `Scan` and `right` is redundant.  For
# example, we can simply pass `add_skipmissing!` to "normal" `foldl`:

ans =  # hide
foldl(add_skipmissing!, eachcol(xs), init=zeros(Int, size(xs, 1)))
#-
@assert ans == [3, 12, 13]  # hide

# However, packaging it as a transducer is sometimes useful as it can
# be composed with other transducers and "bottom" reducing function.
# For example, vectorized version of `cumsum` can easily obtained by
# composing it with `append!` (and then `reshape` after `mapfoldl`):

result = mapfoldl(
    xf_sum_columns(xs[:, 1]),
    Completing(append!),
    eachcol(xs);
    init = Int[])
ans =  # hide
reshape(result, (size(xs, 1), :))
#-
@assert ans == [      # hide
    0  0   1   3      # hide
    3  7  12  12      # hide
    0  6  13  13      # hide
]                     # hide
end  # if VERSION >= v"1.1-"  #src

# Note that we need [`Completing`](@ref) here since `append!` does not
# have the unary method used for [`complete`](@ref
# Transducers.complete) protocol.

# ## Argmax
#
# Another useful operation to do ignoring missing values is
# `argmax`/`argmin`.  It can be implemented using `Enumerate() |>
# Filter(!(ismissing ∘ last))` (see also `add_skipmissing!` above)
# composed with [`ScanEmit`](@ref).  We first need to define a
# function to be called by `ScanEmit`:

##                     ,--- current state
##                     |
##                     |              ,-- input
##                     |              |
function argmax_step((argmax, max), (index, value))
    argmax, max = value > max ? (index, value) : (argmax, max)
    return argmax, (argmax, max)
    ##       \        \
    ##        \        \__ next state
    ##         \
    ##          \__ output
end
nothing  # hide

# This function is passed to `ScanEmit` with the initial state:

xf_argmax = Enumerate() |> Filter(!(ismissing ∘ last)) |>
    ScanEmit(argmax_step, (0, typemin(Int)))
##                          |
##                   initial state
nothing  # hide

# As [`ScanEmit`](@ref) is one of the most complex (and powerful)
# transducer, it may require some comments on how above code works:
#
# * The state `(argmax, max)` is initialized to `(0, typemin(Int))` in
#   `xf_argmax`.  This is the first value passed to the first argument
#   `(argmax, max)` of `argmax_step`.
#
# * The upstream transducer `Enumerate()` provides `(index,
#   value)`-pair which becomes the input (the second argument) of
#   `argmax_step`.
#
# * Function `argmax_step` must return a pair.  The first item becomes
#   the output of `ScanEmit`.  In this case that's the index of the
#   largest item seen so far.
#
# * The second item in the returned pair is fed back to `argmax_step`
#   in the next call.
#
# We have the argmax function by extracting the last output of
# `xf_argmax`:

mapfoldl(xf_argmax, right, [1, 3, missing, 2])

# Side note: We use `typemin(Int)` as the initial value of `max` for
# simplicity.  In practice, it should be
# `typemin(eltype(input_array))`.  A more generic solution is to
# special-case the first invocation by using a singleton like
# `nothing`.  Julia can handle small `Union` type such as this (see
# the next section).  Another solution is to use `Init(>)` from
# InitialValues.jl.

# ## Extrema
#
# Transducer `xf_argmax` in the previous section only outputs the
# index of the maximum element so far.  To output the maximum element
# as well, we can simply use [`Scan`](@ref).  Also, while we are at
# it, let's support both argmax and argmin.  To this end, we
# parametrize the function passed to `Scan` by the comparison function
# `>` and `<`.  Another problem with `xf_argmax` is that it does not
# handle non-`Int` input types.  To properly handle different input
# types, we initialize `Scan`'s state with `nothing` and special-case
# the first invocation to return the input as-is.  Following function
# `argext_step` takes the function `>` or `<` and return a function
# appropriate for `Scan.`

argext_step(should_update) =
    (old, (index, value)) ->
        if old === nothing || should_update(old[2], value)
            (index, value)
        else
            old
        end
nothing  # hide

xf_scanext(should_update) = Scan(argext_step(should_update), nothing)
nothing  # hide

# Passing `<` gives us the argmax transducer:

@time begin  #src
ans = # hide
mapfoldl(
    Enumerate() |> OfType(Tuple{Integer, Number}) |> xf_scanext(<),
    right, [1.0, 3.0, missing, 2.0])
end  #src
#-
@assert ans === (2, 3.0) # hide
@show ans  #src

# We now have transducers `xf_scanext(<)` and `xf_scanext(>)` for
# argmax and argmin, respectively.  We can compute them concurrently
# by [`Zip`](@ref)'ing them together:

xf_fullextrema = Enumerate() |> OfType(Tuple{Integer, Number}) |>
    Zip(xf_scanext(>), xf_scanext(<))

@time begin  #src
ans = # hide
mapfoldl(xf_fullextrema, right, [1.0, 3.0, -1.0, missing, 2.0])
end  #src
#-
@assert ans === ((3, -1.0), (2, 3.0))  # hide
@show ans  #src

# This transducer produces a tuple `((argmin, min), (argmax, max))`.
# To output only indices, append an appropriate `Map`:

xf_argextrema =
    xf_fullextrema |> Map() do ((argmin, min), (argmax, max))
        (argmin, argmax)
    end

@time begin  #src
ans = # hide
mapfoldl(xf_argextrema, right, [1.0, 3.0, -1.0, missing, 2.0])
end  #src
#-
@assert ans === (3, 2)  # hide
@show ans  #src


text:
module TestCSetDataStructures
using Test

using Catlab: @present, generator
using Catlab.Theories: compose, id
using Catlab.CSetDataStructures

# Discrete dynamical systems
############################

@present TheoryDDS(FreeSchema) begin
  X::Ob
  Φ::Hom(X,X)
end

const AbstractDDS = AbstractCSetType(TheoryDDS)
const DDS = CSetType(TheoryDDS, index=[:Φ])
@test DDS <: AbstractDDS
@test DDS <: CSet

dds = DDS()
@test keys(tables(dds)) == (:X,)
@test keys(dds.indices) == (:Φ,)
@test nparts(dds, :X) == 0
@test add_part!(dds, :X) == 1
@test nparts(dds, :X) == 1
@test incident(dds, 1, :Φ) == []

set_subpart!(dds, 1, :Φ, 1)
@test subpart(dds, 1, :Φ) == 1
@test incident(dds, 1, :Φ) == [1]

@test add_part!(dds, :X, Φ=1) == 2
@test add_part!(dds, :X, Φ=1) == 3
@test subpart(dds, :Φ) == [1,1,1]
@test subpart(dds, [2,3], :Φ) == [1,1]
@test incident(dds, 1, :Φ) == [1,2,3]

@test has_part(dds, :X)
@test !has_part(dds, :nonpart)
@test has_part(dds, :X, 3)
@test !has_part(dds, :X, 4)
@test has_part(dds, :X, 1:5) == [true, true, true, false, false]

@test has_subpart(dds, :Φ)
@test !has_subpart(dds, :nonsubpart)
@test_throws ArgumentError subpart(dds, 1, :nonsubpart)
@test_throws ArgumentError incident(dds, 1, :nonsuppart)
@test_throws ArgumentError set_subpart!(dds, 1, :nonsubpart, 1)

# Deletion.
dds = DDS()
add_parts!(dds, :X, 3, Φ=[2,3,3])
rem_part!(dds, :X, 2)
@test nparts(dds, :X) == 2
@test subpart(dds, :Φ) == [0,2]
@test incident(dds, 1, :Φ) == []
@test incident(dds, 2, :Φ) == [2]
rem_part!(dds, :X, 2)
@test nparts(dds, :X) == 1
@test subpart(dds, :Φ) == [0]
rem_part!(dds, :X, 1)
@test nparts(dds, :X) == 0

dds = DDS()
add_parts!(dds, :X, 4, Φ=[2,3,3,4])
@test_throws ErrorException rem_parts!(dds, :X, [4,1])
rem_parts!(dds, :X, [1,4])
@test subpart(dds, :Φ) == [1,1]
@test incident(dds, 1, :Φ) == [1,2]
@test incident(dds, 2, :Φ) == []

# Pretty printing.
dds = DDS()
add_parts!(dds, :X, 3, Φ=[2,3,3])
s = sprint(show, dds)
@test startswith(s, "CSet")
@test occursin("X = 1:3", s)
@test occursin("Φ : X → X = ", s)

s = sprint(show, MIME"text/plain"(), dds)
@test startswith(s, "CSet")
@test occursin("X = 1:3", s)
@test occursin("│ X │", s)

s = sprint(show, MIME"text/html"(), dds)
@test startswith(s, "
")
@test occursin("
", s)
@test endswith(rstrip(s), "
")

# Special case of pretty print: empty table.
empty_dds = DDS()
@test !isempty(sprint(show, empty_dds))
@test !isempty(sprint(show, MIME"text/plain"(), empty_dds))
@test !isempty(sprint(show, MIME"text/html"(), empty_dds))

# Error handling when adding parts.
dds = DDS()
add_parts!(dds, :X, 3, Φ=[1,1,1])
@test_throws AssertionError add_part!(dds, :X, Φ=5)
@test nparts(dds, :X) == 3
@test subpart(dds, :Φ) == [1,1,1]
@test_throws AssertionError add_parts!(dds, :X, 2, Φ=[3,6])
@test nparts(dds, :X) == 3
@test incident(dds, 3, :Φ) == []

# Incidence without indexing.
UnindexedDDS = CSetType(TheoryDDS)
dds = UnindexedDDS()
add_parts!(dds, :X, 4, Φ=[3,3,4,4])
@test isempty(keys(dds.indices))
@test incident(dds, 3, :Φ) == [1,2]
@test incident(dds, 4, :Φ) == [3,4]

# Custom table type supported by Tables.jl, here a named tuple of arrays.
dds = DDS(table_type=NamedTuple)
@test tables(dds).X isa NamedTuple{(:Φ,)}
add_parts!(dds, :X, 3, Φ=[1,1,1])
@test nparts(dds, :X) == 3
@test subpart(dds, :Φ) == [1,1,1]
rem_part!(dds, :X, 2)
@test nparts(dds, :X) == 2

dds′ = copy(dds)
@test tables(dds′).X isa NamedTuple{(:Φ,)}
@test dds′ == dds
empty_dds = typeof(dds)()
@test tables(empty_dds).X isa NamedTuple{(:Φ,)}

# Dendrograms
#############

# Strictly speaking, this data structure is not a single dendrogram but a forest
# of dendrograms (whose roots are the elements fixed under the `parent` map) and
# in order to be valid dendrograms, there must be no nontrivial cycles and the
# `height` map must satisfy `compose(parent, height) ≥ height` pointwise.

@present TheoryDendrogram(FreeSchema) begin
  X::Ob
  R::Data
  parent::Hom(X,X)
  height::Attr(X,R)
end

const AbstractDendrogram = AbstractACSetType(TheoryDendrogram)
const Dendrogram = ACSetType(TheoryDendrogram, index=[:parent])
@test Dendrogram <: AbstractDendrogram
@test Dendrogram <: ACSet
@test Dendrogram{Real} <: AbstractDendrogram{Real}
@test_throws ErrorException CSetType(TheoryDendrogram)

d = Dendrogram{Int}()
add_parts!(d, :X, 3, height=0)
add_parts!(d, :X, 2, height=[10,20])
set_subpart!(d, 1:3, :parent, 4)
set_subpart!(d, [4,5], :parent, 5)

@test nparts(d, :X) == 5
@test subpart(d, 1:3, :parent) isa SubArray{Int,1}
@test subpart(d, 1:3, :parent) == [4,4,4]
@test subpart(d, 4, :parent) == 5
@test subpart(d, :, :parent) == [4,4,4,5,5]
@test incident(d, 4, :parent) == [1,2,3]
@test incident(d, 4:5, :parent) isa SubArray{Vector{Int},1}
@test incident(d, 4:5, :parent) == [[1,2,3], [4,5]]
@test has_subpart(d, :height)
@test subpart(d, [1,2,3], :height) == [0,0,0]
@test subpart(d, 4, :height) == 10
@test subpart(d, :, :height) == [0,0,0,10,20]

# Chained accessors.
@test subpart(d, 3, [:parent, :parent]) == 5
@test subpart(d, 3, [:parent, :height]) == 10
@test incident(d, 5, [:parent, :parent]) == [1,2,3,4,5]
@test incident(d, 10, [:parent, :height]) == [1,2,3]

X, parent, height = TheoryDendrogram[[:X, :parent, :height]]
@test subpart(d, 3, parent) == 4
@test subpart(d, 3, compose(parent, height)) == 10
@test subpart(d, 3, id(X)) == 3
@test incident(d, 10, compose(parent, height)) == [1,2,3]
@test subpart(d, parent) == [4,4,4,5,5]
@test subpart(d, id(X)) == 1:5
@test subpart(d, compose(parent, height)) == [10,10,10,20,20]

# Indexing syntax.
@test d[3, :parent] == 4
@test d[3, [:parent, :height]] == 10
@test d[3:5, [:parent, :height]] == [10,20,20]
@test d[:, :parent] == [4,4,4,5,5]
d2 = copy(d)
d2[1, :parent] = 1
d2[2:3, :parent] = 2:3
@test d2[1:3, :parent] == 1:3
d2[1:3, :parent] = 4
@test d2 == d

# Copying parts.
d2 = Dendrogram{Int}()
copy_parts!(d2, d, X=[4,5])
@test nparts(d2, :X) == 2
@test subpart(d2, [1,2], :parent) == [2,2]
@test subpart(d2, [1,2], :height) == [10,20]

du = disjoint_union(d, d2)
@test nparts(du, :X) == 7
@test subpart(du, :parent) == [4,4,4,5,5,7,7]
@test subpart(du, :height) == [0,0,0,10,20,10,20]

# Pretty printing of data attributes.
s = sprint(show, d)
@test startswith(s, "ACSet")
@test occursin("R = Int64", s)
@test occursin("height : X → R = ", s)

s = sprint(show, MIME"text/plain"(), d)
@test startswith(s, "ACSet")

# Allow type inheritance for data attributes.
d_abs = Dendrogram{Number}()
add_parts!(d_abs, :X, 2, height=[10.0, 4])
@test subpart(d_abs, :height) == [10.0, 4]

# Dendrograms with leaves
#------------------------

@present TheoryLDendrogram <: TheoryDendrogram begin
  L::Ob
  leafparent::Hom(L,X)
end

const LDendrogram = ACSetType(TheoryLDendrogram, index=[:parent, :leafparent])

# Copying between C-sets and C′-sets with C != C′.
ld = LDendrogram{Int}()
copy_parts!(ld, d)
@test nparts(ld, :L) == 0
@test subpart(ld, :parent) == subpart(d, :parent)

add_parts!(ld, :L, 3, leafparent=[2,3,4])
@test subpart(ld, :leafparent) == [2,3,4]
d′ = Dendrogram{Int}()
copy_parts!(d′, ld)
@test d′ == d

# Labeled sets
##############

# The simplest example of a C-set with a data attribute, to test data indexing.

@present TheoryLabeledSet(FreeSchema) begin
  X::Ob
  Label::Data
  label::Attr(X,Label)
end

# Labeled sets with index
#------------------------

const IndexedLabeledSet = ACSetType(TheoryLabeledSet, index=[:label])

lset = IndexedLabeledSet{Symbol}()
@test keys(lset.indices) == (:label,)
add_parts!(lset, :X, 2, label=[:foo, :bar])
@test subpart(lset, :, :label) == [:foo, :bar]
@test incident(lset, :foo, :label) == [1]
@test isempty(incident(lset, :nonkey, :label))

add_part!(lset, :X, label=:foo)
@test incident(lset, :foo, :label) == [1,3]
set_subpart!(lset, 1, :label, :baz)
@test subpart(lset, 1, :label) == :baz
@test incident(lset, [:foo,:baz], :label) == [[3],[1]]
set_subpart!(lset, 3, :label, :biz)
@test incident(lset, :foo, :label) == []

# Deletion with indexed data attribute.
lset = IndexedLabeledSet{Symbol}()
add_parts!(lset, :X, 3, label=[:foo, :foo, :bar])
rem_part!(lset, :X, 1)
@test subpart(lset, :label) == [:bar, :foo]
@test incident(lset, [:foo, :bar], :label) == [[2], [1]]

# Deletion with unitialized data attribute.
lset = IndexedLabeledSet{Symbol}()
add_part!(lset, :X)
rem_part!(lset, :X, 1)
@test nparts(lset, :X) == 0

# Pretty-printing with unitialized data attribute.
lset = IndexedLabeledSet{Symbol}()
add_part!(lset, :X)
@test occursin("#undef", sprint(show, lset))
@test occursin("#undef", sprint(show, MIME"text/plain"(), lset))
@test occursin("#undef", sprint(show, MIME"text/html"(), lset))

# Labeled sets with unique index
#-------------------------------

const UniqueIndexedLabeledSet = ACSetType(TheoryLabeledSet,
                                          unique_index=[:label])

lset = UniqueIndexedLabeledSet{Symbol}()
@test keys(lset.indices) == (:label,)
add_parts!(lset, :X, 2, label=[:foo, :bar])
@test subpart(lset, :, :label) == [:foo, :bar]
@test incident(lset, :foo, :label) == 1
@test incident(lset, [:foo,:bar], :label) == [1,2]
@test incident(lset, :nonkey, :label) == 0

set_subpart!(lset, 1, :label, :baz)
@test subpart(lset, 1, :label) == :baz
@test incident(lset, :baz, :label) == 1
@test incident(lset, :foo, :label) == 0

@test_throws ErrorException set_subpart!(lset, 1, :label, :bar)

# @acset macro
##############

@present TheoryDecGraph(FreeSchema) begin
  E::Ob
  V::Ob
  src::Hom(E,V)
  tgt::Hom(E,V)

  X::Data
  dec::Attr(E,X)
end

const DecGraph = ACSetType(TheoryDecGraph, index=[:src,:tgt])

g = @acset DecGraph{String} begin
  V = 4
  E = 4
  src = [1,2,3,4]
  tgt = [2,3,4,1]
  dec = ["a","b","c","d"]
end

@test nparts(g,:V) == 4
@test subpart(g,:,:src) == [1,2,3,4]
@test incident(g,1,:src) == [1]

# Test mapping
#-------------

f(s::String) = Int(s[1])

h1 = map(g, dec = f)
h2 = map(g, X = f)

@test h1 == h1

@test subpart(h1,:src) == subpart(g,:src)
@test typeof(h1).parameters[3] == Tuple{Int}
@test subpart(h1,:dec) == f.(["a","b","c","d"])

@present TheoryLabelledDecGraph <: TheoryDecGraph begin
  label::Attr(V,X)
end

const LabelledDecGraph = ACSetType(TheoryLabelledDecGraph, index=[:src,:tgt])

g = @acset LabelledDecGraph{String} begin
  V = 4
  E = 4

  src = [1,2,3,4]
  tgt = [2,3,4,1]

  dec = ["a","b","c","d"]
  label = ["w", "x", "y", "z"]
end

h1 = map(g, X = f)
@test subpart(h1,:label) == f.(["w","x","y","z"])

h2 = map(g, dec = f, label = i -> 3)
@test subpart(h2,:dec) == f.(["a","b","c","d"])
@test subpart(h2,:label) == [3,3,3,3]

@test_throws Any map(g, dec = f)

end


text:
function handle_step(obj::ODEFixedStepHandler, arg0::ODEStateAndDerivative, arg1::jboolean)
    return jcall(obj, "handleStep", void, (ODEStateAndDerivative, jboolean), arg0, arg1)
end

function init(obj::ODEFixedStepHandler, arg0::ODEStateAndDerivative, arg1::jdouble)
    return jcall(obj, "init", void, (ODEStateAndDerivative, jdouble), arg0, arg1)
end


text:
export DSSAIntegrator

mutable struct DSSAIntegrator{F,uType,tType,P,S,CB,SA,OPT,TS,chanType, chanS} <: DiffEqBase.DEIntegrator{SSAStepper,Nothing,uType,tType}
    f::F
    u::uType
    t::tType
    tprev::tType
    p::P
    sol::S
    i::Int
    tstop::tType
    cb::CB
    saveat::SA
    save_everystep::Bool
    save_end::Bool
    cur_saveat::Int
    opts::OPT
    tstops::TS
    tstops_idx::Int
    u_modified::Bool
    keep_stepping::Bool  # false if should terminate a simulation
    de_chan::chanType
    chan_sol::chanS
    delayjumpsets::DelayJumpSet
    save_delay_channel::Bool
end

mutable struct DSSASolution{uType,uType2,DType,tType,rateType,P,A,IType,DE,chanS}
    u::uType
    u_analytic::uType2
    errors::DType
    t::tType
    k::rateType
    prob::P
    alg::A
    interp::IType
    dense::Bool
    tslocation::Int
    destats::DE
    retcode::Symbol
    channel::chanS
end
function DSSASolution(ode_sol::DiffEqBase.ODESolution, channel)
    @unpack u, u_analytic, errors, t, k, prob, alg, interp, dense, tslocation, destats, retcode = ode_sol
    DSSASolution{typeof(u), typeof(u_analytic), typeof(errors), typeof(t), typeof(k), typeof(prob), typeof(alg), typeof(interp), typeof(destats), typeof(channel)}(u, u_analytic, errors, t, k, prob, alg, interp, dense, tslocation, destats, retcode, channel)
end

function Base.show(io::IO, m::MIME"text/plain", A::DSSASolution)
    println(io,string("retcode: ",A.retcode))
    println(io,string("Interpolation: "),DiffEqBase.interp_summary(A.interp))
    print(io,"t: ")
    show(io,m,A.t)
    println(io)
    print(io,"u: ")
    show(io,m,A.u)
    println(io)
    print(io,"channel: ")
    show(io,m,A.channel)
    println(io,"\n===\nUse sol.u to check the state variable and sol.channel to check the delay channel solution.\n===")
end

Base.@propagate_inbounds Base.getindex(A::DSSASolution, i::Int) = [A.u[i], A.channel[i]]
Base.@propagate_inbounds Base.getindex(A::DSSASolution, i::Int, ::Colon) = [A.u[j][i] for j in 1:length(A.t)]

function (A::DSSASolution)(s::Symbol,i::Int)
    if s == :channel
        @assert i <= length(A.channel[1])
        return [A.channel[j][i] for j in 1:length(A.t)]
    elseif s ==:u
        @assert i <= length(A.u[1])
        return [A.u[j][i] for j in 1:length(A.t)]
    end
end

(integrator::DSSAIntegrator)(t) = copy(integrator.u)
(integrator::DSSAIntegrator)(out,t) = (out .= integrator.u)

function DiffEqBase.u_modified!(integrator::DSSAIntegrator,bool::Bool)
    integrator.u_modified = bool
end

function DiffEqBase.__solve(djump_prob::DelayJumpProblem,
                    alg::SSAStepper;
                            kwargs...)
    integrator = init(djump_prob,alg;kwargs...)
    solve!(integrator)
    if integrator.save_delay_channel
        DSSASolution(integrator.sol, integrator.chan_sol)
    else
        integrator.sol
    end
end

## Initiate the Problem

function DiffEqBase.__init(djump_prob::DelayJumpProblem,
                         alg::SSAStepper;
                         save_start = true,
                         save_end = true,
                         seed = nothing,
                         alias_jump = Threads.threadid() == 1,
                         saveat = nothing,
                         callback = nothing,
                         tstops = eltype(djump_prob.prob.tspan)[],
                         numsteps_hint=100,
                         save_delay_channel=false)
    if !(djump_prob.prob isa DiscreteProblem)
        error("SSAStepper only supports DiscreteProblems.")
    end
    @assert isempty(djump_prob.jump_callback.continuous_callbacks) # check if continuous_callbacks is empty
    ## Set up the cb which encodes the aggregator information
    if alias_jump
      cb = djump_prob.jump_callback.discrete_callbacks[end]
      if seed !== nothing
          Random.seed!(cb.condition.rng,seed)
      end
    else
      cb = deepcopy(djump_prob.jump_callback.discrete_callbacks[end])
      if seed === nothing
          Random.seed!(cb.condition.rng,seed_multiplier()*rand(UInt64))
      else
          Random.seed!(cb.condition.rng,seed)
      end
    end

    opts = (callback = CallbackSet(callback),)
    prob = djump_prob.prob
    de_chan0 = convert(Vector{Vector{typeof(prob.tspan[1])}},djump_prob.de_chan0)
    if save_start
        t = [prob.tspan[1]]
        u = [copy(prob.u0)]
        chan_sol = [deepcopy(de_chan0)] # DelaySSA: build chan_sol
    else
        t = typeof(prob.tspan[1])[]
        u = typeof(prob.u0)[]
        chan_sol = typeof(de_chan0)[] # DelaySSA: build chan_sol
    end

    sol = DiffEqBase.build_solution(prob,alg,t,u, dense=false,
    calculate_error = false,  destats = DiffEqBase.DEStats(0),
                         interp = DiffEqBase.ConstantInterpolation(t,u))
    save_everystep = any(cb.save_positions)

    if typeof(saveat) <: Number
        _saveat = prob.tspan[1]:saveat:prob.tspan[2]
    else
        _saveat = saveat
    end

   if _saveat !== nothing && !isempty(_saveat) && _saveat[1] == prob.tspan[1]
       cur_saveat = 2
   else
       cur_saveat = 1
   end

   if _saveat !== nothing && !isempty(_saveat)
     sizehint!(u,length(_saveat)+1)
     sizehint!(t,length(_saveat)+1)
   elseif save_everystep
     sizehint!(u,numsteps_hint)
     sizehint!(t,numsteps_hint)
   else
     sizehint!(u,save_start+save_end)
     sizehint!(t,save_start+save_end)
   end

    integrator = DSSAIntegrator(prob.f,copy(prob.u0),prob.tspan[1],prob.tspan[1],prob.p, sol,1,prob.tspan[1],cb,_saveat,save_everystep,save_end,cur_saveat,opts,tstops,1,false,true,deepcopy(de_chan0),chan_sol,djump_prob.delayjumpsets, save_delay_channel)

    cb.initialize(cb,integrator.u,prob.tspan[1],integrator)
    DiffEqBase.initialize!(opts.callback,integrator.u,prob.tspan[1],integrator)
    integrator
end

function DiffEqBase.add_tstop!(integrator::DSSAIntegrator,tstop)
    if tstop > integrator.t
        future_tstops = @view integrator.tstops[integrator.tstops_idx:end]
        insert_index = integrator.tstops_idx + searchsortedfirst(future_tstops, tstop) - 1
        Base.insert!(integrator.tstops, insert_index, tstop)
    end
end

function DiffEqBase.solve!(integrator::DSSAIntegrator)

    end_time = integrator.sol.prob.tspan[2]
    while should_continue_solve(integrator) # It stops before adding a tstop over
        step!(integrator)
    end
    last_t = copy(integrator.t)
    integrator.t = end_time

    if typeof(integrator.cb.affect!) <: DelayDirectJumpAggregation
        saveat_function_direct_method!(integrator, last_t)
    else
        saveat_function!(integrator, last_t)
    end

    if integrator.save_end && integrator.sol.t[end] != end_time
        saveat_end_function!(integrator, last_t)
    end
    DiffEqBase.finalize!(integrator.opts.callback, integrator.u, integrator.t, integrator)
end

function saveat_end_function!(integrator, prev_t)
    # save last t
    end_time = integrator.sol.prob.tspan[2]
    push!(integrator.sol.t,end_time)

    # save last u
    t_final_gap = end_time - prev_t
    if typeof(integrator.cb.affect!) <: DelayDirectJumpAggregation
        T1, T2 = create_Tstruct(integrator.de_chan)
        update_state_final_jump!(integrator.cb.affect!, integrator, t_final_gap, T1, T2)
    end
    push!(integrator.sol.u,copy(integrator.u))

    # save last de_chan
    if integrator.save_delay_channel
        last_chan = deepcopy(integrator.de_chan)
        shift_delay_channel!(last_chan, t_final_gap)
        update_delay_channel!(last_chan)
        push!(integrator.chan_sol,last_chan)
    end
end

"""
    function saveat_function!(integrator, prev_t)

Note that this function does not change u and de_chan, but only takes the snapshots of u and de_chan accordingly
"""
function saveat_function!(integrator, prev_t)
    if integrator.saveat !== nothing && !isempty(integrator.saveat)
        # Split to help prediction
        last_saved_t = prev_t
        prev_de_chan = deepcopy(integrator.de_chan)
        while integrator.cur_saveat < length(integrator.saveat) && (integrator.saveat[integrator.cur_saveat] < integrator.t)
            tgap = integrator.saveat[integrator.cur_saveat] - last_saved_t
            push!(integrator.sol.t,integrator.saveat[integrator.cur_saveat])
            push!(integrator.sol.u,copy(integrator.u))
            if integrator.save_delay_channel
                shift_delay_channel!(prev_de_chan, tgap)
                push!(integrator.chan_sol,deepcopy(prev_de_chan))
            end
            last_saved_t = integrator.saveat[integrator.cur_saveat]
            integrator.cur_saveat += 1
        end
    end
end

"""
    function saveat_function_direct_method!(integrator, prev_t)

Note that this function does not change u and de_chan, but only takes the snapshots of u and de_chan accordingly
"""
function saveat_function_direct_method!(integrator, prev_t)
    # Special to Direct method
    if integrator.saveat !== nothing && !isempty(integrator.saveat)
        T1, T2 = create_Tstruct(integrator.de_chan)
        # Split to help prediction
        last_saved_t = prev_t
        prev_de_chan = deepcopy(integrator.de_chan)
        while integrator.cur_saveat < length(integrator.saveat) && (integrator.saveat[integrator.cur_saveat] < integrator.t)

            tgap = integrator.saveat[integrator.cur_saveat] - last_saved_t
            push!(integrator.sol.t,integrator.saveat[integrator.cur_saveat])

            # Special to Direct method
            update_state_final_jump!(integrator.cb.affect!, integrator, tgap, T1, T2)
            push!(integrator.sol.u,copy(integrator.u))
            if integrator.save_delay_channel
                shift_delay_channel!(prev_de_chan, tgap)
                # Special for Direct method
                update_delay_channel!(prev_de_chan) #
                push!(integrator.chan_sol,deepcopy(prev_de_chan))
            end
            last_saved_t = integrator.saveat[integrator.cur_saveat]
            integrator.cur_saveat += 1
        end
    end
end

# The Jump aggregators should not register the next jump through add_tstop! for SSAIntegrator
# such that we can achieve maximum performance
@inline function register_next_jump_time!(integrator::DSSAIntegrator, p::AbstractDSSAJumpAggregator, t)
    integrator.tstop = p.next_jump_time
    nothing
end

## Delay SSA modify
function DiffEqBase.step!(integrator::DSSAIntegrator)
    integrator.tprev = integrator.t
    next_jump_time = integrator.tstop > integrator.t ? integrator.tstop : typemax(integrator.tstop)

    doaffect = false
    if !isempty(integrator.tstops) &&
        integrator.tstops_idx <= length(integrator.tstops) &&
        integrator.tstops[integrator.tstops_idx] < next_jump_time

        integrator.t = integrator.tstops[integrator.tstops_idx]
        integrator.tstops_idx += 1
    else
        integrator.t = integrator.tstop
        doaffect = true # delay effect until after saveat
    end

    if typeof(integrator.cb.affect!) <: DelayDirectJumpAggregation
        saveat_function_direct_method!(integrator, copy(integrator.tprev))
    else
        saveat_function!(integrator, copy(integrator.tprev))
    end

    # FP error means the new time may equal the old if the next jump time is
    # sufficiently small, hence we add this check to execute jumps until
    # this is no longer true.
    while integrator.t == integrator.tstop
        doaffect && integrator.cb.affect!(integrator)
    end

    if !(typeof(integrator.opts.callback.discrete_callbacks)<:Tuple{})
        discrete_modified, saved_in_cb = DiffEqBase.apply_discrete_callback!(integrator,integrator.opts.callback.discrete_callbacks...)
    else
        saved_in_cb = false
    end

    !saved_in_cb && savevalues!(integrator)

    nothing
end

function DiffEqBase.savevalues!(integrator::DSSAIntegrator,force=false)
    saved, savedexactly = false, false

    # No saveat in here since it would only use previous values,
    # so in the specific case of SSAStepper it's already handled

    if integrator.save_everystep || force
        saved = true
        savedexactly = true
        push!(integrator.sol.t,integrator.t)
        push!(integrator.sol.u,copy(integrator.u))
        if integrator.save_delay_channel
            push!(integrator.chan_sol,deepcopy(integrator.de_chan))
        end
    end

    saved, savedexactly
end

function should_continue_solve(integrator::DSSAIntegrator)
    end_time = integrator.sol.prob.tspan[2]

    # we continue the solve if there is a tstop between now and end_time
    has_tstop = !isempty(integrator.tstops) &&
        integrator.tstops_idx <= length(integrator.tstops) &&
        integrator.tstops[integrator.tstops_idx] < end_time

    # we continue the solve if there will be a jump between now and end_time
    has_jump = integrator.t < integrator.tstop < end_time

    integrator.keep_stepping && (has_jump || has_tstop)
end

num_constant_rate_jumps(aggregator::AbstractDSSAJumpAggregator) = length(aggregator.rates)


text:
@testset "data formatting" begin
    @test true
end


text:
module GridPortionsTests

using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry

domain = (-1,1,-1,1)
partition = (10,10)
oldgrid = CartesianGrid(domain,partition)

const R = 0.7

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

oldcell_to_coods = get_cell_coordinates(oldgrid)

oldcell_to_mask = lazy_map(is_in,oldcell_to_coods)

grid = GridPortion(oldgrid,oldcell_to_mask)
test_grid(grid)

topo = GridTopology(grid)
test_grid_topology(topo)

labels = FaceLabeling(topo)

model = DiscreteModel(grid,topo,labels)
test_discrete_model(model)

#using Gridap.Visualization

#writevtk(model,"model")

#
#write_vtk_file(
#  grid,
#  "grid",
#  celldata=["oldcell"=>grid.cell_to_oldcell],
#  nodaldata=["oldnode"=>grid.node_to_oldnode])

end # module


text:
function pbtf(
    y::Vector{T},                 # response
    x::Vector{T},                 # predictor
    k::Int;                       # order of desired polynomial, only 1 or 2 is acceptable
    thinning::Bool=false,         # whether to apply the thinning technique
    nbins::Int=100,               # number of bins if thinning is applied
    λ::Union{T,Nothing}=nothing,  # the Moreau envelope parameter
    s::T=1e-3,                    # IG(s,r) prior for σ²
    r::T=1e-3,
    s₂::Union{T,Nothing}=nothing, # β'(n-k,s₂) prior for α
    nsample::Int=3000,            # number of MCMC sample iterations
    nthin::Int=5,                 # thinning the chain after MCMC sampling
    verbose::Bool=true            # whether to display progress bar
    ) where T<: Real
    problem = pbtfProblem(y,x,k,thinning=thinning,nbins=nbins,λ=λ,s=s,r=r,s₂=s₂)
    X = [ones(T,problem.n) problem.xgrid]
    q₀ = [problem.TM*X*(X\problem.ybar);0;0]
    if verbose
        elapse_time = @elapsed results = mcmc_with_warmup(GLOBAL_RNG, problem, nsample, initialization=(ϵ=0.01,q=q₀),
        warmup_stages=default_warmup_stages(;stepsize_search=nothing),reporter=ProgressMeterReport())
        println(summarize_tree_statistics(results.tree_statistics))
    else
        elapse_time = @elapsed results = mcmc_with_warmup(GLOBAL_RNG, problem, nsample, initialization=(ϵ=0.01,q=q₀),
        warmup_stages=default_warmup_stages(;stepsize_search=nothing),reporter=NoProgressReport())
    end
    results_matrix = reshape([(results.chain...)...], length(results.chain[1]), length(results.chain))
    results_matrix[1:problem.n,:] = problem.TM \ results_matrix[1:problem.n,:]
    results_matrix[end-1:end,:] = exp.(results_matrix[end-1:end,:])
    chn = Chains(transpose(results_matrix),["β[" .* string.(1:problem.n) .*"]"; "σ²";"α"],thin=nthin)
    result_summary = DataFrame(summarize(chn))
    result_quantile = DataFrame(quantile(chn))
    return result_summary,result_quantile,problem,elapse_time
end

function pbsrtf(
    y::Vector{T},                # response
    x::Vector{T},                # predictor
    k::Int,                      # order of desired polynomial
    restriction::String; # shape restriction
    lb::Union{Vector{T},Nothing}=nothing, #lower bound for β
    ub::Union{Vector{T},Nothing}=nothing, #upper bound for β
    scale::Bool=true,            # whether to scale x and y to the range of 0-10
    thinning::Bool=false,        # whether to apply the thinning technique
    nbins::Int=100,              # number of bins if thinning is applied
    λ::Union{T,Nothing}=nothing, # Moreau envelope parameter
    s::T=1e-3,                   # IG(s,r) prior for σ²
    r::T=1e-3,
    μ::Union{T,Nothing}=nothing, # exp(μ) prior for α
    nsample::Int=1000,           # number of MCMC sample iterations
    nthin::Int=1,                # thinning the chain after MCMC sampling
    verbose::Bool=true           # whether to display progress bar
    ) where T<: Real
    xmin = ymin = 0.0
    xscale = yscale = 1.0
    if scale
        xmin = minimum(x)
        ymin = minimum(y)
        xscale = (maximum(x)-minimum(x))/10.0
        yscale = (maximum(y)-minimum(y))/10.0
        y = (y .- ymin)./yscale
        x = (x .- xmin)./xscale
    end
    problem = pbsrtfProblem(y,x,k,restriction,lb=lb,ub=ub,λ=λ,s=s,r=r,μ=μ,thinning=thinning,nbins=nbins)
    if verbose
        elapse_time = @elapsed results = mcmc_with_warmup(GLOBAL_RNG, problem, nsample, initialization=(ϵ=0.01,),
        warmup_stages=default_warmup_stages(;stepsize_search=nothing,
        init_steps=50,middle_steps=50,doubling_stages=3,terminating_steps=50),reporter=ProgressMeterReport())
        println(summarize_tree_statistics(results.tree_statistics))
    else
        elapse_time = @elapsed results = mcmc_with_warmup(GLOBAL_RNG, problem, nsample, initialization=(ϵ=0.01,),
        warmup_stages=default_warmup_stages(;stepsize_search=nothing,
        init_steps=50,middle_steps=50,doubling_stages=3,terminating_steps=50),reporter=NoProgressReport())
    end
    results_matrix = reshape([(results.chain...)...], length(results.chain[1]), length(results.chain))
    results_matrix[end-1:end,:] = exp.(results_matrix[end-1:end,:])
    problem.x = problem.x .*xscale .+ xmin
    problem.y = problem.y .*yscale .+ ymin
    problem.xgrid = problem.xgrid .*xscale .+ xmin
    problem.ybar = problem.ybar .*yscale .+ ymin
    results_matrix[1:end-2,:] .*= yscale
    results_matrix[1:end-2,:] .+= ymin
    results_matrix[end-1,:] .*= yscale^2
    chn = Chains(transpose(results_matrix),["β[" .* string.(1:problem.n) .*"]"; "σ²";"α"],thin=nthin)
    result_summary = DataFrame(summarize(chn))
    result_quantile = DataFrame(quantile(chn))
    return result_summary,result_quantile,problem,elapse_time
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	load_co2_cap(setup::Dict, path::AbstractString, inputs_co2::Dict)

Function for reading input parameters related to CO$_2$ emissions cap constraints
"""
function load_co2_cap(setup::Dict, path::AbstractString, inputs_co2::Dict)
    # Definition of Cap requirements by zone (as Max Mtons)
    inputs_co2["dfCO2Cap"] = DataFrame(CSV.File(joinpath(path, "CO2_cap.csv"), header = true), copycols = true)
    cap = count(s -> startswith(String(s), "CO_2_Cap_Zone"), names(inputs_co2["dfCO2Cap"]))
    first_col = findall(s -> s == "CO_2_Cap_Zone_1", names(inputs_co2["dfCO2Cap"]))[1]
    last_col = findall(s -> s == "CO_2_Cap_Zone_$cap", names(inputs_co2["dfCO2Cap"]))[1]

    inputs_co2["dfCO2CapZones"] = Matrix{Float64}(inputs_co2["dfCO2Cap"][:, first_col:last_col])
    inputs_co2["NCO2Cap"] = cap

    # Emission limits
    if setup["CO2Cap"] == 1
        #  CO2 emissions cap in mass
        first_col = findall(s -> s == "CO_2_Max_Mtons_1", names(inputs_co2["dfCO2Cap"]))[1]
        last_col = findall(s -> s == "CO_2_Max_Mtons_$cap", names(inputs_co2["dfCO2Cap"]))[1]
        # note the default inputs is in million tons; when scaled, the constraint unit is kton, when not scaled, the constraint unit is ton
        if setup["ParameterScale"] == 1
            inputs_co2["dfMaxCO2"] = Matrix{Float64}(inputs_co2["dfCO2Cap"][:, first_col:last_col]) * (1e6) / ModelScalingFactor
        else
            inputs_co2["dfMaxCO2"] = Matrix{Float64}(inputs_co2["dfCO2Cap"][:, first_col:last_col]) * (1e6)
        end

    elseif (setup["CO2Cap"] == 2 || setup["CO2Cap"] == 3)
        #  CO2 emissions rate applied per MWh
        first_col = findall(s -> s == "CO_2_Max_tons_MWh_1", names(inputs_co2["dfCO2Cap"]))[1]
        last_col = findall(s -> s == "CO_2_Max_tons_MWh_$cap", names(inputs_co2["dfCO2Cap"]))[1]
		# when scaled, the constraint unit is kton, thus the emission rate should be in kton/GWh = ton/MWh, when not scaled, the constraint unit is ton/MWh
        if setup["ParameterScale"] == 1
            inputs_co2["dfMaxCO2Rate"] = Matrix{Float64}(inputs_co2["dfCO2Cap"][:, first_col:last_col])
        else
            inputs_co2["dfMaxCO2Rate"] = Matrix{Float64}(inputs_co2["dfCO2Cap"][:, first_col:last_col])
        end
    end

    println("CO2_cap.csv Successfully Read!")

    return inputs_co2
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
N = 10
#################################
# ConstrainedIndexValueResampling
#################################

# Truncate each of the indices for x at 0.8 their standard deviation around the mean
constraints_x_inds = TruncateStd(0.8)

# Truncate each of the indices for y at 1.5 their standard deviation around the mean
constraints_y_inds = TruncateStd(1.5)

# Truncate each of the values of x at the 20th percentile range
constraints_x_vals = [TruncateQuantiles(0.4, 0.6) for i = 1:N];

# Truncate each of the values of x at the 80th percentile range
constraints_y_vals = [TruncateQuantiles(0.1, 0.9) for i = 1:N];

cs_x = (constraints_x_inds, constraints_x_vals)
cs_y = (constraints_y_inds, constraints_y_vals)

resampling_idxval_scheme = ConstrainedIndexValueResampling(cs_x, cs_y, cs_x)

@test length(resampling_idxval_scheme) == 3

# There should be two
@test length(resampling_idxval_scheme[1]) == 2
@test length(resampling_idxval_scheme[2]) == 2
@test resampling_idxval_scheme[1][1] isa SamplingConstraint
@test resampling_idxval_scheme[1][2] isa Vector{<:SamplingConstraint}
@test length(resampling_idxval_scheme[1][2]) == N
@test resampling_idxval_scheme[2][1] isa SamplingConstraint
@test resampling_idxval_scheme[2][2] isa Vector{<:SamplingConstraint}
@test length(resampling_idxval_scheme[2][2]) == N
@test resampling_idxval_scheme.n == 1

resampling_idxval_scheme = ConstrainedIndexValueResampling(105, cs_x, cs_y)
@test resampling_idxval_scheme.n == 105


text:
###################################################
## visualize.jl
## visualize network (made to be extended)
###################################################

"""
    `NetworkVisualizer` : abstract type that represents visualization in a network
    has to implement ands initialize:
    - attribute `network::Network`
    has to implement, automatically initialized
    - attribute `window::RenderWindow`
    - attribute `view::View`
    - attribute `nodes::Vector{CircleShape}`
    - attribute `roads::Dict{Tuple{Int,Int},Line}`
    - attribute `nodeRadius::Float64` (to scale things)
    - attribute `colors::VizColors` to get all the colors of the visualization

    can implement
    - method `visualInit` => initialize things
    - method `visualEvent`  => receive one event (may be 0 or several per frame)
    - method `visualUpdate` => called each frame, to update and draw objects
    - method `visualRedraw` => called when redrawing
"""
abstract type NetworkVisualizer end

"""
    `visualInit` initialize things
"""
function visualInit(v::NetworkVisualizer)
end

"""
    `visualEvent` => called each frame, is given the events
"""
function visualEvent(v::NetworkVisualizer, event::Event)
end

"""
    `visualStartUpdate` => called each frame, to update and draw objects, is given frame-time
    (before other objects drawn)
"""
function visualStartUpdate(v::NetworkVisualizer,frameTime::Float64)
end

"""
    `visualEndUpdate` => called each frame, to update and draw objects, is given frame-time
     (after other objects drawn)
"""
function visualEndUpdate(v::NetworkVisualizer,frameTime::Float64)
end

"""
    `visualRedraw` => called when change in nodeRadius
"""
function visualRedraw(v::NetworkVisualizer)
end

function Base.show(io::IO, viz::NetworkVisualizer)
    typeName = split(string(typeof(viz)),".")[end]
    println(io,"Network Visualizer: $(typeName)")
end

"""
    `copyVisualData` : copy `ref` visual data into `v`
"""
function copyVisualData(ref::NetworkVisualizer,v::NetworkVisualizer)
    v.window = ref.window
    v.nodes  = ref.nodes
    v.roads  = ref.roads
    v.colors = ref.colors
end

"""
    `visualize`: visualize a network
    - specifics of the visualization are given by NetworkVisualizer object
    - if given network, calls `NodeInfo` visualizer is chosen
"""
function visualize(v::NetworkVisualizer)
    v.nodeRadius = 10.

    #create nodes
    v.nodes = CircleShape[CircleShape() for i in 1:length(v.network.nodes)]

    #create roads
    v.roads = Dict{Tuple{Int,Int},Line}()
    for ((o,d),r) in v.network.roads
        road = Line(Vector2f(0.,0.),Vector2f(1000.,0.))
        v.roads[o,d] = road
    end

    # Defines the window, an event listener, and view
    window_w, window_h = 1200,1200
    v.window = RenderWindow("Network Visualization", window_w, window_h)
    set_vsync_enabled(v.window, true)
    event = Event()

    # Set up the initial view
    minX, maxX, minY, maxY = boundingBox(Tuple{Float64,Float64}[(n.x,n.y) for n in v.network.nodes])
    # Do the Y-axis transformation
    minY, maxY = -maxY, -minY
    networkLength = max(maxX-minX, maxY-minY)
    viewWidth = max(maxX-minX, (maxY-minY)*window_w/window_h)
    viewHeigth = max(maxY-minY, (maxX-minX)*window_h/window_w)
    v.view = View(Vector2f((minX+maxX)/2,(minY+maxY)/2), Vector2f(viewWidth, viewHeigth))
    zoomLevel = 1.0
    hideNodes = true
    # init visualizer
    visualInit(v)
    redraw!(v)
    clock = Clock()
    # gc_enable(false)
    while isopen(v.window)
        frameTime = Float64(as_seconds(restart(clock)))
        while pollevent(v.window, event)
            if get_type(event) == EventType.CLOSED
                close(v.window)
            end
            if get_type(event) == EventType.RESIZED
                size = get_size(event)
                window_w, window_h = size.width, size.height
                viewWidth = max(maxX-minX, (maxY-minY)*window_w/window_h)
                viewHeigth = max(maxY-minY, (maxX-minX)*window_h/window_w)
                set_size(v.view, Vector2f(viewWidth, viewHeigth))
                zoom(v.view, zoomLevel)
            end
            if get_type(event) == EventType.KEY_PRESSED
                k = get_key(event).key_code
                if k == KeyCode.ESCAPE || k == KeyCode.Q
                    close(v.window)
                elseif k == KeyCode.A
                    v.nodeRadius *= 1.3
                    redraw!(v)
                elseif k == KeyCode.S
                    v.nodeRadius /= 1.3
                    redraw!(v)
                elseif k == KeyCode.D
                    hideNodes = !hideNodes
                end
            end
            # additional event processing
            visualEvent(v,event)
        end
		if is_key_pressed(KeyCode.LEFT)
			move(v.view, Vector2f(-networkLength/2*frameTime*zoomLevel,0.))
		end
        if is_key_pressed(KeyCode.RIGHT)
			move(v.view, Vector2f(networkLength/2*frameTime*zoomLevel,0.))
		end
        if is_key_pressed(KeyCode.UP)
			move(v.view, Vector2f(0.,-networkLength/2*frameTime*zoomLevel))
		end
        if is_key_pressed(KeyCode.DOWN)
			move(v.view, Vector2f(0.,networkLength/2*frameTime*zoomLevel))
		end
        if is_key_pressed(KeyCode.Z)
            zoom(v.view, 0.6^frameTime)
            zoomLevel = get_size(v.view).x/viewWidth
		end
		if is_key_pressed(KeyCode.X)
			zoom(v.view, 1/(0.6^frameTime))
            zoomLevel = get_size(v.view).x/viewWidth
		end
        set_view(v.window,v.view)
        clear(v.window, SFML.Color(210,210,210))
        # additional updates
        visualStartUpdate(v, frameTime)
        for road in values(v.roads)
            draw(v.window,road)
        end
        if !hideNodes
            for node in v.nodes
                draw(v.window,node)
            end
        end

        # additional updates
        visualEndUpdate(v, frameTime)

        display(v.window)
    end
    # gc_enable()
end

visualize(n::Network) = visualize(NetworkViz(n))
visualize(r::RoutingPaths) = visualize(RoutingViz(r))

"""
    `redraw!` is a helper function that updates all the coordinates
"""
function redraw!(v::NetworkVisualizer)
    n = v.network
    # create position vectors (inverse y axe for plotting)
    positions = Vector2f[Vector2f(node.x,-node.y) for node in n.nodes]
    for i in 1:length(positions)
       positions[i] = Vector2f(n.nodes[i].x,-n.nodes[i].y)
    end

    #positions nodes
    for (i, no) in enumerate(v.nodes)
        set_radius(no, v.nodeRadius)
        set_position(no, positions[i] - Vector2f(v.nodeRadius,v.nodeRadius))
        set_fillcolor(no, nodeColor(v.colors, v.network.nodes[i]))
    end

    #positions roads
    #SFML bug with lines! have to recreate them...
    for ((o,d),r) in v.network.roads
        x = positions[o].y - positions[d].y
        y = positions[d].x - positions[o].x
        l = sqrt(x^2+y^2); x = x/l; y = y/l
        offset = Vector2f(x*v.nodeRadius*3./5.,y*v.nodeRadius*3./5.)
        road = Line(positions[o]+offset,positions[d]+offset, 4*v.nodeRadius/5.)
        set_fillcolor(road,roadColor(v.colors, r))
        v.roads[o,d] = road
    end
    visualRedraw(v)
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
using Mimi

@defcomp RegionAggregatorSum begin

    inputregions = Index()
    outputregions = Index()

    input_output_mapping = Parameter{String}(index=[inputregions]) # one element per input region containing it's corresponding output region
    input_output_mapping_int = Variable{Int}(index=[inputregions]) # internally computed for speed up

    input_region_names = Parameter{Vector{String}}(index=inputregions)
    output_region_names = Parameter{Vector{String}}(index=outputregions)

    input = Parameter(index=[time, inputregions])
    output = Variable(index=[time, outputregions])

    function init(p,v,d)
        idxs = indexin(p.input_output_mapping, p.output_region_names)
        !isnothing(findfirst(i -> isnothing(i), idxs)) ? error("All provided region names in the RegionAggregatorSum's input_output_mapping Parameter must exist in the output_region_names Parameter.") : nothing
        v.input_output_mapping_int[:] = idxs
    end

    function run_timestep(p,v,d,t)
        v.output[t, :] .= 0.

        for i in d.inputregions
            v.output[t, v.input_output_mapping_int[i]] += p.input[t,i]
        end
    end
end


text:
"""
"""
struct Composite_map{T}
    _nres1::T
    _nres2::T
    function Composite_map(nres1::T, nres2::T) where {T<:Representation}
        new{T}(nres1, nres2)
    end
end

#is_active(it::Composite_map) = is_active(map_to_SAT(it._nres1)) * is_active(map_to_SAT(it._nres2))

export Composite_map

max_stars_repo_name: leikanger/NRES

max_issues_repo_name: leikanger/NRES

max_forks_repo_name: leikanger/NRES


text:
const PLUGINS = first.(splitext.([
	GMAT_jll.libgmatephempropagator_path,
	GMAT_jll.libgmatformation_path,
	GMAT_jll.libgmatproductionpropagators_path,
	GMAT_jll.libgmatthrustfile_path,
	GMAT_jll.libgmatcinterface_path,
	GMAT_jll.libgmatestimation_path,
	GMAT_jll.libgmatfunction_path,
	GMAT_jll.libgmatsavecommand_path,
	GMAT_jll.libgmatdatainterface_path,
	GMAT_jll.libgmateventlocator_path,
	GMAT_jll.libgmatnewparameters_path,
	GMAT_jll.libgmatscripttools_path,
	GMAT_jll.libgmatyukonoptimizer_path,
	GMAT_jll.libgmatekf_path,
	GMAT_jll.libgmatextrapropagators_path,
	GMAT_jll.libgmatpolyhedrongravity_path,
	GMAT_jll.libgmatstation_path,
]))

function startup_file(data_path, output_path)
	return """
#-------------------------------------------------------------------------------
# General Mission Analysis Tool (GMAT) startup file
#-------------------------------------------------------------------------------
# Comment lines start with #
#
# Path/file naming conventions:
#   - Path name properties end with _PATH and can be referenced elsewhere in the
#     file
#   - File name properties end with _FILE
#   - Path/file names are case sensitive on Linux and OS X
#   - Relative paths are relative to the GMAT executable or application bundle
#     (on OS X)
#
# You can add potential and texture files by following the naming convention:
#   - Potential file properties should begin with celestial body name and end
#     with _POT_FILE
#   - Texture file properties should begin with celestial body name and end
#     with _TEXTURE_FILE
#
# If same _FILE is specified multiple times, the last entry will be used.
#
# You can have multiple GMAT_FUNCTION_PATH, MATLAB_FUNCTION_PATH, or
# PYTHON_MODULE_PATH lines. GMAT will search the paths in top-to-bottom order
# and use the first function file found.
#
# In order for a plugin to work, the path to the plugin library for Windows
# (.dll), Linux (.so) or OS X (.dylib) must be specified in a PLUGIN line below.
# The plugin name should be specified without the file extension. Plugins are
# loaded in the order listed.
#
#===============================================================================

#-----------------------------------------------------------
# Runtime options
#-----------------------------------------------------------
#MATLAB_MODE             = SINGLE
#MATLAB_MODE             = SHARED
#MATLAB_MODE             = NO_MATLAB
#WRITE_GMAT_KEYWORD      = ON
#WRITE_PERSONALIZATION_FILE = ON
#NO_SPLASH               = TRUE
#ECHO_COMMANDS           = TRUE

#-----------------------------------------------------------
# Plugins
#-----------------------------------------------------------
$(join("PLUGIN = " .* PLUGINS, "\n"))

#-----------------------------------------------------------
# Directory/file paths
#-----------------------------------------------------------

# Occasional issues have been observed on Windows in which
# GMAT cannot find certain data files, such as icons and DE
# files, at runtime. If this occurs, set ROOT_PATH to the
# absolute path to your top-level GMAT directory.

ROOT_PATH                = $(pwd())

#-----------------------------------------------------------
OUTPUT_PATH              = $output_path
LOG_FILE                 = OUTPUT_PATH/GmatLog.txt
SCREENSHOT_FILE          = OUTPUT_PATH/GmatScreenShot
#-----------------------------------------------------------
MEASUREMENT_PATH         = OUTPUT_PATH/
#-----------------------------------------------------------
VEHICLE_EPHEM_PATH       = OUTPUT_PATH/
#-----------------------------------------------------------
GMAT_INCLUDE_PATH        = ROOT_PATH/userincludes/
#-----------------------------------------------------------
GMAT_FUNCTION_PATH       = ROOT_PATH/userfunctions/gmat
#-----------------------------------------------------------
MATLAB_FUNCTION_PATH     = ROOT_PATH/matlab
MATLAB_FUNCTION_PATH     = ROOT_PATH/userfunctions/matlab
#-----------------------------------------------------------
PYTHON_MODULE_PATH       = ROOT_PATH/userfunctions/python
#-----------------------------------------------------------
DATA_PATH                = $data_path
FILE_UPDATE_PATH         = $data_path
#-----------------------------------------------------------
PLANETARY_EPHEM_SPK_PATH = DATA_PATH/planetary_ephem/spk/
PLANETARY_SPK_FILE       = PLANETARY_EPHEM_SPK_PATH/DE421AllPlanets.bsp
#-----------------------------------------------------------
PLANETARY_EPHEM_DE_PATH  = DATA_PATH/planetary_ephem/de/
DE405_FILE               = PLANETARY_EPHEM_DE_PATH/leDE1941.405
DE421_FILE               = PLANETARY_EPHEM_DE_PATH/leDE1900.421
DE424_FILE               = PLANETARY_EPHEM_DE_PATH/leDE18002100.424
#-----------------------------------------------------------
VEHICLE_EPHEM_SPK_PATH   = DATA_PATH/vehicle/ephem/spk/
VEHICLE_EPHEM_CCSDS_PATH = DATA_PATH/vehicle/ephem/ccsds/
#-----------------------------------------------------------
PLANETARY_COEFF_PATH     = DATA_PATH/planetary_coeff/
EOP_FILE                 = PLANETARY_COEFF_PATH/eopc04_08.62-now
NUTATION_COEFF_FILE      = PLANETARY_COEFF_PATH/NUTATION.DAT
PLANETARY_PCK_FILE       = PLANETARY_COEFF_PATH/SPICEPlanetaryConstantsKernel.tpc
EARTH_LATEST_PCK_FILE  = PLANETARY_COEFF_PATH/earth_latest_high_prec.bpc
EARTH_PCK_PREDICTED_FILE = PLANETARY_COEFF_PATH/SPICEEarthPredictedKernel.bpc
EARTH_PCK_CURRENT_FILE   = PLANETARY_COEFF_PATH/SPICEEarthCurrentKernel.bpc
LUNA_PCK_CURRENT_FILE    = PLANETARY_COEFF_PATH/SPICELunaCurrentKernel.bpc
LUNA_FRAME_KERNEL_FILE   = PLANETARY_COEFF_PATH/SPICELunaFrameKernel.tf
#-----------------------------------------------------------
TIME_PATH                = DATA_PATH/time/
LEAP_SECS_FILE           = TIME_PATH/tai-utc.dat
LSK_FILE                 = TIME_PATH/SPICELeapSecondKernel.tls
#-----------------------------------------------------------
EARTH_POT_PATH           = DATA_PATH/gravity/earth/
LUNA_POT_PATH            = DATA_PATH/gravity/luna/
MARS_POT_PATH            = DATA_PATH/gravity/mars/
VENUS_POT_PATH           = DATA_PATH/gravity/venus/
OTHER_POT_PATH           = DATA_PATH/gravity/other/
#-----------------------------------------------------------
EGM96_FILE               = EARTH_POT_PATH/EGM96.cof
JGM2_FILE                = EARTH_POT_PATH/JGM2.cof
JGM3_FILE                = EARTH_POT_PATH/JGM3.cof
MARS50C_FILE             = MARS_POT_PATH/Mars50c.cof
MGNP180U_FILE            = VENUS_POT_PATH/MGNP180U.cof
LP165P_FILE              = LUNA_POT_PATH/LP165P.cof
#-----------------------------------------------------------
GUI_CONFIG_PATH          = DATA_PATH/gui_config/
PERSONALIZATION_FILE     = GUI_CONFIG_PATH/MyGmat.ini
#-----------------------------------------------------------
ICON_PATH                = DATA_PATH/graphics/icons/
MAIN_ICON_FILE           = ICON_PATH/GMATWin32.ico
#-----------------------------------------------------------
SPLASH_PATH              = DATA_PATH/graphics/splash/
SPLASH_FILE              = SPLASH_PATH/GMATSplashScreenBETA.png
#-----------------------------------------------------------
TEXTURE_PATH             = DATA_PATH/graphics/texture/
SUN_TEXTURE_FILE         = TEXTURE_PATH/Sun.jpg
MERCURY_TEXTURE_FILE     = TEXTURE_PATH/Mercury_JPLCaltech.jpg
VENUS_TEXTURE_FILE       = TEXTURE_PATH/Venus_BjornJonsson.jpg
EARTH_TEXTURE_FILE       = TEXTURE_PATH/ModifiedBlueMarble.jpg
LUNA_TEXTURE_FILE        = TEXTURE_PATH/Moon_HermesCelestiaMotherlode.jpg
MARS_TEXTURE_FILE        = TEXTURE_PATH/Mars_JPLCaltechUSGS.jpg
JUPITER_TEXTURE_FILE     = TEXTURE_PATH/Jupiter_HermesCelestiaMotherlode.jpg
SATURN_TEXTURE_FILE      = TEXTURE_PATH/Saturn_gradiusCelestiaMotherlode.jpg
URANUS_TEXTURE_FILE      = TEXTURE_PATH/Uranus_JPLCaltech.jpg
NEPTUNE_TEXTURE_FILE     = TEXTURE_PATH/Neptune_BjornJonsson.jpg
PLUTO_TEXTURE_FILE       = TEXTURE_PATH/Pluto_JPLCaltech.jpg
#-----------------------------------------------------------
BODY_3D_MODEL_PATH       = DATA_PATH/graphics/models/
#-----------------------------------------------------------
STAR_PATH                = DATA_PATH/graphics/stars/
STAR_FILE                = STAR_PATH/inp_StarCatalog.txt
CONSTELLATION_FILE       = STAR_PATH/inp_Constellation.txt
BORDER_FILE              = STAR_PATH/inp_Border.txt
#-----------------------------------------------------------
VEHICLE_MODEL_PATH       = DATA_PATH/vehicle/models/
SPACECRAFT_MODEL_FILE    = VEHICLE_MODEL_PATH/aura.3ds
#-----------------------------------------------------------
SPAD_PATH                = DATA_PATH/vehicle/spad/
SPAD_SRP_FILE            = SPAD_PATH/SphericalModel.spo
#-----------------------------------------------------------
ATMOSPHERE_PATH          = DATA_PATH/atmosphere/earth
CSSI_FLUX_FILE           = ATMOSPHERE_PATH/SpaceWeather-All-v1.2.txt
SCHATTEN_FILE            = ATMOSPHERE_PATH/SchattenPredict.txt
MARINI_TROPO_FILE        = ATMOSPHERE_PATH/marini.dat
#-----------------------------------------------------------
HELP_FILE                = ROOT_PATH/docs/help/help.chm
#-----------------------------------------------------------
ICRF_FILE                = DATA_PATH/icrf/ICRF_Table.txt
#-----------------------------------------------------------

#-----------------------------------------------------------
# Debug options
#-----------------------------------------------------------
#RUN_MODE                 = TESTING
#RUN_MODE                 = TESTING_NO_PLOTS
#RUN_MODE                 = EXIT_AFTER_RUN
#DEBUG_MATLAB             = ON
#DEBUG_PARAMETERS         = ON
#DEBUG_FILE_PATH          = ON
#HIDE_SAVEMISSION         = TRUE
#PLOT_MODE                = TILE
"""
end


text:
@testset "Quadratic and hermitian forms" begin
  @time include("QuadForm/Basic.jl")
  @time include("QuadForm/Lattices.jl")
  @time include("QuadForm/Genus.jl")
  @time include("QuadForm/Morphism.jl")
  @time include("QuadForm/Enumeration.jl")
  @time include("QuadForm/MassQuad.jl")
  @time include("QuadForm/MassHerm.jl")
  @time include("QuadForm/Quad.jl")
  @time include("QuadForm/QuadBin.jl")
end


text:
using Test

using StaticArrays
using ExactPredicates
using ExactPredicates.Codegen
import Base: complex

import ExactPredicates: incircle_dbg, orient_dbg

p2(x, y=0.0) = (float(x), float(y))
p3(x, y=0, z=0) = (float(x), float(y), float(z))

@testset "easy" begin
    @test incircle(p2(0), p2(1), p2(0,1), p2(1/2,1/2)) == 1
    @test incircle(p2(0), p2(1), p2(0,1), p2(1,1)) == 0
    @test orient(p2(0.0), p2(1.0), p2(0.0,1.0)) == 1
    @test orient(p3(1.0), p3(0.0,1.0), p3(0,0,1.0), p3(0.0)) == 1
    @test insphere(p3(1.0), p3(0.0,1.0), p3(0,0,1.0), p3(0.0), p3(.5,.5,.5)) == 1
    @test closest(p2(0), p2(1), p2(0,1)) == 1
    @test closest(p3(0), p3(1), p3(0,1)) == 1
end

@testset "easy, with retcode" begin
    @test incircle_dbg(p2(0), p2(1), p2(0,1), p2(1/2,1/2)) == (1, Codegen.fastfp_flt)
    @test incircle_dbg(p2(0), p2(1), p2(0,1), p2(1,1)) == (0, Codegen.interval_flt)
    @test orient_dbg(p2(0.0), p2(1.0), p2(0.0,1.0)) == (1, Codegen.fastfp_flt)
end

@testset "pertubations" begin
    a = p2(0.0)
    b = p2(1.0)
    c = p2(0.0, 1.0)

    @test incircle_dbg(a, b, c, a) == (0, Codegen.zerotest_flt)
    @test incircle_dbg(a, b, c, b) == (0, Codegen.interval_flt)
    @test incircle_dbg(a, b, c, c) == (0, Codegen.interval_flt)
    @test incircle_dbg(a, b, b, c) == (0, Codegen.interval_flt)
    @test incircle(a, b, c, p2(nextfloat(0.0))) == 1
    @test incircle(a, b, c, p2(prevfloat(0.0))) == -1

    a = p2(2.0)
    b = p2(1.0, 1.0)
    c = p2(1.0, - 1.0)

    @test incircle(a, b, c, p2(0.0)) == 0
    @test incircle(a, b, c, p2(nextfloat(0.0))) == 1
    @test incircle(a, b, c, p2(prevfloat(0.0))) == -1
end

@testset "exactly collinear" begin
    for i in 1:10
        rpts = rand(1:2^26, 4)
        pts = (rand(1:2^26) + rand(1:2^26)*im)*rpts
        fpts = convert(Vector{Complex{Float64}}, pts)
        @test incircle_dbg(reinterpret(NTuple{2, Float64}, fpts)...) == (0, Codegen.exact_flt)
    end
end

@testset "exactly collinear (small integers)" begin
    for i in 1:10
        rpts = rand(1:2^5, 4)
        pts = (rand(1:2^5) + rand(1:2^5)*im)*rpts
        fpts = convert(Vector{Complex{Float64}}, pts)
        res, flt = incircle_dbg(reinterpret(NTuple{2, Float64}, fpts)...)
        @test res == 0
        @test flt ∈ [Codegen.zerotest_flt, Codegen.interval_flt]
    end
end

@testset "orient consistency" begin
    for i in 1:10
        a,b,c = randn(ComplexF64)*randn(Float64, 3)
        @test orient(a, b, c) == orient(b, c, a) == orient(c, a, b)
    end
end

@testset "fast zero detection" begin
    a, b, c, d = [reim(p) for p in randn(ComplexF64, 4)]
    @test orient_dbg(a, b, b) == (0, Codegen.zerotest_flt)
    @test orient_dbg(a, b, a) == (0, Codegen.zerotest_flt)
end

@testset "incircle consistency" begin
    for i in 1:10
        a,b,c,d = randn(ComplexF64)*[exp(x*im) for x in randn(Float64, 4)] .+ randn(ComplexF64)
        @test incircle(a,b,c,d) == -incircle(b,c,d,a) == incircle(c,d,a,b) == -incircle(d,a,b,c) == incircle(b,a,d,c) ==
            -incircle(a,c,b,d)
    end
end

@testset "orient3 consistency" begin
    for i in 1:10
        u, v, w = randn(SVector{3, SVector{3, Float64}})
        a = randn()*u + randn()*v + randn()*w
        @test orient(u, v, w, a) == -orient(v, w, a, u) == orient(w, a, u, v) == -orient(a, u, v, w) == orient(v, u, a, w)
    end
end

@testset "underflow incircle" begin
    p = randn(SVector{4, SVector{2, Float64}})

    s = Set()
    while all(abs(c) > floatmin(Float64) for r in p for c in r)
        push!(s, incircle(p...))
        p /= 2
    end

    @test length(s) == 1
end

@testset "overflow incircle" begin
    p = randn(SVector{4, SVector{2, Float64}})

    s = Set()
    while all(isfinite(c) for r in p for c in r)
        push!(s, incircle(p...))
        p *= 2
    end

    @test length(s) == 1
end

@testset "closest point 2d" begin
    a, b, p = (-264.2397483183316, -147.68399381193498),
    (-264.23974834572186, -147.68399391802012),
    (-264.2397478247797, -147.68399318319462)

    @test closest(a, b, p) == 1

    a, b, p =  (-754.9415749355259, 315.57575309957866),
    (-754.9415751576867, 315.5757530485543),
    (-754.9415749822347, 315.5757530017328)

    @test closest(a, b, p) == 1
end

import ExactPredicates: coord
struct Point{T}
    x :: T
    y :: T
end
coord(p :: Point) = (p.x, p.y)

@testset "genericity" begin
    @test incircle(Point(0.0,0.0), Point(1.0,0.0), Point(0.0, 1.0), Point(.5, 0.0)) == 1
end

@testset "parallelorder" begin
    for _ in 1:10
        a, b, c, d = randn(ComplexF64, 4)
        @test orient(a, b, c) == parallelorder(a, b, a, c)
    end
end

@testset "meet" begin
    for _ in 1:10
        t = [exp(im*2*pi*a) for a in sort(rand(4))]
        @test meet(t[1], t[2], t[3], t[4]) == -1
        @test meet(t[1], t[3], t[2], t[4]) == 1
        @test meet(t[1], t[3], t[3], t[4]) == 0
        @test meet(t[1], t[3], t[3], t[1]) == 0
    end

    @test meet(p2(0), p2(1), p2(3), p2(4)) == -1
    @test meet(p2(0), p2(3), p2(1), p2(4)) == 0
    @test meet(p2(1), p2(4), p2(1), p2(2)) == 0
end

@testset "rotation" begin
    angles = 3*pi*2*rand(30)
    sort!(angles)
    pts = [exp(im*α) for α in angles]
    @test rotation(pts) == 3
end

@testset "intersectorder" begin
    a, b, pa, pb, qa, qb = (complex(0.0), complex(1.0),
                            complex(0.0, 3.0), complex(2.0, 1.0),
                            complex(0.0, 2.0), complex(2.0, 1.0))

    @test intersectorder(a, b, pa, pb, qa, qb) == -1
    @test intersectorder(a, b, pb, pa, qa, qb) == 1

    for _ in 1:10
        a, b, pa, pb, qa, qb = randn(ComplexF64, 6)
        @test intersectorder(a, b, pa, pb, qa, qb) == intersectorder(qa, qb, a, b, pa, pb)
    end
end


text:
#Demo 5: Logy plots
#-------------------------------------------------------------------------------
using InspectDR
using Colors

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

default_line = line(color=blue, width=3)

#Input data
#-------------------------------------------------------------------------------
x = collect(1.0:1.0:100.0)
b = 2.0 #Base
k = 0.2 #Exponential coefficient

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Log-Y Tests")
mplot.layout[:ncolumns] = 1

plot = add(mplot, InspectDR.Plot2D(:lin, [:lin, :log10],
		title="y = $b^($(k)x)", xlabel="x", ylabels=["y", "y"])
	)
	strip_liny, strip_logy = plot.strips

	#Define extents & scales:
	strip_logy.grid = InspectDR.GridRect(vmajor=true, hmajor=true, hminor=true)

let wfrm #HIDEWARN_0.7
for i in 1:length(plot.strips)
	wfrm = add(plot, x, b.^(k*x), strip=i)
		wfrm.line = default_line
end
end

gplot = display(InspectDR.GtkDisplay(), mplot)

:DONE


text:
using ApproxFun, Base.Test

## PDEs

d=Interval()^2

x=Fun(identity,d[1]);y=Fun(identity,d[2])

#dirichlet(d) is u[-1,:],u[1,:],u[:,-1],u[:,1]
A=[dirichlet(d);lap(d)]
f=[real(exp(-1+1.0im*y));
                        real(exp(1+1im*y));
                        real(exp(x-1im));
                        real(exp(x+1im))]
S=schurfact(A,100)
@time S=schurfact(A,100)
u=S\f
u=S\f
@time u=S\f;
println("Laplace: should be ~0.014, 0.01")

d=Interval()^2
S=schurfact([neumann(d);lap(d)+100I],100)
@time S=schurfact([neumann(d);lap(d)+100I],100)
u=S\ones(4)
u=S\ones(4)
@time u=S\ones(4)
println("Neumann Helmholtz: should be ~0.016, 0.016")

dx=Interval(0.,1.);dt=Interval(0.0,0.54)
d=dx*dt

V=Fun(x->x^2,dx)

Dt=Derivative(d,[0,1]);Dx=Derivative(d,[1,0])

ϵ=0.1

u0=Fun(x->exp(-25*(x-.5)^2)*exp(-1.0im/(5*ϵ)*log(2cosh(5*(x-.5)))),dx)
L=1im*ϵ*Dt+.5*ϵ^2*Dx^2-V⊗1

PO=discretize([timedirichlet(d);L],50)
@time PO=discretize([timedirichlet(d);L],50)
u=PO\u0
@time    u=PO\u0

println("Schrodinger: should be ~0.013,0.015")