text:
include("global_curved.jl")
import PGFPlots

let
  method = :trace

  # SBP interior order
  SBPp   = 6

  # mesh file side set type to actually boundary condition type
  bc_map = [BC_DIRICHLET, BC_DIRICHLET, BC_NEUMANN, BC_NEUMANN,
            BC_JUMP_INTERFACE]
  (verts, EToV, EToF, FToB, EToDomain) = read_inp_2d("meshes/square_circle.inp";
                                                     bc_map = bc_map)
  # EToV defines the element by its vertices
  # EToF defines element by its four faces, in global face number
  # FToB defines whether face is Dirichlet (1), Neumann (2), interior jump (7)
  #      or just an interior interface (0)
  # EToDomain is 1 if element is inside circle; 2 otherwise

  # number of elements and faces
  (nelems, nfaces) = (size(EToV, 2), size(FToB, 1))
  @show (nelems, nfaces)

  # This is needed to fix up points that should be on the boundary of the
  # circle, but the mesh didn't quite put them there
  for v in 1:size(verts, 2)
    x, y = verts[1, v], verts[2, v]
    if abs(hypot(x,y) - 1) < 1e-5
      Q = atan(y, x)
      verts[1, v], verts[2, v] = cos(Q), sin(Q)
    end
  end

  # Plot the original connectivity before mesh warping
  # plot_connectivity(verts, EToV)

  # This is the base mesh size in each dimension
  N1 = N0 = 17

  # EToN0 is the base mesh size (e.g., before refinement)
  EToN0 = zeros(Int64, 2, nelems)
  EToN0[1, :] .= N0
  EToN0[2, :] .= N1

  @assert typeof(EToV) == Array{Int, 2} && size(EToV) == (4, nelems)
  @assert typeof(EToF) == Array{Int, 2} && size(EToF) == (4, nelems)
  @assert maximum(maximum(EToF)) == nfaces

  # Determine secondary arrays
  # FToE : Unique Global Face to Element Number
  #        (the i'th column of this stores the element numbers that share the
  #        global face number i)
  # FToLF: Unique Global Face to Element local face number
  #        (the i'th column of this stores the element local face numbers that
  #        shares the global face number i)
  # EToO : Element to Unique Global Faces Orientation
  #        (the i'th column of this stores the whether the element and global
  #        face are oriented in the same way in physical memory or need to be
  #        rotated)
  # EToS : Element to Unique Global Face Side
  #        (the i'th column of this stores whether an element face is on the
  #        plus side or minus side of the global face)
  (FToE, FToLF, EToO, EToS) = connectivityarrays(EToV, EToF)

  pgf_axis = PGFPlots.Axis(style="width=5cm, height=5cm, ticks=none",
                           xlabel=PGFPlots.L"$x$",
                           ylabel=PGFPlots.L"$y$",
                           xmin = -2, xmax = 2,
                           ymin = -2, ymax = 2)
  for f in 1:nfaces
    if FToB[f] != BC_JUMP_INTERFACE
      (e, lf) = FToE[1,f], FToLF[1,f]
      if lf == 1
        v1, v2 = EToV[1, e], EToV[3, e]
      elseif lf == 2
        v1, v2 = EToV[2, e], EToV[4, e]
      elseif lf == 3
        v1, v2 = EToV[1, e], EToV[2, e]
      else
        v1, v2 = EToV[3, e], EToV[4, e]
      end
      x = verts[1, [v1 v2]][:]
      y = verts[2, [v1 v2]][:]
      push!(pgf_axis, PGFPlots.Linear(x, y, style="no marks, solid, black"))
    end
  end
  push!(pgf_axis, PGFPlots.Circle(0,0,1, style = "very thick, red"))
  PGFPlots.save("square_circle.tikz", pgf_axis)

  # Exact solution
  Lx = maximum(verts[1,:])
  Ly = maximum(abs.(verts[2,:]))

  c = exp(1)/(1 + exp(1))

  vinside(x,y,e) = begin
    r = sqrt.(x .^ 2 + y .^ 2)
    theta = atan.(y,x)
    return c .* (1 .- exp.(-1 .* r .^ 2)) .* r .* sin.(theta)
  end

  voutside(x,y,e) = begin
    r = sqrt.(x .^ 2 + y .^ 2)
    theta = atan.(y,x)
    return (r .- 1) .^ 2 .* cos.(theta) .+ (r .- 1) .* sin.(theta)
  end

  vinside_x(x,y,e) = begin
    r = sqrt.(x .^ 2 + y .^ 2)
    theta = atan.(y,x)
    dtheta_dx = -1 .* sin.(theta) ./ r
    dr_dx = cos.(theta)
    dv_dr = c * (2 .* r .^ 2 .* exp.(-1 .* r .^ 2) .+ 1 .- exp.(-1 .* r .^ 2)) .* sin.(theta)
    dv_dtheta = c .* (1 .- exp.(-1 .* r .^ 2)) .* r .* cos.(theta)
    return dv_dr .* dr_dx + dv_dtheta .* dtheta_dx
  end

  vinside_y(x,y,e) = begin
    r = sqrt.(x .^ 2 + y .^ 2)
    theta = atan.(y,x)
    dtheta_dy = cos.(theta) ./ r
    dr_dy = sin.(theta)
    dv_dr = c * (2 .* r .^ 2 .* exp.(-1 .* r .^ 2) .+ 1 .- exp.(-1 .* r .^ 2)) .* sin.(theta)
    dv_dtheta = c .* (1 .- exp.(-1 .* r .^ 2)) .* r .* cos.(theta)
    return dv_dr .* dr_dy + dv_dtheta .* dtheta_dy
  end

  voutside_x(x,y,e) = begin
    r = sqrt.(x .^ 2 + y .^ 2)
    theta = atan.(y,x)
    dtheta_dx = -1 .* sin.(theta) ./ r
    dr_dx = cos.(theta)
    dv_dr =  2 .* (r .- 1) .* cos.(theta) .+ sin.(theta)
    dv_dtheta = -1 .* (r .- 1) .^ 2 .* sin.(theta) .+ (r .- 1) .* cos.(theta)
    return dv_dr .* dr_dx + dv_dtheta .* dtheta_dx
  end

  voutside_y(x,y,e) = begin
    r = sqrt.(x .^ 2 + y .^ 2)
    theta = atan.(y,x)
    dtheta_dy = cos.(theta) ./ r
    dr_dy = sin.(theta)
    dv_dr = 2 .* (r .- 1) .* cos.(theta) .+ sin.(theta)
    dv_dtheta = -1 .*  (r .- 1) .^ 2 .* sin.(theta) .+ (r .- 1) .* cos.(theta)
    return dv_dr .* dr_dy + dv_dtheta .* dtheta_dy
  end

  polar_laplace(x,y,e) = begin #u_rr + (1/r)*u_r + (1/r^2)*u_theta,theta
    r = sqrt.(x .^ 2 + y .^ 2)
    theta = atan.(y,x)
    if EToDomain[e] == 1
      u_r = c .* (2 .* r .^ 2 .* exp.(-1 .* r .^ 2) .+ 1 .- exp.(-1 .* r .^ 2)) .* sin.(theta)
      u_rr = c .* exp.(-1 .* r .^ 2) .* (6 .* r .- 4 .* r .^ 3) .* sin.(theta)
      return u_rr .+ (1 ./ r) .* u_r .- (c ./ r .^ 2) .* (1 .- exp.(-1 .* r .^ 2)) .* r .* sin.(theta)
      return dv_dr .* dr_dy .+ (cos.(r .- 1) .* cos.(theta)) .* (cos.(theta))
    elseif EToDomain[e] == 2
      return 2 * cos.(theta) .+ (1 ./ r) .* (2 .* (r .- 1) .* cos.(theta) .+ sin.(theta)) .+ (1 ./ r .^ 2) .* (-1 .* (r .- 1) .^ 2 .* cos.(theta) .- (r .- 1) .* sin.(theta))
    else
      error("invalid block")
    end
  end

  vex(x,y,e) = begin
    if EToDomain[e] == 1
      return vinside(x,y,e)
    elseif EToDomain[e] == 2
      return voutside(x,y,e)
    else
      error("invalid block")
    end
  end

  vex_x(x,y,e) = begin
    if EToDomain[e] == 1
      return vinside_x(x,y,e)
    elseif EToDomain[e] == 2
      return voutside_x(x,y,e)
    else
      error("invalid block")
    end
  end

  vex_y(x,y,e) = begin
    if EToDomain[e] == 1
      return vinside_y(x,y,e)
    elseif EToDomain[e] == 2
      return voutside_y(x,y,e)
    else
      error("invalid block")
    end
  end

  ϵ = zeros(4)
  τϵ = copy(ϵ)

  for lvl = 1:length(ϵ)
    # Set up the local grid dimensions
    Nr = EToN0[1, :] * (2^(lvl-1))
    Ns = EToN0[2, :] * (2^(lvl-1))

    #
    # Build the local volume operators
    #

    # Dictionary to store the operators
    OPTYPE = typeof(locoperator(2, 16, 16))
    lop = Dict{Int64, OPTYPE}()

    # Loop over blocks and create local operators
    for e = 1:nelems
      # Get the element corners
      (x1, x2, x3, x4) = verts[1, EToV[:, e]]
      (y1, y2, y3, y4) = verts[2, EToV[:, e]]

      # Initialize the block transformations as transfinite between the corners
      ex = [(α) -> x1 * (1 .- α) / 2 + x3 * (1 .+ α) / 2,
            (α) -> x2 * (1 .- α) / 2 + x4 * (1 .+ α) / 2,
            (α) -> x1 * (1 .- α) / 2 + x2 * (1 .+ α) / 2,
            (α) -> x3 * (1 .- α) / 2 + x4 * (1 .+ α) / 2]
      exα = [(α) -> -x1 / 2 + x3 / 2,
             (α) -> -x2 / 2 + x4 / 2,
             (α) -> -x1 / 2 + x2 / 2,
             (α) -> -x3 / 2 + x4 / 2]
      ey = [(α) -> y1 * (1 .- α) / 2 + y3 * (1 .+ α) / 2,
            (α) -> y2 * (1 .- α) / 2 + y4 * (1 .+ α) / 2,
            (α) -> y1 * (1 .- α) / 2 + y2 * (1 .+ α) / 2,
            (α) -> y3 * (1 .- α) / 2 + y4 * (1 .+ α) / 2]
      eyα = [(α) -> -y1 / 2 + y3 / 2,
             (α) -> -y2 / 2 + y4 / 2,
             (α) -> -y1 / 2 + y2 / 2,
             (α) -> -y3 / 2 + y4 / 2]

      # For blocks on the circle, put in the curved edge transform
      if FToB[EToF[1, e]] == BC_JUMP_INTERFACE
        error("curved face 1 not implemented yet")
      end
      if FToB[EToF[2, e]] == BC_JUMP_INTERFACE
        error("curved face 2 not implemented yet")
      end
      if FToB[EToF[3, e]] == BC_JUMP_INTERFACE
        Q1 = atan(y1, x1)
        Q2 = atan(y2, x2)
        if !(-π/2 < Q1 - Q2 < π/2)
          Q2 -= sign(Q2) * 2 * π
        end
        ex[3] = (α) -> cos.(Q1 * (1 .- α) / 2 + Q2 * (1 .+ α) / 2)
        ey[3] = (α) -> sin.(Q1 * (1 .- α) / 2 + Q2 * (1 .+ α) / 2)
        β3 = (Q2 - Q1) / 2
        exα[3] = (α) -> -β3 .* sin.(Q1 * (1 .- α) / 2 + Q2 * (1 .+ α) / 2)
        eyα[3] = (α) -> +β3 .* cos.(Q1 * (1 .- α) / 2 + Q2 * (1 .+ α) / 2)
      end
      if FToB[EToF[4, e]] == BC_JUMP_INTERFACE
        Q3 = atan(y3, x3)
        Q4 = atan(y4, x4)
        if !(-π/2 < Q3 - Q4 < π/2)
          error("curved face 4 angle correction not implemented yet")
        end
        ex[4] = (α) -> cos.(Q3 * (1 .- α) / 2 + Q4 * (1 .+ α) / 2)
        ey[4] = (α) -> sin.(Q3 * (1 .- α) / 2 + Q4 * (1 .+ α) / 2)
        β4 = (Q4 - Q3) / 2
        exα[4] = (α) -> -β4 .* sin.(Q3 * (1 .- α) / 2 + Q4 * (1 .+ α) / 2)
        eyα[4] = (α) -> +β4 .* cos.(Q3 * (1 .- α) / 2 + Q4 * (1 .+ α) / 2)
      end

      # Create the volume transform as the transfinite blending of the edge
      # transformations
      xt(r,s) = transfinite_blend(ex[1], ex[2], ex[3], ex[4],
                                  exα[1], exα[2], exα[3], exα[4],
                                  r, s)
      yt(r,s) = transfinite_blend(ey[1], ey[2], ey[3], ey[4],
                                  eyα[1], eyα[2], eyα[3], eyα[4],
                                  r, s)

      metrics = create_metrics(SBPp, Nr[e], Ns[e], xt, yt)

      # Build local operators
      lop[e] = locoperator(SBPp, Nr[e], Ns[e], metrics, FToB[EToF[:, e]])
    end

    # If this is the first mesh level plot the mesh
    lvl == 1 && plot_blocks(lop)

    #
    # Do some assemble of the global volume operators
    #
    (M, FbarT, D, vstarts, FToλstarts) =
    LocalGlobalOperators(lop, Nr, Ns, FToB, FToE, FToLF, EToO, EToS,
                         (x) -> cholesky(Symmetric(x)))
    @show lvl
    locfactors = M.F

    # Get a unique array indexes for the face to jumps map
    FToδstarts = bcstarts(FToB, FToE, FToLF, BC_JUMP_INTERFACE, Nr, Ns)

    # Compute the number of volume, trace (λ), and jump (δ) points
    VNp = vstarts[nelems+1]-1
    λNp = FToλstarts[nfaces+1]-1
    δNp = FToδstarts[nfaces+1]-1
    # @show (VNp, λNp, δNp)

    # Build the (sparse) λ matrix using the schur complement and factor
    B = assembleλmatrix(FToλstarts, vstarts, EToF, FToB, locfactors, D, FbarT)
    BF = cholesky(Symmetric(B))

    (bλ, λ, gδ) = (zeros(λNp), zeros(λNp), zeros(λNp))
    (Δ, u, g) = (zeros(VNp), zeros(VNp), zeros(VNp))
    δ = zeros(δNp)

    for f = 1:nfaces
      if FToB[f] == BC_JUMP_INTERFACE
        (e1, e2) = FToE[:, f]
        (lf1, lf2) = FToLF[:, f]
        (xf, yf) = lop[e1].facecoord
        @views δ[FToδstarts[f]:(FToδstarts[f+1]-1)] =
        vex(xf[lf1], yf[lf1], e2) - vex(xf[lf1], yf[lf1], e1)
        #@show δ[FToδstarts[f]:(FToδstarts[f+1]-1)]
      end
    end

    bc_Dirichlet = (lf, x, y, e, δ) -> vex(x, y, e)
    bc_Neumann   = (lf, x, y, nx, ny, e, δ) -> (nx .* vex_x(x, y, e)
                                                + ny .* vex_y(x, y, e))
    in_jump      = (lf, x, y, e, δ) -> begin
      f = EToF[lf, e]
      if EToS[lf, e] == 1
        if EToO[lf, e]
          return -δ[FToδstarts[f]:(FToδstarts[f+1]-1)]
        else
          error("shouldn't get here")
        end
      else
        if EToO[lf, e]
          return  δ[FToδstarts[f]:(FToδstarts[f+1]-1)]
        else
          return  δ[(FToδstarts[f+1]-1):-1:FToδstarts[f]]
        end
      end
    end

    for e = 1:nelems
      gδe = ntuple(4) do lf
        f = EToF[lf, e]
        if EToO[lf, e]
          return @view gδ[FToλstarts[f]:(FToλstarts[f+1]-1)]
        else
          return  @view gδ[(FToλstarts[f+1]-1):-1:FToλstarts[f]]
        end
      end
      locbcarray!((@view g[vstarts[e]:vstarts[e+1]-1]), gδe, lop[e],
                  FToB[EToF[:,e]], bc_Dirichlet, bc_Neumann, in_jump, (e, δ))

      source = (x, y, e) -> (-polar_laplace(x, y, e))
      locsourcearray!((@view g[vstarts[e]:vstarts[e+1]-1]), source, lop[e], e)
    end

    if method == :monolithic
      # Solve the monolithic sytem
      M = blockdiag(ntuple(i->lop[i].M̃, length(lop))...)
      A = [M FbarT'; FbarT Diagonal(D)]
      uλ = A \ [g;gδ]
      u[:] .= uλ[1:length(u)]
    elseif method == :trace
      # Solve for the trace variables then compute displacements
      LocalToGLobalRHS!(bλ, g, gδ, u, locfactors, FbarT, vstarts)
      λ[:] = BF \ bλ

      u[:] = -FbarT' * λ
      u[:] .= g .+ u
      for e = 1:nelems
        F = locfactors[e]
        @views u[vstarts[e]:(vstarts[e+1]-1)] = F \ u[vstarts[e]:(vstarts[e+1]-1)]
        #=
        ldiv!((@view u[vstarts[e]:(vstarts[e+1]-1)]), F,
              (@view u[vstarts[e]:(vstarts[e+1]-1)]))
        =#
      end
    elseif method == :displacements
      # Solve for the discplacement directly
      M = blockdiag(ntuple(i->lop[i].M̃, length(lop))...)
      C = M - FbarT' * Diagonal(1 ./ D) * FbarT
      u[:] = C \ (g - FbarT' * (gδ ./ D))
    end

    for e = 1:nelems
      (x, y) = lop[e].coord
      JH = lop[e].JH
      #@show vex(x[:],y[:],e)' * JH * vex(x[:],y[:],e)
      @views Δ[vstarts[e]:(vstarts[e+1]-1)] = (u[vstarts[e]:(vstarts[e+1]-1)] -
                                               vex(x[:], y[:], e))
      ϵ[lvl] += (Δ[vstarts[e]:(vstarts[e+1]-1)]' * JH * Δ[vstarts[e]:(vstarts[e+1]-1)])
    end

    for f = 1:nfaces
      e1 = FToE[1, f]
      lf1 = FToLF[1, f]
      if FToB[f] == BC_JUMP_INTERFACE
        (xf, yf) = (lop[e1].facecoord[1][lf1], lop[e1].facecoord[2][lf1])
        (nx, ny) = (lop[e1].nx[lf1], lop[e1].ny[lf1])
        τex = nx .* vex_x(xf, yf, e1) + ny .* vex_y(xf, yf, e1)
        HsJ = lop[e1].Hf[lf1] * Diagonal(lop[e1].sJ[lf1])

        # Compute the shear-traction
        λrng = FToλstarts[f]:(FToλstarts[f+1]-1)
        δrng = FToδstarts[f]:(FToδstarts[f+1]-1)
        urng = vstarts[e1]:(vstarts[e1+1]-1)
        τ = computetraction(lop[e1], lf1, u[urng], λ[λrng], δ[δrng])
        Δτ = τ - τex
        τϵ[lvl] += Δτ' * HsJ * Δτ
      end
    end
    ϵ[lvl] = sqrt(ϵ[lvl])
    τϵ[lvl] = sqrt(τϵ[lvl])
    @show (lvl, ϵ[lvl], τϵ[lvl])
  end
  println((log.(ϵ[1:end-1]) - log.(ϵ[2:end])) / log(2))
  println((log.(τϵ[1:end-1]) - log.(τϵ[2:end])) / log(2))

  nothing
end


text:
using Documenter, CategoricalArrays

makedocs(
    modules = [CategoricalArrays],
    format = :html,
    sitename = "CategoricalArrays",
    pages = Any[
        "Overview" => "index.md",
        "Using CategoricalArrays" => "using.md",
        "Implementation details" => "implementation.md",
        "Index" => "functionindex.md"
        ]
    )

deploydocs(
    repo = "github.com/JuliaData/CategoricalArrays.jl.git",
    target = "build",
    julia  = "nightly",
    osname = "linux",
    deps = nothing,
    make = nothing
)

ext: jl


text:
# we write a function that computes x factorial iteratively
function simpleFactorial(x)

	# we check that x is positive
	if x < 0
		error("factorials of non-negative numbers only")
	end

	# we initialize y to 1
	y = 1

	# we loop through all integers from 1 to x and multiply them onto y
	for k = 1:x
		y *= k
	end

	# we return y
 	return y
end


text:
module SampleData

using Celeste: Model, DeterministicVI
import Celeste: Infer
import Celeste: Synthetic
import Celeste.SDSSIO: RunCamcolField, load_field_images, PlainFITSStrategy

using Distributions
using StaticArrays
import WCS, FITSIO, DataFrames

export empty_model_params, datadir,
       sample_ce, perturb_params,
       sample_star_fluxes, sample_galaxy_fluxes,
       gen_sample_star_dataset, gen_sample_galaxy_dataset,
       gen_two_body_dataset, gen_three_body_dataset, gen_n_body_dataset,
       make_elbo_args, true_star_init

const datadir = joinpath(Pkg.dir("Celeste"), "test", "data")

const sample_star_fluxes = [
    4.451805E+03,1.491065E+03,2.264545E+03,2.027004E+03,1.846822E+04]
const sample_galaxy_fluxes = [
    1.377666E+01, 5.635334E+01, 1.258656E+02,
    1.884264E+02, 2.351820E+02] * 100  # 1x wasn't bright enough

# A world coordinate system where the world and pixel coordinates are the same.
const wcs_id = WCS.WCSTransform(2,
                    cd = Float64[1 0; 0 1],
                    ctype = ["none", "none"],
                    crpix = Float64[1, 1],
                    crval = Float64[1, 1]);

# Ensure that test images are available.
wd = pwd()
cd(datadir)
run(`make RUN=3900 CAMCOL=6 FIELD=269`)
cd(wd)

const sample_rcf = RunCamcolField(3900, 6, 269)
const sample_images = load_field_images(PlainFITSStrategy(datadir),sample_rcf)

"""
Turn a images and vector of catalog entries into elbo arguments
that can be used with Celeste.
"""
function make_elbo_args(images::Vector{Image},
                        catalog::Vector{CatalogEntry};
                        active_source=-1,
                        patch_radius_pix::Float64=NaN,
                        include_kl=true)
    patches = Infer.get_sky_patches(images,
                                    catalog,
                                    radius_override_pix=patch_radius_pix)
    S = length(catalog)
    active_sources = active_source > 0 ? [active_source] :
                                          S <= 3 ? collect(1:S) : [1,2,3]
    ElboArgs(images, patches, active_sources; include_kl=include_kl)
end

function empty_model_params(S::Int)
    vp = [DeterministicVI.generic_init_source([ 0., 0. ]) for s in 1:S]
    ElboArgs(Image[],
             vp,
             Matrix{SkyPatch}(S, 0),
             collect(1:S))
end

function sample_ce(pos, is_star::Bool)
    CatalogEntry(pos, is_star, sample_star_fluxes, sample_galaxy_fluxes,
        0.1, 0.7, pi/4, 4.0)
end

# for testing away from the truth, where derivatives != 0
function perturb_params(vp)
    for vs in vp
        vs[ids.is_star] = [ 0.4, 0.6 ]
        vs[ids.pos[1]] += .8
        vs[ids.pos[2]] -= .7
        vs[ids.flux_loc] -= log(10)
        vs[ids.flux_scale] *= 25.
        vs[ids.gal_fracdev] += 0.05
        vs[ids.gal_ab] += 0.05
        vs[ids.gal_angle] += pi/10
        vs[ids.gal_scale] *= 1.2
        vs[ids.color_mean] += 0.5
        vs[ids.color_var] =  1e-1
    end
end

function gen_sample_star_dataset(; perturb=true)
    srand(1)
    images0 = deepcopy(sample_images)
    for b in 1:5
        images0[b].H, images0[b].W = 20, 23
        images0[b].wcs = wcs_id
    end
    catalog = [sample_ce([10.1, 12.2], true),]
    images = Synthetic.gen_blob(images0, catalog)
    ea = make_elbo_args(images, catalog)

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

function gen_sample_galaxy_dataset(; perturb=true, include_kl=true)
    srand(1)
    images0 = deepcopy(sample_images)
    for b in 1:5
        images0[b].H, images0[b].W = 20, 23
        images0[b].wcs = wcs_id
    end
    catalog = [sample_ce([8.5, 9.6], false),]
    images = Synthetic.gen_blob(images0, catalog)
    ea = make_elbo_args(images, catalog; include_kl=include_kl)

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

function gen_two_body_dataset(; perturb=true)
    # A small two-body dataset for quick unit testing.  These objects
    # will be too close to be identifiable.

    srand(1)
    images0 = deepcopy(sample_images)
    for b in 1:5
        images0[b].H, images0[b].W = 20, 23
        images0[b].wcs = wcs_id
    end
    catalog = [
        sample_ce([4.5, 3.6], false),
        sample_ce([10.1, 12.1], true)
    ]
    images = Synthetic.gen_blob(images0, catalog)
    ea = make_elbo_args(images, catalog)

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

function gen_three_body_dataset(; perturb=true)
    srand(1)
    images0 = deepcopy(sample_images)
    for b in 1:5
        images0[b].H, images0[b].W = 112, 238
        images0[b].wcs = wcs_id
    end
    catalog = [
        sample_ce([4.5, 3.6], false),
        sample_ce([60.1, 82.2], true),
        sample_ce([71.3, 100.4], false),
    ];
    images = Synthetic.gen_blob(images0, catalog);
    ea = make_elbo_args(images, catalog);

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

"""
Generate a large dataset with S randomly placed bodies and non-constant
background.
"""
function gen_n_body_dataset(
        S::Int; patch_pixel_radius=20., seed=NaN, perturb=true)
    if !isnan(seed)
        srand(seed)
    end

    images0 = deepcopy(sample_images)
    img_size_h = 900
    img_size_w = 1000
    for b in 1:5
        images0[b].H, images0[b].W = img_size_h, img_size_w
    end

    fluxes = [4.451805E+03,1.491065E+03,2.264545E+03,2.027004E+03,1.846822E+04]

    locations = rand(2, S) .* [img_size_h, img_size_w]
    world_locations = WCS.pix_to_world(images0[3].wcs, locations)

    catalog = CatalogEntry[CatalogEntry(world_locations[:, s], true,
            fluxes, fluxes, 0.1, 0.7, pi/4, 4.0) for s in 1:S];

    images = Synthetic.gen_blob(images0, catalog);

    # Make non-constant background.
    for b=1:5
        images[b].nelec_per_nmgy = fill(images[b].nelec_per_nmgy[1], images[b].H)
        images[b].sky = SkyIntensity(fill(images[b].sky[1,1], images[b].H, images[b].W),
                                     collect(1:images[b].H), collect(1:images[b].W),
                                     ones(images[b].H))
    end

    ea = make_elbo_args(
        images, catalog, patch_radius_pix=patch_pixel_radius)

    vp = Vector{Float64}[DeterministicVI.catalog_init_source(ce) for ce in catalog]
    if perturb
        perturb_params(vp)
    end

    ea, vp, catalog
end

function true_star_init()
    ea, vp, catalog = gen_sample_star_dataset(perturb=false)

    vp[1][ids.is_star] = [ 1.0 - 1e-4, 1e-4 ]
    vp[1][ids.flux_scale] = 1e-4
    vp[1][ids.flux_loc] = log(sample_star_fluxes[3]) - 0.5 * vp[1][ids.flux_scale]
    vp[1][ids.color_var] = 1e-4

    ea, vp, catalog
end

end # End module


text:
"""
    tex_layout(mathexpr::TeXExpr, fontset)

Recursively determine the layout of the math expression represented the given
TeXExpr for the given font set.
"""
function tex_layout(expr, fontset)
    head = expr.head
    args = [expr.args...]
    shrink = 0.6

    try
        if head == :accent
            # TODO
        elseif head == :decorated
            core, sub, super = tex_layout.(args, Ref(fontset))

            core_width = advance(core)

            return Group(
                [core, sub, super],
                Point2f0[
                    (0, 0),
                    (core_width, -0.2),
                    (core_width, xheight(core) - 0.5 * descender(super))],
                [1, shrink, shrink]
            )
        elseif head == :delimited
            elements = tex_layout.(args, Ref(fontset))
            left, content, right = elements

            height = inkheight(content)
            left_scale = max(1, height / inkheight(left))
            right_scale = max(1, height / inkheight(right))
            scales = [left_scale, 1, right_scale]

            dxs = advance.(elements) .* scales
            xs = [0, cumsum(dxs[1:end-1])...]

            # TODO Height calculation for the parenthesis looks wrong
            # TODO Check what the algorithm should be there
            # Center the delimiters in the middle of the bot and top baselines ?
            return Group(elements,
                Point2f0[
                    (xs[1], -bottominkbound(left) + bottominkbound(content)),
                    (xs[2], 0),
                    (xs[3], -bottominkbound(right) + bottominkbound(content))
                ],
                scales
            )
        elseif head == :font
            # TODO
        elseif head == :frac
            numerator = tex_layout(args[1], fontset)
            denominator = tex_layout(args[2], fontset)

            # extend fraction line by half an xheight
            xh = xheight(fontset)
            w = max(inkwidth(numerator), inkwidth(denominator)) + xh/2

            # fixed width fraction line
            lw = thickness(fontset)

            line = HLine(w, lw)
            y0 = xh/2 - lw/2

            # horizontal center align for numerator and denominator
            x1 = (w-inkwidth(numerator))/2
            x2 = (w-inkwidth(denominator))/2

            ytop    = y0 + xh/2 - bottominkbound(numerator)
            ybottom = y0 - xh/2 - topinkbound(denominator)

            return Group(
                [line, numerator, denominator],
                Point2f0[(0,y0), (x1, ytop), (x2, ybottom)]
            )
        elseif head == :function
            name = args[1]
            elements = TeXChar.(collect(name), Ref(fontset), Ref(:function))
            return horizontal_layout(elements)
        elseif head == :group || head == :expr
            elements = tex_layout.(args, Ref(fontset))
            return horizontal_layout(elements)
        elseif head == :integral
            pad = 0.2
            sub, super = tex_layout.(args[2:3], Ref(fontset))

            topint = TeXChar('⌠', fontset, :symbol, raw"\inttop")
            botint = TeXChar('⌡', fontset, :symbol, raw"\intbottom")

            top = Group([topint, super],
                Point2f0[
                    (0, 0),
                    (rightinkbound(topint) + pad, topinkbound(topint) - xheight(super))
                ],
                [1, shrink])
            bottom = Group([botint, sub],
                Point2f0[
                    (0, 0),
                    (rightinkbound(botint) + pad, bottominkbound(botint))
                ],
                [1, shrink])

            return Group(
                [top, bottom],
                Point2f0[
                    (0, xheight(fontset)/2),
                    (0, xheight(fontset)/2 - inkheight(botint) - bottominkbound(botint))
                ]
            )
        elseif head == :space
            return Space(args[1])
        elseif head == :spaced
            sym = tex_layout(args[1], fontset)
            return horizontal_layout([Space(0.2), sym, Space(0.2)])
        elseif head == :sqrt
            content = tex_layout(args[1], fontset)
            sqrt = TeXChar('√', fontset, :symbol, raw"\sqrt")

            relpad = 0.15

            h = inkheight(content)
            ypad = relpad * h
            h += 2ypad

            if h > inkheight(sqrt)
                sqrt = TeXChar('⎷', fontset, :symbol, raw"\sqrtbottom")
            end

            h = max(inkheight(sqrt), h)

            # The root symbol must be manually placed
            y0 = bottominkbound(content) - bottominkbound(sqrt) - ypad/2
            y = y0 + bottominkbound(sqrt) + h
            xpad = advance(sqrt) - inkwidth(sqrt)
            w =  inkwidth(content) + 2xpad

            lw = thickness(fontset)
            hline = HLine(w, lw)
            vline = VLine(inkheight(sqrt) - h, lw)

            return Group(
                [sqrt, vline, hline, content],
                Point2f0[
                    (0, y0),
                    (rightinkbound(sqrt) - lw/2, y),
                    (rightinkbound(sqrt) - lw/2, y - lw/2),
                    (advance(sqrt), 0)
                ]
            )
        elseif head == :symbol
            char, command = args
            return TeXChar(char, fontset, :symbol, command)
        elseif head == :underover
            core, sub, super = tex_layout.(args, Ref(fontset))

            mid = hmid(core)
            dxsub = mid - hmid(sub) * shrink
            dxsuper = mid - hmid(super) * shrink

            under_offset = bottominkbound(core) - (ascender(sub) - xheight(sub)/2) * shrink
            over_offset = topinkbound(core) - descender(super)

            # The leftmost element must have x = 0
            x0 = -min(0, dxsub, dxsuper)

            return Group(
                [core, sub, super],
                Point2f0[
                    (x0, 0),
                    (x0 + dxsub, under_offset),
                    (x0 + dxsuper, over_offset)
                ],
                [1, shrink, shrink]
            )
        end
    catch
        # TODO Better error
        rethrow()
        @error "Error while processing expr"
    end

    @error "Unsupported head $(head) in expr:\n$expr"
end

tex_layout(::Nothing, fontset) = Space(0)

function tex_layout(char::Char, fontset)
    # TODO Move that to parser ?
    if char in "0123456789"
        char_type = :digit
    elseif char in ".,:;!"
        char_type = :punctuation
    elseif char in "[]()+-*/"
        char_type = :symbol
    else
        char_type = :variable
    end
    return TeXChar(char, fontset, char_type)
end

"""
    horizontal_layout(elements)

Layout the elements horizontally, like normal text.
"""
function horizontal_layout(elements)
    dxs = advance.(elements)
    xs = [0, cumsum(dxs[1:end-1])...]

    return Group(elements, Point2f0.(xs, 0))
end

function layout_text(string, fontset)
    isempty(string) && return Space(0)

    elements = TeXChar.(collect(string), Ref(fontset), Ref(:text))
    return horizontal_layout(elements)
end

"""
    unravel(element::TeXElement, pos, scale)

Flatten the layouted TeXElement and produce a single list of base element with
their associated absolute position and scale.
"""
function unravel(group::Group, parent_pos=Point2f0(0), parent_scale=1.0f0)
    positions = [parent_pos .+ pos for pos in parent_scale .* group.positions]
    scales = group.scales .* parent_scale
    elements = []

    for (elem, pos, scale) in zip(group.elements, positions, scales)
        push!(elements, unravel(elem, pos, scale)...)
    end

    return elements
end

unravel(char::ScaledChar, pos, scale) = unravel(char.char, pos, scale*char.scale)
unravel(::Space, pos, scale) = []
unravel(element, pos, scale) = [(element, pos, scale)]

"""
    generate_tex_elements(str)

Create a list of tuple `(texelement, position, scale)` from a string
of LaTeX math mode code. The elements' positions and scales are such as to
approximatively reproduce the LaTeX output.

The elments are of one of the following types
    - `TeXChar` a (unicode) character, in a specific font.
    - `HLine` a horizontal line.
    - `VLine` a vertical line.
"""
function generate_tex_elements(str, fontset=FontSet())
    expr = texparse(str)
    layout = tex_layout(expr, fontset)
    return unravel(layout)
end

# Still hacky as hell
function generate_tex_elements(str::LaTeXString, fontset=FontSet())
    parts = String.(split(str, raw"$"))
    groups = Vector{TeXElement}(undef, length(parts))
    texts = parts[1:2:end]
    maths = parts[2:2:end]

    groups[1:2:end] = layout_text.(texts, Ref(fontset))
    groups[2:2:end] = tex_layout.(texparse.(maths), Ref(fontset))

    return unravel(horizontal_layout(groups))
end


text:
function PlotDipoles(x_dipoles,y_dipoles,Lx,Ly)

    hold on

    %Plot Dipoles
    xd = x_dipoles(:,1);
    xd2 = x_dipoles(:,2);
    yd = y_dipoles(:,1);
    yd2 = y_dipoles(:,2);
    plot([xd xd2],[yd yd2],'.g','Markersize',10)

    for jj = 1:length(x_dipoles(:,1))
      if (xd2(jj) - xd(jj) > Lx/2) && abs(yd2(jj) - yd(jj)) <= Ly/2
        plot([xd2(jj)-Lx xd(jj)],[ yd2(jj) yd(jj)],'-g','Linewidth',2)
        plot([xd2(jj)  xd(jj)+Lx],[ yd2(jj) yd(jj)],'-g','Linewidth',2)
      elseif (xd2(jj) - xd(jj) < -Lx/2) && abs(yd2(jj) - yd(jj)) <= Ly/2
        plot([xd(jj)-Lx,  xd2(jj)],[ yd(jj), yd2(jj)],'-g','Linewidth',2)
        plot([xd(jj),  xd2(jj)+Lx],[ yd(jj) yd2(jj)],'-g','Linewidth',2)
      elseif (yd2(jj) - yd(jj) > Ly/2) && abs(xd2(jj) - xd(jj)) <= Lx/2
        plot([xd2(jj) xd(jj)],[yd2(jj)-Ly yd(jj)],'-g','Linewidth',2)
        plot([xd2(jj) xd(jj)],[yd2(jj) yd(jj)+Ly],'-g','Linewidth',2)
      elseif (yd2(jj) - yd(jj) < -Ly/2) && abs(xd2(jj) - xd(jj)) <= Lx/2
        plot([xd(jj),  xd2(jj)],[ yd(jj)-Ly, yd2(jj)],'-g','Linewidth',2)
        plot([xd(jj),  xd2(jj)],[ yd(jj) yd2(jj)+Ly],'-g','Linewidth',2)
      elseif (xd2(jj) - xd(jj) > Lx/2) && (yd2(jj) - yd(jj) > Ly/2)
        plot([xd2(jj) xd(jj)+Lx],[yd2(jj) yd(jj)+Ly],'-g','Linewidth',2)
        plot([xd2(jj)-Lx xd(jj)],[yd2(jj)-Ly yd(jj)],'-g','Linewidth',2)
      elseif (xd2(jj) - xd(jj) < -Lx/2) && (yd2(jj) - yd(jj) < -Ly/2)
        plot([xd2(jj) xd(jj)-Lx],[yd2(jj) yd(jj)-Ly],'-g','Linewidth',2)
        plot([xd2(jj)+Lx xd(jj)],[yd2(jj)+Ly yd(jj)],'-g','Linewidth',2)
      elseif (xd2(jj) - xd(jj) < -Lx/2) && (yd2(jj) -yd(jj) > Ly/2)
        plot([xd2(jj) xd(jj)-Lx],[yd2(jj) yd(jj)+Ly],'-g','Linewidth',2)
        plot([xd2(jj)+Lx xd(jj)],[yd2(jj)-Ly yd(jj)],'-g','Linewidth',2)
      elseif (xd2(jj) - xd(jj) > Lx/2) && (yd2(jj) -yd(jj) < -Ly/2)
        plot([xd2(jj) xd(jj)+Lx],[yd2(jj) yd(jj)-Ly],'-g','Linewidth',2)
        plot([xd2(jj)-Lx xd(jj)],[yd2(jj)+Ly yd(jj)],'-g','Linewidth',2)
      else
        plot([ xd(jj) xd2(jj)],[ yd(jj) yd2(jj)],'-g','Linewidth',2)

      end
    end

    xlim([-Lx/2 Lx/2])
    ylim([-Ly/2 Ly/2])
    box on
    set(gca,'plotboxaspectRatio',[Lx/Ly 1 1])

return


text:
window = (x=(-1.2,1.5), y=(5.4,6.2))
xy = [(0.171146, 6.18128), (0.825792, 6.1237), (0.900563, 6.13717),
    (-0.543602, 5.82919), (-0.629489, 5.40847), (0.864488, 5.59834),
    (0.0382095, 5.83211), (-0.573306, 6.16146), (-0.118548, 5.74562),
    (0.137276, 5.72065), (-0.387608, 5.45271), (0.709483, 5.68043),
    (-0.769612, 6.06163), (0.523729, 6.07474), (0.389607, 5.90159),
    (-0.559866, 6.02546), (1.31894, 5.84089), (-1.03956, 5.54941),
    (0.765956, 5.61121), (-0.246564, 5.52229), (0.860113, 5.52206),
    (0.576003, 6.12809), (-0.179952, 5.4613), (-1.14176, 6.01474),
    (0.689314, 5.58003), (1.49035, 6.00802), (0.363703, 5.40644),
    (-0.422678, 5.50465), (0.819573, 5.98695), (-0.08627, 6.05889)]
ref = [0.0, 4.3653322423953256e-5, 0.00017530848427371737, 0.0003954137396170454,
    0.0007023905710382694, 0.0010941330363524182, 0.0015666274222397858, 0.0021217332156657775,
    0.0027613239509604037, 0.003484672246775289, 0.0042913167206244385, 0.00517892470218595,
    0.00614875508238788, 0.007202760029804356, 0.008341160543495763, 0.009564503911010491,
    0.01087409358071667, 0.012269344698312201, 0.013750792125568512, 0.015320031173485349,
    0.01697293702895064, 0.018708352664438044, 0.0205289300602286, 0.0224348705029499,
    0.024428426110674306, 0.026509091137849383, 0.028677567773440216, 0.030931700826614672,
    0.03326768504771371, 0.03568943050206452, 0.0381917913940748, 0.04077314873895044,
    0.043419613944554625, 0.046125242533531674, 0.048900285364247065, 0.05174871917671808,
    0.05467257486607324, 0.05766438503293714, 0.06072045899357337, 0.06380869016929558,
    0.06693279226313487, 0.07010646532094111, 0.07332940871453009, 0.07659523992013828,
    0.07991266880424841, 0.08328269017007739, 0.08668307271402931, 0.09012249976969666,
    0.09360888067572704, 0.0971464403870087, 0.10073595750064934, 0.10436252694482595,
    0.10802437936502629, 0.11172157658825133, 0.11545979228092684, 0.11924955042305829,
    0.12309190992784913, 0.12698781546849924, 0.13093984596269448, 0.13494881691918492,
    0.13901412290209092, 0.14313832311172603, 0.14730187769406045, 0.15151181894848642,
    0.15577136564643812, 0.16008148010916545, 0.16444270743772182, 0.16885444111056125,
    0.17332000496449818, 0.17782367112728947, 0.18233519069991422, 0.18686459427444946,
    0.19141130085872982, 0.19598651251784016, 0.2005942384331224, 0.2052296419300934,
    0.20986440774184667, 0.21451658061013812, 0.2191921993089041, 0.22389011449437113,
    0.22859796856603665, 0.233321490175584, 0.23805854733009113, 0.24281739177584483,
    0.24760222733044834, 0.25241471041712726, 0.25725635303125516, 0.2621285633981171,
    0.2670314569532928, 0.27196404263915064, 0.2769194451978795, 0.28185621923197546,
    0.28677862873022897, 0.29170394795948806, 0.296639130861576, 0.3015888196792571,
    0.30655339227950207, 0.3115365099534235, 0.31653820593547866, 0.3215489326308548,
    0.3265528852257139, 0.33156035051241517, 0.3365763311491947, 0.34160407159663797,
    0.34664654645486126, 0.3517080323709918, 0.3567913155315914, 0.3618840544208992,
    0.3669715079865965, 0.37206580814734935, 0.37716060851673616, 0.38223965941515536,
    0.3873160020065506, 0.3923963874015237, 0.3974826481353848, 0.4025759181684342,
    0.4076840068459192, 0.4128089505607022, 0.41795047702792987, 0.4231107072721887,
    0.42829029452455913, 0.43349001884597405, 0.43870952913851835, 0.4439497426336315,
    0.44921155495659426, 0.4544956322902398, 0.45980184876528274, 0.4651322243064496,
    0.47048584836835483, 0.47585767445454086, 0.48120171305696635, 0.4865118725229576,
    0.4918088221886212, 0.4970974541522134, 0.5023820886054363, 0.5076654508620471,
    0.5129484758793467, 0.5182341208828642, 0.5235235399928748, 0.5288185540103508,
    0.5341199082311753, 0.5394267209091863, 0.5447222623341462, 0.550014491055092,
    0.5553086026585501, 0.5606068831140114, 0.5659102102335019, 0.5712204773380845,
    0.5765405320214001, 0.5818682021236018, 0.5872029706731531, 0.5925370793926497,
    0.597874518789654, 0.6032179903279467, 0.6085673763674821, 0.6139262675819204,
    0.6192715296607954, 0.6245868561243617, 0.6298842253115768, 0.635157089302571,
    0.6404167200697142, 0.6456516896333067, 0.6508498986099556, 0.6560242456690262,
    0.6611374893276698, 0.6662060215622143, 0.6712418835868752, 0.6762499414534344,
    0.6812314854721027, 0.6861885008878573, 0.6911227211197069, 0.696034368600978,
    0.700926829085495, 0.7058009550031068, 0.7106248478748473, 0.7153905520579944,
    0.7201155075415395, 0.7248072954285946, 0.7294672938522346, 0.7340980395065325,
    0.7386991352512237, 0.7432582943104953, 0.7477719163148111, 0.7522419705225791,
    0.7566727187886603, 0.7610648401995753, 0.7654218756658812, 0.7697421483832405,
    0.7739856495669473, 0.778184506387009, 0.7823436834703668, 0.7864598207658212,
    0.7905345343205491, 0.7945746459809095, 0.7985832393164524, 0.8025430530865227,
    0.8064171223231932, 0.8102328402692378, 0.8139989877398341, 0.8177104836416438,
    0.8213730377071582, 0.8249878569708711, 0.8285574417979189, 0.832085355982109,
    0.8355720068056933, 0.8390199341387028, 0.8424303607083435, 0.8458029036641665,
    0.8491392263048995, 0.8524396287622842, 0.8557033083682104, 0.8589410390162966,
    0.8621536821033504, 0.865340601334038, 0.8685036668628742, 0.8716416392240589,
    0.8747529132788527, 0.8778165889508787, 0.8807810162479123, 0.8836276129663937,
    0.8864137653591784, 0.8891544555771238, 0.8918512759569056, 0.8945039460548612,
    0.8971066495898357, 0.8996630797540144, 0.9021779793636634, 0.9046640238455361,
    0.9071234433997776, 0.9095589221085179, 0.9119701262004352, 0.9143569373502427,
    0.9167205830294206, 0.9190624170457723, 0.9213817664563575, 0.9236795960482577,
    0.9259557174574908, 0.9282113217176574, 0.9304456554011953, 0.9326587737902291,
    0.9348529649639792, 0.9370262524052397, 0.9391817122886855, 0.9413185588741494,
    0.9434366007076981, 0.9455364521281238, 0.9476183397424743, 0.94968309837709,
    0.9517297171437048, 0.9537165413323675, 0.9556599963451246]
r = 0:0.001:0.25
pp = PointPattern(xy, window)
est = Fest(pp, r)
@test maximum(abs, (est - ref)) < 0.01
est = Fest(pp, r, 1000)
@test maximum(abs, (est - ref)) < 0.001


text:
Base.@pure valadd(::Val{N1},::Val{N2}) where {N1,N2} = Val{N1+N2}()
Base.@pure valsub(::Val{N1},::Val{N2}) where {N1,N2} = N1 > N2 ? Val{N1-N2}() : Val{0}()

tpermute(t::NTuple{N,T}, p::NTuple{N,Int}) where {T,N} = tselect(t, p)
tselect(t::Tuple, ind::NTuple{N,Int}) where {N} = ntuple(n->t[ind[n]], Val(N))

@inline argtail2(a, b, c...) = c
@inline tail2(x::Tuple) = argtail2(x...)

function permutation2swaps(perm)
    p = collect(perm)
    swaps = Vector{Int}()
    N = length(p)
    for k = 1:N-1
        append!(swaps, p[k]-1:-1:k)
        for l = k+1:N
            if p[l] < p[k]
                p[l] += 1
            end
        end
        p[k] = k
    end
    return swaps
end


text:
using Test, Random, AutoGrad, Statistics
using Knet.Ops21: pool
using Knet.KnetArrays: KnetArray
using CUDA: CUDA, CuArray

@testset "ops21/pool" begin

    function pooltest(
        ;
        atype = Array,
        nd = 4,
        op = maximum,
        window = 2,
        padding = 0,
        stride = window,
        propagateNaN = false,
        includePadding = false,
        channelmajor = false,
        alpha = 1,
        gcheck_test = true,
        approx_test = true,
    )
        xd = rand(5:10, nd)
        x = randn(xd...)
        kw = (; op, window, padding, stride, propagateNaN, includePadding, channelmajor, alpha)
        ax = atype(x)
        px = Param(ax)
        r1 = ((isa(ax,Array) || !approx_test) ? true :
              isapprox(pool(x; kw...), Array(pool(ax; kw...))))
        atol,rtol = (eltype(px) == Float64 ? (0.01, 0.05) : (0.1, 0.5))
        r2 = !gcheck_test ? true : @gcheck pool(px; kw...) (;rtol,atol)
        r1 && r2
    end

    @test pooltest(; )
    @test pooltest(; nd=3)
    @test pooltest(; nd=5)
    @test pooltest(; op=mean)
    @test pooltest(; window=3)
    @test pooltest(; window=typemax(Int))
    @test pooltest(; padding=1)
    @test pooltest(; stride=3)
    @test pooltest(; propagateNaN=true)
    @test pooltest(; includePadding=true, op=mean, padding=1)
    @test pooltest(; includePadding=false, op=mean, padding=1)
    @test pooltest(; channelmajor=true)
    @test pooltest(; alpha=2)

    if CUDA.functional()

        @test pooltest(; atype=CuArray)
        @test pooltest(; atype=CuArray, nd=3)
        @test pooltest(; atype=CuArray, nd=5)
        @test pooltest(; atype=CuArray, op=mean)
        @test pooltest(; atype=CuArray, window=3)
        @test pooltest(; atype=CuArray, window=typemax(Int))
        @test pooltest(; atype=CuArray, padding=1)
        @test pooltest(; atype=CuArray, stride=3)
        @test pooltest(; atype=CuArray, propagateNaN=true)
        @test pooltest(; atype=CuArray, includePadding=true, op=mean, padding=1)
        @test_skip pooltest(; atype=CuArray, includePadding=false, op=mean, padding=1)
        @test pooltest(; atype=CuArray, channelmajor=true)
        @test pooltest(; atype=CuArray, alpha=2)

        @test pooltest(; atype=KnetArray)
        @test pooltest(; atype=KnetArray, nd=3)
        @test pooltest(; atype=KnetArray, nd=5)
        @test pooltest(; atype=KnetArray, op=mean)
        @test pooltest(; atype=KnetArray, window=3)
        @test pooltest(; atype=KnetArray, window=typemax(Int))
        @test pooltest(; atype=KnetArray, padding=1)
        @test pooltest(; atype=KnetArray, stride=3)
        @test pooltest(; atype=KnetArray, propagateNaN=true)
        @test pooltest(; atype=KnetArray, includePadding=true, op=mean, padding=1)
        @test_skip pooltest(; atype=KnetArray, includePadding=false, op=mean, padding=1)
        @test pooltest(; atype=KnetArray, channelmajor=true)
        @test pooltest(; atype=KnetArray, alpha=2)

    end
end


text:
# Practical Machine learning
# Decision Trees based learning - Random Forests example
# Chapter 5

typealias TabularData Union{AbstractMatrix,DataFrame}

type Example{T<:TabularData}
    x::T  # tabular data
    y::AbstractVector
    n_labels::Int
    n_features::Int
    sample_weight::Vector{Float64}

    function Example(x::T, y::AbstractVector, sample_weight::Vector{Float64})
        n_labels = length(unique(y))
        n_features = size(x, 2)
        new(x, y, n_labels, n_features, sample_weight)
    end

    Example(x::T, y::AbstractVector) = Example{T}(x, y, ones(Float64, size(x, 1)))
end


text:
# have to import this or @deprecate doesn't work
import Base: getindex, setindex!, delete!, haskey, get, copy!, getproperty, setproperty!

@deprecate(ImageMeta(data::AbstractArray{T,N}, props::AbstractDict{String,Any}) where {T,N},
           ImageMeta(data, to_dict(props)))

@deprecate(ImageMeta(data::AbstractArray{T,N}, props::Dict{String,Any}) where {T,N},
           ImageMeta(data, to_dict(props)))

@deprecate(ImageMeta(data::AbstractArray, props::Dict{<:AbstractString}),
           ImageMeta(data, to_dict(props)))

@deprecate(copy!(imgdest::ImageMeta, imgsrc::ImageMeta, prop1::AbstractString, props::AbstractString...),
           copy!(imgdest, imgsrc, Symbol(prop1), Symbol.(props)...))

@deprecate(delete!(img::ImageMeta, propname::AbstractString),
           delete!(img, Symbol(propname)))

@deprecate(haskey(img::ImageMeta, k::AbstractString),
           hasproperty(img, Symbol(k)))

@deprecate(hasproperty(img::ImageMeta, k::AbstractString),
           hasproperty(img, Symbol(k)))

@deprecate(get(img::ImageMeta, k::AbstractString, default),
           get(img, Symbol(k), default))

@deprecate(getindex(img::ImageMeta, propname::AbstractString),
           getproperty(img, Symbol(propname)))

@deprecate(getproperty(img::ImageMeta, propname::AbstractString),
           getproperty(img, Symbol(propname)))

@deprecate(setproperty!(img::ImageMeta, propname::AbstractString, val),
           setproperty!(img::ImageMeta, Symbol(propname), val))

@deprecate(setindex!(img::ImageMeta, X, propname::AbstractString),
           setproperty!(img, Symbol(propname), X))

@deprecate(data(img::ImageMeta), arraydata(img))

function to_dict(dold::AbstractDict{String})
    dnew = Dict{Symbol,Any}()
    for (k, v) in dold
        dnew[Symbol(k)] = v
    end
    return dnew
end


text:
#
# clock and event routines
#

"""
    setUnit!(sim::Clock, new::FreeUnits)

set a clock to a new time unit in `Unitful`. If necessary convert
current clock times to the new unit.

# Arguments
- `sim::Clock`
- `new::FreeUnits`: new is one of `ms`, `s`, `minute` or `hr` or another Unitful
    `Time` unit.

# Examples
```jldoctest
julia> using Simulate, Unitful

julia> import Unitful: Time, s, minute, hr

julia> c = Clock(t0=60) # setup a new clock with t0=60
Clock: state=Simulate.Undefined(), time=60.0, unit=, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=0.0
julia> tau(c) # current time is 60.0 NoUnits
60.0
julia> setUnit!(c, s)  # set clock unit to Unitful.s
60.0 s
julia> tau(c) # current time is now 60.0 s
60.0 s
julia> setUnit!(c, minute)  # set clock unit to Unitful.minute
1.0 minute
julia> tau(c) # current time is now 1.0 minute
1.0 minute
julia> typeof(tau(c))  # tau(c) now returns a time Quantity ...
Quantity{Float64,𝐓,Unitful.FreeUnits{(minute,),𝐓,nothing}}
julia> isa(tau(c), Time)
true
julia> uconvert(s, tau(c)) # ... which can be converted to other time units
60.0 s
julia> tau(c).val  # it has a value of 1.0
1.0
julia> c.time  # internal clock time is set to 1.0 (is a Float64)
1.0
julia> c.unit  # internal clock unit is set to Unitful.minute
minute
```
"""
function setUnit!(sim::Clock, new::FreeUnits)
    if isa(1new, Time)
        if sim.unit == new
            println("clock is already set to $new")
        elseif sim.unit == NoUnits
            sim.unit = new
        else
            old = sim.unit
            sim.unit = new
            fac = uconvert(new, 1*old).val
            sim.time *= fac
            sim.end_time *= fac
            sim.tev *= fac
            sim.Δt *= fac
            sim.tsa *= fac
        end
    else
        sim.unit = NoUnits
    end
    tau(sim)
end

"""
```
𝐶
Clk
```
`𝐶` (𝐶 = \\itC+tab) or `Clk` is the central simulation clock. Since most
simulations work with **one** time, they should use 𝐶 for time keeping.

# Examples

```jldoctest
julia> using Simulate

julia> reset!(𝐶)
"clock reset to t₀=0.0, sampling rate Δt=0.0."
julia> 𝐶  # central clock
Clock: state=Simulate.Idle(), time=0.0, unit=, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=0.0
julia> Clk  # alias
Clock: state=Simulate.Idle(), time=0.0, unit=, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=0.0
julia> 𝐶.time
0.0
```
"""
const 𝐶 = Clk = Clock()

"""
```
tau(sim::Clock=𝐶)
τ(sim::Clock=𝐶)
```
Return the current simulation time (τ = \\tau+tab).

# Examples

```jldoctest
julia> using Simulate

julia> reset!(𝐶)
"clock reset to t₀=0.0, sampling rate Δt=0.0."
julia> tau() # gives the central time
0.0
julia> τ() # alias, gives the central time
0.0
```
"""
tau(sim::Clock=𝐶) = sim.time*sim.unit
const τ = tau

"""
```
sync!(sim::Clock, to::Clock=𝐶)
```
Force a synchronization of two clocks. Change all registered times of
`sim` accordingly. Convert or force sim.unit to to.unit.
"""
function sync!(sim::Clock, to::Clock=𝐶)
    if (sim.unit == NoUnits) | (sim.unit == to.unit)
        fac = 1
    elseif to.unit == NoUnits
        println(stderr, "Warning: deleted time unit without conversion")
        fac = 1
    else
        fac = uconvert(to.unit, 1sim.unit).val
    end
    Δt = to.time - sim.time*fac
    sim.time = sim.time*fac + Δt
    sim.unit = to.unit
    sim.tsa  = sim.tsa*fac + Δt
    sim.tev  = sim.tev*fac + Δt
    sim.end_time = sim.end_time*fac + Δt
    sim.Δt = to.Δt
    evq = PriorityQueue{SimEvent,Float64}()
    for (ev, t) ∈ pairs(sim.events)
        evq[ev] = t*fac + Δt
    end
    sim.events = evq
    sim
end

"""
```
reset!(sim::Clock, Δt::Number=0; t0::Number=0, hard::Bool=true, unit=NoUnits)
```
reset a clock

# Arguments
- `sim::Clock`
- `Δt::Number=0`: time increment
- `t0::Float64=0` or `t0::Time`: start time
- `hard::Bool=true`: time is reset, all scheduled events and sampling are
    deleted. If hard=false, then only time is reset, event and
    sampling times are adjusted accordingly.
- `unit=NoUnits`: the Time unit for the clock after reset. If a `Δt::Time` is
    given, its Time unit goes into the clock Time unit. If only t0::Time is given,
    its Time unit goes into the clock time unit.

# Examples

```jldoctest
julia> using Simulate, Unitful

julia> import Unitful: s

julia> c = Clock(1s, t0=60s)
Clock: state=Simulate.Undefined(), time=60.0, unit=s, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=1.0
julia> reset!(c)
"clock reset to t₀=0.0, sampling rate Δt=0.0."
julia> c
Clock: state=Simulate.Idle(), time=0.0, unit=, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=0.0
```
"""
function reset!(sim::Clock, Δt::Number=0;
                t0::Number=0, hard::Bool=true, unit=NoUnits)
    if  isa(1unit, Time)
        Δt = isa(Δt, Time) ? uconvert(unit, Δt).val : Δt
        t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
    elseif isa(Δt, Time)
        unit = Unitful.unit(1Δt)
        Δt = Δt.val
        t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
    elseif isa(t0, Time)
        unit = Unitful.unit(t0)
        t0 = t0.val
    else
        nothing
    end
    if hard
        sim.state = Idle()
        sim.time = t0
        sim.unit = unit
        sim.tsa = t0
        sim.tev = t0
        sim.end_time = t0
        sim.evcount = 0
        sim.scount = 0
        sim.Δt = Δt
        sim.events = PriorityQueue{SimEvent,Float64}()
        sim.cevents = SimCond[]
        sim.processes = Dict{Any, SimProcess}()
        sim.sexpr = Sample[]
    else
        sync!(sim, Clock(Δt, t0=t0, unit=unit))
    end
    "clock reset to t₀=$(float(t0*unit)), sampling rate Δt=$(float(Δt*unit))."
end

"""
    nextevent(sim::Clock)

Return the next scheduled event.
"""
nextevent(sim::Clock) = peek(sim.events)[1]

"""
    nextevtime(sim::Clock)

Return the internal time (unitless) of next scheduled event.
"""
nextevtime(sim::Clock) = peek(sim.events)[2]

"""
    simExec(ex::Union{SimExpr, Array{SimExpr,1}}, m::Module=Main)

evaluate the event expressions or SimFunctions.

# Return

the evaluated value or a tuple of evaluated values
"""
function simExec(ex::Union{SimExpr, Array{SimExpr,1}}, m::Module=Main)

    function sexec(x::SimExpr)
        if x isa SimFunction
            ret = x.func(x.arg...; x.kw...)
        else
            return Core.eval(m, x)
        end
        yield()  # to an eventually triggered process
        return ret
    end

    return ex isa SimExpr ? sexec(ex) : Tuple([sexec(x) for x in ex])
end

"""
    checktime(sim::Clock, t::Number)::Float64

check `t` given according to clock settings and return value
"""
function checktime(sim::Clock, t::Number)::Float64
    if isa(t, Real)
        return t
    else
        if sim.unit == NoUnits
            println(stderr, "Warning: clock has no time unit, ignoring units")
            return t.val
        else
            return uconvert(sim.unit, t).val
        end
    end
end

"""
```
event!(sim::Clock, ex::Union{SimExpr, Array, Tuple}, t::Number; scope::Module=Main, cycle::Number=0.0)::Float64
event!(ex::Union{SimExpr, Array, Tuple}, t::Number; scope::Module=Main, cycle::Number=0.0)
```
Schedule an event for a given simulation time.

# Arguments
- `sim::Clock`: simulation clock, if no clock is given, the event goes to 𝐶,
- `ex::{SimExpr, Array, Tuple}`: an expression or SimFunction or an array or tuple of them,
- `t::Float64` or `t::Time`: simulation time,
- `scope::Module=Main`: scope for the expressions to be evaluated
- `cycle::Float64=0.0`: repeat cycle time for an event

# returns
Scheduled internal simulation time (unitless) for that event.
May return a time `> t` from repeated applications of `nextfloat(t)`
if there are events scheduled for `t`.

# Examples
```jldoctest
julia> using Simulate, Unitful

julia> import Unitful: s, minute, hr

julia> myfunc(a, b) = a+b
myfunc (generic function with 1 method)
julia> event!(𝐶, SimFunction(myfunc, 1, 2), 1) # a 1st event
1.0
julia> event!(𝐶, SimFunction(myfunc, 2, 3), 1) #  a 2nd event to the same time
1.0000000000000002
julia> event!(𝐶, SimFunction(myfunc, 3, 4), 1s)
Warning: clock has no time unit, ignoring units
1.0000000000000004
julia> setUnit!(𝐶, s)
0.0 s
julia> event!(𝐶, SimFunction(myfunc, 4, 5), 1minute)
60.0
```
"""
function event!(sim::Clock, ex::Union{SimExpr, Array, Tuple}, t::Number;
                scope::Module=Main, cycle::Number=0.0)::Float64
    t = checktime(sim, t)
    cycle = checktime(sim, cycle)
    while any(i->i==t, values(sim.events)) # in case an event at that time exists
        t = nextfloat(float(t))                  # increment scheduled time
    end
    ev = SimEvent(sconvert(ex), scope, t, cycle)
    sim.events[ev] = t
    return t
end
event!( ex::Union{SimExpr, Array, Tuple}, t::Number; scope::Module=Main, cycle::Number=0.0) =
            event!(𝐶, ex, t, scope=scope, cycle=cycle)

"""
```
event!(sim::Clock, ex::Union{SimExpr, Array, Tuple}, T::Timing, t::Number; scope::Module=Main)::Float64
event!(ex::Union{SimExpr, Array, Tuple}, T::Timing, t::Number; scope::Module=Main)
```
Schedule a timed event, that is an event with a timing.

# Arguments
- `sim::Clock`: simulation clock, if no clock is given, the event goes to 𝐶,
- `ex::{SimExpr, Array, Tuple}`: an expression or SimFunction or an array or tuple of them,
- `T::Timing`: a timing, `at`, `after` or `every` (`before` behaves like `at`),
- `t::Float64` or `t::Time`: simulation time,
- `scope::Module=Main`: scope for the expressions to be evaluated

# returns
Scheduled internal simulation time (unitless) for that event.

# Examples
```jldoctest
julia> using Simulate, Unitful

julia> import Unitful: s, minute, hr

julia> setUnit!(𝐶, s)
0.0 s
julia> myfunc(a, b) = a+b
myfunc (generic function with 1 method)
julia> event!(𝐶, SimFunction(myfunc, 5, 6), after, 1hr)
3600.0
```
"""
function event!(sim::Clock, ex::Union{SimExpr, Array, Tuple}, T::Timing, t::Number;
                scope::Module=Main)
    t = checktime(sim, t)
    if T == after
        event!(sim, sconvert(ex), t + sim.time, scope=scope)
    elseif T == every
        event!(sim, sconvert(ex), sim.time, scope=scope, cycle=t)
    else
        event!(sim, sconvert(ex), t, scope=scope)
    end
end
event!( ex::Union{SimExpr, Array, Tuple}, T::Timing, t::Number; scope::Module=Main) =
            event!(𝐶, ex, T, t; scope=scope)

"calculate the scale form a given number"
function scale(n::Number)
    if n > 0
        i = 1.0
        while !(10^i ≤ n < 10^(i+1))
            n < 10^i ? i -= 1 : i += 1
        end
        return 10^i
    else
        return 1
    end
end

"""
```
event!(sim::Clock, ex::Union{SimExpr, Array, Tuple}, cond::Union{SimExpr, Array, Tuple}; scope::Module=Main):
event!(ex::Union{SimExpr, Array, Tuple}, cond::Union{SimExpr, Array, Tuple}; scope::Module=Main)
```
Schedule a conditional event.

It is executed immediately if the conditions are met, else the condition is
checked at each clock tick Δt. After a conditional event is triggered, it is
removed from the clock. If no sampling rate Δt is setup, a default
sampling rate is setup depending on the scale of the remaining simulation time
``Δt = scale(t_r)/100`` or ``0.01`` if ``t_r = 0``.

# Arguments
- `sim::Clock`: simulation clock, if no clock is given, the event goes to 𝐶,
- `ex::{SimExpr, Array, Tuple}`: an expression or SimFunction or an array or tuple of them,
- `cond::{SimExpr, Array, Tuple}`: a condition is an expression or SimFunction
    or an array or tuple of them. It is true if all expressions or SimFunctions
    therein return true.
- `scope::Module=Main`: scope for the expressions to be evaluated
- `cycle::Float64=0.0`: repeat cycle time for an event

# returns
current simulation time `tau(sim)`.

# Examples
```jldoctest
julia> using Simulate, Unitful

julia> import Unitful: s, minute, hr

```
"""
function event!(sim::Clock, ex::Union{SimExpr, Array, Tuple},
                cond::Union{SimExpr, Array, Tuple}; scope::Module=Main)
    if all(simExec(sconvert(cond)))   # all conditions met
        simExec(sconvert(ex))         # execute immediately
    else
        sim.Δt == 0 ? sim.Δt = scale(sim.end_time - sim.time)/100 : nothing
        push!(sim.cevents, SimCond(sconvert(cond), sconvert(ex), scope))
    end
    return tau(sim)
end
event!( ex::Union{SimExpr, Array, Tuple}, cond::Union{SimExpr, Array, Tuple};
        scope::Module=Main) = event!(𝐶, ex, cond, scope=scope)

"""
```
sample_time!(sim::Clock, Δt::Number)
sample_time!(Δt::Number)
```
set the clock's sample rate starting from now (`tau(sim)`).

# Arguments
- `sim::Clock`: if no clock is given, set the sample rate on 𝐶,
- `Δt::Number`: sample rate, time interval for sampling
"""
function sample_time!(sim::Clock, Δt::Number)
    sim.Δt = checktime(sim, Δt)
    sim.tsa = sim.time + sim.Δt
end
sample_time!(Δt::Number) = sample_time!(𝐶, Δt)

"""
```
sample!(sim::Clock, ex::Union{Expr, SimFunction}; scope::Module=Main)
sample!(ex::Union{Expr, SimFunction}; scope::Module=Main)
```
enqueue an expression for sampling.
# Arguments
- `sim::Clock`: if no clock is given, it samples on 𝐶,
- `ex::Union{Expr, SimFunction}`: an expression or function,
- `Δt::Number=sim.Δt`: set the clock's sampling rate, if no Δt is given, it takes
    the current sampling rate, if that is 0, it calculates one,
- `scope::Module=Main`: optional, an evaluation scope for a given expression.
"""
function sample!(sim::Clock, ex::Union{Expr, SimFunction}, Δt::Number=sim.Δt;
                 scope::Module=Main)
    sim.Δt = Δt == 0 ? scale(sim.end_time - sim.time)/100 : Δt
    push!(sim.sexpr, Sample(ex, scope))
end
sample!(ex::Union{Expr, SimFunction}, Δt::Number=𝐶.Δt; scope::Module=Main) =
    sample!(𝐶, ex, Δt, scope=scope)

"""
    step!(sim::Clock, ::Undefined, ::Init)

initialize a clock.
"""
function step!(sim::Clock, ::Undefined, ::Init)
    sim.state = Idle()
end

"""
    step!(sim::Clock, ::Undefined, σ::Union{Step,Run})

if uninitialized, initialize and then Step or Run.
"""
function step!(sim::Clock, ::Undefined, σ::Union{Step,Run})
    step!(sim, sim.state, Init(0))
    step!(sim, sim.state, σ)
end

"""
    setTimes(sim::Clock)

set clock times for next event or sampling action. The internal clock times
`sim.tev` and `sim.tsa` must always be set to be at least `sim.time`.
"""
function setTimes(sim::Clock)
    if length(sim.events) ≥ 1
        sim.tev = nextevtime(sim)
        sim.tsa = sim.Δt > 0 ? sim.time + sim.Δt : sim.time
    else
        sim.tsa = sim.Δt > 0 ? sim.time + sim.Δt : sim.time
        sim.tev = sim.tsa
    end
end

"""
    step!(sim::Clock, ::Union{Idle,Busy,Halted}, ::Step)

step forward to next tick or scheduled event.

At a tick evaluate (1) all sampling expressions, (2) all conditional events, then
(3) if an event is encountered, evaluate the event expression.

The internal clock times `sim.tev` and `sim.tsa` must always be set to be
at least `sim.time`.
"""
function step!(sim::Clock, ::Union{Idle,Busy,Halted}, ::Step)

    function exec_next_event()
        lock(sim.lock)
        sim.time = sim.tev
        unlock(sim.lock)
        ev = dequeue!(sim.events)
        simExec(ev.ex, ev.scope)
        sim.evcount += 1
        if ev.Δt > 0.0  # schedule repeat event
            event!(sim, ev.ex, sim.time + ev.Δt, scope=ev.scope, cycle=ev.Δt)
        end
        sim.tev = length(sim.events) ≥ 1 ? nextevtime(sim) : sim.time
    end

    function exec_next_tick()
        lock(sim.lock)
        sim.time = sim.tsa
        unlock(sim.lock)
        for s ∈ sim.sexpr
            simExec(s.ex, s.scope)
        end
        cond = [all(simExec(c.cond)) for c in sim.cevents]
        if any(cond)
            ex = sim.cevents[cond]
            sim.cevents = sim.cevents[.!cond]
            [simExec(c.ex) for c in ex]
        end
        sim.scount +=1
    end

    if (sim.tev ≤ sim.time) && (length(sim.events) ≥ 1)
        sim.tev = nextevtime(sim)
    end

    if (length(sim.events) ≥ 1) | (sim.Δt > 0)
        if length(sim.events) ≥ 1
            if (sim.Δt > 0)
                if sim.tsa <= sim.tev
                    exec_next_tick()
                    if sim.tsa == sim.tev
                        exec_next_event()
                    end
                    sim.tsa += sim.Δt
                else
                    exec_next_event()
                end
            else
                exec_next_event()
                sim.tsa = sim.time
            end
        else
            exec_next_tick()
            sim.tsa += sim.Δt
            sim.tev = sim.time
        end
    else
        println(stderr, "step!: nothing to evaluate")
    end
#    length(sim.processes) == 0 || yield() # let processes run
    for i in 1:length(sim.processes) # brute force yield
        yield()
    end
end

"""
    step!(sim::Clock, ::Idle, σ::Run)

Run a simulation for a given duration.

The duration is given with `Run(duration)`. Call scheduled events and evaluate
sampling expressions at each tick in that timeframe.
"""
function step!(sim::Clock, ::Idle, σ::Run)
    length(sim.processes) > 0 ? sleep(0.01) : nothing  # let processes startup
    sim.end_time = sim.time + σ.duration
    sim.evcount = 0
    sim.scount = 0
    sim.state = Busy()
    setTimes(sim)
    while any(i->(sim.time < i ≤ sim.end_time), (sim.tsa, sim.tev))
        step!(sim, sim.state, Step())
        if sim.state == Halted()
            return
        end
    end
    tend = sim.end_time

    # catch remaining events
    while (length(sim.events) ≥ 1) && (sim.tev ≤ tend + Base.eps(tend)*10)
        step!(sim, sim.state, Step())
        tend = nextfloat(tend)
    end

    sim.time = sim.end_time
    sim.state = Idle()
    sleep(0.1)
    "run! finished with $(sim.evcount) clock events, $(sim.scount) sample steps, simulation time: $(sim.time)"
end

"""
    step!(sim::Clock, ::Busy, ::Stop)

Stop the clock.
"""
function step!(sim::Clock, ::Busy, ::Stop)
    sim.state = Halted()
    "Halted after $(sim.evcount) events, simulation time: $(sim.time)"
end

"""
    step!(sim::Clock, ::Halted, ::Resume)

Resume a halted clock.
"""
function step!(sim::Clock, ::Halted, ::Resume)
    sim.state = Idle()
    step!(sim, sim.state, Run(sim.end_time - sim.time))
end

"""
    step!(sim::Clock, q::SState, σ::SEvent)

catch all step!-function.
"""
function step!(sim::Clock, q::SState, σ::SEvent)
    println(stderr, "Warning: undefined transition ",
            "$(typeof(sim)), ::$(typeof(q)), ::$(typeof(σ)))\n",
            "maybe, you should reset! the clock!")
end

"""
    run!(sim::Clock, duration::Number)

Run a simulation for a given duration. Call scheduled events and evaluate
sampling expressions at each tick in that timeframe.
"""
run!(sim::Clock, duration::Number) =
                        step!(sim, sim.state, Run(checktime(sim, duration)))

"""
    incr!(sim::Clock)

Take one simulation step, execute the next tick or event.
"""
incr!(sim::Clock) = step!(sim, sim.state, Step())

"""
    stop!(sim::Clock)

Stop a running simulation.
"""
stop!(sim::Clock) = step!(sim, sim.state, Stop())

"""
    resume!(sim::Clock)

Resume a halted simulation.
"""
resume!(sim::Clock) = step!(sim, sim.state, Resume())

"""
    init!(sim::Clock)

initialize a clock.
"""
init!(sim::Clock) = step!(sim, sim.state, Init(""))


text:
# plotting the results of the experiments

import JLD2
import FileIO: load
import BenchmarkTools
using PGFPlotsX, LinearAlgebra, Statistics, DelimitedFiles
using BSE_k_ISDF

push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usepgfplotslibrary{colorbrewer}")

fig_theme = @pgf {
    "cycle list/Dark2-8",
    mark_options = "solid",
    line_width = "1.5pt",
    grid = "major"}

# %% plot optical absorption spectrum

#diamond
Ha_to_eV = 27.211396641308
Erange, absorption_reference = load("diamond/131313_20/optical_absorption_reference.jld2", "Erange", "absorption")

Z_tol = 0.1
N_iter = 150
absorption_isdf = load("diamond/131313_20/optical_absorption_$(Z_tol)_$(N_iter).jld2", "absorption")

fig = @pgf Axis(
    {fig_theme...,
    "no markers",
    title = "Optical Absorption of Diamond",
    xlabel = "\$\\omega\$ [eV]", ylabel = "\$\\varepsilon_2\$",
    xmin = 0, xmax = 25, ymin = 0, ymax = 38,
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true},
    PlotInc({"dashed", line_width="2.5pt"},
        Table(; x = Ha_to_eV .* Erange, y = absorption_reference[:, 1])),
    LegendEntry("Reference"),
    PlotInc({},
        Table(; x = Ha_to_eV .* Erange, y = absorption_isdf)),
    LegendEntry("ISDF"),
)

pgfsave("diamond/131313_20/diamond_optical_absorption_spectrum.tex", fig, include_preamble = false)
pgfsave("diamond/131313_20/diamond_optical_absorption_spectrum.pdf", fig, include_preamble = false)

# graphene
Ha_to_eV = 27.211396641308
Erange, absorption_reference = load("graphene/optical_absorption_reference.jld2", "Erange", "absorption")

Z_tol = 0.1
N_iter = 100
absorption_isdf = load("graphene/optical_absorption_$(Z_tol)_$(N_iter).jld2", "absorption")

fig = @pgf Axis(
    {fig_theme...,
    "no markers",
    title = "Optical Absorption of Graphene",
    xlabel = "\$\\omega\$ [eV]", ylabel = "\$\\varepsilon_2\$",
    xmin = 1, xmax = 5, ymin = 0, ymax = 15,
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true},
    PlotInc({"dashed", line_width="2.5pt"},
        Table(; x = Ha_to_eV .* Erange, y = clamp.(absorption_reference[:, 1], 0, 20))),
    LegendEntry("Reference"),
    PlotInc({},
        Table(; x = Ha_to_eV .* Erange, y = clamp.(absorption_isdf, 0, 20))),
    LegendEntry("ISDF"),
)

pgfsave("graphene/graphene_optical_absorption_spectrum.tex", fig, include_preamble = false)
pgfsave("graphene/graphene_optical_absorption_spectrum.pdf", fig, include_preamble = false)

# %% plot optical absorption spctrum for different N_μ

example_path = "diamond/131313_20/"

Ha_to_eV = 27.211396641308

Erange, absorption_reference = load(example_path * "/optical_absorption_reference.jld2", "Erange", "absorption")

N_iter = 150
Z_tol_vec = [0.5, 0.1]
optical_absorption_lanc = []
for Z_tol in Z_tol_vec
    absorption = load(example_path * "/optical_absorption_$(Z_tol)_$(N_iter).jld2", "absorption")

    push!(optical_absorption_lanc, absorption)
end

fig = @pgf Axis(
    {fig_theme...,
    title = "Optical Absorption of Diamond",
    xlabel = "\$\\omega\$ [eV]", ylabel = "\$\\varepsilon_2\$",
    xmin = 0, xmax = 25, ymin = 0, ymax = 38,
    legend_pos = "north east",
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true},
    PlotInc({"dashed", line_width="3pt"},
        Table(; x = Ha_to_eV .* Erange, y = absorption_reference[:, 1])),
    LegendEntry("Reference"),
    PlotInc({},
        Table(; x = Ha_to_eV .* Erange, y = optical_absorption_lanc[1])),
    LegendEntry("\$Z_{tol} = 0.5\$"),
    PlotInc({},
        Table(; x = Ha_to_eV .* Erange, y = optical_absorption_lanc[2])),
    LegendEntry("\$Z_{tol} = 0.1\$")
)

pgfsave("diamond/131313_20/diamond_optical_absorption_spectrum_isdf.tex", fig, include_preamble = false)
pgfsave("diamond/131313_20/diamond_optical_absorption_spectrum_isdf.pdf", fig, include_preamble = false)

# %% plotting of error in Z

N_μ_vec, errors_Z_vv, errors_Z_cc, errors_Z_vc =  load("diamond/131313_20/errors_Z.jld2", "N_μ_vec", "errors_Z_vv", "errors_Z_cc", "errors_Z_vc")

fig = @pgf Axis(
    {fig_theme...,
    title = "Error in ISDF", xlabel = "\$N_\\mu^{ij}\$",
    legend_entries = {"\$Z^{cc}\$", "\$Z^{vv}\$", "\$Z^{vc}\$"},
    xmin = N_μ_vec[1], xmax = N_μ_vec[end], ymin = 1e-4, ymax = 1e0,
    ymode = "log",
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true},
    PlotInc({"solid", mark = "*"},
        Table(; x = N_μ_vec, y = errors_Z_cc)),
    PlotInc({"dashed", mark = "cube*"},
        Table(; x = N_μ_vec, y = errors_Z_vv)),
    PlotInc({"dotted", mark = "triangle*"},
        Table(; x = N_μ_vec, y = errors_Z_vc))
)

pgfsave("diamond/131313_20/diamond_errors_isdf.tex", fig, include_preamble = false)
pgfsave("diamond/131313_20/diamond_errors_isdf.pdf", fig, include_preamble = false)

# %% plot benchmark results for diamond

N_μ_vv = 70
N_μ_cc = 220
N_μ_vc = 100

N_ks_vec, setup_times, evaluation_times =  load("diamond/benchmark_$(N_μ_vv)_$(N_μ_cc)_$(N_μ_vc).jld2", "N_ks_vec", "setup_times", "evaluation_times")

fig = @pgf LogLogAxis(
    {fig_theme...,
    title = "Run Time Scaling", xlabel = "\$N_k\$", ylabel = "run time [s]",
    legend_pos = "south east",
    legend_style = "{at={(0.6,0.2)},anchor=west}",
    "width = 0.4\\textwidth",
    "height = 0.3\\textwidth",
    scale_only_axis = true},
    PlotInc({mark = "*"},
        Table(; x = prod.(N_ks_vec), y = setup_times)),
    LegendEntry("initial setup"),
    PlotInc({mark = "square*"},
        Table(; x = prod.(N_ks_vec), y = evaluation_times)),
    LegendEntry("matrix vector product"),
    PlotInc({"dashed"},
        Table(; x = prod.(N_ks_vec), y = 2e-2 * prod.(N_ks_vec))),
    LegendEntry("\$O(N_k)\$")
)

pgfsave("diamond/diamond_timings_k.tex", fig, include_preamble = false)
pgfsave("diamond/diamond_timings_k.pdf", fig, include_preamble = false)


text:
# ----------------------------------------------------------------------------------- #
#  estimate_steady_state: estimates the steady-state vector for the dynamic system
#
# Type: GRN-JULIA
# Version: 1.0
#
# Input arguments:
# problem  - dictionary instance (holds model parameters)
#
# Return arguments:
# XSS - steady state vector
# ----------------------------------------------------------------------------------- #
function estimate_steady_state(problem::Dict{String,Any})

    # Get the initial conditions -
    initial_condition_array = problem["initial_condition_array"];

    # calculate the steady-state soln -
    steady_state_prob = SteadyStateProblem(Balances, initial_condition_array, problem)
    steady_state_soln  = solve(steady_state_prob, DynamicSS(AutoTsit5(Rosenbrock23(autodiff=false));abstol=1e-8,reltol=1e-6,tspan=Inf))

    # return -
    return abs.(steady_state_soln.u)
end


text:
struct MMAModel{T, TV<:AbstractVector{T}, TC<:AbstractVector{<:Function}}
    dim::Int
    objective::Function
    ineq_constraints::TC
    box_max::TV
    box_min::TV
    # Trace flags
    store_trace::Bool
    show_trace::Bool
    extended_trace::Bool
    # Stopping criteria
    maxiter::Base.RefValue{Int}
    ftol::Base.RefValue{T}
    xtol::Base.RefValue{T}
    grtol::Base.RefValue{T}
end

dim(m::MMAModel) = m.dim
min(m::MMAModel, i::Integer) = m.box_min[i]
max(m::MMAModel, i::Integer) = m.box_max[i]
min(m::MMAModel)= m.box_max
max(m::MMAModel) = m.box_min
objective(m::MMAModel) = m.objective
constraints(m::MMAModel) = m.ineq_constraints
constraint(m::MMAModel, i::Integer) = m.ineq_constraints[i]
eval_objective(m, x::AbstractVector{T}, ∇g) where {T} = T(m.objective(x, ∇g))
eval_objective(m, x::AbstractVector{T}) where {T} = m.objective(x, T[])
eval_constraint(m, i, x::AbstractVector{T}, ∇g) where {T} = T(constraint(m, i)(x, ∇g))
eval_constraint(m, i, x::AbstractVector{T}) where {T} = constraint(m, i)(x, T[])
ftol(m) = m.ftol[]
xtol(m) = m.xtol[]
grtol(m) = m.grtol[]
ftol!(m, v) = m.ftol[] = v
xtol!(m, v) = m.xtol[] = v
grtol!(m, v) = m.grtol[] = v

MMAModel(dim, objective, args...; kwargs...) = MMAModel(whichdevice(objective), dim, objective, args...; kwargs...)
MMAModel(::CPU, args...; kwargs...) = MMAModel{Float64, Vector{Float64}, Vector{Function}}(args...; kwargs...)
MMAModel(::GPU, args...; kwargs...) = MMAModel{Float64, CuVector{Float64}, Vector{Function}}(args...; kwargs...)

function MMAModel{T, TV, TC}(dim,
                  objective::Function;
                  maxiter = 200,
                  xtol = eps(T),
                  ftol = sqrt(eps(T)),
                  grtol = sqrt(eps(T)),
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false) where {T, TV, TC}

    mins = ninfsof(TV, dim)
    maxs = infsof(TV, dim)
    MMAModel{T, TV, TC}(dim, objective, Function[],
             mins, maxs, store_trace, show_trace, extended_trace,
             Ref(maxiter), Ref(T(ftol)), Ref(T(xtol)), Ref(T(grtol)))
end

# Box constraints
function box!(m::MMAModel, i::Integer, minb::T, maxb::T) where {T}
    if !(1 <= i <= dim(m))
        throw(ArgumentError("box constraint need to applied to an existing variable"))
    end
    m.box_min[i] = minb
    m.box_max[i] = maxb
end

function box!(m::MMAModel, minb::T, maxb::T) where {T}
    nv = dim(m)
    m.box_min[1:nv] .= minb
    m.box_max[1:nv] .= maxb
end

function box!(m::MMAModel, minbs::AbstractVector{T}, maxbs::AbstractVector{T}) where {T}
    if (length(minbs) != dim(m)) || (length(minbs) != dim(m))
        throw(ArgumentError("box constraint vector must have same size as problem dimension"))
    end
    nv = dim(m)
    map!(identity, m.box_min, minbs)
    map!(identity, m.box_max, maxbs)
end

function ineq_constraint!(m::MMAModel, f::Function)
    push!(m.ineq_constraints, f)
end

function ineq_constraint!(m::MMAModel, fs::Vector{Function})
    for f in fs
        push!(m.ineq_constraints, f)
    end
end


text:
using seisimu

vp = [5000.0,5300,5800,6000]
vs = [3500.0,3200,4000,4500]
rho = 2.8
depth = [700,200,400,600]
Horizon = 1000
pkf = 10
T = 2
ext = 10
iflag = 2
dx = nothing
dz = nothing
dt = nothing
nT = nothing

model = initmodel(vp,vs,rho,depth,Horizon,pkf,T)


text:
# # Example: Gaussian with Given Variance
#
# As a model we are given a univariate Gaussian distribution with unknown mean
# and standard deviation σ = 1. We have one data point at 2.0. For this
# toy example, we can compute the posterior over the mean. But to illustrate
# likelihood-free inference, let us assume here that we can only sample from the
# model:
using LikelihoodfreeInference, Distributions, Random
model(x) = randn() .+ x

# LikelihoodfreeInference.jl passes parameter values as vectors to the model,
# even in the one-dimensional case. In our definition of the model we assume
# that `x[1]` is the mean.
#
# ## Approximate the Posterior
# Our first goal is to find the posterior over the mean given observation and
# a Gaussian prior with mean 0 and standard deviation 5.
data = [2.0]
prior = MultivariateNormal([0.], [5.])

# The true posterior is a Gaussian distribution with mean 25/26*2 and standard
# deviation 26/25
trueposterior = pdf.(Normal.(-1:.01:5, 26/25), 25/26*2.0)
using Plots
gr()
figure = plot(-1:.01:5, trueposterior, label = "posterior")

# Now, we will use an adaptive sequential Monte Carlo method:
smc = AdaptiveSMC(prior = prior, K = 10^4)
result = run!(smc, model, data, verbose = true, maxfevals = 10^6);

# As a Monte Carlo Method the result is a list of particles
particles(smc)
# with corresponding weights
weights(smc)

# The mean of the posterior is given by `weights(smc) .* particles(smc)`, which
# is computed by the `mean` function.
mean(smc)

figure = histogram(vcat(particles(smc)...), weights = weights(smc), normalize = true, label = "AdaptiveSMC")
plot!(figure, -1:.01:5, trueposterior, label = "posterior")

# The `result` above also contains these weights and particles and some
# additional information.
keys(result)

# AdaptiveSMC reduced the epsilon parameter adaptively, as we saw in column
# epsilon of the run above. We can plot this sequence.
scatter(cumsum(result.n_sims)[2:end], result.epsilons,
        yscale = :log10, ylabel = "epsilon", xlabel = "number of model evaluations")

# Alternatively, we may want to use KernelABC.
kabc = KernelABC(prior = prior,
                 kernel = Kernel(),
                 delta = 1e-12,
                 K = 10^4)
result = run!(kabc, model, data, maxfevals = 10^4)
mean(kabc)

figure = histogram(vcat(particles(kabc)...), weights = weights(kabc),
                   xlims = (-1, 5), bins = 100,
                   normalize = true, label = "KernelABC")
plot!(figure, -1:.01:5, trueposterior, label = "posterior")

# ## Point Estimates
# Sometimes we just want a point estimate. We will use BayesianOptimization.jl
# here to minimize the `QDLoss`. We know that the true maximum
# likelihood estimate is at mean = 25/26*2 ≈ 1.923
using BayesianOptimization
p = PointEstimator(optimizer = bo([-10.], [10.]), losstype = QDLoss, prior = prior, K = 100)
result = run!(p, model, data, maxfevals = 5*10^4, verbose = false);
result.x

# KernelRecursiveABC is an alternative method that requires often only few model
# evaluations in low and medium dimensional problems
k = KernelRecursiveABC(prior = prior,
                       kernel = Kernel(),
                       kernelx = Kernel(),
                       delta = 1e-2,
                       K = 100)
result = run!(k, model, data, maxfevals = 2*10^3)
result.x

# ## iid Samples
# Let us suppose here that the data consists of multiple independent and
# identically distributed samples.
data_iid = [[2.0], [1.9], [2.8], [2.1]]

# There are two ways to deal with this data. Either we just assume it is one
# four-dimensional vector
data_onevec = vcat(data_iid...)
# and we define the model as
model_iid_onevec(x) = vcat([model(x) for _ in 1:4]...)
smc = AdaptiveSMC(prior = prior, K = 10^4)
result = run!(smc, model_iid_onevec, data_onevec, verbose = true, maxfevals = 10^6);
histogram(vcat(particles(smc)...), weights = weights(smc), normalize = true, label = "AdaptiveSMC")

# Alternatively, we use another distance function:
model_iid(x) = [model(x) for _ in 1:4]
smc = AdaptiveSMC(prior = prior, K = 10^4, distance = energydistance)
result = run!(smc, model_iid, data_iid, verbose = true, maxfevals = 10^6);
histogram(vcat(particles(smc)...), weights = weights(smc), normalize = true, label = "AdaptiveSMC")


text:
include(joinpath(dirname(@__FILE__), "MyPetStore", "MyPetStore.jl"))

include("test_UserApi.jl")
include("test_StoreApi.jl")
include("test_PetApi.jl")

const server = "http://127.0.0.1/v2"
TestUserApi.test_404(server)
TestUserApi.test_userhook(server)
TestUserApi.test_set_methods()

if get(ENV, "STRESS_PETSTORE", "false") == "true"
    TestUserApi.test_parallel(server)
end

TestUserApi.test(server)
TestStoreApi.test(server)
TestPetApi.test(server)


text:
mutable struct CapacityMarket # Assumes a linear demand curve

    # Parameters

    maxprice::Float64 # $/MW/investment period
    midprice::Float64

    maxprice_capacity::Float64 # MW
    midprice_capacity::Float64
    zeroprice_capacity::Float64

    # Variables

    seg0contribution::VariableRef
    seg1contribution::VariableRef
    seg2contribution::VariableRef

    # Expressions

    totalcontribution::ExpressionRef # Capacity clearing level
    seg1price::ExpressionRef
    seg2price::ExpressionRef
    capacitywelfare::QuadExpressionRef # Capacity market welfare

    # Constraints

    min_seg0contribution::GreaterThanConstraintRef
    max_seg0contribution::LessThanConstraintRef

    min_seg1contribution::GreaterThanConstraintRef
    max_seg1contribution::LessThanConstraintRef

    min_seg2contribution::GreaterThanConstraintRef
    max_seg2contribution::LessThanConstraintRef

    max_totalcontribution::LessThanConstraintRef

    function CapacityMarket(
        maxprice::Float64, midprice::Float64,
        maxprice_capacity::Float64, midprice_capacity::Float64,
        zeroprice_capacity::Float64
    )
        @assert maxprice >= midprice >= 0
        @assert 0 <= maxprice_capacity <= midprice_capacity <= zeroprice_capacity
        new(maxprice, midprice,
            maxprice_capacity, midprice_capacity, zeroprice_capacity)
    end

end

function CapacityMarket(capacitypath::String)

    capacitydata = DataFrame!(CSV.File(joinpath(capacitypath, "rules.csv"),
                                       types=scenarios_capacitymarket_param_types))

    maxprice = first(capacitydata.maxprice)
    midprice = first(capacitydata.midprice)
    maxprice_capacity = first(capacitydata.maxprice_capacity)
    midprice_capacity = first(capacitydata.midprice_capacity)
    zeroprice_capacity = first(capacitydata.zeroprice_capacity)

    return CapacityMarket(
        maxprice, midprice,
        maxprice_capacity, midprice_capacity, zeroprice_capacity)

end

function setup!(market::CapacityMarket, m::Model, ops::Operations, s::AbstractScenario)

    # Variables

    market.seg0contribution = @variable(m)
    set_name(market.seg0contribution, "seg0contribution_capacity_$(s.name)")

    market.seg1contribution = @variable(m)
    set_name(market.seg1contribution, "seg1contribution_capacity_$(s.name)")

    market.seg2contribution = @variable(m)
    set_name(market.seg2contribution, "seg2contribution_capacity_$(s.name)")

    # Expressions

    market.totalcontribution = @expression(m,
        market.seg0contribution + market.seg1contribution + market.seg2contribution)

    market.seg1price = @expression(m,
         market.maxprice - (market.maxprice - market.midprice) /
                           (market.midprice_capacity - market.maxprice_capacity) *
                           market.seg1contribution)

    market.seg2price = @expression(m,
         market.midprice - market.midprice /
                           (market.zeroprice_capacity - market.midprice_capacity) *
                           market.seg2contribution)

    market.capacitywelfare = @expression(m,
        market.seg0contribution * market.maxprice +
        market.seg1contribution * (market.maxprice + market.seg1price) / 2 +
        market.seg2contribution * (market.midprice + market.seg2price) / 2)

    # Constraints

    market.min_seg0contribution = @constraint(m,
        market.seg0contribution >= 0)

    market.max_seg0contribution = @constraint(m,
        market.seg0contribution <= market.maxprice_capacity)

    market.min_seg1contribution = @constraint(m,
        market.seg1contribution >= 0)

    market.max_seg1contribution = @constraint(m,
        market.seg1contribution <= market.midprice_capacity - market.maxprice_capacity)

    market.min_seg2contribution = @constraint(m,
        market.seg2contribution >= 0)

    market.max_seg2contribution = @constraint(m,
        market.seg2contribution <= market.zeroprice_capacity - market.midprice_capacity)

    market.max_totalcontribution = @constraint(m,
        market.totalcontribution <= ucap(ops))

end

welfare(x::CapacityMarket) = x.capacitywelfare


text:
mutable struct RegimeParameterDecomposition
    bar::Float64
    hat::Array{Float64,1}
end

struct Block
    size::Int
    equations::BitArray{1}
    vars::BitArray{1}
    states::BitArray{1}
end

struct Meta
    numstates::Int
    numvars::Int
    numshocks::Int
    numregimes::Int
    numparameters::Int
    numstatic::Int
    numexo::Int
end

struct Parameters
    names::Array{String,1}
    values::Array{Float64,2}
    dictionary::Dict{String,Int}
    isswitching::BitArray{1}
    affectsSS::BitArray{1}
    decomposition::Array{RegimeParameterDecomposition,1}
    values_withbar::Array{Float64,2}
    generic_sym::Array{SymPy.Sym,2}
    generic_hats::Array{SymPy.Sym,2}
    bar_sym::Array{SymPy.Sym,1}
    specificregime_sym::Array{SymPy.Sym,2}
end

struct Variables
    names::Array{String,1}
    dictionary::Dict{String,Int}
    isstate::BitArray
    isstatic::BitArray
    isexo::BitArray
    states_indices::Array{Int,1}
    states_sym::Array{SymPy.Sym,1}
    leads_sym::Array{SymPy.Sym,1}
    contemps_sym::Array{SymPy.Sym,1}
end

struct Shocks
    names::Array{String,1}
    dictionary::Dict{String,Int}
    syms::Array{SymPy.Sym,1}
end

struct Equations
    stringversions::Array{String,1}
    isstatic::BitArray{1}
    isexo::BitArray{1}
    syms::Array{SymPy.Sym,1}
end

mutable struct SteadyState # mutable, because I need to wholesale reassign the steady state values all the time
    system::Array{Function,1}
    values::Array{Float64,1}
end

struct TransitionMatrix
    values::Array{Float64,2}
    genericsyms::Array{SymPy.Sym,1}
    ergodic::Array{Float64,1}
    completelygeneric::SymPy.Sym
end

struct PerturbationSystem
    A1::SparseMatrixCSC{Function}
    A2::SparseMatrixCSC{Function}
    A3::SparseMatrixCSC{Function}
    A4::SparseMatrixCSC{Function}
    B1::SparseMatrixCSC{Function}
    B2::SparseMatrixCSC{Function}
    B3::SparseMatrixCSC{Function}
    B4::SparseMatrixCSC{Function}
    B5::SparseMatrixCSC{Function}
    C1_states::SparseMatrixCSC{Function}
    C2_states::SparseMatrixCSC{Function}
    C_shocks::SparseMatrixCSC{Function}
    C_constant::SparseMatrixCSC{Function}
    C_constant_p::SparseMatrixCSC{Function}
    Xs::Array{SymPy.Sym,3}
    Xx::Array{SymPy.Sym,2}
    Xz::Array{SymPy.Sym,2}
    Js::Array{SymPy.Sym,3}
    Jx::Array{SymPy.Sym,2}
    Jz::Array{SymPy.Sym,2}
    isJ::BitArray{1}
    isY::BitArray{1}
    isW::BitArray{1}
    isZ::BitArray{1}
    isX::BitArray{1}
    blocks::BitArray{1}
end

struct RSDSGEModel
    meta::Meta
    parameters::Parameters
    vars::Variables
    shocks::Shocks
    equations::Equations
    steadystate::SteadyState
    transmatrix::TransitionMatrix
    perturbationsystem::PerturbationSystem
end


text:
module Classification

using LinearAlgebra

export SVC, Logistic

include("utils.jl")
include("Logistic.jl")
include("SVC.jl")

end


text:
# contours

using Plots
x = 1:0.5:20
y = 1:0.5:10
f(x, y) = begin
        (3x + y ^ 2) * abs(sin(x) + cos(y))
    end
X = repeat(reshape(x, 1, :), length(y), 1)
Y = repeat(y, 1, length(x))
Z = map(f, X, Y)
p1 = contour(x, y, f, fill = true)
p2 = contour(x, y, Z)
plot(p1, p2)

#Makie

using CairoMakie

f = Figure(backgroundcolor = :honeydew)
Axis(f[1, 1])

xs = LinRange(0, 10, 100)
ys = LinRange(0, 15, 100)
zs = [cos(x) * sin(y) for x in xs, y in ys]

CairoMakie.contour!(zs,levels=-1:0.1:1)

f

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations