text:
using ThreeBodyDecay
using ThreeBodyDecay.PartialWaveFunctions
using Test

let
    ms = ThreeBodyMasses(m1 = 0.938, m2 = 0.49367, m3 = 0.13957, m0 = 2.46867)
    #
    σs = randomPoint(ms)
    #
    cosζ13_for1_ = cosζ(wr(1,3,1), σs, ms^2)
    cosζ21_for1_ = cosζ(wr(2,1,1), σs, ms^2)
    cosζ23_for1_ = cosζ(wr(2,3,1), σs, ms^2)
    #
    ζ13_for1 = acos(cosζ13_for1_)
    ζ21_for1 = acos(cosζ21_for1_)
    ζ23_for1 = acos(cosζ23_for1_)
    # test of angles
    @test ζ23_for1 ≈ ζ21_for1+ζ13_for1
    #
    # sum rules for the combination of d functions
    prod_13 = [sum(
        wignerd(1,ν,λ,cosζ23_for1_) *
        wignerd(1,λ,0,cosζ21_for1_) * minusone()^λ
        for λ=-1:1) for ν=-1:1]
    just_13 = [wignerd(1,ν,0,cosζ13_for1_) for ν=-1:1]
    @test sum(prod_13 .≈ just_13) == 3
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
"""
Representation of a state for the SingleOCcludedCrosswalk POMDP,
depends on ADM VehicleState type
"""

#### State type

 struct SingleOCFState <: FieldVector{6, Float64}
    ego_y::Float64
    ego_v::Float64
    ped_s::Float64
    ped_T::Float64
    ped_theta::Float64
    ped_v::Float64
end

struct SingleOCFPedState <: FieldVector{4, Float64}
    ped_s::Float64
    ped_T::Float64
    ped_theta::Float64
    ped_v::Float64
end

#### Action type

 struct SingleOCFAction <: FieldVector{2, Float64}
    acc::Float64
    lateral_movement::Float64
end

#### Observaton type

const SingleOCFObs = SingleOCFState

const PEDESTRIAN_OFF_KEY = -1

#### POMDP type

@with_kw mutable struct SingleOCFPOMDP <: POMDP{SingleOCFState, SingleOCFAction, SingleOCFObs}
    env::CrosswalkEnv = CrosswalkEnv()
    PED_SAFETY_DISTANCE::Float64 = 1.5
    ego_type::VehicleDef = VehicleDef()
    ped_type::VehicleDef = VehicleDef(AgentClass.PEDESTRIAN, 1.0 + 2.0*PED_SAFETY_DISTANCE, 1. + 2.0*PED_SAFETY_DISTANCE)
    longitudinal_actions::Vector{Float64} = [1.0, 0.0, -1.0, -2.0, -4.0]
    lateral_actions::Vector{Float64} = [0.0] #[1.0, 0.0, -1.0]
    ΔT::Float64 = 0.2
    PED_A_RANGE::Vector{Float64} = LinRange(-1.0, 1.0, 3)
    PED_THETA_NOISE::Vector{Float64} = LinRange(-0.39/2., 0.39/2., 3)

    EGO_Y_MIN::Float64 = -1.
    EGO_Y_MAX::Float64 = 1.
    EGO_Y_RANGE::Vector{Float64} = LinRange(EGO_Y_MIN, EGO_Y_MAX, 5)

    EGO_V_MIN::Float64 = 0.
    EGO_V_MAX::Float64 = 14.
    EGO_V_RANGE::Vector{Float64} = LinRange(EGO_V_MIN, EGO_V_MAX, 29)

    S_MIN::Float64 = -2.
    S_MAX::Float64 = 50.
    S_RANGE::Vector{Float64} = LinRange(S_MIN, S_MAX, 27)

    T_MIN::Float64 = -5.
    T_MAX::Float64 = 5.
    T_RANGE::Vector{Float64} = LinRange(T_MIN, T_MAX, 11)

    PED_V_MIN::Float64 = 0.
    PED_V_MAX::Float64 = 2.
    PED_V_RANGE::Vector{Float64} = LinRange(PED_V_MIN, PED_V_MAX, 5)

    PED_THETA_MIN::Float64 = 1.57-1.57/2
    PED_THETA_MAX::Float64 = 1.57+1.57/2
    PED_THETA_RANGE::Vector{Float64} = [1.57]# LinRange(PED_THETA_MIN, PED_THETA_MAX, 7)

    state_space_grid::GridInterpolations.RectangleGrid = initStateSpace(EGO_Y_RANGE, EGO_V_RANGE, S_RANGE, T_RANGE, PED_THETA_RANGE, PED_V_RANGE)
    state_space::Vector{SingleOCFState} = getStateSpaceVector(state_space_grid)

    state_space_grid_ped::GridInterpolations.RectangleGrid = initStateSpacePed(S_RANGE, T_RANGE, PED_THETA_RANGE, PED_V_RANGE)
    state_space_ped::Vector{SingleOCFPedState} = getStateSpacePedVector(state_space_grid_ped)

    action_space::Vector{SingleOCFAction} = initActionSpace(longitudinal_actions, lateral_actions)

    ego_vehicle::Vehicle = Vehicle(VehicleState(VecSE2(0.0, 0.0, 0.0), 0.0), VehicleDef(), 1)

    desired_velocity::Float64 = 14. #50.0 / 3.6

    PROBABILITY_PEDESTRIAN_BIRTH::Float64 = 0.6

    COLLISION_COST::Float64 = -500.0
    ACTION_LON_COST::Float64 = -1.0
    ACTION_LAT_COST::Float64 = -1.0
    KEEP_VELOCITY_REWARD::Float64 = 5.0
    KEEP_LANE_REWARD::Float64 = 5.0

    γ::Float64 = 0.95

end

### REWARD MODEL ##################################################################################

function POMDPs.reward(pomdp::SingleOCFPOMDP, s::SingleOCFState, action::SingleOCFAction, sp::SingleOCFState)

    r = 0.

    # is there a collision?
    if collision_checker(pomdp,sp)
        r += (pomdp.COLLISION_COST - sp.ego_v * 20)
        return r
    end

    # keep velocity
    if (sp.ego_v <= pomdp.desired_velocity )
        v_factor = (pomdp.desired_velocity-abs(sp.ego_v - pomdp.desired_velocity)) / pomdp.desired_velocity
        r += (v_factor * pomdp.KEEP_VELOCITY_REWARD)
    else
        r -= pomdp.KEEP_VELOCITY_REWARD
    end

   # stay in center of the road
    if ( abs(sp.ego_y) <= pomdp.EGO_Y_MAX )
        l_factor = 1 - abs(sp.ego_y)  / pomdp.EGO_Y_MAX
        r += (l_factor * pomdp.KEEP_LANE_REWARD)
    else
        r -= pomdp.KEEP_LANE_REWARD
    end

    # costs for longitudinal actions
    if action.acc > 0. ||  action.acc < 0.0
        r += pomdp.ACTION_LON_COST * abs(action.acc)
    end

    # costs for lateral actions
    if abs(action.lateral_movement) > 0
        r += pomdp.ACTION_LAT_COST * abs(action.lateral_movement)
    end

    return r

end

### HELPERS

function POMDPs.isterminal(pomdp::SingleOCFPOMDP, s::SingleOCFState)

    if collision_checker(pomdp,s)
        return true
    end

    return false
end

function POMDPs.discount(pomdp::SingleOCFPOMDP)
    return pomdp.γ
end

function AutomotivePOMDPs.collision_checker(pomdp::SingleOCFPOMDP, s::SingleOCFState)

    object_a_def = pomdp.ego_type
    object_b_def = pomdp.ped_type

    center_a = VecSE2(-object_a_def.length/2, s.ego_y, 0.0)
    center_b = VecSE2(s.ped_s, s.ped_T, s.ped_theta)

    # first fast check:
    @fastmath begin
        Δ = sqrt((center_a.x - center_b.x)^2 + (center_a.y - center_b.y)^2)
        r_a = sqrt(object_a_def.length*object_a_def.length/4 + object_a_def.width*object_a_def.width/4)
        r_b = sqrt(object_b_def.length*object_b_def.length/4 + object_b_def.width*object_b_def.width/4)
    end
    if Δ ≤ r_a + r_b
        # fast check is true, run parallel axis theorem
        Pa = AutomotiveDrivingModels.polygon(center_a, object_a_def)
        Pb = AutomotiveDrivingModels.polygon(center_b, object_b_def)
        return AutomotiveDrivingModels.overlap(Pa, Pb)
    end
    return false
end

function initStateSpace(EGO_Y_RANGE, EGO_V_RANGE, S_RANGE, T_RANGE, PED_THETA_RANGE, PED_V_RANGE)

    return RectangleGrid(EGO_Y_RANGE, EGO_V_RANGE, S_RANGE, T_RANGE, PED_THETA_RANGE, PED_V_RANGE)
end

function initStateSpacePed(S_RANGE, T_RANGE, PED_THETA_RANGE, PED_V_RANGE)

    return RectangleGrid(S_RANGE, T_RANGE, PED_THETA_RANGE, PED_V_RANGE)

end

function getStateSpacePedVector(ped_grid)

    state_space_ped = SingleOCFPedState[]

    for i = 1:length(ped_grid)
        s = ind2x(ped_grid,i)
        push!(state_space_ped,SingleOCFPedState(s[1], s[2], s[3], s[4]))
    end
    # add absent state
    push!(state_space_ped,SingleOCFPedState(-10., -10., 0., 0.))
    return state_space_ped
end

function getStateSpaceVector(grid_space)

    state_space = SingleOCFState[]

    for i = 1:length(grid_space)
        s = ind2x(grid_space,i)
        push!(state_space,SingleOCFState(s[1], s[2], s[3], s[4], s[5], s[6]))
    end

    # add absent state
    push!(state_space,SingleOCFState(0., 0., -10., -10., 0., 0.))
    return state_space
end

function initActionSpace(longitudinal_actions, lateral_actions)

  action_space = SingleOCFAction[]

  for lat_a in lateral_actions
    for lon_a in longitudinal_actions
      push!(action_space, SingleOCFAction(lon_a, lat_a))
    end
  end

  return action_space

end

ext: jl


text:
jtype(tag::BSONType)::DataType =
  tag == null ? Nothing :
  tag == boolean ? Bool :
  tag == int32 ? Int32 :
  tag == int64 ? Int64 :
  tag == double ? Float64 :
  error("Unsupported tag $tag")

parse_cstr(io::IO) = readuntil(io, '\0')

function parse_tag(io::IO, tag::BSONType)
  if tag == null
    nothing
  elseif tag == document
    parse_doc(io)
  elseif tag == array
    parse_array(io)
  elseif tag == string
    len = read(io, Int32)-1
    s = String(read(io, len))
    eof = read(io, 1)
    s
  elseif tag == binary
    len = read(io, Int32)
    subtype = read(io, 1)
    read(io, len)
  else
    read(io, jtype(tag))
  end
end

function parse_array(io::IO)::BSONArray
  len = read(io, Int32)
  ps = BSONArray()

  while (tag = read(io, BSONType)) ≠ eof
    # Note that arrays are dicts with the index as the key
    while read(io, UInt8) != 0x00
      nothing
    end
    push!(ps, parse_tag(io::IO, tag))
  end

  ps
end

function parse_doc(io::IO)::BSONDict
  len = read(io, Int32)
  dic = BSONDict()

  while (tag = read(io, BSONType)) ≠ eof
    k = Symbol(parse_cstr(io))
    dic[k] = parse_tag(io::IO, tag)
  end

  dic
end

backrefs!(x, refs) = applychildren!(x -> backrefs!(x, refs), x)

backrefs!(dict::BSONDict, refs) =
  get(dict, :tag, "") == "backref" ? refs[dict[:ref]] :
  invoke(backrefs!, Tuple{Any,Any}, dict, refs)

function backrefs!(dict)
  haskey(dict, :_backrefs) || return dict
  refs = dict[:_backrefs]
  backrefs!(dict, refs)
  delete!(dict, :_backrefs)
  return dict
end

const tags = Dict{Symbol,Function}()

const raise = Dict{Symbol,Function}()

function _raise_recursive(d::AbstractDict, cache, init)
  if haskey(d, :tag) && haskey(tags, Symbol(d[:tag]))
    if Symbol(d[:tag]) in (:ref, :datatype)
      cache[d] = tags[Symbol(d[:tag])](applychildren!(x -> raise_recursive(x, cache, init), d), init)
    else
      cache[d] = tags[Symbol(d[:tag])](applychildren!(x -> raise_recursive(x, cache, init), d))
    end
  else
    cache[d] = d
    applychildren!(x -> raise_recursive(x, cache, init), d)
  end
end

function raise_recursive(d::AbstractDict, cache, init)
  haskey(cache, d) && return cache[d]
  haskey(d, :tag) && haskey(raise, Symbol(d[:tag])) && return raise[Symbol(d[:tag])](d, cache, init)
  _raise_recursive(d::AbstractDict, cache, init)
end

function raise_recursive(v::BSONArray, cache, init)
  cache[v] = v
  applychildren!(x -> raise_recursive(x, cache, init), v)
end

raise_recursive(x, cache, init) = x

raise_recursive(x, init) = raise_recursive(x, IdDict(), init)

parse(io::IO) = backrefs!(parse_doc(io))
parse(path::String) = open(parse, path)

load(x, init=Main) = raise_recursive(parse(x), init)

function roundtrip(x)
  buf = IOBuffer()
  bson(buf, Dict(:data => x))
  load(seek(buf, 0))[:data]
end


text:
function matgen(n)
  tmp = 1.0 / n / n
  [ tmp * (i - j) * (i + j - 2) for i=1:n, j=1:n ]
end

function mul(a, b)
  m,n = size(a)
  q,p = size(b)
  @assert n == q

  # transpose a for cache-friendliness
  aT = zeros(n,m)
  @simd for i = 1:m
    for j = 1:n
      @inbounds aT[j,i] = a[i,j]
    end
  end

  out = zeros(m,p)
  @simd for i = 1:m
    for j = 1:p
      z = 0.0
      for k = 1:n
        @inbounds z += aT[k,i]*b[k,j]
      end
      @inbounds out[i,j] = z
    end
  end
  out
end

function main(n)
  t = time()
  n = round(Int, n / 2 * 2)
  a = matgen(n)
  b = matgen(n)
  c = mul(a, b)
  v = round(Int, n/2) + 1
  println(c[v, v])
  println(time() - t)
end

function test()
  n = 100
  if length(ARGS) >= 1
    n = parse(Int, ARGS[1])
  end

  println(STDERR, "warming")
  main(200)

  println(STDERR, "bench")
  x = @timed main(n)
  println(STDERR, "Elapsed: $(x[2]), Allocated: $(x[3]), GC Time: $(x[4])")
end

test()

max_stars_repo_name: goddva/benchmarks

max_issues_repo_name: goddva/benchmarks

max_forks_repo_name: goddva/benchmarks


text:
function find_nums(limit)
    count = Dict()
    total = 0

    for i in 2:Int(floor((limit-4)^(1/3)))
        for j in 2:Int(floor((limit-i^3)^(1/2)))
            n = i^3 + j^2
            count[n] = get(count, n, 0) + 1
        end
    end

    for (k, v) in count
        if v == 4 && is_palindrome(string(k))
            total += k
        end
    end

    println(total)
end

function is_palindrome(n)
    k = div(length(n), 2)

    for i in 1:k
        if n[i] != n[end-(i-1)]
            return false
        end
    end

    return true
end

find_nums(Int(1e9))

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
# https://github.com/spencerlyon2/PlotlyJS.jl

supportedArgs(::PlotlyJSBackend) = [
    :annotation,
    # :axis,
    :background_color,
    :color_palette,
    :fillrange,
    :fillcolor,
    :fillalpha,
    :foreground_color,
    :group,
    :label,
    :layout,
    :legend,
    :seriescolor, :seriesalpha,
    :linecolor,
    :linestyle,
    :linetype,
    :linewidth,
    :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    :markeralpha,
    :markerstrokewidth,
    :markerstrokecolor,
    :markerstrokestyle,
    :n,
    :bins,
    :nc,
    :nr,
    # :pos,
    # :smooth,
    :show,
    :size,
    :title,
    :windowtitle,
    :x,
    :xlabel,
    :xlims,
    :xticks,
    :y,
    :ylabel,
    :ylims,
    # :yrightlabel,
    :yticks,
    :xscale,
    :yscale,
    :xflip,
    :yflip,
    :z,
    :marker_z,
    :tickfont,
    :guidefont,
    :legendfont,
    :grid,
    :levels,
    :xerror,
    :yerror,
    :ribbon,
    :quiver,
    :orientation,
    :polar,
  ]
supportedAxes(::PlotlyJSBackend) = [:auto, :left]
supportedTypes(::PlotlyJSBackend) = [:none, :line, :path, :scatter, :steppre, :steppost,
                                   :hist2d, :hist, :density, :bar, :contour, :surface, :path3d, :scatter3d,
                                   :pie, :heatmap] #,, :sticks, :hexbin, :hline, :vline]
supportedStyles(::PlotlyJSBackend) = [:auto, :solid, :dash, :dot, :dashdot]
supportedMarkers(::PlotlyJSBackend) = [:none, :auto, :ellipse, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross,
                                     :pentagon, :hexagon, :octagon, :vline, :hline] #vcat(_allMarkers, Shape)
supportedScales(::PlotlyJSBackend) = [:identity, :log10] #, :ln, :log2, :log10, :asinh, :sqrt]
subplotSupported(::PlotlyJSBackend) = true
stringsSupported(::PlotlyJSBackend) = true

# --------------------------------------------------------------------------------------

function _initialize_backend(::PlotlyJSBackend; kw...)
    @eval begin
        import PlotlyJS
        export PlotlyJS
    end

    for (mime, fmt) in PlotlyJS._mimeformats
        # mime == "image/png" && continue  # don't use plotlyjs's writemime for png
        @eval Base.writemime(io::IO, m::MIME{symbol($mime)}, p::Plot{PlotlyJSBackend}) = writemime(io, m, p.o)
    end

    # override IJulia inline display
    if isijulia()
        IJulia.display_dict(plt::AbstractPlot{PlotlyJSBackend}) = IJulia.display_dict(plt.o)
    end
end

# ---------------------------------------------------------------------------

function _create_plot(pkg::PlotlyJSBackend, d::KW)
    # TODO: create the window/canvas/context that is the plot within the backend (call it `o`)
    # TODO: initialize the plot... title, xlabel, bgcolor, etc
    # o = PlotlyJS.Plot(PlotlyJS.GenericTrace[], PlotlyJS.Layout(),
    #                   Base.Random.uuid4(), PlotlyJS.ElectronDisplay())
    # T = isijulia() ? PlotlyJS.JupyterPlot : PlotlyJS.ElectronPlot
    # o = T(PlotlyJS.Plot())
    o = PlotlyJS.plot()

    Plot(o, pkg, 0, d, KW[])
end

function _add_series(::PlotlyJSBackend, plt::Plot, d::KW)
    syncplot = plt.o

    # add to the data array
    pdict = plotly_series(d, plt.plotargs)
    typ = pop!(pdict, :type)
    gt = PlotlyJS.GenericTrace(typ; pdict...)
    PlotlyJS.addtraces!(syncplot, gt)

    push!(plt.seriesargs, d)
    plt
end

# ---------------------------------------------------------------------------

function _add_annotations{X,Y,V}(plt::Plot{PlotlyJSBackend}, anns::AVec{@compat(Tuple{X,Y,V})})
    # set or add to the annotation_list
    if !haskey(plt.plotargs, :annotation_list)
        plt.plotargs[:annotation_list] = Any[]
    end
    append!(plt.plotargs[:annotation_list], anns)
end

# ----------------------------------------------------------------

function _before_update_plot(plt::Plot{PlotlyJSBackend})
end

# TODO: override this to update plot items (title, xlabel, etc) after creation
function _update_plot(plt::Plot{PlotlyJSBackend}, d::KW)
    pdict = plotly_layout(d)
    syncplot = plt.o
    w,h = d[:size]
    PlotlyJS.relayout!(syncplot, pdict, width = w, height = h)
end

function _update_plot_pos_size(plt::AbstractPlot{PlotlyJSBackend}, d::KW)
end

# ----------------------------------------------------------------

# accessors for x/y data

# function getxy(plt::Plot{PlotlyJSBackend}, i::Int)
#   d = plt.seriesargs[i]
#   d[:x], d[:y]
# end

function setxy!{X,Y}(plt::Plot{PlotlyJSBackend}, xy::Tuple{X,Y}, i::Integer)
  d = plt.seriesargs[i]
  ispolar = get(plt.plotargs, :polar, false)
  xsym = ispolar ? :t : :x
  ysym = ispolar ? :r : :y
  d[xsym], d[ysym] = xy
  # TODO: this is likely ineffecient... we should make a call that ONLY changes the plot data
  PlotlyJS.restyle!(plt.o, i, KW(xsym=>(d[xsym],), ysym=>(d[ysym],)))
  plt
end

# ----------------------------------------------------------------

function _create_subplot(subplt::Subplot{PlotlyJSBackend}, isbefore::Bool)
  # TODO: build the underlying Subplot object.  this is where you might layout the panes within a GUI window, for example
  true
end

function _expand_limits(lims, plt::Plot{PlotlyJSBackend}, isx::Bool)
  # TODO: call expand limits for each plot data
end

function _remove_axis(plt::Plot{PlotlyJSBackend}, isx::Bool)
  # TODO: if plot is inner subplot, might need to remove ticks or axis labels
end

# ----------------------------------------------------------------

# function Base.writemime(io::IO, m::MIME"text/html", plt::AbstractPlot{PlotlyJSBackend})
#     Base.writemime(io, m, plt.o)
# end

# function Base.writemime(io::IO, ::MIME"image/png", plt::AbstractPlot{PlotlyJSBackend})
#     println("here!")
#     writemime_png_from_html(io, plt)
# end

function Base.display(::PlotsDisplay, plt::Plot{PlotlyJSBackend})
    display(plt.o)
end

function Base.display(::PlotsDisplay, plt::Subplot{PlotlyJSBackend})
    error()
end


text:
module HyperoctahedralGroup

using FileIO, Permutations, StaticArrays

export get_elements, get_subgroupdict, get_subgroup_properties
export ⊗

const xyp = SA[0 -1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1]
const xym = SA[0 1 0 0; -1 0 0 0; 0 0 1 0; 0 0 0 1]

const xzp = SA[0 0 -1 0; 0 1 0 0; 1 0 0 0; 0 0 0 1]
const xzm = SA[0 0 1 0; 0 1 0 0; -1 0 0 0; 0 0 0 1]

const xwp = SA[0 0 0 -1; 0 1 0 0; 0 0 1 0; 1 0 0 0]
const xwm = SA[0 0 0 1; 0 1 0 0; 0 0 1 0; -1 0 0 0]

const yzp = SA[1 0 0 0; 0 0 -1 0; 0 1 0 0; 0 0 0 1]
const yzm = SA[1 0 0 0; 0 0 1 0; 0 -1 0 0; 0 0 0 1]

const ywp = SA[1 0 0 0; 0 0 0 -1; 0 0 1 0; 0 1 0 0]
const ywm = SA[1 0 0 0; 0 0 0 1; 0 0 1 0; 0 -1 0 0]

const zwp = SA[1 0 0 0; 0 1 0 0; 0 0 0 -1; 0 0 1 0]
const zwm = SA[1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 -1 0]

function get_elements()
    let elements_file = joinpath(@__DIR__,"../deps/elements.jld2")
        if isfile(elements_file)
            elements = load(elements_file,"elements")
            return elements
        else
            error("Build the package first!")
        end
    end
end

function order(e)
    r = e*e
    n = 1
    while r != e
        r = r*e
        n+=1
    end
    return n
end

include(joinpath(@__DIR__, "../deps/groupmatrix.jl"))

function ⊗(a::Integer,b::Integer)
    if 192 >= a > 0 && 192 >= b > 0
        return g[a,b]
    else
        error("$a o $b Not defined.")
    end
end

function insertsorted!(set,element;startindex = 1)
    if !issorted(set)
        error("Array is not sorted.")
    end
    i = startindex
    while  i <= length(set) && element > set[i]
        i += 1
    end
    i > length(set) ? push!(set, element) : insert!(set, i, element)
    return i
end

function complete!(set)
    added = -1
    queue = vec([(a,b) for a in set, b in set])
    while added != 0
        added = 0
        newElements = empty(set)
        for i in length(queue):-1:1
            a,b = queue[i]
            if ⊗(a,b) ∉ set
                push!(newElements,⊗(a,b))
                added = added+1
            end
            pop!(queue)
        end
        newElements = sort(unique(newElements))
        for a in newElements
            for b in set
                push!(queue,(a,b))
                push!(queue,(b,a))
            end
            for b in newElements
                push!(queue,(a,b))
            end
        end

        i = 1
        for a in newElements
            i = insertsorted!(set,a, startindex = i)
        end
    end
    return set
end

function completewith!(set,toadd)
    toadd = setdiff(toadd,set)
    isempty(toadd) && return set
    added = -1
    queue = vcat(vec([(a,b) for a in set, b in toadd]),vec([(b,a) for a in set, b in toadd]),vec([(a,b) for a in toadd, b in toadd]))
    for a in toadd
        insertsorted!(set,a)
    end
    while added != 0
        added = 0
        newElements = empty(set)
        for i in length(queue):-1:1
            a,b = queue[i]
            if ⊗(a,b) ∉ set
                push!(newElements,⊗(a,b))
                added = added+1
            end
            pop!(queue)
        end
        newElements = sort(unique(newElements))
        for a in newElements
            for b in set
                push!(queue,(a,b))
                push!(queue,(b,a))
            end
            for b in newElements
                push!(queue,(a,b))
            end
        end

        i = 1
        for a in newElements
            i = insertsorted!(set,a, startindex = i)
        end
    end
    return set
end

function isgroup(set)
    for a in set
        for b in set
            if ⊗(a,b) ∉ set
                return false
            end
        end
    end
    return true
end

function get_subgroupdict()
    let subgroupdict_file = joinpath(@__DIR__,"../deps/subgroupdict.jld2")
        if isfile(subgroupdict_file)
            subgroupdict = load(subgroupdict_file,"subgroupdict")
            return subgroupdict
        else
            error("Build the package first!")
        end
    end
end

function get_subgroup_properties()
    let subgroup_properties_file = joinpath(@__DIR__,"../deps/subgroup_properties.jld2")
        if isfile(subgroup_properties_file)
            subgroup_properties = load(subgroup_properties_file,"subgroup_properties")
            return subgroup_properties
        else
            error("Build the package first!")
        end
    end
end

end #of module


text:
"""
Gaussian Mixture Models
"""
module GMM

import Distributions
using DocStringExtensions
import SpecialFunctions

function normal_cdf(
    x::Vector{Float64}, mu::Vector{Float64}, sigma::Vector{Float64}
)::Vector{Float64}

    if length(x) != length(mu) && length(x) != 1 && length(mu) != 1
        throw(DimensionMismatch("x: $(x) and mu: $(mu) are not the same length or 1"))
    end
    top = x .- mu
    bottom = max.(sqrt(2) .* sigma, eps(Float64))
    z = top ./ bottom

    return (1 .+ SpecialFunctions.erf.(z)) ./ 2
end
"""
$(TYPEDSIGNATURES)

Bounded Gaussian Mixture Model (BGMM)
Sample from truncated 1-D Gaussian Mixture Model

GMM1 with low, high
"""
function GMM1(
    weights::Vector{Float64}, mus::Vector{Float64}, sigmas::Vector{Float64},
    low::Float64, high::Float64, sample_size::Int
)::Vector{Float64}
    if low > high
        throw(ArgumentError(string(
            "low (", low, ") should not be greater than high ", high
        )))
    end
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas)
        )))
    end
    if low >= high
        throw(ArgumentError("low is greater or equal to high, low: $(low), high: $(high)"))
    end
    samples = Float64[]
    while length(samples) < sample_size
        active = argmax(rand(Distributions.Multinomial(1, weights)))
        draw = rand(Distributions.Normal(mus[active], sigmas[active]))
        if low <= draw < high
            push!(samples, draw)
        end
    end

    return samples
end
"""
$(TYPEDSIGNATURES)
GMM1 without low, high or q
"""
function GMM1(
    weights::Vector{Float64}, mus::Vector{Float64}, sigmas::Vector{Float64},
    sample_size::Int
)::Vector{Float64}

    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end
    active_arr = getindex.(
        permutedims(
            argmax(rand(Distributions.Multinomial(1, weights), sample_size); dims = 1)
        ),
        1
    )
    samples = rand.(Distributions.Normal.(mus[active_arr], sigmas[active_arr]))

    return reshape(samples, sample_size)
end
"""
$(TYPEDSIGNATURES)
GMM1 with low, high and q
"""
function GMM1(
    weights::Vector{Float64}, mus::Vector{Float64}, sigmas::Vector{Float64},
    low::Float64, high::Float64, q::Float64, sample_size::Int
)::Vector{Float64}
    samples = GMM1(weights, mus, sigmas, low, high, sample_size)

    return round.(samples ./ q) .* q
end
"""
$(TYPEDSIGNATURES)
GMM1 with q
"""
function GMM1(
    weights::Vector{Float64}, mus::Vector{Float64}, sigmas::Vector{Float64},
    q::Float64, sample_size::Int
)::Vector{Float64}
    samples = GMM1(weights, mus, sigmas, sample_size)

    return round.(samples ./ q) .* q
end

# GMM1_lpdf

function logprob(
    samples::Vector{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}, low::Float64, high::Float64, q::Float64, p_accept::Float64
)::Vector{Float64}
    prob = zeros(size(samples))
    for (w, mu, sigma) in zip(weights, mus, sigmas)
        ubound = min.(samples .+ (q / 2.0), high)
        lbound = max.(samples .- (q / 2.0), low)
        # -- two-stage addition is slightly more numerically accurate
        inc_amt = w .* normal_cdf(ubound, [mu], [sigma])
        inc_amt .-= w .* normal_cdf(lbound, [mu], [sigma])
        prob += inc_amt
    end

    return log.(prob) .- log(p_accept)
end
function logprob(
    samples::Vector{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}, q::Float64, p_accept::Float64
)::Vector{Float64}
    prob = zeros(size(samples))
    for (w, mu, sigma) in zip(weights, mus, sigmas)
        ubound = samples .+ (q / 2.0)
        lbound = samples .- (q / 2.0)
        # -- two-stage addition is slightly more numerically accurate
        inc_amt = w .* normal_cdf(ubound, [mu], [sigma])
        inc_amt .-= w .* normal_cdf(lbound, [mu], [sigma])
        prob += inc_amt
    end

    return log.(prob) .- log(p_accept)
end

function logsum_rows(x::Matrix{Float64})::Vector{Float64}
    m = maximum(x, dims = 2)

    return (log.(sum(exp.(x .- m); dims = 2)) .+ m)[:]
end

function mahal(
    samples::Vector{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}, p_accept::Float64
)::Vector{Float64}
    dist = reshape(samples, length(samples), 1) .- reshape(mus, 1, length(mus))
    # mahal size is (n_samples, n_components)
    mahal = (dist ./ reshape(max.(sigmas, eps(Float64)), 1, length(sigmas))) .^ 2
    Z = sqrt.(2pi .* (sigmas .^ 2))
    coef = weights ./ Z ./ p_accept

    return logsum_rows(-0.5 .* mahal .+ log.(reshape(coef, 1, length(coef))))
end

"""
$(TYPEDSIGNATURES)
GMM1_lpdf with low, high and q
"""
function GMM1_lpdf(
    samples::Vector{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}, low::Float64, high::Float64, q::Float64
)::Vector{Float64}
    isempty(samples) && return []
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end

    p_accept = sum(weights .* (normal_cdf([high], mus, sigmas) - normal_cdf([low], mus, sigmas)))
    return logprob(samples, weights, mus, sigmas, low, high, q, p_accept)
end
"""
$(TYPEDSIGNATURES)
GMM1_lpdf with q
"""
function GMM1_lpdf(
    samples::Vector{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}, q::Float64
)::Vector{Float64}
    isempty(samples) && return []
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end
    p_accept = 1.0

    return logprob(samples, weights, mus, sigmas, q, p_accept)
end
"""
$(TYPEDSIGNATURES)
GMM1_lpdf with low, high
"""
function GMM1_lpdf(
    samples::Vector{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}, low::Float64, high::Float64
)::Vector{Float64}
    isempty(samples) && return []
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end
    p_accept = sum(
        weights .* (normal_cdf([high], mus, sigmas) - normal_cdf([low], mus, sigmas))
    )

    return mahal(samples, weights, mus, sigmas, p_accept)
end
"""
$(TYPEDSIGNATURES)
GMM1_lpdf without low, high or q
"""
function GMM1_lpdf(
    samples::Vector{Float64}, weights::Vector{Float64}, mus::Vector{Float64},
    sigmas::Vector{Float64}
)::Vector{Float64}
    isempty(samples) && return []
    if !(length(weights) == length(mus) == length(sigmas))
        throw(DimensionMismatch(string(
            "length(weights): ", length(weights),
            " doesn't equal length(mus): ", length(mus),
            " nor length(sigmas): ", length(sigmas),
        )))
    end

    return mahal(samples, weights, mus, sigmas, 1.0)
end

end # module GMM


text:
using DDF

using Bernstein
using ComputedFieldTypes
using DifferentialForms
using LinearAlgebra
using StaticArrays

@testset "Barycentric basis functions D=$D P=$P R=$R" for D in 0:Dmax,
P in (Pr, Dl),
R in 0:D

    # TODO: all P
    P == Pr || continue

    N = D + 1
    S = Float64
    mfd = simplex_manifold(Val(D), S)
    d = R == 0 ? nothing : deriv(Val(P), Val(R - 1), mfd)

    e = Form{D + 1,0}((1,))

    λs = [SVector{N,S}(i.I) / sum(i.I)
          for i in
              CartesianIndex(ntuple(d -> 0, N)):CartesianIndex(ntuple(d -> 1,
                                                                      N))
          if sum(i.I) > 0]
    λs::Vector{SVector{N,S}}
    @assert all(sum(λ) ≈ 1 for λ in λs)

    basis_λ = Continuum.basis_λ

    if D == 0
        # vertices: 1

        e1 = Form{D + 1,1}((1,))

        if R == 0
            bs = [λ -> λ[1] * e]
        end
    elseif D == 1
        # vertices: 1 2
        # edges:    (1,2)
        @test collect(sparse_column_rows(get_simplices(mfd, 0), ID{0}(1))) ==
              [ID{0}(1)]
        @test collect(sparse_column_rows(get_simplices(mfd, 0), ID{0}(2))) ==
              [ID{0}(2)]
        @test collect(sparse_column_rows(get_simplices(mfd, 1), ID{1}(1))) ==
              [ID{0}(1), ID{0}(2)]

        e1 = Form{D + 1,1}((1, 0))
        e2 = Form{D + 1,1}((0, 1))

        if R == 0
            e = Form{D + 1,R}((1,))
            bs = [λ -> λ[1] * e, λ -> λ[2] * e]
        elseif R == 1
            @test d.values == [-1 1]
            e1 = Form{D + 1,R}((1, 0))
            e2 = Form{D + 1,R}((0, 1))
            bs = [λ -> λ[1] * e2 - λ[2] * e1]
        end
    elseif D == 2
        # vertices: 1 2 3
        # edges:    (1,2) (1,3) (2,3)
        # faces:    (1,2,3)
        @test collect(sparse_column_rows(get_simplices(mfd, 0), ID{0}(1))) ==
              [ID{0}(1)]
        @test collect(sparse_column_rows(get_simplices(mfd, 0), ID{0}(2))) ==
              [ID{0}(2)]
        @test collect(sparse_column_rows(get_simplices(mfd, 0), ID{0}(3))) ==
              [ID{0}(3)]
        @test collect(sparse_column_rows(get_simplices(mfd, 1), ID{1}(1))) ==
              [ID{0}(1), ID{0}(2)]
        @test collect(sparse_column_rows(get_simplices(mfd, 1), ID{1}(2))) ==
              [ID{0}(1), ID{0}(3)]
        @test collect(sparse_column_rows(get_simplices(mfd, 1), ID{1}(3))) ==
              [ID{0}(2), ID{0}(3)]
        @test collect(sparse_column_rows(get_simplices(mfd, 2), ID{2}(1))) ==
              [ID{0}(1), ID{0}(2), ID{0}(3)]

        e1 = Form{D + 1,1}((1, 0, 0))
        e2 = Form{D + 1,1}((0, 1, 0))
        e3 = Form{D + 1,1}((0, 0, 1))

        if R == 0
            bs = [λ -> λ[1] * e, λ -> λ[2] * e, λ -> λ[3] * e]
        elseif R == 1
            @test d.values == [
                -1 1 0
                -1 0 1
                0 -1 1
            ]
            bs = [λ -> λ[1] * e2 - λ[2] * e1, λ -> λ[1] * e3 - λ[3] * e1,
                  λ -> λ[2] * e3 - λ[3] * e2]
        elseif R == 2
            @test d.values == [1 -1 1]
            e12 = e1 ∧ e2
            e13 = e1 ∧ e3
            e23 = e2 ∧ e3
            bs = [λ -> λ[1] * e23 - λ[2] * e13 + λ[3] * e12]
        end
    else
        bs = []
    end

    for n in 1:length(bs), λ in λs
        @test basis_λ(Form{D,R}, n, λ) == bs[n](λ)
    end
end

@testset "Cartesian basis functions D=$D P=$P R=$R" for D in 0:Dmax,
P in (Pr, Dl),
R in 0:D

    # TODO: all P
    P == Pr || continue

    S = Float64
    mfd = orthogonal_simplex_manifold(Val(D), S; optimize_mesh=false)

    e = Form{D,0}((1,))

    xs = SVector{D + 1,SVector{D,S}}(get_coords(mfd)[i]
                                     for i in ID{0}(1):ID{0}(D + 1))
    x2λ = cartesian2barycentric_setup(xs)
    dλ2dx = Continuum.dbarycentric2dcartesian_setup(Form{D,R}, xs)

    if D == 0
        pts = [SVector{D,S}()]
    else
        pts = [sum(i.I .* xs) / sum(i.I)
               for i in
                   CartesianIndex(ntuple(d -> 0, D + 1)):CartesianIndex(ntuple(d -> 1,
                                                                               D +
                                                                               1))
               if sum(i.I) > 0]
    end
    pts::Vector{SVector{D,S}}

    basis_x = Continuum.basis_x

    # Test basis functions at various points
    if D == 0
        if R == 0
            bs = [x -> e]
        end
    elseif D == 1
        e1 = Form{D,1}((1,))
        if R == 0
            bs = [x -> (1 - x[1]) * e, x -> x[1] * e]
        elseif R == 1
            bs = [x -> e1]
        end
    elseif D == 2
        e1 = Form{D,1}((1, 0))
        e2 = Form{D,1}((0, 1))
        if R == 0
            e = Form{D,R,S}((1,))
            bs = [x -> (1 - x[1] - x[2]) * e, x -> x[1] * e, x -> x[2] * e]
        elseif R == 1
            bs = [x -> -(x[2] - 1) * e1 + x[1] * e2,
                  x -> x[2] * e1 - (x[1] - 1) * e2, x -> -x[2] * e1 + x[1] * e2]
        elseif R == 2
            e12 = e1 ∧ e2
            bs = [x -> e12]
        end
    else
        bs = []
    end

    for n in 1:length(bs), x in pts
        @test basis_x(Form{D,R}, x2λ, dλ2dx, n, x) ≈ bs[n](x)
    end
end

# @DISABLED @testset "Cartesian basis functions D=$D P=$P R=$R" for D in
#                                                                          0:Dmax,
# P in (Pr, Dl),
# R in 0:D
#
#     # TODO: all P
#     P == Pr || continue
#
#     S = Float64
#     mfd = simplex_manifold(Val(D), S)
#
#     xs = SVector{D + 1,SVector{D,S}}(get_coords(mfd)[i]
#                                      for i in ID{0}(1):ID{0}(D + 1))
#     x2λ = cartesian2barycentric_setup(xs)
#     dλ2dx = Continuum.dbarycentric2dcartesian_setup(Form{D,R}, xs)
#
#     VR = IDVector{R}
#     I0(is...) = map(ID{0}, is)
#
#     for n in axes(get_simplices(mfd, P == Pr ? R : D - R), 2)
#         for i in axes(get_simplices(mfd, 0), 1)
#             x = get_coords(mfd)[i]
#             bx = Continuum.basis_x(Form{D,R}, x2λ, dλ2dx, Int(n), x)
#             if R == 0
#                 @test bx ≈
#                       Form{D,R}((Int(i) == Int(n) ? get_volumes(mfd, R)[n] : 0,))
#             elseif R == D
#                 if !(bx ≈ Form{D,R}((get_volumes(mfd, R)[n],)))
#                     @show D R n i bx
#                     @show get_volumes(mfd, R)
#                 end
#                 @test bx ≈ Form{D,R}((get_volumes(mfd, R)[n],))
#             elseif D == 2 && R == 1
#                 v1, v2 = VR([I0(1, 2), I0(1, 3), I0(2, 3)])[n]
#                 dx = get_coords(mfd)[v2] - get_coords(mfd)[v1]
#                 #TODO WRONG ORDER OF BASIS FUNCTIONS
#                 if !(bx ≈ Form{D,R}((i in (v1, v2)) * dx))
#                     @show D R n i
#                     for m in axes(get_simplices(mfd, P == Pr ? R : D - R), 2)
#                         v1, v2 = VR([I0(1, 2), I0(1, 3), I0(2, 3)])[m]
#                         dx = get_coords(mfd)[v2] - get_coords(mfd)[v1]
#                         @show Form{D,R}((i in (v1, v2)) * dx)
#                     end
#                 end
#                 v1, v2 = VR([I0(1, 2), I0(1, 3), I0(2, 3)])[n]
#                 dx = get_coords(mfd)[v2] - get_coords(mfd)[v1]
#                 @test bx ≈ Form{D,R}((i in (v1, v2)) * dx)
#             elseif D == 3 && R == 1
#                 v1, v2 = VR([I0(1, 2), I0(1, 3), I0(2, 3), I0(1, 4), I0(2, 4),
#                              I0(3, 4)])[n]
#                 dx = get_coords(mfd)[v2] - get_coords(mfd)[v1]
#                 @test bx ≈ Form{D,R}((i in (v1, v2)) * dx)
#             elseif D == 4 && R == 1
#                 v1, v2 = VR([I0(1, 2), I0(1, 3), I0(2, 3), I0(1, 4), I0(2, 4),
#                              I0(3, 4), I0(1, 5), I0(2, 5), I0(3, 5), I0(4, 5)])[n]
#                 dx = get_coords(mfd)[v2] - get_coords(mfd)[v1]
#                 @test bx ≈ Form{D,R}((i in (v1, v2)) * dx)
#             elseif D == 5 && R == 1
#                 v1, v2 = VR([I0(1, 2), I0(1, 3), I0(2, 3), I0(1, 4), I0(2, 4),
#                              I0(3, 4), I0(1, 5), I0(2, 5), I0(3, 5), I0(4, 5),
#                              I0(1, 6), I0(2, 6), I0(3, 6), I0(4, 6), I0(5, 6)])[n]
#                 dx = Form{D,R}(get_coords(mfd)[v2] - get_coords(mfd)[v1])
#                 @test bx ≈ (i in (v1, v2)) * dx
#             elseif D == 3 && R == 2
#                 v1, v2, v3 = VR([I0(1, 2, 3), I0(1, 2, 4), I0(1, 3, 4),
#                                  I0(2, 3, 4)])[n]
#                 dx = Form{D,1}(get_coords(mfd)[v2] - get_coords(mfd)[v1]) ∧
#                      Form{D,1}(get_coords(mfd)[v3] - get_coords(mfd)[v1]) / 2
#                 @test bx ≈ (i in (v1, v2, v3)) * dx
#             elseif D == 4 && R == 2
#                 v1, v2, v3 = VR([I0(1, 2, 3), I0(1, 2, 4), I0(1, 3, 4),
#                                  I0(2, 3, 4), I0(1, 2, 5), I0(1, 3, 5),
#                                  I0(2, 3, 5), I0(1, 4, 5), I0(2, 4, 5),
#                                  I0(3, 4, 5)])[n]
#                 dx = Form{D,1}(get_coords(mfd)[v2] - get_coords(mfd)[v1]) ∧
#                      Form{D,1}(get_coords(mfd)[v3] - get_coords(mfd)[v1]) / 2
#                 @test bx ≈ (i in (v1, v2, v3)) * dx
#             elseif D == 5 && R == 2
#                 v1, v2, v3 = VR([I0(1, 2, 3), I0(1, 2, 4), I0(1, 3, 4),
#                                  I0(2, 3, 4), I0(1, 2, 5), I0(1, 3, 5),
#                                  I0(2, 3, 5), I0(1, 4, 5), I0(2, 4, 5),
#                                  I0(3, 4, 5), I0(1, 2, 6), I0(1, 3, 6),
#                                  I0(2, 3, 6), I0(1, 4, 6), I0(2, 4, 6),
#                                  I0(3, 4, 6), I0(1, 5, 6), I0(2, 5, 6),
#                                  I0(3, 5, 6), I0(4, 5, 6)])[n]
#                 dx = Form{D,1}(get_coords(mfd)[v2] - get_coords(mfd)[v1]) ∧
#                      Form{D,1}(get_coords(mfd)[v3] - get_coords(mfd)[v1]) / 2
#                 @test bx ≈ (i in (v1, v2, v3)) * dx
#             elseif D == 4 && R == 3
#                 v1, v2, v3, v4 = VR([I0(1, 2, 3, 4), I0(1, 2, 3, 5),
#                                      I0(1, 2, 4, 5), I0(1, 3, 4, 5),
#                                      I0(2, 3, 4, 5)])[n]
#                 dx = Form{D,1}(get_coords(mfd)[v2] - get_coords(mfd)[v1]) ∧
#                      Form{D,1}(get_coords(mfd)[v3] - get_coords(mfd)[v1]) ∧
#                      Form{D,1}(get_coords(mfd)[v4] - get_coords(mfd)[v1]) / 6
#                 @test bx ≈ (i in (v1, v2, v3, v4)) * dx
#             end
#         end
#     end
# end

@testset "Evaluate a function D=$D P=$P R=$R" for D in 0:Dmax,
P in (Pr, Dl),
R in 0:D

    # TODO: all P
    P == Pr || continue

    S = Float64
    mfd = simplex_manifold(Val(D), S)

    T = S

    f0 = zero(Fun{D,P,R,D,S,T}, mfd)
    nvalues = nsimplices(mfd, P == Pr ? R : D - R)
    fi = Fun{D,P,R,D,S,T}(mfd, IDVector{R}([T(i) for i in 1:nvalues]))
    fx = Fun{D,P,R,D,S,SVector{D,T}}(mfd, get_coords(mfd, R))
    f1 = Fun{D,P,R,D,S,T}(mfd, IDVector{R}(rand(T, nvalues)))
    f2 = Fun{D,P,R,D,S,T}(mfd, IDVector{R}(rand(T, nvalues)))
    a = rand(T)
    b = rand(T)

    if R == 0
        for (i, x) in enumerate(get_coords(mfd, R))
            @test evaluate(fi, x) ≈ Form{D,R}((T(i),))
            @test evaluate(fx, x) ≈ Form{D,R}((x,))
        end
    end

    xs = [
        collect(get_coords(mfd, R))
        [random_point(Val(D), mfd) for i in 1:10]
    ]
    for x in xs
        @test evaluate(f0, x) == zero(Form{D,R,T})
        @test evaluate(f1 + f2, x) ≈ evaluate(f1, x) + evaluate(f2, x)
        @test evaluate(a * f1, x) ≈ a * evaluate(f1, x)
    end

    if D == 0
        # vertices: 1
        x1 = get_coords(mfd)[ID{0}(1)]
        fi1 = evaluate(fi, x1)
        if R == 0
            @test fi1 ≈ Form{D,0}((1,))
        else
            @assert false
        end
    elseif D == 1
        # vertices: 1 2
        # edges:    (1,2)
        x1 = get_coords(mfd)[ID{0}(1)]
        x2 = get_coords(mfd)[ID{0}(2)]
        fi1 = evaluate(fi, x1)
        fi2 = evaluate(fi, x2)
        if R == 0
            @test fi1 ≈ Form{D,0}((1,))
            @test fi2 ≈ Form{D,0}((2,))
        elseif R == 1
            @test fi1 ≈ Form{D,1}((1,))
            @test fi2 ≈ Form{D,1}((1,))
        else
            @assert false
        end
    elseif D == 2
        # vertices: 1 2 3
        # edges:    (1,2) (1,3) (2,3)
        # faces:    (1,2,3)
        x1 = get_coords(mfd)[ID{0}(1)]
        x2 = get_coords(mfd)[ID{0}(2)]
        x3 = get_coords(mfd)[ID{0}(3)]
        fi1 = evaluate(fi, x1)
        fi2 = evaluate(fi, x2)
        fi3 = evaluate(fi, x3)
        if R == 0
            b1 = Form{D,0}((S(1),))
            b2 = Form{D,0}((S(1),))
            b3 = Form{D,0}((S(1),))
            @test fi1 ≈ 1 * b1
            @test fi2 ≈ 2 * b2
            @test fi3 ≈ 3 * b3
        elseif R == 1
            b1 = Form{D,1}(x2 - x1)
            b2 = Form{D,1}(x3 - x1)
            b3 = Form{D,1}(x3 - x2)
            @test fi1 ≈ 1 * b1 + 2 * b2
            @test fi2 ≈ 1 * b1 + 3 * b3
            @test fi3 ≈ 2 * b2 + 3 * b3
        elseif R == 2
            b1 = Form{D,1}(x2 - x1) ∧ Form{D,1}(x3 - x1) / 2
            @test fi1 ≈ 1 * b1
            @test fi2 ≈ 1 * b1
            @test fi3 ≈ 1 * b1
        else
            @assert false
        end
    end

    if R < D
        dfi = deriv(fi)
        if D == 0
            @assert false
        elseif D == 1
            # vertices: 1 2
            # edges:    (1,2)
            x1 = get_coords(mfd)[ID{0}(1)]
            x2 = get_coords(mfd)[ID{0}(2)]
            fi1 = evaluate(fi, x1)
            fi2 = evaluate(fi, x2)
            dfi1 = evaluate(dfi, x1)
            dfi2 = evaluate(dfi, x2)
            if R == 0
                @test fi1 ≈ Form{D,0}((1,))
                @test fi2 ≈ Form{D,0}((2,))
                @test dfi1 ≈ (-1 + 2) * Form{D,1}(-x1 + x2)
                @test dfi2 ≈ (-1 + 2) * Form{D,1}(-x1 + x2)
            else
                @assert false
            end
        elseif D == 2
            # vertices: 1 2 3
            # edges:    (1,2) (1,3) (2,3)
            # faces:    (1,2,3)
            x1 = get_coords(mfd)[ID{0}(1)]
            x2 = get_coords(mfd)[ID{0}(2)]
            x3 = get_coords(mfd)[ID{0}(3)]
            fi1 = evaluate(fi, x1)
            fi2 = evaluate(fi, x2)
            fi3 = evaluate(fi, x3)
            dfi1 = evaluate(dfi, x1)
            dfi2 = evaluate(dfi, x2)
            dfi3 = evaluate(dfi, x3)
            if R == 0
                b1 = Form{D,0}((S(1),))
                b2 = Form{D,0}((S(1),))
                b3 = Form{D,0}((S(1),))
                @test fi1 ≈ 1 * b1
                @test fi2 ≈ 2 * b2
                @test fi3 ≈ 3 * b3
                db1 = Form{D,1}(x2 - x1)
                db2 = Form{D,1}(x3 - x1)
                db3 = Form{D,1}(x3 - x2)
                @test dfi1 ≈ (-1 + 2) * db1 + (-1 + 3) * db2
                @test dfi2 ≈ (-1 + 2) * db1 + (-2 + 3) * db3
                @test dfi3 ≈ (-1 + 3) * db2 + (-2 + 3) * db3
            elseif R == 1
                b1 = Form{D,1}(x2 - x1)
                b2 = Form{D,1}(x3 - x1)
                b3 = Form{D,1}(x3 - x2)
                @test fi1 ≈ 1 * b1 + 2 * Form{D,1}(x3 - x1)
                @test fi2 ≈ 1 * b1 + 3 * Form{D,1}(x3 - x2)
                @test fi3 ≈ 2 * b2 + 3 * Form{D,1}(x3 - x2)
                db1 = Form{D,1}(x2 - x1) ∧ Form{D,1}(x3 - x1) / 2
                @test dfi1 ≈ (1 - 2 + 3) * db1
                @test dfi2 ≈ (1 - 2 + 3) * db1
                @test dfi3 ≈ (1 - 2 + 3) * db1
            else
                @assert false
            end
        end
    end
end

@testset "Sample a function D=$D P=$P R=$R" for D in 0:Dmax,
P in (Pr, Dl),
R in 0:D

    # TODO: all P
    P == Pr || continue

    # TODO: all R
    R == 0 || continue

    S = Float64
    # mfd = simplex_manifold(Val(D), S)
    # for level in 1:2
    #     mfd = refined_manifold(mfd; optimize_mesh=false)
    # end
    mfd = large_hypercube_manifold(Val(D), S; nelts=ntuple(d -> 2, D),
                                   optimize_mesh=false)

    # nomenclature: fDRN where D=dim, R=rank, N=basis element
    # D = 0, R = 0
    c₀₀₁ = rand(S)
    f₀₀₁(x) = c₀₀₁ * Form{0,0,S}((1,)) # 1
    f₀₀(x) = f₀₀₁(x)
    # D = 1, R = 0
    c₁₀₁ = rand(S)
    f₁₀₁(x) = c₁₀₁ * Form{1,0,S}((1,)) # 1
    c₁₀₂ = rand(S)
    f₁₀₂(x) = c₁₀₂ * x[1] * Form{1,0,S}((1,)) # x
    f₁₀(x) = f₁₀₁(x) + f₁₀₂(x)
    # D = 1, R = 1
    c₁₁₁ = rand(S)
    f₁₁₁(x) = c₁₁₁ * Form{1,1,S}((1,)) # dx
    f₁₁(x) = f₁₁₁(x)
    # D = 2, R = 0
    c₂₀₁ = rand(S)
    f₂₀₁(x) = c₂₀₁ * Form{2,0,S}((1,)) # 1
    c₂₀₂ = rand(S)
    f₂₀₂(x) = c₂₀₂ * x[1] * Form{2,0,S}((1,)) # x
    c₂₀₃ = rand(S)
    f₂₀₃(x) = c₂₀₃ * x[2] * Form{2,0,S}((1,)) # y
    f₂₀(x) = f₂₀₁(x) + f₂₀₂(x) + f₂₀₃(x)
    # D = 2, R = 1
    c₂₁₁ = S(1) # rand(S)
    f₂₁₁(x) = c₂₁₁ * Form{2,1,S}((1, 0)) # dx
    c₂₁₂ = S(0) # rand(S)
    f₂₁₂(x) = c₂₁₂ * Form{2,1,S}((0, 1)) # dy
    c₂₁₃ = S(0) # rand(S)
    f₂₁₃(x) = c₂₁₃ * (-x[2] * Form{2,1,S}((1, 0)) + x[1] * Form{2,1,S}((0, 1))) # - y dx + x dy
    f₂₁(x) = f₂₁₁(x) + f₂₁₂(x) + f₂₁₃(x)
    # D = 2, R = 2
    c₂₂₁ = rand(S)
    f₂₂₁(x) = c₂₂₁ * Form{2,2,S}((1, 1)) # dx ∧ dy
    f₂₂(x) = f₂₂₁(x)
    # D = 3, R = 0
    c₃₀₁ = rand(S)
    f₃₀₁(x) = c₃₀₁ * Form{3,0,S}((1,)) # 1
    c₃₀₂ = rand(S)
    f₃₀₂(x) = c₃₀₂ * x[1] * Form{3,0,S}((1,)) # x
    c₃₀₃ = rand(S)
    f₃₀₃(x) = c₃₀₃ * x[2] * Form{3,0,S}((1,)) # y
    c₃₀₄ = rand(S)
    f₃₀₄(x) = c₃₀₄ * x[3] * Form{3,0,S}((1,)) # z
    f₃₀(x) = f₃₀₁(x) + f₃₀₂(x) + f₃₀₃(x) + f₃₀₄(x)
    # D = 4, R = 0
    c₄₀ₙ = rand(SVector{5,S})
    function f₄₀ₙ(x, n)
        local y = SVector{5,S}(1, x...)
        return c₄₀ₙ[n] * y[n] * Form{4,0,S}((1,))
    end
    f₄₀(x) = sum(f₄₀ₙ(x, n) for n in 1:5)
    # D = 5, R = 0
    c₅₀ₙ = rand(SVector{6,S})
    function f₅₀ₙ(x, n)
        local y = SVector{6,S}(1, x...)
        return c₅₀ₙ[n] * y[n] * Form{5,0,S}((1,))
    end
    f₅₀(x) = sum(f₅₀ₙ(x, n) for n in 1:5)

    function f(x)
        D == 0 && R == 0 && return f₀₀(x)
        D == 1 && R == 0 && return f₁₀(x)
        D == 1 && R == 1 && return f₁₁(x)
        D == 2 && R == 0 && return f₂₀(x)
        D == 2 && R == 1 && return f₂₁(x)
        D == 2 && R == 2 && return f₂₂(x)
        D == 3 && R == 0 && return f₃₀(x)
        D == 4 && R == 0 && return f₄₀(x)
        D == 5 && R == 0 && return f₅₀(x)
        @assert false
    end
    T = S
    f(zero(SVector{D,S}))::Form{D,R,T}
    f̃ = sample(Fun{D,P,R,D,S,T}, f, mfd)

    xs = [
        collect(get_coords(mfd, R))
        [random_point(Val(R), mfd) for i in 1:10]
    ]
    for x in xs
        fx = f(x)::Form{D,R,T}
        f̃x = evaluate(f̃, x)::Form{D,R,T}
        Ex = norm(f̃x - fx)
        @test 1 + Ex ≈ 1
    end

    #TODO if R < D
    #TODO     function df(x)
    #TODO         @assert m !== missing
    #TODO         local r = m
    #TODO         return r::Form{D,R + 1,T}
    #TODO     end
    #TODO
    #TODO     df̃ = sample(Fun{D,P,R + 1,D,S,T}, df, mfd)
    #TODO     df̃::Fun{D,P,R + 1,D,S,T}
    #TODO
    #TODO     d̃f̃ = deriv(f̃₁)
    #TODO     d̃f̃::Fun{D,P,R + 1,D,S,T}
    #TODO
    #TODO     @test norm(d̃f̃ - df̃, Inf) ≤ sqrt(eps(T))
    #TODO
    #TODO     # We cannot evaluate derivatives pointwise since there are not
    #TODO     # enough degrees of freedom to describe all form components in
    #TODO     # all directions. The discrete function has discontinuous
    #TODO     # normals across interfaces.
    #TODO end
    #TODO
    #TODO if R > 0
    #TODO     function δf(x)
    #TODO         @assert m′ !== missing
    #TODO         local r = bitsign(R) * bitsign(R * (D - R)) * m′
    #TODO         return r::Form{D,R - 1,T}
    #TODO     end
    #TODO
    #TODO     δf̃ = sample(Fun{D,P,R - 1,D,S,T}, δf, mfd)
    #TODO     δf̃::Fun{D,P,R - 1,D,S,T}
    #TODO
    #TODO     δ̃f̃ = coderiv(f̃)
    #TODO     δ̃f̃::Fun{D,P,R - 1,D,S,T}
    #TODO
    #TODO     if !(norm(δ̃f̃ - δf̃, Inf) ≤ sqrt(eps(T)))
    #TODO         @show mfd
    #TODO         @show D P R
    #TODO         @show get_coords(mfd, 0) get_coords(mfd, R)
    #TODO         @show f̃ δ̃f̃ δf̃
    #TODO     end
    #TODO     @test norm(δ̃f̃ - δf̃, Inf) ≤ sqrt(eps(T))
    #TODO end
end

@testset "Project a function D=$D P=$P R=$R" for D in 0:Dmax,
P in (Pr, Dl),
R in 0:D

    # TODO: all P
    P == Pr || continue

    S = Float64
    mfd = simplex_manifold(Val(D), S)
    # for level in 1:2
    #     mfd = refined_manifold(mfd; optimize_mesh=false)
    # end
    # mfd = large_hypercube_manifold(Val(D), S; nelts=ntuple(d -> 2, D),
    #                                optimize_mesh=false)

    # m = R + 1 > D ? missing : rand(Form{D,R + 1,S})
    # m′ = R - 1 < 0 ? missing : rand(Form{D,R - 1,S})
    # b = rand(Form{D,R,S})
    # function f(x)
    #     r = b
    #     if m !== missing
    #         r += m ⋅ Form{D,1,S}(x)
    #     end
    #     if m′ !== missing
    #         r += ⋆(⋆m′ ⋅ Form{D,1,S}(x))
    #     end
    #     return r::Form{D,R,T}
    # end

    # nomenclature: fDRN where D=dim, R=rank, N=basis element
    # D = 0, R = 0
    c₀₀₁ = rand(S)
    f₀₀₁(x) = c₀₀₁ * Form{0,0,S}((1,)) # 1
    f₀₀(x) = f₀₀₁(x)
    # D = 1, R = 0
    c₁₀₁ = rand(S)
    f₁₀₁(x) = c₁₀₁ * Form{1,0,S}((1,)) # 1
    c₁₀₂ = rand(S)
    f₁₀₂(x) = c₁₀₂ * x[1] * Form{1,0,S}((1,)) # x
    f₁₀(x) = f₁₀₁(x) + f₁₀₂(x)
    # D = 1, R = 1
    c₁₁₁ = rand(S)
    f₁₁₁(x) = c₁₁₁ * Form{1,1,S}((1,)) # dx
    f₁₁(x) = f₁₁₁(x)
    # D = 2, R = 0
    c₂₀₁ = rand(S)
    f₂₀₁(x) = c₂₀₁ * Form{2,0,S}((1,)) # 1
    c₂₀₂ = rand(S)
    f₂₀₂(x) = c₂₀₂ * x[1] * Form{2,0,S}((1,)) # x
    c₂₀₃ = rand(S)
    f₂₀₃(x) = c₂₀₃ * x[2] * Form{2,0,S}((1,)) # y
    f₂₀(x) = f₂₀₁(x) + f₂₀₂(x) + f₂₀₃(x)
    # D = 2, R = 1
    c₂₁₁ = S(1) # rand(S)
    f₂₁₁(x) = c₂₁₁ * Form{2,1,S}((1, 0)) # dx
    c₂₁₂ = S(0) # rand(S)
    f₂₁₂(x) = c₂₁₂ * Form{2,1,S}((0, 1)) # dy
    c₂₁₃ = S(0) # rand(S)
    f₂₁₃(x) = c₂₁₃ * (-x[2] * Form{2,1,S}((1, 0)) + x[1] * Form{2,1,S}((0, 1))) # - y dx + x dy
    f₂₁(x) = f₂₁₁(x) + f₂₁₂(x) + f₂₁₃(x)
    # D = 2, R = 2
    c₂₂₁ = rand(S)
    f₂₂₁(x) = c₂₂₁ * Form{2,2,S}((1, 1)) # dx ∧ dy
    f₂₂(x) = f₂₂₁(x)

    function f(x)
        D == 0 && R == 0 && return f₀₀(x)
        D == 1 && R == 0 && return f₁₀(x)
        D == 1 && R == 1 && return f₁₁(x)
        D == 2 && R == 0 && return f₂₀(x)
        D == 2 && R == 1 && return f₂₁(x)
        D == 2 && R == 2 && return f₂₂(x)
        D == 3 && R == 0 && return f₃₀(x)
        D == 4 && R == 0 && return f₄₀(x)
        D == 5 && R == 0 && return f₅₀(x)
        @assert false
    end
    T = S
    f(zero(SVector{D,S}))::Form{D,R,T}
    f̃ = project(Fun{D,P,R,D,S,T}, f, mfd)

    xs = [
        collect(get_coords(mfd, R))
        [random_point(Val(R), mfd) for i in 1:10]
    ]
    for x in xs
        fx = f(x)::Form{D,R,T}
        f̃x = evaluate(f̃, x)::Form{D,R,T}
        Ex = norm(f̃x - fx)
        if !(1 + Ex ≈ 1)
            @show D R get_coords(mfd, R) f̃ x fx f̃x Ex
        end
        @test 1 + Ex ≈ 1
    end

    # TODO: check derivatives etc; translate code from sampling
end


text:
#=
day16:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-16
=#

using AoC, Test

const UINT_BITS = sizeof(UInt) * 8
const UINT_ONES = typemax(UInt)

mutable struct BitStream{T}
    values::Vector{T}
    byteOffset::UInt
    bitOffset::UInt
end
BitStream(values::Vector{T}) where T = BitStream{T}(values,one(UInt),zero(UInt))

function read!(stream::BitStream{T}, n)::UInt where T
    elmask = typemax(unsigned(T))
    elbits = count_ones(elmask)
    stream
    result::UInt = 0
    while n > 0
        if (stream.byteOffset > length(stream.values))
            throw(BoundsError(stream.values, stream.byteOffset))
        end
        remainingBitsInByte = elbits - stream.bitOffset
        bitsToRead = min(remainingBitsInByte, n)

        mask = elmask >> (elbits - bitsToRead)
        bits = (stream.values[stream.byteOffset] >> (remainingBitsInByte - bitsToRead)) & mask

        result = (result << bitsToRead) | bits

        stream.bitOffset += bitsToRead
        n -= bitsToRead
        if stream.bitOffset >= elbits
            stream.bitOffset = 0
            stream.byteOffset += 1
        end
    end

    return result
end
testStream = BitStream([0b10101110_00100000_00000000_00000000])
@test read!(testStream, 3) == 0b101
@test read!(testStream, 3) == 0b011
@test read!(testStream, 5) == 0b10001
testStream = BitStream(UInt8[0b10101110, 0b00100000, 0b00000000, 0b00000000])
@test read!(testStream, 3) == 0b101
@test read!(testStream, 3) == 0b011
@test read!(testStream, 5) == 0b10001

function eof(stream::BitStream)
    return stream.byteOffset > length(stream.values)
end

function lineToStream(line)
    charsPerWord = sizeof(UInt) * 2
    length(line), length(line) % charsPerWord, line
    if length(line) % charsPerWord != 0
        line = rpad(line, length(line) + charsPerWord - (length(line) % charsPerWord), '0')
    end
    return BitStream([line[I:I+charsPerWord-1] for I in 1:charsPerWord:length(line)] |> ll -> parse.(UInt, ll; base=16))
end
@test lineToStream("ABCDEF12").values == UInt[0xabcdef12_00000000]

abstract type AbstractPacket end

struct ValuePacket <: AbstractPacket
    version::UInt
    value::UInt64
end

struct OperatorPacket <: AbstractPacket
    version::UInt
    operator
    packets::Vector{AbstractPacket}
end

function readLiteralValue!(stream)
    bitsRead = 0
    result = 0
    finished = false
    while !finished
        nextValue = read!(stream, 5)
        bitsRead += 5
        result = (result << 4) | (nextValue & 0b1111)
        if nextValue & 0b10000 == 0
            finished = true
        end
    end
    return (result, bitsRead)
end

function readPackets!(stream; maxLength=typemax(Int), maxPackets=typemax(Int))
    bitsRead = 0
    packetsRead = 0
    packets = AbstractPacket[]
    while bitsRead < maxLength && packetsRead < maxPackets && !eof(stream)
        (packet, bitsReadThisTime) = readPacket!(stream)
        push!(packets, packet)
        bitsRead += bitsReadThisTime
        packetsRead += 1
    end
    return (packets, bitsRead)
end

sumOperator(values) = +(0, values...)
productOperator(values) = *(1, values...)
minimumOperator(values) = min(values...)
maximumOperator(values) = max(values...)
greaterThanOperator(values) = values[1] > values[2]
lessThanOperator(values) = values[1] < values[2]
equalToOperator(values) = values[1] == values[2]

function readPacket!(stream)
    version = read!(stream, 3)
    typeId = read!(stream, 3)
    bitsRead = 6

    if typeId == 4
        (value, bitsReadThisTime) = readLiteralValue!(stream)
        bitsRead += bitsReadThisTime
        return (ValuePacket(version, value), bitsRead)
    else
        # operator packet
        lengthTypeId = read!(stream, 1)
        bitsRead += 1
        if lengthTypeId == 0
            contentLength = read!(stream, 15)
            bitsRead += 15
            (packets, bitsReadThisTime) = readPackets!(stream; maxLength = contentLength)
        else
            packetCount = read!(stream, 11)
            bitsRead += 11
            (packets, bitsReadThisTime) = readPackets!(stream; maxPackets = packetCount)
        end
        bitsRead += bitsReadThisTime

        if typeId == 0
            operator = sumOperator
        elseif typeId == 1
            operator = productOperator
        elseif typeId == 2
            operator = minimumOperator
        elseif typeId == 3
            operator = maximumOperator
        elseif typeId == 5
            operator = greaterThanOperator
        elseif typeId == 6
            operator = lessThanOperator
        elseif typeId == 7
            operator = equalToOperator
        end

        return (OperatorPacket(version, operator, packets), bitsRead)
    end
end

totalValue(packet::ValuePacket) = packet.value
totalValue(packet::OperatorPacket) = packet.operator(totalValue.(packet.packets))

totalVersions(packet::ValuePacket) = packet.version
totalVersions(packet::OperatorPacket) = packet.version +(totalVersions.(packet.packets)...)

function part1(line)
    (packet, bitsRead) = readPacket!(lineToStream(line))

    return signed(totalVersions(packet))
end

function part2(line)
    (packet, bitsRead) = readPacket!(lineToStream(line))

    return signed(totalValue(packet))
end

@test part1("D2FE28") == 6
@test part1("8A004A801A8002F478") == 16
@test part1("620080001611562C8802118E34") == 12
@test part1("C0015000016115A2E0802F182340") == 23
@test part1("A0016C880162017C3686B18A3D4780") == 31

@test part2("C200B40A82") == 3
@test part2("04005AC33890") == 54
@test part2("880086C3E88112") == 7
@test part2("CE00C43D881120") == 9
@test part2("D8005AC2A8F0") == 1
@test part2("F600BC2D8F") == 0
@test part2("9C005AC2F8F0") == 0
@test part2("9C0141080250320F1802104A08") == 1

@show lines(16)[1] |> l -> @time part1(l)
@show lines(16)[1] |> l -> @time part2(l)


text:
using GPUBenchmarks, BenchmarkTools, FileIO

benchmark_files = [
    "blackscholes",
    "PDE",
    "poincare",
    "juliaset",
    "mapreduce"
]

GPUBenchmarks.new_run()

for device in GPUBenchmarks.devices()
    name, typ = GPUBenchmarks.init(device)
    @show name typ
end

device = first(GPUBenchmarks.devices())
nothing
file = "juliaset"
device = :cuarrays

for file in benchmark_files
    # for device in GPUBenchmarks.devices()
        @run_julia (JULIA_NUM_THREADS = 8, "-O3", file, device) begin
            using GPUBenchmarks, BenchmarkTools, FileIO
            bench_mod = include(GPUBenchmarks.dir("benchmark", file * ".jl"))
            println("Benchmarking $file $device")
            result = bench_mod.execute(device)
            println("Benchmarking done for $device")
            GPUBenchmarks.append_data!(result)
        end
    # end
end

GPUBenchmarks.last_time_stamp()


text:
#=
Test function

Properties
----------

- There are 9 steady states:
    - Stable
        - (-2.39,2.35)    (2.3,-2.34) (-2.18,-2.12)   (2.06,2.12)
    - Unstable
        - (-2.27,-0.23)   (2.17,0.22) (-0.12,-2.23)   (0.33,2.22) (0.1,0.01)

Reference
---------
1. Source code from https://static-content.springer.com/esm/art%3A10.1038%2Fs41598-017-15889-2/MediaObjects/41598_2017_15889_MOESM2_ESM.rar
=#

module SampleDE

"""
Reference
---------
- `ODE.m`
"""
function ODE!(du,u,p,t)
    du[1]= -1 + 9 *u[1] - 2*u[1]^3 + 9 *u[2] - 2*u[2]^3;
    du[2]=  1 - 11*u[1] + 2*u[1]^3 + 11*u[2] - 2*u[2]^3;
end

function ODE(u,p,t)
    du = deepcopy(u)
    return ODE!(du,u,p,t)
end

u0 = [1.0, 2.0]
p = [2.0, 3.0]

"""
Reference
---------
- `exact_potential.m`
"""
function exact_potential(points, RefPointNum)
    pot = zeros(length(points), 1)
    for i in 1:length(points)
        x = points[i][1]; y = points[i][2];
        pot[i] = -5*(x^2+y^2)+(x^4+y^4)/2+x*y+x;
    end
    pot = pot .- pot[RefPointNum]
    return pot
end

"""
Reference
---------
- `Readme.pdf`
"""
ODE_features = Dict(
    "StableStates" => [ (-2.39,2.35),(2.3,-2.34), (-2.18,-2.12),   (2.06,2.12)  ],
    "UnstableStates" => [(-2.27,-0.23) ,  (2.17,0.22), (-0.12,-2.23),   (0.33,2.22), (0.1,0.01)],
    "SteadyStates" => [(-2.39,2.35),(2.3,-2.34), (-2.18,-2.12),   (2.06,2.12), (-2.27,-0.23) ,  (2.17,0.22), (-0.12,-2.23),   (0.33,2.22) ,(0.1,0.01)]
)

function gradient(dt,x)

    x = reshape(x, 2, 3)

    F=[ [(dt*(1/dt - 6*x[1, 1]^2 + 9)*(9*x[1, 1] + 9*x[2, 1] + (x[1, 1] - x[1, 2])/dt - 2*x[1, 1]^3 - 2*x[2, 1]^3 - 1))/2 + (dt*(6*x[1, 1]^2 - 11)*(11*x[2, 1] - 11*x[1, 1] + (x[2, 1] - x[2, 2])/dt + 2*x[1, 1]^3 - 2*x[2, 1]^3 + 1))/2,x[1, 1]^3 - (9*x[2, 1])/2 - (x[1, 1] - x[1, 2])/(2*dt) - (9*x[1, 1])/2 + x[2, 1]^3 + (dt*(1/dt - 6*x[1, 2]^2 + 9)*(9*x[1, 2] + 9*x[2, 2] + (x[1, 2] - x[1, 3])/dt - 2*x[1, 2]^3 - 2*x[2, 2]^3 - 1))/2 + (dt*(6*x[1, 2]^2 - 11)*(11*x[2, 2] - 11*x[1, 2] + (x[2, 2] - x[2, 3])/dt + 2*x[1, 2]^3 - 2*x[2, 2]^3 + 1))/2 + 1/2,x[1, 2]^3 - (9*x[2, 2])/2 - (x[1, 2] - x[1, 3])/(2*dt) - (9*x[1, 2])/2 + x[2, 2]^3 + 1/2 ],[
    (dt*(1/dt - 6*x[2, 1]^2 + 11)*(11*x[2, 1] - 11*x[1, 1] + (x[2, 1] - x[2, 2])/dt + 2*x[1, 1]^3 - 2*x[2, 1]^3 + 1))/2 - (dt*(6*x[2, 1]^2 - 9)*(9*x[1, 1] + 9*x[2, 1] + (x[1, 1] - x[1, 2])/dt - 2*x[1, 1]^3 - 2*x[2, 1]^3 - 1))/2,(11*x[1, 1])/2 - (11*x[2, 1])/2 - (x[2, 1] - x[2, 2])/(2*dt) - x[1, 1]^3 + x[2, 1]^3 + (dt*(1/dt - 6*x[2, 2]^2 + 11)*(11*x[2, 2] - 11*x[1, 2] + (x[2, 2] - x[2, 3])/dt + 2*x[1, 2]^3 - 2*x[2, 2]^3 + 1))/2 - (dt*(6*x[2, 2]^2 - 9)*(9*x[1, 2] + 9*x[2, 2] + (x[1, 2] - x[1, 3])/dt - 2*x[1, 2]^3 - 2*x[2, 2]^3 - 1))/2 - 1/2,(11*x[1, 2])/2 - (11*x[2, 2])/2 - (x[2, 2] - x[2, 3])/(2*dt) - x[1, 2]^3 + x[2, 2]^3 - 1/2]]
    return hcat(F...)'
end

end


text:
function rdict(args...)
    DefaultOrderedDict{String,Any,typeof(rdict)}(rdict, OrderedDict{String,Any}(args...))
end

const POD_TEMPLATE =
    rdict("apiVersion" => "v1",
          "kind" => "Pod",
          "metadata" => rdict(),
          "spec" => rdict("restartPolicy" => "Never",
                          "containers" => []))

struct KubeError <: Exception
    msg::String
end

KubeError(io::IO) = KubeError(String(take!(io)))

Base.showerror(io::IO, e::KubeError) = print(io, "KubeError: ", e.msg)

"""
    get_pod(name::AbstractString) -> AbstractDict

Retrieve details about the specified pod as a JSON-compatible dictionary. If their is no
pod with the given `name` then a `$KubeError` will be raised.
"""
function get_pod(name::AbstractString)
    kubectl_cmd = `$(kubectl()) get pod/$name -o json`
    err = IOBuffer()
    out = read(pipeline(ignorestatus(kubectl_cmd), stderr=err), String)

    err.size > 0 && throw(KubeError(err))
    return JSON.parse(out; dicttype=OrderedDict)
end

"""
    create_pod(manifest::AbstractDict) -> String

Create a pod based upon the JSON-compatible `manifest`. Returns the name of the pod created.
"""
function create_pod(manifest::AbstractDict)
    # As `kubectl create` can create any resource we'll restrict this function to only
    # creating "Pod" resources.
    kind = manifest["kind"]
    if kind != "Pod"
        throw(ArgumentError("Manifest expected to be of kind \"Pod\" and not \"$kind\""))
    end

    kubectl_cmd = `$(kubectl()) create -f -`
    out = IOBuffer()
    err = IOBuffer()
    open(pipeline(ignorestatus(kubectl_cmd), stdout=out, stderr=err), "w") do p
        write(p.in, JSON.json(manifest))
    end

    err.size > 0 && throw(KubeError(err))

    # Extract the pod name from the output. Needed when using "generateName".
    out_str = String(take!(out))
    m = match(r"pod/(?<name>.*?) created", out_str)
    if m !== nothing
        return m[:name]
    else
        error("Unable to determine the pod name from: \"$out_str\"")
    end
end

"""
    label_pod(name::AbstractString, label::Pair) -> Nothing

Create or replace a metadata label for a given pod `name`. Requires the "patch" permission.
"""
function label_pod(name::AbstractString, label::Pair)
    kubectl_cmd = `$(kubectl()) label --overwrite pod/$name $(join(label, '='))`
    err = IOBuffer()
    run(pipeline(ignorestatus(kubectl_cmd), stdout=devnull, stderr=err))

    err.size > 0 && throw(KubeError(err))
    return nothing
end

"""
    delete_pod(name::AbstractString) -> Nothing

Delete the pod with the given `name`.
"""
function delete_pod(name::AbstractString; wait::Bool=true)
    kubectl_cmd = `$(kubectl()) delete pod/$name --wait=$wait`
    err = IOBuffer()
    run(pipeline(ignorestatus(kubectl_cmd), stdout=devnull, stderr=err))

    err.size > 0 && throw(KubeError(err))
    return nothing
end

"""
    wait_for_running_pod(name::AbstractString; timeout::Real) -> AbstractDict

Wait for the pod with the given `name` to reach the "Running" phase. If the phase is reached
before the `timeout` then the pod details will be returned, otherwise a `TimeoutException`
will be raised.
"""
function wait_for_running_pod(name::AbstractString; timeout::Real)
    pod = nothing

    # https://kubernetes.io/docs/concepts/workloads/pods/pod-lifecycle/#pod-phase
    result = timedwait(timeout; pollint=1) do
        pod = @mock get_pod(name)
        pod["status"]["phase"] != "Pending"
    end

    if result === :ok
        phase = pod !== nothing ? pod["status"]["phase"] : nothing

        if phase == "Running"
            return pod
        else
            error("Worker pod $name expected to in phase \"Running\" and instead is \"$phase\"")
        end
    else
        msg = "timed out after waiting for worker $name to start for $timeout seconds, " *
            "with status:\n" * JSON.json(pod["status"], 4)
        throw(TimeoutException(msg))
    end
end

"""
    exec_pod(name::AbstractString, cmd::Cmd) -> Nothing

Execute `cmd` on the `name`d pod. If the command fails a `KubeError` will be raised.
"""
function exec_pod(name::AbstractString, cmd::Cmd)
    # When an `exec` failure occurs there is some useful information in stdout and stderr
    kubectl_cmd = `$(kubectl()) exec pod/$name -- $cmd`
    err = IOBuffer()
    p = run(pipeline(ignorestatus(kubectl_cmd), stdout=err, stderr=err))

    !success(p) && throw(KubeError(err))
    return nothing
end

"""
    pod_status(name::AbstractString) -> Pair
    pod_status(resource::AbstractDict) -> Pair

Determine the current pod state and, if applicable, the associated reason.
"""
function pod_status end

function pod_status(resource::AbstractDict)
    kind = resource["kind"]
    if kind != "Pod"
        throw(ArgumentError("Manifest expected to be of kind \"Pod\" and not \"$kind\""))
    end

    container_statuses = resource["status"]["containerStatuses"]

    if length(container_statuses) != 1
        throw(ArgumentError("Currently on pods with a single container are supported"))
    end

    container_status = first(container_statuses)
    states = keys(container_status["state"])

    if length(states) != 1
        json = JSON.json("state" => container_status["state"])
        throw(ArgumentError("Currently only a single state is supported:\n$json"))
    end

    state = first(states)
    return state => get(container_status["state"][state], "reason", nothing)
end

pod_status(name::AbstractString) = pod_status(get_pod(name))

"""
    worker_pod_spec(pod=POD_TEMPLATE; kwargs...)

Generate a pod specification for a Julia worker inside a container.
"""
function worker_pod_spec(pod::AbstractDict=POD_TEMPLATE; kwargs...)
    return worker_pod_spec!(deepcopy(pod); kwargs...)
end

function worker_pod_spec!(pod::AbstractDict;
                          manager_name::AbstractString,
                          image::AbstractString,
                          cmd::Cmd,
                          cpu=DEFAULT_WORKER_CPU,
                          memory=DEFAULT_WORKER_MEMORY,
                          service_account_name=nothing)
    pod["metadata"]["generateName"] = "$(manager_name)-worker-"

    # Set a label with the manager name to support easy termination of all workers
    pod["metadata"]["labels"]["manager"] = manager_name

    worker_container =
        rdict("name" => "worker",
              "image" => image,
              "command" => collect(cmd),
              "resources" => rdict("requests" => rdict("cpu" => cpu,
                                                       "memory" => memory),
                                   "limits"   => rdict("cpu" => cpu,
                                                       "memory" => memory)))

    push!(pod["spec"]["containers"], worker_container)

    if service_account_name !== nothing
        pod["spec"]["serviceAccountName"] = service_account_name
    end

    return pod
end

"""
    isk8s() -> Bool

Predicate for testing if the current process is running within a Kubernetes (K8s) pod.
"""
function isk8s()
    in_kubepod = false
    @mock(isfile("/proc/self/cgroup")) || return in_kubepod
    @mock open("/proc/self/cgroup") do fp
        while !eof(fp)
            line = chomp(readline(fp))
            path_name = split(line, ':')[3]
            if startswith(path_name, "/kubepods/")
                in_kubepod = true
                break
            end
        end
    end
    return in_kubepod
end


text:
module V
const NAMES = [
  "B"
  "M"
  "BM"
  "C"
  "CB"
  "CM"
  "Cp"
  "CpB"
  "CpM"
  "CBM"
  "CpBM"
  #(* "TAK1 module*)
  "TAK1"
  "TAK1p"
  "TAK1C"
  "TAK1pC"
  #(* "IKK module*)
  "IKK"
  "IKKC"
  "IKKp"
  "IKKpC"
  "IKKppC"
  "IKKpp"
  "IKKi"
  #(* NFkB module*) c:cytoplasm n:nucleus
  "NFKBc"
  "NFKBn"
  "IKBac"
  "IKBbc"
  "IKBec"
  "IKBan"
  "IKBbn"
  "IKBen"
  "NFKBIKBac"
  "NFKBIKBbc"
  "NFKBIKBec"
  "NFKBIKBan"
  "NFKBIKBbn"
  "NFKBIKBen"
  "IKKIKBac"
  "IKKIKBbc"
  "IKKIKBec"
  "IKKNFKBIKBac"
  "IKKNFKBIKBbc"
  "IKKNFKBIKBec"
  "mRNAac"
  "mRNAbc"
  "mRNAec"
  "mRNAa20c"
  "A20c"
]

for (idx,name) in enumerate(NAMES)
  eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)

end  # module

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim


text:
export gibbssamplescopes!

"""
    gibbssamplescopes!(spn::InfiniteSumNode, x::AbstractMatrix{T})

Apply Gibbs sampling of y_d assignments for each dimension at each components sub-SPN.
"""
function gibbssamplescopes!(spn::InfiniteSumNode,
                            x::AbstractMatrix{T}) where {T<:Real}
    for child in spn.children
        gibbssamplescopes!(child, x)
    end
end

"""
    gibbssamplescopes!(spn::AbstractRegionGraphNode, x::AbstractMatrix{T})

Gibbs sample all y_d assignments for each dimension and update all scopes in-place.
"""
function gibbssamplescopes!(spn::AbstractRegionGraphNode,
                            x::AbstractMatrix{T}) where {T<:Real}
    (N,D) = size(x)
    nodes = topoligicalorder(spn);
    partitionnodes = filter(n -> n isa PartitionGraphNode, nodes)
    dims = collect(1:D)

    shuffle!(dims)

    out = AxisArray(Vector{Vector}(undef, length(nodes)), map(id, nodes))
    y = AxisArray(zeros(Int, length(partitionnodes)), map(n -> id(n), partitionnodes))

    for d in dims

        # Remove dimension from scopes.
        @inbounds for n in nodes
            setscope!(n, d, false)
        end

        # Bottom-up path.
        for n in nodes
            logmllh!(n, x, out, y, d)
        end

        # Top-down assignment.
        _gibbsassignscope!(spn, y, d)

    end
end

function _gibbsassignscope!(n::RegionGraphNode,
                           y::AxisArray, d::Int)
    setscope!(n, d, true)
    for child in children(n)
        _gibbsassignscope!(child, y, d)
    end
end

function _gibbsassignscope!(n::PartitionGraphNode,
                           y::AxisArray, d::Int)
    setscope!(n, d, true)
    _gibbsassignscope!(children(n)[y[id(n)]], y, d)
end

function _gibbsassignscope!(n::FactorizedDistributionGraphNode,
                           y::AxisArray, d::Int)
    setscope!(n, d, true)
end

function _gibbsassignscope!(n::FactorizedMixtureGraphNode,
                           y::AxisArray, d::Int)
    setscope!(n, d, true)
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using Thermodynamics.TemperatureProfiles
using Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates
using StaticArrays
using Test

using CLIMAParameters
using CLIMAParameters.Atmos.SubgridScale: C_smag
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

### Citation
# [Schar2002](@cite)

# ## [Initial Conditions]
function init_schar!(problem, bl, state, aux, localgeo, t)
    ## Problem float-type
    FT = eltype(state)

    (x, y, z) = localgeo.coord

    ## Unpack constant parameters
    param_set = parameter_set(bl)
    R_gas::FT = R_d(param_set)
    c_p::FT = cp_d(param_set)
    c_v::FT = cv_d(param_set)
    p0::FT = MSLP(param_set)
    _grav::FT = grav(param_set)
    γ::FT = c_p / c_v

    c::FT = c_v / R_gas
    c2::FT = R_gas / c_p

    Tiso::FT = 250.0
    θ0::FT = Tiso

    ## Calculate the Brunt-Vaisaila frequency for an isothermal field
    ## Hydrostatic background state
    Brunt::FT = _grav / sqrt(c_p * Tiso)
    Brunt2::FT = Brunt * Brunt
    g2::FT = _grav * _grav

    π_exner::FT = exp(-_grav * z / (c_p * Tiso))
    θ::FT = θ0 * exp(Brunt2 * z / _grav)
    ρ::FT = p0 / (R_gas * θ) * (π_exner)^c

    ## Compute perturbed thermodynamic state:
    T = θ * π_exner
    e_int = internal_energy(param_set, T)
    ts = PhaseDry(param_set, e_int, ρ)

    ## Initial velocity
    z₁::FT = 4000
    z₂::FT = 5000
    u₀::FT = 10
    zscale = (z - z₁) / (z₂ - z₁)
    if z₂ <= z
        u = FT(1)
    elseif z₁ <= z < z₂
        u = (sinpi(zscale / 2))^2
    elseif z <= z₁
        u = FT(0)
    end
    u *= u₀

    ## Initial scalar anomaly profile
    ## Equivalent to a nondiffusive tracer
    Ax::FT = 25000
    Az::FT = 3000
    x₀::FT = 25000
    z₀::FT = 9000
    r = ((x - x₀) / Ax)^2 + ((z - z₀) / Az)^2
    if r <= 1
        χ = (cospi(r / 2))^2
    else
        χ = 0
    end

    ## State (prognostic) variable assignment
    e_kin = FT(1 / 2) * u^2                               # kinetic energy
    e_pot = gravitational_potential(bl.orientation, aux)# potential energy
    ρe_tot = ρ * total_energy(e_kin, e_pot, ts)         # total energy

    state.ρ = ρ
    state.ρu = SVector{3, FT}(ρ * u, 0, 0)
    state.energy.ρe = ρe_tot
    state.tracers.ρχ = ρ * SVector{1, FT}(χ)
end

# Define a `setmax` method
function setmax(f, xmax, ymax, zmax)
    function setmaxima(xin, yin, zin)
        return f(xin, yin, zin; xmax = xmax, ymax = ymax, zmax = zmax)
    end
    return setmaxima
end

function warp_schar(
    xin,
    yin,
    zin;
    xmax = 150000.0,
    ymax = 5000.0,
    zmax = 25000.0,
)
    FT = eltype(xin)
    a::FT = 25000 ## Half-width parameter [m]
    r = sqrt(xin^2 + yin^2)
    h₀::FT = 3000 ## Peak height [m]
    λ::FT = 8000 ## Wavelength
    h_star =
        abs(xin - xmax / 2) <= a ? h₀ * (cospi((xin - xmax / 2) / 2a))^2 : FT(0)
    h = h_star * (cospi((xin - xmax / 2) / λ))^2
    x, y, z = xin, yin, zin + h * (zmax - zin) / zmax
    return x, y, z
end

function config_schar(FT, N, resolution, xmax, ymax, zmax)
    u_relaxation = SVector(FT(10), FT(0), FT(0))

    ## Wave damping coefficient (1/s)
    sponge_ampz = FT(0.5)

    ## Vertical level where the absorbing layer starts
    z_s = FT(20000.0)

    ## Pass the sponge parameters to the sponge calculator
    rayleigh_sponge =
        RayleighSponge{FT}(zmax, z_s, sponge_ampz, u_relaxation, 2)

    ## Setup the source terms for this problem:
    source = (Gravity(), rayleigh_sponge)

    ## Define the reference state:
    T_virt = FT(250)
    temp_profile_ref = IsothermalProfile(param_set, T_virt)
    ref_state = HydrostaticState(temp_profile_ref)
    nothing # hide

    # Define a warping function to build an analytic topography:

    _C_smag = FT(0.21)
    _δχ = SVector{1, FT}(0)

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = Vreman(_C_smag),
        moisture = DryModel(),
        tracers = NTracers{1, FT}(_δχ),
    )

    model = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        init_state_prognostic = init_schar!,
        source = source,
    )

    config = ClimateMachine.AtmosLESConfiguration(
        "ScharScalarAdvection",  # Problem title [String]
        N,                       # Polynomial order [Int]
        resolution,              # (Δx, Δy, Δz) effective resolution [m]
        xmax,                    # Domain maximum size [m]
        ymax,                    # Domain maximum size [m]
        zmax,                    # Domain maximum size [m]
        param_set,               # Parameter set.
        init_schar!,             # Function specifying initial condition
        model = model,           # Model type
        meshwarp = setmax(warp_schar, xmax, ymax, zmax),
    )

    return config
end

# Define a `main` method (entry point)
function main()

    FT = Float64

    ## Define the polynomial order and effective grid spacings:
    N = 4

    ## Define the domain size and spatial resolution
    xmax = FT(150000)
    ymax = FT(2500)
    zmax = FT(25000)
    Δx = FT(500)
    Δy = FT(500)
    Δz = FT(500)
    resolution = (Δx, Δy, Δz)

    t0 = FT(0)
    timeend = FT(10000)

    ## Define the max Courant for the time time integrator (ode_solver).
    ## The default value is 1.7 for LSRK144:
    CFL = FT(1.5)

    ## Assign configurations so they can be passed to the `invoke!` function
    driver_config = config_schar(FT, N, resolution, xmax, ymax, zmax)

    ## Define the time integrator:
    ## We chose an explicit single-rate LSRK144 for this problem
    ode_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_solver_type = ode_solver_type,
        init_on_cpu = true,
        Courant_number = CFL,
    )

    ## State Conservation Callback
    # State variable
    Q = solver_config.Q
    # Volume geometry information
    vgeo = driver_config.grid.vgeo
    M = vgeo[:, Grids._M, :]
    # Unpack prognostic vars
    ρχ₀ = Q[:, 6, :]
    # DG variable sums
    Σρχ₀ = sum(ρχ₀ .* M)
    cb_check_tracer = GenericCallbacks.EveryXSimulationSteps(1000) do
        Q = solver_config.Q
        δρχ = (sum(Q[:, 6, :] .* M) .- Σρχ₀) ./ Σρχ₀
        @show (abs(δρχ))
        nothing
    end

    ## Set up the spectral filter to remove the solutions spurious modes
    ## Define the order of the exponential filter: use 32 or 64 for this problem.
    ## The larger the value, the less dissipation you get:
    filterorder = 64
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            state_auxiliary = solver_config.dg.state_auxiliary,
        )
        nothing
    end
    ## End exponential filter

    ## Invoke solver (calls `solve!` function for time-integrator),
    ## pass the driver, solver and diagnostic config information.
    result = ClimateMachine.invoke!(
        solver_config;
        user_callbacks = (cbfilter, cb_check_tracer),
        check_euclidean_distance = true,
    )

    ## Check that the solution norm is reasonable.
    @test isapprox(result, FT(1); atol = 1.5e-4)
end

# Call `main`
main()


text:
using MeasureTheory
using Base.Iterators
using Statistics

using Random
rng = Random.Xoshiro(3)

x = Chain(Normal()) do xj Normal(μ=xj) end
xobs = rand(rng, x)
y = For(xobs) do xj Poisson(logλ=xj) end
yobs = rand(rng, y)
xv = take(xobs, 10) |> collect
yv = take(yobs, 10) |> collect

take(xobs.parent, 10) |> collect
take(yobs.parent, 10) |> collect

exp.(xv)
yv

# using Plots

# plt = scatter(normcdf.(xv, 1, yv), label=false)

# for j in 1:10
#     xobs = rand(rng, x)
#     yobs = rand(rng, y)
#     xv = take(xobs, 100) |> collect;
#     yv = take(yobs, 100) |> collect;
#     plt = scatter!(plt, normcdf.(xv, 1, yv), label=false)
# end
# plt

using Soss

m = @model begin
    x ~ Chain(Normal()) do xj Normal(μ=xj) end
    y ~ For(xobs) do xj Poisson(logλ=xj) end
end

truth = rand(rng, m())

xobs = take(truth.x, 10) |> collect
yobs = take(truth.y, 10) |> collect

logdensity(m(), (x=xobs, y=yobs))


text:
#!/usr/bin/env julia

using Base: thispatch, thisminor, nextpatch, nextminor
using Base.Pkg.Reqs: Reqs, Requirement
using Base.Pkg.Types
using Base.LinAlg: checksquare
using Base.Random: UUID
using SHA

## Computing UUID5 values from (namespace, key) pairs ##

function uuid5(namespace::UUID, key::String)
    data = [reinterpret(UInt8, [namespace.value]); Vector{UInt8}(key)]
    u = reinterpret(UInt128, sha1(data)[1:16])[1]
    u &= 0xffffffffffff0fff3fffffffffffffff
    u |= 0x00000000000050008000000000000000
    return UUID(u)
end
uuid5(namespace::UUID, key::AbstractString) = uuid5(namespace, String(key))

const uuid_dns = UUID(0x6ba7b810_9dad_11d1_80b4_00c04fd430c8)
const uuid_julia = uuid5(uuid_dns, "julialang.org")
const uuid_package = uuid5(uuid_julia, "package")
const uuid_registry = uuid5(uuid_julia, "registry")

## Loading data into various data structures ##

struct Require
    versions::VersionInterval
    systems::Vector{Symbol}
end

struct Version
    sha1::String
    requires::Dict{String,Require}
end

struct Package
    uuid::UUID
    url::String
    versions::Dict{VersionNumber,Version}
end

Require(versions::VersionInterval) = Require(versions, Symbol[])
Version(sha1::String) = Version(sha1, Dict{String,Require}())

function load_requires(path::String)
    requires = Dict{String,Require}()
    requires["julia"] = Require(VersionInterval())
    isfile(path) || return requires
    for r in filter!(r->r isa Requirement, Reqs.read(path))
        # @assert length(r.versions.intervals) == 1
        new = haskey(requires, r.package)
        versions, systems = r.versions.intervals[1], r.system
        if haskey(requires, r.package)
            versions = versions ∩ requires[r.package].versions
            systems  = systems  ∪ requires[r.package].systems
        end
        requires[r.package] = Require(versions, systems)
    end
    return requires
end

function load_versions(dir::String)
    versions = Dict{VersionNumber,Version}()
    isdir(dir) || return versions
    for ver in readdir(dir)
        path = joinpath(dir, ver)
        sha1 = joinpath(path, "sha1")
        isfile(sha1) || continue
        requires = load_requires(joinpath(path, "requires"))
        versions[VersionNumber(ver)] = Version(readchomp(sha1), requires)
    end
    return versions
end

function load_packages(dir::String)
    pkgs = Dict{String,Package}()
    for pkg in readdir(dir)
        path = joinpath(dir, pkg)
        url = joinpath(path, "url")
        versions = joinpath(path, "versions")
        isfile(url) || continue
        pkgs[pkg] = Package(
            uuid5(uuid_package, pkg),
            readchomp(url),
            load_versions(versions),
        )
    end
    pkgs["julia"] = Package(
        uuid5(uuid_package, "julia"),
        "https://github.com/JuliaLang/julia.git",
        Dict(
            v"0.1.1"  => Version("bddb70374fdbb2256782398c13430fa4d8b55d0b"),
            v"0.1.2"  => Version("1187040d027210cc466f0b6a6d54118fd692cf2d"),
            v"0.2.0"  => Version("68e4f269049b3a77eee1e185e619c3513034b3a6"),
            v"0.2.1"  => Version("7c5dea64f55afb7d9a7224a337d3c0bd5a707c50"),
            v"0.3.0"  => Version("768187890c2709bf2ff06818f40e1cdc79bd44b0"),
            v"0.3.1"  => Version("c03f413bbdb46c00033f4eaad402995cfe3b7be5"),
            v"0.3.2"  => Version("8227746b95146c2921f83d2ae5f37ecd146592d8"),
            v"0.3.3"  => Version("b24213b893ace6ca601722f6cdaac15f2b034b7b"),
            v"0.3.4"  => Version("33920266908e950936ef3f8503810a2d47333741"),
            v"0.3.5"  => Version("a05f87b79ad62beb033817fdfdefa270c9557aaf"),
            v"0.3.6"  => Version("0c24dca65c031820b91721139f0291068086955c"),
            v"0.3.7"  => Version("cb9bcae93a32b42cec02585c387396ff11836aed"),
            v"0.3.8"  => Version("79599ada444fcc63cb08aed64b4ff6a415bb4d29"),
            v"0.3.9"  => Version("31efe690bea1fd8f4e44692e205fb72d34f50ad3"),
            v"0.3.10" => Version("c8ceeefcc1dc25953a644622a895a3adcbc80dad"),
            v"0.3.11" => Version("483dbf5279efe29f94cf35d42de0ab0a673fb34e"),
            v"0.3.12" => Version("80aa77986ec84067f4696439f10fc77c93703bc7"),
            v"0.4.0"  => Version("0ff703b40afddf9b705bd6a06d3a59cb4c089ea5"),
            v"0.4.1"  => Version("cbe1bee3a8f27ba4f349556857bf4615ee5fa68f"),
            v"0.4.2"  => Version("bb73f3489d837e3339fce2c1aab283d3b2e97a4c"),
            v"0.4.3"  => Version("a2f713dea5ac6320d8dcf2835ac4a37ea751af05"),
            v"0.4.4"  => Version("ae683afcc4cc824dcf68251569a6cf8760e4fb6f"),
            v"0.4.5"  => Version("2ac304dfba75fad148d4070ef4f8a2e400c305bb"),
            v"0.4.6"  => Version("2e358ce975029ec97aba5994c17d4a2169c3b085"),
            v"0.4.7"  => Version("ae26b25d43317d7dd3ca05f60b70677aab9c0e08"),
            v"0.5.0"  => Version("3c9d75391c72d7c32eea75ff187ce77b2d5effc8"),
            v"0.5.1"  => Version("6445c82d0060dbe82b88436f0f4371a4ee64d918"),
            v"0.5.2"  => Version("f4c6c9d4bbbd9587d84494e314f692c15ff1f9c0"),
        ),
    )
    return pkgs
end

@eval julia_versions() = $([VersionNumber(0,m) for m=1:5])
julia_versions(f::Function) = filter(f, julia_versions())
julia_versions(vi::VersionInterval) = julia_versions(v->v in vi)

macro clean(ex) :(x = $(esc(ex)); $(esc(:clean)) &= x; x) end

function prune!(pkgs::Associative{String,Package})
    # remove unsatisfiable versions
    while true
        clean = true
        filter!(pkgs) do pkg, p
            filter!(p.versions) do ver, v
                @clean ver == thispatch(ver) > v"0.0.0" &&
                all(v.requires) do kv
                    req, r = kv
                    haskey(pkgs, req) &&
                    any(w->w in r.versions, keys(pkgs[req].versions))
                end
            end
            @clean !isempty(p.versions)
        end
        clean && break
    end
    return pkgs
end

## Load package data ##

const pkgs = load_packages(Pkg.dir("METADATA"))
prune!(pkgs)


text:
module JSONStat

using DataStructures, DataFrames

export readjsondataset, readjsonbundle, writejsondataset, labelframe

optkeys = DataStructures.OrderedDict(
		"dataset" => [
			"version"; "class"; "role"; "label"; "updated"; "source";
			"extension"; "href";"link"; "note"; "error"
			],
		"dimid" => [
			"label"
			],
		"category" => [
			"label"; "child"; "coordinates"; "unit"; "note"
			]
	)

function dimframe(id, dimjsonstat)
	dimid = dimjsonstat[id]
	dimcat = dimid["category"]
	dimdict = DataStructures.OrderedDict()
	dimdict["category"] = DataStructures.OrderedDict()
	if haskey(dimcat, "index")
		dimdata = dimcat["index"]
	else
		dimdata = dimcat["label"]
	end
	for key in keys(dimid)
		if key in optkeys["dimid"]
			dimdict[key] = dimid[key]
		end
	end
	for key in keys(dimcat)
		if key in optkeys["category"]
			dimdict["category"][key] = dimcat[key]
		end
	end
	dimframe = DataFrame()
	if typeof(dimdata)<:DataStructures.OrderedDict
		dimframe[!, Symbol(id)] = collect(keys(dimdata))
	else
		dimframe[!, Symbol(id)] = dimdata
	end
	dimdict["dimframe"] = dimframe
	id => dimdict
end

function readjsondataset(datasetjsonstat)
	if haskey(datasetjsonstat, "version")
		idarr = datasetjsonstat["id"]
	else
		idarr = datasetjsonstat["dimension"]["id"]
	end
	dimension = DataStructures.OrderedDict(map((x)->
		dimframe(x, datasetjsonstat["dimension"]), idarr))
	dimkeys = collect(keys(dimension))
	datasetframe = dimension[dimkeys[1]]["dimframe"]
	for i in 2:size(dimkeys)[1]
		datasetframe = crossjoin(datasetframe, dimension[dimkeys[i]]["dimframe"])
	end
	datasetframe[!, :value] = datasetjsonstat["value"]
	datasetdict = DataStructures.OrderedDict()
	for key in keys(datasetjsonstat)
		if key in optkeys["dataset"]
			datasetdict[key] = datasetjsonstat[key]
		end
	end
	datasetdict["dimension"] = dimension
	datasetdict["datasetframe"] = datasetframe
	datasetdict
end

function labelframe(id, datasetdict)
	labeldict = datasetdict["dimension"][id]["category"]["label"]
	labelframe = DataFrame()
	labelframe[!, Symbol(id)] = collect(keys(labeldict))
	labelframe[!, Symbol("$(id)_label")] = collect(values(labeldict))
	labelframe
end

function readjsonbundle(bundlejsonstat)
	bundledict = DataStructures.OrderedDict()
	for key in keys(bundlejsonstat)
		if haskey(bundlejsonstat[key], "value")
			bundledict[key] = readjsondataset(bundlejsonstat[key])
		end
	end
	bundledict
end

function writedimid(id, datasetjsonstat)
	datasetjsonstat["id"] = vcat(datasetjsonstat["id"], id)
	dimid = datasetjsonstat["dimension"][id]
	dimcat = dimid["category"]
	dimframe = dimid["dimframe"]
	dimsize = size(dimframe)[1]
	datasetjsonstat["size"] = vcat(datasetjsonstat["size"], dimsize)
	dimzip = zip(dimframe[Symbol(id)], collect(0:dimsize-1))
	dimcat["index"] = DataStructures.OrderedDict(dimzip)
	delete!(dimid, "dimframe")
end

function writejsondataset(datasetdict)
	datasetjsonstat = deepcopy(datasetdict)
	if !(haskey(datasetjsonstat, "version"))
		datasetjsonstat["version"] = "2.0"
	end
	datasetjsonstat["id"] = []
	datasetjsonstat["size"] = []
	datasetjsonstat["value"] = datasetjsonstat["datasetframe"][!, :value]
	delete!(datasetjsonstat, "datasetframe")
	idarr = collect(keys(datasetjsonstat["dimension"]))
	foreach((x)->writedimid(x, datasetjsonstat), idarr)
	datasetjsonstat
end

end # module


text:
using HashDicts: HashDict

d = HashDict([("A", 1), ("B", 2)])
d2 = HashDict("A"=>1, "B"=>2)

for (k, v) in [("hi", 2), ("", 0)]
    global d
    d[k] = v
    @assert d[k] == v
end


text:
function NDA_GAIL_JS(;π,
                       S,
                       𝒟_demo,
                       𝒟_nda,
                       Vc::ContinuousNetwork, # value network for estimating cost
                       γ::Float32,
                       λ_gae::Float32 = 0.95f0,
                       αr::Float32=0.5f0,
                       normalize_demo::Bool=true,
                       D::ContinuousNetwork,
                       Dnda::ContinuousNetwork,
                       solver=LagrangePPO,
                       gan_loss::Crux.GANLoss=GAN_BCELoss(),
                       d_opt::NamedTuple=(;),
                       d_opt_nda::NamedTuple=(;),
                       log::NamedTuple=(;),
                       kwargs...)

    d_opt = TrainingParams(;loss = Crux.GAIL_D_loss(gan_loss), name="discriminator_", d_opt...)
    d_opt_nda = TrainingParams(;loss = Crux.GAIL_D_loss(gan_loss), name="nda_discriminator_", d_opt_nda...)
    if normalize_demo
        𝒟_demo = normalize!(deepcopy(𝒟_demo), S, action_space(π))
        𝒟_nda = normalize!(deepcopy(𝒟_nda), S, action_space(π))
    end
    𝒟_demo = 𝒟_demo |> device(π)
    𝒟_nda = 𝒟_nda |> device(π)

    function GAIL_callback(𝒟; info=Dict(), 𝒮)
        batch_train!(D, d_opt, (;), 𝒟_demo, deepcopy(𝒟), info=info)
        batch_train!(Dnda, d_opt_nda, (;), 𝒟_nda, deepcopy(𝒟), info=info)

        # Set the reward
        D_out = value(D, 𝒟[:a], 𝒟[:s]) # This is swapped because a->x and s->y and the convention for GANs is D(x,y)
        r = αr * logσ.(D_out) .- (1f0 - αr) * logcompσ.(D_out)
        ignore() do
            minval, maxval = extrema(D_out)
            info["disc_reward"] = mean(r)
        end
        𝒟[:r] .= r

        # Set the cost
        D_out_nda = value(Dnda, 𝒟[:a], 𝒟[:s])
        r_nda = αr * logσ.(D_out_nda) .- (1f0 - αr) * logcompσ.(D_out_nda)
        c = max.(0, r_nda .- r)
        # c = max.(0, σ.(D_out_nda) .- σ.(D_out))
        ignore() do
            info["disc_nda_cost"] = sum(c) / sum(𝒟[:episode_end])
        end
        𝒟[:cost] .= c

        eps = episodes(𝒟)
        for eprange in eps
            ep = eprange[1]:eprange[2]
            fill_gae!(𝒟, ep, 𝒮.agent.π, λ_gae, γ)
            fill_returns!(𝒟, ep, γ)

            fill_gae!(𝒟, ep, Vc, λ_gae, γ, source=:cost, target=:cost_advantage)
            fill_returns!(𝒟, ep, γ, source=:cost, target=:cost_return)
        end
        𝒟[:advantage] .= whiten(𝒟[:advantage])
        𝒟[:cost_advantage] .= whiten(𝒟[:cost_advantage])

    end
    solver(;π=π, S=S, post_batch_callback=GAIL_callback, Vc=Vc, λ_gae=λ_gae, log=(dir="log/ndagail", period=500, log...), kwargs...)
end


text:
using SimulatedNeuralMoments
using Test

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

max_stars_repo_name: mcreel/SimulatedNeuralMoments

max_issues_repo_name: mcreel/SimulatedNeuralMoments


text:
using BLG_Fraunhofers
using Test

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


text:
module jlmie
    using SpecialFunctions

    include("general_functions.jl")
    include("mie_main.jl")
    export jlmie_mx, jlmie_abcd, jlmie_pt, jlmie_S12
    include("efficiency.jl")
    export jlmie_Qsca, jlmie_Qsca_a_n, jlmie_Qsca_b_n
    include("farfield.jl")
    export jlmie_Isff, jlmie_Isff_EM_n
end

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
ALL_PARAMETERS[:concept] = [(:no_overlap, 1)]

# NOTE - Please use clear_results=false (default) for the real experiments
# icn_benchmark(; clear_results=true)

icn_benchmark()