text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

module equivRotations3D

using StaticArrays
using LinearAlgebra: norm, rank, svd, Diagonal
using ACE.SphericalHarmonics: index_y
using Combinatorics: permutations

export ClebschGordan, IRot3DCoeffs, ERot3DCoeffs, ri_basis, rpi_basis

"""
`ClebschGordan: ` storing precomputed Clebsch-Gordan coefficients; see
`?clebschgordan` for the convention that is use.
"""
struct ClebschGordan{T}
	vals::Dict{Tuple{Int, Int, Int, Int, Int, Int}, T}
end

abstract type Rot3DCoeffs{T} end

"""
`ERot3DCoeffs: ` storing recursively precomputed matrix-valued coefficients for a
rotation-equivariant basis.
"""
struct ERot3DCoeffs{T} <: Rot3DCoeffs{T}
   vals::Vector{Dict}
   cg::ClebschGordan{T}
end

"""
`IRot3DCoeffs: ` storing recursively precomputed coefficients for a
rotation-invariant basis.
"""
struct IRot3DCoeffs{T} <: Rot3DCoeffs{T}
   vals::Vector{Dict}
   cg::ClebschGordan{T}
end

# -----------------------------------
# iterating over an m collection
# -----------------------------------

_mvec(::CartesianIndex{0}) = SVector(Int(0))

_mvec(mpre::CartesianIndex) = SVector(Tuple(mpre)..., - sum(Tuple(mpre)))

struct MRange{N, T2}
   ll::SVector{N, Int}
   cartrg::T2
end

Base.length(mr::MRange) = sum(_->1, _mrange(mr.ll))

"""
Given an l-vector `ll` iterate over all combinations of `mm` vectors  of
the same length such that `sum(mm) == 0`
"""
_mrange(ll) = MRange(ll, Iterators.Stateful(
                     CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)-1))))

function Base.iterate(mr::MRange{1}, args...)
   if isempty(mr.cartrg)
      return nothing
   end
   while !isempty(mr.cartrg)
      popfirst!(mr.cartrg)
   end
   return SVector{1, Int}(0), nothing
end

function Base.iterate(mr::MRange, args...)
   while true
      if isempty(mr.cartrg)
         return nothing
      end
      mpre = popfirst!(mr.cartrg)
      if abs(sum(mpre.I)) <= mr.ll[end]
         return _mvec(mpre), nothing
      end
   end
   error("we should never be here")
end

"""
Given an l-vector `ll` iterate over all combinations of `mm` vectors  of
the same length such that `abs(sum(mm)) <= 1`
"""

struct MRange1{N, T2}
   ll::SVector{N, Int}
   cartrg::T2
end

Base.length(mr::MRange1) = sum(_->1, _mrange1(mr.ll))

_mrange1(ll) = MRange1(ll, Iterators.Stateful(
               filter((x) -> abs(sum(x))<= 1, Tuple.(CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)))))
                     ))

#Iterators.product(1:3, _mrange1(ll))

function Base.iterate(mr::MRange1, args...)
   while true
      if isempty(mr.cartrg)
         return nothing
      end
      mpre = popfirst!(mr.cartrg)
      return SVector(mpre), nothing
   end
end

# ----------------------------------------------------------------------
#     ClebschGordan code
# ----------------------------------------------------------------------

cg_conditions(j1,m1, j2,m2, J,M) =
	cg_l_condition(j1, j2, J)   &&
	cg_m_condition(m1, m2, M)   &&
	(abs(m1) <= j1) && (abs(m2) <= j2) && (abs(M) <= J)

cg_l_condition(j1, j2, J) = (abs(j1-j2) <= J <= j1 + j2)

cg_m_condition(m1, m2, M) = (M == m1 + m2)

"""
`clebschgordan(j1, m1, j2, m2, J, M, T=Float64)` :

A reference implementation of Clebsch-Gordon coefficients based on

https://hal.inria.fr/hal-01851097/document
Equation (4-6)

This heavily uses BigInt and BigFloat and should therefore not be employed
for performance critical tasks, but only precomputation.

The ordering of parameters corresponds to the following convention:
```
clebschgordan(j1, m1, j2, m2, J, M) = C_{j1m1j2m2}^{JM}
```
where
```
   D_{m1k1}^{l1} D_{m2k2}^{l2}}
	=
	∑_j  C_{l1m1l2m2}^{j(m1+m2)} C_{l1k1l2k2}^{j2(k1+k2)} D_{(m1+m2)(k1+k2)}^{j}
```
"""
function clebschgordan(j1, m1, j2, m2, J, M, T=Float64)
	if !cg_conditions(j1, m1, j2, m2, J, M)
		return zero(T)
	end

   N = (2*J+1) *
       factorial(big(j1+m1)) * factorial(big(j1-m1)) *
       factorial(big(j2+m2)) * factorial(big(j2-m2)) *
       factorial(big(J+M)) * factorial(big(J-M)) /
       factorial(big( j1+j2-J)) /
       factorial(big( j1-j2+J)) /
       factorial(big(-j1+j2+J)) /
       factorial(big(j1+j2+J+1))

   G = big(0)
   # 0 ≦ k ≦ j1+j2-J
   # 0 ≤ j1-m1-k ≤ j1-j2+J   <=>   j2-J-m1 ≤ k ≤ j1-m1
   # 0 ≤ j2+m2-k ≤ -j1+j2+J  <=>   j1-J+m2 ≤ k ≤ j2+m2
   lb = (0, j2-J-m1, j1-J+m2)
   ub = (j1+j2-J, j1-m1, j2+m2)
   for k in maximum(lb):minimum(ub)
      bk = big(k)
      G += (-1)^k *
           binomial(big( j1+j2-J), big(k)) *
           binomial(big( j1-j2+J), big(j1-m1-k)) *
           binomial(big(-j1+j2+J), big(j2+m2-k))
   end

   return T(sqrt(N) * G)
end

ClebschGordan(T=Float64) =
	ClebschGordan{T}(Dict{Tuple{Int,Int,Int,Int,Int,Int}, T}())

_cg_key(j1, m1, j2, m2, J, M) = (j1, m1, j2, m2, J, M)
	# Int.((index_y(j1,m1), index_y(j2,m2), index_y(J,M)))

function (cg::ClebschGordan{T})(j1, m1, j2, m2, J, M) where {T}
	if !cg_conditions(j1,m1, j2,m2, J,M)
		return zero(T)
	end
	key = _cg_key(j1, m1, j2, m2, J, M)
	if haskey(cg.vals, key)
		return cg.vals[key]
	end
	val = clebschgordan(j1, m1, j2, m2, J, M, T)
	cg.vals[key] = val
	return val
end

# ----------------------------------------------------------------------
#     IRot3DCoeffs code: generalized cg coefficients
#
#  Note: in this section kk is a tuple of m-values, it is not
#        related to the k index in the 1-p basis (or radial basis)
# ----------------------------------------------------------------------

dicttype(A::Rot3DCoeffs,N::Integer) = dicttype(A::Rot3DCoeffs,Val(N))

dicttype(A::IRot3DCoeffs,::Val{N}) where {N} =
   Dict{Tuple{SVector{N,Int}, SVector{N,Int}, SVector{N,Int}}, Float64}

dicttype(A::ERot3DCoeffs,::Val{N}) where {N} =
      Dict{Tuple{SVector{N,Int}, SVector{N,Int}, SVector{N,Int}}, SMatrix{3,3,ComplexF64,9}}

IRot3DCoeffs(T=Float64) = IRot3DCoeffs(Dict[], ClebschGordan(T))
ERot3DCoeffs(T=Float64) = ERot3DCoeffs(Dict[], ClebschGordan(T))

function get_vals(A::Rot3DCoeffs, valN::Val{N}) where {N}
	if length(A.vals) < N
		for n = length(A.vals)+1:N
			push!(A.vals, dicttype(A,n)())
		end
	end
   	return A.vals[N]::dicttype(A,valN)
end

_key(ll::StaticVector{N}, mm::StaticVector{N}, kk::StaticVector{N}) where {N} =
      (SVector{N, Int}(ll), SVector{N, Int}(mm), SVector{N, Int}(kk))

function (A::IRot3DCoeffs{T})(ll::StaticVector{N},
                            mm::StaticVector{N},
                            kk::StaticVector{N}) where {T, N}

   if       sum(mm) != 0 ||
            sum(kk) != 0 ||
            !all(abs.(mm) .<= ll) ||
            !all(abs.(kk) .<= ll)
      return T(0)
   end
   vals = get_vals(A, Val(N))  # this should infer the type!
   key = _key(ll, mm, kk)
   if haskey(vals, key)
      val  = vals[key]
   else
      val = _compute_val(A, key...)
      vals[key] = val
   end
   return val
end

# the recursion has two steps so we need to define the
# coupling coefficients for N = 1, 2
# TODO: actually this seems false; it is only one recursion step, and a bit
#       or reshuffling should allow us to get rid of the {N = 2} case.

function (A::IRot3DCoeffs{T})(ll::StaticVector{1},
                            mm::StaticVector{1},
                            kk::StaticVector{1}) where {T}
   if ll[1] == mm[1] == kk[1] == 0
      return T(8 * pi^2)
   else
      return T(0)
   end
end

function (A::ERot3DCoeffs{T})(ll::StaticVector{N},
                            mm::StaticVector{N},
                            kk::StaticVector{N}) where {T, N}
   if       abs(sum(mm)) > 1 ||
            abs(sum(kk)) > 1 ||
            !all(abs.(mm) .<= ll) ||
            !all(abs.(kk) .<= ll)
      return get0val(A)
   end
   vals = get_vals(A, Val(N))  # this should infer the type!
   key = _key(ll, mm, kk)
   if haskey(vals, key)
      val  = vals[key]
   else
      val = _compute_val(A, key...)
      vals[key] = val
   end
   return val
end

"""
rmatrices[(m,mu)] = int_{SO(3)} Q^T D^1_{mu,m}(Q) dQ
"""

rmatrices = Dict(
(-1,-1) => SMatrix{3, 3, ComplexF64, 9}(1/6, 1im/6, 0, -1im/6, 1/6, 0, 0, 0, 0),
(-1,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, 1/(3*sqrt(2)), 1im/(3*sqrt(2)), 0),
(-1,1) => SMatrix{3, 3, ComplexF64, 9}(-1/6, -1im/6, 0, -1im/6, 1/6, 0, 0, 0, 0),
(0,-1) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 1/(3*sqrt(2)), 0, 0, -1im/(3*sqrt(2)), 0, 0, 0),
(0,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, 0, 0, 1/3),
(0,1) => SMatrix{3, 3, ComplexF64, 9}(0, 0, -1/(3*sqrt(2)), 0, 0, -1im/(3*sqrt(2)), 0, 0, 0),
(1,-1) => SMatrix{3, 3, ComplexF64, 9}(-1/6, 1im/6, 0, 1im/6, 1/6, 0, 0, 0, 0),
(1,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, -1/(3*sqrt(2)), 1im/(3*sqrt(2)), 0),
(1,1) => SMatrix{3, 3, ComplexF64, 9}(1/6, -1im/6, 0, 1im/6, 1/6, 0, 0, 0, 0)
)

function erot_dot(j1,m1,mu1, j2,m2,mu2)
	"""
	computes < E_{m1 mu1} e_j1 , E_{m2 mu2} e_j2  >
	"""
   E1 = rmatrices(m1,mu1)
   E2 = rmatrices(m2,mu2)
   return dot(E1[:,j1],E2[:,j2])
end

function (A::ERot3DCoeffs{T})(ll::StaticVector{1},
                            mm::StaticVector{1},
                            kk::StaticVector{1}) where {T}
   if ll[1] == 1 && abs(mm[1]) <= 1 && abs(kk[1]) <= 1
      return  rmatrices[(mm[1],kk[1])]
   else
      return get0val(A)
   end
end

#function (A::IRot3DCoeffs{T})(ll::StaticVector{2},
#                            mm::StaticVector{2},
#                            kk::StaticVector{2}) where {T}
#   if ll[1] != ll[2] || sum(mm) != 0 || sum(kk) != 0
#      return T(0)
#   else
#      return T( 8 * pi^2 / (2*ll[1]+1) * (-1)^(mm[1]-kk[1]) )
#   end
#end

# next comes the recursion step for N ≧ 3

get0val(A::IRot3DCoeffs{T}) where T  = T(0)
get0val(A::ERot3DCoeffs{T}) where T = SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, 0, 0, 0)

function _compute_val(A::Rot3DCoeffs{T}, ll::StaticVector{N},
                                        mm::StaticVector{N},
                                        kk::StaticVector{N}) where {T, N}

   val = get0val(A)
   llp = ll[1:N-2]
   mmp = mm[1:N-2]
   kkp = kk[1:N-2]
   for j = abs(ll[N-1]-ll[N]):(ll[N-1]+ll[N])
      if abs(kk[N-1]+kk[N]) > j || abs(mm[N-1]+mm[N]) > j
         continue
      end
		cgk = try
			A.cg(ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
		catch
			@show (ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
			get0val(A)
		end
		cgm = A.cg(ll[N-1], mm[N-1], ll[N], mm[N], j, mm[N-1]+mm[N])
		if cgk * cgm  != 0
			val += cgk * cgm * A( SVector(llp..., j),
								       SVector(mmp..., mm[N-1]+mm[N]),
								       SVector(kkp..., kk[N-1]+kk[N]) )
		end
   end
   return val
end

# ----------------------------------------------------------------------
#   construction of a possible set of generalised CG coefficient;
#   numerically via SVD
# ----------------------------------------------------------------------

function ri_basis(A::IRot3DCoeffs{T}, ll::SVector; ordered=false) where {T}
	CC = compute_Al(A, ll, Val(ordered))
	svdC = svd(CC)
	rk = rank(Diagonal(svdC.S))
	return svdC.U[:, 1:rk]'
end

# unordered
function compute_Al(A::IRot3DCoeffs{T}, ll::SVector, ::Val{false}) where {T}
	len = length(_mrange(ll))
   CC = zeros(T, len, len)
   for (im, mm) in enumerate(_mrange(ll)), (ik, kk) in enumerate(_mrange(ll))
      CC[ik, im] = A(ll, mm, kk)
   end
   return CC
end

function re_basis(A::ERot3DCoeffs{T}, ll::SVector; ordered=false) where {T}
	GG = compute_gl(A, ll, Val(ordered))
	svdG = svd(GG)
	rk = rank(Diagonal(svdC.S))
	return svdC.U[:, 1:rk]'
end

function compute_gl(A::ERot3DCoeffs{T}, ll::SVector, ::Val{false}) where {T}
	len = length(_mrange1(ll))
   GG = zeros(T, len, len) # = Gramian of rotational eqiuvariant span
   # kk = mprime
   for (im, (jm,mm)) in enumerate(Iterators.product(1:3, _mrange1(ll))), (ik, (jk,kk)) in enumerate(Iterators.product(1:3, _mrange1(ll)))
	   for mu in _mrange1(ll)
	  		GG[im,ik] += dot(A(ll, mm, mu)[:,im],A(ll, kk, mu)[:,ik])
		end
   end
   return GG
end

# # ordered; TODO: check this out, clean it up and test it!!!
# function compute_Al(A::IRot3DCoeffs{T}, ll::SVector, ::Val{true}) where {T}
# 	num_mm_sorted = sum(mm -> issorted(mm), _mrange(ll))
# 	# @show num_mm_sorted
# 	num_mm = length(_mrange(ll))
#    CC = zeros(T, num_mm, num_mm_sorted)
# 	im = 0
#    for mm in _mrange(ll)
# 		if issorted(mm) # -> make this sorted relative to ll!!!
# 			im += 1
# 			for (ik, kk) in enumerate(_mrange(ll))
# 		      CC[ik, im] = A(ll, mm, kk)
# 			end
# 		end
# 	end
#    return CC
# end
#
#
# # two utility functions which are probably never used!
#
# compute_Al(ll::SVector{N}; ordered = false) where {N} =
# 		compute_Al(IRot3DCoeffs(N, sum(ll)), ll; ordered=ordered)
#
# compute_Al(A::IRot3DCoeffs, ll::SVector{N}; ordered = false) where {N} =
# 		compute_Al(A, ll, Val(ordered))

# TODO: this could use some documentation

rpi_basis(A::IRot3DCoeffs, zz, nn, ll) =
			rpi_basis(A, SVector(zz...), SVector(nn...), SVector(ll...))

function rpi_basis(A::IRot3DCoeffs,
						 zz::SVector{N},
						 nn::SVector{N, Int},
						 ll::SVector{N, Int}) where {N}
	Uri = ri_basis(A, ll)
	Mri = collect( _mrange(ll) )   # rows...
	G = _gramian(zz, nn, ll, Uri, Mri)
   S = svd(G)
   rk = rank(G; rtol =  1e-7)
	Urpi = S.U[:, 1:rk]'
	return Diagonal(sqrt.(S.S[1:rk])) * Urpi * Uri, Mri
end

function _gramian(zz, nn, ll, Uri, Mri)
   N = length(nn)
   nri = size(Uri, 1)
   @assert size(Uri, 1) == nri
   G = zeros(nri, nri)
   for σ in permutations(1:N)
      if (zz[σ] != zz) || (nn[σ] != nn) || (ll[σ] != ll); continue; end
      for (iU1, mm1) in enumerate(Mri), (iU2, mm2) in enumerate(Mri)
         if mm1[σ] == mm2
            for i1 = 1:nri, i2 = 1:nri
               G[i1, i2] += conj(Uri[i1, iU1]) * Uri[i2, iU2]
            end
         end
      end
   end
   return G
end

end


text:
using Dates
function save_results(sp)
    timestamp = string(now())
    #filename = "optimization_results/results"*timestamp
    filename = "results"*timestamp

    s = scenarios(sp)
    n = length(s)
    n_recourse = length(optimal_recourse_decision(sp,1))
    recourse_results = zeros(n, n_recourse+2)
    for i in 1:n
        recourse_results[i,1] = s[i].data[:t_xi]
        recourse_results[i,2] = s[i].data[:s_xi]*s[i].data[:F_xi]
        recourse_results[i,3:end] = optimal_recourse_decision(sp,i)[:]
    end
    #=scens = zeros(n,2)
    scens[:,1] = [s[i].data[:t_xi] for i in 1:n]
    scens[:,2] = [s[i].data[:s_xi]*s[i].data[:F_xi] for i in 1:n]
    #scens[:,3] = [s[i].data[:t_xi]]
    df = DataFrame(scens, :auto)
    CSV.write(filename*"scenarios.csv")
    for i in 1:n
        ord = optimal_recourse_decision(sp, i)
        df = DataFrame(reshape(ord, (length(ord), 1)), :auto)
        CSV.write(filename*"scenario$i.csv", df |> DataFrame)
    end=#
    od = optimal_decision(sp)
    df = DataFrame(reshape(od, (length(od), 1)), :auto)
    CSV.write(filename*".csv", df |> DataFrame)
    df_recourse = DataFrame(reshape(recourse_results, (n, n_recourse+2)),:auto)
    CSV.write(filename*"recourse_results.csv", df_recourse |> DataFrame)
end

ext: jl


text:
using Random, Optim

Base.@kwdef struct Config
    fuel_dist::Float64 = 2.e6
    seed::Int = 1234
    rng::MersenneTwister = MersenneTwister(seed)
    optimizer::Optim.AbstractOptimizer = NelderMead()
    models::Int = 20
end

struct BookKeep
    mask::BitArray # keeps visited locations
    betas::Vector{Float64} # optimization weights
    pointwise::Vector{Float64}
    function BookKeep(dim, feature_count)
        new(trues(dim), Vector{Float64}(undef, feature_count), Vector{Float64}(undef, dim))
    end
end

# finds most desirable location (local maximum)
function masked_weighted_argmax(mask, scores)
    id = 1
    m = floatmax(Float64) * -1
    @inbounds for i ∈ eachindex(mask)
        if mask[i] && scores[i] > m
            id = i
            m = scores[i]
        end
    end
    return id
end

function find_betas!(coord, locinfo, config, bookkeep, feat, betas; output=false)
    beer_count, traveled = 0, 0.
    bookkeep.mask .= true

    calc_scores!(feat, betas)
    id = masked_weighted_argmax(bookkeep.mask, feat.scores)
    beer_count += size(locinfo.beers[id], 1)
    traveled += bookkeep.pointwise[id]
    bookkeep.mask[id] = false
    output && println("$(Int(traveled÷1000))km $(locinfo.locations[id, :name]) +$(size(locinfo.beers[id], 1))")

    while true
        calc_scores!(feat, betas, id)
        id_new = masked_weighted_argmax(bookkeep.mask, feat.scores)
        bookkeep.pointwise[id_new] + coord.pairwise[id, id_new] + traveled > config.fuel_dist && break
        beer_count += size(locinfo.beers[id], 1)
        traveled += coord.pairwise[id, id_new]
        id = id_new
        bookkeep.mask[id] = false
        output && println("$(Int(traveled÷1000))km $(locinfo.locations[id, :name]) +$(size(locinfo.beers[id], 1))")
    end

    return -beer_count # only beer count is optimized since task is ambiguous
end

function run_optimization(coord, locinfo, config, bookkeep, feat, output=true)
    minimum = 1.
    minimizer = []
    for _ ∈ 1:config.models # few models are run with different betas initializations
        rand!(config.rng, bookkeep.betas)
        res = optimize(betas -> find_betas!(coord, locinfo, config, bookkeep, feat, betas), bookkeep.betas, config.optimizer)
        if Optim.minimum(res) < minimum
            minimum = Optim.minimum(res)
            minimizer = Optim.minimizer(res)
        end
    end

    find_betas!(coord, locinfo, config, bookkeep, feat, minimizer; output=output)
    output && println("Beer count: $(Int(-minimum))")
    return (Int(-minimum))
end


text:
type Serial{A<:AbstractArray, B<:AbstractArray} <: Container{A,B}
    output :: B
    gradInput :: A
    modules :: Array{Layer, 1}
end

Serial(T::Type, input_size, output_size) =
    Serial(Array(T, output_size),
           Array(T, input_size),
           Array(Layer, 0))

function add_layer!{A, B}(self :: Serial{A, B}, layer :: Layer)
    if length(self.modules) == 0
        @assert(shape(self.gradInput) == shape(layer.gradInput),
                "Layer shape mismatch")
    else
        @assert(shape(last(self.modules).output) == shape(layer.gradInput),
                "Layer shape mismatch")
    end
    push!(self.modules, layer)
end

function validate{A, B}(self :: Serial{A, B})
    local input_shape = shape(self.gradInput)

    for layer in self.modules
        @assert(shape(layer.gradInput) == input_shape,
                "Layer shape mismatch")
        input_shape = shape(layer.output)
    end

    @assert(shape(self.output) == input_shape, "Layer shape mismatch")
end

function updateOutput!{A, B}(self :: Serial{A, B}, input :: A)
    for i = 1:length(self.modules)
        input = updateOutput!(self.modules[i], input)
    end
    self.output = input
    return self.output
end

function updateGradInput!{A, B}(self :: Serial{A, B},
                                input :: A,
                                gradOutput :: B)
    if !isempty(self.modules)
        for i = length(self.modules):-1:2
            gradOutput = updateGradInput!(self.modules[i],
                                          self.modules[i-1].output,
                                          gradOutput)
        end

        gradOutput = updateGradInput!(self.modules[1], input, gradOutput)
    end

    self.gradInput = gradOutput
    return self.gradInput
end

function accGradParameters!{A, B}(self :: Serial{A, B},
                                  input :: A,
                                  gradOutput :: B,
                                  scale :: eltype(A))
    if isempty(self.modules)
        return
    end

    for i = length(self.modules):-1:2
        layer = self.modules[i]
        accGradParameters!(layer,
                           self.modules[i-1].output,
                           gradOutput,
                           scale)
        gradOutput = layer.gradInput
    end

    accGradParameters!(self.modules[1], input, gradOutput, scale)
end

export Serial


text:
#=
HDP.jl

Adham Beyki, odinay@gmail.com
27/10/2015

TODO:
Check  CRF sampler
=#

###################################################
## Hierarchical Dirichlet Process Mixture Models ##
###################################################
type HDP{T} <: hMixtureModel
    component::T
    KK::Int
    gg::Float64
    g1::Float64
    g2::Float64
    aa::Float64
    a1::Float64
    a2::Float64

    HDP{T}(c::T,
        KK::Int,
        gg::Float64, g1::Float64, g2::Float64,
        aa::Float64, a1::Float64, a2::Float64) = new(c, KK, gg, g1, g2, aa, a1, a2)
end

HDP{T}(c::T, KK::Int, gg::Real, g1::Real, g2::Real, aa::Real, a1::Real, a2::Real) = HDP{typeof(c)}(c, KK,
    convert(Float64, gg), convert(Float64, g1), convert(Float64, g2),
    convert(Float64, aa), convert(Float64, a1), convert(Float64, a2))

function Base.show(io::IO, hdp::HDP)
  println(io, "Hierarchical Dirichlet Process Mixture Model with $(hdp.KK) $(typeof(hdp.component)) components")
end

type CRFSample
    tji::Vector{Vector{Int}}
    njt::Vector{Vector{Int}}
    kjt::Vector{Vector{Int}}
    zz::Vector{Vector{Int}}

    CRFSample(tji::Vector{Vector{Int}}, njt::Vector{Vector{Int}}, kjt::Vector{Vector{Int}}, zz::Vector{Vector{Int}}) = new(tji, njt, kjt, zz)
end

count_active_clusters(nn::Matrix{Int}) = length(find(x -> x>0, sum(nn, 1)))

function storesample{T}(
    hdp::HDP{T},
    KK_list::Vector{Int},
    KK_dict::Dict{Int, Vector{Vector{Int}}},
    sample_n::Int,
    filename::ASCIIString)

    if endswith(filename, "_")
        dummy_filename = string(filename, sample_n, ".h5")
    else
        dummy_filename = string(filename, "_", sample_n, ".h5")
    end

    KK_dict_keys = collect(keys(KK_dict))
    n_groups = length(KK_dict[KK_dict_keys[1]])
    n_group_j = Array(Int, n_groups)
    for jj = 1:n_groups
        n_group_j[jj] = length(KK_dict[KK_dict_keys[1]][jj])
    end

    println("storing on disk...")
    HDF5.h5open(dummy_filename, "w") do file
        for kk in KK_dict_keys
            zz_flat = Int[]
            for jj=1:n_groups
                append!(zz_flat, KK_dict[kk][jj])
            end
            HDF5.write(file, "zz/$(kk)", zz_flat)
        end
        HDF5.write(file, "sample_n", sample_n)
        HDF5.write(file, "n_groups", n_groups)
        HDF5.write(file, "n_group_j", n_group_j)
    end
end

function sample_hyperparam!(hdp::HDP, n_group_j::Vector{Int}, m::Int)
    #= NOTE
    resampling the group level concentration parameter α0 using auxiliary variables
    w and s, Eq. 50 Teh etal 04 the Gamma distribution in Eq.50 is expressed using
    shape and rate. We have to rescale them to shape and scale representation to be
    able to use standard random gamma functions in Julia to draw from it. Also:
    Gamma(a, 1/b) = Gamma(a) / b
    =#

    n_groups = length(n_group_j)
    w = zeros(Float64, n_groups)
    for jj = 1:n_groups
        w[jj] = rand(Distributions.Beta(hdp.aa+1, n_group_j[jj]))
    end

    pp = n_group_j / hdp.aa
    s = zeros(Int, n_groups)
    for jj = 1:n_groups
        s[jj] = Int(rand() <= pp[jj] / (pp[jj] + 1.0))
    end

    aa_shape = hdp.a1 + m - sum(s)
    aa_rate  = hdp.a2 - sum(log(w))
    hdp.aa = rand(Distributions.Gamma(aa_shape)) / aa_rate

    # resampling the top level concentration parameter γ, Escobar and West 95
    eta_pp = rand(Distributions.Beta(hdp.gg+1, m))
    rr = (hdp.g1 + hdp.KK - 1) / (m*(hdp.g2 - log(eta_pp)))
    pi_eta = rr / (1.0 + rr)

    if rand() < pi_eta
        hdp.gg = rand(Distributions.Gamma(hdp.g1 + hdp.KK)) / (hdp.g2-log(eta_pp))
    else
        hdp.gg = rand(Distributions.Gamma(hdp.g1+hdp.KK-1)) / (hdp.g2-log(eta_pp))
    end
end # sample_hyperparam
function sample_hyperparam!(hdp::HDP, n_group_j::Vector{Int}, m::Int, n_internals::Int)
    for iteration = 1:n_internals
        sample_hyperparam!(hdp, n_group_j, m)
    end
end # sample_hyperparam

function collapsed_gibbs_sampler{T1, T2}(
    hdp::HDP{T1},
    xx::Vector{Vector{T2}},
    zz::Vector{Vector{Int}},
    n_burnins::Int, n_lags::Int, n_samples::Int,
    sample_hyperparam::Bool=true, n_internals::Int=10,
    store_every::Int=100, results_path="", filename::ASCIIString="HDP_results_")

    n_iterations    = n_burnins + (n_samples)*(n_lags+1)
    n_groups        = length(xx)
    n_group_j       = zeros(Int, n_groups)
    nn              = zeros(Int, n_groups, hdp.KK)
    log_likelihood  = 0.0

    for jj = 1:n_groups
        n_group_j[jj] = length(zz[jj])
    end

    # constructing components
    components = Array(typeof(hdp.component), hdp.KK)
    for kk = 1:hdp.KK
        components[kk] = deepcopy(hdp.component)
    end

    KK_list = Int[]
    KK_dict = Dict{Int, Vector{Vector{Int}}}()

    snumbers_file = string(Pkg.dir(), "\\BNP\\src\\StirlingNums_10K.mat")
    snumbers_data = MAT.matread(snumbers_file)
    snumbers = snumbers_data["snumbersNormalizedSparse"]

    #############################################################
    #                  Initializing the model                   #
    #############################################################
    tic()
    print_with_color(:red, "\nInitializing the model\n")

    for jj = 1:n_groups
        for ii = 1:n_group_j[jj]
            kk = zz[jj][ii]
            additem!(components[kk], xx[jj][ii])
            nn[jj, kk] += 1
            log_likelihood += loglikelihood(components[kk], xx[jj][ii])
        end
    end
    push!(KK_list, hdp.KK)

    my_beta = ones(Float64, hdp.KK+1) / (hdp.KK+1)

    elapsed_time = toq()

    ##############################################################
    #                   starting the MCMC chain                  #
    ##############################################################
    for iteration = 1:n_iterations

        if iteration < n_burnins
            print_with_color(:blue, "Burning... ")
        end
        println(@sprintf("iteration: %d, KK=%d, KK mode=%d, gg=%.2f, aa=%.2f, time=%.2f, likelihood=%.2f",
            iteration, hdp.KK, indmax(hist(KK_list, .5:maximum(KK_list)+0.5)[2]),
            hdp.gg, hdp.aa, elapsed_time, log_likelihood))
        log_likelihood = 0.0

        tic()

        # make sure all the components are active
        for kk = 1:hdp.KK
            if sum(nn[:, kk]) == 0
                println("\tcomponent $kk has become inactive")
                nn = del_column(nn, kk)
                splice!(components, kk)

                for ll in 1:n_groups
                    idx = find(x -> x>kk, zz[ll])
                    zz[ll][idx] -= 1
                end

                my_beta[hdp.KK+1] += my_beta[kk]
                splice!(my_beta, kk)
                hdp.KK -= 1
            end
        end

        #=
        Now
        1. resample zz[jj][ii]
            for all jj, ii
                1. remove xx[jj][ii] from its cluster
                2. resample zz[jj][ii]
                3. add x[jj][ii] to the resampled cluster
            end
        2. resample β
        3. resample hyperparams γ, α0
        =#

        # sample zz[jj][ii]
        for jj = randperm(n_groups)
            for ii = randperm(n_group_j[jj])

                ##   1.1   ##
                #  removing the observation
                kk = zz[jj][ii]
                delitem!(components[kk], xx[jj][ii])
                nn[jj, kk] -= 1

                # If as the result of removing xx[jj][ii]
                # phi_k becomes inactive, remove the atom
                if sum(nn[:, kk]) == 0
                    println("\tcomponent $kk has become inactive")
                    nn = del_column(nn, kk)
                    splice!(components, kk)

                    # shifting zz > kk back
                    for ll in 1:n_groups
                        idx = find(x -> x>kk, zz[ll])
                        zz[ll][idx] -= 1
                    end

                    # removing the stick kk
                    my_beta[hdp.KK+1] += my_beta[kk]
                    splice!(my_beta, kk)
                    hdp.KK -= 1
                end

                ##   1.2   ##
                # resample kk
                pp = zeros(Float64, hdp.KK+1)
                for ll = 1:hdp.KK
                    pp[ll] = log(nn[jj, ll] + hdp.aa * my_beta[ll]) + logpredictive(components[ll], xx[jj][ii])
                end
                pp[hdp.KK+1] = log(hdp.aa * my_beta[hdp.KK+1]) + logpredictive(hdp.component, xx[jj][ii])

                lognormalize!(pp)
                kk = sample(pp)

                # if kk is new, we instantiate a new cluster
                if kk == hdp.KK+1
                    println("\tcomponent $(kk) activated.")
                    push!(components, deepcopy(hdp.component))
                    nn = add_column(nn)

                    b = rand(Distributions.Beta(1, hdp.gg))
                    b_new = my_beta[hdp.KK+1]
                    my_beta[hdp.KK+1] = b * b_new
                    push!(my_beta, (1-b)*b_new)
                    hdp.KK += 1
                end

                ##   1.3   ##
                # adding the data point to the resampled cluster
                zz[jj][ii] = kk
                additem!(components[kk], xx[jj][ii])
                nn[jj, kk] += 1
                log_likelihood += loglikelihood(components[kk], xx[jj][ii])

            end # n_group_j
        end # n_groups

        ##   2  ##
        # resampling β vector using auxiliary variable method
        # Eq. 40, 41 Teh etal 2004
        M = zeros(Int, n_groups, hdp.KK)
        for hh in 1:n_internals
            for jj = 1:n_groups
                for kk = 1:hdp.KK
                    if nn[jj, kk] == 0
                        M[jj, kk] = 0
                    else
                        rr = zeros(Float64, nn[jj, kk])
                        for mm = 1:nn[jj, kk]
                            rr[mm] = log(snumbers[nn[jj, kk], mm]) + mm*log(hdp.aa * my_beta[kk])
                        end
                        lognormalize!(rr)
                        M[jj, kk] = sample(rr)
                    end
                end # kk
            end # n_groups
            MM = convert(Vector{Float64}, vec(sum(M, 1)))
            push!(MM, hdp.gg)
            my_beta = rand(Distributions.Dirichlet(MM))

            ##   3   ##
            # resampling hyper-parameters γ and α0
            if sample_hyperparam
                m = sum(M)
                sample_hyperparam!(hdp, n_group_j, m)
            end
        end # n_internals
        elapsed_time = toq()

        # save the sample
        if (iteration-n_burnins) % (n_lags+1) == 0 &&  iteration > n_burnins
            sample_n = convert(Int, (iteration-n_burnins)/(n_lags+1))

            push!(KK_list, hdp.KK)
            KK_dict[hdp.KK] = deepcopy(zz)

            if sample_n % store_every == 0
                normalized_filename = normpath(joinpath(results_path, filename))
                storesample(hdp, KK_list, KK_dict, sample_n, normalized_filename)
            end
        end
    end # iteration

    sample_n = convert(Int, (n_iterations-n_burnins)/(n_lags+1))
    normalized_filename = normpath(joinpath(results_path, filename))
    storesample(hdp, KK_list, KK_dict, sample_n, normalized_filename)

    KK_dict_flat = Dict{Int, Vector{Int}}()
    KK_dict_keys = collect(keys(KK_dict))
    for kk in KK_dict_keys
        zz_flat = Int[]
        for jj=1:n_groups
            append!(zz_flat, KK_dict[kk][jj])
        end
        KK_dict_flat[kk] = zz_flat
    end

    KK_list, KK_dict_flat

end # collapsed_gibbs_sampler!

function CRF_gibbs_sampler{T1, T2}(
    hdp::HDP{T1},
    xx::Vector{Vector{T2}},
    zz::Vector{Vector{Int}},
    n_burnins::Int, n_lags::Int, n_samples::Int,
    sample_hyperparam::Bool=true, n_internals::Int=10,
    store_every::Int=100, filename::ASCIIString="HDP_results_",
    KK_list::Vector{Int}=Int[],
    KK_dict::Dict{Int, Vector{Vector{Int}}}=Dict{Int, Vector{Vector{Int}}}())

    # constructing the components
    components = Array(typeof(hdp.component), hdp.KK)
    for kk = 1:hdp.KK
        components[kk] = deepcopy(hdp.component)
    end

    n_iterations = n_burnins + (n_samples)*(n_lags+1)
    n_groups     = length(xx)
    n_group_j    = Int[length(zz[jj]) for jj = 1:n_groups]
    tji          = Array(Vector{Int}, n_groups)     # tji[jj][ii] represents the table that customer ii in restaurant jj sits at
    njt          = Array(Vector{Int}, n_groups)     # njt[jj][tbl] represents the number of customers sitting at table tbl in restaurant jj
    kjt          = Array(Vector{Int}, n_groups)     # kjt[jj][tbl] represents the dish being served at table tbl at restaurant jj
    mm           = zeros(Int, hdp.KK)                # mm[kk] represents the number of times dish kk is ordered

    if length(KK_list) == 0
        n_sampels_old = 0
        KK_list = zeros(Int, n_samples)
        # KK_dict = Dict{Int, CRFSample}()
        KK_dict = Dict{Int, Vector{Vector{Int}}}()
    else
        n_sampels_old = length(KK_list)
        KK_list = vcat(KK_list, zeros(Int, n_samples))
    end

    ##############################
    #   Initializing the model   #
    ##############################
    tic()

    print_with_color(:red, "\nInitializing the model\n")
    log_likelihood = 0.0

    for jj = 1:n_groups

        kjt[jj] = unique(zz[jj])        # table settings in group jj
        n_tbls  = length(kjt[jj])       # number of tables in group jj

        tji[jj] = zeros(Int, n_group_j[jj])
        njt[jj] = zeros(Int, n_tbls)

        for ii = 1:n_group_j[jj]
            kk = zz[jj][ii]
            additem!(components[kk], xx[jj][ii])
            log_likelihood += loglikelihood(components[kk], xx[jj][ii])

            tbl = findfirst(kjt[jj], kk)
            tji[jj][ii] = tbl
            njt[jj][tbl] += 1
        end
    end

    for jj = 1:n_groups
        for kk in kjt[jj]
            mm[kk] += 1
        end
    end
    KK_list[1] = hdp.KK
    elapsed_time = toq()

    ###################################
    #     starting the MCMC chain     #
    ###################################
    for iteration = 1:n_iterations

        if iteration < n_burnins
            print_with_color(:blue, "Burning... ")
        end

        println(@sprintf("iteration: %d, KK=%d, KK mode=%d, gg=%.2f, aa=%.2f, time=%.2f, likelihood=%.2f",
            iteration, hdp.KK, indmax(hist(KK_list, .5:maximum(KK_list)+0.5)[2]),
            hdp.gg, hdp.aa, elapsed_time, log_likelihood))

        log_likelihood = 0.0

        tic()

        ######################################
        # 1. resampling tables for customers #
        ######################################
        for jj = 1:n_groups

            n_tbls = length(kjt[jj])

            for ii = 1:n_group_j[jj]

                # remove the effect of the data point
                kk = zz[jj][ii]
                tbl = tji[jj][ii]
                delitem!(components[kk], xx[jj][ii])
                njt[jj][tbl] -= 1

                # if by removing xx[jj][ii] table tji[jj][ii] becomes empty
                # remove the table as well
                if njt[jj][tbl] == 0

                    n_tbls -= 1
                    mm[kk] -= 1

                    splice!(kjt[jj], tbl)
                    splice!(njt[jj], tbl)

                    idx = find(x -> x>tbl, tji[jj])
                    tji[jj][idx] -= 1

                    # if tji[jj][ii] (which we have removed it) is the only table which is
                    # serving dish kk and by removing it no other table serves dish kk,
                    # remove dish kk from menu
                    if mm[kk] == 0

                        hdp.KK -= 1
                        splice!(mm, kk)
                        splice!(components, kk)

                        for ll = 1:n_groups
                            idx = find(x -> x>kk, zz[ll])
                            zz[ll][idx] -= 1

                            idx = find(x -> x>kk, kjt[ll])
                            kjt[ll][idx] -= 1
                        end
                    end # if mm[kk]
                end # if njt

                # resample table tji[jj][ii]
                pp = zeros(Float64, n_tbls+1)
                for tbl = 1:n_tbls
                    kk = kjt[jj][tbl]
                    pp[tbl] = log(njt[jj][tbl]) + logpredictive(components[kk], xx[jj][ii])
                end
                pp[n_tbls+1] = log(hdp.aa) + logpredictive(hdp.component, xx[jj][ii])

                lognormalize!(pp)
                tbl = sample(pp)

                # if a new table is selected, draw the dish for the new table
                if tbl == n_tbls+1

                    n_tbls += 1
                    push!(kjt[jj], 0)
                    push!(njt[jj], 0)

                    pp = zeros(Float64, hdp.KK+1)
                    for kk = 1:hdp.KK
                        pp[kk] = log(mm[kk]) + logpredictive(components[kk], xx[jj][ii])
                    end
                    pp[hdp.KK+1] = log(hdp.gg) + logpredictive(hdp.component, xx[jj][ii])

                    lognormalize!(pp)
                    kk = sample(pp)

                    if kk == hdp.KK+1
                        hdp.KK += 1
                        push!(components, deepcopy(hdp.component))
                        push!(mm, 0)
                    end

                    kjt[jj][tbl] = kk
                    mm[kk] += 1
                end

                kk = kjt[jj][tbl]
                tji[jj][ii] = tbl
                njt[jj][tbl] += 1
                additem!(components[kk], xx[jj][ii])
                zz[jj][ii] = kk
            end # n_group_j[jj]
        end # n_groups

        ########################
        # 2. resampling dishes #
        ########################
        for jj = 1:n_groups

            # number of tables in restaurant jj
            n_tbls = length(kjt[jj])

            for tbl = 1:n_tbls

                kk = kjt[jj][tbl]
                mm[kk] -= 1

                tidx = find(x -> x==tbl, tji[jj])

                # if tbl is the only table serving dish kk, by removing it we need to
                # remove the dish from the menu

                if mm[kk] == 0
                    hdp.KK -= 1
                    splice!(mm, kk)
                    splice!(components, kk)

                    for ll = 1:n_groups
                        idx = find(x -> x>kk, zz[ll])
                        zz[ll][idx] -= 1

                        idx = find(x -> x>kk, kjt[ll])
                        kjt[ll][idx] -= 1
                   end

               else
                    # otherwise we need to remove the customers who were sitting at table tbl in
                    # jj, having dish kk, one by one
                    for ll in tidx
                       delitem!(components[kk], xx[jj][ll])
                    end
               end

                # resampling the dish
                pp = zeros(Float64, hdp.KK+1)
                for kk = 1:hdp.KK
                    pp[kk] = log(mm[kk])
                    for ll in tidx
                        pp[kk] += logpredictive(components[kk], xx[jj][ll])
                    end
                end
                pp[hdp.KK+1] += log(hdp.gg)
                for ll in tidx
                    pp[hdp.KK+1] += logpredictive(hdp.component, xx[jj][ll])
                end

                lognormalize!(pp)
                kk = sample(pp)

                if kk == hdp.KK+1
                    hdp.KK += 1
                    push!(components, deepcopy(hdp.component))
                    push!(mm, 0)
                end

                kjt[jj][tbl] = kk
                for ll in tidx
                    additem!(components[kk], xx[jj][ll])
                end
                mm[kk] += 1
                zz[jj][tidx] = kk
            end # tbl
        end # n_groups

        # resampling hyper-parameters
        if sample_hyperparam
            nn = zeros(Int, n_groups, hdp.KK)
            for jj = 1:n_groups
                for ii = 1:n_group_j[jj]
                    nn[jj, zz[jj][ii]] += 1
                end
            end
            m = sum([length(kjt[jj]) for jj=1:n_groups])
            sample_hyperparam!(hdp, n_group_j, m, n_internals)
        end

        for jj = 1:n_groups
            for ii = 1:n_group_j[jj]
            log_likelihood += logpredictive(components[zz[jj][ii]], xx[jj][ii])
            end
        end

        elapsed_time = toq()

        # save the sample
        if (iteration-n_burnins) % (n_lags+1) == 0 &&  iteration > n_burnins
            sample_n = n_sampels_old + convert(Int, (iteration-n_burnins)/(n_lags+1))
            KK_list[sample_n] = hdp.KK
            # crf_sample = CRFSample(deepcopy(tji), deepcopy(njt), deepcopy(kjt), deepcopy(zz))
            # KK_dict[hdp.KK] = crf_sample
            KK_dict[hdp.KK] = deepcopy(zz)
            if sample_n % store_every == 0
                storesample(hdp, KK_list, KK_dict, n_burnins, n_lags, sample_n, filename)
            end
        end

    end # iteration

    sample_n = n_sampels_old + convert(Int, (n_iterations-n_burnins)/(n_lags+1))
    KK_list[sample_n] = hdp.KK
    # crf_sample = CRFSample(deepcopy(tji), deepcopy(njt), deepcopy(kjt), deepcopy(zz))
    # KK_dict[hdp.KK] = crf_sample
    KK_dict[hdp.KK] = deepcopy(zz)
    storesample(hdp, KK_list, KK_dict, n_burnins, n_lags, sample_n, filename)

    KK_list, KK_dict
end # CRF_gibbs_sampler

function truncated_gibbs_sampler{T1, T2}(
    hdp::HDP{T1},
    KK_truncation::Int,
    xx::Vector{Vector{T2}},
    zz::Vector{Vector{Int}},
    n_burnins::Int, n_lags::Int, n_samples::Int,
    sample_hyperparam::Bool=true, n_internals::Int=10,
    store_every::Int=250, results_path="", filename::ASCIIString="HDP_results_")

    # used when sampling the beta variables for pi_j
    BETA_THRESHOLD = 0.0001

    n_iterations    = n_burnins + (n_samples)*(n_lags+1)
    n_groups        = length(xx)
    n_group_j       = zeros(Int, n_groups)
    nn              = zeros(Int, n_groups, KK_truncation)
    pp              = zeros(Float64, KK_truncation)
    log_likelihood  = 0.0

    for jj = 1:n_groups
        n_group_j[jj] = length(zz[jj])
    end

    components = Array(typeof(hdp.component), KK_truncation)
    for kk = 1:KK_truncation
        components[kk] = deepcopy(hdp.component)
    end

    KK_list = Int[]
    KK_dict = Dict{Int, Vector{Vector{Int}}}()

    beta_tilde = zeros(Float64, KK_truncation)
    pi_tilde   = zeros(Float64, n_groups, KK_truncation)

    snumbers_file = string(Pkg.dir(), "\\BNP\\src\\StirlingNums_10K.mat")
    snumbers_data = MAT.matread(snumbers_file)
    snumbers = snumbers_data["snumbersNormalizedSparse"]

    ###################################################################
    #                      Initializing the model                     #
    ###################################################################

    tic()
    print_with_color(:red, "\nInitializing the model\n")

    # 1
    # adding the data points
    for jj = 1:n_groups
        for ii = 1:n_group_j[jj]
            kk = zz[jj][ii]
            additem!(components[kk], xx[jj][ii])
            nn[jj, kk] += 1
            log_likelihood += loglikelihood(components[kk], xx[jj][ii])
        end
    end
    push!(KK_list, count_active_clusters(nn))

    # 2
    # constructing my_beta
    for kk = 1:KK_truncation-1
        gamma1 = 1 + sum(nn[:, kk])
        gamma2 = hdp.gg + sum(nn[:, kk+1:KK_truncation])
        beta_tilde[kk] = gamma1 / (gamma1 + gamma2)
    end
    beta_tilde[KK_truncation] = 1.0
    my_beta = stick_breaking(beta_tilde)

    # 3
    # constructing beta variables for pi_j
    for jj = 1:n_groups
        for kk = 1:KK_truncation-1
            gamma1 = hdp.aa * my_beta[kk] + nn[jj, kk]
            gamma2 = hdp.aa * (1 - sum(my_beta[1:kk])) + sum(nn[jj, kk+1:KK_truncation])
            pi_tilde[jj, kk] = gamma1 / (gamma1 + gamma2)
        end
        pi_tilde[jj, KK_truncation] = 1.0
    end

    elapsed_time = toq()

    ###################################################################
    #                    starting the MCMC chain                      #
    ###################################################################

    for iteration = 1:n_iterations

        # Verbose
        if iteration < n_burnins
            print_with_color(:blue, "Burning... ")
        end
        println(@sprintf("iteration: %d, KK=%d, KK mode=%d, gg=%.2f, aa=%.2f, time=%.2f, likelihood=%.2f",
            iteration, hdp.KK, indmax(hist(KK_list, .5:maximum(KK_list)+0.5)[2]), hdp.gg, hdp.aa,
            elapsed_time, log_likelihood))
        log_likelihood = 0.0

        tic()

        # A
        # resample zz
        for jj = 1:n_groups
            for ii = 1:n_group_j[jj]

                # 1
                # remove the data point
                kk = zz[jj][ii]
                delitem!(components[kk], xx[jj][ii])
                nn[jj, kk] -= 1

                # 2
                # resample for zz[jj][ii]
                pp = zeros(Float64, KK_truncation)
                for kk = 1:KK_truncation
                    pp[kk] = log(pi_tilde[jj, kk]) + sum(log(1-pi_tilde[jj, 1:kk-1])) + logpredictive(components[kk], xx[jj][ii])
                end
                lognormalize!(pp)
                kk = sample(pp)

                # 3
                # add the data point to the newly sampled cluster
                zz[jj][ii] = kk
                nn[jj, kk] += 1
                additem!(components[kk], xx[jj][ii])
                log_likelihood += loglikelihood(components[kk], xx[jj][ii])

            end # ii
        end # jj
        hdp.KK = count_active_clusters(nn)

        # B
        # resample beta and pi

        # resample my_beta
        for kk = 1:KK_truncation-1
            gamma1 = 1 + sum(nn[:, kk])
            gamma2 = hdp.gg + sum(nn[:, kk+1:KK_truncation])
            beta_tilde[kk] = rand(Distributions.Beta(gamma1, gamma2))
        end
        beta_tilde[KK_truncation] = 1.0
        my_beta = stick_breaking(beta_tilde)

        # resample pi_tilde
        for jj = 1:n_groups
            for kk = 1:KK_truncation-1
                gamma1 = hdp.aa * my_beta[kk] + nn[jj, kk]
                gamma2 = hdp.aa * (1 - sum(my_beta[1:kk])) + sum(nn[jj, kk+1:KK_truncation])

                if gamma1 < BETA_THRESHOLD || gamma2 < BETA_THRESHOLD
                    pi_tilde[jj, kk] = 0.0
                else
                    pi_tilde[jj, kk] = rand(Distributions.Beta(gamma1, gamma2))
                end

            end
            pi_tilde[jj, KK_truncation] = 1.0
        end

        # C
        # resample hyper parameters
        if sample_hyperparam
            M = zeros(Int, n_groups, KK_truncation)
            for hh in 1:n_internals
                for jj = 1:n_groups
                    for kk = 1:KK_truncation
                        if sum(nn[:, kk]) != 0
                            if nn[jj, kk] == 0
                                M[jj, kk] = 0
                            else
                                rr = zeros(Float64, nn[jj, kk])
                                for mm = 1:nn[jj, kk]
                                    rr[mm] = log(snumbers[nn[jj, kk], mm]) + mm*log(hdp.aa * my_beta[kk])
                                end
                                lognormalize!(rr)
                                M[jj, kk] = sample(rr)
                            end
                        end
                    end # kk
                end # n_groups

                m = sum(M)
                sample_hyperparam!(hdp, n_group_j, m)
            end # n_internals
        end # sample_hyperparam

        elapsed_time = toq()

        # save the sample
        if (iteration-n_burnins) % (n_lags+1) == 0 &&  iteration > n_burnins
            sample_n = convert(Int, (iteration-n_burnins)/(n_lags+1))

            push!(KK_list, hdp.KK)
            KK_dict[hdp.KK] = deepcopy(zz)

            if sample_n % store_every == 0
                normalized_filename = normpath(joinpath(results_path, filename))
                storesample(hdp, KK_list, KK_dict, sample_n, normalized_filename)
            end
        end
    end # iteration

    sample_n = convert(Int, (n_iterations-n_burnins)/(n_lags+1))
    normalized_filename = normpath(joinpath(results_path, filename))
    storesample(hdp, KK_list, KK_dict, sample_n, normalized_filename)

    # before returning KK_dict, flatten zz to make it consistent
    # with the stored results
    KK_dict_flat = Dict{Int, Vector{Int}}()
    KK_dict_keys = collect(keys(KK_dict))
    for kk in KK_dict_keys
        zz_flat = Int[]
        for jj=1:n_groups
            append!(zz_flat, KK_dict[kk][jj])
        end
        KK_dict_flat[kk] = zz_flat
    end

    KK_list, KK_dict_flat
end

function posterior{T1, T2}(
    hdp::HDP{T1},
    xx::Vector{Vector{T2}},
    KK_dict::Dict{Int, Vector{Vector{Int}}},
    KK::Int)

    n_groups  = length(xx)
    n_group_j = [length(xx[jj]) for jj = 1:n_groups]

    components = Array(typeof(hdp.component), KK)
    for kk = 1:KK
        components[kk] = deepcopy(hdp.component)
    end

    nn = zeros(Int, n_groups, KK)
    zz = KK_dict[KK]

    for jj = 1:n_groups
        for ii =     1:n_group_j[jj]
            kk = zz[jj][ii]
            additem!(components[kk], xx[jj][ii])
            nn[jj, kk] += 1
        end
    end
    pij = nn + hdp.aa
    pij = pij ./ sum(pij, 1)

    pos_components = Array(typeof(posterior(hdp.component)), KK)
    for kk = 1:KK
        pos_components[kk] = posterior(components[kk])
    end

    pos_components, KK_dict[KK], nn, pij
end

function posterior{T1, T2}(
    hdp::HDP{T1},
    xx::Vector{Vector{T2}},
    KK_dict::Dict{Int, CRFSample},
    KK::Int,
    join_tables::Bool=false)

    n_groups  = length(xx)
    n_group_j = [length(xx[jj]) for jj = 1:n_groups]
    nn = zeros(Int, n_groups, KK)
    mm = zeros(Int, KK)

    components = Array(typeof(hdp.component), KK)
    for kk = 1:KK
        components[kk] = deepcopy(hdp.component)
    end

    tji = KK_dict[KK].tji
    njt = KK_dict[KK].njt
    kjt = KK_dict[KK].kjt
    zz  = KK_dict[KK].zz

    for jj = 1:n_groups
        for ii = 1:n_group_j[jj]
            kk = zz[jj][ii]
            additem!(components[kk], xx[jj][ii])
            nn[jj, kk] += 1
        end
    end

    pij = nn + hdp.aa
    pij = pij ./ sum(pij, 1)

    if join_tables
        for jj = 1:n_groups

            kk_list = unique(kjt[jj])

            if length(kk_list) != length(kjt[jj])
                new_kjt = zeros(Int, length(kk_list))
                new_njt = zeros(Int, length(kk_list))
                ll = 1

                for unq_kk in kk_list
                    idx = find(x -> x==unq_kk, kjt[jj])
                    new_kjt[ll] = unq_kk
                    new_njt[ll] = sum(njt[jj][idx])

                    for idid in idx
                        tbl_idx = find(x -> x==idid, tji[jj])
                        tji[jj][tbl_idx] = ll
                    end

                    ll += 1
                end

                kjt[jj] = new_kjt
                njt[jj] = new_njt
            end
        end
    end

    mm = zeros(Int, KK)
    for jj = 1:n_groups
        for kk in kjt[jj]
            mm[kk] += 1
        end
    end

    [posterior(components[kk]) for kk =1:KK], tji, njt, kjt, zz, nn, mm, pij

end

function posterior{T1, T2}(
    hdp::HDP{T1},
    xx::Vector{Vector{T2}},
    zz_flat::Vector{Int},
    KK::Int)

    n_groups = length(xx)
    nn = zeros(Int, n_groups, KK)
    n_group_j = Array(Int, n_groups)

    for jj = 1:n_groups
        n_group_j[jj] = length(xx[jj])
    end

    n_group_cs = cumsum(n_group_j)
    zz = Array(Vector{Int}, n_groups)
    zz[1] = zz_flat[1:n_group_cs[1]]
    for jj = 2:n_groups
        zz[jj] = zz_flat[n_group_cs[jj-1]+1:n_group_cs[jj]]
    end

    components = Array(typeof(hdp.component), KK)
    for kk = 1:KK
        components[kk] = deepcopy(hdp.component)
    end

    zz_unq = sort(unique(zz_flat))
    for jj = 1:n_groups
        for ii = 1:n_group_j[jj]
            kk = find(x -> x==zz[jj][ii], zz_unq)[1]
            additem!(components[kk], xx[jj][ii])
            nn[jj, kk] += 1
        end
    end
    pij = nn + hdp.aa
    pij = pij ./ sum(pij, 1)

    pos_components = Array(typeof(posterior(hdp.component)), KK)
    for kk = 1:KK
        pos_components[kk] = posterior(components[kk])
    end

    pos_components, nn, pij
end


text:
function minimize!(objective::ObjWrapper, x0, scheme, options=OptOptions())
    minimize!(objective, (x0, nothing), (scheme, NWI()), options)
end
function minimize!(objective::ObjWrapper, x0, approach::Tuple{<:Any, <:TRSPSolver}, options=OptOptions())
    minimize!(objective, (x0, nothing), approach, options)
end
function minimize!(objective::ObjWrapper, state0::Tuple, scheme, options=OptOptions())
    minimize!(objective, state0, (scheme, NWI()), options)
end
function minimize!(objective::ObjWrapper, state0::Tuple, approach::Tuple{<:Any, <:TRSPSolver}, options=OptOptions())
    tr_minimize!(objective, state0, approach, options)
end
include("optimize/inplace_loop.jl")
include("optimize/outofplace_loop.jl")


text:
@testset "1022.sum-of-root-to-leaf-binary-numbers.jl" begin
    @test sum_root_to_leaf(TreeNode{Int}([1,0,1,0,1,0,1])) == 22
    @test sum_root_to_leaf(TreeNode(1)) == 1
    @test sum_root_to_leaf(TreeNode(0)) == 0
    @test sum_root_to_leaf(TreeNode{Int}([1, 1])) == 3
end


text:
"""
builds the constructor which will organize cell features
eventually there will be spherocylindrical features
"""
struct Cell{T<:Real}
    p::Vector{T} #position of the central point of the cell
    v::Vector{T} #velocity
    d::Int64 # calculated death (aka division) of the cell
	o::Int64 # orientation, [0:180], of line segment of cell center +/- portion of length with respect to vertical axis of photo
    v2::Vector{T} # velocity half step
    a::Vector{T} # acceleration
	birth::Int64
	motherID::Int64
	sisterID::Int64
	daughterID::Vector{Int64}
end

Cell(p::Vector{T}, v::Vector{T}, d::Int64, o::Int64; birth=0, motherID=0, sisterID=0, daughterID=zeros(Int64,2)) where T<:Real = Cell(p,v,d,o,zeros(T,2),zeros(T,2),birth,motherID,sisterID,daughterID)

function makeCellArray(n, sparsity_level, μ)
	"""
	Draws starting positions & velocities from Normal Distr and initiates cells.
	"""
	d = Normal(0,sparsity_level)
	pos = rand(d, (2,n))
	vel = rand(Normal(), (2,n))
	pos .*= sparsity_level
	vel ./= 1000.

	t_division = rand(DiscreteUniform(100,μ),n)
	#randomkick = rand(DiscreteUniform(4,μ),n)
	orientation = rand(DiscreteUniform(0,180),n)

	cellarray = Cell[]
	for i in 1:n
		boop = [Cell(pos[:,i],vel[:,i],t_division[i],orientation[i])]
		append!(cellarray,boop)
	end
	return cellarray
end

max_stars_repo_name: andi-bergeson/simgym

max_issues_repo_name: andi-bergeson/simgym

max_forks_repo_name: andi-bergeson/simgym


text:
@testset "Bus Constructors" begin
    tBus = Bus(nothing)
    tLoadZones = LoadZones(nothing)

    bus = Bus(1, "test", PowerSystems.SLACK::BusType, 0.0, 0.0, (min=0.0, max=0.0), nothing)
    @test PowerSystems.get_bustype(bus) == PowerSystems.REF::BusType

    @test_throws(PowerSystems.DataFormatError,
                 Bus(1, "test", PowerSystems.ISOLATED::BusType, 0.0, 0.0,
                     (min=0.0, max=0.0), nothing))
end

@testset "Generation Constructors" begin
    tThreePartCost = ThreePartCost(nothing)
    @test tThreePartCost isa PowerSystemType
    tTwoPartCost = TwoPartCost(nothing)
    @test tTwoPartCost isa PowerSystemType
    tTechThermal = TechThermal(nothing)
    @test tTechThermal isa PowerSystemType
    tThermalGen = ThermalStandard(nothing)
    @test tThermalGen isa PowerSystems.Component
    tTechHydro = TechHydro(nothing)
    @test tTechHydro isa PowerSystemType
    tHydroFix = HydroFix(nothing)
    @test tHydroFix isa PowerSystems.Component
    tHydroDispatch = HydroDispatch(nothing)
    @test tHydroDispatch isa PowerSystems.Component
    tHydroStorage = HydroStorage(nothing)
    @test tHydroStorage isa PowerSystems.Component
    tTechRenewable = TechRenewable(nothing)
    @test tTechRenewable isa PowerSystemType
    tRenewableFix = RenewableFix(nothing)
    @test tRenewableFix isa PowerSystems.Component
    tRenewableDispatch = RenewableDispatch(nothing)
    @test tRenewableDispatch isa PowerSystems.Component
    tRenewableDispatch = RenewableDispatch(nothing)
    @test tRenewableDispatch isa PowerSystems.Component
end

@testset "Storage Constructors" begin
    tStorage = GenericBattery(nothing)
    @test tStorage isa PowerSystems.Component
end

@testset "Load Constructors" begin
    tPowerLoad = PowerLoad(nothing)
    @test tPowerLoad isa PowerSystems.Component
    tPowerLoadPF = PowerLoadPF(nothing)
    @test tPowerLoadPF isa PowerSystems.Component
    tPowerLoad = PowerLoad("init", true, Bus(nothing), nothing, 0.0, 0.0, 0.0, 0.0)
    @test tPowerLoad isa PowerSystems.Component
    tPowerLoadPF = PowerLoadPF("init", true, Bus(nothing), nothing, 0.0, 0.0, 1.0)
    @test tPowerLoadPF isa PowerSystems.Component
    tLoad = InterruptibleLoad(nothing)
    @test tLoad isa PowerSystems.Component
end

@testset "Branch Constructors" begin
    tLine = Line(nothing)
    @test tLine isa PowerSystems.Component
    tMonitoredLine = MonitoredLine(nothing)
    @test tMonitoredLine isa PowerSystems.Component
    tHVDCLine = HVDCLine(nothing)
    @test tHVDCLine isa PowerSystems.Component
    tVSCDCLine = VSCDCLine(nothing)
    @test tVSCDCLine isa PowerSystems.Component
    tTransformer2W = Transformer2W(nothing)
    @test tTransformer2W isa PowerSystems.Component
    tTapTransformer = TapTransformer(nothing)
    @test tTapTransformer isa PowerSystems.Component
    tPhaseShiftingTransformer = PhaseShiftingTransformer(nothing)
    @test tPhaseShiftingTransformer isa PowerSystems.Component
end

@testset "Service Constructors" begin
    #tProportionalReserve = ProportionalReserve(nothing)
    #@test tProportionalReserve isa PowerSystems.Service
    tStaticReserve = StaticReserve(nothing)
    @test tStaticReserve isa PowerSystems.Service
end

@testset "Forecast Constructors" begin
    tg = RenewableFix(nothing)
    forecast_data = PowerSystems.TimeSeries.TimeArray([DateTime("01-01-01"), DateTime("01-01-01")+Hour(1)], [1.0, 1.0])
    #Deterministic Tests
    tDeterministicForecast = PSY.Deterministic(tg,"scalingfactor", Hour(1),DateTime("01-01-01"),24)
    @test tDeterministicForecast isa PowerSystems.Forecast
    tDeterministicForecast = PSY.Deterministic(tg,"scalingfactor", forecast_data)
    @test tDeterministicForecast isa PowerSystems.Forecast
    #Probabilistic Tests
    tProbabilisticForecast = PSY.Probabilistic(tg,"scalingfactor", Hour(1), DateTime("01-01-01"),[0.5, 0.5], 24)
    @test  tProbabilisticForecast isa PowerSystems.Forecast
    tProbabilisticForecast = PSY.Probabilistic(tg,"scalingfactor", [1.0], forecast_data)
    @test  tProbabilisticForecast isa PowerSystems.Forecast
    #Scenario Tests
    tScenarioForecast = PSY.ScenarioBased(tg, "scalingfactor", Hour(1), DateTime("01-01-01"), 2, 24)
    @test  tScenarioForecast isa PowerSystems.Forecast
    tScenarioForecast = PSY.ScenarioBased(tg,"scalingfactor",forecast_data)
    @test  tScenarioForecast isa PowerSystems.Forecast
end


text:
# add more julia processes
using Distributed
nprocs() < 2 ? addprocs(2-nprocs()) : nothing

# tell Julia that you want to use these modules/namespaces
using RoME, Distributions
## Inter-operating visualization packages for Caesar/RoME/IncrementalInference exist
# using RoMEPlotting
using Compose
# Using this for trend plotting.
using Gadfly
# Using for capturing data.
using DataFrames

numVariables = 300
solveEveryNVariables = 20
lagLength = 30

# Standard Hexagonal example for totalIterations - solve every iterationsPerSolve iterations.
function runHexagonalExample(fg::AbstractDFG, totalIterations::Int, iterationsPerSolve::Int)::DataFrame
    # Add the first pose :x0
    addVariable!(fg, :x0, Pose2)

    # Add at a fixed location PriorPose2 to pin :x0 to a starting location
    addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3), 0.01*Matrix{Float64}(LinearAlgebra.I, 3,3))))

    # Add a landmark l1
    addVariable!(fg, :l1, Point2, tags=["LANDMARK"])

    # Drive around in a hexagon a number of times
    solveTimes = DataFrame(GraphSize = [], TimeBuildBayesTree = [], TimeSolveGraph = [])
    for i in 0:totalIterations
        psym = Symbol("x$i")
        nsym = Symbol("x$(i+1)")
        @info "Adding pose $nsym..."
        addVariable!(fg, nsym, Pose2)
        pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal( [0.1;0.1;0.1].^2 ) )))
        @info "Adding odometry factor between $psym -> $nsym..."
        addFactor!(fg, [psym;nsym], pp )

        if i % 6 == 0
            @info "Creating factor between $psym and l1..."
            p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
            addFactor!(fg, [psym; :l1], p2br)
        end
        if i % iterationsPerSolve == 0 && i != 0
            @info "Performing inference!"
            if getSolverParams(fg).isfixedlag
                @info "Quasi fixed-lag is enabled (a feature currently in testing)!"
                fifoFreeze!(fg)
            end
            tBuild = @timed tree = wipeBuildNewTree!(fg)
            tInfer = @timed inferOverTree!(fg, tree, N=100)
            graphSize = length([ls(fg)[1]..., ls(fg)[2]...])
            push!(solveTimes, (graphSize, tBuild[2], tInfer[2]))
        end
    end
    return solveTimes
end

# start with an empty factor graph object
fg = initfg()
# DO NOT enable fixed-lag operation
solverTimesForBatch = runHexagonalExample(fg, numVariables, solveEveryNVariables)
# savejld(fg, "batchSolve.jld2")

#############################
##### ---- Fixed lag mode ---
#############################

fgFixedLag = initfg()
fgFixedLag.solverParams.isfixedlag = true
fgFixedLag.solverParams.qfl = lagLength

solverTimesFixedLag = runHexagonalExample(fgFixedLag, numVariables, solveEveryNVariables)
# savejld(fgFixedLag, "fixedSolve.jld2")

#### Visualization

# Plot the many iterations to see that it succeeded.
# Batch
# drawPosesLandms(fg)

# Fixed lag
# drawPosesLandms(fgFixedLag)

#### PLOTTING of Results

using Gadfly
using Colors
using CSV

# Make a clean dataset
rename!(solverTimesForBatch, :TimeBuildBayesTree => :Batch_BayedBuild, :TimeSolveGraph => :Batch_SolveGraph);
rename!(solverTimesFixedLag, :TimeBuildBayesTree => :FixedLag_BayedBuild, :TimeSolveGraph => :FixedLag_SolveGraph);
timingMerged = DataFrames.join(solverTimesForBatch, solverTimesFixedLag, on=:GraphSize)
CSV.write("timing_comparison.csv", timingMerged)

# PP = []
# push!(PP, Gadfly.layer(x=timingMerged[:GraphSize], y=timingMerged[:FixedLag_SolveGraph], Geom.path, Theme(default_color=colorant"green"))[1]);
# push!(PP, Gadfly.layer(x=timingMerged[:GraphSize], y=timingMerged[:Batch_SolveGraph], Geom.path, Theme(default_color=colorant"magenta"))[1]);
#
# plt = Gadfly.plot(PP...,
#     Guide.title("Solving Time vs. Iteration for Fixed-Lag Operation"),
#     Guide.xlabel("Solving Iteration"),
#     Guide.ylabel("Solving Time (seconds)"),
#     Guide.manual_color_key("Legend", ["fixed", "batch"], ["green", "magenta"]))
# Gadfly.draw(PNG("results_comparison.png", 12cm, 15cm), plt)

#### EXPORT the data
using JLD2, FileIO
@save "results.jld2" timingMerged


text:
# Copyright 2022 RelationalAI, Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# The primary API level interface to the RAI REST API. These entry points
# provide convient Juilia language bindings to the corresponding opeartions,
# but a minimum of additional functionality. The purpose is to present service
# functionality as direclty as possible, but in a way that is natural for the
# Julia language.

import JSON3
import Arrow
using Base.Threads: @spawn

using Mocking: Mocking, @mock  # For unit testing, by mocking API server responses

const PATH_DATABASE = "/database"
const PATH_ENGINE = "/compute"
const PATH_OAUTH_CLIENTS = "/oauth-clients"
const PATH_TRANSACTION = "/transaction"
const PATH_ASYNC_TRANSACTIONS = "/transactions"
const PATH_USERS = "/users"
const ARROW_CONTENT_TYPE = "application/vnd.apache.arrow.stream"

struct HTTPError <: Exception
    status_code::Int
    status_text::String
    details::Union{String,Nothing}
    HTTPError(status_code) = new(status_code, HTTP.statustext(status_code), nothing)
    HTTPError(status_code, details) = new(status_code, HTTP.statustext(status_code), details)
end

function Base.show(io::IO, e::HTTPError)
    println(io, "$(e.status_code) $(e.status_text)")
    isnothing(e.details) && return
    try
        println(io, JSON3.read(e.details))
    catch
        println(io, e.details)
    end
end

# Returns a `Dict` constructed from the given pairs, filtering out pairs where
# the value is `nothing`.
function _filter(pairs::Pair...)
    d = Dict((k => v) for (k, v) in pairs if !isnothing(v))
    return length(d) > 0 ? d : nothing
end

# Returns a URL constructed from context settings and the given path.
function _mkurl(ctx::Context, path)
    return "$(ctx.scheme)://$(ctx.host):$(ctx.port)$path"
end

function _print_request(method, path, query, body)
    println("$method $path")
    !isnothing(query) && for (k, v) in query
        println("$k: $v")
    end
    !isnothing(body) && println(String(body))
end

function _request(ctx::Context, method, path; query = nothing, body = UInt8[], kw...)
    # _print_request(method, path, query, body);
    try
        rsp = request(ctx, method, _mkurl(ctx, path); query = query, body = body, kw...)
        return JSON3.read(rsp.body)
    catch e
        if e isa HTTP.ExceptionRequest.StatusError
            throw(HTTPError(e.status, String(e.response.body)))
        else
            rethrow()
        end
    end
end

_delete(ctx::Context, path; body = nothing, kw...) =
    _request(ctx, "DELETE", path; body = body, kw...)

_get(ctx::Context, path; query = nothing, kw...) =
    _request(ctx, "GET", path; query = query, kw...)

_patch(ctx::Context, path; body = nothing, kw...) =
    _request(ctx, "PATCH", path; body = body, kw...)

_post(ctx::Context, path; body = nothing, kw...) =
    _request(ctx, "POST", path; body = body, kw...)

_put(ctx::Context, path; body = nothing, kw...) =
    _request(ctx, "PUT", path; body = body, kw...)

function create_engine(ctx::Context, engine::AbstractString; size = nothing, kw...)
    isnothing(size) && (size = "XS")
    data = Dict("region" => ctx.region, "name" => engine, "size" => size)
    return _put(ctx, PATH_ENGINE; body = JSON3.write(data), kw...)
end

function create_oauth_client(ctx::Context, name::AbstractString, permissions; kv...)
    isnothing(permissions) && (permissions = [])
    data = Dict("name" => name, "permissions" => permissions)
    return _post(ctx, PATH_OAUTH_CLIENTS; body = JSON3.write(data), kv...)
end

function create_user(ctx::Context, email::AbstractString, roles = nothing; kw...)
    isnothing(roles) && (roles = [])
    data = Dict("email" => email, "roles" => roles)
    return _post(ctx, PATH_USERS; body = JSON3.write(data), kw...)
end

function delete_database(ctx::Context, database::AbstractString; kw...)
    data = Dict("name" => database)
    return _delete(ctx, PATH_DATABASE; body = JSON3.write(data), kw...)
end

function delete_engine(ctx::Context, engine::AbstractString; kw...)
    data = Dict("name" => engine)
    return _delete(ctx, PATH_ENGINE; body = JSON3.write(data), kw...)
end

function delete_model(ctx::Context, database::AbstractString, engine::AbstractString, model::AbstractString; kw...)
    tx = Transaction(ctx.region, database, engine, "OPEN"; readonly = false)
    actions = [_make_delete_models_action([model])]
    return _post(ctx, PATH_TRANSACTION; query = query(tx), body = body(tx, actions...), kw...)
end

function delete_oauth_client(ctx::Context, id::AbstractString; kw...)
    return _delete(ctx, joinpath(PATH_OAUTH_CLIENTS, id); kw...)
end

function delete_user(ctx::Context, userid::AbstractString; kw...)
    return _delete(ctx, joinpath(PATH_USERS, userid); kw...)
end

function disable_user(ctx::Context, userid::AbstractString; kw...)
    return update_user(ctx, userid; status = "INACTIVE", kw...)
end

function enable_user(ctx::Context, userid::AbstractString; kw...)
    return update_user(ctx, userid; status = "ACTIVE", kw...)
end

function get_engine(ctx::Context, engine::AbstractString; kw...)
    query = Dict("name" => engine, "deleted_on" => "")
    rsp = _get(ctx, PATH_ENGINE; query = query, kw...).computes
    length(rsp) == 0 && throw(HTTPError(404))
    return rsp[1]
end

function get_database(ctx::Context, database::AbstractString; kw...)
    query = Dict("name" => database)
    rsp = _get(ctx, PATH_DATABASE; query = query, kw...).databases
    length(rsp) == 0 && throw(HTTPError(404))
    return rsp[1]
end

# todo: move to rel query
function get_model(ctx::Context, database::AbstractString, engine::AbstractString, name::AbstractString; kw...)
    models = _list_models(ctx, database, engine; kw...)
    for model in models
        model["name"] == name && return model["value"]
    end
    throw(HTTPError(404))
end

function get_oauth_client(ctx::Context, id::AbstractString; kw...)
    return _get(ctx, joinpath(PATH_OAUTH_CLIENTS, id); kw...).client
end

# Returns the user with the given email.
function find_user(ctx::Context, email::AbstractString; kw...)
    rsp = list_users(ctx)
    for item in rsp
        item["email"] == email && return item
    end
    return nothing
end

function get_user(ctx::Context, userid::AbstractString; kw...)
    return _get(ctx, joinpath(PATH_USERS, userid); kw...).user
end

function list_databases(ctx::Context; state = nothing, kw...)
    query = _filter("state" => state)
    return _get(ctx, PATH_DATABASE; query = query, kw...).databases
end

function list_engines(ctx::Context; state = nothing, kw...)
    query = _filter("state" => state)
    return _get(ctx, PATH_ENGINE; query = query, kw...).computes
end

function list_oauth_clients(ctx::Context; kw...)
    return _get(ctx, PATH_OAUTH_CLIENTS; kw...).clients
end

function list_users(ctx::Context; kw...)
    return _get(ctx, PATH_USERS; kw...).users
end

function update_user(ctx::Context, userid::AbstractString; status = nothing, roles = nothing, kw...)
    data = _filter("status" => status, "roles" => roles)
    return _patch(ctx, joinpath(PATH_USERS, userid); body = JSON3.write(data), kw...)
end

"""
    Transaction

Represents the required and optional parameters associated with a transaction
request.
"""
mutable struct Transaction
    region::String
    database::String
    engine::Union{String,Nothing}
    mode::String
    abort::Bool
    nowait_durable::Bool
    readonly::Bool
    source::Union{String,Nothing}
    version::Int

    Transaction(region, database, engine, mode; source = nothing, readonly = false) =
        new(region, database, engine, mode, false, false, readonly, source, 0)
end

# Returns the serialized request body for the given transaction.
function body(tx::Transaction, actions...)::String
    data = _filter(
        "type" => "Transaction",
        "dbname" => tx.database,
        "mode" => !isnothing(tx.mode) ? tx.mode : "OPEN",
        "computeName" => tx.engine,
        "source_dbname" => tx.source,
        "abort" => tx.abort,
        "readonly" => tx.readonly,
        "nowaite_durable" => tx.nowait_durable,
        "version" => tx.version,
        "actions" => _make_actions(actions...))
    return JSON3.write(data)
end

# Returns the request query params for the given transaction.
function query(tx::Transaction)
    return _filter(
        "dbname" => tx.database,
        "compute_name" => tx.engine,
        "open_mode" => tx.mode,
        "region" => tx.region,
        "source_dbname" => tx.source)
end

function _create_mode(source, overwrite)
    if !isnothing(source)
        return overwrite ? "CLONE_OVERWRITE" : "CLONE"
    else
        return overwrite ? "CREATE_OVERWRITE" : "CREATE"
    end
end

# Wraps each of the given actions with a LabeledAction.
function _make_actions(actions...)
    result = []
    for (i, action) in enumerate(actions)
        item = Dict(
            "name" => "action$i",
            "type" => "LabeledAction",
            "action" => action)
        push!(result, item)
    end
    return result
end

function _make_delete_models_action(models::Vector)
    return Dict(
        "type" => "ModifyWorkspaceAction",
        "delete_source" => models)
end

function _make_load_model_action(name, model)
    return Dict(
        "type" => "InstallAction",
        "sources" => [_make_query_source(name, model)])
end

function _make_list_models_action()
    return Dict("type" => "ListSourceAction")
end

function _make_list_edb_action()
    return Dict("type" => "ListEdbAction")
end

function _make_query_action(source, ::Nothing)
    return Dict(
        "type" => "QueryAction",
        "source" => _make_query_source("query", source),
        "persist" => [], # todo: remove
        "inputs" => [],
        "outputs" => [])
end

function _make_query_action(source, inputs::Dict)
    return Dict(
        "type" => "QueryAction",
        "source" => _make_query_source("query", source),
        "persist" => [], # todo: remove
        "inputs" => [_make_query_action_input(k, v) for (k, v) in inputs],
        "outputs" => [])
end

function _make_query_action_input(name, value)
    return Dict(
        "type" => "Relation",
        "columns" => [[value]],
        "rel_key" => _make_relkey(name, _reltype(value)))
end

function _make_relkey(name, key)
    return Dict(
        "type" => "RelKey",
        "name" => name,
        "keys" => [key],
        "values" => [])
end

function _make_query_source(name, model)
    return Dict(
        "type" => "Source",
        "name" => name,
        "path" => "",
        "value" => model)
end

function _reltype(_::AbstractString)
    return "RAI_VariableSizeStrings.VariableSizeString"
end

function create_database(ctx::Context, database::AbstractString, engine::AbstractString; source = nothing, overwrite = false, kw...)
    mode = _create_mode(source, overwrite)
    tx = Transaction(ctx.region, database, engine, mode; source = source)
    return _post(ctx, PATH_TRANSACTION; query = query(tx), body = body(tx), kw...)
end

# Execute the given query string, using any given optioanl query inputs.
# todo: consider create_transaction
# todo: consider create_transaction to better align with future transaction
#   resource model
function exec_v1(ctx::Context, database::AbstractString, engine::AbstractString, source; inputs = nothing, readonly = false, kw...)
    source isa IO && (source = read(source, String))
    tx = Transaction(ctx.region, database, engine, "OPEN"; readonly = readonly)
    data = body(tx, _make_query_action(source, inputs))
    return _post(ctx, PATH_TRANSACTION; query = query(tx), body = data, kw...)
end
# todo: when we have async transactions, add a variation that dispatches and
#   waits .. consider creating two entry points for readonly and readwrite.

"""
    exec(ctx, "database", "engine", "query source"; kwargs...)

Synchronously execute a provided query string in the supplied database. This function
creates a Transaction using the supplied engine, and then polls the Transaction until it has
completed, and returns a Dict holding the Transaction resource and its results and metadata.

## Keyword Arguments:
- `readonly = false`: If true, this is a "read-only query", and the effects of this
   transaction will not be committed to the database.
- `inputs`: Optional dictionary of input pairs, mapping a relation name to a value.
   (Deprecated - the format of the inputs Dict will change in upcoming releases.)

# Examples:
```julia
julia> exec(ctx, "my_database", "my_engine", "2 + 2")
Dict{String, Any} with 4 entries:
  "metadata"    => Union{}[]
  "problems"    => Union{}[]
  "results"     => Pair{String, Arrow.Table}["/:output/Int64"=>Arrow.Table with 1 r…
  "transaction" => {…

julia> exec(ctx, "my_database", "my_engine", \"""
           def insert:my_relation = 1, 2, 3
           \""",
           readonly = false,
       )
Dict{String, Any} with 4 entries:
  "metadata"    => Union{}[]
  "problems"    => Union{}[]
  "results"     => Any[]
  "transaction" => {…
```
"""
function exec(ctx::Context, database::AbstractString, engine::AbstractString, source; inputs = nothing, readonly = false, kw...)
    txn = exec_async(ctx, database, engine, source; inputs=inputs, readonly=readonly, kw...)
    try
        backoff = Base.ExponentialBackOff(
                n = typemax(Int),
                first_delay = 0.5,
                factor = 1.1,
                max_delay = 120,  # 2 min
            )
        for duration in backoff
            transaction_is_done(txn) && break

            txn = get_transaction(ctx, transaction_id(txn))
            sleep(duration)
        end
        if haskey(txn, "results")
            return txn
        else
            id = transaction_id(txn)
            t = @spawn get_transaction(ctx, id)
            m = @spawn get_transaction_metadata(ctx, id)
            p = @spawn get_transaction_problems(ctx, id)
            r = @spawn get_transaction_results(ctx, id)
            return Dict(
                "transaction" => fetch(t),
                "metadata" => fetch(m),
                "problems" => fetch(p),
                "results" => fetch(r),
            )
        end
    catch
        @info "TXN" txn
        # Always print out the transaction id so that users can still get the txn ID even
        # if there's an error during polling (such as an InterruptException).
        #@info """Exception while polling for transaction:\n"id": $(repr(transaction_id(txn)))"""
        rethrow()
    end
end

function exec_async(ctx::Context, database::AbstractString, engine::AbstractString, source; inputs = nothing, readonly = false, kw...)
    source isa IO && (source = read(source, String))
    tx_body = Dict(
        "dbname" => database,
        "engine_name" => engine,
        "query" => source,
        #"nowait_durable" => self.nowait_durable, # TODO: currently unsupported
        "readonly" => readonly,
        # "sync_mode" => "async"
    )
    if inputs !== nothing
        tx_body["v1_inputs"] = [_make_query_action_input(k, v) for (k, v) in inputs]
    end
    body = JSON3.write(tx_body)
    path = _mkurl(ctx, PATH_ASYNC_TRANSACTIONS)
    rsp = @mock request(ctx, "POST", path; body = body, kw...)
    return _parse_response(rsp)
end

function _parse_response(rsp)
    content_type = HTTP.header(rsp, "Content-Type")
    if lowercase(content_type) == "application/json"
        content = HTTP.body(rsp)
        # async mode
        return Dict(
            "transaction" => JSON3.read(content),
        )
    elseif occursin("multipart/form-data", lowercase(content_type))
        # sync mode
        return _parse_multipart_fastpath_sync_response(rsp)
    else
        error("Unknown response content-type, for response:\n$(rsp)")
    end
end

function get_transaction(ctx::Context, id::AbstractString; kw...)
    path = PATH_ASYNC_TRANSACTIONS * "/$id"
    rsp = _get(ctx, path; kw...)
    return rsp.transaction
end

function transaction_is_done(txn)
    if haskey(txn, "transaction")
        txn = txn["transaction"]
    end
    return txn["state"] ∈ ("COMPLETED", "ABORTED")
end

function transaction_id(txn)
    if haskey(txn, "transaction")
        txn = txn["transaction"]
    end
    return txn["id"]
end

function get_transaction_metadata(ctx::Context, id::AbstractString; kw...)
    path = PATH_ASYNC_TRANSACTIONS * "/$id/metadata"
    rsp = _get(ctx, path; kw...)
    return rsp
end

function get_transaction_problems(ctx::Context, id::AbstractString; kw...)
    path = PATH_ASYNC_TRANSACTIONS * "/$id/problems"
    rsp = _get(ctx, path; kw...)
    return rsp
end

function get_transaction_results(ctx::Context, id::AbstractString; kw...)
    path = PATH_ASYNC_TRANSACTIONS * "/$id/results"
    path = _mkurl(ctx, path)
    rsp = @mock request(ctx, "GET", path; kw...)
    content_type = HTTP.header(rsp, "Content-Type")
    if !occursin("multipart/form-data", content_type)
        throw(HTTPError(400, "Unexpected response content-type for rsp:\n$rsp"))
    end
    return _parse_multipart_results_response(rsp)
end

function _parse_multipart_fastpath_sync_response(msg)
    # TODO: in-place conversion to Arrow without copying the bytes.
    #   ... HTTP.parse_multipart_form() copies the bytes into IOBuffers.
    parts = _parse_multipart_form(msg)
    @assert parts[1].name == "transaction"
    @assert parts[2].name == "metadata"

    problems_idx = findfirst(p->p.name == "problems", parts)
    problems = JSON3.read(parts[problems_idx])
    results = _extract_multipart_results_response(parts)

    return Dict(
        "transaction" => JSON3.read(parts[1]),
        "metadata" => JSON3.read(parts[2]),
        "problems" => problems,
        "results" => results,
    )
end

function _parse_multipart_results_response(msg)
    # TODO: in-place conversion to Arrow without copying the bytes.
    #   ... HTTP.parse_multipart_form() copies the bytes into IOBuffers.
    parts = _parse_multipart_form(msg)
    return _extract_multipart_results_response(parts)
end
function _extract_multipart_results_response(parts)
    return [
        (part.name => Arrow.Table(part.data)) for part in parts
            if part.contenttype == ARROW_CONTENT_TYPE
    ]
end

function list_edbs(ctx::Context, database::AbstractString, engine::AbstractString; kw...)
    tx = Transaction(ctx.region, database, engine, "OPEN"; readonly = true)
    data = body(tx, _make_list_edb_action())
    rsp = _post(ctx, PATH_TRANSACTION; query = query(tx), body = data, kw...)
    length(rsp.actions) == 0 && return []
    return rsp.actions[1].result.rels
end

function _list_models(ctx::Context, database::AbstractString, engine::AbstractString; kw...)
    tx = Transaction(ctx.region, database, engine, "OPEN"; readonly = true)
    data = body(tx, _make_list_models_action())
    rsp = _post(ctx, PATH_TRANSACTION; query = query(tx), body = data, kw...).actions
    length(rsp) == 0 && return []
    return rsp[1].result.sources
end

function list_models(ctx::Context, database::AbstractString, engine::AbstractString; kw...)
    models = _list_models(ctx, database, engine; kw...)
    return [model["name"] for model in models]
end

function _gen_literal(value)
    return "$value"
end

function _gen_literal(value::Dict)
    items = ["$(_gen_literal(v)),$(_gen_literal(k))" for (k, v) in value]
    return "{" + join(items, ";") + "}"
end

function _gen_literal(value::String)
    s = replace(value, "'" => "\\'")
    return "'$s'"
end

_gen_literal(value::Symbol) = ":$value"

function _gen_literal(value::Vector)
    items = [_gen_literal(item) for item in value]
    return "{" + join(items, ",") + "}"
end

function _gen_config(name, value)
    isnothing(value) && return ""
    return "def config:syntax:$name=$(_gen_literal(value))"
end

function _gen_config(syntax::Dict)
    length(syntax) == 0 && return ""
    items = [_gen_config(k, v) for (k, v) in syntax if !isnothing(v)]
    return join(items, '\n') * '\n'
end

_gen_config(::Nothing) = ""

_read_data(d::AbstractString) = d
_read_data(d::IO) = read(d, String)

# todo: need to uniquify config and data so it doesn't conflict with those
#   names if they already exist in the database.
# todo: add support for config:path
function load_csv(
    ctx::Context, database, engine, relation, data;
    delim = nothing, header = nothing, header_row = nothing,
    escapechar = nothing, quotechar = nothing, kw...
)
    inputs = Dict("data" => _read_data(data))
    syntax = Dict{String,String}()
    syntax = _filter(
        "delim" => delim,
        "header" => header,
        "header_row" => header_row,
        "escapechar" => escapechar,
        "quotechar" => quotechar)
    source = _gen_config(syntax)
    source *= """def config:data = data
                 def insert:$relation = load_csv[config]"""
    return exec(ctx, database, engine, source; inputs = inputs, readonly = false, kw...)
end

# todo: need to uniquify config and data
# todo: add support for config:path
# todo: data should be string or io
function load_json(ctx::Context, database::AbstractString, engine::AbstractString, relation::AbstractString, data; kw...)
    inputs = Dict("data" => _read_data(data))
    source = """def config:data = data\n
                def insert:$relation = load_json[config]"""
    return exec(ctx, database, engine, source; inputs = inputs, readonly = false, kw...)
end

function load_model(ctx::Context, database::AbstractString, engine::AbstractString, models::Dict; kw...)
    tx = Transaction(ctx.region, database, engine, "OPEN"; readonly = false)
    actions = [_make_load_model_action(name, model) for (name, model) in models]
    return _post(ctx, PATH_TRANSACTION; query = query(tx), body = body(tx, actions...), kw...)
end

# --- utils -------------------------
# Patch for older versions of HTTP package that don't support parsing multipart responses:
if hasmethod(HTTP.MultiPartParsing.parse_multipart_form, (HTTP.Response,))
    # Available as of HTTP v0.9.18:
    _parse_multipart_form = HTTP.MultiPartParsing.parse_multipart_form
else
    # This function is copied directly from this PR: https://github.com/JuliaWeb/HTTP.jl/pull/817
    function _parse_multipart_form(msg::HTTP.Message)
        # parse boundary from Content-Type
        m = match(r"multipart/form-data; boundary=(.*)$", msg["Content-Type"])
        m === nothing && return nothing

        boundary_delimiter = m[1]

        # [RFC2046 5.1.1](https://tools.ietf.org/html/rfc2046#section-5.1.1)
        length(boundary_delimiter) > 70 && error("boundary delimiter must not be greater than 70 characters")

        return HTTP.MultiPartParsing.parse_multipart_body(HTTP.payload(msg), boundary_delimiter)
    end
end
# -----------------------------------

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
### SIMULATION FUNCTIONS ###

"""
    simulate_step(filter::AbstractFilter, x::AbstractVector, u::AbstractVector, rng::AbstractRNG=Random.GLOBAL_RNG)

Run a step of simulation starting at state x, taking action u, and using the
motion and measurement equations specified by the filter.
"""
function simulate_step(filter::AbstractFilter, x::AbstractVector{<:Number},
                        u::AbstractVector{<:Number}, rng::AbstractRNG=Random.GLOBAL_RNG)

    xn = predict(filter.d, x, u, rng)

    yn = measure(filter.o, xn, u, rng)

    return xn, yn
end

"""
    simulation(filter::AbstractFilter, b0::GaussianBelief,
                action_sequence::Vector{AbstractVector}})

Run a simulation to get positions and measurements. Samples starting point from
GaussianBelief b0, the runs action_sequence with additive gaussian noise all
specified by AbstractFilter filter to return a simulated state and measurement
history.
"""
function simulation(filter::AbstractFilter, b0::GaussianBelief,
                    action_sequence::Vector{<:AbstractArray},
                    rng::AbstractRNG = Random.GLOBAL_RNG)

    # make initial state
    s0 = rand(rng, b0)

    # simulate action sequence
    state_history = [s0]
    measurement_history = Vector{AbstractVector{typeof(s0[1])}}()
    for u in action_sequence
        xn, yn = simulate_step(filter, state_history[end], u)
        push!(state_history, xn)
        push!(measurement_history, yn)
    end

    ## TODO: Is it worth separating out start state from state_history
    return state_history, measurement_history
end

"""
    run_filter(filter::AbstractFilter, b0::GaussianBelief, action_history::Vector{AbstractVector},
            measurement_history::Vector{AbstractVector})

Given an initial belief b0, matched-size arrays for action and measurement
histories and a filter, update the beliefs using the filter, and return a
vector of all beliefs.
"""
function run_filter(filter::AbstractFilter, b0::GaussianBelief, action_history::Vector{A},
            measurement_history::Vector{B}) where {A<:AbstractVector, B<:AbstractVector}

        # assert matching action and measurement sizes
        @assert length(action_history) == length(measurement_history)

        # initialize belief vector
        beliefs = [b0]

        # iterate through and update beliefs
        for (u, y) in zip(action_history, measurement_history)
            bn = update(filter, beliefs[end], u, y)
            push!(beliefs, bn)
        end

        return beliefs
end

"""
    unpack(belief_history::Vector{<:GaussianBelief};
        dims::Vector{Int}=[])

Given a history of beliefs, return an unpacked (time steps, state dim)-sized array of
predicted means and a (time steps, state dim, state dim)-sized array of
covariances. One can optionally specify dimensions indices dims to output
reduced state information.
"""
function unpack(belief_history::Vector{<:GaussianBelief};
    dims::Vector{Int}=Vector{Int}())

    # set default to condense all dimensions
    if length(dims) == 0
        dims = collect(1:length(belief_history[1].μ))
    end

    # set output sizes
    μ = zeros(typeof(belief_history[1].μ[1]),
        length(belief_history), length(dims))
    Σ = zeros(typeof(belief_history[1].μ[1]),
        length(belief_history), length(dims), length(dims))

    # iterate over belief_history and place elements appropriately
    for (i, belief) in enumerate(belief_history)
        μ[i,:] = belief.μ[dims]
        Σ[i,:,:] = belief.Σ[dims,dims]
    end

    return μ, Σ
end
### TODO: add in-place update!, predict!, and measure! functions


text:
export trouve_intervalle
function trouve_intervalle(h :: AbstractNLPModel,
                           t₀ :: Float64,
                           inc0 :: Float64;
                           verbose :: Bool = false)

        iter = 1
        h₀ = obj(h, [t₀])[1]
        g₀ = grad(h, [t₀])[1]
        sd = -sign(g₀)
        inc = inc0
        t₁ = t₀ + sd * inc

        h₁ = obj(h, [t₁])[1]
        g₁ = grad(h, [t₁])[1]

        verbose &&
            @printf("iter t₀        g₀        h₀         t₁        ")
        verbose &&
            @printf("g₀        h₀\n")
        verbose &&
            @printf("%4d %7.2e %7.2e  %7.2e  %7.2e  %7.2e  %7.2e \n",
                iter,  t₀, g₀, h₀, t₁, g₀, h₀)

        while (g₁ * sd < 0.0) && (h₁ < h₀) && (iter < 30)
          inc = inc * 4
          t₀ = t₁
          h₀ = h₁
          g₀ = g₁
          t₁ = t₀ + sd * inc
          h₁ = obj(h, [t₁])[1]
          g₁ = grad(h, [t₁])[1]
          verbose && @printf("%4d %7.2e %7.2e  %7.2e  %7.2e  %7.2e  %7.2e \n",
                              iter,  t₀, g₀, h₀, t₁, g₀, h₀)
          iter = iter + 1
        end

        while (g₁ * sd < 0.0) & (iter < 30)
          tₘ = (t₁ + t₀) / 2
          hₘ = obj(h, [tₘ])[1]
          gₘ = grad(h, [tₘ])[1]
          if gₘ * sd>0
            t₁ = tₘ
            g₁ = gₘ
            h₁ = hₘ
          else
            if hₘ < h₀
              t₀ = tₘ
              h₀ = hₘ
              g₀ = gₘ
            else
              t₁ = tₘ
              h₁ = hₘ
              g₁ = gₘ
            end
          end
          verbose &&
            @printf("%4d %7.2e %7.2e  %7.2e  %7.2e  %7.2e  %7.2e \n",
                    iter,  t₀, g₀, h₀, t₁, g₀, h₀)
          iter = iter + 1
        end

        a = min(t₀, t₁)
        b = max(t₀, t₁)

        return (a, b)

end


text:
# Julia wrapper for header: /opt/ffmpeg/include/libswscale/swscale.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    #swscale_version,
    swscale_configuration,
    swscale_license,
    sws_getCoefficients,
    sws_isSupportedInput,
    sws_isSupportedOutput,
    sws_isSupportedEndiannessConversion,
    sws_alloc_context,
    sws_init_context,
    sws_freeContext,
    sws_getContext,
    sws_scale,
    sws_setColorspaceDetails,
    sws_getColorspaceDetails,
    sws_allocVec,
    sws_getGaussianVec,
    sws_getConstVec,
    sws_getIdentityVec,
    sws_scaleVec,
    sws_normalizeVec,
    sws_convVec,
    sws_addVec,
    sws_subVec,
    sws_shiftVec,
    sws_cloneVec,
    sws_printVec2,
    sws_freeVec,
    sws_getDefaultFilter,
    sws_freeFilter,
    sws_getCachedContext,
    sws_convertPalette8ToPacked32,
    sws_convertPalette8ToPacked24,
    sws_get_class

function swscale_version()
    ccall((:swscale_version,libswscale),Uint32,())
end

function swscale_configuration()
    ccall((:swscale_configuration,libswscale),Ptr{Uint8},())
end

function swscale_license()
    ccall((:swscale_license,libswscale),Ptr{Uint8},())
end

function sws_getCoefficients(colorspace::Integer)
    ccall((:sws_getCoefficients,libswscale),Ptr{Cint},(Cint,),colorspace)
end

function sws_isSupportedInput(pix_fmt::AVPixelFormat)
    ccall((:sws_isSupportedInput,libswscale),Cint,(AVPixelFormat,),pix_fmt)
end

function sws_isSupportedOutput(pix_fmt::AVPixelFormat)
    ccall((:sws_isSupportedOutput,libswscale),Cint,(AVPixelFormat,),pix_fmt)
end

function sws_isSupportedEndiannessConversion(pix_fmt::AVPixelFormat)
    ccall((:sws_isSupportedEndiannessConversion,libswscale),Cint,(AVPixelFormat,),pix_fmt)
end

function sws_alloc_context()
    ccall((:sws_alloc_context,libswscale),Ptr{SwsContext},())
end

function sws_init_context(sws_context,srcFilter,dstFilter)
    ccall((:sws_init_context,libswscale),Cint,(Ptr{SwsContext},Ptr{SwsFilter},Ptr{SwsFilter}),sws_context,srcFilter,dstFilter)
end

function sws_freeContext(swsContext)
    ccall((:sws_freeContext,libswscale),Void,(Ptr{SwsContext},),swsContext)
end

function sws_getContext(srcW::Integer,srcH::Integer,srcFormat::AVPixelFormat,dstW::Integer,dstH::Integer,dstFormat::AVPixelFormat,flags::Integer,srcFilter,dstFilter,param)
    ccall((:sws_getContext,libswscale),Ptr{SwsContext},(Cint,Cint,AVPixelFormat,Cint,Cint,AVPixelFormat,Cint,Ptr{SwsFilter},Ptr{SwsFilter},Ptr{Cdouble}),srcW,srcH,srcFormat,dstW,dstH,dstFormat,flags,srcFilter,dstFilter,param)
end

function sws_scale(c,srcSlice,srcStride,srcSliceY::Integer,srcSliceH::Integer,dst,dstStride)
    ccall((:sws_scale,libswscale),Cint,(Ptr{SwsContext},Ptr{Ptr{Uint8}},Ptr{Cint},Cint,Cint,Ptr{Ptr{Uint8}},Ptr{Cint}),c,srcSlice,srcStride,srcSliceY,srcSliceH,dst,dstStride)
end

function sws_setColorspaceDetails(c,inv_table,srcRange::Integer,table,dstRange::Integer,brightness::Integer,contrast::Integer,saturation::Integer)
    ccall((:sws_setColorspaceDetails,libswscale),Cint,(Ptr{SwsContext},Ptr{Cint},Cint,Ptr{Cint},Cint,Cint,Cint,Cint),c,inv_table,srcRange,table,dstRange,brightness,contrast,saturation)
end

function sws_getColorspaceDetails(c,inv_table,srcRange,table,dstRange,brightness,contrast,saturation)
    ccall((:sws_getColorspaceDetails,libswscale),Cint,(Ptr{SwsContext},Ptr{Ptr{Cint}},Ptr{Cint},Ptr{Ptr{Cint}},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),c,inv_table,srcRange,table,dstRange,brightness,contrast,saturation)
end

function sws_allocVec(length::Integer)
    ccall((:sws_allocVec,libswscale),Ptr{SwsVector},(Cint,),length)
end

function sws_getGaussianVec(variance::Cdouble,quality::Cdouble)
    ccall((:sws_getGaussianVec,libswscale),Ptr{SwsVector},(Cdouble,Cdouble),variance,quality)
end

function sws_getConstVec(c::Cdouble,length::Integer)
    ccall((:sws_getConstVec,libswscale),Ptr{SwsVector},(Cdouble,Cint),c,length)
end

function sws_getIdentityVec()
    ccall((:sws_getIdentityVec,libswscale),Ptr{SwsVector},())
end

function sws_scaleVec(a,scalar::Cdouble)
    ccall((:sws_scaleVec,libswscale),Void,(Ptr{SwsVector},Cdouble),a,scalar)
end

function sws_normalizeVec(a,height::Cdouble)
    ccall((:sws_normalizeVec,libswscale),Void,(Ptr{SwsVector},Cdouble),a,height)
end

function sws_convVec(a,b)
    ccall((:sws_convVec,libswscale),Void,(Ptr{SwsVector},Ptr{SwsVector}),a,b)
end

function sws_addVec(a,b)
    ccall((:sws_addVec,libswscale),Void,(Ptr{SwsVector},Ptr{SwsVector}),a,b)
end

function sws_subVec(a,b)
    ccall((:sws_subVec,libswscale),Void,(Ptr{SwsVector},Ptr{SwsVector}),a,b)
end

function sws_shiftVec(a,shift::Integer)
    ccall((:sws_shiftVec,libswscale),Void,(Ptr{SwsVector},Cint),a,shift)
end

function sws_cloneVec(a)
    ccall((:sws_cloneVec,libswscale),Ptr{SwsVector},(Ptr{SwsVector},),a)
end

function sws_printVec2(a,log_ctx,log_level::Integer)
    ccall((:sws_printVec2,libswscale),Void,(Ptr{SwsVector},Ptr{AVClass},Cint),a,log_ctx,log_level)
end

function sws_freeVec(a)
    ccall((:sws_freeVec,libswscale),Void,(Ptr{SwsVector},),a)
end

function sws_getDefaultFilter(lumaGBlur::Cfloat,chromaGBlur::Cfloat,lumaSharpen::Cfloat,chromaSharpen::Cfloat,chromaHShift::Cfloat,chromaVShift::Cfloat,verbose::Integer)
    ccall((:sws_getDefaultFilter,libswscale),Ptr{SwsFilter},(Cfloat,Cfloat,Cfloat,Cfloat,Cfloat,Cfloat,Cint),lumaGBlur,chromaGBlur,lumaSharpen,chromaSharpen,chromaHShift,chromaVShift,verbose)
end

function sws_freeFilter(filter)
    ccall((:sws_freeFilter,libswscale),Void,(Ptr{SwsFilter},),filter)
end

function sws_getCachedContext(context,srcW::Integer,srcH::Integer,srcFormat::AVPixelFormat,dstW::Integer,dstH::Integer,dstFormat::AVPixelFormat,flags::Integer,srcFilter,dstFilter,param)
    ccall((:sws_getCachedContext,libswscale),Ptr{SwsContext},(Ptr{SwsContext},Cint,Cint,AVPixelFormat,Cint,Cint,AVPixelFormat,Cint,Ptr{SwsFilter},Ptr{SwsFilter},Ptr{Cdouble}),context,srcW,srcH,srcFormat,dstW,dstH,dstFormat,flags,srcFilter,dstFilter,param)
end

function sws_convertPalette8ToPacked32(src,dst,num_pixels::Integer,palette)
    ccall((:sws_convertPalette8ToPacked32,libswscale),Void,(Ptr{Uint8},Ptr{Uint8},Cint,Ptr{Uint8}),src,dst,num_pixels,palette)
end

function sws_convertPalette8ToPacked24(src,dst,num_pixels::Integer,palette)
    ccall((:sws_convertPalette8ToPacked24,libswscale),Void,(Ptr{Uint8},Ptr{Uint8},Cint,Ptr{Uint8}),src,dst,num_pixels,palette)
end

function sws_get_class()
    ccall((:sws_get_class,libswscale),Ptr{AVClass},())
end