text:
__precompile__()

module Baseline
using Flux
abstract type AbstractBaseline end

struct Linear{DT<:AbstractFloat} <: AbstractBaseline
   nfeat::Int
   T::Int
   numT::Int
   reg_coeff::DT
   coeff::Vector{DT} # states x N

   # workspace
   featnd::Array{DT,3}
   A::Matrix{DT}
   bline::Vector{DT}
   function Linear{DT}(reg::AbstractFloat, ns::Integer, T::Integer, numT::Integer) where DT<:AbstractFloat
      nfeat = 2*ns+3+1
      featnd=zeros(DT, nfeat, T, numT)
      A = zeros(DT, nfeat, nfeat)
      bline = zeros(DT, T*numT)
      return new(nfeat, T, numT, reg, zeros(DT, nfeat),
                 featnd, A, bline)
   end
end

struct Quadratic{DT<:AbstractFloat} <: AbstractBaseline
   nfeat::Int
   T::Int
   numT::Int
   reg_coeff::DT
   coeff::Vector{DT} # states x N

   # workspace
   featnd::Array{DT,3}
   A::Matrix{DT}
   bline::Vector{DT}
   function Quadratic{DT}(reg::AbstractFloat, ns::Integer, T::Integer, numT::Integer) where DT<:AbstractFloat
      nfeat = ns+convert(Int,ns*(ns+1)/2)+ns+ns+4+1
      featnd=zeros(DT, nfeat, T, numT)
      gettimefeatures!(featnd)

      A = zeros(DT, nfeat, nfeat)
      bline = zeros(DT, T*numT)
      return new(nfeat, T, numT, reg, zeros(DT, nfeat),
                 featnd, A, bline)
   end
end

struct NN{DT<:AbstractFloat} <: AbstractBaseline
   nhidden::Int
   nfeat::Int
   T::Int
   numT::Int

   loss::Function
   opt::Function

   #workspace
   featnd::Array{DT,2}
   m::Any
   #w::Array{Any}
   bline::Vector{DT}
   function NN{DT}(ns::Integer, T::Integer, numT::Integer, nhidden::Integer) where DT<:AbstractFloat
      nfeat = 1*ns+3+1
      featnd=zeros(DT, nfeat, T*numT)
      m = Chain(Dense(nfeat, nhidden, tanh),
                Dense(nhidden, nhidden, tanh),
                Dense(nhidden, 1)) |> gpu
      bline = zeros(DT, T*numT)
      return new(nhidden, nfeat, T, numT,
                 (a, b) -> sum((m(a) - b).^2)/length(b),
                 ADAM(params(m)), #, 0.0001),
                 featnd, m, bline)
   end
end

######################################################################## common

function gettimefeatures!{DT<:AbstractFloat}(feat::Array{DT,3})
   T    = size(feat, 2)
   numT = size(feat, 3)
   al = convert(Array{DT}, linspace(1, T, T)') / DT(T)
   for t=1:numT
      feat[end-4,:,t] = al
      feat[end-3,:,t] = al.^2
      feat[end-2,:,t] = al.^3
      feat[end-1,:,t] = al.^4
      feat[end,:,t]   = DT(1.0)
   end
end

function getlinearfeatures(obs::Matrix{Float64})
   T = size(obs,2)
   # TODO add clipping
   o = max.(obs, -10.0)
   o = min.(o, 10.0)
   al = linspace(1, T, T)' / T
   feat = vcat(o, o.^2, al, al.^2, al.^3, ones(1, T))
end

function getquadfeatures!{DT<:AbstractFloat}(feat::Array{DT,3},
                                             obs::Matrix{DT},
                                             t::Integer)
   n = size(obs, 1)
   T = size(feat, 2)

   feat[1:n,:,t] = obs[:,(t-1)*T+1:t*T] #implicite conversion between 32, 64
   feat[1:n,:,t] = max.(feat[1:n,:,t], DT(-10.0))
   feat[1:n,:,t] = min.(feat[1:n,:,t], DT(10.0))
   feat[1:n,:,t] /= DT(10.0)

   qfsize = convert(Int, n*(n+1)/2)
   qf = view(feat, (n+1):(n+qfsize) , :, t) #Matrix{Float64}(qfsize, T)
   k = 1
   @inbounds for i=1:n
      for j=i:n
         for x=1:T
            qf[k,x] = feat[i,x,t].*feat[j,x,t]
         end
         k += 1
      end
   end

   i0 = n+qfsize
   i1 = i0+n
   @inbounds for j = 1:n
      for i=1:T
         feat[i0+j,i,t] = sin(pi*feat[j,i,t])
         feat[i1+j,i,t] = cos(pi*feat[j,i,t])
      end
   end
end

function getNNfeatures(obs::Matrix{Float64})
   T = size(obs,2)
   # TODO add clipping
   o = max.(obs, -10.0)
   o = min.(o, 10.0)
   # al = linspace(1, T, T)' / 100.0
   al = linspace(1, T, T)' / 1000.0
   # feat = vcat(o, o.^2, al, al.^2, al.^3, ones(1, T))
   feat = vcat(o, al, al.^2, al.^3, ones(1, T))
end

function prefit!{DT<:AbstractFloat}(b::AbstractBaseline, returns::Vector{DT})
   featmat = reshape(b.featnd, b.nfeat, b.numT*b.T)
   #b.A[:] = (featmat*featmat')
   A_mul_Bt!(b.A, featmat, featmat)

   return featmat*returns
end

# returns: T*numT
function fit!{DT<:AbstractFloat}(b::AbstractBaseline, returns::Vector{DT})

   target = prefit!(b, returns) # sets A matrix, returns target

   for i=1:b.nfeat # set diag
      b.A[i,i] += b.reg_coeff
   end
   #b.coeff[:] = b.A\(target)
   if isposdef(b.A)
      A_ldiv_B!(b.coeff, cholfact!(b.A), target) # needs A to be factorized
   else
      A_ldiv_B!(b.coeff, lufact!(convert(Array{Float64}, b.A)), target) # needs A to be factorized
   end
end

function fit!{DT<:AbstractFloat}(b::NN, returns::Vector{DT})
   nsamples = length(returns)
   const nbatch = 50

   X = b.featnd |> gpu
   Y = returns |> gpu

   rndidx = Flux.chunk(shuffle(1:nsamples), nbatch)

   datas = [ (X[:,i], Y[i]') for i in rndidx ]

   opt = ADAM(params(b.m))
   evalcb() = @show(b.loss(b.featnd, returns'))

   BLAS.set_num_threads(Threads.nthreads())
   for i=1:3
      Flux.train!(b.loss, datas, b.opt)#, cb = Flux.throttle(evalcb, 1))
   end

end

######################################################################## linear

function predict!{T<:AbstractFloat}(b::Linear, obs::Matrix{T})
   nthread = min.(b.numT, Threads.nthreads()) # cant have more threads than data
   Threads.@threads for tid=1:nthread
      #for tid=1:nthread
      thread_range = Distributed.splitrange(b.numT, nthread)[tid]

      for t=thread_range
         b.featnd[:,:,t] = getlinearfeatures(obs[:,(t-1)*b.T+1:t*b.T])
      end
   end
   b.bline[:] = reshape(b.featnd, b.nfeat, b.numT*b.T)'*b.coeff
end

##################################################################### quadratic
# TODO this is hacky and this whole module is a disgrace
function predict{DT<:AbstractFloat}(b::Quadratic, obs::Matrix{DT})
   nsamples = size(obs, 2)
   featnd=zeros(DT, b.nfeat, 1, nsamples)
   gettimefeatures!(featnd)
   for t=1:nsamples
      getquadfeatures!(featnd, obs, t)
   end
   values = reshape(featnd, b.nfeat, nsamples)'*b.coeff
   return values
end

function predict!{DT<:AbstractFloat}(b::Quadratic, obs::Matrix{DT})
   nthread = min.(b.numT, Threads.nthreads()) # cant have more threads than data
   Threads.@threads for tid=1:nthread
      #for tid=1:nthread
      thread_range = Distributed.splitrange(b.numT, nthread)[tid]

      for t=thread_range
         getquadfeatures!(b.featnd, obs, t)
      end
   end
   b.bline[:] = reshape(b.featnd, b.nfeat, b.numT*b.T)'*b.coeff
end

##################################################################### Neural Network

function predict!{T<:AbstractFloat}(b::NN, obs::Matrix{T})
   b.featnd[:,:] = getNNfeatures(obs)
   b.bline .= b.m(b.featnd).data[1,:]
end

end # module


text:
# We `import` methods from Nodes so we can extend them.
# The import must occur *before* the extensions/methods defined
# in splash and game scenes.
# import .Ranger.Nodes:
#     transition, get_replacement, visit

using .Ranger

module RangerGame
using ..Ranger
using ..Ranger.Rendering
using ..Ranger.Events
using ..Ranger.Geometry
using ..Ranger.Nodes
using ..Ranger.Nodes.Scenes
using ..Ranger.Nodes.Custom
using ..Ranger.Animation
using ..Ranger.Nodes.Filters

orange = Rendering.Orange()
white = Rendering.White()
darkgray = Rendering.DarkGray()
lightgray = Rendering.LightGray()
yellow = Rendering.Yellow()
blue = Rendering.Blue()
lightblue = Rendering.LightNavyBlue()
red = Rendering.Red()
lime = Rendering.Lime()
olive = Rendering.Olive()
lightpurple = Rendering.LightPurple()
peach = Rendering.Peach()

include("orbit_system_node.jl")
include("scene_boot.jl")
include("splash_scene.jl")

include("ship_node.jl")
include("game_layer.jl")

include("game_scene.jl")

export go

function build(world::Ranger.World)
    println("Building: ", world.title)

    game = GameScene(world, "GameScene")
    build(game, world)

    splash = SplashScene(world, "SplashScene", game)
    build(splash, world)
    splash.transitioning.pause_for = 0.25 * 1000.0

    boot = SceneBoot(world, "SceneBoot", splash)

    Ranger.Engine.push(boot)

    Nodes.print_tree(game)

    println("Built");
end

function go()
    world = Ranger.Engine.initialize("Template 0", build)

    if world ≠ nothing
        Ranger.Engine.run(world);
    end
end

end # Module ----------------------------------------------------------

ext: jl

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
""" Categories of C-sets and attributed C-sets.
"""
module CSets
export ACSetTransformation, CSetTransformation, SubACSet, SubCSet,
  components, force, is_natural, homomorphism, homomorphisms, is_homomorphic,
  isomorphism, isomorphisms, is_isomorphic,
  generate_json_acset, parse_json_acset, read_json_acset, write_json_acset

using Base.Meta: quot
using AutoHashEquals
using JSON
using Reexport
using StaticArrays: SVector

@reexport using ...CSetDataStructures
using ...GAT, ..FreeDiagrams, ..Limits, ..Subobjects, ..Sets, ..FinSets
import ..Limits: limit, colimit, universal
import ..Subobjects: Subobject, SubobjectBiHeytingAlgebra,
  implies, ⟹, subtract, \, negate, ¬, non, ~
import ..FinSets: FinSet, FinFunction, FinDomFunction, force, as_predicate
using ...Theories: Category, CatDesc, AttrDesc, ob, hom, attr, adom, acodom
import ...Theories: dom, codom, compose, ⋅, id,
  meet, ∧, join, ∨, top, ⊤, bottom, ⊥

# FinSets interop
#################

""" Create `FinSet` for part of attributed C-set.
"""
FinSet(X::ACSet, type::Symbol) = FinSet{Int,Int}(nparts(X, type))

""" Create `FinFunction` for part or subpart of attributed C-set.

The subpart must be of kind `Ob`. For indexed subparts, the index is included.
"""
FinFunction(X::ACSet, name::Symbol) = fin_function(X, Val{name})

@generated function fin_function(X::ACSet{CD,AD,Ts,Idxed},
    ::Type{Val{name}}) where {CD,AD,Ts,Idxed,name}
  if name ∈ ob(CD)
    quote
      FinFunction(identity, FinSet(X, $(QuoteNode(name))))
    end
  elseif name ∈ hom(CD)
    quote
      FinFunction(subpart(X, $(QuoteNode(name))),
                  FinSet(X, $(QuoteNode(codom(CD, name)))),
                  index=$(name ∈ Idxed ? :(X.indices.$name) : false))
    end
  else
    throw(ArgumentError("$(repr(name)) not in $(ob(CD)) or $(hom(CD))"))
  end
end

""" Create `FinDomFunction` for part or subpart of attributed C-set.

The codomain is always of type `TypeSet`, regardless of whether the subpart is
of kind `Ob` or `Data`. For indexed subparts, the index is included.
"""
FinDomFunction(X::ACSet, name::Symbol) = fin_dom_function(X, Val{name})

@generated function fin_dom_function(X::ACSet{CD,AD,Ts,Idxed},
    ::Type{Val{name}}) where {CD,AD,Ts,Idxed,name}
  if name ∈ ob(CD)
    quote
      n = nparts(X, $(QuoteNode(name)))
      FinDomFunction(1:n, FinSet(n), TypeSet{Int}())
    end
  elseif name ∈ hom(CD) || name ∈ attr(AD)
    quote
      FinDomFunction(subpart(X, $(QuoteNode(name))),
                     index=$(name ∈ Idxed ? :(X.indices.$name) : false))
    end
  else
    throw(ArgumentError(
      "$(repr(name)) not in $(ob(CD)), $(hom(CD)), or $(attr(AD))"))
  end
end

# C-set transformations
#######################

""" Transformation between attributed C-sets.

A homomorphism of C-sets is a natural transformation: a transformation between
functors C → Set satisfying the naturality axiom for all morphisms in C. This
struct records the data of a transformation; it does not enforce naturality, but
see [`is_natural`](@ref).

A C-set transformation has a component for every object in C. When C-sets have
attributes, the data types are assumed to be fixed. Thus, the naturality axiom
for data attributes is a commutative triangle, rather than a commutative square.
"""
@auto_hash_equals struct ACSetTransformation{CD <: CatDesc, AD <: AttrDesc{CD},
    Comp <: NamedTuple, Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}}
  components::Comp
  dom::Dom
  codom::Codom

  function ACSetTransformation{CD,AD}(components::NamedTuple, X::Dom, Y::Codom) where
      {Ob, CD <: CatDesc{Ob}, AD <: AttrDesc{CD},
       Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}}
    @assert keys(components) ⊆ Ob
    coerced_components = NamedTuple{Ob}(
      coerce_component(ob, get(components, ob) do; Int[] end, X, Y)
      for ob in Ob)
    new{CD,AD,typeof(coerced_components),Dom,Codom}(coerced_components, X, Y)
  end
end

function coerce_component(ob::Symbol, f::FinFunction{Int,Int}, X, Y)
  length(dom(f)) == nparts(X,ob) || error("Domain error in component $ob")
  length(codom(f)) == nparts(Y,ob) || error("Codomain error in component $ob")
  return f
end
function coerce_component(ob::Symbol, f, X, Y)::FinFunction{Int,Int}
  FinFunction(f, nparts(X,ob), nparts(Y,ob))
end

ACSetTransformation(components, X::Dom, Y::Codom) where
    {CD, AD, Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}} =
  ACSetTransformation{CD,AD}(components, X, Y)
ACSetTransformation(X::Dom, Y::Codom; components...) where
    {CD, AD, Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}} =
  ACSetTransformation{CD,AD}((; components...), X, Y)

const CSetTransformation{CD, Comp,
                         Dom <: AbstractCSet{CD}, Codom <: AbstractCSet{CD}} =
  ACSetTransformation{CD,AttrDesc{CD,(),(),(),()},Comp,Dom,Codom}

CSetTransformation(components, X::AbstractCSet, Y::AbstractCSet) =
  ACSetTransformation(components, X, Y)
CSetTransformation(X::AbstractCSet, Y::AbstractCSet; components...) =
  ACSetTransformation(X, Y; components...)

components(α::ACSetTransformation) = α.components
Base.getindex(α::ACSetTransformation, ob) = α.components[ob]

""" Is the transformation between C-sets a natural transformation?

Uses the fact that to check whether a transformation is natural, it suffices to
check the naturality equation on a generating set of morphisms.
"""
function is_natural(α::ACSetTransformation{CD,AD}) where {CD,AD}
  X, Y = dom(α), codom(α)
  for (f, c, d) in zip(hom(CD), dom(CD), codom(CD))
    Xf, Yf, α_c, α_d = subpart(X,f), subpart(Y,f), α[c], α[d]
    all(Yf[α_c(i)] == α_d(Xf[i]) for i in eachindex(Xf)) || return false
  end
  for (f, c) in zip(attr(AD), adom(AD))
    Xf, Yf, α_c = subpart(X,f), subpart(Y,f), α[c]
    all(Yf[α_c(i)] == Xf[i] for i in eachindex(Xf)) || return false
  end
  return true
end

map_components(f, α::ACSetTransformation) =
  ACSetTransformation(map(f, components(α)), dom(α), codom(α))
force(α::ACSetTransformation) = map_components(force, α)

# Finding C-set transformations
###############################

""" Find a homomorphism between two attributed ``C``-sets.

Returns `nothing` if no homomorphism exists. For many categories ``C``, the
``C``-set homomorphism problem is NP-complete and thus this procedure generally
runs in exponential time. It works best when the domain object is small.

This procedure uses the classic backtracking search algorithm for a
combinatorial constraint satisfaction problem (CSP). As is well known, the
homomorphism problem for relational databases is reducible to CSP. Since the
C-set homomorphism problem is "the same" as the database homomorphism problem
(insofar as attributed C-sets are "the same" as relational databases), it is
also reducible to CSP. Backtracking search for CSP is described in many computer
science textbooks, such as (Russell & Norvig 2010, *Artificial Intelligence*,
Third Ed., Chapter 6: Constraint satisfaction problems, esp. Algorithm 6.5). In
our implementation, the search tree is ordered using the popular heuristic of
"minimum remaining values" (MRV), also known as "most constrained variable."

To restrict to *monomorphisms*, or homomorphisms whose components are all
injective functions, set the keyword argument `monic=true`. To restrict only
certain components to be injective or bijective, use `monic=[...]` or
`iso=[...]`. For example, setting `monic=[:V]` for a graph homomorphism ensures
that the vertex map is injective but imposes no constraints on the edge map.

To restrict the homomorphism to a given partial assignment, set the keyword
argument `initial`. For example, to fix the first source vertex to the third
target vertex in a graph homomorphism, set `initial=(V=Dict(1 => 3),)`.

See also: [`homomorphisms`](@ref), [`isomorphism`](@ref).
"""
function homomorphism(X::AbstractACSet, Y::AbstractACSet; kw...)
  result = nothing
  homomorphisms(X, Y; kw...) do α
    result = α; return true
  end
  result
end

""" Find all homomorphisms between two attributed ``C``-sets.

This function is at least as expensive as [`homomorphism`](@ref) and when no
homomorphisms exist, it is exactly as expensive.
"""
function homomorphisms(X::AbstractACSet{CD,AD}, Y::AbstractACSet{CD,AD};
                       kw...) where {CD,AD}
  results = ACSetTransformation{CD,AD}[]
  homomorphisms(X, Y; kw...) do α
    push!(results, map_components(deepcopy, α)); return false
  end
  results
end
homomorphisms(f, X::AbstractACSet, Y::AbstractACSet;
              monic=false, iso=false, initial=(;)) =
  backtracking_search(f, X, Y, monic=monic, iso=iso, initial=initial)

""" Is the first attributed ``C``-set homomorphic to the second?

A convenience function based on [`homomorphism`](@ref).
"""
is_homomorphic(X::AbstractACSet, Y::AbstractACSet; kw...) =
  !isnothing(homomorphism(X, Y; kw...))

""" Find an isomorphism between two attributed ``C``-sets, if one exists.

See [`homomorphism`](@ref) for more information about the algorithms involved.
"""
isomorphism(X::AbstractACSet, Y::AbstractACSet; initial=(;)) =
  homomorphism(X, Y, iso=true, initial=initial)

""" Find all isomorphisms between two attributed ``C``-sets.

This function is at least as expensive as [`isomorphism`](@ref) and when no
homomorphisms exist, it is exactly as expensive.
"""
isomorphisms(X::AbstractACSet, Y::AbstractACSet; initial=(;)) =
  homomorphisms(X, Y, iso=true, initial=initial)
isomorphisms(f, X::AbstractACSet, Y::AbstractACSet; initial=(;)) =
  homomorphisms(f, X, Y, iso=true, initial=initial)

""" Are the two attributed ``C``-sets isomorphic?

A convenience function based on [`isomorphism`](@ref).
"""
is_isomorphic(X::AbstractACSet, Y::AbstractACSet; kw...) =
  !isnothing(isomorphism(X, Y; kw...))

""" Internal state for backtracking search for ACSet homomorphisms.
"""
struct BacktrackingState{CD <: CatDesc, AD <: AttrDesc{CD},
    Assign <: NamedTuple, PartialAssign <: NamedTuple,
    Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}}
  """ The current assignment, a partially-defined homomorphism of ACSets. """
  assignment::Assign
  """ Depth in search tree at which assignments were made. """
  assignment_depth::Assign
  """ Inverse assignment for monic components or if finding a monomorphism. """
  inv_assignment::PartialAssign
  """ Domain ACSet: the "variables" in the CSP. """
  dom::Dom
  """ Codomain ACSet: the "values" in the CSP. """
  codom::Codom
end

function backtracking_search(f, X::AbstractACSet{CD}, Y::AbstractACSet{CD};
                             monic=false, iso=false, initial=(;)) where {Ob, CD<:CatDesc{Ob}}
  # Fail early if no monic/isos exist on cardinality grounds.
  if iso isa Bool
    iso = iso ? Ob : ()
  end
  for c in iso
    nparts(X,c) == nparts(Y,c) || return false
  end
  if monic isa Bool
    monic = monic ? Ob : ()
  end
  # Injections between finite sets are bijections, so reduce to that case.
  monic = unique([iso..., monic...])
  for c in monic
    nparts(X,c) <= nparts(Y,c) || return false
  end

  # Initialize state variables for search.
  assignment = NamedTuple{Ob}(zeros(Int, nparts(X, c)) for c in Ob)
  assignment_depth = map(copy, assignment)
  inv_assignment = NamedTuple{Ob}(
    (c in monic ? zeros(Int, nparts(Y, c)) : nothing) for c in Ob)
  state = BacktrackingState(assignment, assignment_depth, inv_assignment, X, Y)

  # Make any initial assignments, failing immediately if inconsistent.
  for (c, c_assignments) in pairs(initial)
    for (x, y) in partial_assignments(c_assignments)
      assign_elem!(state, 0, Val{c}, x, y) || return false
    end
  end

  # Start the main recursion for backtracking search.
  backtracking_search(f, state, 1)
end

function backtracking_search(f, state::BacktrackingState, depth::Int)
  # Choose the next unassigned element.
  mrv, mrv_elem = find_mrv_elem(state, depth)
  if isnothing(mrv_elem)
    # No unassigned elements remain, so we have a complete assignment.
    return f(ACSetTransformation(state.assignment, state.dom, state.codom))
  elseif mrv == 0
    # An element has no allowable assignment, so we must backtrack.
    return false
  end
  c, x = mrv_elem

  # Attempt all assignments of the chosen element.
  Y = state.codom
  for y in parts(Y, c)
    assign_elem!(state, depth, Val{c}, x, y) &&
      backtracking_search(f, state, depth + 1) &&
      return true
    unassign_elem!(state, depth, Val{c}, x)
  end
  return false
end

""" Find an unassigned element having the minimum remaining values (MRV).
"""
function find_mrv_elem(state::BacktrackingState{CD}, depth) where CD
  mrv, mrv_elem = Inf, nothing
  Y = state.codom
  for c in ob(CD), (x, y) in enumerate(state.assignment[c])
    y == 0 || continue
    n = count(can_assign_elem(state, depth, Val{c}, x, y) for y in parts(Y, c))
    if n < mrv
      mrv, mrv_elem = n, (c, x)
    end
  end
  (mrv, mrv_elem)
end

""" Check whether element (c,x) can be assigned to (c,y) in current assignment.
"""
function can_assign_elem(state::BacktrackingState, depth,
                         ::Type{Val{c}}, x, y) where c
  # Although this method is nonmutating overall, we must temporarily mutate the
  # backtracking state, for several reasons. First, an assignment can be a
  # consistent at each individual subpart but not consistent for all subparts
  # simultaneously (consider trying to assign a self-loop to an edge with
  # distinct vertices). Moreover, in schemas with non-trivial endomorphisms, we
  # must keep track of which elements we have visited to avoid looping forever.
  ok = assign_elem!(state, depth, Val{c}, x, y)
  unassign_elem!(state, depth, Val{c}, x)
  return ok
end

""" Attempt to assign element (c,x) to (c,y) in the current assignment.

Returns whether the assignment succeeded. Note that the backtracking state can
be mutated even when the assignment fails.
"""
@generated function assign_elem!(state::BacktrackingState{CD,AD}, depth,
                                 ::Type{Val{c}}, x, y) where {CD, AD, c}
  quote
    y′ = state.assignment.$c[x]
    y′ == y && return true  # If x is already assigned to y, return immediately.
    y′ == 0 || return false # Otherwise, x must be unassigned.
    if !isnothing(state.inv_assignment.$c) && state.inv_assignment.$c[y] != 0
      # Also, y must unassigned in the inverse assignment.
      return false
    end

    # Check attributes first to fail as quickly as possible.
    X, Y = state.dom, state.codom
    $(map(out_attr(AD, c)) do f
        :(subpart(X,x,$(quot(f))) == subpart(Y,y,$(quot(f))) || return false)
      end...)

    # Make the assignment and recursively assign subparts.
    state.assignment.$c[x] = y
    state.assignment_depth.$c[x] = depth
    if !isnothing(state.inv_assignment.$c)
      state.inv_assignment.$c[y] = x
    end
    $(map(out_hom(CD, c)) do (f, d)
        :(assign_elem!(state, depth, Val{$(quot(d))}, subpart(X,x,$(quot(f))),
                       subpart(Y,y,$(quot(f)))) || return false)
      end...)
    return true
  end
end

""" Unassign the element (c,x) in the current assignment.
"""
@generated function unassign_elem!(state::BacktrackingState{CD}, depth,
                                   ::Type{Val{c}}, x) where {CD, c}
  quote
    state.assignment.$c[x] == 0 && return
    assign_depth = state.assignment_depth.$c[x]
    @assert assign_depth <= depth
    if assign_depth == depth
      X = state.dom
      if !isnothing(state.inv_assignment.$c)
        y = state.assignment.$c[x]
        state.inv_assignment.$c[y] = 0
      end
      state.assignment.$c[x] = 0
      state.assignment_depth.$c[x] = 0
      $(map(out_hom(CD, c)) do (f, d)
          :(unassign_elem!(state, depth, Val{$(quot(d))},
                           subpart(X,x,$(quot(f)))))
        end...)
    end
  end
end

""" Get assignment pairs from partially specified component of C-set morphism.
"""
partial_assignments(x::AbstractDict) = pairs(x)
partial_assignments(x::AbstractVector) =
  ((i,y) for (i,y) in enumerate(x) if !isnothing(y) && y > 0)

# FIXME: Should these accessors go elsewhere?
in_hom(CD, c) = [dom(CD,f) => f for f in hom(CD) if codom(CD,f) == c]
out_hom(CD, c) = [f => codom(CD,f) for f in hom(CD) if dom(CD,f) == c]
out_attr(AD, c) = [f for f in attr(AD) if dom(AD, f) == c]

# Category of C-sets
####################

@instance Category{ACSet, ACSetTransformation} begin
  dom(α::ACSetTransformation) = α.dom
  codom(α::ACSetTransformation) = α.codom

  id(X::ACSet) = ACSetTransformation(map(id, fin_sets(X)), X, X)

  function compose(α::ACSetTransformation, β::ACSetTransformation)
    # Question: Should we incur cost of checking that codom(β) == dom(α)?
    ACSetTransformation(map(compose, components(α), components(β)),
                        dom(α), codom(β))
  end
end

fin_sets(X::ACSet) = map(table -> FinSet(length(table)), tables(X))

@cartesian_monoidal_instance CSet CSetTransformation
@cocartesian_monoidal_instance ACSet ACSetTransformation

# Limits and colimits
#####################

""" Limit of C-sets that stores the pointwise limits in FinSet.
"""
struct CSetLimit{Ob <: AbstractCSet, Diagram, Cone <: Multispan{Ob},
                 Limits <: NamedTuple} <: AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
  limits::Limits
end

""" Colimit of attributed C-sets that stores the pointwise colimits in FinSet.
"""
struct ACSetColimit{Ob <: AbstractACSet, Diagram, Cocone <: Multicospan{Ob},
                    Colimits <: NamedTuple} <: AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
  colimits::Colimits
end

# Compute limits and colimits of C-sets by reducing to those in FinSet using the
# "pointwise" formula for (co)limits in functor categories.

function limit(diagram::AbstractFreeDiagram{ACS}) where
    {CD <: CatDesc, ACS <: AbstractCSet{CD}}
  limits = map(limit, unpack_diagram(diagram))
  Xs = cone_objects(diagram)
  Y = ACS()
  for (c, lim) in pairs(limits)
    add_parts!(Y, c, length(ob(lim)))
  end
  for (f, c, d) in zip(hom(CD), dom(CD), codom(CD))
    Yfs = map(legs(limits[c]), Xs) do π, X
      compose(π, FinFunction(subpart(X, f), nparts(X, d)))
    end
    Yf = universal(limits[d], Multispan(ob(limits[c]), Yfs))
    set_subpart!(Y, f, collect(Yf))
  end
  πs = pack_components(map(legs, limits), map(X -> Y, Xs), Xs)
  CSetLimit(diagram, Multispan(Y, πs), limits)
end

function universal(lim::CSetLimit, cone::Multispan)
  components = map(universal, lim.limits, unpack_diagram(cone))
  CSetTransformation(components, apex(cone), ob(lim))
end

function colimit(diagram::AbstractFreeDiagram{ACS}) where
    {CD <: CatDesc, AD <: AttrDesc{CD}, Ts, ACS <: AbstractACSet{CD,AD,Ts}}
  # Colimit of C-set without attributes.
  colimits = map(colimit, unpack_diagram(diagram))
  Xs = cocone_objects(diagram)
  Y = ACS()
  for (c, colim) in pairs(colimits)
    add_parts!(Y, c, length(ob(colim)))
  end
  for (f, c, d) in zip(hom(CD), dom(CD), codom(CD))
    Yfs = map(legs(colimits[d]), Xs) do ι, X
      compose(FinFunction(subpart(X, f), nparts(X, d)), ι)
    end
    Yf = universal(colimits[c], Multicospan(ob(colimits[d]), Yfs))
    set_subpart!(Y, f, collect(Yf))
  end
  ιs = pack_components(map(legs, colimits), Xs, map(X -> Y, Xs))

  # Set data attributes by canonical inclusion from attributes in diagram.
  for (attr, c, d) in zip(attr(AD), adom(AD), acodom(AD))
    T = Ts.parameters[d]
    data = Vector{Union{Some{T},Nothing}}(nothing, nparts(Y, c))
    for (ι, X) in zip(ιs, Xs)
      for i in parts(X, c)
        j = ι[c](i)
        if isnothing(data[j])
          data[j] = Some(subpart(X, i, attr))
        else
          val1, val2 = subpart(X, i, attr), something(data[j])
          val1 == val2 || error(
            "ACSet colimit does not exist: $attr attributes $val1 != $val2")
        end
      end
    end
    set_subpart!(Y, attr, map(something, data))
  end

  ACSetColimit(diagram, Multicospan(Y, ιs), colimits)
end

function universal(colim::ACSetColimit, cocone::Multicospan)
  components = map(universal, colim.colimits, unpack_diagram(cocone))
  ACSetTransformation(components, ob(colim), apex(cocone))
end

""" Diagram in C-Set → named tuple of diagrams in FinSet
"""
unpack_diagram(discrete::DiscreteDiagram{<:AbstractACSet}) =
  map(DiscreteDiagram, unpack_sets(ob(discrete)))
unpack_diagram(span::Multispan{<:AbstractACSet}) =
  map(Multispan, fin_sets(apex(span)), unpack_components(legs(span)))
unpack_diagram(cospan::Multicospan{<:AbstractACSet}) =
  map(Multicospan, fin_sets(apex(cospan)), unpack_components(legs(cospan)))
unpack_diagram(para::ParallelMorphisms{<:AbstractACSet}) =
  map(ParallelMorphisms, unpack_components(hom(para)))

function unpack_diagram(d::Union{FreeDiagram{ACS},BipartiteFreeDiagram{ACS}}) where
    {Ob, CD <: CatDesc{Ob}, ACS <: AbstractACSet{CD}}
  NamedTuple{Ob}([ map(d, Ob=X -> FinSet(X, ob), Hom=α -> α[ob]) for ob in Ob ])
end

""" Vector of C-sets → named tuple of vectors of FinSets
"""
unpack_sets(Xs::AbstractVector{<:AbstractACSet{CD}}) where
    {Ob, CD <: CatDesc{Ob}} =
  NamedTuple{Ob}([ map(X -> FinSet(X, ob), Xs) for ob in Ob ])

""" Vector of C-set transformations → named tuple of vectors of FinFunctions
"""
unpack_components(αs::AbstractVector{<:ACSetTransformation{CD}}) where
    {Ob, CD <: CatDesc{Ob}} =
  NamedTuple{Ob}([ map(α -> α[ob], αs) for ob in Ob ])

""" Named tuple of vectors of FinFunctions → vector of C-set transformations
"""
function pack_components(fs::NamedTuple{Ob}, doms, codoms) where Ob
  components = map((x...) -> NamedTuple{Ob}(x), fs...) # XXX: Is there a better way?
  map(ACSetTransformation, components, doms, codoms)
end

""" Objects in diagram that will have explicit legs in limit cone.

Encodes common conventions such as, when taking a pullback of a cospan, not
explicitly including a cone leg for the cospan apex since it can be computed
from the other legs.

FIXME: Should this function be part of the official limits interface?
"""
cone_objects(diagram) = ob(diagram)
cone_objects(diagram::BipartiteFreeDiagram) = ob₁(diagram)
cone_objects(cospan::Multicospan) = feet(cospan)
cone_objects(para::ParallelMorphisms) = SVector(dom(para))

""" Objects in diagram that will have explicit legs in colimit cocone.
"""
cocone_objects(diagram) = ob(diagram)
cocone_objects(diagram::BipartiteFreeDiagram) = ob₂(diagram)
cocone_objects(span::Multispan) = feet(span)
cocone_objects(para::ParallelMorphisms) = SVector(codom(para))

# Sub-C-sets
############

const SubCSet{CD} = Subobject{<:AbstractCSet{CD}}
const SubACSet{CD,AD} = Subobject{<:AbstractACSet{CD,AD}}

components(A::SubACSet) = map(Subobject, components(hom(A)))
force(A::SubACSet) = Subobject(force(hom(A)))

""" Construct subobject of C-set from components of inclusion or predicate.

This function constructs a subobject from the components of the inclusion
morphism, given as a named tuple or as keyword arguments. The components can
also be specified as predicates (boolean vectors).
"""
function Subobject(X::T, components) where T <: AbstractACSet
  U = T()
  components = map(coerce_subob_component, components)
  copy_parts!(U, X, components)
  Subobject(ACSetTransformation(components, U, X))
end
Subobject(X::AbstractACSet; components...) = Subobject(X, (; components...))

coerce_subob_component(f::FinFunction{Int}) = collect(f)
coerce_subob_component(A::SubFinSet{Int}) = collect(hom(A))
coerce_subob_component(f::AbstractVector{Int}) = f
coerce_subob_component(pred::Union{AbstractVector{Bool},BitVector}) =
  findall(pred)

@instance SubobjectBiHeytingAlgebra{AbstractACSet,SubACSet} begin
  @import ob
  meet(A::SubACSet, B::SubACSet) = meet(A, B, SubOpBoolean())
  join(A::SubACSet, B::SubACSet) = join(A, B, SubOpBoolean())
  top(X::AbstractACSet) = top(X, SubOpWithLimits())
  bottom(X::AbstractACSet) = bottom(X, SubOpWithLimits())

  implies(A::SubACSet, B::SubACSet) = implies(A, B, SubOpBoolean())
  subtract(A::SubACSet, B::SubACSet) = subtract(A, B, SubOpBoolean())
  negate(A::SubACSet) = implies(A, bottom(ob(A)), SubOpBoolean())
  non(A::SubACSet) = subtract(top(ob(A)), A, SubOpBoolean())
end

function meet(A::SubACSet, B::SubACSet, ::SubOpBoolean)
  Subobject(common_ob(A, B), map(components(A), components(B)) do A₀, B₀
    meet(A₀, B₀, SubOpBoolean())
  end)
end
function join(A::SubACSet, B::SubACSet, ::SubOpBoolean)
  Subobject(common_ob(A, B), map(components(A), components(B)) do A₀, B₀
    join(A₀, B₀, SubOpBoolean())
  end)
end
top(X::AbstractACSet, ::SubOpBoolean) =
  Subobject(X, map(X₀ -> top(X₀, SubOpBoolean()), fin_sets(X)))
bottom(X::AbstractACSet, ::SubOpBoolean) =
  Subobject(X, map(X₀ -> bottom(X₀, SubOpBoolean()), fin_sets(X)))

""" Implication of sub-C-sets.

By (Reyes et al 2004, Proposition 9.1.5), the implication ``A ⟹ B`` for two
sub-``C``-sets ``A,B ↪ X`` is given by

``x ∈ (A ⟹ B)(c) iff ∀f: c → c′, x⋅f ∈ A(c′) ⟹ x⋅f ∈ B(c′)``

for all ``c ∈ C`` and ``x ∈ X(c)``. By the definition of implication and De
Morgan's law in classical logic, this is equivalent to

``x ∉ (A ⟹ B)(c) iff ∃f: c → c′, x⋅f ∈ A(c′) ∧ x⋅f ∉ B(c′)``.

In this form, we can clearly see the duality to formula and algorithm for
subtraction of sub-C-sets ([`subtract`](@ref)).
"""
function implies(A::SubACSet{CD}, B::SubACSet{CD}, ::SubOpBoolean) where CD
  X = common_ob(A, B)
  A, B = map(as_predicate, components(A)), map(as_predicate, components(B))
  D = map(X₀ -> trues(length(X₀)), fin_sets(X))

  function unset!(c, x)
    D[c][x] = false
    for (c′,x′) in all_incident(X, Val{c}, x)
      if D[c′][x′]; unset!(c′,x′) end
    end
  end

  for c in ob(CD), x in parts(X,c)
    if D[c][x] && A[c][x] && !B[c][x]; unset!(c,x) end
  end
  Subobject(X, D)
end

""" Subtraction of sub-C-sets.

By (Reyes et al 2004, Sec 9.1, pp. 144), the subtraction ``A ∖ B`` for two
sub-``C``-sets ``A,B ↪ X`` is given by

``x ∈ (A ∖ B)(c) iff ∃f: c′ → c, ∃x′ ∈ f⁻¹⋅x, x′ ∈ A(c′) ∧ x′ ∉ B(c′)``

for all ``c ∈ C`` and ``x ∈ X(c)``. Compare with [`implies`](@ref).
"""
function subtract(A::SubACSet{CD}, B::SubACSet{CD}, ::SubOpBoolean) where CD
  X = common_ob(A, B)
  A, B = map(as_predicate, components(A)), map(as_predicate, components(B))
  D = map(X₀ -> falses(length(X₀)), fin_sets(X))

  function set!(c, x)
    D[c][x] = true
    for (c′,x′) in all_subparts(X, Val{c}, x)
      if !D[c′][x′]; set!(c′,x′) end
    end
  end

  for c in ob(CD), x in parts(X,c)
    if !D[c][x] && A[c][x] && !B[c][x]; set!(c,x) end
  end
  Subobject(X, D)
end

function common_ob(A::Subobject, B::Subobject)
  (X = ob(A)) == ob(B) ||
    error("Subobjects have different base objects: $(ob(A)) != $(ob(B))")
  return X
end

# FIXME: Should these two accessors go elsewhere?

@generated function all_subparts(X::AbstractACSet{CD},
                                 ::Type{Val{c}}, x::Int) where {CD,c}
  Expr(:tuple, map(out_hom(CD, c)) do (f,c′)
    :($(quot(c′)), subpart(X,x,$(quot(f))))
  end...)
end

@generated function all_incident(X::AbstractACSet{CD},
                                 ::Type{Val{c}}, x::Int) where {CD,c}
  Expr(:call, GlobalRef(Iterators, :flatten),
    Expr(:tuple, map(in_hom(CD, c)) do (c′,f)
      :(Tuple{Symbol,Int}[ ($(quot(c′)),x′) for x′ in incident(X,x,$(quot(f))) ])
    end...))
end

# Serialization
###############

""" Serialize an ACSet object to a JSON string
"""
function generate_json_acset(x::T) where T <: AbstractACSet
  JSON.json(x.tables)
end

""" Deserialize a dictionary from a parsed JSON string to an object of the given ACSet type
"""
function parse_json_acset(type::Type{T}, input::Dict) where T <: AbstractACSet
  out = type()
  for (k,v) ∈ input
    add_parts!(out, Symbol(k), length(v))
  end
  for l ∈ values(input)
    for (i, j) ∈ enumerate(l)
      for (k,v) ∈ j
        vtype = eltype(out[Symbol(k)])
        if !(typeof(v) <: vtype)
          v = vtype(v)
        end
        set_subpart!(out, i, Symbol(k), v)
      end
    end
  end
  out
end

""" Deserialize a JSON string to an object of the given ACSet type
"""
function parse_json_acset(type::Type{T}, input::String) where T <: AbstractACSet
  parse_json_acset(type, JSON.parse(input))
end

""" Read a JSON file to an object of the given ACSet type
"""
function read_json_acset(type::Type{T}, input::String) where T <: AbstractACSet
  parse_json_acset(type, open(f->read(f, String), input))
end

""" Serialize an ACSet object to a JSON file
"""
function write_json_acset(x::T, fname::AbstractString) where T <: AbstractACSet
  open(string(fname), "w") do f
    write(f, generate_json_acset(x))
  end
end

end


text:
function ChainRulesCore.rrule(
    config::RuleConfig{>:HasReverseMode},
    ::typeof(__foldl__),
    rf,
    init,
    xs::AbstractArray,
)
    rf_fake = ReduceSomePullbackXF(config)'(rf)
    pbtype = Core.Compiler.return_type(
        unreduced_or_throw ∘ foldl_nocomplete,
        Tuple{Typeof(rf_fake),Typeof(init),Typeof(xs)},
    )

    # TODO: Use StructArrays to store `next` pullbacks more efficiently?
    pullbacks = similar(xs, pbtype)
    rf_rr = ArrayPullbackXF(config)'(rf)
    acc, complete_pullback = __foldl__(rf_rr, init, zip(xs, referenceable(pullbacks)))

    function foldl_pullback(ȳ)
        _, r̄f, ācc = complete_pullback(ȳ)
        # TODO: Lazily allocate `x̄s` for
        x̄s = similar(xs, float(eltype(xs)))
        @inline function rf_pb((r̄f, ācc, allnotan), i)
            @inbounds pb = pullbacks[i]
            _, r̄f′, ācc′, x̄ = pb(ācc)
            @inbounds x̄s[i] = x̄
            allnotan &= x̄ isa NoTangent
            return (add!!(r̄f, r̄f′), ācc′, allnotan)
        end
        r̄f, ācc, allnotan = __foldl__(
            Completing(rf_pb),
            (r̄f, ācc, true),
            reverse(eachindex(x̄s, pullbacks)),
        )
        # TODO: Define `reverse(zip(arrays...))` or something equivalent
        # ācc = __foldl__(..., reverse(zip(referenceable(x̄s), pullbacks)))
        return NoTangent(), r̄f, ācc, allnotan ? NoTangent() : x̄s
    end
    return acc, foldl_pullback
end

"""
    ArrayPullbackXF(config::RuleConfig{>:HasReverseMode})

Pullback of reducing function (aka loop body) as a transducer. The "Array"
prefix indicates that the sotrage of the pullbacks are allocated externally as
an array.
"""
struct ArrayPullbackXF{Config} <: Transducer
    config::Config
end

@inline function Transducers.next(rf::R_{ArrayPullbackXF}, acc, (x, r))
    acc′, pb = @inlinecall rrule_via_ad(xform(rf).config, next, inner(rf), acc, x)
    r[] = pb
    return acc′
end

@inline Transducers.complete(rf::R_{ArrayPullbackXF}, acc) =
    rrule_via_ad(xform(rf).config, complete, inner(rf), acc)

#=
"""
    GenericPullbackXF

Generic fallback implementation (but how to construct the output of
`foldl_pullback`?).
"""
struct GenericPullbackXF{Config} <: Transducer
    config::Config
end

@inline function Transducers.next(rf::R_{GenericPullbackXF}, (pullbacks, acc), x)
    acc′, pb = @inlinecall rrule_via_ad(xform(rf).config, next, inner(rf), acc, x)
    push!(pullbacks, pb)
    return pullbacks, acc′
end

@inline Transducers.complete(rf::R_{GenericPullbackXF}, acc) =
    rrule_via_ad(xform(rf).config, complete, inner(rf), acc)
=#

"""
    ReduceSomePullbackXF(config::RuleConfig{>:HasReverseMode})

Transducer that reduces to a pullback of `next` at unspecified point. Used only
for inference.
"""
struct ReduceSomePullbackXF{Config} <: Transducer
    config::Config
end

const _SOME_BOOL = Ref(false)

Transducers.next(rf::R_{ReduceSomePullbackXF}, acc, x) =
    if _SOME_BOOL[]
        _, pb = rrule_via_ad(xform(rf).config, next, inner(rf), acc, x)
        reduced(pb)
    else
        next(inner(rf), acc, x)
    end

Transducers.complete(rf::R_{ReduceSomePullbackXF}, _) = nothing

unreduced_or_throw(x::Reduced) = unreduced(x)
unreduced_or_throw(::Any) = error("unreachable")


text:
abstract type AbstractAcquisition end

setparams!(a, model) = nothing
function acquisitionfunction(a, model)
    x -> begin
        μ, σ² = mean_var(model, x)
        a.(μ, σ²)
    end
end

"""
The probability of improvement measures the probability that a point `x` leads
to an improvement upon an incumbent target `τ`. For Gaussian distributions it is
given by

    Φ[(μ(x) - τ)/σ(x)]

where `Φ` is the standard normal cumulative distribution function and `μ(x)`, `σ(x)`
are mean and standard deviation of the distribution at point `x`.
"""
mutable struct ProbabilityOfImprovement <: AbstractAcquisition
    τ::Float64
end
@inline function (a::ProbabilityOfImprovement)(μ, σ²)
    σ² == 0 && return float(μ > a.τ)
    normal_cdf(μ - a.τ, σ²)
end
ProbabilityOfImprovement(; τ = -Inf) = ProbabilityOfImprovement(τ)

"""
The expected improvement measures the expected improvement `x - τ` of a point `x`
upon an incumbent target `τ`. For Gaussian distributions it is given by

    (μ(x) - τ) * ϕ[(μ(x) - τ)/σ(x)] + σ(x) * Φ[(μ(x) - τ)/σ(x)]

where `ϕ` is the standard normal distribution function and `Φ` is the standard
normal cumulative function, and `μ(x)`, `σ(x)` are mean and standard deviation
of the distribution at point `x`.
"""
mutable struct ExpectedImprovement <: AbstractAcquisition
    τ::Float64
end
ExpectedImprovement(; τ = -Inf) = ExpectedImprovement(τ)
function setparams!(a::Union{ExpectedImprovement,ProbabilityOfImprovement}, model)
    a.τ = max(maxy(model), a.τ)
end
@inline function (a::ExpectedImprovement)(μ, σ²)
    σ² == 0 && return μ > a.τ ? μ - a.τ : 0.
    abs(μ - a.τ) * normal_cdf(μ - a.τ, σ²) + √σ² * normal_pdf(μ - a.τ, σ²)
end

abstract type BetaScaling end
"""
Scales `βt` of `UpperConfidenceBound` as

    βt = √(2 * log(t^(D/2 + 2) * π^2/(3δ)))

where `t` is the number of observations, `D` is the dimensionality of the input
data points and δ is a small constant (default δ = 0.1).

See Brochu E., Cora V. M., de Freitas N. (2010), "A Tutorial on Bayesian
Optimization of Expensive Cost Functions, with Application to Active User
Modeling and Hierarchical Reinforcement Learning", https://arxiv.org/abs/1012.2599v1
page 16.
"""
struct BrochuBetaScaling <: BetaScaling
    δ::Float64
end
"""
Applies no scaling to `βt` of `UpperConfidenceBound`.
"""
struct NoBetaScaling <: BetaScaling end
"""
For Gaussian distributions the upper confidence bound at `x` is given by
    μ(x) + βt * σ(x)

where `βt` is a fixed parameter in the case of `NoBetaScaling` or an observation
size dependent parameter in the case of e.g. `BrochuBetaScaling`.
"""
mutable struct UpperConfidenceBound{S} <: AbstractAcquisition
    scaling::S
    βt::Float64
end
"""
    UpperConfidenceBound(; scaling = BrochuBetaScaling(.1), βt = 1)
"""
UpperConfidenceBound(; scaling = BrochuBetaScaling(.1), βt = 1.) = UpperConfidenceBound(scaling, βt)
function setparams!(a::UpperConfidenceBound{BrochuBetaScaling}, model)
    D, nobs = dims(model)
    nobs == 0 && (nobs = 1)
    a.βt = sqrt(2*log(nobs^(D/2 + 2)*π^2/(3*a.scaling.δ)))
end
(a::UpperConfidenceBound)(μ, σ²) = μ + a.βt * √σ²

"""
The acquisition function associated with `ThompsonSamplingSimple` draws
independent samples for each input `x` a function value from the model. Together
with a gradient-free optimization method this leads to proposal points that
might be similarly distributed as the maxima of true Thompson samples from GPs.
True Thompson samples from a GP are simply functions from a GP. Maximizing these
samples can be tricky, see e.g. http://hildobijl.com/Downloads/GPRT.pdf
chapter 6.
"""
struct ThompsonSamplingSimple <: AbstractAcquisition end
acquisitionfunction(a::ThompsonSamplingSimple, model) = x -> myrand(model, x)

struct MaxMean <: AbstractAcquisition end
acquisitionfunction(a::MaxMean, model) = x -> mean_var(model, x)[1]

"""
The mutual information measures the amount of information gained by querying at
x. The parameter γ̂ gives a lower bound for the information on f from the queries
{x}. For a Gaussian this is
    γ̂ = ∑σ²(x)
and the mutual information at x is
    μ(x) + √(α)*(√(σ²(x)+γ̂) - √(γ̂))

where `μ(x)`, `σ(x)` are mean and standard deviation
of the distribution at point `x`.

See Contal E., Perchet V., Vayatis N. (2014), "Gaussian Process Optimization
with Mutual Information" http://proceedings.mlr.press/v32/contal14.pdf
"""
mutable struct MutualInformation <: AbstractAcquisition
    sqrtα::Float64
    γ̂::Float64
end
MutualInformation(; α = 1.0, γ̂ = 0.0) = MutualInformation(sqrt(α), γ̂)
function setparams!(a::MutualInformation, model)
    D, nobs = dims(model)
    if iszero(nobs)
        a.γ̂ = 0.0
    else
        last_x = @view model.x[:, end]
        μ, σ2 = mean_var(model, last_x)
        a.γ̂ += σ2
    end
end
(a::MutualInformation)(μ, σ²) = μ + a.sqrtα * (sqrt(σ² + a.γ̂) - sqrt(a.γ̂))

# TODO see
# https://github.com/HildoBijl/GPRT/blob/7166548b8587201fabc671a0647aac2ff96f3555/Chapter6/Chapter6.m#L723
# and corresponding thesis page 173
# mutable struct ThompsonSampling{K} <: AbstractAcquisition
# end

# TODO
# mutable struct EntropySearch <: AbstractAcquisition
# end

# TODO
# mutable struct PredictiveEntropySearch <: AbstractAcquisition
# end


text:
@inline function generic_qkv_attention(mixingf, scoref, q, k, v, args...)
    return unwrap_collapse(mixing(mixingf, v, scoref, q, k, args...))
end

@inline function generic_multihead_qkv_attention(mixingf, scoref, head, q, k, v, args...)
    hq = CollapsedDimArray((move_head_dim_out ∘ split_head)(head, q), static(2), static(ndims(q)))
    hk = CollapsedDimArray((move_head_dim_out ∘ split_head)(head, k), static(2), static(ndims(k)))
    hv = CollapsedDimArray((move_head_dim_out ∘ split_head)(head, v), static(2), static(ndims(v)))
    a = generic_qkv_attention(mixingf, scoref, hq, hk, hv, args...)
    return (merge_head ∘ move_head_dim_in)(a)
end

@inline naive_qkv_attention(q, k, v, mask=nothing) = unwrap_collapse(weighted_sum_mixing(
    normalized_score(NNlib.softmax, masked_score, GenericAttenMaskOp(), mask, scaled_dot_product_score, q, k),
    v))

@inline function multihead_qkv_attention(head, q, k, v, mask=nothing)
    hq = CollapsedDimArray((move_head_dim_out ∘ split_head)(head, q), static(2), static(ndims(q)))
    hk = CollapsedDimArray((move_head_dim_out ∘ split_head)(head, k), static(2), static(ndims(k)))
    hv = CollapsedDimArray((move_head_dim_out ∘ split_head)(head, v), static(2), static(ndims(v)))
    a = naive_qkv_attention(hq, hk, hv, mask)
    return (merge_head ∘ move_head_dim_in)(a)
end

generic_attention(mixingf, scoref, q::AbstractArray, k::AbstractArray, v::AbstractArray, args...) = generic_qkv_attention(mixingf, scoref, q, k, v)
generic_attention(mixingf, scoref, head::Integer, q::AbstractArray, k::AbstractArray, v::AbstractArray, args...) = generic_multihead_qkv_attention(mixingf, scoref, head, q, k, v)
attention(head::Integer, q::AbstractArray, k::AbstractArray, v::AbstractArray, mask=nothing) = multihead_qkv_attention(head, q, k, v, mask)
attention(q::AbstractArray, k::AbstractArray, v::AbstractArray, mask=nothing) = naive_qkv_attention(q, k, v, mask)


text:
using Spark

function basic_disc_fault(master::Spark.Master)
    rdd = Spark.input(master, "RDDA.txt", "int_reader")
    master.workers[1].active = false
    mapped_rdd = Spark.map(master, rdd, "double_map")
    @assert mapped_rdd == false
    # recovery
    new_rdd = Spark.recover(master, rdd)
    collection = Spark.collect(master, new_rdd)
    new_mapped_rdd = Spark.map(master, new_rdd, "double_map")
    @assert Spark.count(master, new_rdd) == Spark.count(master, new_mapped_rdd)
    for kv in collection
        values = Spark.lookup(master, new_mapped_rdd, kv[1])
        @assert 2*(kv[2][1]) == values[1]
    end
    master.workers[1].active = true
end

function join_disc_fault(master::Spark.Master)
    rdd_a = Spark.input(master, "RDDA.txt", "int_reader")
    rdd_b = Spark.input(master, "RDDA.txt", "int_reader")
    # first join
    joined_rdd = Spark.join(master, rdd_a, rdd_b)
    # first collection
    joined_collection = Spark.collect(master, joined_rdd)
    # failure
    master.workers[1].active = false
    # first map
    mapped_rdd = Spark.map(master, joined_rdd, "double_map")
    @assert mapped_rdd == false
    # recover
    new_joined_rdd = Spark.recover(master, joined_rdd)
    # second collect
    joined_collection = Spark.collect(master, new_joined_rdd)
    @assert new_joined_rdd != false
    new_mapped_rdd = Spark.map(master, new_joined_rdd, "double_map")
    @assert new_mapped_rdd != false
    @assert Spark.count(master, new_joined_rdd) == 10
    for kv in joined_collection
        @assert length(kv[2]) == 2
    end
    master.workers[1].active = true
end


text:
mutable struct AllSubsetRegressionResult <: ModelSelection.ModelSelectionResult

    datanames
    modelavg_datanames

    data
    bestresult_data
    modelavg_data

    fixedvariables
    outsample
    criteria
    modelavg
    ttest
    residualtest
    orderresults
    nobs

    function AllSubsetRegressionResult(
            datanames,
            modelavg_datanames,
            fixedvariables,
            outsample,
            criteria,
            modelavg,
            ttest,
            residualtest,
            orderresults
        )

        new(datanames, modelavg_datanames, nothing, nothing, nothing, fixedvariables, outsample, criteria, modelavg, ttest, residualtest, orderresults, 0)
    end
end


text:
#= Most of the code in this file is based on code by Professor Jean Michel (reused here with his
permission). See the following paper for his context and for his code for the C version of GAP3:

Jean Michel. The development version of the CHEVIE package of GAP3. J. Algebra, 435:308–336, 2015.
=#

#=
The aim of the code below is to test the main conjecture of the author's paper on minimally
dominant elements (see arXiv:2110.09266).

The input is as follows:
. A Weyl group W
. A list L of pairs of words (in GAP notation) representing two elements in a conjugacy class of W,
    one is of minimal length and the other one is minimally dominant. The pairs in this list are
    separated by the string "end cc".
. A prime number p

Such a list L can be computed with the other file in this repository.
Running test_main_conjecture(W,L,p) then verifies whether the corresponding Bruhat cells intersect
the same conjugacy classes in the associated (simply connected) reductive group, over an
algebraically closed field of characteristic p.

Input example:
> W = coxgroup(:B,3)
> L = [[2, 1, 3, 2, 1, 3, 2, 3] ,
       [1, 2, 1, 2] ,
       "end cc"  ,
       [1, 2, 1, 3, 2, 1, 3, 2] ,
       [1, 3] ,
       "end cc"  ,
       [2, 1, 2, 3] ,
       [2, 1] ,
       "end cc"  ,
       [1, 3, 2, 1, 3, 2] ,
       [3, 2] ,
       "end cc"  ,
       [2, 1, 3, 2, 1, 3, 2] ,
       [3] ,
       "end cc"  ,
       [3, 2, 1, 2, 3] ,
       [1] ,
       "end cc"]
> p = 7
> test_main_conjecture(W, L, p)
=#

using Gapjm     # https://github.com/jmichel7/Gapjm.jl
using LinearAlgebra

function f(W, wrd, q)
    permutedims(
        toM(
            DLLefschetz(Tbasis(hecke(W, q))(W(wrd...))).v[
                UnipotentCharacters(W).harishChandra[1][:charNumbers]]))
end

function g(W, q, p)
    uc = UnipotentClasses(W, p)
    t = ICCTable(uc, 1; q)
    triv = map(1:length(uc)) do i
        findfirst(==([i, charinfo(uc.classes[i].Au)[:positionId]]),
            uc.springerseries[1][:locsys])
    end
    permutedims(t.L[triv, triv] * t.scalar[triv, :] * diagm(q .^ t.dimBu))
end

function vdash(W, wrd, p)
    q = Pol()^2
    m = f(W, wrd, q) * g(W, q, p)
    map(j -> axes(m, 2)[(!iszero).(m[j, :])], axes(m, 1))
end

function test_main_conjecture(W, L, p)
    indices_of_conjugacy_classes = []
    count_cc = 0
    for wrd in L
        if wrd == "end cc"
            count_cc += 1
            if length(Set(indices_of_conjugacy_classes)) != 1
                print("Counterexample at conjugacy class ", count_cc, " in L\n")
            end
            indices_of_conjugacy_classes = []
        else
            append!(indices_of_conjugacy_classes, vdash(W, wrd, p))
        end
    end
end

max_stars_repo_name: waamm/minimally-dominant

max_issues_repo_name: waamm/minimally-dominant

max_forks_repo_name: waamm/minimally-dominant


text:
@testset "ModelParameters" begin
    @test mP_1.U ≈ 1.1
    @test mP_1.μ ≈ 1.2
    @test mP_1.β ≈ 1.3
    @test mP_1.n ≈ 1.4
end

@testset "SimulationParameters" begin
    @test sP_1.n_iω == 1
    @test sP_1.n_iν == 2
    @test sP_1.shift == false
    @test sP_1.tc_type_f == :nothing
    @test sP_1.tc_type_b == :nothing
    @test sP_1.λc_type == :nothing
    @test sP_1.ωsum_type == :common
    @test sP_1.λ_rhs == :native
    @test sP_1.fullChi == false
    @test sP_1.χFillType == LadderDGA.zero_χ_fill
    @test sP_1.bosonic_tail_coeffs == [0,1,2,3]
    @test sP_1.fermionic_tail_coeffs == [0,1,2,3]
    @test sP_1.usable_prct_reduction == 0.1
end


text:
__precompile__()

module Xpress

    fn = joinpath(dirname(@__FILE__),"..","deps","deps.jl")
    if isfile(fn)
        include(fn)
    else
        error("Xpress not properly installed. Please run Pkg.build(\"Xpress\")")
    end

    ### imports

    using Compat
    using Compat.SparseArrays
    using Compat.LinearAlgebra

    import Base.show, Base.copy

    ### exports
    export

    ## xprs_env
    free_env,

    ## xprs_params
    getparam, setparam!, setparams!,

    ## xprs_model
    set_sense!, update_model!, reset_model!,
    read_model, write_model,
    copy,

    ## xprs_attrs
    model_name, model_sense, model_type,
    num_vars, num_constrs, num_sos, num_qconstrs,
    num_cnzs, num_qnzs, num_qcnzs,
    is_qp, is_qcp, is_mip,

    lowerbounds, upperbounds, objcoeffs, set_objcoeffs!,

    ## xprs_vars
    XPRS_CONTINUOUS, XPRS_BINARY, XPRS_INTEGER,
    add_var!, add_vars!, add_cvar!, add_cvars!,
    add_bvar!, add_bvars!, add_ivar!, add_ivars!,
    del_vars!,

    ## xprs_constrs
    add_constr!, add_constrs!, add_constrs_t!,
    add_rangeconstr!, add_rangeconstrs!, add_rangeconstrs_t!,
    get_constrmatrix, add_sos!, del_constrs!, chg_coeffs!,

    ## xprs_quad
    add_qpterms!, add_qconstr!,

    ## higher level
    xpress_model,

    ## xprs_solve
    optimize, computeIIS, get_solution,
    get_status, OptimInfo, get_optiminfo, get_objval

    ### include source files

    include("xprs_full_defines.jl")

    include("xprs_common.jl")
    include("xprs_env.jl")

    include("xprs_model.jl")
    include("xprs_params.jl")
    include("xprs_attrs.jl")

    include("xprs_vars.jl")

    include("xprs_constrs.jl")
    include("xprs_quad.jl")
    include("xprs_highlevel.jl")

    include("xprs_solve.jl")
    include("xprs_callbacks.jl")
    include("xprs_iis.jl")

    include("XpressSolverInterface.jl")
    include("MOIWrapper.jl")

    # license checker
    include("xprs_userlic.jl")

    function __init__()

        # some lics require special check
        userlic()

        # start Xpress with XPRSinit
        Env()
    end
end


text:
import LinearAlgebra: I
import RungeKutta: @big, _big, istriustrict, istrilstrict

@testset "$(rpad("Utility Functions",80))" begin

    @test _big(1)    == BigFloat(1)
    @test _big(1//1) == BigFloat(1)
    @test _big(1.0)  == BigFloat(1)
    @test _big("1")  == BigFloat(1)
    @test _big(:(1)) == BigFloat(1)

    x = @big [1  2.0  1//3]
    y = big.([1  2.0  1//3])
    z = [big(1)  big(2.0)  big(1//3)]

    @test x != y
    @test x == z

    for n in 1:5
        A = zeros(n,n)
        B = zeros(n,n)

        for i in 1:n
            for j in i+1:n
                A[i,j] = rand()
                B[j,i] = rand()
            end
        end

        @test istriustrict(A)
        @test istrilstrict(B)

        @test !istriustrict(A .+ Matrix(I, n, n))
        @test !istrilstrict(B .+ Matrix(I, n, n))

        @test !istriustrict(rand(n,n))
        @test !istrilstrict(rand(n,n))
    end

end


text:
# Declaring variables and using data types in Julia

# 5 basic types of data:
# Integers, Floating point numbers, booleans, strings, and characters

# Variable names typically are lowercase. Use underscores if the name
# is too long to be easily read.
first = 1
second = 2.0
stringVar = "This is a string"
initial = 'J'
boolVal = true

println(first)
println(second)
println(stringVar)
println(initial)
println(boolVal)

# Variables can be redefined and change type
first = 1.0
println(first)

# Strings have to be in double quotes. This will error: anotherString = 'hello world'

# Constant values are declared with const and are all uppercase
# Note: this is only useful in the global scope
const CONSTANT = 42

# Constants of the same type *can* be reassigned, but with a warning
CONSTANT = 43

# Constants of different types cannot be reassigned, this is an error
# Cannot do: CONSTANT = 43.0

# Type Annotation identifies a variable as a particular type
function testFunc()
    x::String = "some string"
end

testFunc()

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
export abm_data_exploration

"""
    abm_data_exploration(model::ABM, agent_step!, model_step!, params=Dict(); kwargs...)
Open an interactive application for exploring an agent based model and
the impact of changing parameters on the time evolution. Requires `Agents`.

The application evolves an ABM interactively and plots its evolution, while allowing
changing any of the model parameters interactively and also
showing the evolution of collected data over time (if any are asked for, see below).
The agent based model is plotted and animated exactly as in [`abm_play`](@ref),
and the arguments `model, agent_step!, model_step!` are propagated there as-is.

Calling `abm_data_exploration` returns: `figure, agent_df, model_df`. So you can save the
figure, but you can also access the collected data (if any).

## Interaction
Besides the basic time evolution interaction of [`abm_play`](@ref), additional
functionality here allows changing model parameters in real time, based on the provided
fourth argument `params`. This is a dictionary which decides which
parameters of the model will be configurable from the interactive application.
Each entry of `params` is a pair of `Symbol` to an `AbstractVector`, and provides a range
of possible values for the parameter named after the given symbol (see example online).
Changing a value in the parameter slides is only updated into the actual model when
pressing the "update" button.

The "reset" button resets the model to its original agent and space state but it updates
it to the currently selected parameter values. A red vertical line is displayed in the
data plots when resetting, for visual guidance.

## Keywords

* `ac, am, as, scheduler, offset, aspect, scatterkwargs`: propagated to [`abm_plot`](@ref).
* `adata, mdata`: Same as the keyword arguments of `Agents.run!`, and decide which data of the
  model/agents will be collected and plotted below the interactive plot. Notice that
  data collection can only occur on plotted steps (and thus steps not plotted due to
  "spu" are also not data-collected).
* `alabels, mlabels`: If data are collected from agents or the model with `adata, mdata`,
  the corresponding plots have a y-label named after the collected data. Instead, you can
  give `alabels, mlabels` (vectors of strings with exactly same length as `adata, mdata`),
  and these labels will be used instead.
* `when = true`: When to perform data collection, as in `Agents.run!`.
* `spu = 1:100`: Values that the "spu" slider will obtain.
"""
function abm_data_exploration(
        model, agent_step!, model_step!, params = Dict();
        mdata = nothing,
        adata = nothing,
        alabels = nothing,
        mlabels = nothing,
        when = true,
        spu = 1:100,
        kwargs...
    )

    # preinitialize a bunch of stuff
    model0 = deepcopy(model)
    modelobs = Observable(model) # only useful for resetting
    ac = get(kwargs, :ac, JULIADYNAMICS_COLORS[1])
    as = get(kwargs, :as, 10)
    am = get(kwargs, :am, :circle)
    scheduler = get(kwargs, :scheduler, model.scheduler)
    offset = get(kwargs, :offset, nothing)

    !isnothing(adata) && @assert adata[1] isa Tuple "Only aggregated agent data are allowed."
    !isnothing(alabels) && @assert length(alabels) == length(adata)
    !isnothing(mlabels) && @assert length(mlabels) == length(mdata)
    df_agent = Agents.init_agent_dataframe(model, adata)
    df_model = Agents.init_model_dataframe(model, mdata)
    L = (isnothing(adata) ? 0 : size(df_agent)[2]-1) + (isnothing(mdata) ? 0 : size(df_model)[2]-1)
    s = 0 # current step
    P = length(params)

    # Initialize main layout and abm axis
    figure = Figure(; resolution = (1600, 800), backgroundcolor = DEFAULT_BG)
    abmax = figure[1,1][1,1] = if dimensionality(model) == 3
        Axis3(figure)
    else
        Axis(figure)
    end

    # Initialize the ABM plot stuff
    abmstepper = abm_init_stepper_and_plot!(abmax, figure, model; kwargs...)
    speed, slep, step, run, reset, update = abm_controls_play!(figure, model, spu, true)

    # Initialize parameter controls & data plots
    datalayout = figure[:, 2] = GridLayout(tellheight = false)
    slidervals = abm_param_controls!(figure, datalayout, model, params, L)
    if L > 0
        N = Observable([0]) # steps that data are recorded at.
        data, axs = init_abm_data_plots!(
            figure, datalayout, model, df_agent, df_model, adata, mdata, N, alabels, mlabels
        )
    end

    # Clicking the step button
    on(step) do clicks
        n = speed[]
        Agents.step!(abmstepper, model, agent_step!, model_step!, n)
        if L > 0
            s += n
            if L > 0 && Agents.should_we_collect(s, model, when) # update collected data
                push!(N.val, s)
                update_abm_data_plots!(data, axs, model, df_agent, df_model, adata, mdata, N)
            end
        end
    end

    # Clicking the run button
    isrunning = Observable(false)
    on(run) do clicks; isrunning[] = !isrunning[]; end
    on(run) do clicks
        @async while isrunning[]
            n = speed[]
            model = modelobs[] # this is useful only for the reset button
            Agents.step!(abmstepper, model, agent_step!, model_step!, n)
            if L > 0
                s += n
                if L > 0 && Agents.should_we_collect(s, model, when) # update collected data
                    push!(N.val, s)
                    update_abm_data_plots!(data, axs, model, df_agent, df_model, adata, mdata, N)
                end
            end
            slep[] == 0 ? yield() : sleep(slep[])
            isopen(figure.scene) || break # crucial, ensures computations stop if closed window.
        end
    end

    # Clicking the reset button
    on(reset) do clicks
        modelobs[] = deepcopy(model0)
        Agents.step!(abmstepper, model, agent_step!, model_step!, 0)
        L > 0 && add_reset_line!(axs, s)
        update[] = update[] + 1 # also trigger parameter updates
    end

    # Clicking the update button
    on(update) do clicks
        model = modelobs[]
        update_abm_properties!(model, slidervals)
    end

    display(figure)
    return figure, df_agent, df_model
end

function abm_param_controls!(figure, datalayout, model, params, L)
    slidervals = Dict{Symbol, Observable}()
    for (i, (l, vals)) in enumerate(params)
        startvalue = has_key(model.properties, l) ?
            get_value(model.properties, l) : vals[1]
        sll = labelslider!(figure, string(l), vals; sliderkw = Dict(:startvalue => startvalue))
        slidervals[l] = sll.slider.value # directly add the observable
        datalayout[i+L, :] = sll.layout
    end
    return slidervals
end

function init_abm_data_plots!(
        figure, datalayout, model, df_agent, df_model, adata, mdata, N, alabels, mlabels
    )
    Agents.collect_agent_data!(df_agent, model, adata, 0)
    Agents.collect_model_data!(df_model, model, mdata, 0)
    La = isnothing(adata) ? 0 : size(df_agent)[2]-1
    Lm = isnothing(mdata) ? 0 : size(df_model)[2]-1
    data, axs = [], []

    # Plot all quantities
    # TODO: make scatter+line plot 1.
    for i in 1:La
        x = Agents.aggname(adata[i])
        val = Observable([df_agent[end, x]])
        push!(data, val)
        ax = datalayout[i, :] = Axis(figure)
        push!(axs, ax)
        ax.ylabel = isnothing(alabels) ? x : alabels[i]
        c = JULIADYNAMICS_COLORS[mod1(i, 3)]
        lines!(ax, N, val, color = c)
        scatter!(
            ax, N, val; marker = MARKER, markersize = 5Makie.px,
            color = c, strokewidth = 0.5
        )
    end
    for i in 1:Lm
        x = Agents.aggname(mdata[i])
        val = Observable([df_model[end, x]])
        push!(data, val)
        ax = datalayout[i+La, :] = Axis(figure)
        push!(axs, ax)
        ax.ylabel = isnothing(mlabels) ? x : mlabels[i]
        c = JULIADYNAMICS_COLORS[mod1(i+La, 3)]
        lines!(ax, N, val, color = c)
        scatter!(ax, N, val, marker = MARKER, markersize = 5Makie.px, color = c,
                 strokewidth = 0.5)
    end
    if La+Lm > 1
        for ax in @view(axs[1:end-1]); hidexdecorations!(ax, grid = false); end
    end
    axs[end].xlabel = "step"
    return data, axs
end

function update_abm_data_plots!(data, axs, model, df_agent, df_model, adata, mdata, N)
    Agents.collect_agent_data!(df_agent, model, adata, N[][end])
    Agents.collect_model_data!(df_model, model, mdata, N[][end])
    La = isnothing(adata) ? 0 : size(df_agent)[2]-1
    Lm = isnothing(mdata) ? 0 : size(df_model)[2]-1

    for i in 1:La
        o = data[i]
        x = Agents.aggname(adata[i])
        val = df_agent[end, x]
        push!(o[], val)
        o[] = o[] #update plot
    end
    for i in 1:Lm
        o = data[i+La]
        x = Agents.aggname(mdata[i])
        val = df_model[end, x]
        push!(o[], val)
        o[] = o[] #update plot
    end
    # TODO: Maybe we can optimize this by setting the limits ourselves (we have the data)
    for ax in axs; autolimits!(ax); end
end

function update_abm_properties!(model, slidervals)
    for l in keys(slidervals)
        v = slidervals[l][]
        if has_key(model.properties, l)
            set_value!(model.properties, l, v)
        else
            throw(KeyError("$l"))
        end
    end
end

function add_reset_line!(axs, s)
    for ax in axs
        # vline!(ax, s; color = "#c41818")
        vlines!(ax, [s], color = "#c41818")
    end
end


text:
@testset "string macro" begin
    @test i"abc" == "abc"
    @test i"abc" == i"abc"

    @test i"a\na\na\na" == join(fill("a", 4), "\n")
    @test addr_eq(i"a\na\na\na", intern(join(fill("a", 4), "\n")))

    @test i"a $(2*54) c" == "a 108 c"
    @test addr_eq(i"a $(2*54) c", i"a 108 c")

    x = "cruel"
    @test i"hello $x world" == "hello cruel world"
    @test addr_eq(i"hello $x world", i"hello cruel world")
end


text:
"""
Add values to extras
"""
function addextras(data, outputtype, filename, path)
    data.extras[ModelSelection.generate_extra_key(OUTPUT_EXTRAKEY, data.extras)] = Dict(
        :outputtype => outputtype,
        :filename => filename,
        :path => path
    )
    return data
end

function get_array_details(arr)
    dict = Dict(arr)
    Dict(
        ((length(arr) == 1) ? "var" : "vars") => Dict(
            "names" => map(string, collect(keys(dict))),
            "values" => collect(values(dict))
        ),
        ((length(arr) == 1) ? "vars" : "var") => false
    )
end

function get_array_simple_details(arr)
    Dict(
        ((length(arr) == 1) ? "var" : "vars") => Dict(
            "names" => map(string, arr)
        ),
        ((length(arr) == 1) ? "vars" : "var") => false
    )
end


text:
# functional derivative calculation

import PDESolver._evalFunctionalDeriv_q

"""
  Computes a 3D array of hte derivative of a functional wrt eqn.q.

  The derivative is evaluated at the state in eqn.q_vec.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * func: AbstractBoundaryFunctional to evaluate

  **Inputs/Outputs**

   * func_deriv_arr: array to hold derivative of function wrt eqn.q, same
                     size as equation.q

  **Options Keys**

  This function is not compatible with `precompute_q_bndry` = false
"""
function _evalFunctionalDeriv_q(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol}, opts,
                           func::AbstractBoundaryFunctional,
                           func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol}

  if mesh.isDG
    boundaryinterpolate!(mesh.sbpface, mesh.bndryfaces, eqn.q, eqn.q_bndry)
  end

  calcFunctionalDeriv(mesh, sbp, eqn, opts, func, func_deriv_arr)

  return nothing
end

@doc """
### EulerEquationMod. calcFunctionalDeriv

Computes a 3D array of the derivative of a functional w.r.t eqn.q on all
mesh nodes.

**Inputs**

*  `mesh` : Abstract DG mesh type
*  `sbp`  : Summation-By-parts operator
*  `eqn`  : Euler equation object
*  `opts` : Options dictionary
*  `func` : Functional object of super-type AbstractBoundaryFunctional
                      that is needed for computing the adjoint vector.
                      Depending on the functional being computed, a different
                      method based on functional type may be needed to be
                      defined.
*  `func_deriv_arr` : 3D array that stores the derivative of the functional
                      w.r.t. eqn.q. The array is the same size as eqn.q

**Outputs**

*  None

"""->
function calcFunctionalDeriv(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol}, opts,
                           func::AbstractBoundaryFunctional,
                           func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol}

  integrand = zeros(eqn.q_bndry)  #TODO: only allocate enough space for
                                  #      the boundaries the functional is using
  node_info = Array{Int}(3)

  # Populate integrand
  for itr = 1:length(func.bcnums)
    bcnum = func.bcnums[itr]

    start_index = mesh.bndry_offsets[bcnum]
    end_index = mesh.bndry_offsets[bcnum+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i

    nfaces = length(bndry_facenums)
#    q2 = zeros(Tsol, mesh.numDofPerNode)
    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      global_facenum = idx_range[i]
      for j = 1:mesh.sbpface.numnodes
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        q = sview(eqn.q_bndry, :, j, global_facenum)
#        convertToConservative(eqn.params, q, q2)
        aux_vars = ro_sview(eqn.aux_vars_bndry, :, j, global_facenum)
        x = ro_sview(mesh.coords_bndry, :, j, global_facenum)
        nrm = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
        node_info[1] = itr; node_info[2] = j; node_info[3] = i
        integrand_i = sview(integrand, :, j, global_facenum)

        calcIntegrandDeriv(opts, eqn.params, q, aux_vars, nrm, integrand_i, node_info,
                           func)
      end  # End for j = 1:mesh.sbpface.numnodes
    end    # End for i = 1:nfaces
  end      # End for itr = 1:length(functional_edges)

  fill!(func_deriv_arr, 0)
  boundaryintegrate!(mesh.sbpface, mesh.bndryfaces, integrand, func_deriv_arr)

  return nothing
end  # End function calcFunctionalDeriv

"""
  Method for LiftCoefficient, which delegates most of the computation to the
  Lift functional.
"""
function calcFunctionalDeriv(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol}, opts,
                           func::AeroCoefficients,
                           func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol}

  calcFunctionalDeriv(mesh, sbp, eqn, opts, func.func, func_deriv_arr)

  Ma = eqn.params.Ma
  fac = 0.5*eqn.params.rho_free*Ma*Ma

  scale!(func_deriv_arr, 1./fac)

  return nothing
end

"""
  Method for SolutionDeviation
"""
function calcFunctionalDeriv(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol}, opts,
                           func::SolutionDeviation,
                           func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol}

  val_bar = 1
  u_bar = zeros(Tsol, mesh.numDofPerNode)
  u_bar_bar = zeros(Tsol, mesh.numDofPerNode)
  for i=1:mesh.numEl

    # compute u_bar
    vol = zero(Tmsh)  # area/volume of the element
    fill!(u_bar, 0)
    for j=1:mesh.numNodesPerElement
      w_j = sbp.w[j]/mesh.jac[j, i]
      vol += w_j
      for k=1:mesh.numDofPerNode
        u_bar[k] += eqn.q[k, j, i]*w_j
      end
    end

    for k=1:mesh.numDofPerNode
      u_bar[k] /= vol
    end
#=
    # compute \int || u - u_bar ||
    for j=1:mesh.numNodesPerElement
      w_j = sbp.w[j]/mesh.jac[j, i]
      for k=1:mesh.numDofPerNode
        delta_u = eqn.q[k, j, i] - u_bar[k]
        val += delta_u*w_j*delta_u
      end
    end
=#
    #----------------------------------
    # reverse sweep
    fill!(u_bar_bar, 0)
    for j=1:mesh.numNodesPerElement
      w_j = sbp.w[j]/mesh.jac[j, i]
      for k=1:mesh.numDofPerNode
        delta_u = eqn.q[k, j, i] - u_bar[k]
        #val += delta_u*w_j*delta_u

        delta_u_bar = 2*delta_u*w_j*val_bar
        func_deriv_arr[k, j, i] += delta_u_bar
        u_bar_bar[k]       -= delta_u_bar
      end
    end

    vol_bar = zero(Tmsh)
    for k=1:mesh.numDofPerNode
      u_bar[k] *= vol  # restore primal value
      u_bar_bar[k] /= vol
    end

    for j=1:mesh.numNodesPerElement
      w_j = sbp.w[j]/mesh.jac[j, i]
      for k=1:mesh.numDofPerNode
        func_deriv_arr[k, j, i] += u_bar_bar[k]*w_j
      end
    end

  end  # end i

  return nothing
end

#------------------------------------------------------------------------------
# generic fallback: complex step it

function calcIntegrandDeriv(opts, params::ParamType{2},
          q::AbstractArray{Tsol,1},
          aux_vars::AbstractArray{Tres, 1},
          nrm::AbstractArray{Tmsh},
          integrand_deriv::AbstractArray{Tsol, 1}, node_info,
          func::Tfunc) where {Tsol, Tres, Tmsh, Tfunc<:AbstractBoundaryFunctional}

  h = 1e-20
  pert = Complex128(0, h)
  for i=1:length(q)
    q[i] += pert
    val = calcBoundaryFunctionalIntegrand(params, q, aux_vars, nrm, node_info,
                                          func)
    integrand_deriv[i] = imag(val)/h
    q[i] -= pert
  end

  return nothing
end


text:
module ModiaResult

const path = dirname(dirname(@__FILE__))

include("AbstractInterface.jl")
include("NoPlot.jl")
include("SilentNoPlot.jl")
include("Utilities.jl")
include("OverloadedMethods.jl")
#include("CompareResults.jl")
include("ResultViews.jl")

end # module

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
"""
# module BisectPy

# Examples

```jldoctest
julia>
```
"""
module BisectPy

using Reexport: @reexport

include("Bisect.jl")
@reexport using .Bisect

include("Insort.jl")
@reexport using .Insort

include("Find.jl")

end # module


text:
using Test, SearchLight, Car

@testset "Car unit tests" begin

  ### Your tests here
  @test 1 == 1

end;


text:
@testset "ilds.jl" begin
    @testset "expandIlds!()" begin
        ### Checking status ###
        # :NodeLimitStop
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.numberOfNodes = 1
        toCall = Stack{Function}()
        @test SeaPearl.expandIlds!(toCall,0, model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(),nothing) == :NodeLimitStop
        @test isempty(toCall)

        # :SolutionLimitStop
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.numberOfSolutions = 0

        toCall = Stack{Function}()
        @test SeaPearl.expandIlds!(toCall,0, model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(),nothing) == :SolutionLimitStop
        @test isempty(toCall)

        #:TimeLimitStop
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.searchingTime = 0

        toCall = Stack{Function}()
        SeaPearl.tic()
        @test SeaPearl.expandIlds!(toCall,0, model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(),nothing) == :TimeLimitStop
        @test isempty(toCall)

        # :Infeasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(3, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.expandIlds!(toCall,0, model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(),nothing) == :Infeasible
        @test isempty(toCall)

        # :Feasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 2, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.expandIlds!(toCall,0, model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(),nothing) == :FoundSolution
        @test isempty(toCall)

        ### Checking stack ###
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.expandIlds!(toCall,0, model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(),nothing) == :Feasible
        @test length(toCall) == 3

        @test pop!(toCall)(model) == :SavingState
        @test length(model.trailer.prior) == 1 # saveState!()

        @test pop!(toCall)(model) == :FoundSolution
        @test length(model.statistics.solutions) == 1 # Found a solution

        @test pop!(toCall)(model) == :BackTracking
        @test length(model.trailer.prior) == 0 # restoreState!()
    end

    @testset "initroot(::ILDSearch)" begin
        # :Feasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 2, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.initroot!(toCall, SeaPearl.ILDSearch(0), model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :FoundSolution
        @test isempty(toCall)

    end

    @testset "search!(::ILDSearch)" begin
        ### Checking status ###
        # :LimitStop
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.numberOfNodes = 1
        x = SeaPearl.IntVar(1, 1, "x", trailer)
        SeaPearl.addVariable!(model, x)

        @test SeaPearl.search!(model, SeaPearl.ILDSearch(0), SeaPearl.MinDomainVariableSelection()) == :NodeLimitStop

        # :SolutionLimitStop
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.numberOfSolutions = 0

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        @test SeaPearl.search!(model, SeaPearl.ILDSearch(0), SeaPearl.MinDomainVariableSelection()) == :SolutionLimitStop

        # :Infeasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(3, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        @test SeaPearl.search!(model, SeaPearl.ILDSearch(0), SeaPearl.MinDomainVariableSelection()) == :Infeasible

        # :Optimal
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 2, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        @test SeaPearl.search!(model, SeaPearl.ILDSearch(0), SeaPearl.MinDomainVariableSelection()) == :NonOptimal
        @test length(model.statistics.solutions) == 1

        ### Checking more complex statistics.solutions ###
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        @test SeaPearl.search!(model, SeaPearl.ILDSearch(1), SeaPearl.MinDomainVariableSelection()) == :NonOptimal
        @test length(model.statistics.solutions) == 2
        @test model.statistics.solutions[2] == Dict("x" => 2,"y" => 2)
        @test model.statistics.solutions[1] == Dict("x" => 3,"y" => 3)

    end

    @testset "search!() with a BasicHeuristic" begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        @test SeaPearl.search!(model, SeaPearl.ILDSearch(1), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :NonOptimal
        @test model.statistics.solutions[2] == Dict("x" => 2,"y" => 2)
        @test model.statistics.solutions[1] == Dict("x" => 3,"y" => 3)

        SeaPearl.empty!(model)

        x = SeaPearl.IntVar(2, 3, "x", model.trailer)
        y = SeaPearl.IntVar(2, 3, "y", model.trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, model.trailer))

        my_heuristic(x::SeaPearl.IntVar; cpmodel=nothing) = minimum(x.domain)
        @test SeaPearl.search!(model, SeaPearl.ILDSearch(1), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(my_heuristic)) == :NonOptimal
        @test model.statistics.solutions[2] == Dict("x" => 3,"y" => 3)

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x1 = SeaPearl.IntVar(1, 2, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 3, "x2", trailer)
        x3 = SeaPearl.IntVar(2, 4, "x3", trailer)
        SeaPearl.addVariable!(model, x1)
        SeaPearl.addVariable!(model, x2)
        SeaPearl.addVariable!(model, x3)

        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x1, x2, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x2, x3, trailer))
        @test SeaPearl.search!(model, SeaPearl.ILDSearch(3), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :NonOptimal

        #in this specific finding order with the given heuristic and iLDS method.
        @test model.statistics.solutions[1] == Dict("x1" => 2,"x2" => 3,"x3" => 4)
        @test model.statistics.solutions[2] == Dict("x1" => 2,"x2" => 3,"x3" => 2)
        @test model.statistics.solutions[3] == Dict("x1" => 2,"x2" => 1,"x3" => 4)
        @test model.statistics.solutions[4] == Dict("x1" => 1,"x2" => 3,"x3" => 4)
        @test model.statistics.solutions[5] == Dict("x1" => 2,"x2" => 1,"x3" => 3)
        @test model.statistics.solutions[6] == Dict("x1" => 1,"x2" => 3,"x3" => 2)
        @test model.statistics.solutions[7] == Dict("x1" => 1,"x2" => 2,"x3" => 4)
        @test model.statistics.solutions[8] == Dict("x1" => 2,"x2" => 1,"x3" => 2)
        @test model.statistics.solutions[9] == Dict("x1" => 1,"x2" => 2,"x3" => 3)

    end
    @testset "search!() with a BasicHeuristic - out_solver being true" begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        @test SeaPearl.search!(model, SeaPearl.ILDSearch(1), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(); out_solver=true) == :FoundSolution
        @test length(model.statistics.solutions) == 1
        @test model.statistics.solutions[1] == Dict("x" => 3,"y" => 3)

        SeaPearl.empty!(model)

        x = SeaPearl.IntVar(2, 3, "x", model.trailer)
        y = SeaPearl.IntVar(2, 3, "y", model.trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, model.trailer))

        my_heuristic(x::SeaPearl.IntVar; cpmodel=nothing) = minimum(x.domain)
        @test SeaPearl.search!(model, SeaPearl.ILDSearch(1), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(my_heuristic); out_solver=true) == :FoundSolution
        @test length(model.statistics.solutions) == 1
        @test model.statistics.solutions[1] == Dict("x" => 2,"y" => 2)

    end

    @testset "search!() with a LearnedHeuristic I" begin

        approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
        target_approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
        gnnlayers = 1
        approximator_model = SeaPearl.CPNN(
            graphChain = Flux.Chain(
                SeaPearl.GraphConv(3 => 64, Flux.leakyrelu),
                [approximator_GNN for i = 1:gnnlayers]...
            ),
            nodeChain = Flux.Chain(
                Flux.Dense(64, 32, Flux.leakyrelu),
                Flux.Dense(32, 32, Flux.leakyrelu),
                Flux.Dense(32, 16, Flux.leakyrelu),
            ),
            outputChain = Flux.Dense(16, 4),
        )
        target_approximator_model = SeaPearl.CPNN(
            graphChain = Flux.Chain(
                SeaPearl.GraphConv(3 => 64, Flux.leakyrelu),
                [target_approximator_GNN for i = 1:gnnlayers]...
            ),
            nodeChain = Flux.Chain(
                Flux.Dense(64, 32, Flux.leakyrelu),
                Flux.Dense(32, 32, Flux.leakyrelu),
                Flux.Dense(32, 16, Flux.leakyrelu),
            ),
            outputChain = Flux.Dense(16, 4),
        )

        agent = RL.Agent(
            policy = RL.QBasedPolicy(
                learner = RL.DQNLearner(
                    approximator = RL.NeuralNetworkApproximator(
                        model = approximator_model,
                        optimizer = ADAM(0.001f0)
                    ),
                    target_approximator = RL.NeuralNetworkApproximator(
                        model = target_approximator_model,
                        optimizer = ADAM(0.001f0)
                    ),
                    loss_func = Flux.Losses.huber_loss,
                    stack_size = nothing,
                    γ = 0.99f0,
                    batch_size = 32,
                    update_horizon = 1,
                    min_replay_history = 1,
                    update_freq = 1,
                    target_update_freq = 100,
                ),
                explorer = RL.EpsilonGreedyExplorer(
                    ϵ_stable = 0.01,
                    kind = :exp,
                    ϵ_init = 1.0,
                    warmup_steps = 0,
                    decay_steps = 500,
                    step = 1,
                    is_break_tie = false,
                    #is_training = true,
                    rng = MersenneTwister(33)
                )
            ),
            trajectory = RL.CircularArraySLARTTrajectory(
                capacity = 500,
                state = SeaPearl.DefaultTrajectoryState[] => (),
                action = Int => (),
                legal_actions_mask = Vector{Bool} => (4, ),
            )
        )

        # define the value selection
        valueSelection = SeaPearl.LearnedHeuristic(agent)
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x1 = SeaPearl.IntVar(2, 3, "x1", trailer)
        x2 = SeaPearl.IntVar(2, 3, "x2", trailer)
        x3 = SeaPearl.IntVar(3, 4, "x3", trailer)
        x4 = SeaPearl.IntVar(2, 5, "x4", trailer)
        SeaPearl.addVariable!(model, x1)
        SeaPearl.addVariable!(model, x2)
        SeaPearl.addVariable!(model, x3)
        SeaPearl.addVariable!(model, x4)

        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x1, x2, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x2, x3, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x3, x4, trailer))

        # define the variable selection
        variableSelection = SeaPearl.MinDomainVariableSelection()

        # launch the search
        SeaPearl.search!(model, SeaPearl.ILDSearch(6), variableSelection, valueSelection)

        possible_solutions = [
            Dict("x1" => 2, "x2" => 3, "x3" => 4, "x4" => 2),
            Dict("x1" => 2, "x2" => 3, "x3" => 4, "x4" => 3),
            Dict("x1" => 2, "x2" => 3, "x3" => 4, "x4" => 5),
            Dict("x1" => 3, "x2" => 2, "x3" => 4, "x4" => 2),
            Dict("x1" => 3, "x2" => 2, "x3" => 4, "x4" => 3),
            Dict("x1" => 3, "x2" => 2, "x3" => 4, "x4" => 5),
            Dict("x1" => 3, "x2" => 2, "x3" => 3, "x4" => 2),
            Dict("x1" => 3, "x2" => 2, "x3" => 3, "x4" => 4),
            Dict("x1" => 3, "x2" => 2, "x3" => 3, "x4" => 5)
        ]

        for solution in model.statistics.solutions
            @test solution in possible_solutions
        end

    end

    @testset "search!() with a LearnedHeuristic out of the solver" begin

        approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
        target_approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
        gnnlayers = 1
        approximator_model = SeaPearl.CPNN(
            graphChain = Flux.Chain(
                SeaPearl.GraphConv(3 => 64, Flux.leakyrelu),
                [approximator_GNN for i = 1:gnnlayers]...
            ),
            nodeChain = Flux.Chain(
                Flux.Dense(64, 32, Flux.leakyrelu),
                Flux.Dense(32, 32, Flux.leakyrelu),
                Flux.Dense(32, 16, Flux.leakyrelu),
            ),
            outputChain = Flux.Dense(16, 4),
        )
        target_approximator_model = SeaPearl.CPNN(
            graphChain = Flux.Chain(
                SeaPearl.GraphConv(3 => 64, Flux.leakyrelu),
                [target_approximator_GNN for i = 1:gnnlayers]...
            ),
            nodeChain = Flux.Chain(
                Flux.Dense(64, 32, Flux.leakyrelu),
                Flux.Dense(32, 32, Flux.leakyrelu),
                Flux.Dense(32, 16, Flux.leakyrelu),
            ),
            outputChain = Flux.Dense(16, 4),
        )

        agent = RL.Agent(
            policy = RL.QBasedPolicy(
                learner = RL.DQNLearner(
                    approximator = RL.NeuralNetworkApproximator(
                        model = approximator_model,
                        optimizer = ADAM(0.001f0)
                    ),
                    target_approximator = RL.NeuralNetworkApproximator(
                        model = target_approximator_model,
                        optimizer = ADAM(0.001f0)
                    ),
                    loss_func = Flux.Losses.huber_loss,
                    stack_size = nothing,
                    γ = 0.99f0,
                    batch_size = 32,
                    update_horizon = 1,
                    min_replay_history = 1,
                    update_freq = 1,
                    target_update_freq = 100,
                ),
                explorer = RL.EpsilonGreedyExplorer(
                    ϵ_stable = 0.01,
                    kind = :exp,
                    ϵ_init = 1.0,
                    warmup_steps = 0,
                    decay_steps = 500,
                    step = 1,
                    is_break_tie = false,
                    #is_training = true,
                    rng = MersenneTwister(33)
                )
            ),
            trajectory = RL.CircularArraySLARTTrajectory(
                capacity = 500,
                state = SeaPearl.DefaultTrajectoryState[] => (),
                action = Int => (),
                legal_actions_mask = Vector{Bool} => (4, ),
            )
        )

        # define the value selection
        valueSelection = SeaPearl.LearnedHeuristic(agent)

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x1 = SeaPearl.IntVar(1, 2, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 2, "x2", trailer)
        x3 = SeaPearl.IntVar(2, 3, "x3", trailer)
        x4 = SeaPearl.IntVar(1, 4, "x4", trailer)
        SeaPearl.addVariable!(model, x1)
        SeaPearl.addVariable!(model, x2)
        SeaPearl.addVariable!(model, x3)
        SeaPearl.addVariable!(model, x4)

        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x1, x2, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x2, x3, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x3, x4, trailer))

        # define the variable selection
        variableSelection = SeaPearl.MinDomainVariableSelection()

        # launch the search
        SeaPearl.search!(model, SeaPearl.ILDSearch(6), variableSelection, valueSelection; out_solver=true)

        possible_solutions = [
            Dict("x1" => 1, "x2" => 2, "x3" => 3, "x4" => 1),
            Dict("x1" => 1, "x2" => 2, "x3" => 3, "x4" => 2),
            Dict("x1" => 1, "x2" => 2, "x3" => 3, "x4" => 4),
            Dict("x1" => 2, "x2" => 1, "x3" => 3, "x4" => 1),
            Dict("x1" => 2, "x2" => 1, "x3" => 3, "x4" => 2),
            Dict("x1" => 2, "x2" => 1, "x3" => 3, "x4" => 4),
            Dict("x1" => 2, "x2" => 1, "x3" => 2, "x4" => 1),
            Dict("x1" => 2, "x2" => 1, "x3" => 2, "x4" => 3),
            Dict("x1" => 2, "x2" => 1, "x3" => 2, "x4" => 4)
        ]

        for solution in model.statistics.solutions
            @test solution in possible_solutions
        end

    end

end


text:
###
### logistic and logit
###

logistic(x::Real) = inv(one(x) + exp(-x))

function logistic_logjac(x::Real)
    mx = -abs(x)
    mx - 2*log1p(exp(mx))
end

logit(x::Real) = log(x / (one(x) - x))

logit_logjac(y) = -log(y) - log1p(-y)

###
### calculations
###

"""
    $SIGNATURES

Number of elements (strictly) above the diagonal in an ``n×n`` matrix.
"""
unit_triangular_dimension(n::Int) = n * (n-1) ÷ 2

###
### type calculations
###

"""
    $(SIGNATURES)

Extend element type of argument so that it is closed under the algebra used by this package.

Pessimistic default for non-real types.
"""
function extended_eltype(::Type{S}) where S
    T = eltype(S)
    T <: Real ? typeof(√(one(T))) : Any
end

extended_eltype(x::T) where T = extended_eltype(T)

####
#### random values
####

"Shared part of docstrings for keyword arguments of or passed to [`random_reals`](@ref)."
const _RANDOM_REALS_KWARGS_DOC = """
A standard multivaritate normal or Cauchy is used, depending on `cauchy`, then scaled with
`scale`. `rng` is the random number generator used.
"""

_random_reals_scale(rng::AbstractRNG, scale::Real, cauchy::Bool) =
    cauchy ? scale / abs2(randn(rng)) : scale * 1.0

"""
$(SIGNATURES)

Random real number.

$(_RANDOM_REALS_KWARGS_DOC)
"""
random_real(; scale::Real = 1, cauchy::Bool = false, rng::AbstractRNG = GLOBAL_RNG) =
    randn(rng) * _random_reals_scale(rng, scale, cauchy)

"""
$(SIGNATURES)

Random vector in ``ℝⁿ``.

$(_RANDOM_REALS_KWARGS_DOC)
"""
random_reals(n::Integer; scale::Real = 1, cauchy::Bool = false,
             rng::AbstractRNG = GLOBAL_RNG) =
                 randn(rng, n) .* _random_reals_scale(rng, scale, cauchy)


text:
"""
    eigDAG_Distance(𝚽, Q, numEigs; edge_weights = 1)

compute DAG distances between pairwise graph Laplacian eigenvectors.

# Input Arguments
- `𝚽::Matrix{Float64}`: matrix of graph Laplacian eigenvectors, 𝜙ⱼ₋₁ (j = 1,...,size(𝚽,1)).
- `Q::Matrix{Float64}`: incidence matrix of the graph.
- `numEigs::Int64`: number of eigenvectors considered.
- `edge_weight::Any`: default value is 1, stands for unweighted graph
    (i.e., all edge weights equal to 1). For weighted graph, edge_weight is the
    weights vector, which stores the affinity weight of each edge.

# Output Argument
- `dis::Matrix{Float64}`: a numEigs x numEigs distance matrix, dis[i,j] = d_DAG(𝜙ᵢ₋₁, 𝜙ⱼ₋₁).

"""
function eigDAG_Distance(𝚽, Q, numEigs; edge_weight = 1)
    dis = zeros(numEigs, numEigs)
    abs_∇𝚽 = abs.(Q' * 𝚽)
    for i = 1:numEigs, j = i+1:numEigs
        dis[i,j] = (edge_weight == 1) ? norm(abs_∇𝚽[:,i]-abs_∇𝚽[:,j],2) : norm((abs_∇𝚽[:,i]-abs_∇𝚽[:,j]).*sqrt.(edge_weight),2)
    end
    return dis + dis'
end

function eigDAG_Distance_normalized(𝚽,Q,numEigs; edge_weight = 1)
    dis = zeros(numEigs,numEigs)
    abs_∇𝚽 = abs.(Q' * 𝚽)
    for i = 1:numEigs, j = i+1:numEigs
        dis[i,j] = (edge_weight == 1) ? norm(abs_∇𝚽[:,i]-abs_∇𝚽[:,j],2) : norm((abs_∇𝚽[:,i]-abs_∇𝚽[:,j]).*sqrt.(edge_weight),2)
        dis[i,j] /= norm(𝚽[:,i] .* 𝚽[:,j],2)
    end
    return dis + dis'
end