text:
using Catlab.CategoricalAlgebra
using Catlab.Present
using Catlab.Theories
using Catlab.WiringDiagrams
using Catlab.Programs
using Catlab.Programs.RelationalPrograms: parse_relation_context, parse_relation_call

using Catlab.CategoricalAlgebra.FinSets
using Catlab.CategoricalAlgebra.CSetDataStructures: AttributedCSet
using DataStructures: IntDisjointSets, find_root!, in_same_set
using Dates

struct Schema
  arities::Dict{Symbol, Int}
end

ids = x -> Symbol("id_$x")
arg = (x, y) -> Symbol("$(x)_$y")

"""
We can slightly complicate this by splitting V up into an object for each color
and arities can be a list of symbols.
"""
function schema_to_csettype(sch::Schema)
  pres = Presentation(FreeSchema)

  Id, Name = [FreeSchema.Data{:generator}([x], []) for x in [:Id, :Name]]
  add_generator!(pres, Id)
  add_generator!(pres, Name)

  relobs = [Ob(FreeSchema, s) for s in keys(sch.arities)]
  vob = Ob(FreeSchema, :V)

  add_generator!(pres, vob)
  add_generator!(pres, FreeSchema.Attr{:generator}([:id_V, vob, Id], [vob, Id]))
  add_generator!(pres, FreeSchema.Attr{:generator}([:sort, vob, Name], [vob, Name]))

  for x in relobs
    add_generator!(pres, x)
    add_generator!(pres, FreeSchema.Attr{:generator}(
      [ids(x.args[1]), x, Id], [x, Id]))
  end

  for (n, (s, a)) in enumerate(sch.arities)
    relob = relobs[n]
    for src_sym in 1:a
      add_generator!(pres, Hom(arg(s,src_sym), relob, vob))
    end
  end

  cset = ACSetType(pres, unique_index=vcat([:id_V],
                   map(ids, collect(keys(sch.arities)))))
  return ACSetType(pres){Int, Symbol}, cset{Int, Symbol}
end

"""
Dependency: When the body pattern is seen, the result pattern is asserted to
exist and eqs hold between subsets of the body variables.
"""
struct Dep
  name :: Symbol
  body :: ACSet
  res :: ACSet
  eqs :: Set{Set{Int}}
end

function Dep(name::Symbol, body::ACSet,res::ACSet)::Dep
  return Dep(name, body,res,Set{Set{Int}}())
end

mutable struct Sequence{State,Params}
  state::State
  params::Params
  func::Function
end

"""
Give a sequence of rules to run. The last iteration of rules will be repeated.
"""
function seq(s::Vector{Vector{Symbol}})::Function
  !isempty(s) || error("s cannot be empty")
  function f(Σ, i, _, _, _)
    ds = s[min(i, end)]
    return isempty(ds) ? Σ : [getdep(Σ, d) for d in ds]
  end
  return f
end

allseq = seq([Symbol[]])

function sizes(x::ACSet{CD}) where {CD}
  [o=>nparts(x,o) for o in ob(CD)]
end

getdep(Σ::Set{Dep}, name::Symbol)::Dep = only([d for d in Σ if d.name==name])

function maxes(x::ACSet{CD})::Dict{Symbol, Int} where {CD}
  Dict([o => max0(x[ids(o)]) for o in ob(CD)])
end

function merge_maxes(d1::Dict{Symbol, Int},
                     d2::Dict{Symbol, Int})::Dict{Symbol, Int}
  Dict([o => max(d1[o], d2[o]) for o in keys(d1)])
end

Base.isempty(dI::ACSet) = sum(nparts(dI, o)
  for o in filter(!=(:V), ob(typeof(dI).parameters[1]))) == 0

"""Adds content of y to x"""
function addstuff!(X::ACSet{CD}, Y::ACSet{CD})::Nothing where {CD}
  xids = Set(X[:id_V])
  xmap = Int[]
  for (vi, vsrt) in Y.tables[:V]
    if vi ∈ xids
      nxt = only(incident(X, vi, :id_V))
    else
      nxt = add_part!(X, :V, id_V=vi, sort=vsrt)
    end
    push!(xmap, nxt)
  end
  mxs= maxes(X)
  for o in filter(!=(:V), ob(CD))
    oseen = Set([collect(r)[1:end-1] for r in X.tables[o]])
    for r in Y.tables[o]
      cr = collect(r)
      newx = xmap[cr[1:end-1]]
      if newx ∉ oseen
        push!(oseen, newx)
        push!(newx, cr[end]+mxs[o])
        add_part!(X, o; Dict(zip(keys(r), newx))...)
      end
    end
  end
  return nothing
end

function Base.union(x::ACSet{CD}, y::ACSet{CD})::ACSet{CD} where {CD}
  if nparts(x, :V) == 0 return y
  elseif nparts(y,:V) ==0 return x
  end

  overlap = AttributedCSet{typeof(x).parameters...}()
  o_ids = Dict([let s = ids(o); o=>collect(x[s] ∩ y[s]) end
                for o in ob(CD)])
  o_indexes = filter(x->!isempty(x[2]), [o => vcat(incident(x, v, ids(o))...)
               for (o, v) in collect(o_ids)])
  if !isempty(o_indexes)
    copy_parts!(overlap, x; Dict(o_indexes)...)
  end
  inexact = collect(filter(!=(:id_V), map(ids, ob(CD))))
  fx = homomorphism(overlap,x)
  fy = homomorphism(overlap,y, inexact=inexact)
  !(fx===nothing) || error("overlap $overlap x $x")
  !(fy===nothing) || error("overlap $overlap y $y")
  res = apex(pushout(fx,fy))
  trim!(res)
  return res
end

function unionfind_to_dict(d::IntDisjointSets, idvec::Vector{Int})::Dict{Int,Int}
  dic, mins = Dict{Int, Int}(), Set{Int}()
  for i in 1:length(d)
    iid = idvec[i]
    for m in mins
      if in_same_set(d, i, m)
        dic[iid] = idvec[m]
        break
      end
    end
    if !haskey(dic, iid)
      push!(mins, i)
    end
  end
  return dic
end
"""
Apply a mapping to a DB instance. The mapping refers to V id's. {a ↦ b}
collapses the vertex identified by `a` into the vertex identified by `b`.
If `b` is not found in the DB instance, a vertex gets created for it.
"""
function appdict(d::Dict{Int,Int}, I::ACSet{CD})::ACSet{CD} where {CD}
  I_ = deepcopy(I)
  seenV = Set{Int}()
  for o in filter(!=(:V), ob(CD))
    del_list, seen = Int[], Set()
    for (i, r) in enumerate(I.tables[o])
      for (k, v) in filter(x->x[1]!=ids(o), collect(zip(keys(r), r)))
        v_id = I_[:id_V][v]
        new_v_id = get(d, v_id, v_id)
        new_v_ = incident(I_, new_v_id, :id_V)
        new_v = (isempty(new_v_)
                ? add_part!(I_, :V, id_V=new_v_id, sort=I_[:sort][v])
                : only(new_v_))
        push!(seenV, new_v)
        set_subpart!(I_, i, k, new_v)
      end
      new_r = collect(I_.tables[o][i])[1:end-1]
      if new_r in seen
        push!(del_list, i)
      else
        push!(seen, new_r)
      end
    end
    rem_parts!(I_, o, del_list)
  end
  rem_parts!(I_, :V, sort(collect(setdiff(parts(I_, :V), seenV))))
  return I_
end

function appdict!(d::Dict{Int,Int}, I::ACSet{CD})::Nothing where {CD}
  seenV = Set{Int}()
  for o in filter(!=(:V), ob(CD))
    del_list, seen = Int[], Set()
    for (i, r) in enumerate(I.tables[o])
      for (k, v) in filter(x->x[1]!=ids(o), collect(zip(keys(r), r)))
        v_id = I[:id_V][v]
        new_v_id = get(d, v_id, v_id)
        new_v_ = incident(I, new_v_id, :id_V)
        new_v = (isempty(new_v_)
                ? add_part!(I, :V, id_V=new_v_id, sort=I[:sort][v])
                : only(new_v_))
        push!(seenV, new_v)
        set_subpart!(I, i, k, new_v)
      end
      new_r = collect(I.tables[o][i])[1:end-1]
      if new_r in seen
        push!(del_list, i)
      else
        push!(seen, new_r)
      end
    end
    rem_parts!(I, o, del_list)
  end
  rem_parts!(I, :V, sort(collect(setdiff(parts(I, :V), seenV))))
  return nothing
end

appdict(d::Dict{Int,Int}, v::Vector{Int})::Vector{Int} = [get(d, x, x) for x in v]

max0 = x -> isempty(x) ? 0 : maximum(x)

"""Offset IDs componentwise"""
function offset(I::ACSet{CD}, maxes::Dict{Symbol, Int})::ACSet{CD} where {CD}
  res = deepcopy(I)
  for o in ob(CD)
    atr, off = ids(o), maxes[o]
    for i in sort(res[atr], rev=true)
      set_subpart!(res, only(incident(res, i, atr)), atr, i + off)
    end
  end
  return res
end

"""Creates a new version of add_on with fresh ids
I - a base that will have y added to it
mapping - a subset of node IDs in add_on are mapped to IDs in I
add_on - a pattern that will be added to I
"""
function fresh(mapping::Dict{Int,Int},
               add_on::ACSet{CD}, mxs::Dict{Symbol, Int})::ACSet where {CD}
  new_add_on1 = offset(add_on, mxs)
  μ = Dict([k+mxs[:V]=>v for (k, v) in pairs(mapping)])
  new_add_on2 = appdict(μ, new_add_on1)
  trim!(new_add_on2)
  return new_add_on2
end

"""Remove nodes that appear in no tuples. Remove duplicate tuples."""
function trim!(I::ACSet{CD})::Nothing where {CD}
  keep_v = Set{Int}()
  for o in filter(!=(:V), ob(CD))
    del_list, seen = Int[], Set()
    for (i, r) in enumerate(I.tables[o])
      vs = collect(r)[1:end-1]
      if vs in seen
        push!(del_list, i)
      else
        push!(seen, vs)
        union!(keep_v, vs)
      end
    end
    rem_parts!(I, o, del_list)
  end
  del_v = sort(collect(setdiff(parts(I, :V), keep_v)))
  rem_parts!(I, :V, del_v)
  return nothing
end

# function Base.iterate(IΣ::restricted_chase)::Tuple{Tuple{ACSet, Dict{Int,Int}}, chase_state}
#   iterate(IΣ, chase_state(IΣ.I, IΣ.I, 1))
# end

# function Base.iterate(IΣ::core_chase)::Tuple{Tuple{ACSet, Dict{Int,Int}}, core_chase_state}
#   iterate(IΣ, core_chase_state(IΣ.I, 1, Dict([k.name=>Dict([o=>Set{Int}()
#      for o in IΣ.It.body.body.parameters[1].parameters[1]]) for k in IΣ.Σ]), false))
# end

"""
Terminate determines how many chase iterations to do.
Sequence determines what dependencies to trigger each iteration
"""
function core_chase(I::ACSet{CD}, Σ::Set{Dep}, Select::Function=allseq,
                    Terminate::Function=(_,_,_,_)->false, maxi::Int=3) where {CD}
  m, seen = Dict{Int,Int}(), Dict{Symbol, Dict{Symbol, Set{Int}}}([
    d.name=>Dict([o=>Set{Int}() for o in ob(CD)]) for d in Σ])
  for i in 1:maxi
    println("--------Iteration $i--------")
    μ, N = IntDisjointSets(nparts(I, :V)), ACSet{typeof(I).parameters...}()
    for τ in Select(Σ, i, I, N, μ)
      core_apply_trigger!(I, τ, seen[τ.name], N, μ)
    end
    m = unionfind_to_dict(μ, I[:id_V])
    appdict!(m, I)
    appdict!(m, N)
    addstuff!(I, N)
    # I = findcore(I; loose=true)
    if Terminate(I, N, m, i)
      return (I, m)
    end
  end
  return (I, m)
end

# function Base.iterate(IΣ::restricted_chase, state::Union{Nothing,chase_state}
#                       )::Union{Nothing, Tuple{Tuple{ACSet, Dict{Int,Int}}, chase_state}}
#   if isempty(state.dI) # Termination condition
#     return nothing
#   else # Normal unpacking of state
#     I, ∆I, counter, s = state.I, state.dI, state.counter, IΣ.S
#     N, μ = IΣ.It() , IntDisjointSets(nparts(I, :V))
#   end
#   ∆I = I # until we figure out how to do triggerwise-∆I
#   println("counter $counter")
#   while true
#     τ = select(s, IΣ.Σ, counter, I, ∆I, N, μ)
#     if τ===nothing break end
#     N, μ =  apply_trigger(I, τ, ∆I, N, μ)
#   end
#   m = unionfind_to_dict(μ, I[:id_V])
#   I, ∆I = update_inst(I, ∆I, N, m)
#   return (I, m), chase_state(I, ∆I, counter+1) # output, new_state
# end
# Base.IteratorSize(::restricted_chase) = Base.IsInfinite()
# Base.IteratorSize(::core_chase) = Base.IsInfinite()

apply_trigger(I::ACSet, τ::Dep)=apply_trigger(
  I, τ, I, ACSet{typeof(I).parameters...}(), IntDisjointSets(nparts(I, :V)))

function update_inst(I::ACSet, ∆I::ACSet, N::ACSet, μ::Dict{Int,Int})
  ∆I = appdict(μ, N) # TODO understand what line 14 means
  # prune!(∆I, setdiff(Set(N[:id_V]), Set(I[:id_V])))
  I = union(appdict(μ, I), ∆I)
  return I, ∆I
end

# function apply_trigger(I::ACSet, τ::Dep, ∆I::ACSet, N::ACSet, μ::IntDisjointSets)
#   print("applying $(τ.name)")
#   q, ps = tgd_relation(τ.body,  Dict{Int,Int}())
#   trigs = (isempty(τ.body) ? [NamedTuple()] : query(I, q, ps))
#   println("$(length(trigs)) potential triggers")
#   for (i, trig) in enumerate(trigs)
#     print("$i ")
#     #println("potential trigger in I $trig")
#     # make sure trigger at least somewhat pertains to the 'unknown' subset of I
#     d = zip([parse(Int, string(k)[2:end]) for k in keys(trig)], trig)
#     d2 = Dict{Int,Int}([k=>I[:id_V][v] for (k,v) in d])

#     # ACTUALLY we should check that any of the relations intersect,
#     # not the vertices (?) Would require change to output of query. TODO
#     if isempty(τ.body) || !isempty(intersect(values(d2), ∆I[:id_V]))
#       # print("[relevant trigger]")
#       # Quotient by equivalences
#       for eqset in τ.eqs
#         eqids = [only(incident(I, d2[x], :id_V)) for x in eqset]
#         # println("unioning $eqids")
#         for (x, y) in zip(eqids, eqids[2:end])
#           union!(μ, x, y)
#         end
#       end

#       # Add to graph if needed
#       if nparts(τ.res, :V) > 0
#         q, ps = tgd_relation(τ.res, d2)
#         uI = appdict(unionfind_to_dict(μ, I[:id_V]), I ∪ N)
#         if isempty(query(uI, q, ps))  # THIS SHOULD JUST BE A IS_HOMOMORPHIC CHECK
#           # println("active tgd trigger w/ $d2")
#           newstuff = fresh(d2, τ.res, merge_maxes(maxes(N), maxes(I)))
#           # println("newstuff $newstuff")
#           N = N ∪ newstuff
#         end
#       end
#     end
#   end
#   println("")
#   return N, μ
# end

function core_apply_trigger!(I::ACSet{CD}, τ::Dep, seen::Dict{Symbol, Set{Int}},
                             N::ACSet, μ::IntDisjointSets)::Nothing where {CD}
  print("applying $(τ.name) on I = $(sizes(I)) and N = $(sizes(N))")
  strt = now()
  q, ps = tgd_relation(τ.body, Dict{Int,Int}())
  trigs = (isempty(τ.body) ? [NamedTuple()] : query(I, q, ps))
  println("$(length(trigs)) potential triggers (query time $(now()-strt))")
  strt = now()
  for (i, trig) in enumerate(trigs)
    #print("$i ")
    # make sure trigger at least somewhat pertains to the 'unknown' subset of I
    unseen = false
    d2 = Dict{Int,Int}()

    for (k, v) in zip(keys(trig), trig)
      ostr, istr = split(string(k), "__")
      o, i = Symbol(ostr), parse(Int, istr)
      idv = I[ids(o)][v]
      if o == :V
        d2[i] = I[:id_V][v]
        if !unseen
          for vid in I[:id_V]
            if in_same_set(μ, idv, vid) && idv ∉ seen[:V]
              unseen = true
              break
            end
          end
        end
      elseif !unseen && idv ∉ seen[o]
        unseen = true
      end
    end
    if isempty(τ.body) || unseen
      # print(" [relevant] ")
      # Quotient by equivalences
      for eqset in τ.eqs
        eqids = [only(incident(I, d2[x], :id_V)) for x in eqset]
        # print(" unioning $eqids ")
        for (x, y) in zip(eqids, eqids[2:end])
          union!(μ, x, y)
        end
      end
      #ttest = appdict(d2, τ.res)
      resids = Set(τ.res[:id_V])
      init = (V=Dict([(first(incident(τ.res, k, :id_V)) =>
                       first(incident(I, v, :id_V)))
                      for (k,v) in pairs(d2) if k ∈ resids]),)
      inex = collect(map(ids, ob(CD)))
      if !is_homomorphic(τ.res, I, initial=init, inexact=inex)
        # Add to graph if needed
        if nparts(τ.res, :V) > 0
          newstuff = fresh(d2, τ.res, merge_maxes(maxes(N), maxes(I)))
          # println("new V ids $(sort(collect(setdiff(newstuff[:id_V], I[:id_V]))))")
          addstuff!(N, newstuff)
        end
      end
    end
  end
  println("time: $(now()-strt)")
  return nothing
end

"""
Given a DB inst as a pattern, create a conjunctive query that returns its
matches. Optionally constrain particular vertices to have particular IDs.
"""
function tgd_relation(X::ACSet{CD}, bound::Dict{Int,Int})::Tuple{UndirectedWiringDiagram, NamedTuple} where {CD}
  vars = Dict([o=>[Symbol("$(o)__$i") for i in X[ids(o)]] for o in ob(CD)])
  typs = ["V(_id=$i, sort=s$i)" for i in vars[:V]]
  params = Set(["s$i=:$n" for (i, n) in zip(vars[:V], X[:sort])])
  bodstr = vcat(["begin"], typs)
  for (k, v) in collect(bound)
    bound_ind = findfirst(==(k), X[:id_V])
    if !(bound_ind===nothing)
      push!(bodstr, "V(_id=$(vars[bound_ind]),id_V=i$v)")
      push!(params, ("i$v=$v"))
    end
  end
  for o in filter(!=(:V), ob(CD))
    for row in X.tables[o]
      is = collect(enumerate(collect(row)))
      _, oid = pop!(is)
      push!(bodstr, "$o(_id=$(o)__$oid,$(join(["$(arg(o,i))=$(vars[:V][v])"
                                for (i, v) in is],",")))")
    end
  end
  push!(bodstr, "end")
  exstr = "($(join(["$v=$v" for vs in values(vars) for v in vs],",") ))"
  ctxstr = "($(join(vcat(
    ["s$v::Symbol" for v in vars[:V]],
    ["$v::$k" for (k, vs) in pairs(vars) for v in vs],
    collect(Set(["i$v::Int" for v in values(bound)]))
    ),",")))"
  pstr = "($(join(params,",")),)"
  ex  = Meta.parse(exstr)
  ctx = Meta.parse(ctxstr)
  hed = Expr(:where, ex, ctx)
  bod = Meta.parse(join(bodstr, "\n"))
  # println("ex $exstr\n ctx $ctxstr\n bod $bodstr\n params $pstr")
  ps = isempty(bound) ? NamedTuple() : eval(Meta.parse(pstr))
  res = parse_relation_diagram(hed, bod)
  # Remove isolated junctions
  iso = setdiff(parts(res, :Junction), res[:junction])
  isempty(iso ∩ res[:outer_junction]) || error("$res")
  rem_parts!(res, :Junction, iso)
  return res, ps
end

function findcore(I::ACSet{CD}; loose::Bool=false)::ACSet{CD} where {CD}
  println("STARTING CORE SEARCH")
  strt = now()
  inexact = collect(filter(!=(:id_V), map(ids, ob(CD))))

  hs = [h.components[:V].func for h in (loose
          ? [homomorphism(I, I, inexact=inexact)]
          : homomorphisms(I, I, inexact=inexact))]
  minh = sort(hs, by=x->length(Set(x)))[1]
  μ = Dict(zip(I[:id_V], I[:id_V][minh]))
  res = appdict(μ, I)
  println(map(sizes, [I, res]))
  println("core dt = $(now()-strt)")
  return res
end


text:
export TensorSpace, ⊗, ProductSpace, factor, factors, nfactors

#  SV is a tuple of d spaces
abstract type AbstractProductSpace{SV,DD,RR} <: Space{DD,RR} end

spacetype(::AbstractProductSpace{SV},k) where {SV} = SV.parameters[k]

##### Tensorizer
# This gives the map from coefficients to the
# tensor entry of a tensor product of d spaces
# findfirst is overriden to get efficient inverse
# blocklengths is a tuple of block lengths, e.g., Chebyshev()^2
# would be Tensorizer((1:∞,1:∞))
# ConstantSpace() ⊗ Chebyshev()
# would be Tensorizer((1:1,1:∞))
# and Chebyshev() ⊗ ArraySpace([Chebyshev(),Chebyshev()])
# would be Tensorizer((1:∞,2:2:∞))

struct Tensorizer{DMS<:Tuple}
    blocks::DMS
end

const TrivialTensorizer{d} = Tensorizer{NTuple{d,Ones{Int,1,Tuple{OneToInf{Int}}}}}

Base.eltype(a::Tensorizer) = NTuple{length(a.blocks),Int}
Base.eltype(::Tensorizer{NTuple{d,T}}) where {d,T} = NTuple{d,Int}
dimensions(a::Tensorizer) = map(sum,a.blocks)
Base.length(a::Tensorizer) = mapreduce(sum,*,a.blocks)

# (blockrow,blockcol), (subrow,subcol), (rowshift,colshift), (numblockrows,numblockcols), (itemssofar, length)
start(a::Tensorizer{Tuple{AA,BB}}) where {AA,BB} = (1,1), (1,1), (0,0), (a.blocks[1][1],a.blocks[2][1]), (0,length(a))

function next(a::Tensorizer{Tuple{AA,BB}}, ((K,J), (k,j), (rsh,csh), (n,m), (i,tot))) where {AA,BB}
    ret = k+rsh,j+csh
    if k==n && j==m  # end of block
        if J == 1 || K == length(a.blocks[1])   # end of new block
            B = K+J # next block
            J = min(B, length(a.blocks[2]))::Int  # don't go past new block
            K = B-J+1   # K+J-1 == B
        else
            K,J = K+1,J-1
        end
        k = j = 1
        if i+1 < tot # not done yet
            n,m = a.blocks[1][K], a.blocks[2][J]
            rsh,csh = sum(a.blocks[1][1:K-1]), sum(a.blocks[2][1:J-1])
        end
    elseif k==n
        k  = 1
        j += 1
    else
        k += 1
    end
    ret, ((K,J), (k,j), (rsh,csh), (n,m), (i+1,tot))
end

done(a::Tensorizer, ((K,J), (k,j), (rsh,csh), (n,m), (i,tot))) = i ≥ tot

iterate(a::Tensorizer) = next(a, start(a))
function iterate(a::Tensorizer, st)
    done(a,st) && return nothing
    next(a, st)
end

cache(a::Tensorizer) = CachedIterator(a)

function Base.findfirst(::TrivialTensorizer{2},kj::Tuple{Int,Int})
    k,j=kj
    if k > 0 && j > 0
        n=k+j-2
        (n*(n+1))÷2+k
    else
        0
    end
end

function Base.findfirst(sp::Tensorizer{Tuple{<:AbstractFill{S},<:AbstractFill{T}}},kj::Tuple{Int,Int}) where {S,T}
    k,j=kj

    if k > 0 && j > 0
        a,b = getindex_value(sp.blocks[1]),getindex_value(sp.blocks[2])
        kb1,kr = fldmod(k-1,a)
        jb1,jr = fldmod(j-1,b)
        nb=kb1+jb1
        a*b*(nb*(nb+1)÷2+kb1)+a*jr+kr+1
    else
        0
    end
end

# which block of the tensor
# equivalent to sum of indices -1

# block(it::Tensorizer,k) = Block(sum(it[k])-length(it.blocks)+1)
block(ci::CachedIterator{T,TrivialTensorizer{2}},k::Int) where {T} =
    Block(k == 0 ? 0 : sum(ci[k])-length(ci.iterator.blocks)+1)

block(::TrivialTensorizer{2},n::Int) =
    Block(floor(Integer,sqrt(2n) + 1/2))

block(sp::Tensorizer{<:Tuple{<:AbstractFill{S},<:AbstractFill{T}}},n::Int) where {S,T} =
    Block(floor(Integer,sqrt(2floor(Integer,(n-1)/(getindex_value(sp.blocks[1])*getindex_value(sp.blocks[2])))+1) + 1/2))
_cumsum(x) = cumsum(x)
_cumsum(x::Number) = x
block(sp::Tensorizer,k::Int) = Block(findfirst(x->x≥k, _cumsum(blocklengths(sp))))
block(sp::CachedIterator,k::Int) = block(sp.iterator,k)

blocklength(it,k) = blocklengths(it)[k]
blocklength(it,k::Block) = blocklength(it,k.n[1])
blocklength(it,k::BlockRange) = blocklength(it,Int.(k))

blocklengths(::TrivialTensorizer{2}) = 1:∞

blocklengths(it::Tensorizer) = tensorblocklengths(it.blocks...)
blocklengths(it::CachedIterator) = blocklengths(it.iterator)

function getindex(it::TrivialTensorizer{2},n::Integer)
    m=block(it,n)
    p=findfirst(it,(1,m))
    j=1+n-p
    j,m-j+1
end

# could be cleaned up using blocks
function getindex(it::Tensorizer{<:Tuple{<:AbstractFill{S},<:AbstractFill{T}}},n::Integer) where {S,T}
    a,b = getindex_value(it.blocks[1]),getindex_value(it.blocks[2])
    nb1,nr = fldmod(n-1,a*b) # nb1 = "nb" - 1, i.e. using zero-base
    m1=block(it,n).n[1]-1
    pb1=fld(findfirst(it,(1,b*m1+1))-1,a*b)
    jb1=nb1-pb1
    kr1,jr1 = fldmod(nr,a)
    b*jb1+jr1+1,a*(m1-jb1)+kr1+1
end

blockstart(it,K)::Int = K==1 ? 1 : sum(blocklengths(it)[1:K-1])+1
blockstop(it,::Infinity) = ∞
_K_sum(bl::AbstractVector, K) = sum(bl[1:K])
_K_sum(bl::Integer, K) = bl
blockstop(it, K)::Int = _K_sum(blocklengths(it), K)

blockstart(it,K::Block) = blockstart(it,K.n[1])
blockstop(it,K::Block) = blockstop(it,K.n[1])

blockrange(it,K) = blockstart(it,K):blockstop(it,K)
blockrange(it,K::BlockRange) = blockstart(it,first(K)):blockstop(it,last(K))

# convert from block, subblock to tensor
subblock2tensor(rt::TrivialTensorizer{2},K,k) =
    (k,K.n[1]-k+1)

subblock2tensor(rt::CachedIterator{II,TrivialTensorizer{2}},K,k) where {II} =
    (k,K.n[1]-k+1)

subblock2tensor(rt::CachedIterator,K,k) = rt[blockstart(rt,K)+k-1]

# tensorblocklengths gives calculates the block sizes of each tensor product
#  Tensor product degrees are taken to be the sum of the degrees
#  a degree is which block you are in

tensorblocklengths(a) = a   # a single block is not modified
tensorblocklengths(a, b) = conv(a,b)
tensorblocklengths(a,b,c,d...) = tensorblocklengths(tensorblocklengths(a,b),c,d...)

# TensorSpace
# represents the tensor product of several subspaces
"""
    TensorSpace(a::Space,b::Space)

represents a tensor product of two 1D spaces `a` and `b`.
The coefficients are interlaced in lexigraphical order.

For example, consider
```julia
Fourier()*Chebyshev()  # returns TensorSpace(Fourier(),Chebyshev())
```
This represents functions on `[-π,π) x [-1,1]`, using the Fourier basis for the first argument
and Chebyshev basis for the second argument, that is, `φ_k(x)T_j(y)`, where
```
φ_0(x) = 1,
φ_1(x) = sin x,
φ_2(x) = cos x,
φ_3(x) = sin 2x,
φ_4(x) = cos 2x
…
```
By Choosing `(k,j)` appropriately, we obtain a single basis:
```
φ_0(x)T_0(y) (= 1),
φ_0(x)T_1(y) (= y),
φ_1(x)T_0(y) (= sin x),
φ_0(x)T_2(y), …
```
"""
struct TensorSpace{SV,D,R} <:AbstractProductSpace{SV,D,R}
    spaces::SV
end

tensorizer(sp::TensorSpace) = Tensorizer(map(blocklengths,sp.spaces))
blocklengths(S::TensorSpace) = tensorblocklengths(map(blocklengths,S.spaces)...)

# the evaluation is *, so the type will be the same as *
# However, this fails for some any types
tensor_eval_type(a,b) = Base.promote_op(*,a,b)
tensor_eval_type(::Type{Vector{Any}},::Type{Vector{Any}}) = Vector{Any}
tensor_eval_type(::Type{Vector{Any}},_) = Vector{Any}
tensor_eval_type(_,::Type{Vector{Any}}) = Vector{Any}

TensorSpace(sp::Tuple) =
    TensorSpace{typeof(sp),typeof(mapreduce(domain,×,sp)),
                mapreduce(rangetype,(a,b)->tensor_eval_type(a,b),sp)}(sp)

dimension(sp::TensorSpace) = mapreduce(dimension,*,sp.spaces)

for OP in (:spacescompatible,:(==))
    @eval $OP(A::TensorSpace{SV,D,R},B::TensorSpace{SV,D,R}) where {SV,D,R} =
        all(Bool[$OP(A.spaces[k],B.spaces[k]) for k=1:length(A.spaces)])
end

canonicalspace(T::TensorSpace) = TensorSpace(map(canonicalspace,T.spaces))

TensorSpace(A...) = TensorSpace(tuple(A...))
TensorSpace(A::ProductDomain) = TensorSpace(tuple(map(Space,A.domains)...))
⊗(A::TensorSpace,B::TensorSpace) = TensorSpace(A.spaces...,B.spaces...)
⊗(A::TensorSpace,B::Space) = TensorSpace(A.spaces...,B)
⊗(A::Space,B::TensorSpace) = TensorSpace(A,B.spaces...)
⊗(A::Space,B::Space) = TensorSpace(A,B)

domain(f::TensorSpace) = ×(domain.(f.spaces)...)
Space(sp::ProductDomain) = TensorSpace(sp)

setdomain(sp::TensorSpace, d::ProductDomain) = TensorSpace(setdomain.(factors(sp), factors(d)))

*(A::Space, B::Space) = A⊗B
^(A::Space, p::Integer) = p == 1 ? A : A*A^(p-1)

## TODO: generalize
components(sp::TensorSpace{Tuple{S1,S2}}) where {S1<:Space{D,R},S2} where {D,R<:AbstractArray} =
    [s ⊗ sp.spaces[2] for s in components(sp.spaces[1])]

components(sp::TensorSpace{Tuple{S1,S2}}) where {S1,S2<:Space{D,R}} where {D,R<:AbstractArray} =
    [sp.spaces[1] ⊗ s for s in components(sp.spaces[2])]

Base.size(sp::TensorSpace{Tuple{S1,S2}}) where {S1<:Space{D,R},S2} where {D,R<:AbstractArray} =
    size(sp.spaces[1])

Base.size(sp::TensorSpace{Tuple{S1,S2}}) where {S1,S2<:Space{D,R}} where {D,R<:AbstractArray} =
    size(sp.spaces[2])

# TODO: Generalize to higher dimensions
getindex(sp::TensorSpace{Tuple{S1,S2}},k::Integer) where {S1<:Space{D,R},S2} where {D,R<:AbstractArray} =
    sp.spaces[1][k] ⊗ sp.spaces[2]

getindex(sp::TensorSpace{Tuple{S1,S2}},k::Integer) where {S1,S2<:Space{D,R}} where {D,R<:AbstractArray} =
    sp.spaces[1] ⊗ sp.spaces[2][k]

length(sp::TensorSpace{Tuple{S1,S2}}) where {S1<:Space{D,R},S2} where {D,R<:AbstractArray} =
    length(sp.spaces[1])

length(sp::TensorSpace{Tuple{S1,S2}}) where {S1,S2<:Space{D,R}} where {D,R<:AbstractArray} =
    length(sp.spaces[2])

iterate(sp::TensorSpace{Tuple{S1,S2}},k...) where {S1<:Space{D,R},S2} where {D,R<:AbstractArray} =
    iterate(components(sp),k...)

iterate(sp::TensorSpace{Tuple{S1,S2}},k...) where {S1,S2<:Space{D,R}} where {D,R<:AbstractArray} =
    iterate(components(sp),k...)

# every column is in the same space for a TensorSpace
# TODO: remove
columnspace(S::TensorSpace,_) = S.spaces[1]

struct ProductSpace{S<:Space,V<:Space,D,R} <: AbstractProductSpace{Tuple{S,V},D,R}
    spacesx::Vector{S}
    spacey::V
end

ProductSpace(spacesx::Vector,spacey) =
    ProductSpace{eltype(spacesx),typeof(spacey),typeof(mapreduce(domain,×,sp)),
                mapreduce(s->eltype(domain(s)),promote_type,sp)}(spacesx,spacey)

# TODO: This is a weird definition
⊗(A::Vector{S},B::Space) where {S<:Space} = ProductSpace(A,B)
domain(f::ProductSpace) = domain(f.spacesx[1])×domain(f.spacesy)

nfactors(d::AbstractProductSpace) = length(d.spaces)
factors(d::AbstractProductSpace) = d.spaces
factor(d::AbstractProductSpace,k) = factors(d)[k]

isambiguous(A::TensorSpace) = isambiguous(A.spaces[1]) || isambiguous(A.spaces[2])

Base.transpose(d::TensorSpace) = TensorSpace(d.spaces[2],d.spaces[1])

## Transforms

for (plan, plan!, Typ) in ((:plan_transform, :plan_transform!, :TransformPlan),
                           (:plan_itransform, :plan_itransform!, :ITransformPlan))
    @eval begin
        $plan!(S::TensorSpace, M::AbstractMatrix) = $Typ(S,(($plan(S.spaces[1],size(M,1)),size(M,1)),
                                                             ($plan(S.spaces[2],size(M,2)),size(M,2))),
                                                             Val{true})

        function *(T::$Typ{<:Any,<:TensorSpace,true}, M::AbstractMatrix)
            n=size(M,1)

            for k=1:size(M,2)
                M[:,k]=T.plan[1][1]*M[:,k]
            end
            for k=1:n
                M[k,:]=T.plan[2][1]*M[k,:]
            end
            M
        end

        function *(T::$Typ{TT,SS,false},v::AbstractVector) where {SS<:TensorSpace,TT}
            P = $Typ(T.space,T.plan,Val{true})
            P*AbstractVector{rangetype(SS)}(v)
        end
    end
end

function plan_transform(sp::TensorSpace, ::Type{T}, n::Integer) where {T}
    NM=n
    if isfinite(dimension(sp.spaces[1])) && isfinite(dimension(sp.spaces[2]))
        N,M=dimension(sp.spaces[1]),dimension(sp.spaces[2])
    elseif isfinite(dimension(sp.spaces[1]))
        N=dimension(sp.spaces[1])
        M=NM÷N
    elseif isfinite(dimension(sp.spaces[2]))
        M=dimension(sp.spaces[2])
        N=NM÷M
    else
        N=M=round(Int,sqrt(n))
    end

    TransformPlan(sp,((plan_transform(sp.spaces[1],T,N),N),
                    (plan_transform(sp.spaces[2],T,M),M)),
                Val{false})
end

function plan_transform!(sp::TensorSpace, ::Type{T}, n::Integer) where {T}
    P = plan_transform(sp, T, n)
    TransformPlan(sp, P.plan, Val{true})
end

plan_transform(sp::TensorSpace, v::AbstractVector) = plan_transform(sp,eltype(v),length(v))
plan_transform!(sp::TensorSpace, v::AbstractVector) = plan_transform!(sp,eltype(v),length(v))

function plan_itransform(sp::TensorSpace, v::AbstractVector{T}) where {T}
    N,M = size(totensor(sp, v)) # wasteful
    ITransformPlan(sp,((plan_itransform(sp.spaces[1],T,N),N),
                    (plan_itransform(sp.spaces[2],T,M),M)),
                Val{false})
end

function *(T::TransformPlan{TT,<:TensorSpace,true},v::AbstractVector) where TT # need where TT
    N,M = T.plan[1][2],T.plan[2][2]
    V=reshape(v,N,M)
    fromtensor(T.space,T*V)
end

*(T::ITransformPlan{TT,<:TensorSpace,true},v::AbstractVector) where TT  =
    vec(T*totensor(T.space,v))

## points

points(d::Union{EuclideanDomain{2},BivariateSpace},n,m) = points(d,n,m,1),points(d,n,m,2)

function points(d::BivariateSpace,n,m,k)
    ptsx=points(columnspace(d,1),n)
    ptst=points(factor(d,2),m)

    promote_type(eltype(ptsx),eltype(ptst))[fromcanonical(d,x,t)[k] for x in ptsx, t in ptst]
end

##  Fun routines

fromtensor(S::Space,M::AbstractMatrix) = fromtensor(tensorizer(S),M)
totensor(S::Space,M::AbstractVector) = totensor(tensorizer(S),M)

# we only copy upper triangular of coefficients
function fromtensor(it::Tensorizer,M::AbstractMatrix)
    n,m=size(M)
    ret=zeros(eltype(M),blockstop(it,max(n,m)))
    k = 1
    for (K,J) in it
        if k > length(ret)
            break
        end
        if K ≤ n && J ≤ m
            ret[k] = M[K,J]
        end
        k += 1
    end
    ret
end

function totensor(it::Tensorizer,M::AbstractVector)
    n=length(M)
    B=block(it,n)
    ds = dimensions(it)

    ret=zeros(eltype(M),sum(it.blocks[1][1:min(B.n[1],length(it.blocks[1]))]),
                        sum(it.blocks[2][1:min(B.n[1],length(it.blocks[2]))]))
    k=1
    for (K,J) in it
        if k > n
            break
        end
        ret[K,J] = M[k]
        k += 1
    end
    ret
end

for OP in (:block,:blockstart,:blockstop)
    @eval begin
        $OP(s::TensorSpace, ::Infinity) = ∞
        $OP(s::TensorSpace, M::Block) = $OP(tensorizer(s),M)
        $OP(s::TensorSpace, M) = $OP(tensorizer(s),M)
    end
end

function points(sp::TensorSpace,n)
    pts=Array{float(eltype(domain(sp)))}(undef,0)
    a,b = sp.spaces
    if isfinite(dimension(a)) && isfinite(dimension(b))
        N,M=dimension(a),dimension(b)
    elseif isfinite(dimension(a))
        N=dimension(a)
        M=n÷N
    elseif isfinite(dimension(b))
        M=dimension(b)
        N=n÷M
    else
        N=M=round(Int,sqrt(n))
    end

    for y in points(b,M),
        x in points(a,N)
        push!(pts,Vec(x...,y...))
    end
    pts
end

itransform(sp::TensorSpace,cfs) = vec(itransform!(sp,coefficientmatrix(Fun(sp,cfs))))

evaluate(f::AbstractVector,S::AbstractProductSpace,x) = ProductFun(totensor(S,f),S)(x...)
evaluate(f::AbstractVector,S::AbstractProductSpace,x,y) = ProductFun(totensor(S,f),S)(x,y)

coefficientmatrix(f::Fun{<:AbstractProductSpace}) = totensor(space(f),f.coefficients)

#TODO: Implement
# function ∂(d::TensorSpace{<:IntervalOrSegment{Float64}})
#     @assert length(d.spaces) ==2
#     PiecewiseSpace([d[1].a+im*d[2],d[1].b+im*d[2],d[1]+im*d[2].a,d[1]+im*d[2].b])
# end

union_rule(a::TensorSpace,b::TensorSpace) = TensorSpace(map(union,a.spaces,b.spaces))

## Convert from 1D to 2D

# function isconvertible{T,TT}(sp::Space{Segment{Vec{2,TT}},<:Real},ts::TensorSpace)
#     d1 = domain(sp)
#     d2 = domain(ts)
#     if d2
#     length(ts.spaces) == 2 &&
#     ((domain(ts)[1] == Point(0.0) && isconvertible(sp,ts.spaces[2])) ||
#      (domain(ts)[2] == Point(0.0) && isconvertible(sp,ts.spaces[1])))
#  end

isconvertible(sp::UnivariateSpace,ts::TensorSpace{SV,D,R}) where {SV,D<:EuclideanDomain{2},R} = length(ts.spaces) == 2 &&
    ((domain(ts)[1] == Point(0.0) && isconvertible(sp,ts.spaces[2])) ||
     (domain(ts)[2] == Point(0.0) && isconvertible(sp,ts.spaces[1])))

# coefficients(f::AbstractVector,sp::ConstantSpace,ts::TensorSpace{SV,D,R}) where {SV,D<:EuclideanDomain{2},R} =
#     f[1]*ones(ts).coefficients

#
# function coefficients(f::AbstractVector,sp::Space{IntervalOrSegment{Vec{2,TT}}},ts::TensorSpace{Tuple{S,V},D,R}) where {S,V<:ConstantSpace,D<:EuclideanDomain{2},R,TT} where {T<:Number}
#     a = domain(sp)
#     b = domain(ts)
#     # make sure we are the same domain. This will be replaced by isisomorphic
#     @assert first(a) ≈ Vec(first(factor(b,1)),factor(b,2).x) &&
#         last(a) ≈ Vec(last(factor(b,1)),factor(b,2).x)
#
#     coefficients(f,sp,setdomain(factor(ts,1),a))
# end

function coefficients(f::AbstractVector,sp::UnivariateSpace,ts::TensorSpace{SV,D,R}) where {SV,D<:EuclideanDomain{2},R}
    @assert length(ts.spaces) == 2

    if factor(domain(ts),1) == Point(0.0)
        coefficients(f,sp,ts.spaces[2])
    elseif factor(domain(ts),2) == Point(0.0)
        coefficients(f,sp,ts.spaces[1])
    else
        error("Cannot convert coefficients from $sp to $ts")
    end
end

function isconvertible(sp::Space{Segment{Vec{2,TT}}},ts::TensorSpace{SV,D,R}) where {TT,SV,D<:EuclideanDomain{2},R}
    d1 = domain(sp)
    d2 = domain(ts)
    if length(ts.spaces) ≠ 2
        return false
    end
    if d1.a[2] ≈ d1.b[2]
        isa(factor(d2,2),Point) && factor(d2,2).x ≈ d1.a[2] &&
            isconvertible(setdomain(sp,Segment(d1.a[1],d1.b[1])),ts[1])
    elseif d1.a[1] ≈ d1.b[1]
        isa(factor(d2,1),Point) && factor(d2,1).x ≈ d1.a[1] &&
            isconvertible(setdomain(sp,Segment(d1.a[2],d1.b[2])),ts[2])
    else
        return false
    end
end

function coefficients(f::AbstractVector,sp::Space{Segment{Vec{2,TT}}},
                            ts::TensorSpace{SV,D,R}) where {TT,SV,D<:EuclideanDomain{2},R}
    @assert length(ts.spaces) == 2
    d1 = domain(sp)
    d2 = domain(ts)
    if d1.a[2] ≈ d1.b[2]
        coefficients(f,setdomain(sp,Segment(d1.a[1],d1.b[1])),factor(ts,1))
    elseif d1.a[1] ≈ d1.b[1]
        coefficients(f,setdomain(sp,Segment(d1.a[2],d1.b[2])),factor(ts,2))
    else
        error("Cannot convert coefficients from $sp to $ts")
    end
end

Fun(::typeof(identity), S::TensorSpace) = Fun(xyz->collect(xyz),S)

ext: jl


text:
struct GraphSVMModel
    svm::LIBSVM.SVM
    kernel::KernelFunctions.Kernel
    graphs::AbstractVector
end

function svmtrain(graphs::AbstractVector{<:AbstractGraph}, labels, kernel::KernelFunctions.Kernel; kwargs...)

    X = kernelmatrix(kernel, graphs)
    svm = svmtrain(X, labels, kernel=Kernel.Precomputed; kwargs...)

    return GraphSVMModel(svm, kernel, graphs)
end

function svmpredict(model::GraphSVMModel, unpredicted_graphs::AbstractVector{<:AbstractGraph})

    graphs = model.graphs
    kernel = model.kernel

    # TODO might only be necessary to do the calculations for support vectors
    X = kernelmatrix(kernel, graphs, unpredicted_graphs)
    return svmpredict(model.svm, X)[1] # for simplicity return only the labels for now
end


text:
# functions related to chi-square distribution

import .RFunctions:
    chisqpdf,
    chisqlogpdf,
    chisqcdf,
    chisqccdf,
    chisqlogcdf,
    chisqlogccdf,
    chisqinvcdf,
    chisqinvccdf,
    chisqinvlogcdf,
    chisqinvlogccdf


text:
include("ToricDivisors.jl")

######################
# 1: Attributes of ToricVarieties
######################

function affine_open_covering( v::toric_variety )
    gap_cover = GAP.Globals.AffineOpenCovering( v.GapToricVariety )
    return [ toric_variety( v ) for v in gap_cover ]
end
export affine_open_covering

struct cox_ring
           GapCoxRing::GapObj
end
export cox_ring

function cox_ring( v::toric_variety )
    gap_ring = GAP.Globals.CoxRing( v.GapToricVariety )
    return cox_ring( gap_ring )
end
export cox_ring

function list_of_variables_of_cox_ring(v::toric_variety)
    vars = GAP.Globals.ListOfVariablesOfCoxRing(v.GapToricVariety)
    return Vector{String}(vars)
end
export list_of_variables_of_cox_ring

struct class_group
           GapClassGroup::GapObj
end
export class_group

function class_group( v::toric_variety )
    gap_class_group = GAP.Globals.ClassGroup( v.GapToricVariety )
    return class_group( gap_class_group )
end
export class_group

struct torus_invariant_divisor_group
           GapTorusInvariantDivisorGroup::GapObj
end
export torus_invariant_divisor_group

function torus_invariant_divisor_group( v::toric_variety )
    gap_TorusInvariantDivisorGroup = GAP.Globals.TorusInvariantDivisorGroup( v.GapToricVariety )
    return torus_invariant_divisor_group( gap_TorusInvariantDivisorGroup )
end
export torus_invariant_divisor_group

struct map_from_character_to_principal_divisor
           GapMapFromCharacterToPrincipalDivisor::GapObj
end
export map_from_character_to_principal_divisor

function map_from_character_to_principal_divisor( v::toric_variety )
    gap_MapFromCharacterToPrincipalDivisor = GAP.Globals.MapFromCharacterToPrincipalDivisor( v.GapToricVariety )
    return map_from_character_to_principal_divisor( gap_MapFromCharacterToPrincipalDivisor )
end
export map_from_character_to_principal_divisor

struct map_from_weil_divisors_to_class_group
           GapMapFromWeilDivisorsToClassGroup::GapObj
end
export map_from_weil_divisors_to_class_group

function map_from_weil_divisors_to_class_group( v::toric_variety )
    gap_MapFromWeilDivisorsToClassGroup = GAP.Globals.MapFromWeilDivisorsToClassGroup( v.GapToricVariety )
    return map_from_weil_divisors_to_class_group( gap_MapFromWeilDivisorsToClassGroup )
end
export map_from_weil_divisors_to_class_group

function dimension( v::toric_variety )
    return GAP.Globals.Dimension(v.GapToricVariety)::Int
end
export dimension

function dimension_of_torusfactor( v::toric_variety )
    return GAP.Globals.DimensionOfTorusfactor( v.GapToricVariety )::Int
end
export dimension_of_torusfactor

struct coordinate_ring_of_torus
           GapCoordinateRingOfTorus::GapObj
end
export coordinate_ring_of_torus

function coordinate_ring_of_torus( v::toric_variety )
    gap_CoordinateRingOfTorus = GAP.Globals.CoordinateRingOfTorus( v.GapToricVariety )
    return coordinate_ring_of_torus( gap_CoordinateRingOfTorus )
end
export coordinate_ring_of_torus

function list_of_variables_of_coordinate_ring_of_torus(v::toric_variety)
    vars = GAP.Globals.ListOfVariablesOfCoordinateRingOfTorus(v.GapToricVariety)
    return Vector{String}(vars)
end
export list_of_variables_of_coordinate_ring_of_torus

function is_product_of(v::toric_variety)
    factors = GAP.Globals.IsProductOf(v.GapToricVariety)
    return [ toric_variety( f ) for f in factors ]
end
export is_product_of

struct character_lattice
           GapCharacterLattice::GapObj
end
export character_lattice

function character_lattice( v::toric_variety )
    gap_CharacterLattice = GAP.Globals.CharacterLattice( v.GapToricVariety )
    return character_lattice( gap_CharacterLattice )
end
export character_lattice

function torus_invariant_prime_divisors( v::toric_variety )
    divisors = GAP.Globals.TorusInvariantPrimeDivisors( v.GapToricVariety )
    return [ toric_divisor( d ) for d in divisors ]
end
export torus_invariant_prime_divisors

struct irrelevant_ideal
           GapIrrelevantIdeal::GapObj
end
export irrelevant_ideal

function irrelevant_ideal( v::toric_variety )
    gap_IrrelevantIdeal = GAP.Globals.IrrelevantIdeal( v.GapToricVariety )
    return irrelevant_ideal( gap_IrrelevantIdeal )
end
export irrelevant_ideal

struct stanley_reisner_ideal
           GapSRIdeal::GapObj
end
export stanley_reisner_ideal

function stanley_reisner_ideal( v::toric_variety )
    gap_SRIdeal = GAP.Globals.SRIdeal( v.GapToricVariety )
    return stanley_reisner_ideal( gap_SRIdeal )
end
export stanley_reisner_ideal

struct morphism_from_cox_variety
           GapMorphismFromCoxVariety::GapObj
end
export morphism_from_cox_variety

function morphism_from_cox_variety( v::toric_variety )
    gap_MorphismFromCoxVariety = GAP.Globals.MorphismFromCoxVariety( v.GapToricVariety )
    return morphism_from_cox_variety( gap_MorphismFromCoxVariety )
end
export morphism_from_cox_variety

function cox_variety( v::toric_variety )
    gap_CoxVariety = GAP.Globals.CoxVariety( v.GapToricVariety )
    return toric_variety( gap_CoxVariety )
end
export cox_variety

struct fan
           gap_fan::GapObj
           rays::Vector{Vector{Int}}
           cones::Vector{Vector{Int}}
end
export fan

function fan_of_variety( v::toric_variety )
    # collect data
    gap_fan = GAP.Globals.FanOfVariety( v.GapToricVariety )
    rays = Vector{Vector{Int}}( GAP.Globals.RayGenerators( gap_fan ) )
    cones = Vector{Vector{Int}}( GAP.Globals.RaysInMaximalCones( gap_fan ) )
    cones = [ findall( x -> x == 1, c ) for c in cones ]

    # return the fan
    return fan( gap_fan, rays, cones )
end
export fan_of_variety

function fan( v::toric_variety )
    return fan_of_variety( v )
end
export fan

struct cartier_torus_invariant_divisor_group
           GapCartierTorusInvariantDivisorGroup::GapObj
end
export cartier_torus_invariant_divisor_group

function cartier_torus_invariant_divisor_group( v::toric_variety )
    gap_CartierTorusInvariantDivisorGroup = GAP.Globals.CartierTorusInvariantDivisorGroup( v.GapToricVariety )
    return cartier_torus_invariant_divisor_group( gap_CartierTorusInvariantDivisorGroup )
end
export cartier_torus_invariant_divisor_group

struct picard_group
           GapPicardGroup::GapObj
end
export picard_group

function picard_group( v::toric_variety )
    gap_PicardGroup = GAP.Globals.PicardGroup( v.GapToricVariety )
    return picard_group( gap_PicardGroup )
end
export picard_group

function name_of_variety( v::toric_variety )
    if ! ( Bool( GAP.Globals.HasNameOfVariety( v.GapToricVariety ) ) )
            return "No name set for this variety"
    end

    return String( GAP.Globals.NameOfVariety( v.GapToricVariety ) )
end
export name_of_variety

function set_name_of_variety( v::toric_variety, s::String )
    GAP.Globals.SetNameOfVariety( v.GapToricVariety, GapObj( s ) )
    return true
end
export set_name_of_variety

struct zariski_cotangent_sheaf
           GapZariskiCotangentSheaf::GapObj
end
export zariski_cotangent_sheaf

function zariski_cotangent_sheaf( v::toric_variety )
    gap_ZariskiCotangentSheaf = GAP.Globals.ZariskiCotangentSheaf( v.GapToricVariety )
    return zariski_cotangent_sheaf( gap_ZariskiCotangentSheaf )
end
export zariski_cotangent_sheaf

struct cotangent_sheaf
           GapCotangentSheaf::GapObj
end
export cotangent_sheaf

function cotangent_sheaf( v::toric_variety )
    gap_CotangentSheaf = GAP.Globals.CotangentSheaf( v.GapToricVariety )
    return cotangent_sheaf( gap_CotangentSheaf )
end
export cotangent_sheaf

function euler_characteristic( v::toric_variety )
    return GAP.Globals.EulerCharacteristic( v.GapToricVariety )::Int
end
export euler_characteristic

#struct underlying_sheaf
#          GapUnderlyingSheaf::GapObj
#end
#export underlying_sheaf

#function underlying_sheaf( v::toric_variety )
#   gap_Underlying = GAP.Globals.UnderlyingSheaf( v.GapToricVariety )
#    return underlying_sheaf( gap_UnderlyingSheaf )
#end
#export underlying_sheaf

######################
# 2: Methods of ToricVarieties
######################

function coordinate_ring_of_torus( v::toric_variety, names::Vector{String} )
    gap_names = [ GapObj( names[ i ] ) for i in 1 : size(names)[1] ]
    gap_names = GapObj( gap_names )
    gap_CoordinateRingOfTorus = GAP.Globals.CoordinateRingOfTorus( v.GapToricVariety, gap_names )
    return coordinate_ring_of_torus( gap_CoordinateRingOfTorus )
end
export coordinate_ring_of_torus

function cox_ring( v::toric_variety, names::String )
    gap_names = GapObj( names )
    gap_CoxRing = GAP.Globals.CoxRing( v.GapToricVariety, gap_names )
    return cox_ring( gap_CoxRing )
end
export cox_ring

function Base.:*( v::toric_variety, w::toric_variety )
    gap_ToricVariety = v.GapToricVariety * w.GapToricVariety
    return toric_variety( gap_ToricVariety )
end
export *

#=

#! @Description
#!  Computes the rational function corresponding to the character grid element elem or to the list of
#!  integers elem. This computation needs to know the coordinate ring of the torus of the variety vari. By
#!  default this ring is introduced with variables x1 to xn where n is the dimension of the variety. If
#!  different variables should be used, then CoordinateRingOfTorus has to be set accordingly before calling this method.
#! @Returns a homalg element
#! @Arguments elem, vari
DeclareOperation( "CharacterToRationalFunction",
                  [ IsList, IsToricVariety ] );

#! @Description
#!  Returns a list of the currently defined Divisors of the toric variety.
#! @Returns a list
#! @Arguments vari
DeclareOperation( "WeilDivisorsOfVariety",
                  [ IsToricVariety ] );

#! @Description
#!
#! @Arguments vari
DeclareOperation( "Factors",
                  [ IsToricVariety ] );

#! @Description
#!
#! @Arguments vari, p
DeclareOperation( "BlowUpOnIthMinimalTorusOrbit",
                  [ IsToricVariety, IsInt ] );

#! @Description
#!
DeclareGlobalFunction( "ZariskiCotangentSheafViaEulerSequence" );

#! @Description
#!
DeclareGlobalFunction( "ZariskiCotangentSheafViaPoincareResidueMap" );

#! @Description
#!
#! @Arguments vari, p
DeclareOperation( "ithBettiNumber",
                  [ IsToricVariety, IsInt ] );

#! @Description
#!
#! @Arguments vari, p
DeclareOperation( "NrOfqRationalPoints",
                  [ IsToricVariety, IsInt ] );

=#

max_stars_repo_name: fingolfin/JToric

max_issues_repo_name: fingolfin/JToric

max_forks_repo_name: fingolfin/JToric


text:
#=
    General helping functions
=#
logprior(m::AbstractBQModel) = Base.Fix1(logpdf, p_0(m))
prior(m::AbstractBQModel) = exp ∘ logprior(m)

integrand(m::AbstractBQModel) = exp ∘ logintegrand(m)
logjoint(m::AbstractBQModel) = x -> logprior(m)(x) + logintegrand(m)(x)

"""
    BayesModel(prior, logintegrand) <: AbstractBQModel

Model inheriting from `AbstractMCMC.AbstractModel`.
`prior` should be a multivariate distribution from `Distributions.jl`
at the moment `prior` has to be a `MvNormal` but this will improved in a later version
`logintegrand` should be the log of the function to integrate.
"""
struct BayesModel{Tp,Ti} <: AbstractBQModel{Tp,Ti}
    prior::Tp
    logintegrand::Ti
end

p_0(m::BayesModel{<:AbstractMvNormal}) = m.prior
p_0(m::BayesModel) = MvNormal(ones(dim(m.prior)))
logintegrand(m::BayesModel{<:AbstractMvNormal}) = m.logintegrand
function logintegrand(m::AbstractBQModel)
    return function reweightedlogintegrand(x)
        return m.logintegrand(x) + logpdf(m.prior, x) - logprior(m)(x)
    end
end


text:
abstract type DataSet{DS} end

getproperty(ds::DS, k::Symbol) where {DS<:DataSet{<:DataSet}} =
    hasfield(DS, k) ? getfield(ds, k) : getproperty(getfield(ds, :_super), k)
setproperty!(ds::DS, k::Symbol, v) where {DS<:DataSet{<:DataSet}} =
    hasfield(DS, k) ? setfield!(ds, k, v) : setproperty!(getfield(ds, :_super), k, v)
propertynames(ds::DS) where {DS′<:DataSet, DS<:DataSet{DS′}} =
    union(fieldnames(DS), fieldnames(DS′))

function new_dataset(::Type{DS}; kwargs...) where {DS′<:DataSet, DS<:DataSet{DS′}}
    kw  = filter(((k,_),)->  k in fieldnames(DS),  kwargs)
    kw′ = filter(((k,_),)->!(k in fieldnames(DS)), kwargs)
    DS(_super=DS′(;kw′...); kw...)
end

copy(ds::DS) where {DS<:DataSet} =
    DS(((k==:_super ? copy(v) : v) for (k,v) in pairs(fields(ds)))...)

# needed until fix to https://github.com/FluxML/Zygote.jl/issues/685
Zygote.grad_mut(ds::DataSet) = Ref{Any}((;(propertynames(ds) .=> nothing)...))

# Stores variables needed to construct the posterior
@kwdef mutable struct BaseDataSet <: DataSet{Nothing}
    d                # data
    Cϕ               # ϕ covariance
    Cf               # unlensed field covariance
    Cf̃ = nothing     # lensed field covariance (not always needed)
    Cn               # noise covariance
    Cn̂ = Cn          # approximate noise covariance, diagonal in same basis as Cf
    M  = 1           # user mask
    M̂  = M           # approximate user mask, diagonal in same basis as Cf
    B  = 1           # beam and instrumental transfer functions
    B̂  = B           # approximate beam and instrumental transfer functions, diagonal in same basis as Cf
    D  = 1           # mixing matrix for mixed parametrization
    G  = 1           # reparametrization for ϕ
    P  = 1           # pixelization operator (if estimating field on higher res than data)
    L  = LenseFlow   # lensing operator, possibly cached for memory reuse
end

function subblock(ds::DS, block) where {DS<:DataSet}
    DS(map(collect(pairs(fields(ds)))) do (k,v)
        @match (k,v) begin
            ((:Cϕ || :G || :L), v)              => v
            (_, L::Union{Nothing,FuncOp,Real})  => L
            (_, L)                              => getindex(L,block)
        end
    end...)
end

function (ds::DataSet)(θ::NamedTuple)
    DS = typeof(ds)
    DS(map(fieldvalues(ds)) do v
        (v isa Union{ParamDependentOp,DataSet}) ? v(θ) : v
    end...)
end
(ds::DataSet)(;θ...) = ds((;θ...))

adapt_structure(to, ds::DS) where {DS <: DataSet} = DS(adapt(to, fieldvalues(ds))...)

@doc doc"""
    resimulate(ds::DataSet; [f, ϕ, n])

Make a new DataSet with the data replaced by a simulation, potentially given a
fixed f, ϕ, or n, if any are provided.

Returns a named tuple of `(ds, f, ϕ, n, f̃)`
"""
resimulate(ds::DataSet; kwargs...) = resimulate!(copy(ds); kwargs...)

@doc doc"""
    resimulate!(ds::DataSet; [f, ϕ, n])

Replace the data in this DataSet in-place with a simulation, potentially given a
fixed f, ϕ, or n, if any are provided.

Returns a named tuple of `(ds, f, ϕ, n, f̃)`
"""
function resimulate!(
    ds::DataSet;
    f=nothing, ϕ=nothing, n=nothing, f̃=nothing,
    rng=global_rng_for(ds.d), seed=nothing
)

    @unpack M,P,B,L,Cϕ,Cf,Cn,d = ds()
    D = batchsize(d)

    if (f̃ == nothing)
        if (ϕ == nothing)
            ϕ = simulate(batch(Cϕ,D), rng=rng, seed=seed)
        end
        if (f == nothing)
            f = simulate(batch(Cf,D), rng=rng, seed=(seed==nothing ? nothing : seed+1))
        end
        f̃ = L(ϕ)*f
    else
        f = ϕ = nothing
    end
    if (n == nothing)
        n = simulate(batch(Cn,D), rng=rng, seed=(seed==nothing ? nothing : seed+2))
    end

    ds.d = d = M*P*B*f̃ + n

    @namedtuple(ds,f,ϕ,n,f̃,d)

end

@doc doc"""
    load_sim

Create a `BaseDataSet` object with some simulated data, returing the DataSet
and simulated truths. E.g.

```julia
@unpack f,ϕ,ds = load_sim(;
    θpix  = 2,
    Nside = 128,
    pol   = :I,
    T     = Float32
);
```
For rectangular maps, Nside expects (Ny, Nx), i.e. (Nrow, Ncol).
If Nside isa Int, the code interprets it as a square map.
"""
function load_sim(;

    # basic configuration
    θpix,
    θpix_data = θpix,
    Nside,
    pol,
    T = Float32,
    storage = Array,
    Nbatch = nothing,

    # noise parameters, or set Cℓn or even Cn directly
    μKarcminT = 3,
    ℓknee = 100,
    αknee = 3,
    Cℓn = nothing,
    Cn = nothing,

    # beam parameters, or set B directly
    beamFWHM = 0,
    B = nothing, B̂ = nothing,

    # mask parameters, or set M directly
    pixel_mask_kwargs = nothing,
    bandpass_mask = LowPass(3000),
    M = nothing, M̂ = nothing,

    # theory
    Cℓ = nothing,
    fiducial_θ = NamedTuple(),
    rfid = nothing,

    rng = global_rng_for(storage),
    seed = nothing,
    D = nothing,
    G = nothing,
    Nϕ_fac = 2,
    L = LenseFlow,
    ∂mode = fourier∂

)

    # the biggest ℓ on the 2D fourier grid
    ℓmax = round(Int,ceil(√2*fieldinfo(Flat(θpix=θpix,Nside=Nside)).nyq)+1)

    # CMB Cℓs
    if (rfid != nothing)
        @warn "`rfid` will be removed in a future version. Use `fiducial_θ=(r=...,)` instead."
        fiducial_θ = merge(fiducial_θ,(r=rfid,))
    end
    Aϕ₀ = get(fiducial_θ, :Aϕ, 1)
    fiducial_θ = Base.structdiff(fiducial_θ, NamedTuple{(:Aϕ,)}) # remove Aϕ key if present
    if (Cℓ == nothing)
        Cℓ = camb(;fiducial_θ..., ℓmax=ℓmax)
    else
        if !isempty(fiducial_θ)
            error("Can't pass both `Cℓ` and `fiducial_θ` parameters which affect `Cℓ`, choose one or the other.")
        elseif maximum(Cℓ.total.TT.ℓ) < ℓmax
            error("ℓmax of `Cℓ` argument should be higher than $ℓmax for this configuration.")
        end
    end
    r₀ = Cℓ.params.r

    # noise Cℓs (these are non-debeamed, hence beamFWHM=0 below; the beam comes in via the B operator)
    if (Cℓn == nothing)
        Cℓn = noiseCℓs(μKarcminT=μKarcminT, beamFWHM=0, ℓknee=ℓknee, αknee=αknee, ℓmax=ℓmax)
    end

    # some things which depend on whether we chose :I, :P, or :IP
    pol = Symbol(pol)
    S,ks,F,F̂,nF = @match pol begin
        :I  => (S0,  (:TT,),            FlatMap,    FlatFourier,    1)
        :P  => (S2,  (:EE,:BB),         FlatQUMap,  FlatEBFourier,  2)
        :IP => (S02, (:TT,:EE,:BB,:TE), FlatIQUMap, FlatIEBFourier, 3)
        _   => throw(ArgumentError("`pol` should be one of :I, :P, or :IP"))
    end

    # pixelization
    Pix = Flat(Nside=Nside, θpix=θpix, ∂mode=∂mode)
    if (θpix_data == θpix)
        Pix_data = Pix
        P = 1
    else
        Pix_data = Flat(Nside=Nside.÷(θpix_data÷θpix), θpix=θpix_data, ∂mode=∂mode)
        P = FuncOp(
            op  = f -> ud_grade(f, θpix_data, deconv_pixwin=false, anti_aliasing=false),
            opᴴ = f -> ud_grade(f, θpix,      deconv_pixwin=false, anti_aliasing=false)
        )
    end

    # covariances
    Cϕ₀ = adapt(storage, Cℓ_to_Cov(Pix,      T, S0, (Cℓ.total.ϕϕ)))
    Cfs = adapt(storage, Cℓ_to_Cov(Pix,      T, S,  (Cℓ.unlensed_scalar[k] for k in ks)...))
    Cft = adapt(storage, Cℓ_to_Cov(Pix,      T, S,  (Cℓ.tensor[k]          for k in ks)...))
    Cf̃  = adapt(storage, Cℓ_to_Cov(Pix,      T, S,  (Cℓ.total[k]           for k in ks)...))
    Cn̂  = adapt(storage, Cℓ_to_Cov(Pix_data, T, S,  (Cℓn[k]                for k in ks)...))
    if (Cn == nothing); Cn = Cn̂; end
    Cf = ParamDependentOp((;r=r₀,   _...)->(Cfs + T(r/r₀)*Cft))
    Cϕ = ParamDependentOp((;Aϕ=Aϕ₀, _...)->(T(Aϕ) * Cϕ₀))

    # data mask
    if (M == nothing)
        Mfourier = adapt(storage, Cℓ_to_Cov(Pix_data, T, S, ((k==:TE ? 0 : 1) * bandpass_mask.diag.Wℓ for k in ks)...; units=1))
        if (pixel_mask_kwargs != nothing)
            Mpix = adapt(storage, Diagonal(F{Pix_data}(repeated(T.(make_mask(Nside.÷(θpix_data÷θpix),θpix_data; pixel_mask_kwargs...).Ix),nF)...)))
        else
            Mpix = 1
        end
        M = Mfourier * Mpix
        if (M̂ == nothing)
            M̂ = Mfourier
        end
    else
        if (M̂ == nothing)
            M̂ = M
        end
    end
    if (M̂ isa DiagOp{<:BandPass})
        M̂ = Diagonal(M̂ * one(diag(Cf)))
    end

    # beam
    if (B == nothing)
        B = adapt(storage, Cℓ_to_Cov(Pix, T, S, ((k==:TE ? 0 : 1) * sqrt(beamCℓs(beamFWHM=beamFWHM)) for k=ks)..., units=1))
    end
    if (B̂ == nothing)
        B̂ = B
    end

    # creating lensing operator cache
    Lϕ = alloc_cache(L(diag(Cϕ)),diag(Cf))

    # put everything in DataSet
    ds = BaseDataSet(;@namedtuple(d=nothing, Cn, Cn̂, Cf, Cf̃, Cϕ, M, M̂, B, B̂, D, P, L=Lϕ)...)

    # simulate data
    @unpack ds,f,f̃,ϕ,n = resimulate(ds, rng=rng, seed=seed)

    # with the DataSet created, we now more conveniently create the mixing matrices D and G
    if (G == nothing)
        Nϕ = quadratic_estimate(ds,(pol in (:P,:IP) ? :EB : :TT)).Nϕ / Nϕ_fac
        G₀ = @. nan2zero(sqrt(1 + 2/($Cϕ()/Nϕ)))
        ds.G = ParamDependentOp((;Aϕ=Aϕ₀, _...)->(pinv(G₀) * sqrt(I + 2 * Nϕ * pinv(Cϕ(Aϕ=Aϕ)))))
    end
    if (D == nothing)
        σ²len = T(deg2rad(5/60)^2)
        ds.D = ParamDependentOp(
            function (;r=r₀, _...)
                Cfr = Cf(;r=r)
                sqrt(Diagonal(diag(Cfr) .+ σ²len .+ 2*diag(Cn̂)) * pinv(Cfr))
            end,
        )
    end

    if Nbatch != nothing
        ds.d = identity.(batch(ds.d, Nbatch))
        ds.L = alloc_cache(L(identity.(batch(zero(diag(Cϕ)), Nbatch))), ds.d)
    end

    return adapt(storage, @namedtuple(f, f̃, ϕ, n, ds, ds₀=ds(), T, P=Pix, Cℓ, L))

end

@deprecate load_sim_dataset(args...; kwargs...) load_sim(args...; kwargs...)


text:
push!(LOAD_PATH, "./jl")
using RBT
include("../tpch.jl")

setdb(tpch)

function RBT.query(state, expr; params...)
    params_signature = ([(name, typeof(param)) for (name, param) in params]...)
    q = RBT.prepare(state, expr; params_signature...)
    println(repr(q))
    println(repr(q.pipe))
    r = RBT.execute(q; params...)
    #display(r)
    #println()
    return r
end

@query(region)
@query(nation)
@query(customer)
@query(supplier)
@query(part)
@query(partsupp)
@query(order)
@query(lineitem)

# Q1
@query(
    lineitem
    :filter(shipdate <= date("1998-12-01") - DELTA*days)
    :group(returnflag, linestatus)
    :select(
        returnflag,
        linestatus,
        sum_qty => sum(lineitem.quantity),
        sum_base_price => sum(lineitem.extendedprice),
        sum_disc_price => sum(lineitem.(extendedprice*(1-discount))),
        sum_charge => sum(lineitem.(extendedprice*(1-discount)*(1+tax))),
        avg_qty => mean(lineitem.quantity),
        avg_price => mean(lineitem.extendedprice),
        avg_disc => mean(lineitem.discount),
        count_order => count(lineitem)),
    DELTA=90)

# Q2
@query(
    partsupp
    :filter(
        (part.size == SIZE) &
        contains(part.type_, TYPE) &
        (supplier.nation.region.name == REGION))
    :filter(supplycost == min(and_around(part).supplycost))
    :sort(
        supplier.acctbal:desc,
        supplier.nation.name,
        part.name,
        part.id)
    :select(
        supplier.acctbal,
        supplier.name,
        supplier.nation,
        part.id,
        part.mfgr,
        supplier.address,
        supplier.phone,
        supplier.comment)
    :take(100),
    SIZE=8, TYPE="BRASS", REGION="EUROPE")

# Q3
@query(
    order
    :filter(
        (customer.mktsegment == SEGMENT) &
        (orderdate < DATE))
    :select(
        id,
        revenue =>
            lineitem:filter(shipdate > DATE).(extendedprice*(1-discount)):sum,
        orderdate,
        shippriority)
    :sort(revenue:desc, orderdate)
    :take(10),
    SEGMENT="BUILDING", DATE=Date("1995-03-15"))

# Q4
@query(
    order
    :filter(
        DATE <= orderdate < DATE + 3_months &&
        any(lineitem.(commitdate < receiptdate)))
    :group(orderpriority)
    :select(
        orderpriority,
        order_count => count(order)),
    DATE=Date("1993-07-01"))

# Q5
@query(
    lineitem
    :filter(
        DATE <= order.orderdate < DATE + 1_year &&
        supplier.nation == order.customer.nation &&
        supplier.nation.region.name == REGION)
    :group(supplier.nation)
    :select(
        nation,
        revenue => sum(lineitem.(extendedprice*(1-discount))))
    :sort(revenue:desc),
    REGION="ASIA", DATE=Date("1994-01-01"))

# Q6
@query(
    lineitem
    :filter(
        DATE <= shipdate < DATE + 1_year &&
        DISCOUNT-0.011 <= discount <= DISCOUNT+0.011 &&
        quantity < QUANTITY)
    :group
    :select(
        revenue => sum(lineitem.(extendedprice*discount))),
    DATE=Date("1994-01-01"), DISCOUNT=0.06, QUANTITY=24)

# Q7
@query(
    lineitem
    :define(
        supp_nation => supplier.nation.name,
        cust_nation => order.customer.nation.name)
    :filter(
        date("1995-01-01") <= shipdate <= date("1996-12-31") && (
            (supp_nation == NATION1 && cust_nation == NATION2) ||
            (supp_nation == NATION2 && cust_nation == NATION1)))
    :group(
        supp_nation,
        cust_nation,
        year => year(shipdate))
    :select(
        supp_nation,
        cust_nation,
        year,
        revenue => sum(lineitem.(extendedprice*(1-discount)))),
    NATION1="FRANCE", NATION2="GERMANY")

# Q8
@query(
    lineitem
    :define(
        year => year(order.orderdate),
        supp_nation => supplier.nation.name,
        volume => extendedprice*(1-discount))
    :filter(
        part.type_ == TYPE &&
        order.customer.nation.region.name == REGION &&
        1995 <= year <= 1996)
    :group(year)
    :select(
        year,
        mkt_share =>
            sum(lineitem:filter(supp_nation == NATION).volume) / sum(lineitem.volume)),
    NATION="CANADA", REGION="AMERICA", TYPE="ECONOMY ANODIZED STEEL")

# Q9
@query(
    lineitem
    :define(
        amount =>
            extendedprice*(1-discount) - partsupp.supplycost*quantity)
    :filter(contains(part.name, COLOR))
    :group(
        nation => supplier.nation.name,
        year => year(order.orderdate))
    :sort(nation, year:desc)
    :select(
        nation,
        year,
        sum_profit => sum(lineitem.amount)),
    COLOR="green")

# Q10
@query(
    customer
    :define(
        returned =>
            order:filter(DATE <= orderdate < DATE+3_months).lineitem:filter(returnflag == "R"))
    :select(
        id,
        name,
        revenue => sum(returned.(extendedprice*(1-discount))),
        acctbal,
        nation,
        address,
        phone,
        comment)
    :sort(revenue:desc)
    :take(20),
    DATE=Date("1993-10-01"))

# Q11
@query(
    partsupp
    :filter(supplier.nation.name == NATION)
    :group(part)
    :define(value => sum(partsupp.(supplycost * availqty)))
    :filter(value > sum(and_around.value)*FRACTION)
    :select(part.id, value)
    :sort(value:desc),
    NATION="GERMANY", FRACTION=0.0001)

# Q12
@query(
    lineitem
    :define(
        high => order.orderpriority in ["1-URGENT", "2-HIGH"])
    :filter(
        shipmode in SHIPMODES &&
        shipdate < commitdate < receiptdate &&
        DATE <= receiptdate < DATE + 1_year)
    :group(shipmode)
    :select(
        shipmode,
        high_line_count => count(lineitem:filter(high)),
        low_line_count => count(lineitem:filter(!high))),
    SHIPMODES=["MAIL", "SHIP"], DATE=Date("1994-01-01"))

# Q13
@query(
    customer
    :group(
        c_count => count(order:filter(!contains(comment, WORDS))))
    :select(
        c_count,
        custdist => count(customer))
    :sort(
        custdist:desc,
        c_count:desc),
    WORDS=r".*special.*requests.*")

# Q14
@query(
    lineitem
    :define(
        volume => extendedprice*(1-discount),
        promo => startswith(part.type_, "PROMO"))
    :filter(DATE <= shipdate < DATE + 1_month)
    :group
    :select(
        promo_revenue =>
            100*sum(lineitem:filter(promo).volume)/sum(lineitem.volume)),
    DATE=Date("1995-09-01"))

# Q15
@query(
    supplier
    :define(
        total_revenue =>
            sum(lineitem:filter(DATE <= shipdate < DATE+3_months).(extendedprice*(1-discount))))
    :first(total_revenue)
    :select(
        id,
        name,
        address,
        phone,
        total_revenue),
    DATE=Date("1996-01-01"))

# Q16
@query(
    partsupp
    :filter(
        part.brand != BRAND &&
        !startswith(part.type_, TYPE) &&
        part.size in SIZES &&
        !contains(supplier.comment, r".*Customer.*Complaints.*"))
    :group(
        part.brand,
        part.type_,
        part.size)
    :select(
        brand,
        type_,
        size,
        supplier_cnt => count(unique(partsupp.supplier)))
    :sort(supplier_cnt:desc),
    BRAND="Brand#45", TYPE="MEDIUM POLISHED", SIZES=[49, 14, 23, 45, 19, 3, 36, 9])

# Q17
@query(
    lineitem
    :filter(part.brand == BRAND && part.container == CONTAINER)
    :filter(quantity < 0.2 * mean(and_around(part).quantity))
    :group
    :select(
        avg_yearly => sum(lineitem.extendedprice)/7),
    BRAND="Brand#23", CONTAINER="MED JAR")

# Q18
@query(
    order
    :select(
        customer.name,
        custkey => customer.id,
        orderkey => id,
        orderdate,
        totalprice,
        quantity => sum(lineitem.quantity))
    :filter(quantity > QUANTITY)
    :sort(
        totalprice:desc,
        orderdate)
    :take(100),
    QUANTITY=300)

# Q19
@query(
    lineitem
    :filter(
        (
            part.brand == BRAND1 &&
            part.container in ["SM CASE", "SM BOX", "SM PACK", "SM PKG"] &&
            QUANTITY1 <= quantity <= QUANTITY1+10 &&
            1 <= part.size <= 5 &&
            shipmode in ["AIR", "AIR REG"] &&
            shipinstruct == "DELIVER IN PERSON"
        )
        ||
        (
            part.brand == BRAND2 &&
            part.container in ["MED BAG", "MED BOX", "MED PKG", "MED PACK"] &&
            QUANTITY2 <= quantity <= QUANTITY2+10 &&
            1 <= part.size <= 10 &&
            shipmode in ["AIR", "AIR REG"] &&
            shipinstruct == "DELIVER IN PERSON"
        )
        ||
        (
            part.brand == BRAND3 &&
            part.container in ["LG CASE", "LG BOX", "LG PACK", "LG PKG"] &&
            QUANTITY3 <= quantity <= QUANTITY3+10 &&
            1 <= part.size <= 15 &&
            shipmode in ["AIR", "AIR REG"] &&
            shipinstruct == "DELIVER IN PERSON"
        ))
    :group
    :select(
        revenue => sum(lineitem.(extendedprice*(1-discount)))),
    QUANTITY1=1, QUANTITY2=10, QUANTITY3=20,
    BRAND1="Brand#53", BRAND2="Brand#35", BRAND3="Brand#32")

# Q20
@query(
    partsupp
    :define(
        expectedqty =>
            sum(lineitem:filter(DATE <= shipdate < DATE + 1_year).quantity))
    :filter(
        supplier.nation.name == NATION &&
        startswith(part.name, COLOR) &&
        availqty > 0.5*expectedqty)
    :group(supplier)
    :select(
        supplier.name,
        supplier.address)
    :sort(name),
    COLOR="forest", DATE=Date("1994-01-01"), NATION="CANADA")

# Q21
@query(
    lineitem
    :filter((
        order.orderstatus == "F" &&
        supplier.nation.name == NATION &&
        receiptdate > commitdate &&
        any(order.lineitem.supplier != SUPP) &&
        !any(order.lineitem.(supplier != SUPP && receiptdate > commitdate)))
        :given(SUPP => supplier))
    :group(supplier)
    :select(
        supplier.name,
        numwait => count(lineitem))
    :sort(
        numwait:desc,
        name)
    :take(100),
    NATION="SAUDI ARABIA")

# Q22
@query(
    customer
    :define(
        cntrycode => substr(phone, 1, 2))
    :filter(
        cntrycode in IS &&
        acctbal > 0_usd)
    :filter(
        !exists(order) &&
        acctbal > mean(and_around.acctbal))
    :group(cntrycode)
    :select(
        cntrycode,
        numcust => count(customer),
        totacctbal => sum(customer.acctbal)),
    IS=["13", "31", "23", "29", "30", "18", "17"])

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct PolicyDefinitionProperties <: SwaggerModel
    policyType::Any # spec type: Union{ Nothing, String } # spec name: policyType
    displayName::Any # spec type: Union{ Nothing, String } # spec name: displayName
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    policyRule::Any # spec type: Union{ Nothing, Any } # spec name: policyRule

    function PolicyDefinitionProperties(;policyType=nothing, displayName=nothing, description=nothing, policyRule=nothing)
        o = new()
        validate_property(PolicyDefinitionProperties, Symbol("policyType"), policyType)
        setfield!(o, Symbol("policyType"), policyType)
        validate_property(PolicyDefinitionProperties, Symbol("displayName"), displayName)
        setfield!(o, Symbol("displayName"), displayName)
        validate_property(PolicyDefinitionProperties, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(PolicyDefinitionProperties, Symbol("policyRule"), policyRule)
        setfield!(o, Symbol("policyRule"), policyRule)
        o
    end
end # type PolicyDefinitionProperties

const _property_map_PolicyDefinitionProperties = Dict{Symbol,Symbol}(Symbol("policyType")=>Symbol("policyType"), Symbol("displayName")=>Symbol("displayName"), Symbol("description")=>Symbol("description"), Symbol("policyRule")=>Symbol("policyRule"))
const _property_types_PolicyDefinitionProperties = Dict{Symbol,String}(Symbol("policyType")=>"String", Symbol("displayName")=>"String", Symbol("description")=>"String", Symbol("policyRule")=>"Any")
Base.propertynames(::Type{ PolicyDefinitionProperties }) = collect(keys(_property_map_PolicyDefinitionProperties))
Swagger.property_type(::Type{ PolicyDefinitionProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_PolicyDefinitionProperties[name]))}
Swagger.field_name(::Type{ PolicyDefinitionProperties }, property_name::Symbol) =  _property_map_PolicyDefinitionProperties[property_name]

const _allowed_PolicyDefinitionProperties_policyType = ["NotSpecified", "BuiltIn", "Custom"]

function check_required(o::PolicyDefinitionProperties)
    true
end

function validate_property(::Type{ PolicyDefinitionProperties }, name::Symbol, val)
    if name === Symbol("policyType")
        Swagger.validate_param(name, "PolicyDefinitionProperties", :enum, val, _allowed_PolicyDefinitionProperties_policyType)
    end
end


text:
"""
    module GeneralizedGammaDistribution

This module provides the generalized gamma distribution.

Constructors: GeneralizedGamma

Parameterizations:

Three parameterizations are supported. They are represented by types `RflexsurvParams`, `StacyParams`, and `WikipediaParams`.
`RflexsurvParams` is the default parameterization.

These types may be given as the first argument to the constructor `GeneralizedGamma`, and to `params`.
"""
module GeneralizedGammaDistribution
import Distributions

export GeneralizedGamma
export RflexsurvParams, StacyParams, WikipediaParams

# Use parameterization of R flexsurv implementation
"""
    GeneralizedGamma(μ=1.0, σ=1.0, Q=1.0)

The *Generalized gamma distribution* has probability density
function (in the Wikipedia parameterization)

```math
f(x; a, d, p) = \\frac{(p/a^d) x^{d-1} e^{-(x/a)^p}}{\\Gamma(d/p)}
\\quad x > 0
```

```julia
params(RflexsurvParams, d)  # Get the parameters, i.e. (μ, σ, Q)
params(d)                   # Get the parameters, i.e. (μ, σ, Q)
params(WikipediaParams, d)  # Get the parameters, (a, d, p)
params(StacyParams, d)      # Get the parameters, (b, d, p)
```

External links

* [Generalized gamma distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalized_gamma_distribution)
"""
struct GeneralizedGamma{T<:Real} <: Distributions.ContinuousUnivariateDistribution
    μ::T
    σ::T
    Q::T
    gamma_dist::Distributions.Gamma  # gamma distribution used to compute random samples
end

"""
    RflexsurvParams

Parameterization of the generalized gamma distribution used in the R flexsurve implementation.
This is the default parameterization
"""
struct RflexsurvParams # <: AbstractParameterization
end

"""
    WikipediaParams

Parameterization of the generalized gamma distribution found in the Wikipedia article
"""
struct WikipediaParams # <: AbstractParameterization
end

"""
    StacyParams

The Stacy Parameterization of the generalized gamma distribution. This parameterization was
used in the older dgengamma.orig implementation.
"""
struct StacyParams # <: AbstractParameterization
end

"""
    GeneralizedGamma(::Type{RflexsurvParams}, μ, σ, Q)
    GeneralizedGamma(μ, σ, Q)

The `generalized gamma distribution` parameterized as in the R `flexsurv`
implementation[^1].

[^1]: Prentice, R. L. (1974). A log gamma model and its maximum likelihood estimation. Biometrika 61(3):539-544
"""
function GeneralizedGamma(::Type{RflexsurvParams}, μ=1.0, σ=1.0, Q=1.0)
    gamma_dist = Distributions.Gamma(1 / Q^2, 1)
    return GeneralizedGamma(μ, σ, Q, gamma_dist)
end
GeneralizedGamma(μ::Real=1.0, σ::Real=1.0, Q::Real=1.0) = GeneralizedGamma(RflexsurvParams, μ, σ, Q)

# Parameterization from wikipedia page
# PDF is: \frac{p/a^d}{\Gamma(d/p)} x^{d-1}e^{-(x/a)^p}
# The R page for gengamma flexsurv says that Q < 0 is allowed.
# The wikipedia page does not have abs(d) and abs(p), but this works, because the -i π cancels
# Must have d > 0, p > 0  or  d < 0, p < 0
"""
    GeneralizedGamma(::Type{WikipediaParams}, a=1.0, d=1.0, p=1.0)

The `generalized gamma distribution` parameterized as on the Wikipedia page.
Requires either `d > 0, p > 0` or `d < 0, p < 0`.
"""
function GeneralizedGamma(::Type{WikipediaParams}, a=1.0, d=1.0, p=1.0)
    μ = log(a) + (log(abs(d)) - log(abs(p))) / p
    σ = 1 / sqrt(p * d)
    Q = sqrt(p / d) * sign(p)
    return GeneralizedGamma(μ, σ, Q)
end

"""
    GeneralizedGamma(::Type{StacyParams}, b=1.0, d=1.0, p=1.0)

The `generalized gamma distribution` with the Stacy parameterization. This follows
the older R convention `dgengamma.orig` [^2].

[^2]: Stacy, E. W. (1962). A generalization of the gamma distribution. Annals of Mathematical Statistics 33:1187-92.
"""
function GeneralizedGamma(::Type{StacyParams}, b=1.0, d=1.0, p=1.0)
    a = b^(-1 / p)
    return GeneralizedGamma(WikipediaParams, a, d, p)
end

@inline Distributions.partype(d::GeneralizedGamma{T}) where {T<:Real} = T

# Use the algorithm given in the documentation pages for the R package flexsurv.
function Distributions.rand(p::GeneralizedGamma)
    Qs = p.Q^2
    gamma_deviate = rand(p.gamma_dist)  # only saves 10 or so percent time.
    w = log(Qs * gamma_deviate) / p.Q
    return exp(p.μ + p.σ * w)
end

#  mean in wikipedia parameterization is
#  a * gamma((d+1)/p)/gamma(d/p)
function Distributions.mean(p::GeneralizedGamma)
    Qs = p.Q^2
    iQs = 1 / Qs
    a = Qs^(p.σ / p.Q) * exp(p.μ)
    return a * Distributions.gamma(iQs + p.σ / p.Q) / Distributions.gamma(iQs)
end

function Distributions.pdf(p::GeneralizedGamma, x::Real)
    (d, p, a) = Distributions.params(WikipediaParams, p)
    return (p / a^d)/ Distributions.gamma(d / p) * x^(d - 1) * exp(-(x / a)^p)
end

function Distributions.mode(p::GeneralizedGamma{T}) where T
    (d, p, a) = Distributions.params(WikipediaParams, p)
    return d <= 1 ? zero(T) : a * ((d - 1) / p)^(1/p)
end

Distributions.support(::GeneralizedGamma{T}) where {T} = Distributions.RealInterval(zero(T), convert(T, Inf))

function Distributions.var(dt::GeneralizedGamma)
    (d, p, a) = Distributions.params(WikipediaParams, dt)
    d2 = Distributions.gamma(d / p)
    return a^2 * ((Distributions.gamma((d + 2) / p) / d2) -
                  (Distributions.gamma((d + 1) / p) / d2)^2)
end

# How do we include the lower incomplete gamma function ? How is it done in Distributions ?
# It is buried somewhere in the layers of code.
# function Distributions.cdf(p::GeneralizedGamma, x::Real)
#     (d, p, a) = Distributions.params(WikipediaParams, p)
#     num = lower_incomplete_gamma(d/p, (x/a)^p)
#     return num /  Distributions.gamma(d / p)
# end

"""
    params(::Type{WikipediaParams}, dt::GeneralizedGamma)

Return parameters of the `generalized Gamma distribution` following
the convention on the Wikipedia page.
"""
function Distributions.params(::Type{WikipediaParams}, dt::GeneralizedGamma)
    d = 1 / (dt.σ * dt.Q)
    p = (dt.Q) / (dt.σ)
    a = abs(dt.Q)^(2 / p) * exp(dt.μ)
    return (a, d, p)
end

"""
   params(::Type{StacyParams}, dt::GeneralizedGamma)

Return parameters of the `generalized Gamma distribution` following the parameterization
of the older R implementation `dgengamma.orig` [^2].

[^2]: Stacy, E. W. (1962). A generalization of the gamma distribution. Annals of Mathematical Statistics 33:1187-92.
"""
function Distributions.params(::Type{StacyParams}, dt::GeneralizedGamma)
    (a, d, p) = Distributions.params(WikipediaParams, dt)
    b = a^(-p)
    return (b, d, p)
end

"""
   params(::Type{RflexsurvParams}, dt::GeneralizedGamma)

Return parameters of the `generalized Gamma distribution` following
the convention of the R `flexsurv` function.
"""
function Distributions.params(::Type{RflexsurvParams}, d::GeneralizedGamma)
    return (d.μ, d.σ, d.Q)
end

Distributions.params(d::GeneralizedGamma) = Distributions.params(RflexsurvParams, d)

end # module


text:
import StaticLint: hasfile, canloadfile, loadfile, setfile, getfile, getsymbolserver
import StaticLint: getpath, setpath, getroot, setroot, getcst, setcst, scopepass, getserver, setserver
hasfile(server::LanguageServerInstance, path::String) = hasdocument(server, URI2(filepath2uri(path)))
canloadfile(server::LanguageServerInstance, path::String) = isfile(path)
function loadfile(server::LanguageServerInstance, path::String)
    source = read(path, String)
    uri = filepath2uri(path)
    doc = Document(uri, source, true, server)
    StaticLint.setfile(server, path, doc)
end
setfile(server::LanguageServerInstance, path::String, x::Document) = setdocument!(server, URI2(filepath2uri(path)), x)
getfile(server::LanguageServerInstance, path::String) = getdocument(server, URI2(filepath2uri(path)))
getsymbolserver(server::LanguageServerInstance) = server.symbol_store

getpath(d::Document) = d.path
function setpath(d::Document, path::String)
    d.path = path
    return d
end

getroot(d::Document) = d.root
function setroot(d::Document, root::Document)
    d.root = root
    return d
end

getcst(d::Document) = d.cst
function setcst(d::Document, cst::EXPR)
    d.cst = cst
    return d
end

getserver(file::Document) = file.server
function setserver(file::Document, server::LanguageServerInstance)
    file.server = server
    return file
end


text:
using Test
using LinearAlgebra
using OrdinaryDiffEq, DiffEqDevTools

# Since non-FSAL methods are FSAL'ed artificially in OrdinaryDiffEq.jl,
# we need to exclude this last FSAL stage for the comparison with
# tableaus.
function coefficients_as_in_tableau(A, b, c, tab)
  if size(A) == size(tab.A)
    AA = A
    bb = b
    cc = c
  else
    AA = A[1:end-1, 1:end-1]
    bb = b[1:end-1]
    cc = c[1:end-1]
  end
  AA, bb, cc
end

function ode_tableau_tests(T)
  let erk = Heun()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructHeun()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = Euler()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructEuler()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = RK4()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructRK4()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = SSPRK22()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructSSPRK22()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = SSPRK33()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructSSPRK33()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = SSPRK432()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructSSPRK43()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = SSPRK104()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructSSPRK104()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = Tsit5()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructTsitouras5()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = BS3()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructBogakiShampine3()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = BS5()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructBogakiShampine5()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = Vern7()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructVerner7()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end

  let erk = Vern8()
    A, b, c = deduce_Butcher_tableau(erk)
    tab = constructVerner8()
    AA, bb, cc = coefficients_as_in_tableau(A, b, c, tab)
    @test AA ≈ tab.A
    @test bb ≈ tab.α
    @test cc ≈ tab.c
  end
end

ode_tableau_tests(Float32)
ode_tableau_tests(Float64)


text:
############################################################################################
#Markov Kernel
markov_latent = Int32.(rand(_rng, Distributions.Categorical(3), N))
markov_data = randn(_rng, Float16, size(markov_latent, 1))
markov_param = (;
    μ=Param([-.1, 0.0, 0.1], fill(Normal(0.0, 10), 3)),
    σ=Param([5.0, 3.0, 2.0], fill(truncated(Normal(3.0, 10.0), 0.0, 10.0), 3)),
    p=Param(
        [[0.2, 0.6, 0.2], [0.2, 0.6, 0.2], [0.2, 0.6, 0.2]], [Dirichlet(3, 3) for i in 1:3]
    ),
    latent=Param(markov_latent, Fixed()),
)
struct HMM <: ModelName end
hmm = ModelWrapper(HMM(), markov_param)
markov_objective = Objective(hmm, markov_data, :latent)

function ModelWrappers.generate(_rng::Random.AbstractRNG, objective::Objective{<:ModelWrapper{HMM}})
    @unpack model, data = objective
    @unpack μ, σ, latent = model.val
    return rand(_rng, Normal(Float32(2.), Float32(3.)))
end

function ModelWrappers.predict(_rng::Random.AbstractRNG, objective::Objective{<:ModelWrapper{HMM}})
    @unpack model, data = objective
    @unpack μ, σ, latent = model.val
    #NOTE: Check if correct prediction is applied even if standard predict is overloaded
	return Int(3) #rand(_rng, Normal(μ[latent[end]], σ[latent[end]]))
end

function get_dynamics(model::ModelWrapper{<:HMM}, θ)
    @unpack μ, σ, p = θ
    dynamicsᵉ = [Normal(μ[iter], σ[iter]) for iter in eachindex(μ)]
    dynamicsˢ = [Categorical(p[iter]) for iter in eachindex(μ)]
    return dynamicsᵉ, dynamicsˢ
end
function BaytesFilters.dynamics(objective::Objective{<:ModelWrapper{<:HMM}})
    @unpack model, data = objective
    dynamicsᵉ, dynamicsˢ = get_dynamics(model, model.val)

    initialˢ = Categorical(fill(1 / length(dynamicsᵉ), length(dynamicsᵉ)))
    transition(particles, iter) = dynamicsˢ[particles[iter - 1]]
    evidence(particles, iter) = dynamicsᵉ[particles[iter]]

    return Markov(initialˢ, transition, evidence)
end
dynamics(markov_objective)

############################################################################################
"Forward Filter HMM"
function filter_forward(objective::Objective{<:ModelWrapper{HMM}})
    @unpack model, data = objective
## Map Parameter to observation and state probabilities
    @unpack p = model.val
    dynamicsᵉ, _ = get_dynamics(model, model.val)
#!NOTE: Working straight with parameter instead of distributions here as easier to implement
    dynamicsˢ           = transpose( reduce(hcat, p ) ) #[ Categorical(p[iter]) for iter in eachindex(μ) ]
    initialˢ            = get_stationary( dynamicsˢ )
    structure = BaytesCore.ByRows()
## Assign Log likelihood
    ℓℒᵈᵃᵗᵃ = zeros(Float64, size(data,1), size(dynamicsˢ, 1) )
    Base.Threads.@threads for state in Base.OneTo( size(ℓℒᵈᵃᵗᵃ, 2) )
    for iter in Base.OneTo( size(ℓℒᵈᵃᵗᵃ, 1) )
            ℓℒᵈᵃᵗᵃ[iter, state] += logpdf( dynamicsᵉ[state], grab(data, iter, structure) )
        end
    end
## Initialize
    Nstates = size(dynamicsˢ, 1)
    α = zeros( size(ℓℒᵈᵃᵗᵃ) ) # α represents P( sₜ | e₁:ₜ ) for each t, which is numerically more stable than classical forward probabilities p(sₜ, e₁:ₜ)
    c = vec_maximum( @view(ℓℒᵈᵃᵗᵃ[1,:]) ) # Used for stable logsum() calculations for incremental log likelihood addition, i.e.: log( exp(x) + exp(y) ) == x + log(1 + exp(y -x ) ) for y >= x.
    ℓℒ = 0.0 #Log likelihood container
## Calculate initial probabilities
    for stateₜ in Base.OneTo(Nstates)
         α[1,stateₜ] += initialˢ[stateₜ] * exp(ℓℒᵈᵃᵗᵃ[1,stateₜ]-c) #Calculate initial p(s₁, e₁) ∝ P(s₁) * p(e₁ | s₁)
    end
## Normalize P( s₁, e₁ ) and return normalizing constant P(e₁), which can be used to calculate ℓℒ = P(e₁:ₜ) = p(e₁) ∏_k p(eₖ | e₁:ₖ₋₁)
    norm = sum( @view(α[1,:]) )
    α[1,:] ./= norm
    ℓℒ += log(norm) + c# log(norm) = p(eₜ | e₁:ₜ₋₁), c is constant that is added back after removing on top for numerical stability
## Loop through sequence
    for t = 2:size(α, 1)
        c = vec_maximum( @view(ℓℒᵈᵃᵗᵃ[t,:]) ) # Used for stable logsum() calculations, i.e.: log( exp(x) + exp(y) ) == x + log(1 + exp(y -x ) ) for y >= x.
        ## Calculate ∑_sₜ₋₁ P( sₜ | sₜ₋₁) * P(sₜ₋₁ | e₁:ₜ₋₁) - we sum over sₜ₋₁ - states are per row
        for stateₜ in Base.OneTo(Nstates)
            for stateₜ₋₁ in Base.OneTo(Nstates)
                α[t,stateₜ] += dynamicsˢ[stateₜ₋₁, stateₜ] * α[t-1,stateₜ₋₁] # * for both log and non-log version as inside log( sum(probability(...)))
            end
            ## Then multiply with ℒᵈᵃᵗᵃ corrector P( eₜ | sₜ )
            α[t,stateₜ] *= exp(ℓℒᵈᵃᵗᵃ[t,stateₜ]-c) # - c for higher numerical stability, will be added back to likelihood increment below
        end
        ## Normalize α and obtain P( eₜ | e₁:ₜ₋₁)
        norm = sum( @view(α[t,:]) )
        α[t,:] ./= norm
        ## Add normalizing constant p(eₜ | e₁:ₜ₋₁) to likelihood term
        ℓℒ += log(norm)+c
    end
    return (α, ℓℒ)
end
filter_forward(markov_objective)

############################################################################################
function ModelWrappers.simulate(rng::Random.AbstractRNG, model::ModelWrapper{HMM}; Nsamples = 1000)
    dynamicsᵉ, dynamicsˢ = get_dynamics(model, model.val)
    latentⁿᵉʷ = Vector{eltype(model.val.latent)}(undef, Nsamples)
    observedⁿᵉʷ = zeros(Nsamples)

    stateₜ = rand(rng, dynamicsˢ[ rand(1:length(dynamicsˢ) ) ] )
    latentⁿᵉʷ[1] = stateₜ
    observedⁿᵉʷ[1] = rand(rng, dynamicsᵉ[stateₜ])

    for iter in 2:size(observedⁿᵉʷ,1)
            stateₜ = rand(rng, dynamicsˢ[ stateₜ ]) #stateₜ for t-1 overwritten
            latentⁿᵉʷ[iter] = stateₜ
            observedⁿᵉʷ[iter] = rand(rng, dynamicsᵉ[stateₜ])
    end
    return observedⁿᵉʷ, latentⁿᵉʷ
end


text:
using Multibreak

# This file serves the dual purpose of both testing the package and
# being a tutorial to the functionality.

# The first example showcases the semantics of multiple
# `break`/`continue` in ordinary Julia code. Effectively only the
# first `break`/`continue` is executed and additional ones are simply
# unreachable code.
@testset MBTestSet "standard Julia" begin
    out = []
    for i = 1:5
        for j = 1:3
            if i == 1 && j == 1
                continue
            elseif i == 2 && j == 2
                break
            elseif i == 3 && j == 2
                break; continue
            elseif i == 4 && j == 2
                break; break
            end
            push!(out, (i, j))
        end
        push!(out, i)
    end

    @test out == [(1, 2), (1, 3), 1,
                  (2, 1), 2,
                  (3, 1), 3,
                  (4, 1), 4,
                  (5, 1), (5, 2), (5, 3), 5]
end

# With the `@multibreak` macro, the semantics of multiple
# `break`/`continue` is changed to a sequence of `break`/`continue`
# statements outwards through the enclosing loops. Once a `continue`
# is reached, the process ends. Thus there can be any number of
# `break`, optionally followed by a `continue`. Obviously there cannot
# be more `break`/`continue` than nested loop levels, however.
@testset MBTestSet "multibreak macro" begin
    out = []
    @multibreak begin
        for i = 1:5
            for j = 1:3
                if i == 1 && j == 1
                    continue
                elseif i == 2 && j == 2
                    break
                elseif i == 3 && j == 2
                    break; continue
                elseif i == 4 && j == 2
                    break; break
                end
                push!(out, (i, j))
            end
            push!(out, i)
        end
    end

    @test out == [(1, 2), (1, 3), 1,
                  (2, 1), 2,
                  (3, 1),
                  (4, 1)]
end

# This shows how the previous example can be coded with the `goto` and
# `label` macros instead of the `@multibreak` macro. This is in fact
# how the `@multibreak` macro works internally and should make the
# semantics clear.
@testset MBTestSet "goto and label" begin
    out = []
    for i = 1:5
        for j = 1:3
            if i == 1 && j == 1
                @goto loop2continue
            elseif i == 2 && j == 2
                @goto loop2break
            elseif i == 3 && j == 2
                @goto loop1continue
            elseif i == 4 && j == 2
                @goto loop1break
            end
            push!(out, (i, j))
            @label loop2continue
        end
        @label loop2break
        push!(out, i)
        @label loop1continue
    end
    @label loop1break

    @test out == [(1, 2), (1, 3), 1,
                  (2, 1), 2,
                  (3, 1),
                  (4, 1)]
end

# `@multibreak` can be applied directly to the outer loop statement
# instead of to an enclosing block as was done earlier. If the nested
# loop is the only thing in the block, this is only a stylistic
# difference.
@testset MBTestSet "multibreak applied to outer for" begin
    out = []
    @multibreak for i = 1:5
        for j = 1:3
            if i == 1 && j == 1
                continue
            elseif i == 2 && j == 2
                break
            elseif i == 3 && j == 2
                break; continue
            elseif i == 4 && j == 2
                break; break
            end
            push!(out, (i, j))
        end
        push!(out, i)
    end

    @test out == [(1, 2), (1, 3), 1,
                  (2, 1), 2,
                  (3, 1),
                  (4, 1)]
end

# If preferred, `@multibreak` can be applied to an entire function.
@testset MBTestSet "multibreak applied to function" begin
    @multibreak function f()
        out = []
        for i = 1:5
            for j = 1:3
                if i == 1 && j == 1
                    continue
                elseif i == 2 && j == 2
                    break
                elseif i == 3 && j == 2
                    break; continue
                elseif i == 4 && j == 2
                    break; break
                end
                push!(out, (i, j))
            end
            push!(out, i)
        end
        return out
    end

    @test f() == [(1, 2), (1, 3), 1,
                  (2, 1), 2,
                  (3, 1),
                  (4, 1)]
end

# When `@multibreak` is applied to an entire function, multiple
# `break`/`continue` is available in all loops in the function.
@testset MBTestSet "multibreak applied to function with two loops" begin
    @multibreak function f()
        m = 0
        for i = 1:5
            for j = 1:3
                m += 1
                if i + j > 5
                    break; continue
                end
            end
        end

        n = 0
        for i = 1:5
            for j = 1:3
                n += 1
                if i + j > 5
                    break; break
                end
            end
        end

        return m, n
    end

    @test f() == (12, 9)
end

# You are not allowed to break out of more loops than there are.
@testset MBTestSet "too many break" begin
    @multibreak function f()
        for i = 1:5
            break; break
        end
    end

    @test_throws ErrorException f()
end

# Neither continue after breaking out of all the loops.
@testset MBTestSet "too many break and continue" begin
    @multibreak function f()
        for i = 1:5
            for j = 1:3
                if i + j > 5
                    break; break; continue
                end
            end
        end
    end

    @test_throws ErrorException f()
end

# It doesn't make sense to break out of a loop after having continued it.
@testset MBTestSet "continue before break" begin
    @multibreak function f()
        for i = 1:5
            for j = 1:3
                if i + j > 5
                    continue; break
                end
            end
        end
    end

    @test_throws ErrorException f()
end

# Neither can you continue repeatedly.
@testset MBTestSet "multiple continue" begin
    @multibreak function f()
        for i = 1:5
            for j = 1:3
                for k = 1:2
                    if i + j + k > 5
                        break; continue; continue
                    end
                end
            end
        end
    end

    @test_throws ErrorException f()
end

# For clarity, breaking cannot pass function call boundaries.
@testset MBTestSet "cannot break through function barrier" begin
    @multibreak function f()
        for i = 1:5
            for j = 1:3
                break; break; break
            end
        end
    end

    # Calling `f` gives an error rather than breaking the `k` loop.
    for k = 1:2
        @test_throws ErrorException f()
    end
end

# Note, in sufficiently simple cases, nested loops can be written as a
# single for statement with multiple variables. In that case a single
# `break` breaks out of all the loops whereas a single `continue`
# continues the inner loop. If that is enough, there is no need to
# involve the `@multibreak` macro.
@testset MBTestSet "single loop with multiple variables" begin
    n = 0
    for i = 1:5, j = 1:3
        n += 1
        if i + j > 5
            break
        end
    end
    @test n == 9

    n = 0
    for i = 1:5, j = 1:3
        n += 1
        if i + j > 5
            continue
        end
    end
    @test n == 15
end

# Naturally for loops with multiple variables can in turn be nested,
# in which case the `@multibreak` macro can be utilized.
@testset MBTestSet "nested loops with multiple variables" begin
    n = 0
    @multibreak begin
        for i = 1:3, j = 1:3
            for k = 1:3, l = 1:3
                n += 1
                if i + j + k + l > 6
                    break; break
                end
                if i + j + k + l > 5
                    break; continue
                end
            end
        end
    end
    @test n == 10
end

# So there's also a kind of loop called `while`, huh? We had better try
# that as well.
@testset MBTestSet "nested while loops" begin
    n = 0

    @multibreak begin
        while true
            while true
                n += 1
                break; break
            end
        end
    end
    @test n == 1

    @multibreak begin
        i = 0
        while i <= 5
            i += 1
            j = 1
            while j <= 3
                n += 1
                if i + j > 4
                    break; continue
                end
                j += 1
            end
            i += 1
        end
    end
    @test n == 9
end

# `for` and `while` loops can be mixed.
@testset MBTestSet "mixed for and while loops" begin
    n = 0

    @multibreak begin
        i = 1
        while i <= 3
            i += 1
            for j = i:5
                n += 1
                if i + j > 7
                    break; break
                end
                if i + j > 5
                    break; continue
                end
            end
        end
    end
    @test n == 5
end

# A word of warning. The `@multibreak` macro also transforms single
# `break` and `continue` to `@goto` and `@label`. A side effect of the
# implementation is that dead code can come alive. The macro could be
# refined to place the `@goto` at the position of the first
# `break`/`continue` rather than at the end of the block but it's not
# really worth the added complexity. Just don't do this in code where
# you need the `@multibreak`.
@testset MBTestSet "zombie code" begin
    I_am_a_zombie = false

    @multibreak begin
        while true
            break
            I_am_a_zombie = true
        end
    end
    @test I_am_a_zombie
end


text:
using Documenter, SARSOP

makedocs(
    modules = [SARSOP],
    format = Documenter.HTML(),
    sitename = "SARSOP.jl"
)

deploydocs(
    repo = "github.com/JuliaPOMDP/SARSOP.jl.git",
)


text:
module MouseButtons

export MouseButton

@enum MouseButton::UInt8 begin
  LEFT = 1
  MIDDLE = 2
  RIGHT = 3
  WHEEL_UP = 4
  WHEEL_DOWN = 5
end

end

using .MouseButtons


text:
export Normalize

"""
$(TYPEDEF)

Represents a Normalization operation.
"""
@auto_hash_equals struct Normalize <: Layer
    mean::Array{Real,1}
    std::Array{Real,1}
end

function Base.show(io::IO, p::Normalize)
    print(io, "Normalize(means: $(p.mean), stds: $(p.std))")
end

function apply(p::Normalize, x::Array{<:JuMPReal})
    padded_shape = (ones(Int, ndims(x) - 1)..., length(p.mean))
    m = reshape(p.mean, padded_shape)
    s = reshape(p.std, padded_shape)
    output = (x .- m) ./ s
    return output
end

(p::Normalize)(x::Array{<:JuMPReal}) = apply(p, x)


text:
module AVLTrees

import Base: iterate, haskey, getkey, getindex, setindex!, length,
            eltype, isempty, insert!, popfirst!, insert!, delete!
            print, show, firstindex, pop!, popfirst!

include("node.jl")
include("tree.jl")
include("set.jl")

export AVLTree, AVLSet, findkey
end # module


text:
#=
    %FILE%

  Created by %USER% on %YEAR%/%MONTH%/%DAY%.
  Copyright (c) %YEAR%年 %USER% All rights reserved.
=#

%HERE%


text:
module C
c = 3
end

module C2
module SubC
c2 = 3.5
end
end


text:
module TestConfig

import ..repl_output
using AMGX, Defer, Test
using AMGX: Config, AMGXException

@testset "Config" begin
    @scope @testset "Dict" begin
        cfg = @! Config()
        @test occursin("uninitialized", repl_output(cfg))
        d = Dict("max_levels" => 10)
        AMGX.create!(cfg, d)
        str = sprint((io, cfg) -> show(io, MIME("text/plain"), cfg), cfg)
        @test !occursin("uninitialized", repl_output(cfg))
        @test occursin("@", repl_output(cfg))

        d = Dict("max_lovels" => 10)
        @test_throws AMGXException("Incorrect amgx configuration provided.") Config(d)
    end
end

end # module