text:
using FactCheck, Images, Color, FixedPointNumbers

facts("IO") do
    workdir = joinpath(tempdir(), "Images")
    if !isdir(workdir)
        mkdir(workdir)
    end

    context("Gray png") do
        a = rand(2,2)
        aa = convert(Array{Ufixed8}, a)
        fn = joinpath(workdir, "2by2.png")
        Images.imwrite(a, fn)
        b = Images.imread(fn)
        @fact convert(Array, b) --> aa
        Images.imwrite(aa, fn)
        b = Images.imread(fn)
        @fact convert(Array, b) --> aa
        aaimg = Images.grayim(aa)
        open(fn, "w") do file
            writemime(file, MIME("image/png"), aaimg, minpixels=0)
        end
        b = Images.imread(fn)
        @fact b --> aaimg
        aa = convert(Array{Ufixed16}, a)
        Images.imwrite(aa, fn)
        b = Images.imread(fn)
        @fact convert(Array, b) --> aa
        aa = Ufixed12[0.6 0.2;
                      1.4 0.8]
        open(fn, "w") do file
            writemime(file, MIME("image/png"), Images.grayim(aa), minpixels=0)
        end
        b = Images.imread(fn)
        @fact Images.data(b) --> Ufixed8[0.6 0.2;
                                        1.0 0.8]
    end

    context("Color") do
        fn = joinpath(workdir, "2by2.png")
        img = Images.colorim(rand(3,2,2))
        img24 = convert(Images.Image{RGB24}, img)
        Images.imwrite(img24, fn)
        b = Images.imread(fn)
        imgrgb8 = convert(Images.Image{RGB{Ufixed8}}, img)
        @fact Images.data(imgrgb8) --> Images.data(b)
    end

    context("Writemime's use of restrict") do
        abig = Images.grayim(rand(Uint8, 1024, 1023))
        fn = joinpath(workdir, "big.png")
        open(fn, "w") do file
            writemime(file, MIME("image/png"), abig, maxpixels=10^6)
        end
        b = Images.imread(fn)
        @fact Images.data(b) --> convert(Array{Ufixed8,2}, Images.data(Images.restrict(abig, (1,2))))
    end

    context("More writemime tests") do
        a = Images.colorim(rand(Uint8, 3, 2, 2))
        fn = joinpath(workdir, "2by2.png")
        open(fn, "w") do file
            writemime(file, MIME("image/png"), a, minpixels=0)
        end
        b = Images.imread(fn)
        @fact Images.data(b) --> Images.data(a)

        abig = Images.colorim(rand(Uint8, 3, 1021, 1026))
        fn = joinpath(workdir, "big.png")
        open(fn, "w") do file
            writemime(file, MIME("image/png"), abig, maxpixels=10^6)
        end
        b = Images.imread(fn)
        @fact Images.data(b) --> convert(Array{RGB{Ufixed8},2}, Images.data(Images.restrict(abig, (1,2))))

        # Issue #269
        abig = Images.colorim(rand(Uint16, 3, 1024, 1023))
        open(fn, "w") do file
            writemime(file, MIME("image/png"), abig, maxpixels=10^6)
        end
        b = Images.imread(fn)
        @fact Images.data(b) --> convert(Array{RGB{Ufixed8},2}, Images.data(Images.restrict(abig, (1,2))))
    end

    context("Colormap usage") do
        datafloat = reshape(linspace(0.5, 1.5, 6), 2, 3)
        dataint = round(Uint8, 254*(datafloat .- 0.5) .+ 1)  # ranges from 1 to 255
        # build our colormap
        b = RGB(0,0,1)
        w = RGB(1,1,1)
        r = RGB(1,0,0)
        cmaprgb = Array(RGB{Float64}, 255)
        f = linspace(0,1,128)
        cmaprgb[1:128] = [(1-x)*b + x*w for x in f]
        cmaprgb[129:end] = [(1-x)*w + x*r for x in f[2:end]]
        img = Images.ImageCmap(dataint, cmaprgb)
        Images.imwrite(img,joinpath(workdir,"cmap.jpg"))
    end

    context("Alpha") do
        c = reinterpret(Images.BGRA{Ufixed8}, [0xf0884422]'')
        fn = joinpath(workdir, "alpha.png")
    Images.imwrite(c, fn)
        C = Images.imread(fn)
        # @test C[1] == c[1]  # disabled because Travis has a weird, old copy of ImageMagick for which this fails (see #261)
        Images.imwrite(reinterpret(ARGB32, [0xf0884422]''), fn)
        D = Images.imread(fn)
        # @test D[1] == c[1]
    end

    context("3D TIFF (issue #307)") do
        A = Images.grayim(rand(0x00:0xff, 2, 2, 4))
        fn = joinpath(workdir, "3d.tif")
        Images.imwrite(A, fn)
        B = Images.imread(fn)
        @fact A --> B
    end

    context("Clamping (issue #256)") do
        A = grayim(rand(2,2))
        A[1,1] = -0.4
        fn = joinpath(workdir, "2by2.png")
        @fact_throws InexactError Images.imwrite(A, fn)
        Images.imwrite(A, fn, mapi=Images.mapinfo(Images.Clamp, A))
        B = Images.imread(fn)
        A[1,1] = 0
        @fact B --> map(Ufixed8, A)
    end

    @unix_only context("Reading from a stream (issue #312)") do
        fn = joinpath(workdir, "2by2.png")
        io = open(fn)
        img = Images.imread(io, Images.ImageMagick)
        close(io)
        @fact isa(img, Images.Image) --> true
    end
end


text:
# functions that interact with the accelerator
export IterActivation, AccuracyActivation, ImmediateActivation
const CA = COSMOAccelerators

"Activate accelerator after `start_iter` iterations of the main algorithm."
struct IterActivation <: AbstractActivationReason
  start_iter::Int64
  function IterActivation(start_iter::Int64)
    start_iter < 2 && ArgumentError("start_iter has to be at least 2.")
    new(start_iter)
  end
end

"Activate accelerator after accuracy of main algorithm <= `start_accuracy`."
struct AccuracyActivation <: AbstractActivationReason
  start_accuracy::Float64
  function AccuracyActivation(start_accuracy::Float64)
    start_accuracy >= 0. && ArgumentError("start_accuracy has to be a non-negative number.")
    new(start_accuracy)
  end
end

# dispatch on the activation reason
function check_activation!(ws::Workspace, activation_reason::ImmediateActivation, num_iter::Int64)
    if !ws.accelerator_active && num_iter >= 2
        ws.accelerator_active = true
    end
end

function check_activation!(ws::Workspace, activation_reason::IterActivation, num_iter::Int64)
    if !ws.accelerator_active && num_iter >= activation_reason.start_iter
        ws.accelerator_active = true
    end
end
check_activation!(ws::Workspace, activation_reason::AccuracyActivation, num_iter::Int64) = nothing

function check_activation!(ws::Workspace, activation_reason::AccuracyActivation, r::ResultInfo) where {T <: AbstractFloat}
    if !ws.accelerator_active
        tol = activation_reason.start_accuracy

        if isapprox_primal_feasible(r, tol, tol) && isapprox_dual_feasible(r, tol, tol)
            ws.accelerator_active = true
        end
    end
end

check_activation!(ws::Workspace, activation_reason::Union{IterActivation, ImmediateActivation}, r::ResultInfo) where {T <: AbstractFloat}= nothing

"""
	acceleration_pre!(accelerator, ws, num_iter)

A function that the accelerator can use before the nominal ADMM operator step.

In the case of an Anderson Accelerator this is used to calculate an accelerated candidate vector, that overwrites the current iterate.
"""
function acceleration_pre!(accelerator::CA.AbstractAccelerator, ws::Workspace{T}, num_iter::Int64) where {T <: AbstractFloat}
	COSMO.check_activation!(ws, ws.activation_reason, num_iter)
	if ws.accelerator_active
		if ws.settings.verbose_timing
			time_start = time()
			CA.update!(accelerator, ws.vars.w, ws.vars.w_prev, num_iter)
			ws.times.update_time += time() - time_start
			time_start = time()
			CA.accelerate!(ws.vars.w, ws.vars.w_prev, accelerator, num_iter)
			ws.times.accelerate_time += time() - time_start
		else
			CA.update!(accelerator, ws.vars.w, ws.vars.w_prev, num_iter)
			# overwrite w here
			CA.accelerate!(ws.vars.w, ws.vars.w_prev, accelerator, num_iter)
		end
	end
end
acceleration_pre!(accelerator::CA.EmptyAccelerator, ws::Workspace{T}, num_iter::Int64) where {T <: AbstractFloat} = nothing

"""
	acceleration_post!(accelerator, ws, num_iter, safeguarding_iter)

A function that the accelerator can use after the nominal ADMM operator step.

In the case of an Anderson Accelerator this is used to check the quality of the accelerated candidate vector and take measures if the vector is of bad quality.
"""
function acceleration_post!(accelerator::CA.AndersonAccelerator, ws::Workspace{T}, num_iter::Int64) where {T <: AbstractFloat}
	if ws.accelerator_active && CA.was_successful(accelerator)
		if ws.settings.safeguard
			# norm(w_prev - w, 2)
			# use accelerator.f here to get f = (x - g) as w_prev has been overwritten before this function call
			nrm_f = norm(accelerator.f, 2)
			nrm_tol = nrm_f * ws.settings.safeguard_tol

			# compute residual norm of accelerated point
			nrm_f_acc = compute_accelerated_res_norm!(accelerator.f, ws.vars.w, ws.vars.w_prev)

			# Safeguarding check: f(w_acc) = w_prev_acc - w_acc <= τ f(w_prev) = τ * (w_prev - w)
			# if safeguarding check fails, discard the accelerated candidate point and do a normal ADMM step instead
            if nrm_f_acc > nrm_tol
                CA.log!(accelerator, num_iter, :acc_guarded_declined)
				# don't use accelerated candidate point. Reset w = g_last
				reset_accelerated_vector!(ws.vars.w, ws.vars.w_prev, accelerator)
				m, n = ws.p.model_size
				ws.times.proj_time += admm_z!(ws.vars.s, ws.vars.w, ws.p.C, n)
				admm_x!(ws.vars.s, ws.ν, ws.s_tl, ws.ls, ws.sol, ws.vars.w, ws.kkt_solver, ws.p.q, ws.p.b, ws.ρvec, ws.settings.sigma, m, n)
				admm_w!(ws.vars.s, ws.x_tl, ws.s_tl, ws.vars.w, ws.settings.alpha, m, n);
				ws.safeguarding_iter += 1
            else
                CA.log!(accelerator, num_iter, :acc_guarded_accepted)
			end
        else
            CA.log!(accelerator, num_iter, :acc_unguarded)
		end
	end
end

acceleration_post!(accelerator::CA.AbstractAccelerator, ws::Workspace{T}, num_iter::Int64) where {T <: AbstractFloat} = nothing

"Computes the nonexpansive operator residual norm, i.e. f(w_prev) = || w_prev - w ||_2."
function compute_accelerated_res_norm!(f::AbstractVector{T}, w::AbstractVector{T}, w_prev::AbstractVector{T}) where {T <: AbstractFloat}
	@. f = w_prev - w
	return norm(f, 2)
end

"Reset ADMM variables to last non-accelerated point."
reset_accelerated_vector!(w::AbstractVector{T}, w_prev::AbstractVector{T}, aa::CA.AndersonAccelerator{T}) where {T <: AbstractFloat} = reset_accelerated_vector!(w, w_prev, aa.g_last)
function reset_accelerated_vector!(w::AbstractVector{T}, w_prev::AbstractVector{T}, g_last::AbstractVector{T}) where {T <: AbstractFloat}
		@. w_prev = g_last
		@. w = g_last
end

ext: jl


text:
"""
Function object to map points on the cube surface to the sphere of radius r.
We need sides of length 2 centred at the origin using the approach
in http://mathproofs.blogspot.com/2005/07/mapping-cube-to-sphere.html.

Construct arguments:
  r  : sphere radius
Function object arguments:
  xyz: 3D cartesian coordinates of a point on the cube surface
"""
struct MapCubeToSphere{T} <: Function
  radius::T
end

function (map::MapCubeToSphere{T})(xyz) where T
  x,y,z = xyz
  xₛ = x*sqrt(1.0-y^2/2-z^2/2+y^2*z^2/(3.0))
  yₛ = y*sqrt(1.0-z^2/2-x^2/2+x^2*z^2/(3.0))
  zₛ = z*sqrt(1.0-x^2/2-y^2/2+x^2*y^2/(3.0))
  map.radius*Point(xₛ,yₛ,zₛ)
end

function generate_Γface_to_bgface(model)
  panel_to_entity=[23,26,24,25,21,22]
  entity_to_panel=Dict(23=>1,26=>2,24=>3,25=>4,21=>5,22=>6)
  ptr_ncells_panel=zeros(Int64,length(panel_to_entity)+1)
  labels = get_face_labeling(model)
  bgface_to_mask = Gridap.Geometry.get_face_mask(labels,"boundary",2)
  Γface_to_bgface = findall(bgface_to_mask)
  face_to_entity  = labels.d_to_dface_to_entity[3][Γface_to_bgface]
  for entity in face_to_entity
    panel=entity_to_panel[entity]
    ptr_ncells_panel[panel+1]=ptr_ncells_panel[panel+1]+1
  end
  ptr_ncells_panel[1]=1
  for i=1:length(ptr_ncells_panel)-1
    ptr_ncells_panel[i+1]=ptr_ncells_panel[i+1]+ptr_ncells_panel[i]
  end
  Γface_to_bgface_panelwise=similar(Γface_to_bgface)
  for (i,entity) in enumerate(face_to_entity)
    panel=entity_to_panel[entity]
    Γface_to_bgface_panelwise[ptr_ncells_panel[panel]]=Γface_to_bgface[i]
    ptr_ncells_panel[panel]=ptr_ncells_panel[panel]+1
  end
  Γface_to_bgface_panelwise
end

function CubedSphereDiscreteModel(n,order; radius=1)
  function _cell_vector_to_dof_vector!(dof_vector,cell_node_ids, cell_vector)
    cache_cell_node_ids = array_cache(cell_node_ids)
    cache_cell_vector   = array_cache(cell_vector)
    for k=1:length(cell_node_ids)
       current_node_ids = getindex!(cache_cell_node_ids,cell_node_ids,k)
       current_values   = getindex!(cache_cell_vector,cell_vector,k)
       for (i,id) in enumerate(current_node_ids)
        dof_vector[current_node_ids[i]]=current_values[i]
       end
    end
  end

  domain = (-1,1,-1,1,-1,1)
  cells  = (n,n,n)
  model  = CartesianDiscreteModel(domain,cells)

  # Restrict model to cube surface
  Γface_to_bgface=generate_Γface_to_bgface(model)
  cube_surface_model = Gridap.Geometry.BoundaryDiscreteModel(Polytope{2},model,Γface_to_bgface)

  # Generate high-order FE map and ordering
  vector_reffe=ReferenceFE(lagrangian,VectorValue{3,Float64},order)
  V = FESpace(cube_surface_model,vector_reffe; conformity=:H1)
  vh = interpolate(MapCubeToSphere(radius),V)
  scalar_reffe=ReferenceFE(QUAD,lagrangian,Float64,order)
  xref=Gridap.ReferenceFEs.get_node_coordinates(scalar_reffe)
  xrefₖ=Fill(xref,num_cells(cube_surface_model))
  vhx=lazy_map(evaluate,Gridap.CellData.get_data(vh),xrefₖ)
  V = FESpace(cube_surface_model,scalar_reffe; conformity=:H1)
  node_coordinates = Vector{Point{3,Float64}}(undef,num_free_dofs(V))
  cell_node_ids    = get_cell_dof_ids(V)
  _cell_vector_to_dof_vector!(node_coordinates,cell_node_ids,vhx)
  cell_types  = collect(Fill(1,num_cells(cube_surface_model)))
  cell_reffes = [scalar_reffe]

  cube_surface_grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
                                                       Gridap.Arrays.Table(cell_node_ids),
                                                       cell_reffes,
                                                       cell_types,
                                                       Gridap.Geometry.Oriented())
  topology=Gridap.Geometry.get_grid_topology(cube_surface_model)
  face_labeling=Gridap.Geometry.get_face_labeling(cube_surface_model)
  Gridap.Geometry.UnstructuredDiscreteModel(cube_surface_grid,topology,face_labeling)
end

struct AnalyticalMapCubedSphereDiscreteModel{T,B,C} <: Gridap.Geometry.DiscreteModel{2,3}
  cell_map::T
  cubed_sphere_model::B
  trian::C
  function AnalyticalMapCubedSphereDiscreteModel(n;radius=1)
    domain = (-1,1,-1,1,-1,1)
    cells  = (n,n,n)
    model  = CartesianDiscreteModel(domain,cells)

    # Restrict model to cube surface
    Γface_to_bgface=generate_Γface_to_bgface(model)
    cube_surface_model = Gridap.Geometry.BoundaryDiscreteModel(Polytope{2},model,Γface_to_bgface)

    m1=Fill(Gridap.Fields.GenericField(MapCubeToSphere(radius)),num_cells(cube_surface_model))
    m2=get_cell_map(cube_surface_model)
    m=lazy_map(∘,m1,m2)

    cubed_mesh_model=CubedSphereDiscreteModel(n,1)

    # Wrap up BoundaryTriangulation
    btrian=Triangulation(cubed_mesh_model)
    trian=AnalyticalMapCubedSphereTriangulation(m,btrian)

    # Build output object
    T=typeof(m)
    B=typeof(cubed_mesh_model)
    C=typeof(trian)
    GC.gc()
    new{T,B,C}(m,cubed_mesh_model,trian)
  end
end

Gridap.Geometry.get_cell_map(model::AnalyticalMapCubedSphereDiscreteModel) = model.cell_map
Gridap.Geometry.get_grid(model::AnalyticalMapCubedSphereDiscreteModel) = Gridap.Geometry.get_grid(model.cubed_sphere_model)
Gridap.Geometry.get_grid_topology(model::AnalyticalMapCubedSphereDiscreteModel) = Gridap.Geometry.get_grid_topology(model.cubed_sphere_model)
Gridap.Geometry.get_face_labeling(model::AnalyticalMapCubedSphereDiscreteModel) = Gridap.Geometry.get_face_labeling(model.cubed_sphere_model)
Gridap.Geometry.get_triangulation(a::AnalyticalMapCubedSphereDiscreteModel) = a.trian
Gridap.Geometry.Triangulation(a::AnalyticalMapCubedSphereDiscreteModel) = a.trian

function CubedSphereDiscreteModel(n;radius=1)
  AnalyticalMapCubedSphereDiscreteModel(n;radius)
end

const CSDMT = Union{AnalyticalMapCubedSphereDiscreteModel,
                 <:Gridap.Geometry.UnstructuredDiscreteModel{2,3}}

function Gridap.CellData.get_normal_vector(model::CSDMT)
    cell_normal = Gridap.Geometry.get_facet_normal(model)
    Gridap.CellData.GenericCellField(cell_normal,Triangulation(model),ReferenceDomain())
end

function Gridap.Geometry.get_facet_normal(model::CSDMT)
  function _unit_outward_normal(v::Gridap.Fields.MultiValue{Tuple{2,3}},sign_flip::Bool)
    n1 = v[1,2]*v[2,3] - v[1,3]*v[2,2]
    n2 = v[1,3]*v[2,1] - v[1,1]*v[2,3]
    n3 = v[1,1]*v[2,2] - v[1,2]*v[2,1]
    n = VectorValue(n1,n2,n3)
    (-1)^sign_flip*n/norm(n)
  end

  # Get the Jacobian of the cubed sphere mesh
  trian = Triangulation(model)
  map   = get_cell_map(trian)
  Jt    = lazy_map(∇,map)

  # Get the index of the panel for each element
  fl = get_face_labeling(model)
  panel_id  = fl.d_to_dface_to_entity[3]
  sign_flip = [panel_id[i] == 25 || panel_id[i] == 21 || panel_id[i] == 24 for i=1:length(panel_id)]
  fsign_flip = lazy_map(Gridap.Fields.ConstantField,sign_flip)

  lazy_map(Operation(_unit_outward_normal),Jt,fsign_flip)
end


text:
function ln_like_pma(θ_system, pma::ProperMotionAnom)
    ll = 0.0

    for i in eachindex(pma.ra_epoch, pma.dec_epoch)
        pm_ra_star = 0.0
        pm_dec_star = 0.0

        # The model can support multiple planets
        for key in keys(θ_system.planets)
            θ_planet = θ_system.planets[key]

            if θ_planet.mass < 0
                return -Inf
            end

            # TODO: we are creating these from scratch for each observation instead of sharing them
            orbit = construct_elements(θ_system, θ_planet)

            # RA and dec epochs are usually slightly different
            # Note the unit conversion here from jupiter masses to solar masses to
            # make it the same unit as the stellar mass (element.mu)
            pm_ra_star += propmotionanom(orbit, pma.ra_epoch[i], θ_planet.mass*mjup2msol)[1]
            pm_dec_star += propmotionanom(orbit, pma.dec_epoch[i], θ_planet.mass*mjup2msol)[2]
        end

        residx = pm_ra_star - pma.pm_ra[i]
        residy = pm_dec_star - pma.pm_dec[i]
        σ²x = pma.σ_pm_ra[i]^2
        σ²y = pma.σ_pm_dec[i]^2
        χ²x = -0.5residx^2 / σ²x - log(sqrt(2π * σ²x))
        χ²y = -0.5residy^2 / σ²y - log(sqrt(2π * σ²y))

        ll += χ²x + χ²y
    end

    return ll
end

# TODO: image modelling for multi planet systems do not consider how "removing" one planet
# might increase the contrast of another.
# function ln_like_images(θ_system, system)
#     ll = 0.0
#     for key in keys(θ_system.planets)
#         θ_planet = θ_system.planets[key]
#         elements = construct_elements(θ_system, θ_planet)

#         if (elements.a <= 0 ||
#             elements.e < 0 ||
#             elements.plx < 0 ||
#             elements.μ <= 0)
#             ll += NaN
#             continue
#         end

#         ll += ln_like_images_element(elements, θ_planet, system)
#     end

#     # # Connect the flux at each epoch to an overall flux in this band for this planet
#     # # fᵢ = θ_band.epochs
#     # # ll += -1/2 * sum(
#     # #     (fᵢ .- θ_band.f).^2
#     # # ) / (θ_band.σ_f² * mean(fᵢ)^2)

#     # And connect that flux to a modelled Teff and mass
#     # f_model = model_interpolator(θ_planet.Teff, θ_planet.mass)
#     # ll += -1/2 * (f_model - θ_band)^2 /  (θ_planet.σ_f_model² * f_model^2)

#     return ll
# end

function ln_like_images(elements::DirectOrbits.AbstractElements, θ_planet, system)
    images = system.images
    T = eltype(θ_planet)
    ll = zero(T)
    for i in eachindex(images.epoch)

        # Calculate position at this epoch
        ra, dec = kep2cart(elements, images.epoch[i])
        # x must be negated to go from sky coordinates (ra increasing to left) to image coordinates (ra increasing to right).
        x = -ra
        y = dec

        # Get the photometry in this image at that location
        # Note in the following equations, subscript x (ₓ) represents the current position (both x and y)
        f̃ₓ = lookup_coord(images.image[i], (x, y), images.platescale[i])

        # Find the uncertainty in that photometry value (i.e. the contrast)
        r = √(x^2 + y^2)
        σₓ = images.contrast[i](r / images.platescale[i])

        # When we get a position that falls outside of our available
        # data (e.g. under the coronagraph) we cannot say anything
        # about the likelihood. This is equivalent to σₓ→∞ or log likelihood
        # of zero.
        if !isfinite(σₓ) || !isfinite(f̃ₓ)
            continue
        end

        band = images.band[i]

        # Verify the user has specified a prior or model for this band.
        if !hasproperty(θ_planet, band)
            error("No photometry prior for the band $band was specified, and neither was mass.")
        end
        # TODO: verify this is type stable
        f_band = getproperty(θ_planet, band)
        # Direct imaging likelihood.
        # Notes: we are assuming that the different images fed in are not correlated.
        # The general multivariate Gaussian likleihood is exp(-1/2 (x⃗-μ⃗)ᵀ𝚺⁻¹(x⃗-μ⃗)) + √((2π)ᵏ|𝚺|)
        # Because the images are uncorrelated, 𝚺 is diagonal and we can separate the equation
        # into a a product of univariate Gaussian likelihoods or sum of log-likelihoods.
        # That term for each image is given below.

        # Ruffio et al 2017, eqn (31)
        # Mawet et al 2019, eqn (8)

        σₓ² = σₓ^2
        ll += -1 / (2σₓ²) * (f_band^2 - 2f_band * f̃ₓ)
    end

    return ll
end

# Astrometry
function ln_like_astrom(elements, astrom::Astrometry)
    ll = 0.0

    for i in eachindex(astrom.epoch)
        x, y = kep2cart(elements, astrom.epoch[i])
        residx = astrom.ra[i] - x
        residy = astrom.dec[i] - y
        σ²x = astrom.σ_ra[i ]^2
        σ²y = astrom.σ_dec[i]^2
        χ²x = -(1/2)*residx^2 / σ²x - log(sqrt(2π * σ²x))
        χ²y = -(1/2)*residy^2 / σ²y - log(sqrt(2π * σ²y))
        ll += χ²x + χ²y
    end
    return ll
end

# Photometry
function ln_like_phot(photometry, θ_planet)
    ll = 0.0
    for i in eachindex(photometry.band)
        band = photometry.band[i]
        phot_param = getproperty(θ_planet, band)
        phot_meas = photometry.phot[i]
        if !isfinite(phot_param)
            return -Inf
        end
        σ_phot = photometry.σ_phot[i]
        resid = phot_param - phot_meas
        σ² = σ_phot^2
        χ² = -(1/2)*resid^2 / σ² - log(sqrt(2π * σ²))
        ll += χ²
    end
    return ll
end

# Overall log likelihood of the system given the parameters θ_system
function ln_like(θ_system, system::System)
    ll = 0.0
    # Fail fast if we have a negative stellar mass.
    # Users should endeavour to use priors on e.g. stellar mass
    # that are strictly positive.
    if hasproperty(θ_system, :μ) && θ_system.μ <= 0
        return -Inf
    end
    # Go through each planet in the model and add its contribution
    # to the ln-likelihood.
    for (θ_planet, planet) in zip(θ_system.planets, system.planets)

        if !isnothing(planet.photometry)
            ll += ln_like_phot(planet.photometry, θ_planet)
        end

        # We don't construct the elements object if there is no data requiring it.
        # This also means we can model e.g. photometry directly without specifying
        # all the orbital parameters.
        if isnothing(planet.astrometry) && isnothing(system.images)
            continue
        end

        # Like negative stellar mass, users should use priors with supports
        # that do not include these invalid values. But if we see them,
        # give zero likelihood right away instead of an inscrutable error
        # from some code expecting these invariants.
        if θ_planet.a <= 0 || θ_planet.e < 0 || θ_planet.e >= 1
            return -Inf
        end

        # Resolve the combination of system and planet parameters
        # as a KeplerianElements object. This pre-computes
        # some factors used in various calculations.
        kep_elements = construct_elements(θ_system, θ_planet)

        if !isnothing(planet.astrometry)
            ll += ln_like_astrom(kep_elements, planet.astrometry)
        end

        if !isnothing(system.images)
            ll += ln_like_images(kep_elements, θ_planet, system)
        end

    end

    # TODO: PMA is re-calculating some factors used in kep_elements.
    # Should think of a way to integrate it into the loop above
    if !isnothing(system.propermotionanom)
        ll += ln_like_pma(θ_system, system.propermotionanom)
    end

    return ll
end

    # Hierarchical parameters over multiple planets
    # if haskey(system.priors.priors, :σ_i²)
    #     # If the sampler wanders into negative variances, return early to prevent
    #     # taking square roots of negative values later on
    #     if θ.σ_i² < 0
    #         return -Inf
    #     end

    #     # hierarchical priors here
    #     sum_iᵢ = zero(θ.i)
    #     sum_iᵢθi² = zero(θ.i)
    #     for θ_planet in θ.planets
    #         sum_iᵢ += θ_planet.i
    #         sum_iᵢθi² += (θ_planet.i .- θ.i)^2
    #     end
    #     ll += -1/2 * sum_iᵢθi² / θ.σ_i²  - log(sqrt(2π * θ.σ_i²))
    # end
    # if haskey(system.priors.priors, :σ_Ω²)
    #     # If the sampler wanders into negative variances, return early to prevent
    #     # taking square roots of negative values later on
    #     if θ.σ_Ω² < 0
    #         return -Inf
    #     end

    #     # hierarchical priors here
    #     sum_Ωᵢ = zero(θ.Ω)
    #     sum_ΩᵢθΩ² = zero(θ.Ω)
    #     for θ_planet in θ.planets
    #         _, Ωᵢ, _  = get_ωΩτ(θ_system, θ_planet)
    #         sum_Ωᵢ += Ωᵢ
    #         sum_ΩᵢθΩ² += (Ωᵢ .- θ.Ω)^2
    #     end
    #     ll += -1/2 * sum_ΩᵢθΩ² / θ.σ_Ω²  - log(sqrt(2π * θ.σ_Ω²))
    # end

# This is a straight forward implementation that unfortunately is not type stable.
# This is because we are looping over a heterogeneous container
# function make_ln_prior(priors)
#     return function ln_prior(params)
#         lp = zero(first(params))
#         for i in eachindex(params)
#             pd = priors[i]
#             param = params[i]
#             lp += logpdf(pd, param)
#         end
#         return lp
#     end
# end

function make_ln_prior(system::System)

    # This function uses meta-programming to unroll all the code at compile time.
    # This is a big performance win, since it avoids looping over all the different
    # types of distributions that might be specified as priors.
    # Otherwise we would have to loop through an abstract vector and do runtime dispatch!
    # This way all the code gets inlined into a single tight numberical function in most cases.

    i = 0
    prior_evaluations = Expr[]

    # System priors
    for prior_distribution in values(system.priors.priors)
        i += 1
        ex = :(
            lp += $logpdf($prior_distribution, arr[$i])
        )
        push!(prior_evaluations,ex)
    end

    # Planet priors
    for planet in system.planets
        # for prior_distribution in values(planet.priors.priors)
        for (key, prior_distribution) in zip(keys(planet.priors.priors), values(planet.priors.priors))
            i += 1
            # Work around for Beta distributions.
            # Outside of the range [0,1] logpdf returns -Inf.
            # This works fine, but AutoDiff outside this range causes a DomainError.
            if typeof(prior_distribution) <: Beta
                ex = :(
                    lp += 0 <= arr[$i] < 1 ? $logpdf($prior_distribution, arr[$i]) : -Inf
                )
            else
                ex = :(
                    lp += $logpdf($prior_distribution, arr[$i])
                )
            end
            push!(prior_evaluations,ex)
        end
    end

    # Here is the function we return.
    # It maps an array of parameters into our nested named tuple structure
    # Note: eval() would normally work fine here, but sometimes we can hit "world age problemms"
    # The RuntimeGeneratedFunctions package avoids these in all cases.
    return @RuntimeGeneratedFunction(:(function (arr)
        l = $i
        @boundscheck if length(arr) != l
            error("Expected exactly $l elements in array (got $(length(arr)))")
        end
        lp = zero(first(arr))
        # Add contributions from planet priors
        @inbounds begin
           $(prior_evaluations...)
        end
        return lp
    end))
end

# Same as above, but assumes the input to the log prior was sampled
# using transformed distributions from Bijectors.jl
# Uses logpdf_with_trans() instead of logpdf to make the necessary corrections.
function make_ln_prior_transformed(system::System)

    i = 0
    prior_evaluations = Expr[]

    # System priors
    for prior_distribution in values(system.priors.priors)
        i += 1
        ex = :(
            lp += $logpdf_with_trans($prior_distribution, arr[$i], true)
        )
        push!(prior_evaluations,ex)
    end

    # Planet priors
    for planet in system.planets
        # for prior_distribution in values(planet.priors.priors)
        for (key, prior_distribution) in zip(keys(planet.priors.priors), values(planet.priors.priors))
            i += 1
            ex = :(
                lp += $logpdf_with_trans($prior_distribution, arr[$i], true)
            )
            push!(prior_evaluations,ex)
        end
    end

    # Here is the function we return.
    # It maps an array of parameters into our nested named tuple structure
    # Note: eval() would normally work fine here, but sometimes we can hit "world age problemms"
    # The RuntimeGeneratedFunctions package avoids these in all cases.
    return @RuntimeGeneratedFunction(:(function (arr)
        l = $i
        @boundscheck if length(arr) != l
            error("Expected exactly $l elements in array (got $(length(arr)))")
        end
        lp = zero(first(arr))
        # Add unrolled prior evaluations
        @inbounds begin
           $(prior_evaluations...)
        end
        return lp
    end))
end

# Replaces `θ = Bijectors.invlink.(priors_vec, θ_t)` with a type stable
# unrolled version.
function make_Bijector_invlinkvec(priors_vec)

    i = 0
    parameter_transformations = Expr[]

    # System priors
    for prior_distribution in priors_vec
        i += 1
        ex = :(
            theta_out[$i] = $(Bijectors.invlink)($prior_distribution, arr[$i])
        )
        push!(parameter_transformations, ex)
    end

    # Here is the function we return.
    # It maps an array of parameters into our nested named tuple structure
    # Note: eval() would normally work fine here, but sometimes we can hit "world age problemms"
    # The RuntimeGeneratedFunctions package avoids these in all cases.
    return @RuntimeGeneratedFunction(:(function (arr)
        l = $i
        theta_out = @MVector zeros(eltype(arr), l)
        # theta_out = zeros(eltype(arr), l)
        @boundscheck if length(arr) != l
            error("Expected exactly $l elements in array (got $(length(arr)))")
        end
        # Add unrolled parameter transformations to fill theta_out
        @inbounds begin
           $(parameter_transformations...)
        end
        return theta_out
    end))
end

function ln_post(θ, system::System)
    return ln_prior(θ, system) + ln_like(θ, system)
end


text:
# Syntax patterns
function prefix(io, call, argument1)
    print(io, call)
    print(io, ' ')
    print(io, argument1)
end

function parenthesize_prefix(io, call, argument1)
    print(io, call)
    print(io, ' ')
    print(io, '(')
    print(io, argument1)
    print(io, ')')
end

function postfix(io, call, argument1)
    print(io, argument1)
    print(io, ' ')
    print(io, call)
end

function tight_infix(io, call, argument1, argument2)
    print(io, argument1)
    print(io, call)
    print(io, argument2)
end

function parenthesize_infix(io, call, argument1, argument2)
    print(io, argument1)
    print(io, ' ')
    print(io, call)
    print(io, ' ')
    print(io, argument2)
end

function infix(io, call, argument1, argument2, arguments...)
    print(io, argument1)
    print(io, ' ')
    print(io, call)
    print(io, ' ')
    infix(io, call, argument2, arguments...)
end

function infix(io, call, argument1)
    print(io, argument1)
end

function head_call_tail(io, call, argument1, arguments...)
    print(io, argument1)
    print(io, ' ')
    print(io, call)
    print(io, ' ')
    join(io, arguments, ", ")
end

function call_tail_head(io, call, argument1, arguments...)
    print(io, call)
    print(io, ' ')
    join(io, arguments, ", ")
    print(io, ' ')
    print(io, argument1)
end

function when_then_else(io, condition, action, arguments...)
    print(io, "WHEN ")
    print(io, condition)
    print(io, " THEN ")
    print(io, action)
    print(io, ' ')
    when_then_else(io, arguments...)
end

function when_then_else(io, action)
    print(io, "ELSE ")
    print(io, action)
end

function case(io, call, arguments...)
    print(io, call)
    print(io, ' ')
    when_then_else(io, arguments...)
end

# Use the call to determine the syntax pattern
function show(io::IO, sql_expression::SQLExpression)
    call = sql_expression.call
    arguments = sql_expression.arguments

    if call === :.
        tight_infix(io, call, arguments...)
    elseif call === :(=)
        infix(io, call, arguments...)
    elseif call === :<
        infix(io, call, arguments...)
    elseif call === Symbol("<>")
        infix(io, call, arguments...)
    elseif call === Symbol("||")
        infix(io, call, arguments...)
    elseif call === :*
        infix(io, call, arguments...)
    elseif call === :+
        infix(io, call, arguments...)
    elseif call === :%
        infix(io, call, arguments...)
    elseif call === :AND
        infix(io, call, arguments...)
    elseif call === :AS
        parenthesize_infix(io, call, arguments...)
    elseif call === :CASE
        case(io, call, arguments...)
    elseif call === :DESC
        postfix(io, call, arguments...)
    elseif call === :FROM
        parenthesize_prefix(io, call, arguments...)
    elseif call === Symbol("GROUP BY")
        head_call_tail(io, call, arguments...)
    elseif call === :IN
        infix(io, call, arguments...)
    elseif call === Symbol("INNER JOIN")
        infix(io, call, arguments...)
    elseif call === Symbol("IS NOT DISTINCT FROM")
        infix(io, call, arguments...)
    elseif call === Symbol("IS NULL")
        postfix(io, call, arguments...)
    elseif call === :LIKE
        infix(io, call, arguments...)
    elseif call === :LIMIT
        infix(io, call, arguments...)
    elseif call === :OFFSET
        infix(io, call, arguments...)
    elseif call === :ON
        infix(io, call, arguments...)
    elseif call === :OR
        infix(io, call, arguments...)
    elseif call === Symbol("ORDER BY")
        head_call_tail(io, call, arguments...)
    elseif call === :NOT
        prefix(io, call, arguments...)
    elseif call === :SELECT
        call_tail_head(io, call, arguments...)
    elseif call === Symbol("SELECT DISTINCT")
        call_tail_head(io, call, arguments...)
    elseif call === :WHERE
        infix(io, call, arguments...)
    else
        print(io, call)
        print(io, '(')
        join(io, arguments, ", ")
        print(io, ')')
    end
end

# Default to * if no columns are specifcied
function finalize(sql_expression::SQLExpression)
    if sql_expression.call in (:FROM, Symbol("GROUP BY"), Symbol("INNER JOIN"))
        SQLExpression(:SELECT, sql_expression, :*)
    else
        sql_expression
    end
end

function finalize(something)
    something
end


text:
module SymDiff

using LinearAlgebra
using Random
using Test

using Parameters: @with_kw, @unpack
using Setfield: set
using StaticArrays: SVector

using SymEngine: @vars, diff, subs, symbols
import SymEngine

using Bifurcations.ArrayUtils: cast_container, _normalize!
using Bifurcations.BifurcationsBase: contkind, HopfCont, SaddleNodeCont
using Bifurcations.Codim2: DiffEqCodim2Problem, get_augsys_cache, augsys,
    NormalizingAS, BackReferencingAS, ds_state, ds_eigvec
using Bifurcations.Continuations: get_prob_cache, get_u0, residual
import Bifurcations.Continuations: residual_jacobian

∙(a, b) = sum(a .* b)

sjac(f, x) = [diff(f[i], x[j]) for i in 1:length(f), j in 1:length(x)]

@with_kw struct SymDiffCodim2Problem
    prob::DiffEqCodim2Problem
    param1::SymEngine.Basic
    param2::SymEngine.Basic
    x::Vector
    v::Vector
    v0::Union{Vector{SymEngine.Basic}, Nothing}
    w::Union{SymEngine.Basic, Nothing}
    f::Vector
    H::Vector
    dfdx::Matrix
    dHdu::Matrix
end

interleave(x, y) = transpose([x y])[:]

function SymDiffCodim2Problem(prob::DiffEqCodim2Problem)
    @vars param1 param2
    p = prob.de_prob.p
    p = set(p, prob.param_axis1, param1)
    p = set(p, prob.param_axis2, param2)

    n = length(prob.de_prob.u0)
    x = [symbols("x_$i") for i in 1:n]
    # x = cast_container(typeof(prob.de_prob.u0), x)
    x = SVector(x...)
    f = prob.de_prob.f(x, p, prob.de_prob.tspan[1])
    dfdx = sjac(f, x)

    if contkind(prob) isa SaddleNodeCont
        w = nothing
        v = [symbols("v_$i") for i in 1:n]
        if prob.augmented_system isa NormalizingAS
            v0 = v
        else
            v0 = [symbols("v0_$i") for i in 1:n]
        end
        u = vcat(x, v, param1, param2)
        H = vcat(
            f,
            dfdx * v,
            v0 ∙ v - 1,
        )
        if prob.augmented_system isa NormalizingAS
            v0 = nothing
        end
    else
        @assert contkind(prob) isa HopfCont
        @vars w
        vr = [symbols("vr_$i") for i in 1:n]
        vi = [symbols("vi_$i") for i in 1:n]
        v = interleave(vr, vi)
        u = vcat(x, v, w, param1, param2)
        H_head = vcat(
            f,
            interleave(
                dfdx * vr .+ w .* vi,
                dfdx * vi .- w .* vr,
            ),
        )
        if prob.augmented_system isa NormalizingAS
            H = vcat(
                H_head,
                vr ∙ vr + vi ∙ vi - 1,
                vr ∙ vi,
            )
            v0 = nothing
        else
            v0r = [symbols("v0r_$i") for i in 1:n]
            v0i = [symbols("v0i_$i") for i in 1:n]
            v0 = interleave(v0r, v0i)
            H = vcat(
                H_head,
                v0r ∙ vr + v0i ∙ vi - 1,
                v0r ∙ vi - v0i ∙ vr,
            )
        end
    end

    dHdu = sjac(H, u)  # aka J

    return SymDiffCodim2Problem(
        prob,
        param1,
        param2,
        x,
        v,
        v0,
        w,
        f,
        H,
        dfdx,
        dHdu,
    )
end

@with_kw struct SymDiffCodim2Cache{AC}
    prob::SymDiffCodim2Problem
    augsys_cache::AC
end

SymDiffCodim2Cache(prob::DiffEqCodim2Problem) =
    SymDiffCodim2Cache(SymDiffCodim2Problem(prob))
SymDiffCodim2Cache(prob) = SymDiffCodim2Cache(prob, get_augsys_cache(prob.prob))

function all_symbols(prob::SymDiffCodim2Problem)
    @unpack param1, param2, x, v, v0, w = prob
    syms = Any[x, v]
    v0 !== nothing && push!(syms, v0)
    w  !== nothing && push!(syms, w)
    return vcat(syms..., param1, param2)
end

all_symbols(prob_cache::SymDiffCodim2Cache) = all_symbols(prob_cache.prob)

function all_values(u, prob_cache::SymDiffCodim2Cache)
    param1, param2 = u[end-1:end]
    x = ds_state(prob_cache.prob.prob, u)
    v = ds_eigvec(prob_cache.prob.prob, u)
    if contkind(prob_cache.prob.prob) isa SaddleNodeCont
        vals = Any[x, v]
    else
        @assert contkind(prob_cache.prob.prob) isa HopfCont
        vals = Any[x, interleave(real(v), imag(v))]
    end
    if contkind(prob_cache.prob.prob) isa HopfCont
        if augsys(prob_cache.prob.prob) isa BackReferencingAS
            push!(vals, interleave(
                real(prob_cache.augsys_cache.v),
                imag(prob_cache.augsys_cache.v),
            ))
        end
        push!(vals, u[end-2])  # w
    else
        if augsys(prob_cache.prob.prob) isa BackReferencingAS
            push!(vals, prob_cache.augsys_cache.v)
        end
    end
    return vcat(vals..., param1, param2) :: AbstractVector{<: Real}
end

function evalenv(u, prob_cache::SymDiffCodim2Cache)
    syms = all_symbols(prob_cache)
    vals = all_values(u, prob_cache)
    if length(syms) != length(vals)
        error("""Number of symbols and values are different:
        #symbols = $(length(syms))
        #values  = $(length(vals))
        symbols  = $syms
        values   = $vals
        """)
    end
    return zip(syms, vals)
end

evalsym(e, env) = Float64(subs(e, env...))
evalsym(e::AbstractArray, env) = evalsym.(e, Ref(env))

function residual_jacobian(u, prob_cache::SymDiffCodim2Cache)
    env = evalenv(u, prob_cache)
    H = evalsym(prob_cache.prob.H, env)
    J = evalsym(prob_cache.prob.dHdu, env)
    return H, J
end

function _randn(rng, ::Type{Complex{T}}, n) where T
    reals = randn(rng, T, 2n)
    cs = reinterpret(Complex{T}, reals)
    @assert length(cs) == n
    return cs
end

_randn(args...) = randn(args...)

function test_residual_jacobian(prob::DiffEqCodim2Problem;
                                num = 1,
                                rng = MersenneTwister(0))
    sd_cache = SymDiffCodim2Cache(prob)
    prob_cache = get_prob_cache(prob)
    u0 = get_u0(prob)
    for _ in 1:num
        u = typeof(u0)(randn(rng, eltype(u0), length(u0)))
        if prob.augmented_system isa BackReferencingAS
            let v = sd_cache.augsys_cache.v
                v = typeof(v)(_randn(rng, eltype(v), length(v)))
                v = _normalize!(v)
                sd_cache.augsys_cache.v = v
                prob_cache.augsys_cache.v = v
            end
        end

        H1_actual = residual(u, prob_cache)
        H2_actual, J_actual = residual_jacobian(u, prob_cache)
        H_desired, J_desired = residual_jacobian(u, sd_cache)

        @test H1_actual ≈ H_desired
        @test H2_actual ≈ H_desired
        @test J_actual ≈ J_desired
    end
end

end  # module


text:
# A bunch of MCMC probability density function samplers
#
# Refs:
# https://darrenjw.wordpress.com/2010/08/15/metropolis-hastings-mcmc-algorithms/
# https://theclevermachine.wordpress.com/2012/11/04/mcmc-multivariate-distributions-block-wise-component-wise-updates/
#
# https://pymc-devs.github.io/pymc/theory.html

"""
    metropolis(pdf, sample_ppdf, theta0;
               niter=10^5,
               nburnin=niter÷2,
               nthin=1,
               use_progress_meter=true,
               hasblob=false, # whether the pdf also returns a blob
               init_blobs=init_output_vector, # initialize storage for blobs, of form (blob0, nsamples) -> container
               reduce_blob!=(blobs, blob) -> push!(blobs, blob) # add or reduce one blob into blob-storage

Metropolis sampler, the simplest MCMC method. Can only be used with symmetric
proposal distributions (`ppdf`).

Input:

- pdf -- The probability log-density function to sample.

         It can also returns an arbitrary blob of
         something as second output argument.
         `p, blob = pdf(theta)` where `theta` are the parameters.

- sample_ppdf -- draws a sample for the proposal/jump distribution
                 `sample_ppdf(theta)`.  Needs to be symmetric:
                 `ppdf(theta1)==ppdf(theta2)`.  Example:
                 `theta -> c*randn() + theta`
- theta0 -- initial value of parameters (<:AbstractVector)

Optional keyword input:

- niter -- number of steps to take (10^5)
- nburnin -- number of initial steps discarded, aka burn-in (niter/3)
- nthin -- only store every n-th sample (default=1)
- use_progress_meter=true : whether to show a progress meter
- hasblob=false -- whether the pdf also returns a blob
- init_blobs=init_output_vector -- initialize storage for blobs, of form (blob0, nsamples) -> container
- reduce_blob!=(blobs, blob) -> push!(blobs, blob) -- add or reduce one blob into blob-storage

Notes:
- single threaded

Output:

- samples:
  - if typeof(theta0)==Vector then Array{eltype(theta0)}(length(theta0), niter-nburnin)
  - if typeof(theta0)!=Vector then Array{typeof(theta0)}(niter-nburnin)
- accept_ratio: ratio accepted to total steps
- logdensities: the value of the log-density for each sample
- blobs: anything else that the pdf-function returns as second argument
"""
function metropolis(pdf, sample_ppdf, theta0;
                    niter=10^5,
                    nburnin=niter÷2,
                    nthin=1,
                    use_progress_meter=true,
                    hasblob=false, # whether the pdf also returns a blob
                    init_blobs=init_output_vector, # initialize storage for blobs, of form (blob0, nsamples) -> container
                    reduce_blob! =(blobs, blob) -> push!(blobs, blob) # add one blob to blob-storage
                    )
    # initialize
    theta0 = deepcopy(theta0)
    pdf_ = hasblob ? pdf : t -> (pdf(t), nothing)
    p0, blob0 = pdf_(theta0)
    prog = use_progress_meter ? Progress(length((1-nburnin):(niter-nburnin))÷nthin, 1, "Metropolis, niter=$niter: ", 25) : nothing

    # run
    #@inferred _metropolis(pdf_, sample_ppdf, theta0, p0, blob0, niter, nburnin, nthin, prog, hasblob, init_blobs, reduce_blob!)
    _metropolis(pdf_, sample_ppdf, theta0, p0, blob0, niter, nburnin, nthin, prog, hasblob, init_blobs, reduce_blob!)
end

"Makes output arrays"
function init_output_vector(v0s, nsamples)
    vs = typeof(v0s)[];
    sizehint!.(vs, nsamples)
    return vs
end
init_output_vector(v0s::Nothing, nsamples) = nothing

function _metropolis(pdf, sample_ppdf, theta0, p0, blob0, niter, nburnin, nthin, prog, hasblob, init_blobs, reduce_blob!)
    # number of stored samples
    nsamples = (niter-nburnin) ÷ nthin

    # initialize output arrays
    thetas = init_output_vector(theta0, nsamples)
    blobs = init_blobs(blob0, nsamples)
    logdensities = init_output_vector(p0, nsamples)

    naccept = 0
    nn = 1
    for n=(1-nburnin):(niter-nburnin)
        # take a step:
        theta1 = sample_ppdf(theta0)
        p1, blob1 = pdf(theta1)
        # if p1/p0>rand() then accept:
        if  p1-p0 > log(rand())
            theta0 = theta1
            blob0 = blob1
            p0 = p1
            naccept += 1
        end
        # storage
        if rem(n, nthin)==0
            prog!=nothing && ProgressMeter.next!(prog; showvalues =
                                                 [(:accept_ratio, round(naccept/nn, sigdigits=3)),
                                                  (:burnin_phase, n<=0)])
            if n>0 # after burnin
                push!(thetas, theta0)
                hasblob && reduce_blob!(blobs, blob0)
                push!(logdensities, p0)
            end
        end
        nn +=1 # number of steps
        if n==0 # reset counters at end of burnin phase
            naccept=0
            nn=1
        end
    end
    accept_ratio = naccept/(niter-nburnin)
    return thetas, accept_ratio, logdensities, blobs
end

"""
    emcee(pdf, theta0s;
          niter=10^5,
          nburnin=niter÷2,
          nthin=1,
          a_scale=2.0, # step scale parameter.  Probably needn't be adjusted
          use_progress_meter=true,
          hasblob=false,
          init_blobs=init_output_vector, # initialize storage for blobs, of form (blob0, nsamples_walker) -> container
          reduce_blob!=(blobs, blob) -> push!(blobs, blob) # add one blob to blob-storage

The affine invariant MCMC sampler, aka MC hammer in its Python emcee implementation

Input:

- pdf -- The probability log-density function to sample.

         The pdf can also returns an arbitrary blob of something.
         `p, blob = pdf(theta)` where `theta` are the parameters.
- theta0s -- initial parameters, a Vector of length of the number of walkers.
             Consider creating them with `make_theta0s(theta0, ballradius)`

Note that theta0s[1] need to support `.+`, `.*`, and `length`.

Optional key-word input:
- niter -- total number of steps to take (10^5) (==total number of log-density evaluations).
- nburnin -- total number of initial steps discarded, aka burn-in (niter/3)
- nthin -- only store every n-th sample (default=1)
- use_progress_meter=true : whether to show a progress meter
- hasblob=false -- whether the pdf also returns a blob
- init_blobs=init_output_vector -- initialize storage for blobs of one walker, of form `(blob0, nsamples) -> container`
- reduce_blob!=(blobs, blob) -> push!(blobs, blob) -- add or reduce one blob into blob-storage (of one walker)

Output:
- samples: of type & size `Matrix{typeof(theta0)}(niter-nburnin, nwalkers)`
- accept_ratio: ratio of accepted to total steps
- logdensities: the value of the log-density for each sample
- blobs: anything else that the pdf-function returns as second argument

Notes:
- use `squash_walkers(samples)` to concatenate all walkers-chains into one.

Reference:
- Goodman, Jonathan, and Jonathan Weare, 2010
  http://dx.doi.org/10.2140/camcos.2010.5.65
- emcee https://github.com/dfm/emcee
- Foreman-Mackey et al. 2013, emcee: The MCMC hammer,
  https://github.com/dfm/emcee

Example

    thetas, accept_ratio, logdensities = emcee(x->-sum(x.^2), [1.0, 2, 5, 8])
"""
function emcee(pdf, theta0s;
               niter=10^5,
               nburnin=niter÷2,
               nthin=1,
               a_scale=2.0, # step scale parameter.  Probably needn't be adjusted
               use_progress_meter=true,
               hasblob=false,
               init_blobs=init_output_vector, # initialize storage for blobs, of form (blob0, nsamples_walker) -> container
               reduce_blob! =(blobs, blob) -> push!(blobs, blob) # add one blob to blob-storage
               )
    theta0s = deepcopy(theta0s)

    @assert a_scale>1
    nwalkers = length(theta0s)
    @assert iseven(nwalkers) "Use an even number of walkers."
    niter_walker = niter ÷ nwalkers
    nburnin_walker = nburnin ÷ nwalkers
    @assert nwalkers>=length(theta0s[1])+2 "Use more walkers: at least DOF+2, but better many more."

    # initialize
    pdf_ = hasblob ? pdf : t -> (pdf(t), nothing)
    tmp = pdf_.(theta0s)
    p0s, blob0s = getindex.(tmp, 1), getindex.(tmp, 2)

    # initialize progress meter
    prog = use_progress_meter ? Progress(length((1-nburnin_walker):(niter_walker-nburnin_walker))÷nthin, 1, "emcee, niter=$niter, nwalkers=$nwalkers: ", 25) : nothing
    # do the MCMC
    _emcee(pdf_, theta0s, p0s, blob0s, niter_walker, nburnin_walker, nwalkers, nthin, a_scale, prog, hasblob, init_blobs, reduce_blob!)
end

"Makes output arrays"
init_output_emcee(v0s, nwalkers, nsamples_walker, init_fn=init_output_vector) =
    [init_fn(v0s[i], nsamples_walker) for i=1:nwalkers]
init_output_emcee(v0s::Vector{Nothing}, nwalkers, nsamples_walker, init_fn=init_output_vector) = nothing

"g-pdf, see eq. 10 of Foreman-Mackey et al. 2013."
g_pdf(z, a) = 1/a<=z<=a ? 1/sqrt(z) * 1/(2*(sqrt(a)-sqrt(1/a))) : zero(z)

"Inverse cdf of g-pdf, see eq. 10 of Foreman-Mackey et al. 2013."
cdf_g_inv(u, a) = (u*(sqrt(a)-sqrt(1/a)) + sqrt(1/a) )^2

"Sample from g using inverse transform sampling.  a=2.0 is recommended."
sample_g(a) = cdf_g_inv(rand(), a)

function _emcee(pdf, theta0s, p0s, blob0s, niter_walker, nburnin_walker, nwalkers, nthin, a_scale, prog, hasblob, init_blobs, reduce_blob!)
    # number of stored samples of each walker
    nsamples_walker = (niter_walker-nburnin_walker) ÷ nthin

    # initialize output
    thetas = init_output_emcee(theta0s, nwalkers, nsamples_walker)
    blobs = init_output_emcee(blob0s, nwalkers, nsamples_walker, init_blobs)
    logdensities = init_output_emcee(p0s, nwalkers, nsamples_walker)

    # initialization and work arrays:
    naccept = zeros(Int, nwalkers)
    N = length(theta0s[1])
    nn = 1
    for n in (1-nburnin_walker):(niter_walker-nburnin_walker)
        for batch = 1:2
            ncs, ncos = circshift(SVector(1:nwalkers÷2, nwalkers÷2+1:nwalkers), batch-1)
            Threads.@threads for nc in ncs
                # draw a random other walker
                no = rand(ncos)
                # sample g (eq. 10)
                z = sample_g(a_scale)

                # sample a new theta with the stretch move (eq. 7):
                theta1 = theta0s[no] .+ z .* (theta0s[nc] .- theta0s[no])
                # and calculate the theta1 density:
                p1, blob1 = pdf(theta1)

                # if z^(N-1)*p1/p0>rand() then accept:
                if (N-1)*log(z) + p1 - p0s[nc] >= log(rand())
                    theta0s[nc] = theta1
                    p0s[nc] = p1

                    blob0s[nc] = blob1
                    naccept[nc] += 1
                end

                if n>0 && rem(n,nthin)==0 # store theta after burn-in
                    push!(thetas[nc], theta0s[nc])
                    hasblob && reduce_blob!(blobs[nc], blob0s[nc])
                    push!(logdensities[nc], p0s[nc])
                end
            end # for nc =1:ncs
        end # batch = 1:2
        if rem(n,nthin)==0 # store theta after burn-in
            macc = mean(naccept)
            sacc = sqrt(var(naccept)) # std(naccept) is not type-stable!
            outl = sum(abs.(naccept.-macc).>2*sacc)
            prog!=nothing && ProgressMeter.next!(prog;
                                                 showvalues = [(:accept_ratio_mean, round(macc/nn,sigdigits=3)),
                                                               (:accept_ratio_std, round(sacc/nn,sigdigits=3)),
                                                               (:accept_ratio_outliers, outl),
                                                               (:burnin_phase, n<=0)])
        end
        if n==0 # reset after burnin
            fill!(naccept,0)
            nn = 1
        end
        nn +=1
    end # for n=(1-nburnin_walker):(niter_walker-nburnin_walker)
    accept_ratio = [na/(niter_walker-nburnin_walker) for na in naccept]
    return thetas, accept_ratio, logdensities, blobs
end

"""
    make_theta0s(theta0::T, ball_radius::T, pdf, nwalkers;
                      ball_radius_halfing_steps=7, # make the ball smaller this often
                      ntries=100*nwalkers, # how many tries per pall radius
                      hasblob=false) where T

Tries to find theta0s with nonzero density within a ball from theta0.  To be used with
`emcee`.

Example

    pdf = x-> -sum(x.^2)
    nwalkers = 100
    samples = emcee(pdf, make_theta0s([0.0, 0.0], [0.1, 0.1], pdf, nwalkers), nburnin=0, use_progress_meter=false);
"""
function make_theta0s(theta0::T, ball_radius, pdf, nwalkers;
                      ball_radius_halfing_steps=7,
                      ntries=100,
                      hasblob=false) where T
    npara = length(theta0)
    if ball_radius isa Number && !(T<:Number)
        ball_radius = ones(npara) * ball_radius
    end
    @assert length(ball_radius)==npara

    theta0s = T[]

    for i=1:nwalkers
        for k=1:ball_radius_halfing_steps
            j = 0
            ball_radius *= 1/2^(k-1)
            for j=1:ntries
                tmp = if npara==1
                    theta0 .+ randn().*ball_radius
                else
                    theta0 .+ randn(npara).*ball_radius
                end
                if hasblob
                    p0, blob0 = pdf(tmp)
                else
                    p0, blob0 = pdf(tmp), nothing
                end
                if p0>-Inf
                    push!(theta0s, tmp)
                    break
                end
            end
            length(theta0s)==i && break # found suitable thetas for this walker
            j==ntries && k==ball_radius_halfing_steps &&
                error("Could not find suitable initial theta.  PDF is zero in too many places inside ball.")
        end
    end
    return theta0s
end

"""
    squash_walkers(thetas, accept_ratio, blobs, logdensities;
                                                         drop_low_accept_ratio=false,
                                                         drop_fact=2,
                                                         verbose=true,
                                                         order=false, # if true the samples are ordered, this may not work with blobs not stored in vectors
                                                         merge_blobs!=append!) # function to merge one blobs into another blobs (blobs1, blobs1) -> merge blobs2 into blobs1

Puts the samples of all emcee walkers into one vector.

This can also drop walkers which have a low accept ratio (this can happen with
emcee, but I am anything but sure whether this is ok),
by setting drop_low_accept_ratio=true.  drop_fact -> increase to drop more walkers.

Returns:
- thetas
- mean(accept_ratio[walkers2keep])
- log-densities
- blobs
"""
function squash_walkers(thetas, accept_ratio, logdensities=nothing, blobs=nothing;
                        drop_low_accept_ratio=false,
                        drop_fact=2,
                        verbose=true,
                        order=false, # if true the samples are ordered, this may not work with blobs not stored in vectors
                        merge_blobs! =append!) # function to merge one blobs into another blobs (blobs1,blobs2) ->

    nwalkers=length(accept_ratio)
    if drop_low_accept_ratio
        walkers2keep = Int[]
        # this is a bit of a hack, but emcee seems to sometimes have
        # walkers which are stuck in oblivion.  These have a very low
        # acceptance ratio, thus filter them out.
        ma,sa = median(accept_ratio), std(accept_ratio)
        verbose && println("Median accept ratio is $ma, standard deviation is $sa\n")
        for nc=1:nwalkers
            if accept_ratio[nc]<=ma-drop_fact*sa # this 1 is heuristic
                verbose && println("Dropping walker $nc with low accept ratio $(accept_ratio[nc])")
                continue
            end
            push!(walkers2keep, nc)
        end
    else
        walkers2keep = collect(1:nwalkers)
    end

    t = copy(thetas[walkers2keep[1]])
    append!.(Ref(t), thetas[walkers2keep[2:end]])

    if logdensities==nothing
        l = nothing
    else
        l = copy(logdensities[walkers2keep[1]])
        append!.(Ref(l), logdensities[walkers2keep[2:end]])
    end

    if blobs==nothing
        b = nothing
    else
        b = deepcopy(blobs[walkers2keep[1]])
        merge_blobs!.(Ref(b), blobs[walkers2keep[2:end]])
    end

    if order # order such that early samples are early in thetas
        nc = length(walkers2keep)
        ns = length(thetas[1])
        perm = sortperm(vcat([collect(1:ns) for i=1:nc]...))
        if b!=nothing
            b = b[perm]
        end
        if l!=nothing
            l = l[perm]
        end
        t = t[perm]
    end
    return t, mean(accept_ratio[walkers2keep]), l, b
end

## TODO
# """
#          retrace_samples(pdf, thetas_in;
#                          logpdf=true,
#                          blob_reduce! = default_blob_reduce!,
#                          use_progress_meter=true)

# This function allows to re-run the pdf for some, given thetas (`thetas_in`).  This is probably
# a bit niche, but hey.

# This I use when only sampling the prior but then want to get some function-blob evaluations.
# That way the evaluation of an expensive forward function can be avoided.
# """
# function retrace_samples(pdf, thetas_in;
#                          logpdf=true,
#                          blob_reduce! = default_blob_reduce!,
#                          use_progress_meter=true,
#                          )
#     # initialize
#     nwalkers = 0
#     niter = size(thetas_in,2)
#     nburnin = 0
#     nthin = 1
#     pdf_, p0, theta0, blob0, thetas, blobs, nwalkers, pdftype, logposts =
#         _initialize(pdf, thetas_in[:,1], niter, nburnin, logpdf, nwalkers, nthin, blob_reduce!, make_SharedArray=false)

#     prog = use_progress_meter ? Progress(length((1-nburnin):(niter-nburnin))÷nthin, 1, "Retracing samples, niter=$niter: ", 25) : nothing

#     # run
#     _retrace_samples!(p0, thetas_in, blob0, thetas, blobs, pdf_, niter,
#                       nburnin, nthin, pdftype, blob_reduce!, prog, logposts)
# end
# function _retrace_samples!(p0, thetas_in, blob0, thetas, blobs, pdf, niter,
#                            nburnin, nthin, pdftype, blob_reduce!, prog, logposts)
#     @inbounds for ni=1:size(thetas_in,2)
#         theta = thetas_in[:,ni]
#         # take a step:
#         p, blob = pdf(theta)

#         _setindex!(thetas, theta, ni)
#         blob_reduce!(blobs, blob, ni)
#         logposts[ni] = p

#         prog!=nothing && ProgressMeter.next!(prog)
#     end
#     accept_ratio = -1
#     return thetas, accept_ratio, blobs, logposts
# end


text:
module MultiTankMaterialBalance

export read_rates, read_params, process_params!, save_rates!, save_params!
export getvalue, getvalues, getparams!, getparams!!, setparams!, solve!, update_targ!, calc_well_index!
export NonlinearProblem, LinearScaling, SigmoidScaling
export DenseLinearSolver, RecursiveLinearSolver, SparseLinearSolver
export NewtonSolver, FittingSet, TargetFunction, AdjointSolver

using CSV
using DataFrames
using DataFramesMeta
using StringEncodings
using Dates
using Underscores
using LinearAlgebra
using SparseArrays
using UnPack
using Printf
using SuiteSparse
using RecursiveFactorization
using AMD
using Metis
using LoopVectorization

include("types.jl")
include("utils.jl")
include("problem.jl")
include("linalg.jl")
include("solve.jl")
include("newton.jl")
include("scale.jl")
include("parameters.jl")
include("target.jl")
include("gradient.jl")
include("adjoint.jl")
include("data.jl")

end


text:
module ChemistryFeaturization

using SimpleWeightedGraphs

encodable_elements(a::Any) = throw(MethodError(encodable_elements, a))
encode(a::Any, object_to_be_encoded) = throw(MethodError(encode, a))
decode(a::Any, encoded_features) = throw(MethodError(decode, a))
elements(a::Any) = throw(MethodError(elements, a))
output_shape(a::Any) = throw(MethodError(output_shape, a))

include("data.jl")
export Data

include("abstracts/abstracttypes.jl")
export AbstractType

include("codecs/codecs.jl")
export Codec

include("utils/Utils.jl")
export Utils

include("atoms/atoms.jl")
export Atoms

using .Atoms: AtomGraph, visualize, elements
export AtomGraph, visualize

include("features/features.jl")
export FeatureDescriptor
using .FeatureDescriptor
export ElementFeatureDescriptor, SpeciesFeatureDescriptor, BondFeatureDescriptor
export output_shape, get_value

include("featurizations/featurizations.jl")
export Featurization
using .Featurization: GraphNodeFeaturization, encode
export GraphNodeFeaturization, encode

export encodable_elements, decode

include("featurizedatoms.jl")
export FeaturizedAtoms, featurize, decode

end


text:
@testset "846.hand-of-straights.jl" begin
    @test is_n_straight_hand([1, 2, 3, 6, 2, 3, 4, 7, 8], 3) == true
    @test is_n_straight_hand([1, 2, 3, 4, 5], 4) == false
    @test is_n_straight_hand([1, 1, 2, 3], 4) == false
end


text:
using MathOptInterface
const MOI = MathOptInterface

include("scaled_psd_cone_bridge.jl")

MOI.Utilities.@product_of_sets(
    Cones,
    MOI.Zeros,
    MOI.Nonnegatives,
    MOI.SecondOrderCone,
    ScaledPSDCone,
    MOI.ExponentialCone,
    MOI.DualExponentialCone,
    MOI.PowerCone{T},
    MOI.DualPowerCone{T},
)

const OptimizerCache{T} = MOI.Utilities.GenericModel{
    Cdouble,
    MOI.Utilities.ObjectiveContainer{Cdouble},
    MOI.Utilities.VariablesContainer{Cdouble},
    MOI.Utilities.MatrixOfConstraints{
        Cdouble,
        MOI.Utilities.MutableSparseMatrixCSC{
            Cdouble,
            T,
            MOI.Utilities.ZeroBasedIndexing,
        },
        Vector{Cdouble},
        Cones{Cdouble},
    },
}

function _to_sparse(
    ::Type{T},
    A::MOI.Utilities.MutableSparseMatrixCSC{
        Cdouble,
        T,
        MOI.Utilities.ZeroBasedIndexing,
    },
) where {T}
    return -A.nzval, A.rowval, A.colptr
end

mutable struct MOISolution
    primal::Vector{Float64}
    dual::Vector{Float64}
    slack::Vector{Float64}
    ret_val::Int
    raw_status::String
    objective_value::Float64
    dual_objective_value::Float64
    objective_constant::Float64
    solve_time_sec::Float64
    iterations::Int
end

function MOISolution()
    return MOISolution(
        Float64[],
        Float64[],
        Float64[],
        0, # SCS_UNFINISHED
        "",
        NaN,
        NaN,
        NaN,
        0.0,
        0,
    )
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    cones::Union{Nothing,Cones{Cdouble}}
    sol::MOISolution
    silent::Bool
    options::Dict{Symbol,Any}
    function Optimizer(; kwargs...)
        optimizer = new(nothing, MOISolution(), false, Dict{Symbol,Any}())
        for (key, value) in kwargs
            MOI.set(optimizer, MOI.RawOptimizerAttribute(String(key)), value)
        end
        return optimizer
    end
end

function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges)
    return [ScaledPSDConeBridge{Cdouble}]
end

MOI.get(::Optimizer, ::MOI.SolverName) = "SCS"

MOI.is_empty(optimizer::Optimizer) = optimizer.cones === nothing

function MOI.empty!(optimizer::Optimizer)
    optimizer.cones = nothing
    optimizer.sol = MOISolution()
    return
end

###
### MOI.RawOptimizerAttribute
###

function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value)
    return optimizer.options[Symbol(param.name)] = value
end

function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute)
    return optimizer.options[Symbol(param.name)]
end

###
### MOI.Silent
###

MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.silent = value
    return
end

MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

###
### MOI.AbstractModelAttribute
###

function MOI.supports(
    ::Optimizer,
    ::Union{
        MOI.ObjectiveSense,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    },
)
    return true
end

###
### MOI.AbstractVariableAttribute
###

# TODO(odow): support variable primal starts
# function MOI.supports(
#     ::Optimizer,
#     ::MOI.VariablePrimalStart,
#     ::Type{MOI.VariableIndex},
# )
#     return true
# end

###
### MOI.AbstractConstraintAttribute
###

# TODO(odow): support constraint primal and dual starts
# function MOI.supports(
#     ::Optimizer,
#     ::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
#     ::Type{<:MOI.ConstraintIndex},
# )
#     return true
# end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorAffineFunction{Cdouble}},
    ::Type{
        <:Union{
            MOI.Zeros,
            MOI.Nonnegatives,
            MOI.SecondOrderCone,
            ScaledPSDCone,
            MOI.ExponentialCone,
            MOI.DualExponentialCone,
            MOI.PowerCone{Cdouble},
            MOI.DualPowerCone{Cdouble},
        },
    },
)
    return true
end

function _map_sets(f, ::Type{T}, sets, ::Type{S}) where {T,S}
    F = MOI.VectorAffineFunction{Cdouble}
    cis = MOI.get(sets, MOI.ListOfConstraintIndices{F,S}())
    return T[f(MOI.get(sets, MOI.ConstraintSet(), ci)) for ci in cis]
end

function MOI.optimize!(
    dest::Optimizer,
    src::MOI.Utilities.UniversalFallback{OptimizerCache{T}},
) where {T}
    linear_solver = get(dest.options, :linear_solver, DirectSolver)
    if T != scsint_t(linear_solver)
        cache = MOI.Utilities.UniversalFallback(
            OptimizerCache{scsint_t(linear_solver)}(),
        )
        MOI.copy_to(cache, src)
        return MOI.optimize!(dest, cache)
    end
    # The real stuff starts here.
    MOI.empty!(dest)
    index_map = MOI.Utilities.identity_index_map(src)
    Ab = src.model.constraints
    A = Ab.coefficients

    model_attributes = MOI.get(src, MOI.ListOfModelAttributesSet())
    objective_function_attr =
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Cdouble}}()
    for attr in model_attributes
        if attr != MOI.ObjectiveSense() && attr != objective_function_attr
            throw(MOI.UnsupportedAttribute(attr))
        end
    end
    max_sense = false
    if MOI.ObjectiveSense() in model_attributes
        max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    end
    objective_constant = 0.0
    c = zeros(A.n)
    if objective_function_attr in model_attributes
        obj = MOI.get(src, objective_function_attr)
        objective_constant = MOI.constant(obj)
        for term in obj.terms
            c[term.variable.value] += (max_sense ? -1 : 1) * term.coefficient
        end
    end
    # `model.primal` contains the result of the previous optimization.
    # It is used as a warm start if its length is the same, e.g.
    # probably because no variable and/or constraint has been added.
    if A.n != length(dest.sol.primal)
        resize!(dest.sol.primal, A.n)
        fill!(dest.sol.primal, 0.0)
    end
    if A.m != length(dest.sol.dual)
        resize!(dest.sol.dual, A.m)
        fill!(dest.sol.dual, 0.0)
    end
    if A.m != length(dest.sol.slack)
        resize!(dest.sol.slack, A.m)
        fill!(dest.sol.slack, 0.0)
    end
    # Set starting values and throw error for other variable attributes
    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    MOI.Utilities.pass_attributes(dest, src, index_map, vis_src)
    # Set starting values and throw error for other constraint attributes
    for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
        cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
        MOI.Utilities.pass_attributes(dest, src, index_map, cis_src)
    end
    options = copy(dest.options)
    if dest.silent
        options[:verbose] = 0
    end
    dest.cones = deepcopy(Ab.sets)
    sol = scs_solve(
        linear_solver,
        A.m,
        A.n,
        A,
        Ab.constants,
        c,
        Ab.sets.num_rows[1],
        Ab.sets.num_rows[2] - Ab.sets.num_rows[1],
        _map_sets(MOI.dimension, T, Ab, MOI.SecondOrderCone),
        _map_sets(MOI.side_dimension, T, Ab, ScaledPSDCone),
        div(Ab.sets.num_rows[5] - Ab.sets.num_rows[4], 3),
        div(Ab.sets.num_rows[6] - Ab.sets.num_rows[5], 3),
        vcat(
            _map_sets(set -> set.exponent, Float64, Ab, MOI.PowerCone{Cdouble}),
            _map_sets(
                set -> -set.exponent,
                Float64,
                Ab,
                MOI.DualPowerCone{Cdouble},
            ),
        ),
        dest.sol.primal,
        dest.sol.dual,
        dest.sol.slack;
        options...,
    )
    dest.sol = MOISolution(
        sol.x,
        sol.y,
        sol.s,
        sol.ret_val,
        raw_status(sol.info),
        (max_sense ? -1 : 1) * sol.info.pobj,
        (max_sense ? -1 : 1) * sol.info.dobj,
        objective_constant,
        (sol.info.setupTime + sol.info.solveTime) / 1000,
        sol.info.iter,
    )
    return index_map, false
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    linear_solver = get(dest.options, :linear_solver, DirectSolver)
    T = scsint_t(linear_solver)
    cache = MOI.Utilities.UniversalFallback(OptimizerCache{T}())
    index_map = MOI.copy_to(cache, src)
    MOI.optimize!(dest, cache)
    return index_map, false
end

function MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec)
    return optimizer.sol.solve_time_sec
end

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return optimizer.sol.raw_status
end

"""
    ADMMIterations()

The number of ADMM iterations completed during the solve.
"""
struct ADMMIterations <: MOI.AbstractModelAttribute end

MOI.is_set_by_optimize(::ADMMIterations) = true

function MOI.get(optimizer::Optimizer, ::ADMMIterations)
    return optimizer.sol.iterations
end

# Implements getter for result value and statuses
# SCS returns one of the following integers:
# -7 SCS_INFEASIBLE_INACCURATE
# -6 SCS_UNBOUNDED_INACCURATE
# -5 SCS_SIGINT
# -4 SCS_FAILED
# -3 SCS_INDETERMINATE
# -2 SCS_INFEASIBLE  : primal infeasible, dual unbounded
# -1 SCS_UNBOUNDED   : primal unbounded, dual infeasible
#  0 SCS_UNFINISHED  : never returned, used as placeholder
#  1 SCS_SOLVED
#  2 SCS_SOLVED_INACCURATE
function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    s = optimizer.sol.ret_val
    @assert -7 <= s <= 2
    if s == -7
        return MOI.ALMOST_INFEASIBLE
    elseif s == -6
        return MOI.ALMOST_DUAL_INFEASIBLE
    elseif s == 2
        return MOI.ALMOST_OPTIMAL
    elseif s == -5
        return MOI.INTERRUPTED
    elseif s == -4
        return MOI.NUMERICAL_ERROR
    elseif s == -3
        return MOI.SLOW_PROGRESS
    elseif s == -2
        return MOI.INFEASIBLE
    elseif s == -1
        return MOI.DUAL_INFEASIBLE
    elseif s == 1
        return MOI.OPTIMAL
    else
        @assert s == 0
        return MOI.OPTIMIZE_NOT_CALLED
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.objective_value
    if !MOI.Utilities.is_ray(MOI.get(optimizer, MOI.PrimalStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.dual_objective_value
    if !MOI.Utilities.is_ray(MOI.get(optimizer, MOI.DualStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end

function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    elseif optimizer.sol.ret_val in (-3, 1, 2)
        return MOI.FEASIBLE_POINT
    elseif optimizer.sol.ret_val in (-6, -1)
        return MOI.INFEASIBILITY_CERTIFICATE
    end
    return MOI.INFEASIBLE_POINT
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.VariablePrimal,
    vi::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.primal[vi.value]
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.slack[MOI.Utilities.rows(optimizer.cones, ci)]
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    elseif optimizer.sol.ret_val in (-3, 1, 2)
        return MOI.FEASIBLE_POINT
    elseif optimizer.sol.ret_val in (-7, -2)
        return MOI.INFEASIBILITY_CERTIFICATE
    end
    return MOI.INFEASIBLE_POINT
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.dual[MOI.Utilities.rows(optimizer.cones, ci)]
end

MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1


text:
# ---
# title: JuliaRL\_DDPG\_Pendulum
# cover: assets/JuliaRL_DDPG_Pendulum.png
# description: DDPG applied to Pendulum
# date: 2021-05-22
# author: "[Jun Tian](https://github.com/findmyway)"
# ---

#+ tangle=true
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses
using IntervalSets

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:DDPG},
    ::Val{:Pendulum},
    ::Nothing;
    seed = 123,
)
    rng = StableRNG(seed)
    inner_env = PendulumEnv(T = Float32, rng = rng)
    A = action_space(inner_env)
    low = A.left
    high = A.right
    ns = length(state(inner_env))

    env = ActionTransformedEnv(
        inner_env;
        action_mapping = x -> low + (x + 1) * 0.5 * (high - low),
    )
    init = glorot_uniform(rng)

    create_actor() = Chain(
        Dense(ns, 30, relu; init = init),
        Dense(30, 30, relu; init = init),
        Dense(30, 1, tanh; init = init),
    )

    create_critic() = Chain(
        Dense(ns + 1, 30, relu; init = init),
        Dense(30, 30, relu; init = init),
        Dense(30, 1; init = init),
    )

    agent = Agent(
        policy = DDPGPolicy(
            behavior_actor = NeuralNetworkApproximator(
                model = create_actor(),
                optimizer = ADAM(),
            ),
            behavior_critic = NeuralNetworkApproximator(
                model = create_critic(),
                optimizer = ADAM(),
            ),
            target_actor = NeuralNetworkApproximator(
                model = create_actor(),
                optimizer = ADAM(),
            ),
            target_critic = NeuralNetworkApproximator(
                model = create_critic(),
                optimizer = ADAM(),
            ),
            γ = 0.99f0,
            ρ = 0.995f0,
            batch_size = 64,
            start_steps = 1000,
            start_policy = RandomPolicy(-1.0..1.0; rng = rng),
            update_after = 1000,
            update_every = 1,
            act_limit = 1.0,
            act_noise = 0.1,
            rng = rng,
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = 10000,
            state = Vector{Float32} => (ns,),
            action = Float32 => (),
        ),
    )

    stop_condition = StopAfterStep(10_000, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalRewardPerEpisode()
    Experiment(agent, env, stop_condition, hook, "# Play Pendulum with DDPG")
end

#+ tangle=false
using Plots
using Statistics
pyplot() #hide
ex = E`JuliaRL_DDPG_Pendulum`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_DDPG_Pendulum.png") #hide

# ![](assets/JuliaRL_DDPG_Pendulum.png)


text:
using ThreeBodyTB
using Test
using Suppressor

SRCDIR=ThreeBodyTB.SRCDIR

@testset "testing prototypes" begin
    @suppress begin
        include("$SRCDIR/../code_for_dataset_gen/Prototypes.jl")
        t =  oxidation_guess("Na", "Cl")
        tref = [["Na", "Cl", :core_binary],["Na", "Cl", :A1B1], ["Na", "Cl", "pd3te"], ["Na", "Cl", "bcc_5lay"], ["Na", "Cl", "fcc_5lay"]]

        @test t == tref

        t =  oxidation_guess("Na", "Mg")
        tref = [["Na", "Mg", :core_binary],["Na", "Mg", :metals], ["Na", "Mg", "hh_oxygen"], ["Na", "Mg", "hh_oxygen2"], ["Na", "Mg", "hex_oxygen"]]

        @test t == tref

        ThreeBodyTB.set_units(both="eVAng")

    end
end


text:
using PrivateMultiplicativeWeights

function marginals(d=20, order=3, n=1000)
    data_matrix = rand(0:1, d, n)
    data_matrix[3, :] = data_matrix[1, :] .* data_matrix[2, :]
    mwem(Parities(d, order), Tabular(data_matrix), MWParameters(verbose=true))
end

function factored_marginals(d=20, order=3, n=1000)
    data_matrix = rand(0:1, d, n)
    data_matrix[3, :] = data_matrix[1, :] .* data_matrix[2, :]
    mwem(FactorParities(d, order), Tabular(data_matrix), MWParameters(verbose=true))
end

function range_queries()
    histogram = vcat(zeros(100), ones(100), zeros(100))
    data = Histogram(histogram, 100)
    mwem(SeriesRangeQueries(300), data, MWParameters(repetitions=100, verbose=true))
end

print("Elapsed time: ", @elapsed marginals())
print("\n\n")

print("Elapsed time: ", @elapsed factored_marginals())
print("\n\n")

print("Elapsed time: ", @elapsed range_queries())
print("\n\n")


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

struct __enum_ZigBeeMsg_ZigBeeCmdType <: ProtoEnum
    ZCL::Int32
    NETWORK_MGMT::Int32
    __enum_ZigBeeMsg_ZigBeeCmdType() = new(0,1)
end #struct __enum_ZigBeeMsg_ZigBeeCmdType
const ZigBeeMsg_ZigBeeCmdType = __enum_ZigBeeMsg_ZigBeeCmdType()

struct __enum_ZigBeeMsg_ZCLCmd_ZCLCmdType <: ProtoEnum
    ON_OFF::Int32
    LEVEL::Int32
    COLOR_CONTROL::Int32
    IDENTIFY::Int32
    __enum_ZigBeeMsg_ZCLCmd_ZCLCmdType() = new(0,1,2,3)
end #struct __enum_ZigBeeMsg_ZCLCmd_ZCLCmdType
const ZigBeeMsg_ZCLCmd_ZCLCmdType = __enum_ZigBeeMsg_ZCLCmd_ZCLCmdType()

struct __enum_ZigBeeMsg_ZCLCmd_OnOffCmd_ZCLOnOffCmdType <: ProtoEnum
    ON::Int32
    OFF::Int32
    TOGGLE::Int32
    __enum_ZigBeeMsg_ZCLCmd_OnOffCmd_ZCLOnOffCmdType() = new(0,1,2)
end #struct __enum_ZigBeeMsg_ZCLCmd_OnOffCmd_ZCLOnOffCmdType
const ZigBeeMsg_ZCLCmd_OnOffCmd_ZCLOnOffCmdType = __enum_ZigBeeMsg_ZCLCmd_OnOffCmd_ZCLOnOffCmdType()

mutable struct ZigBeeMsg_ZCLCmd_OnOffCmd <: ProtoType
    _type::Int32
    ZigBeeMsg_ZCLCmd_OnOffCmd(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_OnOffCmd

struct __enum_ZigBeeMsg_ZCLCmd_LevelCmd_ZCLLevelCmdType <: ProtoEnum
    MOVE_TO_LEVEL::Int32
    MOVE::Int32
    __enum_ZigBeeMsg_ZCLCmd_LevelCmd_ZCLLevelCmdType() = new(0,1)
end #struct __enum_ZigBeeMsg_ZCLCmd_LevelCmd_ZCLLevelCmdType
const ZigBeeMsg_ZCLCmd_LevelCmd_ZCLLevelCmdType = __enum_ZigBeeMsg_ZCLCmd_LevelCmd_ZCLLevelCmdType()

mutable struct ZigBeeMsg_ZCLCmd_LevelCmd_MoveToLevelCmdParams <: ProtoType
    level::UInt32
    transition_time::UInt32
    ZigBeeMsg_ZCLCmd_LevelCmd_MoveToLevelCmdParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_LevelCmd_MoveToLevelCmdParams

mutable struct ZigBeeMsg_ZCLCmd_LevelCmd_MoveCmdParams <: ProtoType
    mode::UInt32
    rate::UInt32
    ZigBeeMsg_ZCLCmd_LevelCmd_MoveCmdParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_LevelCmd_MoveCmdParams

mutable struct ZigBeeMsg_ZCLCmd_LevelCmd <: ProtoType
    _type::Int32
    move_to_level_params::ZigBeeMsg_ZCLCmd_LevelCmd_MoveToLevelCmdParams
    move_params::ZigBeeMsg_ZCLCmd_LevelCmd_MoveCmdParams
    ZigBeeMsg_ZCLCmd_LevelCmd(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_LevelCmd

struct __enum_ZigBeeMsg_ZCLCmd_ColorControlCmd_ZCLColorControlCmdType <: ProtoEnum
    MOVETOHUE::Int32
    MOVETOSAT::Int32
    MOVETOHUEANDSAT::Int32
    MOVETOCOLORTEMP::Int32
    __enum_ZigBeeMsg_ZCLCmd_ColorControlCmd_ZCLColorControlCmdType() = new(0,1,2,3)
end #struct __enum_ZigBeeMsg_ZCLCmd_ColorControlCmd_ZCLColorControlCmdType
const ZigBeeMsg_ZCLCmd_ColorControlCmd_ZCLColorControlCmdType = __enum_ZigBeeMsg_ZCLCmd_ColorControlCmd_ZCLColorControlCmdType()

struct __enum_ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams_DirectionParam <: ProtoEnum
    SHORTEST_DISTANCE::Int32
    LONGEST_DISTANCE::Int32
    UP::Int32
    DOWN::Int32
    __enum_ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams_DirectionParam() = new(0,1,2,3)
end #struct __enum_ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams_DirectionParam
const ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams_DirectionParam = __enum_ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams_DirectionParam()

mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams <: ProtoType
    hue::UInt32
    direction::Int32
    transition_time::UInt32
    ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams

mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToSatCmdParams <: ProtoType
    saturation::UInt32
    transition_time::UInt32
    ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToSatCmdParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToSatCmdParams

mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueAndSatCmdParams <: ProtoType
    hue::UInt32
    saturation::UInt32
    transition_time::UInt32
    ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueAndSatCmdParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueAndSatCmdParams

mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToColorTempCmdParams <: ProtoType
    color_temperature::UInt32
    transition_time::UInt32
    ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToColorTempCmdParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToColorTempCmdParams

mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd <: ProtoType
    _type::Int32
    movetohue_params::ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams
    movetosat_params::ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToSatCmdParams
    movetohueandsat_params::ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueAndSatCmdParams
    movetocolortemp_params::ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToColorTempCmdParams
    ZigBeeMsg_ZCLCmd_ColorControlCmd(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_ColorControlCmd

struct __enum_ZigBeeMsg_ZCLCmd_IdentifyCmd_ZCLIdentifyCmdType <: ProtoEnum
    IDENTIFY_ON::Int32
    IDENTIFY_OFF::Int32
    __enum_ZigBeeMsg_ZCLCmd_IdentifyCmd_ZCLIdentifyCmdType() = new(0,1)
end #struct __enum_ZigBeeMsg_ZCLCmd_IdentifyCmd_ZCLIdentifyCmdType
const ZigBeeMsg_ZCLCmd_IdentifyCmd_ZCLIdentifyCmdType = __enum_ZigBeeMsg_ZCLCmd_IdentifyCmd_ZCLIdentifyCmdType()

mutable struct ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOnCmdParams <: ProtoType
    endpoint::UInt32
    identify_time::UInt32
    ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOnCmdParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOnCmdParams

mutable struct ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOffCmdParams <: ProtoType
    identify_time::UInt32
    ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOffCmdParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOffCmdParams

mutable struct ZigBeeMsg_ZCLCmd_IdentifyCmd <: ProtoType
    _type::Int32
    identify_on_params::ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOnCmdParams
    identify_off_params::ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOffCmdParams
    ZigBeeMsg_ZCLCmd_IdentifyCmd(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd_IdentifyCmd

mutable struct ZigBeeMsg_ZCLCmd <: ProtoType
    _type::Int32
    onoff_cmd::ZigBeeMsg_ZCLCmd_OnOffCmd
    level_cmd::ZigBeeMsg_ZCLCmd_LevelCmd
    colorcontrol_cmd::ZigBeeMsg_ZCLCmd_ColorControlCmd
    identify_cmd::ZigBeeMsg_ZCLCmd_IdentifyCmd
    node_id::UInt32
    endpoint_index::UInt32
    ZigBeeMsg_ZCLCmd(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_ZCLCmd

struct __enum_ZigBeeMsg_NetworkMgmtCmd_NetworkMgmtCmdTypes <: ProtoEnum
    CREATE_NWK::Int32
    LEAVE_NWK::Int32
    NODE_LEAVE_NWK::Int32
    PERMIT_JOIN::Int32
    DISCOVERY_INFO::Int32
    RESET_PROXY::Int32
    IS_PROXY_ACTIVE::Int32
    NETWORK_STATUS::Int32
    __enum_ZigBeeMsg_NetworkMgmtCmd_NetworkMgmtCmdTypes() = new(0,1,2,3,4,5,6,7)
end #struct __enum_ZigBeeMsg_NetworkMgmtCmd_NetworkMgmtCmdTypes
const ZigBeeMsg_NetworkMgmtCmd_NetworkMgmtCmdTypes = __enum_ZigBeeMsg_NetworkMgmtCmd_NetworkMgmtCmdTypes()

mutable struct ZigBeeMsg_NetworkMgmtCmd_PermitJoinParams <: ProtoType
    time::UInt32
    ZigBeeMsg_NetworkMgmtCmd_PermitJoinParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_NetworkMgmtCmd_PermitJoinParams

mutable struct ZigBeeMsg_NetworkMgmtCmd_NodeLeaveNetParams <: ProtoType
    node_id::UInt32
    ZigBeeMsg_NetworkMgmtCmd_NodeLeaveNetParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_NetworkMgmtCmd_NodeLeaveNetParams

struct __enum_ZigBeeMsg_NetworkMgmtCmd_ClusterDescription_ClusterInOutType <: ProtoEnum
    SERVER_IN::Int32
    CLIENT_OUT::Int32
    __enum_ZigBeeMsg_NetworkMgmtCmd_ClusterDescription_ClusterInOutType() = new(0,1)
end #struct __enum_ZigBeeMsg_NetworkMgmtCmd_ClusterDescription_ClusterInOutType
const ZigBeeMsg_NetworkMgmtCmd_ClusterDescription_ClusterInOutType = __enum_ZigBeeMsg_NetworkMgmtCmd_ClusterDescription_ClusterInOutType()

mutable struct ZigBeeMsg_NetworkMgmtCmd_ClusterDescription <: ProtoType
    _type::Int32
    cluster_id::UInt32
    ZigBeeMsg_NetworkMgmtCmd_ClusterDescription(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_NetworkMgmtCmd_ClusterDescription

mutable struct ZigBeeMsg_NetworkMgmtCmd_EndPointDescription <: ProtoType
    endpoint_index::UInt32
    profile_id::UInt32
    device_id::UInt32
    clusters::Base.Vector{ZigBeeMsg_NetworkMgmtCmd_ClusterDescription}
    ZigBeeMsg_NetworkMgmtCmd_EndPointDescription(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_NetworkMgmtCmd_EndPointDescription

mutable struct ZigBeeMsg_NetworkMgmtCmd_NodeDescription <: ProtoType
    node_id::UInt32
    eui64::UInt64
    endpoints::Base.Vector{ZigBeeMsg_NetworkMgmtCmd_EndPointDescription}
    ZigBeeMsg_NetworkMgmtCmd_NodeDescription(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_NetworkMgmtCmd_NodeDescription

struct __enum_ZigBeeMsg_NetworkMgmtCmd_NetworkStatus_Status <: ProtoEnum
    NO_NETWORK::Int32
    JOINING_NETWORK::Int32
    JOINED_NETWORK::Int32
    JOINED_NETWORK_NO_PARENT::Int32
    LEAVING_NETWORK::Int32
    __enum_ZigBeeMsg_NetworkMgmtCmd_NetworkStatus_Status() = new(0,1,2,3,4)
end #struct __enum_ZigBeeMsg_NetworkMgmtCmd_NetworkStatus_Status
const ZigBeeMsg_NetworkMgmtCmd_NetworkStatus_Status = __enum_ZigBeeMsg_NetworkMgmtCmd_NetworkStatus_Status()

mutable struct ZigBeeMsg_NetworkMgmtCmd_NetworkStatus <: ProtoType
    _type::Int32
    ZigBeeMsg_NetworkMgmtCmd_NetworkStatus(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_NetworkMgmtCmd_NetworkStatus

mutable struct ZigBeeMsg_NetworkMgmtCmd <: ProtoType
    _type::Int32
    permit_join_params::ZigBeeMsg_NetworkMgmtCmd_PermitJoinParams
    is_proxy_active::Bool
    node_leave_params::ZigBeeMsg_NetworkMgmtCmd_NodeLeaveNetParams
    node_info::ZigBeeMsg_NetworkMgmtCmd_NodeDescription
    connected_nodes::Base.Vector{ZigBeeMsg_NetworkMgmtCmd_NodeDescription}
    network_status::ZigBeeMsg_NetworkMgmtCmd_NetworkStatus
    ZigBeeMsg_NetworkMgmtCmd(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg_NetworkMgmtCmd

mutable struct ZigBeeMsg <: ProtoType
    _type::Int32
    zcl_cmd::ZigBeeMsg_ZCLCmd
    network_mgmt_cmd::ZigBeeMsg_NetworkMgmtCmd
    ZigBeeMsg(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZigBeeMsg

struct __enum_ZWaveMsg_ZWaveOperations <: ProtoEnum
    UNDEF::Int32
    SEND::Int32
    ADDNODE::Int32
    REMOVENODE::Int32
    SETDEFAULT::Int32
    LIST::Int32
    __enum_ZWaveMsg_ZWaveOperations() = new(0,1,2,3,4,5)
end #struct __enum_ZWaveMsg_ZWaveOperations
const ZWaveMsg_ZWaveOperations = __enum_ZWaveMsg_ZWaveOperations()

mutable struct ZWaveMsg_ZWaveCommandInfo <: ProtoType
    cmd::Int32
    param::Base.Vector{AbstractString}
    ZWaveMsg_ZWaveCommandInfo(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZWaveMsg_ZWaveCommandInfo

mutable struct ZWaveMsg_ZWaveClassInfo <: ProtoType
    zwave_class::Int32
    command::Base.Vector{ZWaveMsg_ZWaveCommandInfo}
    ZWaveMsg_ZWaveClassInfo(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZWaveMsg_ZWaveClassInfo

mutable struct ZWaveMsg_ZWaveNode <: ProtoType
    service_name::AbstractString
    zwave_class::Base.Vector{ZWaveMsg_ZWaveClassInfo}
    ZWaveMsg_ZWaveNode(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZWaveMsg_ZWaveNode

mutable struct ZWaveMsg_ZWaveCommand <: ProtoType
    zwclass::Int32
    cmd::Int32
    params::Array{UInt8,1}
    ZWaveMsg_ZWaveCommand(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZWaveMsg_ZWaveCommand

mutable struct ZWaveMsg <: ProtoType
    operation::Int32
    service_to_send::AbstractString
    zwave_cmd::ZWaveMsg_ZWaveCommand
    node::Base.Vector{ZWaveMsg_ZWaveNode}
    result::Bool
    ZWaveMsg(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ZWaveMsg

mutable struct LircParams <: ProtoType
    device::AbstractString
    command::AbstractString
    config::AbstractString
    LircParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct LircParams

export ZigBeeMsg_ZigBeeCmdType, ZigBeeMsg_ZCLCmd_ZCLCmdType, ZigBeeMsg_ZCLCmd_OnOffCmd_ZCLOnOffCmdType, ZigBeeMsg_ZCLCmd_OnOffCmd, ZigBeeMsg_ZCLCmd_LevelCmd_ZCLLevelCmdType, ZigBeeMsg_ZCLCmd_LevelCmd_MoveToLevelCmdParams, ZigBeeMsg_ZCLCmd_LevelCmd_MoveCmdParams, ZigBeeMsg_ZCLCmd_LevelCmd, ZigBeeMsg_ZCLCmd_ColorControlCmd_ZCLColorControlCmdType, ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams_DirectionParam, ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueCmdParams, ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToSatCmdParams, ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToHueAndSatCmdParams, ZigBeeMsg_ZCLCmd_ColorControlCmd_MoveToColorTempCmdParams, ZigBeeMsg_ZCLCmd_ColorControlCmd, ZigBeeMsg_ZCLCmd_IdentifyCmd_ZCLIdentifyCmdType, ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOnCmdParams, ZigBeeMsg_ZCLCmd_IdentifyCmd_IdentifyOffCmdParams, ZigBeeMsg_ZCLCmd_IdentifyCmd, ZigBeeMsg_ZCLCmd, ZigBeeMsg_NetworkMgmtCmd_NetworkMgmtCmdTypes, ZigBeeMsg_NetworkMgmtCmd_PermitJoinParams, ZigBeeMsg_NetworkMgmtCmd_NodeLeaveNetParams, ZigBeeMsg_NetworkMgmtCmd_ClusterDescription_ClusterInOutType, ZigBeeMsg_NetworkMgmtCmd_ClusterDescription, ZigBeeMsg_NetworkMgmtCmd_EndPointDescription, ZigBeeMsg_NetworkMgmtCmd_NodeDescription, ZigBeeMsg_NetworkMgmtCmd_NetworkStatus_Status, ZigBeeMsg_NetworkMgmtCmd_NetworkStatus, ZigBeeMsg_NetworkMgmtCmd, ZigBeeMsg, ZWaveMsg_ZWaveOperations, ZWaveMsg_ZWaveCommandInfo, ZWaveMsg_ZWaveClassInfo, ZWaveMsg_ZWaveNode, ZWaveMsg_ZWaveCommand, ZWaveMsg, LircParams


text:
module NonIterableCellMapsTests

using Test
using Gridap

struct Foo <: NonIterableCellMap{Point{2},1,Int,1} end

foo = Foo()

@test isa(foo,CellMap)

try
  iterate(foo)
catch LoadError
end

try
  iterate(foo,2)
catch LoadError
end

end # module