text:
export ambient_space, rank, gram_matrix, inner_product, involution, ishermitian, is_quadratic, is_regular,
       is_local_square, is_isometric, is_rationally_isometric, is_isotropic, quadratic_space,
       hermitian_space, diagonal, invariants, hasse_invariant, witt_invariant, orthogonal_basis, fixed_field,
       restrict_scalars, orthogonal_complement

################################################################################
#
#  Maps
#
################################################################################

mutable struct AbsSpaceMor{D, T} <: Map{D, D, HeckeMap, AbsSpaceMor}
  header::MapHeader{D, D}
  matrix::T

  function AbsSpaceMor(V::D, W::D, B::T) where {D, T}
    z = new{D, T}()
    z.header = MapHeader{D, D}(V, W)
    z.matrix = B
    return z
  end
end

function hom(V::AbsSpace, W::AbsSpace, B::MatElem; check::Bool = false)
  @req base_ring(V) == base_ring(W) "Spaces must have the same base field"
  @req nrows(B) == dim(V) && ncols(B) == dim(W) """
  Dimension mismatch. Matrix ($(nrows(B))x$(ncols(B))) must be of
  dimensions $(dim(V))x$(dim(W)).
  """
  if check
    GV = gram_matrix(V)
    GW = gram_matrix(W)
    fl = B * GW * transpose(_map(B, involution(W))) == GV
    if !fl
      error("Matrix does not define a morphism of spaces")
    end
  end
  return AbsSpaceMor(V, W, B)
end

function image(f::AbsSpaceMor, v::Vector)
  V = domain(f)
  w = matrix(base_ring(V), 1, length(v), v) * f.matrix
  return vec(collect(w))
end

function compose(f::AbsSpaceMor, g::AbsSpaceMor)
  @req codomain(f) === domain(g) "incompatible morphisms"
  return hom(domain(f), codomain(g), f.matrix * g.matrix)
end

################################################################################
#
#  Creation
#
################################################################################

################################################################################
#
#  Basic invariants
#
################################################################################

@doc Markdown.doc"""
    rank(V::AbsSpace) -> Int

Return the rank of the space `V`.
"""
rank(L::AbsSpace) = rank(L.gram)

@doc Markdown.doc"""
    dim(V::AbsSpace) -> Int

Return the dimension of the space `V`.
"""
dim(V::AbsSpace) = nrows(V.gram)

@doc Markdown.doc"""
    gram_matrix(V::AbsSpace) -> MatElem

Return the Gram matrix of the space `V`.
"""
gram_matrix(V::AbsSpace) = V.gram

# Once we have quaternion spaces the following makes more sense

@doc Markdown.doc"""
    base_ring(V::AbsSpace) -> NumField

Return the algebra over which the space `V` is defined.
"""
base_ring(V::AbsSpace) = _base_algebra(V)

@doc Markdown.doc"""
    fixed_field(V::AbsSpace) -> NumField

Return the fixed field of the space `V`.
"""
fixed_field(::AbsSpace)

@doc Markdown.doc"""
    involution(V::AbsSpace) -> NumField

Return the involution of the space `V`.
"""
involution(V::AbsSpace)

################################################################################
#
#  Predicates
#
################################################################################

# TODO: Maybe cache this?
@doc Markdown.doc"""
    is_regular(V::AbsSpace) -> Bool

Return whether the space `V` is regular, that is, if the Gram matrix
has full rank.
"""
function is_regular(V::AbsSpace)
  return rank(V) == dim(V)
end

@doc Markdown.doc"""
    is_quadratic(V::AbsSpace) -> Bool

Return whether the space `V` is quadratic.
"""
is_quadratic(::AbsSpace)

@doc Markdown.doc"""
    ishermitian(V::AbsSpace) -> Bool

Return whether the space `V` is hermitian.
"""
ishermitian(::AbsSpace)

################################################################################
#
#  Determinant and discriminant
#
################################################################################

function det(V::AbsSpace)
  d = det(gram_matrix(V))
  return fixed_field(V)(d)
end

@doc Markdown.doc"""
    det(V::AbsSpace) -> FieldElem

Return the determinant of the space `V` as an element of its fixed field.
"""
det(::AbsSpace)

@doc Markdown.doc"""
    discriminant(V::AbsSpace) -> FieldElem

Return the discriminant of the space `V` as an element of its fixed field.
"""
function discriminant(V::AbsSpace)
  d = det(V)
  n = mod(rank(V), 4)
  if n == 0 || n == 1
    return d
  else
    return -d
  end
end

function _discriminant(G)
  d = det(G)
  n = nrows(G)
  if mod(n, 4) == 0 || mod(n, 4) == 1
    return d
  else
    return -d
  end
end

################################################################################
#
#  Computing inner products
#
################################################################################

@doc Markdown.doc"""
    gram_matrix(V::AbsSpace, M::MatElem) -> MatElem

Return the Gram matrix of the rows of `M` with respect to the Gram matrix of the space `V`.
"""
function gram_matrix(V::AbsSpace{T}, M::MatElem{S}) where {S, T}
  @req ncols(M) == dim(V) "Matrix must have $(dim(V)) columns ($(ncols(M)))"
  if S === elem_type(T)
    return M * gram_matrix(V) * transpose(_map(M, involution(V)))
  else
    Mc = change_base_ring(base_ring(V), M)
    return Mc * gram_matrix(V) * transpose(_map(Mc, involution(V)))
  end
end

@doc Markdown.doc"""
    gram_matrix(V::AbsSpace, S::Vector{Vector}) -> MatElem

Return the Gram matrix of the sequence `S` with respect to the Gram matrix of the space `V`.
"""
function gram_matrix(V::AbsSpace{T}, S::Vector{Vector{U}}) where {T, U}
  m = zero_matrix(base_ring(V), length(S), rank(V))
  for i in 1:length(S)
    if length(S[i]) != rank(V)
      throw(error("Vectors must be of length $(rank(V))"))
    end
    for j in 1:rank(V)
      m[i, j] = S[i][j]
    end
  end
  return gram_matrix(V, m)
end

@doc Markdown.doc"""
    inner_product(V::AbsSpace, v::Vector, w::Vector) -> FieldElem

Return the inner product of `v` and `w` with respect to the bilinear form of the space `V`.
"""
inner_product(V::AbsSpace, v::Vector, w::Vector)

@doc Markdown.doc"""
    inner_product(V::AbsSpace, v::MatElem, w::MatElem) -> MatElem

Shortcut for `v * gram_matrix(V) * adjoint(w)`.
"""
inner_product(V::AbsSpace, v::MatElem, w::MatElem)

_inner_product(L::AbsLat, v, w) = inner_product(ambient_space(L), v, w)

################################################################################
#
#  Diagonalization
#
################################################################################

@doc Markdown.doc"""
    orthogonal_basis(V::AbsSpace) -> MatElem

Return a matrix `M`, such that the rows of `M` form an orthgonal basis of the space `V`.
"""
function orthogonal_basis(V::AbsSpace)
  _, B = _gram_schmidt(gram_matrix(V), involution(V))
  return B
end

@doc Markdown.doc"""
    diagonal(V::AbsSpace) -> Vector{FieldElem}

Return a vector of elements $a_1,\dotsc,a_n$ such that the space `V` is isometric to
the diagonal space $\langle a_1,\dotsc,a_n \rangle$.

The elements are contained in the fixed field of `V`.
"""
diagonal(V::AbsSpace)

################################################################################
#
#  Gram-Schmidt
#
################################################################################

# Clean this up
function _gram_schmidt(M::MatElem, a, nondeg = true)
  F = deepcopy(M)
  K = base_ring(F)
  n = nrows(F)
  S = identity_matrix(K, n)
  T = identity_matrix(K, n)
  okk = is_diagonal(F)
  if !okk
    for i in 1:n
      if iszero(F[i,i])
        zero!(T)
        for i in 1:nrows(T)
          T[i, i] = 1
        end
        ok = 0
        for j in (i + 1):n
          if !iszero(F[j, j])
            ok = j
            break
          end
        end
        if ok != 0 # ok !== nothing
          j = ok
          T[i,i] = 0
          T[j,j] = 0
          T[i,j] = 1
          T[j,i] = 1
        else
          ok = 0
          for j in (i + 1):n
            if !iszero(F[i, j])
              ok = j
              break
            end
          end
          if ok == 0
            if nondeg
              error("Matrix is not of full rank")
            end
          else
            j = ok
            T[i, j] = 1 // (2 * F[j, i])
          end
        end
        if ok == 0
          continue
        end

        #S = T * S
        S = mul!(S, T, S)
        #F = T * F * transpose(_map(T, a))
        F = mul!(F, T, F)
        F = mul!(F, F, transpose(_map(T, a)))
      end

      zero!(T)
      for i in 1:nrows(T)
        T[i, i] = 1
      end

      for j in (i + 1):n
        T[j, i] = divexact(-F[j, i], F[i, i])
      end
      #F = T * F * transpose(_map(T, a))
      F = mul!(F, T, F)
      F = mul!(F, F, transpose(_map(T, a)))
      #S = T * S
      S = mul!(S, T, S)
    end
    @assert is_diagonal(F)
  end
  @hassert :Lattice 1 S * M * transpose(_map(S, a)) == F
  return F, S
end

################################################################################
#
#  Local isometry
#
################################################################################

@doc Markdown.doc"""
    is_isometric(L::AbsSpace, M::AbsSpace, p::Union{InfPlc, NfOrdIdl}) -> Bool

Return whether the spaces `L` and `M` are isometric over the completion at `p`.
"""
is_isometric(L::AbsSpace, M::AbsSpace, p)

################################################################################
#
#  Global isometry
#
################################################################################

@doc Markdown.doc"""
    is_isometric(L::AbsSpace, M::AbsSpace) -> Bool

Return whether the spaces `L` and `M` are isometric.
"""
is_isometric(L::AbsSpace, M::AbsSpace)

################################################################################
#
#  Definitness
#
################################################################################

# Returns 0 if V is not definite
# Returns an element a != 0 such that a * canonical_basis of V has
# positive Gram matrix
function _isdefinite(V::AbsSpace)
  E = base_ring(V)
  K = fixed_field(V)
  if (!is_totally_real(K)) || (ishermitian(V) && !is_totally_complex(E))
    return zero(K)
  end
  D = diagonal(V)
  signs_to_consider = Tuple{InfPlc, Int}[]
  for v in real_places(K)
    S = Int[sign(d, v) for d in D]
    if length(unique(S)) != 1
      return zero(K)
    else
      push!(signs_to_consider, (v, S[1]))
    end
  end
  if length(signs_to_consider) == 1
    return K(signs_to_consider[1][2])
  else
    return element_with_signs(K, signs_to_consider)
  end
end

@doc Markdown.doc"""
    is_positive_definite(V::AbsSpace) -> Bool

Return whether the space `V` is positive definite.
"""
function is_positive_definite(V::AbsSpace)
  E = base_ring(V)
  K = fixed_field(V)
  if (!is_totally_real(K)) || (ishermitian(V) && !is_totally_complex(E))
    return false
  end
  D = diagonal(V)
  for d in D
    if !is_totally_positive(d)
      return false
    end
  end
  return true
end

@doc Markdown.doc"""
    is_negative_definite(V::AbsSpace) -> Bool

Return whether the space `V` is negative definite.
"""
function is_negative_definite(V::AbsSpace)
  E = base_ring(V)
  K = fixed_field(V)
  if (!is_totally_real(K)) || (ishermitian(V) && !is_totally_complex(E))
    return false
  end
  D = diagonal(V)
  for d in D
    if !is_totally_positive(-d)
      return false
    end
  end
  return true
end

@doc Markdown.doc"""
    is_definite(V::AbsSpace) -> Bool

Return whether the space `V` is definite.
"""
function is_definite(V::AbsSpace)
  return is_positive_definite(V) || is_negative_definite(V)
end

################################################################################
#
#  Isotropic
#
################################################################################

@doc Markdown.doc"""
    is_isotropic(V::AbsSpace, p::Union{NfOrdIdl, InfPlc}) -> Bool

Given a space `V` and a place `p` in the fixed field `K` of `V`, return
whether the completion of `V` at `p` is isotropic.
"""
is_isotropic(::AbsSpace, p)

is_isotropic(V::AbsSpace, p::InfPlc) = _isisotropic(V, p)

# this is badly written, no need to compute d
function _isisotropic(D::Vector{fmpq}, p::PosInf)
  n = length(D)
  if n <= 1
    return false
  end
  E = parent(D[1])
  d = reduce(*, D, init = one(E))
  if d == 0
    return true
  elseif n <= 1
    return false
  else
    return length(unique!(fmpq[sign(d) for d in D])) == 2
  end
end

# this is badly written, no need to compute d
function _isisotropic(D::Vector, p::InfPlc)
  n = length(D)
  if n <= 1
    return false
  end
  E = parent(D[1])
  d = reduce(*, D, init = one(E))
  if d == 0
    return true
  elseif n <= 1
    return false
  elseif is_complex(p)
    return true
  else
    return length(unique!(Int[sign(d, p) for d in D])) == 2
  end
end

# this looks wrong
function _isisotropic(V::AbsSpace, p::InfPlc)
  n = rank(V)
  d = det(V)
  E = base_ring(V)
  if d == 0
    return true
  elseif n <= 1
    return false
  elseif is_complex(p)
    return true
  else
    D = diagonal(V)
    return length(unique!(Int[sign(d, p) for d in D])) == 2
  end
end

################################################################################
#
#  Restriction of scalars
#
################################################################################

# TODO: Use absolute_coordinates
function restrict_scalars(V::AbsSpace, K::FlintRationalField,
                                       alpha = one(base_ring(V)))
  E = base_ring(V)
  n = rank(V)
  d = absolute_degree(E)
  B = absolute_basis(E)
  v = elem_type(E)[zero(E) for i in 1:n]
  w = elem_type(E)[zero(E) for i in 1:n]
  G = zero_matrix(FlintQQ, d * n, d * n)
  r = 1
  c = 1
  for i in 1:n
    for bi in 1:length(B)
      v[i] = B[bi]
      c = 1
      for j in 1:n
        for bj in 1:length(B)
          w[j] = B[bj]
          G[r, c] = trace(alpha * inner_product(V, v, w), FlintQQ)
          w[j] = zero(E)
          c = c + 1
        end
      end
      v[i] = zero(E)
      r = r + 1
    end
  end

  return quadratic_space(FlintQQ, G, check = false), VecSpaceRes(E, rank(V))
end

################################################################################
#
#  Orthogonal complement
#
################################################################################

@doc Markdown.doc"""
    orthogonal_complement(V::AbsSpace, M::MatElem)

Given a space `V` and a subspace `W` with basis matrix `M`, returns a basis
matrix of the orthogonal complement of `W` inside `V`.
"""
function orthogonal_complement(V::AbsSpace, M::MatElem)
  N = gram_matrix(V) * _map(transpose(M), involution(V))
  r, K = left_kernel(N)
  @assert r == nrows(K)
  return K
end

################################################################################
#
#  Orthogonal sum
#
################################################################################

function _orthogonal_sum(V::AbsSpace, W::AbsSpace)
  K = base_ring(V)
  G = diagonal_matrix(gram_matrix(V), gram_matrix(W))
  n = dim(V) + dim(W)
  i1 = zero_matrix(K, dim(V), n)
  for i in 1:dim(V)
    i1[i, i] = 1
  end
  i2 = zero_matrix(K, dim(W), n)
  for i in 1:dim(W)
    i2[i, i + dim(V)] = 1
  end
  return G, i1, i2
end

@doc Markdown.doc"""
    orthogonal_sum(V::AbsSpace, W::AbsSpace) -> AbsSpace, AbsSpaceMor, AbsSpaceMor

Given two spaces `V` and `W` of the same kind (either both hermitian or both quadratic)
and defined over the same algebra, return their orthogonal sum $V \oplus W$. It is given with
the two natural embeddings $V \to V\oplus W$ and $W \to V\oplus W$.
"""
orthogonal_sum(V::AbsSpace, W::AbsSpace)

function orthogonal_sum(V::QuadSpace, W::QuadSpace)
  @req base_ring(V) === base_ring(W) "Base algebra must be equal"
  G, i1, i2 = _orthogonal_sum(V, W)
  VplusW = quadratic_space(base_ring(V), G)
  f1 = hom(V, VplusW, i1)
  f2 = hom(W, VplusW, i2)
  return VplusW, f1, f2
end

function orthogonal_sum(V::HermSpace, W::HermSpace)
  @req base_ring(V) === base_ring(W) "Base algebra must be equal"
  G, i1, i2 = _orthogonal_sum(V, W)
  VplusW = hermitian_space(base_ring(V), G)
  f1 = hom(V, VplusW, i1)
  f2 = hom(W, VplusW, i2)
  return VplusW, f1, f2
end

################################################################################
#
#  Embeddings
#
################################################################################

@doc Markdown.doc"""
    is_locally_represented_by(U::T, V::T, p::NfOrdIdl) where T <: AbsSpace -> Bool

Given two spaces `U` and `V` over the same algebra `E`, and a prime ideal `p` in
the maximal order $\mathcal O_K$ of their fixed field `K`, return whether `U` is
represented by `V` locally at `p`, i.e. whether $U_p$ embeds in $V_p$.
"""
is_locally_represented_by(::AbsSpace, ::AbsSpace, p)

@doc Markdown.doc"""
    is_represented_by(U::T, V::T) where T <: AbsSpace -> Bool

Given two spaces `U` and `V` over the same algebra `E`, return whether `U` is
represented by `V`, i.e. whether `U` embeds in `V`.
"""
is_represented_by(::AbsSpace, ::AbsSpace)


text:
replicated = replicate(data, w1)
coll = collapse(replicated'; order=4, standardize=true)
recombined = recombine(coll, w1)

@test all(repl .== replicated)
@test length(recombined) == length(w1)
@test all(recombined .== recombine(collapse(repl'; order=4, standardize=true), w1))

@test_throws DimensionMismatch replicate(data', w1)

ext: jl


text:
using RDatasets, DataFrames, FixedEffectModels, Base.Test

df = dataset("plm", "Cigar")
df[:pState] = pool(df[:State])
df[:pYear] = pool(df[:Year])

for method in [:cholesky, :qr, :lsmr]

	result = reg(Sales ~ Price |> pYear, df, save = true, method = method)
	@test_approx_eq  result.augmentdf[1, :pYear] 164.77833189721005

	result = reg(Sales ~ Price |> pYear + pState, df, save = true, method = method)
	@test_approx_eq_eps  result.augmentdf[1, :pYear] + result.augmentdf[1, :pState]  140.6852 1e-3

	result = reg(Sales ~ Price |> Year&pState, df, save = true, method = method)
	@test_approx_eq_eps result.augmentdf[1, :YearxpState]  1.742779  1e-3

	result = reg(Sales ~ Price |> pState + Year&pState, df, save = true, method = method)
	@test_approx_eq_eps  result.augmentdf[1, :pState]  -91.690635 1e-1

	result = reg(Sales ~ Price |> pState, df, save = true, subset = (df[:State] .<= 30), method = method)

	@test isna(result.augmentdf[1380,:pState])

end
	# add test with IV & weight


text:
export Circulant, Toeplitz

immutable Toeplitz{T} <: AbstractArray{T, 2}
	c :: Vector{T}
end

#XXX Inefficient but works
getindex(T::Toeplitz, i, j) = getindex(full(T), i, j)
isassigned(T::Toeplitz, i, j) = isassigned(full(T), i, j)

size(T::Toeplitz, r::Int) = (r==1 || r==2) ? 1 + div(length(T.c),2) :
    throw(ArgumentError("Invalid dimension $r"))
size(T::Toeplitz) = size(T,1), size(T,2)

function full{T}(To::Toeplitz{T})
	n=size(To, 1)
	M=Array(T, n, n)
	for i=1:n
		M[i:n,i] = To.c[n:2n-i]
		M[1:i-1,i] = To.c[n-i+1:n-1]
	end
	M
end

immutable Circulant{T} <: AbstractArray{T, 2}
	c :: Vector{T}
end

#XXX Inefficient but works
getindex(C::Circulant, i, j) = getindex(full(C), i, j)
isassigned(C::Circulant, i, j) = isassigned(full(C), i, j)

size(C::Circulant, r::Int) = (r==1 || r==2) ? length(C.c) :
    throw(ArgumentError("Invalid dimension $r"))
size(C::Circulant) = size(C,1), size(C,2)

function full{T}(C::Circulant{T})
	n=size(C, 1)
	M=Array(T, n, n)
	for i=1:n
		M[i:n,i] = C.c[1:n-i+1]
		M[1:i-1,i] = C.c[n-i+2:n]
	end
	M
end


text:
module LabeledVectors
export LabeledVector, indexof, labelof, unlabel, label, labels

struct LabeledVector{S,T} <: AbstractVector{T}
  elts::Vector{T}
  labels::Vector{S}
  indexes::Dict{S,Int}
end

flip(t::Tuple{S,T}) where {S,T} = (t[2],t[1])

function LabeledVector(labels::Vector{S}, elts::Vector{T}) where {S,T}
  @assert unique(labels) == labels
  @assert length(labels) == length(elts)
  label_to_idx = Dict{S,Int}(flip.(enumerate(labels)))
  LabeledVector{S,T}(elts, labels, label_to_idx)
end

function LabeledVector(vs::Vector{Pair{S,T}}) where {S,T}
  LabeledVector(first.(vs), last.(vs))
end

function LabeledVector{S,T}() where {S,T}
  LabeledVector(S[], T[])
end

labelof(lv::LabeledVector, i::Int) = lv.labels[i]

indexof(lv::LabeledVector{S}, l::S) where {S} = lv.indexes[l]

labels(lv::LabeledVector) = lv.labels

unlabel(lv::LabeledVector) = lv.elts

label(v::Vector{T}, labels::Vector{S}) where {S,T} = LabeledVector(labels, v)

Base.getindex(lv::LabeledVector, i::Int) = lv.elts[i]
Base.getindex(lv::LabeledVector{S}, l::S) where {S} = lv.elts[indexof(lv,l)]

Base.length(lv::LabeledVector) = length(lv.elts)
Base.size(lv::LabeledVector) = (length(lv),)

function Base.push!(lv::LabeledVector{S,T}, k::S, v::T) where {S,T}
  @assert k ∉ lv.labels
  push!(lv.elts, v)
  push!(lv.labels, k)
  lv.indexes[k] = length(lv.elts)
end

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule TOPCOM_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("TOPCOM")
JLLWrappers.@generate_main_file("TOPCOM", UUID("36f60fef-b880-50dc-9289-4aaecee93cc3"))
end  # module TOPCOM_jll


text:
function initialize_db(db_name::String)
    db = SQLite.DB(joinpath(db_name * ".sqlite"))
    create_party_table!(db)
    create_statement_table!(db)
    create_opinion_table!(db)
    create_results_table!(db)
    create_sequences_table!(db)
    create_parliament_table!(db)
    create_parameters_table!(db)
    return db
end

function create_party_table!(db::SQLite.DB)
    SQLite.execute(db, """
        CREATE TABLE IF NOT EXISTS party
        (
            party_id INTEGER NOT NULL PRIMARY KEY,
            party_shorthand TEXT,
            party_name TEXT
        );
    """)
    return true
end

function create_statement_table!(db::SQLite.DB)
    SQLite.execute(db, """
        CREATE TABLE IF NOT EXISTS statement
        (
            statement_id INTEGER NOT NULL PRIMARY KEY,
            statement_title TEXT,
            statement TEXT
        );
    """)
    return true
end

function create_opinion_table!(db::SQLite.DB)
    SQLite.execute(db, """
        CREATE TABLE IF NOT EXISTS opinion
        (
            party_id INTEGER NOT NULL,
            statement_id INTEGER NOT NULL,
            position INTEGER,
            position_rationale TEXT,
            FOREIGN KEY(party_id) REFERENCES party(party_id),
            FOREIGN KEY(statement_id) REFERENCES statement(statement_id),
            PRIMARY KEY(party_id, statement_id)
        );
    """)
    return true
end

function create_results_table!(db)
    SQLite.execute(db, """
        CREATE TABLE IF NOT EXISTS results
        (
            agent_id INTEGER NOT NULL,
            party TEXT,
            step INTEGER NOT NULL,
            rep INTEGER NOT NULL,
            statement_id INTEGER NOT NULL,
            position REAL,
            seq INTEGER NOT NULL,
            batchname TEXT,
            PRIMARY KEY(agent_id, step, rep, statement_id, seq, batchname)
        );
    """)
    return true
end

function create_sequences_table!(db)
    SQLite.execute(db, """
        CREATE TABLE IF NOT EXISTS sequences
        (
            seq_id INTEGER NOT NULL,
            step INTEGER,
            party_1 TEXT,
            party_2 TEXT,
            batchname TEXT,
            FOREIGN KEY(seq_id, batchname) REFERENCES results(seq, batchname),
            PRIMARY KEY(seq_id, batchname, step)
        );
    """)
    return true
end

function create_parliament_table!(db)
    SQLite.execute(db, """
        CREATE TABLE IF NOT EXISTS parliament
        (
            party_shorthand TEXT,
            batchname TEXT,
            seats INTEGER,
            PRIMARY KEY(party_shorthand, batchname)
        );
    """)
    return true
end

function create_parameters_table!(db)
    SQLite.execute(db, """
        CREATE TABLE IF NOT EXISTS parameters
        (
            batchname TEXT,
            group_size INTEGER,
            parliament_majority INTEGER,
            required_consensus REAL,
            PRIMARY KEY(batchname)
        );
    """)
    return true
end


text:
struct Product{A,B}
    a::A
    b::B
end
@inline function Base.size(p::Product)
    M = @inbounds size(p.a)[1]
    (M, Base.tail(size(p.b))...)
end
@inline function Base.size(p::Product, i::Integer)
    i == 1 && return @inbounds size(p.a)[1]
    @inbounds size(p.b)[i]
end
@inline Base.length(p::Product) = prod(size(p))
@inline Base.broadcastable(p::Product) = p
@inline numdims(A) = ndims(A) # fallback
@inline numdims(::Type{Product{A,B}}) where {A,B} = numdims(B)
@inline Base.ndims(::Type{Product{A,B}}) where {A,B} = numdims(B)
# This numdims nonsense is a hack to avoid type piracy in defining:
@inline numdims(::Type{B}) where {N, S <: Base.Broadcast.AbstractArrayStyle{N}, B <: Base.Broadcast.Broadcasted{S}} = N

Base.Broadcast._broadcast_getindex_eltype(::Product{A,B}) where {T, A <: AbstractVecOrMat{T}, B <: AbstractVecOrMat{T}} = T
function Base.Broadcast._broadcast_getindex_eltype(p::Product)
    promote_type(
        Base.Broadcast._broadcast_getindex_eltype(p.a),
        Base.Broadcast._broadcast_getindex_eltype(p.b)
    )
end

# recursive_eltype(::Type{A}) where {T, A <: AbstractArray{T}} = T
# recursive_eltype(::Type{NTuple{N,T}}) where {N,T<:Union{Float32,Float64}} = T
# recursive_eltype(::Type{Float32}) = Float32
# recursive_eltype(::Type{Float64}) = Float64
# recursive_eltype(::Type{Tuple{T}}) where {T} = T
# recursive_eltype(::Type{Tuple{T1,T2}}) where {T1,T2} = promote_type(recursive_eltype(T1), recursive_eltype(T2))
# recursive_eltype(::Type{Tuple{T1,T2,T3}}) where {T1,T2,T3} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3))
# recursive_eltype(::Type{Tuple{T1,T2,T3,T4}}) where {T1,T2,T3,T4} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3), recursive_eltype(T4))
# recursive_eltype(::Type{Tuple{T1,T2,T3,T4,T5}}) where {T1,T2,T3,T4,T5} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3), recursive_eltype(T4), recursive_eltype(T5))

# function recursive_eltype(::Type{Broadcasted{S,A,F,ARGS}}) where {S,A,F,ARGS}
#     recursive_eltype(ARGS)
# end

"""
    A *ˡ B

A lazy product of `A` and `B`. While functionally identical to `A * B`, this may avoid the
need for intermediate storage for any computations in `A` or `B`.  Example:

    @avx @. a + B *ˡ (c + d')

which is equivalent to

     a .+ B * (c .+ d')

It should only be used inside an `@avx` block, and to materialize the result it cannot be
the final operation.
"""
@inline *ˡ(a::A, b::B) where {A,B} = Product{A,B}(a, b)
@inline Base.Broadcast.broadcasted(::typeof(*ˡ), a::A, b::B) where {A, B} = Product{A,B}(a, b)
# TODO: Need to make this handle A or B being (1 or 2)-D broadcast objects.
function add_broadcast!(
    ls::LoopSet, mC::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(prod::Type{<:Product}), elementbytes::Int
)
    A, B = prod.parameters
    K = gensym(:K)
    mA = gensym(:Aₘₖ)
    mB = gensym(:Bₖₙ)
    pushpreamble!(ls, Expr(:(=), mA, Expr(:(.), bcname, QuoteNode(:a))))
    pushpreamble!(ls, Expr(:(=), mB, Expr(:(.), bcname, QuoteNode(:b))))
    pushpreamble!(ls, Expr(:(=), K, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__,Symbol(@__FILE__)), Expr(:ref, Expr(:call, :size, mB), 1))))
    k = gensym(:k)
    add_loop!(ls, Loop(k, 1, K), k)
    m = loopsyms[1];
    if numdims(B) == 1
        bloopsyms = Symbol[k]
        cloopsyms = Symbol[m]
        reductdeps = Symbol[m, k]
        kvec = bloopsyms
    elseif numdims(B) == 2
        n = loopsyms[2];
        bloopsyms = Symbol[k,n]
        cloopsyms = Symbol[m,n]
        reductdeps = Symbol[m, k, n]
        kvec = Symbol[k]
    else
        throw("B must be a vector or matrix.")
    end
    # load A
    # loadA = add_load!(ls, gensym(:A), productref(A, mA, m, k), elementbytes)
    loadA = add_broadcast!(ls, gensym(:A), mA, Symbol[m,k], A, elementbytes)
    # load B
    loadB = add_broadcast!(ls, gensym(:B), mB, bloopsyms, B, elementbytes)
    # set Cₘₙ = 0
    # setC = add_constant!(ls, zero(promote_type(recursive_eltype(A), recursive_eltype(B))), cloopsyms, mC, elementbytes)
    # targetC will be used for reduce_to_add
    mCt = gensym(mC)
    targetC = add_constant!(ls, gensym(:zero), cloopsyms, mCt, elementbytes, :numericconstant)
    push!(ls.preamble_zeros, (identifier(targetC), IntOrFloat))
    setC = add_constant!(ls, gensym(:zero), cloopsyms, mC, elementbytes, :numericconstant)
    push!(ls.preamble_zeros, (identifier(setC), IntOrFloat))
    setC.reduced_children = kvec
    # compute Cₘₙ += Aₘₖ * Bₖₙ
    reductop = Operation(
        ls, mC, elementbytes, :vfmadd231, compute, reductdeps, kvec, Operation[loadA, loadB, setC]
    )
    reductop = pushop!(ls, reductop, mC)
    reductfinal = Operation(
        ls, mCt, elementbytes, :reduce_to_add, compute, cloopsyms, kvec, Operation[reductop, targetC]
    )
    pushop!(ls, reductfinal, mCt)
end

struct LowDimArray{D,T,N,A<:DenseArray{T,N}} <: DenseArray{T,N}
    data::A
end
@inline Base.pointer(A::LowDimArray) = pointer(A.data)
Base.@propagate_inbounds Base.getindex(A::LowDimArray, i...) = getindex(A.data, i...)
Base.size(A::LowDimArray) = Base.size(A.data)
@generated function VectorizationBase.stridedpointer(A::LowDimArray{D,T,N}) where {D,T,N}
    s = Expr(:call, Expr(:(.), Expr(:(.), :LoopVectorization, QuoteNode(:VectorizationBase)), QuoteNode(:staticmul)), T, Expr(:tuple, [Expr(:ref, :strideA, n) for n ∈ 1+D[1]:N if D[n]]...))
    f = D[1] ? :PackedStridedPointer : :SparseStridedPointer
    Expr(:block, Expr(:meta,:inline), Expr(:(=), :strideA, Expr(:call, :strides, Expr(:(.), :A, QuoteNode(:data)))),
         Expr(:call, Expr(:(.), :VectorizationBase, QuoteNode(f)), Expr(:call, :pointer, Expr(:(.), :A, QuoteNode(:data))), s))
end
function LowDimArray{D}(data::A) where {D,T,N,A <: AbstractArray{T,N}}
    LowDimArray{D,T,N,A}(data)
end
function extract_all_1_array!(ls::LoopSet, bcname::Symbol, N::Int, elementbytes::Int)
    refextract = gensym(bcname)
    ref = Expr(:ref, bcname); append!(ref.args, [1 for n ∈ 1:N])
    pushpreamble!(ls, Expr(:(=), refextract, ref))
    return add_constant!(ls, refextract, elementbytes) # or replace elementbytes with sizeof(T) ? u
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(LDA::Type{<:LowDimArray}), elementbytes::Int
)
    D,T,N::Int,_ = LDA.parameters
    if !any(D)
        return extract_all_1_array!(ls, bcname, N, elementbytes)
    end
    fulldims = Symbol[loopsyms[n] for n ∈ 1:N if D[n]::Bool]
    ref = ArrayReference(bcname, fulldims)
    add_simple_load!(ls, destname, ref, elementbytes, true, false )::Operation
end
function add_broadcast_adjoint_array!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{A}, elementbytes::Int
) where {T,N,A<:AbstractArray{T,N}}
    parent = gensym(:parent)
    pushpreamble!(ls, Expr(:(=), parent, Expr(:call, :parent, bcname)))
    # isone(length(loopsyms)) && return extract_all_1_array!(ls, bcname, N, elementbytes)
    ref = ArrayReference(parent, Symbol[loopsyms[N + 1 - n] for n ∈ 1:N])
    add_simple_load!( ls, destname, ref, elementbytes, true, true )::Operation
end
function add_broadcast_adjoint_array!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{<:AbstractVector}, elementbytes::Int
)
    # isone(length(loopsyms)) && return extract_all_1_array!(ls, bcname, N, elementbytes)
    ref = ArrayReference(bcname, Symbol[loopsyms[2]])
    add_simple_load!( ls, destname, ref, elementbytes, true, true )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Adjoint{T,A}}, elementbytes::Int
) where {T, A <: AbstractArray{T}}
    add_broadcast_adjoint_array!( ls, destname, bcname, loopsyms, A, elementbytes )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Transpose{T,A}}, elementbytes::Int
) where {T, A <: AbstractArray{T}}
    add_broadcast_adjoint_array!( ls, destname, bcname, loopsyms, A, elementbytes )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{<:AbstractArray{T,N}}, elementbytes::Int
) where {T,N}
    add_simple_load!(ls, destname, ArrayReference(bcname, @view(loopsyms[1:N])), elementbytes, true, true)
end
function add_broadcast!(
    ls::LoopSet, ::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{T}, elementbytes::Int
) where {T<:Number}
    add_constant!(ls, bcname, elementbytes) # or replace elementbytes with sizeof(T) ? u
end
function add_broadcast!(
    ls::LoopSet, ::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{Base.RefValue{T}}, elementbytes::Int
) where {T}
    refextract = gensym(bcname)
    pushpreamble!(ls, Expr(:(=), refextract, Expr(:ref, bcname)))
    add_constant!(ls, refextract, elementbytes) # or replace elementbytes with sizeof(T) ? u
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{SubArray{T,N,A,S,B}}, elementbytes::Int
) where {T,N,N2,A<:AbstractArray{T,N2},B,N3,S <: Tuple{Int,Vararg{Any,N3}}}
    inds = Vector{Symbol}(undef, N+1)
    inds[1] = Symbol("##DISCONTIGUOUSSUBARRAY##")
    inds[2:end] .= @view(loopsyms[1:N])
    add_simple_load!(ls, destname, ArrayReference(bcname, inds), elementbytes, true, true)
end
BroadcastedArray{S<:Broadcast.AbstractArrayStyle,F,A} = Broadcasted{S,Nothing,F,A}
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(B::Type{<:BroadcastedArray}),
    elementbytes::Int
)
    S,_,F,A = B.parameters
    instr = get(FUNCTIONSYMBOLS, F) do
        f = gensym(:func)
        pushpreamble!(ls, Expr(:(=), f, Expr(:(.), bcname, QuoteNode(:f))))
        Instruction(bcname, f)
    end
    args = A.parameters
    Nargs = length(args)
    bcargs = Expr(:(.), bcname, QuoteNode(:args))
    # this is the var name in the loop
    parents = Operation[]
    deps = Symbol[]
    # reduceddeps = Symbol[]
    for (i,arg) ∈ enumerate(args)
        argname = gensym(:arg)
        pushpreamble!(ls, Expr(:(=), argname, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__,Symbol(@__FILE__)), Expr(:ref, bcargs, i))))
        # dynamic dispatch
        parent = add_broadcast!(ls, gensym(:temp), argname, loopsyms, arg, elementbytes)::Operation
        push!(parents, parent)
        mergesetdiffv!(deps, loopdependencies(parent), reduceddependencies(parent))
        # if !(isload(parent) || isconstant(parent))# && parent.instruction.instr ∉ (:reduced_add, :reduced_prod, :reduce_to_add, :reduce_to_prod)
            # mergesetv!(reduceddeps, reduceddependencies(parent))
        # end
        # pushparent!(parents, deps, reduceddeps, parent)
    end
    op = Operation(
        length(operations(ls)), destname, elementbytes, instr, compute, deps, NOPARENTS, parents
    )
    pushop!(ls, op, destname)
end

# size of dest determines loops
# function vmaterialize!(
@generated function vmaterialize!(
    dest::StridedArray{T,N}, bc::BC, ::Val{Mod}
) where {T <: NativeTypes, N, BC <: Broadcasted, Mod}
    # we have an N dimensional loop.
    # need to construct the LoopSet
    # @show typeof(dest)
    loopsyms = [gensym(:n) for n ∈ 1:N]
    ls = LoopSet(Mod)
    ls.isbroadcast[] = true
    sizes = Expr(:tuple)
    for (n,itersym) ∈ enumerate(loopsyms)
        Nsym = gensym(:N)
        add_loop!(ls, Loop(itersym, 1, Nsym), itersym)
        push!(sizes.args, Nsym)
    end
    pushpreamble!(ls, Expr(:(=), sizes, Expr(:call, :size, :dest)))
    elementbytes = sizeof(T)
    add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
    add_simple_store!(ls, :dest, ArrayReference(:dest, loopsyms), elementbytes)
    resize!(ls.loop_order, num_loops(ls)) # num_loops may be greater than N, eg Product
    # return ls
    q = lower(ls)
    push!(q.args, :dest)
    pushfirst!(q.args, Expr(:meta,:inline))
    # @show q
    q
     # ls
end
@generated function vmaterialize!(
    dest′::Union{Adjoint{T,A},Transpose{T,A}}, bc::BC, ::Val{Mod}
) where {T <: NativeTypes, N, A <: StridedArray{T,N}, BC <: Broadcasted, Mod}
    # we have an N dimensional loop.
    # need to construct the LoopSet
    loopsyms = [gensym(:n) for n ∈ 1:N]
    ls = LoopSet(Mod)
    ls.isbroadcast[] = true
    pushpreamble!(ls, Expr(:(=), :dest, Expr(:call, :parent, :dest′)))
    sizes = Expr(:tuple)
    for (n,itersym) ∈ enumerate(loopsyms)
        Nsym = gensym(:N)
        add_loop!(ls, Loop(itersym, 1, Nsym), itersym)
        push!(sizes.args, Nsym)
    end
    pushpreamble!(ls, Expr(:(=), sizes, Expr(:call, :size, :dest′)))
    elementbytes = sizeof(T)
    add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
    add_simple_store!(ls, :dest, ArrayReference(:dest, reverse(loopsyms)), elementbytes)
    resize!(ls.loop_order, num_loops(ls)) # num_loops may be greater than N, eg Product
    q = lower(ls)
    push!(q.args, :dest′)
    pushfirst!(q.args, Expr(:meta,:inline))
    q
    # ls
end
function vmaterialize!(
    dest::StridedArray{T,N}, bc::Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(identity),Tuple{T2}}, ::Val{Mod}
) where {T <: NativeTypes, N, T2 <: Number, Mod}
    arg = T(first(bc.args))
    @avx for i ∈ eachindex(dest)
        dest[i] = arg
    end
    dest
end
function vmaterialize!(
    dest′::Union{Adjoint{T,A},Transpose{T,A}}, bc::Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(identity),Tuple{T2}}, ::Val{Mod}
) where {T <: NativeTypes, N, A <: StridedArray{T,N}, T2 <: Number, Mod}
    arg = T(first(bc.args))
    dest = parent(dest′)
    @avx for i ∈ eachindex(dest)
        dest[i] = arg
    end
    dest′
end

@inline function vmaterialize(bc::Broadcasted, ::Val{Mod}) where {Mod}
    ElType = Base.Broadcast.combine_eltypes(bc.f, bc.args)
    vmaterialize!(similar(bc, ElType), bc, Val{Mod}())
end

vmaterialize!(dest, bc, ::Val{mod}) where {mod} = Base.Broadcast.materialize!(dest, bc)


text:
function sqldf_R(query::String)::DataFrame

    # Normalize Query
    query = replace(query,"\n"=>" ")

    # Extract tables in query
    tables = qtables(query)

    # Prepare R
    for t in tables
        Main.eval(Main.Meta.parse("@rput "*t))
    end

    # Retrieve and Return
    rquery = "c7e881694dfe63f2 = sqldf('"*query*"')"
    rquery = replace(rquery,"\"" => "\\\"")

    eval(Meta.parse("R\""*rquery*"\""))
    R"colnames(c7e881694dfe63f2) <- make.unique(colnames(c7e881694dfe63f2),sep=':')"
    res = rcopy(R"c7e881694dfe63f2")

    # Clean Up
    R"rm(c7e881694dfe63f2)"
    for t in tables
        eval(Meta.parse("R"*"\"rm("*t*")\""))
    end
    R"gc()"

    return res
end


text:
# Benchmark sparseKernel.jl
using Meganet, BenchmarkTools

# Setup example from tests
history = Pkg.dir("Meganet")*"/benchmarks/micro/bm_sparseKernel.jld"
nK = [100; 80]
nT = 5000;
A = sprandn(nK[1],nK[2],.5)
Q = sprandn(nnz(A),nT,.7)
K = getSparseKernel(Float64,A,Q)

# JthetaTmv
th  = initTheta(K)
dth = initTheta(K)
nex = 1
Y  = randn(nFeatIn(K),nex)+nex
Z  = randn(nFeatOut(K),nex)-nex

out = JthetaTmv(K, Z, th, Y)
out_old = Meganet.JthetaTmv_old(K, Z, th, Y)
@code_warntype JthetaTmv(K, Z, th, Y)
trial = @benchmark JthetaTmv(K, Z, th, Y)

# Save benchmarks
if true
    Meganet.updatehistory!(history, trial)
    hist = JLD.load(history, "hist")
    judge(hist)
end


text:
using ZipFile

mutable struct NLModel <: NEOSModel
    solver::NEOSSolver
    xmlmodel::String
    last_results::String
    inner::AmplNLWriter.AmplNLMathProgModel
    status::Symbol
end

MPB.LinearQuadraticModel(s::NEOSSolver{S, :NL}) where {S} = NLModel(
    s,
    "",
    "",
    AmplNLWriter.AmplNLMathProgModel("", String[], ""),
    NOTSOLVED
)

function MPB.loadproblem!(m::NLModel, A::AbstractMatrix, x_l, x_u, c, g_l, g_u, sense)
    AmplNLWriter.loadproblem!(AmplNLWriter.AmplNLLinearQuadraticModel(m.inner), A, x_l, x_u, c, g_l, g_u, sense)
end

function neos_writexmlmodel!(m::NLModel)
    # There is no non-linear binary type, only non-linear discrete, so make
    # sure binary vars have bounds in [0, 1]
    for i in 1:m.inner.nvar
        if m.inner.vartypes[i] == :Bin
            if m.inner.x_l[i] < 0
                m.inner.x_l[i] = 0
            end
            if m.inner.x_u[i] > 1
                m.inner.x_u[i] = 1
            end
        end
    end
    io = IOBuffer()
    print(io, "<model>")
    AmplNLWriter.make_var_index!(m.inner)
    AmplNLWriter.make_con_index!(m.inner)
    AmplNLWriter.write_nl_file(io, m.inner)
    print(io, "</model>")
    # Convert the model to NL and add
    m.xmlmodel = replace(m.solver.template, r"<model>.*</model>"is => String(take!(io)))
end

# Wrapper functions
for f in [:getvartype,:getsense,:status,:getsolution,:getobjval,:numvar,:numconstr,:get_solve_result,:get_solve_result_num,:get_solve_message,:get_solve_exitcode,:getsolvetime]
  @eval $f(m::NLModel) = $f(m.inner)
end
for f in [:setvartype!,:setsense!,:setwarmstart!]
  @eval $f(m::NLModel, x) = $f(m.inner, x)
end

function solfilename(job)
    "https://neos-server.org/neos/jobs/$(10_000 * floor(Int, job.number / 10_000))/$(job.number)-$(job.password)-solver-output.zip"
end

MPB.getreducedcosts(m::NLModel) = fill(NaN, numvar(m))
MPB.getconstrduals(m::NLModel) = fill(NaN, numconstr(m))
MPB.getobjbound(m::NLModel) = NaN

function MPB.status(m::NLModel)
    m.status
end

function parseresults!(m::NLModel, job)
    m.status = SOVLERERROR

    # https://neos-server.org/neos/jobs/5710000/5711322-FLWbgxPt-solver-output.zip
    if m.solver.print_level >= 2
        @info("Getting solution file from $(solfilename(job))...")
    end
    try
        res = HTTP.request("GET",solfilename(job))
    catch
        @info("No results found...")
        m.status = NORESULTS
        return
    end

    if res.status != 200
        @error("Error retrieving results for job $(job.number):$(job.password). Response status is $(res.status).")
        m.status = NORESULTS
        return
    end

    if m.solver.print_level >= 2
        @info("Extracting file from .zip")
    end
    io = IOBuffer()
    write(io, res.data)
    z = ZipFile.Reader(io)
    @assert length(z.files) == 1 # there should only be one .sol file in here
    sol = readstring(z.files[1])
    close(io)
    io = IOBuffer()
    write(io, sol)
    if m.solver.print_level >= 2
        @info("Reading results")
    end
    seekstart(io)
    AmplNLWriter.read_results(io, m.inner)
    m.status = SOLVED
end


text:
module TestACSetViews
using Test

using Catlab: @present, generator
using Catlab.Theories: compose, id
using Catlab.Graphs
using Catlab.Graphs.BasicGraphs: TheoryGraph
using Catlab.CSetDataStructures
using Catlab.CategoricalAlgebra.ACSetViews
using Catlab.CategoricalAlgebra.CSets

using Tables

# Basic Indexing and Special Methods
####################################

@present TheoryLocatedGraph <: TheoryGraph begin
  Num::Data
  x::Attr(V,Num)
  y::Attr(V,Num)
end

const LocatedGraph = ACSetType(TheoryLocatedGraph)

g = @acset LocatedGraph{Float64} begin
  V = 5
  E = 4
  src = [1,1,1,1]
  tgt = [2,3,4,5]
  x = [0.,1.,0,-1.,0]
  y = [0.,0.,1.,0.,-1.]
end

edges = ACSetView(g,:E)

upward_edges = @select_where(edges, src.y <= tgt.y)
@test parts(upward_edges) == [1,2,3]

@test upward_edges[1,:src] == 1.
@test upward_edges[1,[:src,:x]] == 0.
@test upward_edges.src == [1,1,1]
@test upward_edges[:,[:tgt,:x]] == [1.,0.,-1.]

@test @compute_prop(upward_edges, (src.x - tgt.x)^2 + (src.y - tgt.y)^2) == [1,1,1]

# Tables.jl Interface
#####################

vertices = ACSetView(g,:V)

@test Tables.istable(vertices)
@test Tables.columnaccess(vertices)
@test Tables.columns(vertices) == vertices
@test Tables.columnnames(vertices) == (:x,:y)
@test Tables.getcolumn(vertices,:x)[2] == 1.
@test Tables.getcolumn(vertices,2,:x) == 1.

# Morphisms
###########

g, h = WeightedGraph{Float64}(2), WeightedGraph{Float64}(4)
add_edge!(g, 1, 2, weight=2.0)
add_edges!(h, [1,2,3], [2,3,4], weight=[1.0,2.0,3.0])
α = ACSetTransformation((V=[2,3], E=[2]), g, h)
αE = ACSetViewMorphism(α, :E)

@test is_natural(αE)
@test αE.f == α.components[:E]

end


text:
struct BruteForce <: HamiltonianSolver end

min_path(T, stops, ::BruteForce) = _min_length_path_exhaustive(T, stops)

function _min_length_path_exhaustive(T::TravelTimes, stops::Vector{<:AbstractStop})
    _min_length_path_exhaustive(T.travel_Times, stops)
end

# sdlve the minimum length hamiltonian path by brute-force search
# computes the minimum path of the given stops, starting from the first stop given
# in `stops` and finishing in any other other stops
function _min_length_path_exhaustive(T::AbstractMatrix, stops::Vector{<:AbstractStop})
    p0 = stops[1]

    # FIXME refactor
    stop_names = [name(s) for s in stops]
    M = view(T, stop_names, stop_names)

    perm = permutations(stops[2:end])

    num_stops = length(stops)

    minvalue = Inf
    local result
    for p in perm
        value = M[name(p0), name(p[1])]
        for i in 1:num_stops-2
            value += M[name(p[i]), name(p[i+1])]
        end

        if value < minvalue
            minvalue = value
            result = vcat(p0, p)
        end
    end

    return (minvalue, result)
end


text:
function normAm(A,m)
    #NORMAM   Estimate of 1-norm of power of matrix.
    #   NORMAM(A,m) estimates norm(A^m,1).
    #   If A has nonnegative elements the estimate is exact.
    #   [C,MV] = NORMAM(A,m) returns the estimate C and the number MV of
    #   matrix-vector products computed involving A or A^*.

    #   Reference: A. H. Al-Mohy and N. J. Higham, A New Scaling and Squaring
    #   Algorithm for the Matrix Exponential, SIAM J. Matrix Anal. Appl. 31(3):
    #   970-989, 2009.

    #   Awad H. Al-Mohy and Nicholas J. Higham, September 7, 2010.

    t = 1; # Number of columns used by NORMEST1.

    n = size(A,1);

    # function afun_power(flag,X)
    #     #AFUN_POWER  Function to evaluate matrix products needed by NORMEST1.
    #
    #     if flag == "dim"
    #         Z = n;
    #     elseif flag == "real"
    #         if isreal(A)
    #             Z = 1
    #         else
    #             Z = 0
    #         end
    #     else
    #         (p,q) = size(X);
    #         if p != n,
    #             error("Dimension mismatch")
    #         end
    #
    #         if flag=="notransp"
    #             for i = 1:m
    #                 X = A*X;
    #             end
    #         elseif flag=="transp"
    #             for i = 1:m
    #                 X = A'*X
    #             end
    #         end
    #
    #         Z = X;
    #
    #     end
    #
    #     return Z
    # end

    if isequal(A,abs(A))
        e = ones(n,1);
        for j=1:m         # for positive matrices only
            e = A'*e;
        end
        c = norm(e,Inf);
        mv = m;
    else
        #(c,_,_,it) = normest1(@afun_power,t);
        #mv = it[2]*t*m;
        c = norm(A^m,1);
        mv = 6*t*m;
    end

    return (c,mv)

end


text:
#!/usr/bin/env julia

print(sum(factorial.(BigInt.(1:3000))))

max_stars_repo_name: robindaumann/fac-bench

max_issues_repo_name: robindaumann/fac-bench

max_forks_repo_name: robindaumann/fac-bench


text:
#####
##### univariate building blocks
#####

"""
$(SIGNATURES)

Represents the sequence of Chebyshev polynomials (of the first kind).
"""
struct Chebyshev end

Base.Broadcast.broadcastable(C::Chebyshev) = Ref(C)

####
#### set length calculations
####

"""
$(SIGNATURES)

Length of a disjoint set with index `k ≥ 1`.
"""
function set_length(::Chebyshev, k::Integer)
    @argcheck k > 0
    if k == 1
        1
    elseif k == 2
        2
    else
        2^(k - 2)
    end
end

"""
$(SIGNATURES)

Indexes for disjoint set with index `k ≥ 1`.
"""
function set_range(::Chebyshev, k::Integer)
    @argcheck k > 0
    if k == 1
        1:1
    elseif k == 2
        2:3
    else
        K = 2^(k - 2)
        start = 2 + K
        start:(start + K - 1)
    end
end

"""
$(SIGNATURES)

Set of additional (Chebyshev-Lobatto) gridpoints for index `k ≥ 1`. Sets are disjoint,
contain elements in increasing order, which for `k ∈ 2:K` contain the Chebyshev-Lobatto
gridpoints

```math
cos(πj / (N+ 1))
```
for ``j=0, …, N+1`` where ``N = 2ᴷ``, in an interleaved pattern. `k == 1` is just `[0.0]`.
"""
function gridpoint_set(::Chebyshev, k::Integer)
    if k == 1
        [0.0]
    elseif k == 2
        [-1.0, 1.0]
    else
        @argcheck k ≥ 1
        N = 2^(k - 1)
        cospi.(((N-1):-2:1) ./ N)
    end
end

function all_gridpoints(::Chebyshev, K::Integer)
    mapreduce(k -> gridpoint_set(Chebyshev(), k), vcat, 1:K)
end

# struct ChebyshevPolynomials{T}
#     x::T
# end

# Base.IteratorEltype(::Type{ChebyshevPolynomials{T}}) where {T} = float(T)

# Base.IteratorSize(::Type{<:ChebyshevPolynomials}) = Base.IsInfinite()

# function Base.iterate(c::ChebyshevPolynomials{T}) where {T}
#     S = float(T)
#     one(S), (one(S) / 2, zero(S))
# end

# function Base.iterate(c::ChebyshevPolynomials, (yp, ypp))
#     y = 2 * c.x * yp - ypp
#     y, (y, yp)
# end

function evaluate!(z::AbstractVector, ::Chebyshev, x::Real)
    N = length(z)
    if N ≥ 1
        ypp = z[1] = 1
        if N ≥ 2
            yp = z[2] = x
            for i in 3:N
                y = 2 * yp * x - ypp
                ypp = yp
                yp = y
                z[i] = y
            end
        end
    end
    z
end

"""
$(SIGNATURES)

Evaluate the first `N` basis functions at `x`, returning a vector.
"""
function evaluate(N::Integer, C::Chebyshev, x::Real)
    evaluate!(Vector{float(eltype(x))}(undef, N), C, x)
end

"""
$(SIGNATURES)

In-place version of [`evaluate`](@ref), returns the first argument.
"""
function evaluate!(Z::AbstractMatrix, ::Chebyshev, x::AbstractVector)
    N = size(Z, 2)
    @argcheck size(Z, 1) == length(x)
    if N ≥ 1
        Z[:, 1] .= 1
        if N ≥ 2
            Z[:, 2] .= x
            for j in 3:N
                @. Z[:, j] = 2 * x * Z[:, j - 1] - Z[:, j - 2]
            end
        end
    end
    Z
end

"""
$(SIGNATURES)

Evaluate the first `N` basis functions at elements of `x`, one row for each element.
"""
function evaluate(N::Integer, C::Chebyshev, x::AbstractVector)
    evaluate!(Matrix{float(eltype(x))}(undef, length(x), N), C, x)
end


text:
using Cxx
using OpenNI2
using CVCore
using CVHighGUI

const ni2 = OpenNI2

isesc(key) = key == 27
genfilename(ext=".png") =
    joinpath(dirname(@__FILE__), string(now(), "_", time_ns(), ext))

# Turn this on if you prefer to use Kinect v2
# assumed that libfreenect2-openni2 driver is installed
use_libfreenect2_openni2_driver = false

always_save = false
show_ir = true
show_depth = true

if use_libfreenect2_openni2_driver
    w, h = 512, 424
else
    w, h = 640>>1, 480>>1
end

function setVideoMode(stream, si, w, h, pxfmt)
    modes = ni2.getSupportedVideoModes(si)
    found = false
    for mode in modes
        if Int(ni2.getResolutionX(mode)) == w &&
                Int(ni2.getResolutionY(mode)) == h &&
                ni2.getPixelFormat(mode) == pxfmt
            println("$w x $h : video mode found")
            ni2.setVideoMode(stream, mode)
            found = true
            break
        end
    end
    !found && error("video mode not found")
end

ni2.initialize()
device = ni2.DevicePtr()
ni2.open(device)
ni2.setDepthColorSyncEnabled(device, false)

di = ni2.getDeviceInfo(device)
@show ni2.getName(di)
@show ni2.getVendor(di)

depth = ni2.VideoStreamPtr()
ni2.create(depth, device, ni2.SENSOR_DEPTH)
!use_libfreenect2_openni2_driver && ni2.setMirroringEnabled(depth, true)
setVideoMode(depth, ni2.getSensorInfo(device, ni2.SENSOR_DEPTH), w, h,
    ni2.PIXEL_FORMAT_DEPTH_1_MM)

ir = ni2.VideoStreamPtr()
ni2.create(ir, device, ni2.SENSOR_IR)
!use_libfreenect2_openni2_driver && ni2.setMirroringEnabled(ir, true)
setVideoMode(ir, ni2.getSensorInfo(device, ni2.SENSOR_IR), w, h,
    ni2.PIXEL_FORMAT_GRAY16)

ni2.setImageRegistrationMode(device, ni2.IMAGE_REGISTRATION_OFF)

foreach(ni2.start, [depth, ir])

frame = ni2.VideoFrameRef()

while true
    readyIndex = ni2.waitForAnyStream([depth, ir])
    if readyIndex == 0
        ni2.readFrame(depth, frame)
        arr = convert(Array{ni2.DepthPixel,2}, frame)
        scaledarr = arr * 1/maximum(arr)
        show_depth && imshow("depth", scaledarr)
    elseif readyIndex == 1
        ni2.readFrame(ir, frame)
        arr = convert(Array{ni2.Grayscale16Pixel,2}, frame)
        scaledarr = arr * 1/maximum(arr)
        show_ir && imshow("ir", scaledarr)
    end

    key = waitKey(delay=1)
    isesc(key) && break

    if key == 's' || always_save
        scale!(scaledarr, 255.)
        ext = readyIndex == 0 ? "_depth.png" : "_ir.png"
        fname = genfilename(ext)
        @show fname
        imwrite(fname, Mat(scaledarr))
    end

    rand() > 0.95 && gc(false)
end

destroyAllWindows()

foreach(ni2.stop, [depth, ir])
foreach(ni2.destroy, [depth, ir])
ni2.close(device)
ni2.shutdown()

# Force deallocate
depth=0;ir=0;device=0;gc()


text:
using NIfTI
using FinEtools

# file = niread("example4d.nii", mmap=false)
# display(file.header)
# voxel_size(file.header)
# size(file)
# vb = VoxelBoxVolume(file.raw[:, :, :, 1], voxel_size(file.header))
# vtkexport("example4d.vtk", vb)

# file = niread("bigger_subset_MW_Ear_R_rotY-thinner-5mat-fin.nii")
# display(file.header)
# voxel_size(file.header)
# size(file)
# vb = VoxelBoxVolume(file.raw[:, :, :, 1], voxel_size(file.header))
# vtkexport("bigger_subset_MW_Ear_R_rotY-thinner-5mat-fin.vtk", vb)

File = "Head_256x256x126-256x256x252"
file = niread(File * ".nii")
display(file.header)
voxel_size(file.header)
size(file)
vb = VoxelBoxVolume(file.raw[:, :, :, 1], voxel_size(file.header))
vtkexport(File * ".vtk", vb)


text:
__precompile__()

module MultiscaleGraphSignalTransforms

include("GraphSignal.jl")
include("GraphPartition.jl")
include("BasisSpecification.jl")
include("common.jl")
include("GHWT.jl")
include("GHWT_2d.jl")
include("GHWT_tf_1d.jl")
include("GHWT_tf_2d.jl")
include("HGLET.jl")
include("GraphSig_Plot.jl")
include("gplot.jl")
include("partition_fiedler.jl")

using Reexport
@reexport using .GraphSignal, .GraphPartition, .BasisSpecification, .GHWT, .GHWT_2d, .GHWT_tf_1d, .GHWT_tf_2d, .HGLET

export dvec2dmatrix, dmatrix2dvec, levlist2levlengths!, bsfull, bs_haar, bs_level, bs_walsh, dvec_Threshold, rs_to_region, GraphSig_Plot, gplot, gplot!, partition_fiedler
export cost_functional, dmatrix_flatten, dmatrix_ldb_flatten

## export functions of NGWP.jl
using LinearAlgebra, SparseArrays, LightGraphs, SimpleWeightedGraphs, Clustering
using JuMP, Clp, Optim, Statistics, QuadGK, Arpack
import Plots: plot, plot!, scatter, scatter!
import StatsBase: crosscor

include("dualgraph.jl")
include("eigDAG_Distance.jl")
include("eigHAD_Distance.jl")
include("eigROT_Distance.jl")
include("eigsROT_Distance.jl")
include("eigTSD_Distance.jl")
include("helpers.jl")
include("LP-HGLET.jl")
include("LP-NGWP.jl")
include("natural_distances.jl")
include("utils.jl")
include("varimax.jl")
include("NGWF.jl")
include("NGWP.jl")
include("PC-NGWP.jl")
include("SunFlowerGraph.jl")
include("VM-NGWP.jl")

export eigDAG_Distance, eigHAD_Distance, eigHAD_Affinity
export eigROT_Distance, ROT_Distance, eigEMD_Distance
export eigsROT_Distance
export eigTSD_Distance, K_functional
export natural_eigdist
export SunFlowerGraph, dualgraph
export pc_ngwp, pairclustering, mgslp
export vm_ngwp, varimax
export lp_ngwp, rising_cutoff, find_pairinds, pair_inds_shadding, lp_ngwp_analysis
export LPHGLET_Synthesis, LPHGLET_Analysis_All, HGLET_dictionary, LPHGLET_dictionary
export HGLET_DST4_Synthesis, LPHGLET_DST4_Synthesis, HGLET_DST4_dictionary, LPHGLET_DST4_dictionary
export HGLET_DST4_Analysis, LPHGLET_DST4_Analysis
export unitary_folding_operator, keep_folding!
export ngwp_analysis, ngwp_bestbasis, NGWP_jkl
export natural_eigdist
export nat_spec_filter, ngwf_all_vectors, rngwf_all_vectors, ngwf_vector, frame_approx, rngwf_lx
export scatter_gplot, scatter_gplot!, stem, stem!
export standardize_eigenvectors!, spike, characteristic, χ, sort_wavelets, transform2D
export getall_expansioncoeffs, approx_error_plot, getall_expansioncoeffs2, approx_error_plot2

end


text:
#################################
# Indices API                   #
#################################
# Get daily index values for a specific location. Index availability varies by location.

# Contains a centralized function to make calls to the API.
include("request.jl")

function get_indexvalues_for_group_of_indices(apikey::String, time::Int64, groupID::String, locationKey::String, language::String, details::Bool)
    # Returns daily index data for a specific group of indices, by location key.

    url = "http://dataservice.accuweather.com/indices/v1/daily/"

    # Add required parameters to the URL
    url = url * string(time) * "day/"
    # Allowed values for 'time': '1', '5', '10' and '15'.

    url = url * locationKey * "/groups/"
    url = url * groupID
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end
    if details == true url = url * "&details=true" end

    return make_API_call(url)
end

function get_indexvalues_for_index(apikey::String, time::Int64, indexID::String, locationKey::String, language::String, details::Bool)
    # Returns daily index data for a specific index, by location key.

    url = "http://dataservice.accuweather.com/indices/v1/daily/"

    # Add required parameters to the URL
    url = url * string(time) * "day/"
    # Allowed values for 'time': '1', '5', '10' and '15'.

    url = url * locationKey * "/"
    url = url * indexID
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end
    if details == true url = url * "&details=true" end

    return make_API_call(url)
end

function get_indexvalues_all(apikey::String, time::Int64, indexID::String, locationKey::String, language::String, details::Bool)
    # Returns daily index data for all indices, by location key.

    url = "http://dataservice.accuweather.com/indices/v1/daily/"

    # Add required parameters to the URL
    url = url * string(time) * "day/"
    # Allowed values for 'time': '1', '5', '10' and '15'.

    url = url * locationKey
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end
    if details == true url = url * "&details=true" end

    return make_API_call(url)
end

function get_indices(apikey::String, language::String)
    # Returns metadata for all daily indices.

    url = "http://dataservice.accuweather.com/indices/v1/daily"

    # Add required parameters to the UR
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end

    return make_API_call(url)
end

function get_indexgroups(apikey::String, language::String)
    # Returns metadata for all index groups.

    url = "http://dataservice.accuweather.com/indices/v1/daily/groups"

    # Add required parameters to the UR
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end

    return make_API_call(url)
end

function get_indices_for_group(apikey::String, groupID::String, language::String)
    # Returns metadata for all indices in a specified group.

    url = "http://dataservice.accuweather.com/indices/v1/daily/groups/"

    # Add required parameters to the URL
    url = url * groupID
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end

    return make_API_call(url)
end

function get_index(apikey::String, indexID::String, language::String)
    # Returns metadata for a specific index type.

    url = "http://dataservice.accuweather.com/indices/v1/daily/"

    # Add required parameters to the URL
    url = url * indexID
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end

    return make_API_call(url)
end


text:
# ellipse to rotated and shifted ellipse
# initial and final landmark positions observed

# just as exp1, but only larger noise on landmarks 9--18,
# and also updating of the initial positions
using Revise

using BridgeLandmarks
const BL = BridgeLandmarks
using RCall
using Random
using Distributions
using DataFrames
using DelimitedFiles
using CSV
using StaticArrays
using LinearAlgebra
using JLD2
using FileIO

Random.seed!(9)

workdir = @__DIR__
cd(workdir)
include(joinpath(BL.dir(),"scripts", "postprocessing.jl"))
outdir = workdir
mkpath(joinpath(outdir, "forward"))

#-------- read data ----------------------------------------------------------
dat = load("data_exp1.jld2")
xobs0 = dat["xobs0"]
xobsT = dat["xobsT"]
n = dat["n"]
nshapes = dat["nshapes"]

################################# start settings #################################
ITER = 500
subsamples = 0:10:ITER

model = [:ms, :ahs][1]
fixinitmomentato0 = false
obs_atzero = true
updatescheme =  [:innov, :mala_mom]

# noise on observations
σobs0 = 0.01
Σobs0 = [σobs0^2 * one(UncF) for i in 1:n]
σobsT = vcat(fill(σobs,9), fill(0.1,n-9))
ΣobsT = [σobsT[i]^2 * one(UncF) for i in 1:n]
Σobs = [Σobs0, ΣobsT]

T = 1.0; dt = 0.01; t = 0.0:dt:T; tt_ =  tc(t,T)

################################# MCMC tuning pars #################################
ρinit = 0.9              # pcN-step
covθprop = [0.04 0. 0.; 0. 0.04 0.; 0. 0. 0.04]
if model==:ms
    δinit = [0.0001, 0.1] # first comp is not used
else
    δinit = [0.001, 0.1] # first comp is not used
end
η(n) = min(0.2, 10/n)  # adaptation rate for adjusting tuning pars
adaptskip = 20  # adapt mcmc tuning pars every adaptskip iters
maxnrpaths = 10 # update at most maxnrpaths Wiener increments at once
tp = tuningpars_mcmc(ρinit, maxnrpaths, δinit,covθprop,η,adaptskip)

################################# initialise P #################################
ainit = mean(norm.([xobs0[i]-xobs0[i-1] for i in 2:n]))/2.0
if model == :ms
    cinit = 0.2
    γinit = 2.0
    P = MarslandShardlow(ainit, cinit, γinit, 0.0, n)
elseif model == :ahs
    cinit = 0.02
    γinit = 0.2
    stdev = 0.75
    nfsinit = construct_nfs(2.5, stdev, γinit)
    P = Landmarks(ainit, cinit, n, 2.5, stdev, nfsinit)
end

################## prior specification with θ = (a, c, γ) ########################
priorθ = product_distribution([Exponential(ainit), Exponential(cinit), Exponential(γinit)])
κ = 100.0
priormom = MvNormalCanon(zeros(d*n), gramkernel(xobs0,P)/κ)

#########################
xobsT = [xobsT]
xinit = State(xobs0, zeros(PointF,P.n))
mT = zeros(PointF, n)

start = time() # to compute elapsed time
    Xsave, parsave, objvals, accpcn, accinfo, δ, ρ, covθprop =
    lm_mcmc(tt_, (xobs0,xobsT), Σobs, mT, P,
              obs_atzero, fixinitmomentato0, ITER, subsamples,
              xinit, tp, priorθ, priormom, updatescheme,
            outdir)
elapsed = time() - start

#----------- post processing -------------------------------------------------
println("Elapsed time: ",round(elapsed/60;digits=2), " minutes")
perc_acc_pcn = mean(accpcn)*100
println("Acceptance percentage pCN step: ", round(perc_acc_pcn;digits=2))
write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir)
write_info(model,ITER, n, tt_, updatescheme, Σobs, tp, ρ, δ, perc_acc_pcn, elapsed, outdir)
write_observations(xobs0, xobsT, n, nshapes, x0,outdir)
write_acc(accinfo, accpcn, nshapes,outdir)
write_params(parsave, 0:ITER, outdir)
write_noisefields(P, outdir)


text:
scatter_defaults() = Dict([
    :strokewidth => 0,
    :color => main_color()
]);

"""
	$(SIGNATURES)

Line graph. Simple wrapper around `lines`.
"""
function scatter_plot(xV, yV :: AbstractVector{F};
    fig = blank_plot(), pos = (1,1), kwargs ...) where F

    args = merge(scatter_defaults(), kwargs);
    ax = fig[pos...] = Axis(fig; args...);
    scatter!(ax, xV, yV; args...);
    return fig, ax
end

function scatter_plot(yV :: AbstractVector{F}; kwargs ...) where F
    return scatter_plot(1 : length(yV), yV; kwargs...);
end

function scatter_plot(xV, yM :: AbstractMatrix{F};
    fig = blank_plot(), pos = (1,1), kwargs...) where F
    nr, nc = size(yM);
    @assert nr == length(xV);
    args = merge(scatter_defaults(), kwargs);
    ax = make_axis(fig, pos; args...);
    for j = 1 : nc
        add_scatter!(ax, xV, yM[:, j]; args..., color = get_colors(j, nc));
    end
    return fig, ax
end

"""
	$(SIGNATURES)

Add line to a plot.
"""
function add_scatter!(ax :: Axis, x, y; kwargs...)
    scatter!(ax, x, y; kwargs...);
end

add_scatter!(p :: Makie.FigureAxisPlot, x, y; kwargs...) =
    add_scatter!(p.axis, x, y; kwargs...);

# -----------


text:
#---------------------------------------
# types for a single EB sample
#---------------------------------------

"""
    EBayesSample{T}

Abstract type representing empirical Bayes samples with realizations of type `T`.
"""
abstract type EBayesSample{T} end

abstract type ContinuousEBayesSample{T} <: EBayesSample{T} end
abstract type DiscreteEBayesSample{T} <: EBayesSample{T} end

"""
    likelihood_distribution(Z::EBayesSample, μ::Number)

Returns the distribution ``p(\\cdot \\mid \\mu)`` of ``Z \\mid \\mu`` (the return type being
a `Distributions.jl` Distribution).

# Examples
```jldoctest
julia> likelihood_distribution(StandardNormalSample(1.0), 2.0)
Normal{Float64}(μ=2.0, σ=1.0)
```
"""
function likelihood_distribution end

"""
    response(Z::EBayesSample{T})

Returns the concrete realization of `Z` as type `T`, thus dropping the information about the
likelihood.

# Examples
```jldoctest
julia> response(StandardNormalSample(1.0))
1.0
```
"""
function response(Z::EBayesSample)
   Z.Z
end

function set_response(Z::EBayesSample, znew=missing)
    Z = @set Z.Z = znew
end

function nuisance_parameter end

function Base.Float64(Z::EBayesSample{<:Number})
    Base.Float64(response(Z))
end

Base.isnan(Z::EBayesSample) = Base.isnan(response(Z))
Base.isfinite(Z::EBayesSample) = Base.isfinite(response(Z))

function Base.isless(a::EBayesSample, b::EBayesSample)
    Base.isless( (response(a), nuisance_parameter(a)),
                 (response(b), nuisance_parameter(b)) )
end

Base.isless(a::EBayesSample, b::Number) = Base.isless(response(a),b)
Base.isless(a::Number, b::EBayesSample) = Base.isless(a, response(b))

"""
	marginalize(Z::EBayesSample, prior::Distribution)

Given a `prior` distribution ``G`` and  `EBayesSample` ``Z``,
return that marginal distribution of ``Z``. Works for `EBayesSample{Missing}``,
i.e., no realization is needed.

# Examples
```jldoctest
julia> marginalize(StandardNormalSample(1.0), Normal(2.0, sqrt(3)))
Normal{Float64}(μ=2.0, σ=1.9999999999999998)
````
"""
function marginalize end

function marginalize(Z::EBayesSample, G::Dirac)
    likelihood_distribution(Z, G.value)
end

function posterior end

Broadcast.broadcastable(Z::EBayesSample) = Ref(Z)

function multiplicity(Zs)
    collect(StatsBase.weights(Zs))
end
function StatsBase.weights(Zs::AbstractVector{<:EBayesSample})
    uweights(length(Zs))
end

summarize_by_default(Zs) = false

# trait
abstract type Skedasticity end
struct Homoskedastic <: Skedasticity end
struct Heteroskedastic <: Skedasticity end

function skedasticity(Zs::AbstractVector{EB}) where {EB<:EBayesSample}
    # probably should check that we have a concrete type..
    #Base.typename(typeof(BinomialSample(2,3)))
    length(unique(nuisance_parameter.(Zs))) == 1 ? Homoskedastic() : Heteroskedastic()
end

# avoid piracy
likelihood(Z::EBayesSample, param) = _pdf(likelihood_distribution(Z, param), response(Z))
loglikelihood(Z::EBayesSample, param) =
    _logpdf(likelihood_distribution(Z, param), response(Z))
loglikelihood(Z::EBayesSample, prior::Distribution) =
    _logpdf(marginalize(Z, prior), response(Z))

function loglikelihood(Zs::AbstractVector{<:EBayesSample}, prior)
    sum(loglikelihood.(Zs, prior))
end

"""
    pdf(prior::Distribution, Z::EBayesSample)

Given a `prior` ``G`` and `EBayesSample` ``Z``, compute the marginal density of `Z`.

# Examples
```jldoctest
julia> Z = StandardNormalSample(1.0)
Z=     1.0 | σ=1.0
julia> prior = Normal(2.0, sqrt(3))
Normal{Float64}(μ=2.0, σ=1.7320508075688772)
julia> pdf(prior, Z)
0.17603266338214976
julia> pdf(Normal(2.0, 2.0), 1.0)
0.17603266338214976
```
"""
pdf(prior::Distribution, Z::EBayesSample) = _pdf(marginalize(Z, prior), response(Z))

"""
    cdf(prior::Distribution, Z::EBayesSample)

Given a `prior` ``G`` and `EBayesSample` ``Z``, evaluate the CDF of the marginal
distribution of ``Z`` at `response(Z)`.
"""
cdf(prior::Distribution, Z::EBayesSample) = _cdf(marginalize(Z, prior), response(Z))
"""
    ccdf(prior::Distribution, Z::EBayesSample)

Given a `prior` ``G`` and `EBayesSample` ``Z``, evaluate the complementary CDF of the marginal
distribution of ``Z`` at `response(Z)`.
"""
ccdf(prior::Distribution, Z::EBayesSample) = _ccdf(marginalize(Z, prior), response(Z))

logpdf(prior::Distribution, Z::EBayesSample) = _logpdf(marginalize(Z, prior), response(Z))

_pdf(dbn, z) = pdf(dbn, z)
_logpdf(dbn, z) = logpdf(dbn, z)

_cdf(dbn, z) = cdf(dbn, z)
_ccdf(dbn, z) = ccdf(dbn, z)
#const EBInterval{T} = Union{[Interval{T, S, R} for S in [Open, Closed, Unbounded], R in [Open, Closed, Unbounded]]...} where T

const EBInterval{T} = Union{
    Interval{T,Unbounded,Closed},
    Interval{T,Open,Unbounded},
    Interval{T,Open,Closed},
} where {T}

const EBIntervalFlipped{T} = Union{
    Interval{T,Closed,Open},
    Interval{T,Unbounded,Open},
    Interval{T,Closed,Unbounded},
} where {T}

function _cdf(dbn, interval::Interval{T,S,Unbounded}) where {T,S}
    one(eltype(dbn))
end

function _cdf(
    dbn::ContinuousUnivariateDistribution,
    interval::Interval{T,S,B},
) where {T,S,B<:Bounded}
    _cdf(dbn, last(interval))
end

#-- Unbounded - Unbounded
function _pdf(dbn, interval::Interval{T,Unbounded,Unbounded}) where {T}
    one(eltype(dbn))
end
function _logpdf(dbn, interval::Interval{T,Unbounded,Unbounded}) where {T}
    zero(eltype(dbn))
end

#-- Unbounded - Closed
function _pdf(dbn, interval::Interval{T,Unbounded,Closed}) where {T}
    cdf(dbn, last(interval))
end
function _logpdf(dbn, interval::Interval{T,Unbounded,Closed}) where {T}
    logcdf(dbn, last(interval))
end

#-- Open - Unbounded
function _pdf(dbn, interval::Interval{T,Open,Unbounded}) where {T}
    ccdf(dbn, first(interval))
end
function _logpdf(dbn, interval::Interval{T,Open,Unbounded}) where {T}
    logccdf(dbn, first(interval))
end

#-- Open - Closed
function _pdf(dbn, interval::Interval{T,Open,Closed}) where {T}
    cdf(dbn, last(interval)) - cdf(dbn, first(interval))
end
function _logpdf(dbn, interval::Interval{T,Open,Closed}) where {T}
    logdiffcdf(dbn, last(interval), first(interval))
end

# In the case of continuous distributions, we can delegate most results to  (Open, Closed)
for f in [:_pdf, :_logpdf]
    @eval begin
        # Closed - Closed and Closed - Open and Open - Open
        function $f(
            dbn::ContinuousUnivariateDistribution,
            interval::Union{
                Interval{T,Closed,Closed},
                Interval{T,Closed,Open},
                Interval{T,Open,Open},
            },
        ) where {T}
            _interval = Interval{T,Open,Closed}(first(interval), last(interval))
            $f(dbn, _interval)
        end
        # Unbounded - Open
        function $f(
            dbn::ContinuousUnivariateDistribution,
            interval::Interval{T,Unbounded,Open},
        ) where {T}
            _interval = Interval{T,Unbounded,Closed}(first(interval), last(interval))
            $f(dbn, _interval)
        end
        # Closed - Unbounded
        function $f(
            dbn::ContinuousUnivariateDistribution,
            interval::Interval{T,Closed,Unbounded},
        ) where {T}
            _interval = Interval{T,Open,Unbounded}(first(interval), last(interval))
            $f(dbn, _interval)
        end
    end
end

#TODO: Handle discrete distributions (less important)
function _pdf(dbn::DiscreteDistribution, interval::Interval{T,Closed,Closed}) where {T}
    cdf(dbn, last(interval)) - cdf(dbn, first(interval)) + pdf(dbn, first(interval))
end

function _cdf(
    dbn::DiscreteDistribution,
    interval::Interval{T,S,Closed},
) where {T,S}
    _cdf(dbn, last(interval))
end

function _pdf(dbn::DiscreteDistribution, interval::Interval{T,Closed,Unbounded}) where {T}
    ccdf(dbn, first(interval)) + pdf(dbn, first(interval))
end

function _logpdf(dbn::DiscreteDistribution, interval::Interval{T,Closed,Unbounded}) where {T}
    log(_pdf(dbn, interval))
end

function _support(d::Distribution)
    distributions_interval_to_interval(support(d))
end

function distributions_interval_to_interval(interval::Distributions.RealInterval)
    _lb = isinf(interval.lb) ? nothing : interval.lb
    _ub = isinf(interval.ub) ? nothing : interval.ub
    Interval(_lb, _ub)
end

struct MultinomialSummary{T,D<:AbstractDict{T,Int}}
    store::D #TODO, use other container
end

MultinomialSummary(vals, cnts) = MultinomialSummary(SortedDict(Dict(vals .=> cnts)))

const VectorOrSummary{T} = Union{AbstractVector{T},MultinomialSummary{T}}

# Does the distinction of the two things below really make sense?
function (Zs_summary::MultinomialSummary)(Z)
    get(Zs_summary.store, Z, zero(Int))
end

function Base.getindex(Zs_summary::MultinomialSummary, i)
    Base.getindex(Zs_summary.store, i)
end

Base.keys(Zs_summary::MultinomialSummary) = Base.keys(Zs_summary.store)
Base.values(Zs_summary::MultinomialSummary) = Base.values(Zs_summary.store)
Base.length(Zs_summary::MultinomialSummary) = Base.length(Zs_summary.store)

# TODO: Move to Dictionaries.jl

function Base.broadcasted(::typeof(response), Zs_summary::MultinomialSummary)
    response.(keys(Zs_summary))
end

function Base.broadcasted(::typeof(pdf), prior, Zs_summary::MultinomialSummary)
    # Should this also return keys?
    pdf.(prior, collect(keys(Zs_summary.store)))
end

function Base.broadcasted(f, Zs_summary::MultinomialSummary)
    # Should this also return keys?
    f.(keys(Zs_summary.store))
end

function StatsBase.weights(Zs_summary::MultinomialSummary)
    fweights(collect(values(Zs_summary.store)))
end

summarize(Zs::AbstractVector) = MultinomialSummary(SortedDict(countmap(Zs)))
summarize(Zs::AbstractVector, ws::StatsBase.AbstractWeights) = MultinomialSummary(SortedDict(countmap(Zs, ws)))

summarize(Zs::MultinomialSummary) = Zs

function skedasticity(Zs_summary::MultinomialSummary)
    all_unique_samples = collect(keys(Zs_summary.store))
    skedasticity(all_unique_samples)
end

function loglikelihood(mult::MultinomialSummary, prior)
    sum([n * loglikelihood(Z, prior) for (Z, n) in mult.store])
end

nobs(Zs_summary::MultinomialSummary) = sum(values(Zs_summary.store))
nobs(Zs::AbstractVector{<:EBayesSample}) = length(Zs)

StatsBase.fit(::Nothing, ::VectorOrSummary) = nothing

# Recall difference between truncating marginal VS truncated likelihood
# Truncated(Normal())
#------------------------------------------
# MarginalTruncated{...} and Truncated{...}
#------------------------------------------
# w Dirac prior these are the same!


text:
# ====================================================================
#                         MODULE DataFitting
# ====================================================================
module DataFitting

using Printf, PrettyTables
using Statistics, Distributions
using DataStructures

import Base.show
import Base.ndims
import Base.size
import Base.length
import Base.getindex
import Base.propertynames
import Base.getproperty
import Base.setproperty!
import Base.append!
import Base.isequal

export Domain, CartesianDomain, Measures,
    FuncWrap, ScalarParam, Smooth,
    Model, getparamvalues, setparamvalues!, resetcounters!,
    componenttype, addcomp!,
    add_dom!, rm_dom!, dom_count,
    addexpr!, replaceexpr!, setflag!,
    evaluate!, fit!, fit, probe,
    test_component

const compsep = "_"

include("Types.jl")
include("components.jl")
include("private.jl")
include("userinterface.jl")
include("show.jl")

end


text:
using Pumas, PumasTutorials, CSV

data = PumasTutorials.tutorial_data("data/nca","multiple_dose_IVbolus_7BLQ_test")
data = CSV.read(data,missingstring="NA")
first(data,10)

timeu = u"hr"
concu = u"mg/L"
amtu  = u"mg"

pop = read_nca(data, id=:ID, time=:time, conc=:DV, amt=:DOSE, ii=24timeu,
    route=:Formulation, occasion=:OCC,timeu=timeu, concu=concu, amtu=amtu,llq=0.4concu)

NCA.auc(pop,auctype=:last,method=:linear)

report = NCAReport(pop)
report = NCA.to_dataframe(report)

names(report)

CSV.write("./tutorials/nca/report_MD_IVbolus_7BLQ.csv", report)


text:
function run_singlethread(num_sims::Integer) :: Float64
    @assert num_sims > 0
    acum = 0

    for _ in 1:num_sims

        # `x` and `y` are scalars
        x = rand(Float64)
        y = rand(Float64)

        r = sqrt(x^2 + y^2)

        if r <= 1.0
            acum += 1
        end
    end

    area = acum / num_sims
    return 4*area
end


text:
using Plots

function main()
    x = 0:0.1:6
    y1 = sin.(x)
    y2 = cos.(x)

    plot(
        x, y1,
        title = "sin & cos",
        label = "sin",
        xlabel = "x",
        ylabel = "y"
    )
    plot!(
        x, y2,
        line = :dash,
        label = "cos"
    )
end

main()


text:
#!/usr/bin/env julia

#######
# CLI #
#######
import Celeste.ArgumentParse
parser = Celeste.ArgumentParse.ArgumentParser()
ArgumentParse.add_argument(
    parser,
    "--ais-output",
    help="Directory containing source-specific AIS sample files",
)
ArgumentParse.add_argument(
    parser,
    "--output-dir",
    help="Directory to output results dataframes",
)
ArgumentParse.add_argument(
    parser,
    "--truth-csv",
    help="True file (coadd, or prior generated file)",
)
ArgumentParse.add_argument(
    parser,
    "--vb-csv",
    help="VB inference file for comparison",
)
ArgumentParse.add_argument(
    parser,
    "--photo-csv",
    help="Photo inference file for comparison"
)
parsed_args = ArgumentParse.parse_args(parser, ARGS)
println(parsed_args)
println("======= scoring MCMC results ============")

#######################
# Analysis Output Dir #
#######################
import JLD, CSV
using DataFrames, StatsBase
import Celeste: AccuracyBenchmark, Model, MCMC

if haskey(parsed_args, "output-dir")
  output_dir = parsed_args["output-dir"]
else
  output_dir = "./mcmc_result_dataframes"
end
println(" ... saving output to ", output_dir)
if !isdir(output_dir)
  mkdir(output_dir)
end

#######################################################################
# collect MCMC files and summarize them into a results dataframe      #
#######################################################################
if haskey(parsed_args, "ais-output")
    ais_output_dir = parsed_args["ais-output"]
else
    ais_output_dir = "ais-output/synthetic-run-10-31"
end
ais_out = filter!(r"\.jld$", readdir(ais_output_dir))
results = [JLD.load(joinpath(ais_output_dir, f))["res"] for f in ais_out];
results = [res for res in results if !res["is_sky_bad"]]
@printf " ... found %d results files in directory %s \n" length(results) ais_output_dir

#####################################################
# Compute prob star/gal (based on model evidence    #
#####################################################
function compute_pstar(res; prior_is_star=.28) #[.28, .72])
    # is_star = [.28, .72] vs [.999, .001] vs [.5, .5]
    type_chain  = zeros(length(res["gal_bootstrap"]))
    lnprob_a    = log(prior_is_star)
    lnprob_nota = log(1.-prior_is_star)
    for n in 1:length(type_chain)
        # normalizing constant is ln p(data | star)
        # so p(star | data) \propto p(data | star) p(star)
        lnprob_star = res["star_bootstrap"][n] + lnprob_a
        lnprob_gal  = res["gal_bootstrap"][n] + lnprob_nota
        lnsum = Model.logsumexp([lnprob_star, lnprob_gal]) # normalize
        type_chain[n] = lnprob_star - lnsum
    end
    ave_pstar = Model.logsumexp(type_chain) - log(length(type_chain))
    return ave_pstar
end

for res in results
  res["ave_pstar"] = compute_pstar(res; prior_is_star=.28)
end

# must correct for off-by-pixel value
star_summary, gal_summary, results_df = MCMC.consolidate_samples(results;
    summarize=median)
#results_df[:dec] -= .396 / 3600.
#results_df[:ra]  -= .397 / 3600.

######################################################################
# Compute the per chain R Hats (TODO put these in the alg itself)    #
######################################################################
function compute_ss_rhat(sdf; num_chains=25, num_samps_per_chain=25)
    chains = []
    si = 1
    for i in 1:num_chains
        endi = si + num_samps_per_chain - 1
        #@printf "chunked %d : %d \n" si endi
        push!(chains, Array(sdf[si:endi,10:10]))
        si += num_samps_per_chain
    end
    return MCMC.potential_scale_reduction_factor(chains)
end

star_psrf = mean(vcat([compute_ss_rhat(res["star_samples"])
                       for res in results]...), 2)
gal_psrf  = mean(vcat([compute_ss_rhat(res["gal_samples"])
                       for res in results]...), 2)

######################################################################
# Load comparison data frames                                        #
######################################################################
if haskey(parsed_args, "truth-csv")
  truth_file = parsed_args["truth-csv"]
else
  truth_file = expanduser("~/Proj/Celeste.jl/benchmark/accuracy/output/prior_6328dda483.csv")
end
println(" ... loaded truth dataframe from ", truth_file)
truth = AccuracyBenchmark.read_catalog(truth_file)

if haskey(parsed_args, "vb-csv")
  vb_file = parsed_args["vb-csv"]
else
  vb_file = expanduser("~/Proj/Celeste.jl/benchmark/accuracy/output/prior_6328dda483_synthetic_6c24c1c9ad_predictions_342302e66a.csv")
end
println(" ... loaded VB comparison DF from ", vb_file)
vbdf = AccuracyBenchmark.read_catalog(vb_file)

if haskey(parsed_args, "photo-csv")
  photo_file = parsed_args["photo-csv"]
else
  println(" ... no photo file passed in --- loading VB as dummy ")
  photo_file = vb_file
end
photodf = AccuracyBenchmark.read_catalog(photo_file)

############################################################################
# Match prediction dataframes first --- throw away the same sources as VB  #
############################################################################
matched_truth, matched_prediction_dfs = AccuracyBenchmark.match_catalogs(
    truth, [vbdf, results_df, photodf]) #, tol=.396/3600*4)
@printf " ... found %d matched sources for comparison \n" size(matched_truth, 1)
truth = matched_truth
vbdf = matched_prediction_dfs[1]
results_df = matched_prediction_dfs[2]
photodf = matched_prediction_dfs[3]
@printf "Found %d nan r fluxes in the VB dataframe\n" sum(isnan.(vbdf[:flux_r_nmgy]))

#bad_idx = isnan.(vbdf[:flux_r_nmgy])
#truth = truth[.~bad_idx,:]
#vbdf = vbdf[.~bad_idx,:]
#results_df = results_df[.~bad_idx,:]
#photodf = photodf[.~bad_idx,:]

############################################
# Compare MCMC results samples to truth    #
############################################
function get_percentile_df(truth, sample_list)

    STDERR_COLUMNS = Set([
        :log_flux_r_stderr,
        :color_ug_stderr,
        :color_gr_stderr,
        :color_ri_stderr,
        :color_iz_stderr,
    ])

    # store these values --- easier that way
    truth[:log_flux_r] = log.(truth[:flux_r_nmgy])

    # go through each truth source, match it to a sample in the results list
    _, _, sample_summary = MCMC.consolidate_samples(sample_list; summarize=mean)
    #sample_summary[:dec] -= .396 / 3600.
    #sample_summary[:ra]  -= .396 / 3600.
    samp_pos = convert(Array, sample_summary[[:ra, :dec]])

    function get_percentile_column(column_name)
        percentiles, names = []  , []
        for i in 1:size(truth, 1)

            # match this truth position to the appropriate set of samples
            pos_i = convert(Array, truth[[:ra, :dec]][i, :])
            dists = sum(broadcast(-, samp_pos, pos_i).^2, 2)
            near_dist, near_idx = findmin(dists)

            # grab samples of appropriate source type
            sdict = results[near_idx]
            pstar = exp(sdict["ave_pstar"])
            samps = (pstar > .5) ? sdict["star_samples"] : sdict["gal_samples"]

            # truth percentile
            true_val = truth[column_name][i]
            true_percentile = sum(true_val .> sort(samps[column_name])) / float(size(samps, 1))

            push!(percentiles, true_percentile)
            push!(names, column_name)
        end
        return percentiles, names
    end

    # assemble dataframe of percentile values
    column_names = [:log_flux_r, :color_ug, :color_gr, :color_ri, :color_iz]
    res = [get_percentile_column(c) for c in column_names]
    pers = vcat([r[1] for r in res]...)
    names = vcat([r[2] for r in res]...)
    return DataFrame([names, pers], [:name, :percentile])
end

###############################################
# Compute scores for comparison data frames   #
###############################################
scores = AccuracyBenchmark.score_predictions(truth, [vbdf, results_df])
println("Comparing vb (first) to MCMC (second)")
println(repr(scores))

# uncertainties
vb_udf = AccuracyBenchmark.get_uncertainty_df(truth, vbdf);
vb_uscore = AccuracyBenchmark.score_uncertainty(vb_udf);

mc_udf = AccuracyBenchmark.get_uncertainty_df(truth, results_df);
mc_uscore = AccuracyBenchmark.score_uncertainty(mc_udf);
mc_perdf = get_percentile_df(truth, results)

println("VB Uncertainty")
println(vb_uscore)
println("MC Uncertainty")
println(mc_uscore)

# output to output dir
CSV.write(joinpath(output_dir, "uscore_vb.csv"), vb_uscore)
CSV.write(joinpath(output_dir, "uscore_mc.csv"), mc_uscore)
CSV.write(joinpath(output_dir, "udf_vb.csv"), vb_udf)
CSV.write(joinpath(output_dir, "udf_mc.csv"), mc_udf)
CSV.write(joinpath(output_dir, "perdf_mc.csv"), mc_perdf)

########################################################################
# Compute bootstrap resampled uncertainties for prediction estimates   #
########################################################################
function score_predictions_bootstrap(truth::DataFrame,
                                     prediction_dfs::Vector{DataFrame};
                                     num_bootstrap=1000,
                                     lo=2.5,
                                     hi=97.5)
    matched_truth, matched_prediction_dfs = AccuracyBenchmark.match_catalogs(truth, prediction_dfs)

    # in each sample, resample matched truth and matched pred dataframe rows
    scoredfs = []
    for bsamp in 1:num_bootstrap
        nsamps = size(matched_truth,1)
        bidx = [rand(1:nsamps) for n in 1:nsamps]
        mtruth = matched_truth[bidx,:]
        mpreds = [m[bidx,:] for m in matched_prediction_dfs]

        error_dfs = [AccuracyBenchmark.get_error_df(mtruth, predictions) for predictions in mpreds]
        @assert length(prediction_dfs) <= 2
        sdf = AccuracyBenchmark.get_scores_df(
            matched_truth,
            error_dfs[1],
            length(error_dfs) > 1 ? Nullable(error_dfs[2]) : Nullable{DataFrame}(),
        )
        sdf[:sample_i] = ones(Int64, size(sdf,1))*bsamp
        push!(scoredfs, sdf)
    end

    # assemble hi and lo values in the bootstrap samples
    sdf = scoredfs[1]
    boot_dict = Dict(:first_lo => [],
                     :first_hi => [])
                     #:field => [])
    if :second in names(sdf)
        boot_dict[:second_lo] = []
        boot_dict[:second_hi] = []
    end

    nrows = size(sdf, 1)
    for r in 1:nrows
        first_vals  = [sdf[r,:first] for sdf in scoredfs]
        push!(boot_dict[:first_lo], percentile(first_vals, lo))
        push!(boot_dict[:first_hi], percentile(first_vals, hi))
        if :second in names(sdf)
            second_vals = [sdf[r,:second] for sdf in scoredfs]
            push!(boot_dict[:second_lo], percentile(second_vals, lo))
            push!(boot_dict[:second_hi], percentile(second_vals, hi))
        end
        #push!(boot_dict[:field], sdf[r,:field])
    end

    scoremean = AccuracyBenchmark.score_predictions(truth, prediction_dfs)
    score_cis = DataFrame(boot_dict)
    return hcat(scoremean, score_cis), vcat(scoredfs...)
end

lo, hi = 2.5, 97.5
scoredf, scoredf_boot = score_predictions_bootstrap(truth, [vbdf, results_df]);
println("Bootstrap scoredf:")
println(scoredf)
CSV.write(joinpath(output_dir, "scoredf.csv"), scoredf)
CSV.write(joinpath(output_dir, "scoredf_boot.csv"), scoredf_boot)

##############################
# look at star ROC curve     #
##############################
#matched_truth, matched_prediction_dfs = AccuracyBenchmark.match_catalogs(
#    truth, [vbdf, results_df, photodf]) #, tol=.396/3600*4)
#@printf " ... found %d matched sources for comparison \n" size(matched_truth, 1)

# look at ROC curve for pstar
true_star = matched_truth[:is_star]
pstar_vb  = matched_prediction_dfs[1][:is_star]
pstar_mc  = matched_prediction_dfs[2][:is_star]
pstar_df = DataFrame(Dict(:true_star => true_star,
                          :pstar_vb => pstar_vb,
                          :pstar_mc => pstar_mc))
CSV.write(joinpath(output_dir, "pstardf.csv"), pstar_df)
CSV.write(joinpath(output_dir, "matched_truth.csv"), matched_truth)
CSV.write(joinpath(output_dir, "matched_vb.csv"), matched_prediction_dfs[1])
CSV.write(joinpath(output_dir, "matched_mc.csv"), matched_prediction_dfs[2])
CSV.write(joinpath(output_dir, "matched_photo.csv"), matched_prediction_dfs[3])