text:
export diageye

using StatsFuns: logistic
using StatsFuns: softmax, softmax!

import LinearAlgebra
import LoopVectorization: @turbo

diageye(::Type{T}, n::Int) where { T <: Real } = Matrix{T}(I, n, n)
diageye(n::Int)                                = diageye(Float64, n)

function normalize_sum(x::Array{Float64,1})
    x ./ sum(x)
end

const sigmoid = logistic

dtanh(x) = 1 - tanh(x)^2
"""
    negate_inplace!(A)

Returns `-A`, modifying and reusing `A` storage if possible.

See also: [`mul_inplace!`](@ref)
"""
function negate_inplace! end

negate_inplace!(A::AbstractArray) = -A
negate_inplace!(A::Array)         = map!(-, A, A)
negate_inplace!(A::Real)          = -A

"""
    mul_inplace!(alpha, A)

Returns `alpha * A`, modifying and reusing `A` storage if possible.

See also: [`negate_inplace!`](@ref)
"""
function mul_inplace! end

mul_inplace!(alpha, A::AbstractArray) = alpha * A
mul_inplace!(alpha, A::Array)         = lmul!(alpha, A)
mul_inplace!(alpha, A::Real)          = alpha * A

"""
    rank1update(A, x)
    rank1update(A, x, y)

Helper function for A + x * y'. Uses optimised BLAS version for AbstractFloats and fallbacks to a generic implementation in case of differentiation
"""
function rank1update end

rank1update(A::AbstractMatrix, x::AbstractVector)                    = rank1update(eltype(A), eltype(x), eltype(x), A, x, x)
rank1update(A::AbstractMatrix, x::AbstractVector, y::AbstractVector) = rank1update(eltype(A), eltype(x), eltype(y), A, x, y)

rank1update(A::Real, x::Real)          = rank1update(A, x, x)
rank1update(A::Real, x::Real, y::Real) = A + x * y

function rank1update(::Type{ T }, ::Type{ T }, ::Type{T}, A::Matrix, x::Vector, y::Vector) where { T <: AbstractFloat }
    return LinearAlgebra.BLAS.ger!(one(T), x, y, copy(A))
end

function rank1update(::Type{ T1 }, ::Type{ T2 }, ::Type{ T3 }, A::AbstractMatrix, x::AbstractVector, y::AbstractVector) where { T1 <: Real, T2 <: Real , T3 <: Real }
    return A + x * y'
end

"""
    mul_trace(A, B)

Computes tr(A * B) wihtout allocating A * B.
"""
function mul_trace end

mul_trace(A::Real, B::Real) = A * B

function mul_trace(A::AbstractMatrix, B::AbstractMatrix)
    sA, sB = size(A), size(B)
    @assert (sA === sB) && (length(sA) === 2) && (first(sA) === last(sA))
    result = zero(promote_type(eltype(A), eltype(B)))
    @turbo for i in 1:first(sA), j in 1:first(sA)
        @inbounds result += A[i, j] * B[j, i]
    end
    return result
end


text:
# ========= GaussSBP approximation types ============
# Note: we define type aliases outside of the @muladd block to avoid Revise breaking when code
# inside the @muladd block is edited. See https://github.com/trixi-framework/Trixi.jl/issues/801
# for more details.

# GaussSBP ApproximationType: e.g., Gauss nodes on quads/hexes
struct GaussSBP end

function tensor_product_quadrature(element_type::Line, r1D, w1D)
  return r1D, w1D
end

function tensor_product_quadrature(element_type::Quad, r1D, w1D)
  sq, rq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D))
  ws, wr = vec.(StartUpDG.NodesAndModes.meshgrid(w1D))
  wq = wr .* ws
  return rq, sq, wq
end

function tensor_product_quadrature(element_type::Hex, r1D, w1D)
  rq, sq, tq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D, r1D, r1D))
  wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(w1D, w1D, w1D))
  wq = wr .* ws .* wt
  return rq, sq, tq, wq
end

# type parameters for `TensorProductFaceOperator`.
abstract type AbstractGaussOperator end
struct Interpolation <: AbstractGaussOperator end
# - `Projection{ScaleByFaceWeights=Static.False()}` corresponds to the operator `Pf = M \ Vf'`,
#   which is used in `VolumeIntegralFluxDifferencing`.
# - `Projection{ScaleByFaceWeights=Static.True()}` corresponds to the quadrature-based lifting
#   operator `LIFT = M \ (Vf' * diagm(rd.wf))`, which is used in `SurfaceIntegralWeakForm`
struct Projection{ScaleByFaceWeights}  <: AbstractGaussOperator end

#   TensorProductGaussFaceOperator{Tmat, Ti}
#
# Data for performing tensor product interpolation from volume nodes to face nodes.
struct TensorProductGaussFaceOperator{NDIMS, OperatorType <: AbstractGaussOperator,
                                      Tmat, Tweights, Tfweights, Tindices}
  interp_matrix_gauss_to_face_1d::Tmat
  inv_volume_weights_1d::Tweights
  face_weights::Tfweights
  face_indices_tensor_product::Tindices
  nnodes_1d::Int
  nfaces::Int
end

# constructor for a 2D operator
function TensorProductGaussFaceOperator(operator::AbstractGaussOperator,
                                        dg::DGMulti{2, Quad, GaussSBP})
  rd = dg.basis

  rq1D, wq1D = StartUpDG.gauss_quad(0, 0, polydeg(dg))
  interp_matrix_gauss_to_face_1d = polynomial_interpolation_matrix(rq1D, [-1; 1])

  nnodes_1d = length(rq1D)

  # Permutation of indices in a tensor product form
  indices = reshape(1:length(rd.rf), nnodes_1d, rd.Nfaces)
  face_indices_tensor_product = zeros(Int, 2, nnodes_1d, ndims(rd.elementType))
  for i in 1:nnodes_1d # loop over nodes in one face
    face_indices_tensor_product[:, i, 1] .= indices[i, 1:2]
    face_indices_tensor_product[:, i, 2] .= indices[i, 3:4]
  end

  T_op = typeof(operator)
  Tm = typeof(interp_matrix_gauss_to_face_1d)
  Tw = typeof(inv.(wq1D))
  Tf = typeof(rd.wf)
  Ti = typeof(face_indices_tensor_product)
  return TensorProductGaussFaceOperator{2, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d,
                                                                 inv.(wq1D), rd.wf,
                                                                 face_indices_tensor_product,
                                                                 nnodes_1d, rd.Nfaces)
end

# constructor for a 3D operator
function TensorProductGaussFaceOperator(operator::AbstractGaussOperator,
                                        dg::DGMulti{3, Hex, GaussSBP})
  rd = dg.basis

  rq1D, wq1D = StartUpDG.gauss_quad(0, 0, polydeg(dg))
  interp_matrix_gauss_to_face_1d = polynomial_interpolation_matrix(rq1D, [-1; 1])

  nnodes_1d = length(rq1D)

  # Permutation of indices in a tensor product form
  indices = reshape(1:length(rd.rf), nnodes_1d, nnodes_1d, rd.Nfaces)
  face_indices_tensor_product = zeros(Int, 2, nnodes_1d, nnodes_1d, ndims(rd.elementType))
  for j in 1:nnodes_1d, i in 1:nnodes_1d # loop over nodes in one face
    face_indices_tensor_product[:, i, j, 1] .= indices[i, j, 1:2]
    face_indices_tensor_product[:, i, j, 2] .= indices[i, j, 3:4]
    face_indices_tensor_product[:, i, j, 3] .= indices[i, j, 5:6]
  end

  T_op = typeof(operator)
  Tm = typeof(interp_matrix_gauss_to_face_1d)
  Tw = typeof(inv.(wq1D))
  Tf = typeof(rd.wf)
  Ti = typeof(face_indices_tensor_product)
  return TensorProductGaussFaceOperator{3, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d,
                                                                 inv.(wq1D), rd.wf,
                                                                 face_indices_tensor_product,
                                                                 nnodes_1d, rd.Nfaces)
end

# specialize behavior of `mul_by!(A::TensorProductGaussFaceOperator)`
@inline function mul_by!(A::TensorProductGaussFaceOperator)
  return (out, x) -> tensor_product_gauss_face_operator!(out, A, x)
end

@inline function tensor_product_gauss_face_operator!(out::AbstractMatrix,
                                                     A::TensorProductGaussFaceOperator,
                                                     x::AbstractMatrix)
  @threaded for col in Base.OneTo(size(out, 2))
    tensor_product_gauss_face_operator!(view(out, :, col), A, view(x, :, col))
  end
end

# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Interpolates values from volume Gauss nodes to face nodes on one element.
@inline function tensor_product_gauss_face_operator!(out::AbstractVector,
                                                     A::TensorProductGaussFaceOperator{2, Interpolation},
                                                     x::AbstractVector)

  @unpack interp_matrix_gauss_to_face_1d, face_indices_tensor_product = A
  @unpack nnodes_1d, nfaces = A

  fill!(out, zero(eltype(out)))

  # for 2D GaussSBP nodes, the indexing is first in x, then in y
  x = reshape(x, nnodes_1d, nnodes_1d)

  # interpolation in the x-direction
  @turbo for i in Base.OneTo(nnodes_1d) # loop over nodes in a face
    index_left  = face_indices_tensor_product[1, i, 1]
    index_right = face_indices_tensor_product[2, i, 1]
    for jj in Base.OneTo(nnodes_1d)      # loop over "line" of volume nodes
      out[index_left]  = out[index_left]  + interp_matrix_gauss_to_face_1d[1, jj] * x[jj, i]
      out[index_right] = out[index_right] + interp_matrix_gauss_to_face_1d[2, jj] * x[jj, i]
    end
  end

  # interpolation in the y-direction
  @turbo for i in Base.OneTo(nnodes_1d) # loop over nodes in a face
    index_left  = face_indices_tensor_product[1, i, 2]
    index_right = face_indices_tensor_product[2, i, 2]
    for jj in Base.OneTo(nnodes_1d)               # loop over "line" of volume nodes
      out[index_left]  = out[index_left]  + interp_matrix_gauss_to_face_1d[1, jj] * x[i, jj]
      out[index_right] = out[index_right] + interp_matrix_gauss_to_face_1d[2, jj] * x[i, jj]
    end
  end
end

# Interpolates values from volume Gauss nodes to face nodes on one element.
@inline function tensor_product_gauss_face_operator!(out::AbstractVector,
                                                     A::TensorProductGaussFaceOperator{3, Interpolation},
                                                     x::AbstractVector)

  @unpack interp_matrix_gauss_to_face_1d, face_indices_tensor_product = A
  @unpack nnodes_1d, nfaces = A

  fill!(out, zero(eltype(out)))

  # for 3D GaussSBP nodes, the indexing is first in y, then x, then z.
  x = reshape(x, nnodes_1d, nnodes_1d, nnodes_1d)

  # interpolation in the y-direction
  @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face
    index_left  = face_indices_tensor_product[1, i, j, 2]
    index_right = face_indices_tensor_product[2, i, j, 2]
    for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes
      out[index_left]  = out[index_left]  + interp_matrix_gauss_to_face_1d[1, jj] * x[jj, i, j]
      out[index_right] = out[index_right] + interp_matrix_gauss_to_face_1d[2, jj] * x[jj, i, j]
    end
  end

  # interpolation in the x-direction
  @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face
    index_left  = face_indices_tensor_product[1, i, j, 1]
    index_right = face_indices_tensor_product[2, i, j, 1]
    for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes
      out[index_left]  = out[index_left]  + interp_matrix_gauss_to_face_1d[1, jj] * x[i, jj, j]
      out[index_right] = out[index_right] + interp_matrix_gauss_to_face_1d[2, jj] * x[i, jj, j]
    end
  end

  # interpolation in the z-direction
  @turbo for i in Base.OneTo(nnodes_1d), j in Base.OneTo(nnodes_1d) # loop over nodes in a face
    index_left  = face_indices_tensor_product[1, i, j, 3]
    index_right = face_indices_tensor_product[2, i, j, 3]
    for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes
      # The ordering (i,j) -> (j,i) needs to be reversed for this last face.
      # This is due to way we define face nodes for Hex() types in StartUpDG.jl.
      out[index_left]  = out[index_left]  + interp_matrix_gauss_to_face_1d[1, jj] * x[j, i, jj]
      out[index_right] = out[index_right] + interp_matrix_gauss_to_face_1d[2, jj] * x[j, i, jj]
    end
  end
end

# Projects face node values to volume Gauss nodes on one element.
@inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector,
                                                     A::TensorProductGaussFaceOperator{2, Projection{ApplyFaceWeights}},
                                                     x::AbstractVector) where {ApplyFaceWeights}

  @unpack interp_matrix_gauss_to_face_1d, face_indices_tensor_product = A
  @unpack inv_volume_weights_1d, nnodes_1d, nfaces = A

  fill!(out_vec, zero(eltype(out_vec)))

  # for 2D GaussSBP nodes, the indexing is first in x, then y
  out = reshape(out_vec, nnodes_1d, nnodes_1d)

  if ApplyFaceWeights == true
    @turbo for i in eachindex(x)
      x[i] = x[i] * A.face_weights[i]
    end
  end

  # interpolation in the x-direction
  @turbo for i in Base.OneTo(nnodes_1d) # loop over face nodes
    index_left  = face_indices_tensor_product[1, i, 1]
    index_right = face_indices_tensor_product[2, i, 1]
    for jj in Base.OneTo(nnodes_1d) # loop over a line of volume nodes
      out[jj, i] = out[jj, i] + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left]
      out[jj, i] = out[jj, i] + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right]
    end
  end

  # interpolation in the y-direction
  @turbo for i in Base.OneTo(nnodes_1d)
    index_left  = face_indices_tensor_product[1, i, 2]
    index_right = face_indices_tensor_product[2, i, 2]
    # loop over a line of volume nodes
    for jj in Base.OneTo(nnodes_1d)
      out[i, jj] = out[i, jj] + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left]
      out[i, jj] = out[i, jj] + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right]
    end
  end

  # apply inv(M)
  @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d)
    out[i, j] = out[i, j] * inv_volume_weights_1d[i] * inv_volume_weights_1d[j]
  end
end

# Interpolates values from volume Gauss nodes to face nodes on one element.
@inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector,
                                                     A::TensorProductGaussFaceOperator{3, Projection{ApplyFaceWeights}},
                                                     x::AbstractVector) where {ApplyFaceWeights}

  @unpack interp_matrix_gauss_to_face_1d, face_indices_tensor_product = A
  @unpack inv_volume_weights_1d, nnodes_1d, nfaces = A

  fill!(out_vec, zero(eltype(out_vec)))

  # for 3D GaussSBP nodes, the indexing is first in y, then x, then z.
  out = reshape(out_vec, nnodes_1d, nnodes_1d, nnodes_1d)

  if ApplyFaceWeights == true
    @turbo for i in eachindex(x)
      x[i] = x[i] * A.face_weights[i]
    end
  end

  # interpolation in the y-direction
  @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face
    index_left  = face_indices_tensor_product[1, i, j, 2]
    index_right = face_indices_tensor_product[2, i, j, 2]
    for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes
      out[jj, i, j] = out[jj, i, j] + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left]
      out[jj, i, j] = out[jj, i, j] + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right]
    end
  end

  # interpolation in the x-direction
  @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face
    index_left  = face_indices_tensor_product[1, i, j, 1]
    index_right = face_indices_tensor_product[2, i, j, 1]
    for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes
      out[i, jj, j] = out[i, jj, j] + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left]
      out[i, jj, j] = out[i, jj, j] + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right]
    end
  end

  # interpolation in the z-direction
  @turbo for i in Base.OneTo(nnodes_1d), j in Base.OneTo(nnodes_1d) # loop over nodes in a face
    index_left  = face_indices_tensor_product[1, i, j, 3]
    index_right = face_indices_tensor_product[2, i, j, 3]
    for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes
      # The ordering (i,j) -> (j,i) needs to be reversed for this last face.
      # This is due to way we define face nodes for Hex() types in StartUpDG.jl.
      out[j, i, jj] = out[j, i, jj] + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left]
      out[j, i, jj] = out[j, i, jj] + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right]
    end
  end

  # apply inv(M)
  @turbo for k in Base.OneTo(nnodes_1d), j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d)
    out[i, j, k] = out[i, j, k] * inv_volume_weights_1d[i] * inv_volume_weights_1d[j] * inv_volume_weights_1d[k]
  end
end

# Specialized constructor for GaussSBP approximation type on quad elements. Restricting to
# VolumeIntegralFluxDifferencing for now since there isn't a way to exploit this structure
# for VolumeIntegralWeakForm yet.
function DGMulti(element_type::Union{Quad, Hex},
                 approximation_type::GaussSBP,
                 volume_integral::VolumeIntegralFluxDifferencing,
                 surface_integral=SurfaceIntegralWeakForm(surface_flux);
                 polydeg::Integer,
                 surface_flux=flux_central,
                 kwargs...)

  # explicitly create tensor product Gauss quadrature rule with polydeg+1 points
  r1D, w1D = StartUpDG.gauss_quad(0, 0, polydeg)
  gauss_rule_vol = tensor_product_quadrature(element_type, r1D, w1D)
  gauss_rule_face = tensor_product_quadrature(StartUpDG.face_type(element_type), r1D, w1D)

  rd = RefElemData(element_type, Polynomial(), polydeg,
                   quad_rule_vol=gauss_rule_vol,
                   quad_rule_face=gauss_rule_face,
                   kwargs...)

  # Since there is no dedicated GaussSBP approximation type implemented in StartUpDG, we simply
  # initialize `rd = RefElemData(...)` with the appropriate quadrature rules and modify the
  # rd.approximationType manually so we can dispatch on the `GaussSBP` type.
  # This uses the Setfield @set macro, which behaves similarly to `Trixi.remake`.
  rd_gauss = @set rd.approximationType = GaussSBP()

  # We will modify the face interpolation operator of rd_gauss later, but want to do so only after
  # the mesh is initialized, since the face interpolation operator is used for that.
  return DG(rd_gauss, nothing #= mortar =#, surface_integral, volume_integral)
end

# For now, this is mostly the same as `create_cache` for DGMultiFluxDiff{<:Polynomial}.
# In the future, we may modify it so that we can specialize additional parts of GaussSBP() solvers.
function create_cache(mesh::VertexMappedMesh, equations,
                      dg::DGMultiFluxDiff{<:GaussSBP, <:Union{Quad, Hex}}, RealT, uEltype)

  rd = dg.basis
  @unpack md = mesh

  cache = invoke(create_cache, Tuple{typeof(mesh), typeof(equations), DGMultiFluxDiff, typeof(RealT), typeof(uEltype)},
                 mesh, equations, dg, RealT, uEltype)

  # for change of basis prior to the volume integral and entropy projection
  r1D, _ = StartUpDG.gauss_lobatto_quad(0, 0, polydeg(dg))
  rq1D, _ = StartUpDG.gauss_quad(0, 0, polydeg(dg))
  interp_matrix_lobatto_to_gauss_1D = polynomial_interpolation_matrix(r1D, rq1D)
  interp_matrix_gauss_to_lobatto_1D = polynomial_interpolation_matrix(rq1D, r1D)
  NDIMS = ndims(rd.elementType)
  interp_matrix_lobatto_to_gauss = SimpleKronecker(NDIMS, interp_matrix_lobatto_to_gauss_1D, uEltype)
  interp_matrix_gauss_to_lobatto = SimpleKronecker(NDIMS, interp_matrix_gauss_to_lobatto_1D, uEltype)
  inv_gauss_weights = inv.(rd.wq)

  # specialized operators to perform tensor product interpolation to faces for Gauss nodes
  interp_matrix_gauss_to_face = TensorProductGaussFaceOperator(Interpolation(), dg)
  Pf = TensorProductGaussFaceOperator(Projection{Static.False()}(), dg)

  # `LIFT` matrix for Gauss nodes - this is equivalent to `Pf` scaled by `diagm(rd.wf)`,
  # where `rd.wf` are Gauss node face quadrature weights.
  gauss_LIFT = TensorProductGaussFaceOperator(Projection{Static.True()}(), dg)

  nvars = nvariables(equations)
  rhs_volume_local_threaded   = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)  for _ in 1:Threads.nthreads()]
  gauss_volume_local_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)  for _ in 1:Threads.nthreads()]

  return (; cache..., Pf, gauss_LIFT, inv_gauss_weights,
         rhs_volume_local_threaded, gauss_volume_local_threaded,
         interp_matrix_lobatto_to_gauss, interp_matrix_gauss_to_lobatto,
         interp_matrix_gauss_to_face)
end

# TODO: DGMulti. Address hard-coding of `entropy2cons!` and `cons2entropy!` for this function.
function entropy_projection!(cache, u, mesh::VertexMappedMesh, equations, dg::DGMultiFluxDiff{<:GaussSBP})

  rd = dg.basis
  @unpack Vq = rd
  @unpack VhP, entropy_var_values, u_values = cache
  @unpack projected_entropy_var_values, entropy_projected_u_values = cache
  @unpack interp_matrix_lobatto_to_gauss, interp_matrix_gauss_to_face = cache

  @threaded for e in eachelement(mesh, dg, cache)
    apply_to_each_field(mul_by!(interp_matrix_lobatto_to_gauss), view(u_values, :, e), view(u, :, e))
  end

  # transform quadrature values to entropy variables
  cons2entropy!(entropy_var_values, u_values, equations)

  volume_indices = Base.OneTo(rd.Nq)
  face_indices = (rd.Nq + 1):(rd.Nq + rd.Nfq)

  # Interpolate volume Gauss nodes to Gauss face nodes (note the layout of
  # `projected_entropy_var_values = [vol pts; face pts]`).
  entropy_var_face_values = view(projected_entropy_var_values, face_indices, :)
  apply_to_each_field(mul_by!(interp_matrix_gauss_to_face), entropy_var_face_values, entropy_var_values)

  # directly copy over volume values (no entropy projection required)
  entropy_projected_volume_values = view(entropy_projected_u_values, volume_indices, :)
  @threaded for i in eachindex(u_values)
    entropy_projected_volume_values[i] = u_values[i]
  end

  # transform entropy to conservative variables on face values
  entropy_projected_face_values = view(entropy_projected_u_values, face_indices, :)
  entropy2cons!(entropy_projected_face_values, entropy_var_face_values, equations)

  return nothing
end

# Assumes cache.flux_face_values is already computed.
# Enables tensor product evaluation of `LIFT::TensorProductGaussFaceOperator`.
function calc_surface_integral!(du, u, surface_integral::SurfaceIntegralWeakForm,
                                mesh::VertexMappedMesh, equations,
                                dg::DGMultiFluxDiff{<:GaussSBP}, cache)
  @unpack gauss_volume_local_threaded, rhs_volume_local_threaded = cache
  @unpack interp_matrix_gauss_to_lobatto, gauss_LIFT = cache

  @threaded for e in eachelement(mesh, dg, cache)

    # applies LIFT matrix, output is stored at Gauss nodes
    gauss_volume_local = gauss_volume_local_threaded[Threads.threadid()]
    apply_to_each_field(mul_by!(gauss_LIFT), gauss_volume_local, view(cache.flux_face_values, :, e))

    # interpolate result back to Lobatto nodes for ease of analysis, visualization
    rhs_volume_local = rhs_volume_local_threaded[Threads.threadid()]
    apply_to_each_field(mul_by!(interp_matrix_gauss_to_lobatto), rhs_volume_local, gauss_volume_local)

    for i in eachindex(rhs_volume_local)
      du[i, e] = du[i, e] + rhs_volume_local[i]
    end
  end
end

function calc_volume_integral!(du, u, mesh::VertexMappedMesh,
                               have_nonconservative_terms::Val{false}, equations,
                               volume_integral, dg::DGMultiFluxDiff{<:GaussSBP},
                               cache)

  @unpack entropy_projected_u_values = cache
  @unpack fluxdiff_local_threaded, rhs_local_threaded, rhs_volume_local_threaded = cache

  # After computing the volume integral, we transform back to Lobatto nodes.
  # This allows us to reuse the other DGMulti routines as-is.
  @unpack interp_matrix_gauss_to_lobatto = cache
  @unpack Ph, Pf, inv_gauss_weights = cache

  rd = dg.basis
  volume_indices = Base.OneTo(rd.Nq)
  face_indices = (rd.Nq + 1):(rd.Nq + rd.Nfq)

  @threaded for e in eachelement(mesh, dg, cache)
    fluxdiff_local = fluxdiff_local_threaded[Threads.threadid()]
    fill!(fluxdiff_local, zero(eltype(fluxdiff_local)))
    u_local = view(entropy_projected_u_values, :, e)

    local_flux_differencing!(fluxdiff_local, u_local, e,
                             have_nonconservative_terms, volume_integral,
                             has_sparse_operators(dg),
                             mesh, equations, dg, cache)

    # convert `fluxdiff_local::Vector{<:SVector}` to `rhs_local::StructArray{<:SVector}`
    # for faster performance when using `apply_to_each_field`.
    rhs_local = rhs_local_threaded[Threads.threadid()]
    for i in Base.OneTo(length(fluxdiff_local))
      rhs_local[i] = fluxdiff_local[i]
    end

    # stores rhs contributions only at Gauss volume nodes
    rhs_volume_local = rhs_volume_local_threaded[Threads.threadid()]

    # Here, we exploit that under a Gauss nodal basis the structure of the projection
    # matrix `Ph = [diagm(1 ./ wq), Pf]` such that `Ph * [u; uf] = (u ./ wq) + Pf * uf`.
    local_volume_flux = view(rhs_local, volume_indices)
    local_face_flux = view(rhs_local, face_indices)

    # initialize rhs_volume_local = Pf * local_face_flux
    apply_to_each_field(mul_by!(Pf), rhs_volume_local, local_face_flux)

    # accumulate volume contributions
    for i in eachindex(rhs_volume_local)
      rhs_volume_local[i] = rhs_volume_local[i] + local_volume_flux[i] * inv_gauss_weights[i]
    end

    # transform rhs back to Lobatto nodes
    apply_to_each_field(mul_by!(interp_matrix_gauss_to_lobatto),
                        view(du, :, e), rhs_volume_local)
  end

end

end # @muladd

ext: jl


text:
module test_julia_comparison

using Test

struct C
    one
    two
end

struct B
    cs::Vector{C}
end

struct A
    bs::Vector{B}
end

for typ in (A, B)
    @eval function Base.:(==)(l::T, r::T) where {T <: $typ}
        for name in fieldnames(T)
            getfield(l, name) != getfield(r, name) && return false
        end
        return true
    end
end

x = [A([B([C(1, 2)])])]
y = [A([B([C(1, 2)])])]

@test x == y

end # module test_julia_comparison


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

"""
    save(term, output, [termopts,] [font,] [size,] [linewidth,] [background,] [handle]) -> nothing

Save current figure (or figure specified by `handle`) using the specified `term`. Optionally,
the font, size, linewidth, and background may be specified as arguments.
"""
function save(; term::String,
                output::String,
                handle::Int = gnuplot_state.current,
                font        = "",
                size        = "",
                linewidth   = 0,
                background  = "",
                saveopts    = config[:saveopts])

    debug("term = $term, output = $output, saveopts = $saveopts", "save")

    # process arguments
    isempty(output) && throw(DomainError("Please specify an output filename."))
    isempty(term) && throw(DomainError("Please specify a file format"))

    h = findfigure(handle)
    h == 0 && throw(DomainError(h, "requested figure does not exist."))
    fig = gnuplot_state.figs[h]

    (term == "pdf" || term == :pdf) && (term = "pdfcairo")
    (term == "png" || term == :png) && (term = "pngcairo")
    (term == "eps" || term == :eps) && (term = "epscairo")

    # create print configuration
    pc = "set term $term "
    if isempty(saveopts)
        !isempty(font) && (pc *= "font '$font' ")
        !isempty(size) && (pc *= "size $size ")
        linewidth > 0 && (pc *= "linewidth $linewidth ")
        !isempty(background) && (pc *= "background '$background' ")
    else
        pc *= saveopts
    end

    # send gnuplot commands
    llplot(fig, printstring=(pc,output))

    return nothing
end


text:
# Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
# [@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

# 12 steps to Navier–Stokes
# =====
# ***

# You see where this is going ... we'll do 2D diffusion now and next we will combine steps 6 and 7 to solve Burgers' equation. So make sure your previous steps work well before continuing.

# Step 7: 2D Diffusion
# ----
# ***

# And here is the 2D-diffusion equation:

# $$\frac{\partial u}{\partial t} = \nu \frac{\partial ^2 u}{\partial x^2} + \nu \frac{\partial ^2 u}{\partial y^2}$$

# You will recall that we came up with a method for discretizing second order derivatives in Step 3, when investigating 1-D diffusion.  We are going to use the same scheme here, with our forward difference in time and two second-order derivatives.

# $$\frac{u_{i,j}^{n+1} - u_{i,j}^n}{\Delta t} = \nu \frac{u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n}{\Delta x^2} + \nu \frac{u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n}{\Delta y^2}$$

# Once again, we reorganize the discretized equation and solve for $u_{i,j}^{n+1}$

# $$
# \begin{split}
# u_{i,j}^{n+1} = u_{i,j}^n &+ \frac{\nu \Delta t}{\Delta x^2}(u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n) \\
# &+ \frac{\nu \Delta t}{\Delta y^2}(u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n)
# \end{split}
# $$

using Plots

# ## Variable declarations
nₓ = 31
ny = 31
nₜ = 17
nᵤ = .05
c = 1
Δx = 2 / (nₓ - 1)
Δy = 2 / (ny - 1)
σ = 0.25
Δt = σ * Δx

rangex = range(start = 0, stop = 2, length = nₓ)
rangey = range(start = 0, stop = 2, length = ny)
indices = CartesianIndices((rangex, rangey))

uₙ = ones(ny, nₓ) ##create a 1xn vector of 1's
u  = ones(ny, nₓ)

# ## Assign initial conditions
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
u[indices] .= 2.0
##set hat function I.C. : v(.5<=x<=1 && .5<=y<=1 ) is 2
v[indices] .= 2.0

# TODO PLOT
fig = pyplot.figure()
ax = fig.gca(projection='3d')
X, Y = numpy.meshgrid(x, y)
surf = ax.plot_surface(X, Y, u, rstride=1, cstride=1, cmap=cm.viridis,
        linewidth=0, antialiased=False)

ax.set_xlim(0, 2)
ax.set_ylim(0, 2)
ax.set_zlim(1, 2.5)

ax.set_xlabel('$x$')
ax.set_ylabel('$y$');

# $$
# \begin{split}
# u_{i,j}^{n+1} = u_{i,j}^n &+ \frac{\nu \Delta t}{\Delta x^2}(u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n) \\
# &+ \frac{\nu \Delta t}{\Delta y^2}(u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n)
# \end{split}
# $$

#Run through nt timesteps
# TODO wrap in a let block maybe?
function diffuse(nt)
    rangex = range(start = 0, stop = 2, length = nₓ)
    rangey = range(start = 0, stop = 2, length = ny)
    indices = CartesianIndices((rangex, rangey))
    u[indices] = 2.0

    for n in 1:nt
        un = copy(n)
        # TODO check indexing
        u[2:end, 2:end] = (un[2:end,2:end] +
                        nu * Δt / Δx^2 *
                        (un[2:end, 3:end] - 2 * un[2:end, 2:end] + un[2:end, 1:end-1]) +
                        nu * Δt / Δy^2 *
                        (un[3:end, 2:end-1] - 2 * un[2:end, 2:end] + un[begin:end-1, 2:end]))
        u[begin, :] = 1.0
        u[end, :] = 1.0
        u[:, begin] = 1.0
        u[:, end] = 1.0
    end
    u
end

# TODO PLOT
    fig = pyplot.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis,
        linewidth=0, antialiased=True)
    ax.set_zlim(1, 2.5)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$');

## TODO Plot these 3
diffuse(10)

diffuse(14)

diffuse(50)

## Learn More

The video lesson that walks you through the details for Steps 5 to 8 is **Video Lesson 6** on You Tube:

# TODO YOUTUBE
# src="https://www.youtube.com/embed/tUg_dE3NXoY"

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
module MCMCDiagnostics

using Statistics: mean, var
using StatsBase: mean_and_var

export ess_factor_estimate, effective_sample_size, potential_scale_reduction

"""
    autocorrelation(x, k, v = var(x))

Estimate of lag-`k` autocorrelation of `x` from a variogram. `v` is the variance
of `x`, used when supplied.

See Gelman et al (2013), section 11.4.
"""
function autocorrelation(x::AbstractVector, k::Integer, v = var(x))
    x1 = @view(x[1:(end-k)])
    x2 = @view(x[(1+k):end])
    V = sum((x1 .- x2).^2) / length(x1)
    1 - V / (2*v)
end

"""
    ess_factor_estimate(x, v = var(x))

Estimate for effective sample size factor.

Return `τ, K` where `τ` is estimated effective sample size / sample size, and
`K` is the last lag used for autocorrelation estimation.

# Notes

See Gelman et al (2013), section 11.4.

`τ` is capped at 1, this is relevant when the sample has large negative
autocorrelation (happens with HMC/NUTS).

Some implementations (eg Stan) use FFT for autocorrelations, which yields the
whole spectrum. In practice, a <50-100 lags are usually sufficient for
reasonable samplers, so the “naive” version may be more efficient.
"""
function ess_factor_estimate(x::AbstractVector, v = var(x))
    N = length(x)
    τ_inv = 1 + 2 * autocorrelation(x, 1, v)
    K = 2
    while K < N - 2
        Δ = autocorrelation(x, K, v) + autocorrelation(x, K + 1, v)
        if Δ < 0
            break
        else
            τ_inv += 2*Δ
            K += 2
        end
    end
    min(1 / τ_inv, one(τ_inv)), K
end

"""
    effective_sample_size(x, v = var(x))

Effective sample size of vector `x`.

Estimated from autocorrelations. See Gelman et al (2013), section 11.4.

When the variance `v` is supplied, it saves some calculation time.
"""
function effective_sample_size(x::AbstractVector, v = var(x))
    τ, _ = ess_factor_estimate(x, v)
    τ * length(x)
end

"""
    potential_scale_reduction(chains...)

Potential scale reduction factor (for possibly ragged chains).

Also known as R̂. Always ≥ 1 by construction, but values much larger than 1 (say
1.05) indicate poor mixing.

Uses formula from Stan Development Team (2017), section 28.3.
"""
function potential_scale_reduction(chains::AbstractVector...)
    mvs = mean_and_var.(chains)
    W = mean(last.(mvs))
    B = var(first.(mvs))
    √(1 + B / W)
end

end # module


text:
# test_lazy.jl

function test_lazy()
    @testset "$(rpad("Shifted sequences",80))" begin
        test_shift() end
    @testset "$(rpad("Reversed sequences",80))" begin
        test_reverse() end
    @testset "$(rpad("Downsampled sequences",80))" begin
        test_downsample() end
    @testset "$(rpad("Upsampled sequences",80))" begin
        test_upsample() end
    @testset "$(rpad("Modulated sequences",80))" begin
        test_modulate() end
    @testset "$(rpad("Convolutions",80))" begin
        test_convolution() end
end

function test_shift()
    a = CompactSequence([1,2,3,4,5,6])
    test_shifted_sequence(shift(a, 2), a)
    test_shifted_sequence(ShiftedSequence(a, 2), a)
end

function test_shifted_sequence(s, f)
    n = sequence_shift(s)

    @test s[4] == f[4-n]

    @test iscompact(s) == iscompact(f)

    if iscompact(f)
        @test first(nzrange(s)) == first(nzrange(f))+n
        @test last(nzrange(s)) == last(nzrange(f))+n
    end

    # Check that an additional shift is added to the existing shift
    if typeof(s) <: ShiftedSequence
        @test sequence_shift(shift(s, 4)) == n+4
    end
end

function test_reverse()
    a = CompactSequence([1,2,3,4,5,6])
    test_reversed_sequence(reverse(a), a)
    test_reversed_sequence(ReversedSequence(a), a)
end

function test_reversed_sequence(s, f)
    @test s[4] == f[-4]
    @test iscompact(s) == iscompact(f)

    if iscompact(f)
        @test first(nzrange(s)) == -last(nzrange(f))
        @test last(nzrange(s)) == -first(nzrange(f))
    end

    if typeof(s) <: ReversedSequence
        @test reverse(s) == f
    end
end

function test_downsample()
    a = CompactSequence([1,2,3,4,5,6])
    test_downsampled_sequence(downsample(a, 2), a, 2)
    test_downsampled_sequence(DownsampledSequence(a, 2), a, 2)
end

function test_downsampled_sequence(s, f, M)
    @test samplefactor(s) == M
    for k in (-2, 0, 1, 4)
        @test s[k] == f[M*k]
    end
end

function test_upsample()
    a = CompactSequence([1,2,3,4,5,6])
    test_upsampled_sequence(upsample(a, 2), a, 2)
    test_upsampled_sequence(UpsampledSequence(a, 2), a, 2)
end

function test_upsampled_sequence(s, f, M)
    @test samplefactor(s) == M
    for k in (-2, 0, 1, 4)
        @test s[M*k] == f[k]
    end
    for m in 1:M-1
        @test s[m] == zero(eltype(s))
    end
end

function test_modulate()
    a = CompactSequence([1.0,2.0,3.0,4.0,5.0,6.0])
    test_modulated_sequence(modulate(a), a)
    test_modulated_sequence(ModulatedSequence(a), a)
end

function test_modulated_sequence(s, f)
    for k in -2:2
        @test s[k] == (-one(eltype(s)) )^k * f[k]
    end

    @test iscompact(s) == iscompact(f)

    if iscompact(f)
        @test first(nzrange(s)) == first(nzrange(f))
        @test last(nzrange(s)) == last(nzrange(f))
    end
end

function test_convolution()
    a = CompactSequence([1.0,2.0,3.0,4.0,5.0,6.0])
    b = CompactSequence(OffsetArray([4.1,2.2,3.3], -2:0))
    @test a*a == convolve(a, a)
    test_convolved_sequences(a*a, a, a)
    test_convolved_sequences(a*b, a, b)
end

function test_convolved_sequences(s, f, g)
    n = 4
    z = zero(eltype(f))
    @test s[n] == sum([f[k]*g[n-k] for k in nzrange(f)])

    @test iscompact(s) == iscompact(f) && iscompact(g)
end


text:
@testset "Joint Life" begin

    @testset "Joint Last Survivor unknown status" begin

        @testset "ALMCR §9.4" begin
            ℓ₁ = [43302,42854,42081,41351,40050]
            ℓ₂ = [47260,47040,46755,46500,46227]

            ps₁ =  ℓ₁ ./ ℓ₁[1]
            qs₁ = [1 - ps₁[t] / ps₁[t - 1] for t in 2:5 ]

            ps₂ =  ℓ₂ ./ ℓ₂[1]
            qs₂ = [1 - ps₂[t] / ps₂[t - 1] for t in 2:5 ]

            m1 = UltimateMortality(qs₁, start_age=65)
            m2 = UltimateMortality(qs₂, start_age=60)

            @test decrement(m1, 65, 66) == 1 - ℓ₁[2] / ℓ₁[1]
            @test decrement(m2, 60, 61) == 1 - ℓ₂[2] / ℓ₂[1]

            l1 = SingleLife(mort = m1, issue_age = 65)
            l2 = SingleLife(mort = m2, issue_age = 60)

            jl = JointLife(lives=(l1, l2), contingency = LastSurvivor(), joint_assumption=Frasier())

            @test isapprox( survival(jl, 2) , 0.9997, atol = 1e-4)

            ins = LifeContingency(jl,Yields.Constant(0.05))
            ins_l1 = LifeContingency(jl.lives[1],Yields.Constant(0.05))
            ins_l2 = LifeContingency(jl.lives[2],Yields.Constant(0.05))
            # problem 9.1.f
            @test isapprox( present_value(AnnuityDue(ins,n=5)), 4.5437, atol = 1e-4)

            @test isapprox( present_value(AnnuityDue(ins)), 4.5437, atol = 1e-4)
        end

        @testset "CIA tables" begin
            m1 = MortalityTables.table("1986-92 CIA – Male Smoker, ANB")
            m2 = MortalityTables.table("1986-92 CIA – Female Nonsmoker, ANB")
            l1 = SingleLife(mort = m1.ultimate, issue_age = 40)
            l2 = SingleLife(mort = m2.ultimate, issue_age = 37)

            jl = JointLife(lives=(l1, l2), contingency=LastSurvivor(), joint_assumption=Frasier())

            @testset "independent lives" begin
                for time in 1:40
                    tpx = survival(l1,time)
                    tpy = survival(l2,time)

                    @test survival(jl, time) == tpx + tpy - tpx * tpy
                end
            end

            q_annual = [0.00000141120,0.00000478349,0.00000921100,0.00001508953,0.00002255325,0.00003179413,0.00004327035,0.00005782855,0.00007524178,0.00009695039,0.00012370408,0.00015606419,0.00019612090,0.00024529620,0.00030502659,0.00037768543,0.00046587280,0.00057325907,0.00070238876,0.00085952824,0.00104786439,0.00127420204,0.00154640464,0.00187374100,0.00226186101,0.00272457482,0.00327193743,0.00391984214,0.00468084563,0.00557778027,0.00662135901,0.00783684326,0.00924475551,0.01086696400,0.01272977389,0.01485150445,0.01727277853,0.02000566587,0.02308982834,0.02654184882,0.03039819991,0.03469299691,0.03944737905,0.04469931521,0.05048294677,0.05683522818,0.06379854263,0.07142629169,0.07976664920,0.08888579420,0.09884411368,0.10971505011,0.12159455247,0.13456410066,0.14873755448,0.16421406941,0.18111392905,0.19954137714,0.21962070363,0.24144744481,0.26515968560,0.29099967654,0.31901974549,0.35111309929,0.39733960046,0.47004540325,0.58255000000,0.74852030000,1.00000000000]
            q_cumulative = [0.00000141120,0.00000619469,0.00001540563,0.00003049492,0.00005304748,0.00008483992,0.00012810660,0.00018592774,0.00026115553,0.00035808060,0.00048174038,0.00063772939,0.00083372522,0.00107881691,0.00138351443,0.00176067732,0.00222572986,0.00279771301,0.00349813669,0.00435465818,0.00539795948,0.00666528343,0.00820138084,0.01005975458,0.01229886182,0.01498992747,0.01821281880,0.02206126957,0.02663884979,0.03206804441,0.03847706939,0.04601237389,0.05483175625,0.06510286552,0.07700389465,0.09071177542,0.10641770955,0.12429441828,0.14451430984,0.16722048170,0.19253547997,0.22054884408,0.25129614928,0.28476269870,0.32086998531,0.35946849466,0.40033347120,0.44316542761,0.48758225561,0.53312891378,0.57927637250,0.62543608637,0.67098101782,0.71525516126,0.75760741222,0.79741168545,0.83410325108,0.86720651682,0.89637071504,0.92139174110,0.94223548231,0.95904493828,0.97211041164,0.98190281145,0.98909354112,0.99422007198,0.99758716905,0.99939322200,1.00000000000]
            @testset "precalced vectors" begin
                for time in 1:40
                    @test isapprox(decrement(jl,time), q_cumulative[time], atol = 1e-6)
                end

                for time in 1:40
                    q′ =  1 - survival(jl, time) / survival(jl, time-1)

                    @test isapprox(q′, q_annual[time], atol = 1e-6)
                end

                for time in 1:40
                    @test isapprox(survival(jl, time), 1 - q_cumulative[time], atol = 1e-6)
                end
            end
        end

    end
end


text:
using Pkg.Artifacts
using Pkg.BinaryPlatforms
using URIs

function tectonic()
    pkgname = "tectonic"
    origin = "https://github.com/tectonic-typesetting/tectonic/releases/download"
    version = v"0.8.0"
    build = 1

    downloads = Dict(
        "$origin/tectonic%40$version/tectonic-$version-x86_64-unknown-linux-musl.tar.gz" => Linux(:x86_64),
        "$origin/tectonic%40$version/tectonic-$version-x86_64-pc-windows-msvc.zip" => Windows(:x86_64),
        "$origin/tectonic%40$version/tectonic-$version-x86_64-apple-darwin.tar.gz" => MacOS(:x86_64),
    )

    build_path = joinpath(@__DIR__, "build")

    ispath(build_path) && rm(build_path; recursive=true, force=true)

    mkpath(build_path)

    artifact_toml = joinpath(@__DIR__, "Artifacts.toml")

    mktempdir() do temp_path
        for (url, platform) in downloads
            file = basename(url)
            @info "downloading" platform url file
            download(url, file)

            @info "unpacking" file
            product_hash = create_artifact() do artifact_dir
                if endswith(file, ".zip")
                    run(`unzip $file -d $artifact_dir`)
                elseif endswith(file, ".tar.gz")
                    run(`tar -xvf $file -C $artifact_dir`)
                else
                    error("cannot unpack that: $file")
                end
                files = readdir(artifact_dir)
            end
            @info "artifact" product_hash

            archive_filename = "$pkgname-$version+$build-$(triplet(platform)).tar.gz"
            download_hash = archive_artifact(product_hash, joinpath(build_path, archive_filename))

            @info "binding" archive_filename
            bind_artifact!(
                artifact_toml,
                "tectonic_bin",
                product_hash,
                platform=platform,
                force=true,
                download_info=Tuple[
                    (
                     "https://github.com/MichaelHatherly/Tectonic.jl/releases/download/tectonic-$(URIs.escapeuri("$(version)+$(build)"))/$archive_filename",
                     download_hash
                    )
                ]
            )
        end
    end
end

function biber()
    pkgname = "biber"
    origin = "https://sourceforge.net/projects/biblatex-biber/files/biblatex-biber"
    version = "2.14"
    build = 1

    downloads = Dict(
        "$origin/$version/binaries/Linux/biber-linux_x86_64.tar.gz" => Linux(:x86_64),
        "$origin/$version/binaries/Windows/biber-MSWIN64.zip" => Windows(:x86_64),
        "$origin/$version/binaries/OSX_Intel/biber-darwin_x86_64.tar.gz" => MacOS(:x86_64),
    )

    build_path = joinpath(@__DIR__, "build")

    ispath(build_path) && rm(build_path; recursive=true, force=true)

    mkpath(build_path)

    artifact_toml = joinpath(@__DIR__, "Artifacts.toml")

    mktempdir() do temp_path
        for (url, platform) in downloads
            file = basename(url)
            @info "downloading" platform url file
            download(url, file)

            @info "unpacking" file
            product_hash = create_artifact() do artifact_dir
                if endswith(file, ".zip")
                    run(`unzip $file -d $artifact_dir`)
                elseif endswith(file, ".tar.gz")
                    run(`tar -xvf $file -C $artifact_dir`)
                else
                    error("cannot unpack that: $file")
                end
                files = readdir(artifact_dir)
            end
            @info "artifact" product_hash

            archive_filename = "$pkgname-$version+$build-$(triplet(platform)).tar.gz"
            download_hash = archive_artifact(product_hash, joinpath(build_path, archive_filename))

            @info "binding" archive_filename
            bind_artifact!(
                artifact_toml,
                "biber_bin",
                product_hash,
                platform=platform,
                force=true,
                download_info=Tuple[
                    (
                     "https://github.com/MichaelHatherly/Tectonic.jl/releases/download/$(URIParser.escape("$(version)+$(build)"))/$archive_filename",
                     download_hash
                    )
                ]
            )
        end
    end
end


text:
"""

We assume that the initial data y0 is given at tspan[1], and that
tspan[end] is the last integration time.

"""
function ode{T<:Number}(F, y0, tspan::AbstractVector{T}, stepper::AbstractStepper;
                        jacobian::Function  = (t,y)->fdjacobian(F, t, y),
                        # we need these options explicitly for the dtinit
                        reltol::T   = eps(T)^T(1//3)/10,
                        abstol::T   = eps(T)^T(1//2)/10,
                        initstep::T = dtinit(F, y0, tspan, reltol, abstol; order=order(stepper)),
                        kargs...)

    # TODO: any ideas on how we could improve the interface so that we
    # don't have to use the ugly call to dtinit as a default?

    t0 = tspan[1]

    # construct a solver
    ode  = explicit_ineff(t0,y0,F,jac=jacobian)

    opts = Options{T}(;
                      tspan    = tspan,
                      reltol   = reltol,
                      abstol   = abstol,
                      initstep = initstep,
                      kargs...)
    solver = solve(ode,stepper,opts)

    # handle different directions of time integration
    if issorted(tspan)
        # do nothing, we are already set with the solver
        solution = collect(dense(solver))
    elseif issorted(reverse(tspan))
        # Reverse the time direction if necessary.  dense() only works
        # for positive time direction.

        # TODO: still ugly but slightly less bandaid-like then the
        # previous solution
        solution = map(ty->(2*t0-ty[1],ty[2]),collect(dense(reverse_time(solver))))
    else
        warn("Unsorted output times are not supported")
        return ([t0],[y0])
    end

    n = length(solution)

    # convert a list of pairs to a pair of arrays
    # TODO: leave it out as a list of pairs?
    tn = Array(T,n)
    yn = Array(typeof(y0),n)

    for (n,(t,y)) in enumerate(solution)
        tn[n] = t
        yn[n] = isa(y0,Number) ? y[1] : y
    end

    return (tn,yn)
end

ode23s{T}(F,y0,t0::Vector{T};kargs...)        = ode(F,y0,t0,ModifiedRosenbrockStepper{T}(); kargs...)
ode1{T}(F,y0,t0::Vector{T};kargs...)          = ode(F,y0,t0,RKStepperFixed{T}(bt_feuler); kargs...)
ode2_midpoint{T}(F,y0,t0::Vector{T};kargs...) = ode(F,y0,t0,RKStepperFixed{T}(bt_midpoint); kargs...)
ode2_heun{T}(F,y0,t0::Vector{T};kargs...)     = ode(F,y0,t0,RKStepperFixed{T}(bt_heun); kargs...)
ode4{T}(F,y0,t0::Vector{T};kargs...)          = ode(F,y0,t0,RKStepperFixed{T}(bt_rk4); kargs...)
ode21{T}(F,y0,t0::Vector{T};kargs...)         = ode(F,y0,t0,RKStepperAdaptive{T}(bt_rk21); kargs...)
ode23{T}(F,y0,t0::Vector{T};kargs...)         = ode(F,y0,t0,RKStepperAdaptive{T}(bt_rk23); kargs...)
ode45_fe{T}(F,y0,t0::Vector{T};kargs...)      = ode(F,y0,t0,RKStepperAdaptive{T}(bt_rk45); kargs...)
ode45_dp{T}(F,y0,t0::Vector{T};kargs...)      = ode(F,y0,t0,RKStepperAdaptive{T}(bt_dopri5); kargs...)
ode78{T}(F,y0,t0::Vector{T};kargs...)         = ode(F,y0,t0,RKStepperAdaptive{T}(bt_feh78); kargs...)

const ode45 = ode45_dp

"""

A nasty hack to convert a solver with negative time direction into a
solver with positive time direction.  This is necessary as negative
time direction is not supported by steppers (including the dense
output).  This only works for ExplicitODE.

"""
function reverse_time(sol::Solver)
    ode, options, stepper = sol.ode, sol.options, sol.stepper

    t0 = ode.t0
    y0 = ode.y0

    # TODO: improve the implementation
    function F_reverse!(t,y,dy)
        ode.F!(2*t0-t,y,dy)
        dy[:]=-dy
    end

    # TODO: is that how the jacobian changes?
    function jac_reverse!(t,y,J)
        ode.jac!(2*t0-t,y,J)
        J[:]=-J
    end

    # ExplicitODE is immutable
    ode_reversed = ExplicitODE(t0,y0,F_reverse!,jac_reverse!)
    stopevent = options.stopevent

    # TODO: we are modifying options here, should we construct new
    # options insted?
    options.tstop     = 2*t0-options.tstop
    options.tspan     = reverse(2*t0.-options.tspan)
    options.stopevent = (t,y)->stopevent(2*t0-t,y)
    return solve(ode_reversed,stepper,options)
end


text:
# TODO:
# (1) Is the mixed attribute handled correctly? Should it propogate 'into' groups, sequence etc.?

function transform_xsd_ast(ast::ASTNode, startelement)

  elementglobaldefs = Dict{AbstractString,ASTNode}()
  substitutiongroups = Dict{AbstractString,Array{AbstractString}}()
  substitutiongroupdefs = Dict{AbstractString,ASTNode}()
  attributeglobaldefs = Dict{AbstractString,ASTNode}()
  groupglobaldefs = Dict{AbstractString,ASTNode}()
  attributegroupglobaldefs = Dict{AbstractString,ASTNode}()
  typeglobaldefs = Dict{AbstractString,ASTNode}()

  process_xsd_element_nodes(ast, elementglobaldefs, substitutiongroups)
  process_xsd_attribute_nodes(ast, attributeglobaldefs)
  process_xsd_group_nodes(ast, groupglobaldefs)
  process_xsd_attributegroup_nodes(ast, attributegroupglobaldefs)
  process_xsd_simpletype_nodes(ast, typeglobaldefs)
  process_xsd_complextype_nodes(ast, typeglobaldefs)

  process_xsd_mixed_arg(ast)

  process_xsd_occurs_args(ast)
  process_xsd_use_arg(ast)

  process_xsd_compositor_nodes(ast)
  process_xsd_extension_nodes(ast)
  process_xsd_restriction_nodes(ast)
  process_xsd_union_nodes(ast)
  process_xsd_list_nodes(ast)

  process_xsd_any_nodes(ast)
  process_xsd_anyattribute_nodes(ast)

  process_xsd_substitution_groups(ast, substitutiongroups, substitutiongroupdefs)
  process_xsd_element_refs(ast, elementglobaldefs, substitutiongroupdefs)
  process_xsd_attribute_refs(ast, attributeglobaldefs)
  process_xsd_group_refs(ast, groupglobaldefs)
  process_xsd_attributegroup_refs(ast, attributegroupglobaldefs)
  process_xsd_type_refs(ast, ast, typeglobaldefs)

  if haskey(elementglobaldefs,startelement)
    push!(ast.refs, elementglobaldefs[startelement])
  else
    error("start element $(startelement) is not globally defined")
  end

end

# TODO:
#   fixed element/attribute content

# TODO - currently LightXML doesn't allow us to process namespace definitions, so
# we make the assumption that xs: is the prefix for XML Schema, and any other prefix
# is that of the target namespace

process_xsd_is_builtin_type(typename) = startswith(typename,"xs:") # TODO parameterise xs namespace

# assume any prefix other than xs: is that of the target namespace, remove it
function process_xsd_handle_namespace(fullname)
  shortname = fullname
  if !process_xsd_is_builtin_type(fullname)
    splitnames = split(fullname,":")
    shortname = splitnames[end]
  end
  shortname
end

function process_xsd_element_nodes(parentnode::ASTNode, elementglobaldefs, substitutiongroups)

  for node in parentnode.children

    if node.func == :element

      if haskey(node.args, :ref)

        node.func = :elementRef

      else

        if parentnode.func == :xsd

          elementname = node.args[:name]
          elementglobaldefs[elementname] = node

          if !haskey(substitutiongroups, elementname)
            substitutiongroups[elementname] = AbstractString[]
          end
          if !(haskey(node.args, :abstract) && (node.args[:abstract]=="true"))
            push!(substitutiongroups[elementname], elementname)
          end

          if haskey(node.args, :substitutionGroup)

            substitutiongroupname = process_xsd_handle_namespace(node.args[:substitutionGroup])

            if !haskey(substitutiongroups, substitutiongroupname )
              substitutiongroups[substitutiongroupname] = AbstractString[]
            end
            push!(substitutiongroups[substitutiongroupname], elementname)

            if substitutiongroupname != elementname
              # if no type specification, then must be the same type as the head of the substitution group
              if !haskey(node.args, :type) && isempty(node.children)
                elementsubnode = ASTNode(:elementSub)
                elementsubnode.args[:ref] = substitutiongroupname
                push!(node.children, elementsubnode)
              end
            end

          end

        end

        if haskey(node.args, :type)
          typerefnode = ASTNode(:typeRef)
          typerefnode.args[:ref] = node.args[:type]
          insert!(node.children, 1, typerefnode)
        end

        if haskey(node.args, :fixed)
          warn(node, "currently do not handle fixed attribute")
        end

      end

    end

    process_xsd_element_nodes(node, elementglobaldefs, substitutiongroups)

  end

end

function process_xsd_attribute_nodes(parentnode::ASTNode, attributeglobaldefs)

  for node in parentnode.children

    if node.func == :attribute

      if haskey(node.args, :ref)

        node.func = :attributeRef

      else

        if parentnode.func == :xsd
          attributeglobaldefs[node.args[:name]] = node
        end

        if haskey(node.args, :type)
          typerefnode = ASTNode(:typeRef)
          typerefnode.args[:ref] = node.args[:type]
          insert!(node.children, 1, typerefnode)
        end

        if isempty(node.children)
          warn(node, "No type found for attribute, so assuming xs:string")
          typerefnode = ASTNode(:typeRef)
          typerefnode.args[:ref] = "xs:string"
          insert!(node.children, 1, typerefnode)
        end

        if haskey(node.args, :fixed)
          warn(node, "currently do not handle fixed attribute")
        end

      end

    end

    process_xsd_attribute_nodes(node, attributeglobaldefs)

  end

end

function process_xsd_group_nodes(parentnode::ASTNode, groupglobaldefs)

  for node in parentnode.children

    if node.func == :group

      if parentnode.func == :xsd

        node.func = :group
        # TODO could check for name existence and uniqueness - but currently assume valid XSD
        groupglobaldefs[node.args[:name]] = node

      else

        node.func = :groupRef

      end

    end

    process_xsd_group_nodes(node, groupglobaldefs)

  end

end

function process_xsd_attributegroup_nodes(parentnode::ASTNode, attributegroupglobaldefs)

  for node in parentnode.children

    if node.func == :attributeGroup

      if parentnode.func == :xsd

        node.func = :attributeGroup
        # TODO could check for name existence and uniqueness - but currently assume valid XSD
        attributegroupglobaldefs[node.args[:name]] = node

      else

        node.func = :attributeGroupRef

      end

    end

    process_xsd_attributegroup_nodes(node, attributegroupglobaldefs)

  end

end

function process_xsd_simpletype_nodes(parentnode::ASTNode, typeglobaldefs)

  for node in parentnode.children

    if node.func == :simpleType

      if parentnode.func == :xsd

        # TODO could check for name existence and uniqueness - but currently assume valid XSD
        typeglobaldefs[node.args[:name]] = node

      end

    end

    process_xsd_simpletype_nodes(node, typeglobaldefs)

  end

end

function process_xsd_complextype_nodes(parentnode::ASTNode, typeglobaldefs)

  for node in parentnode.children

    if node.func == :complexType

      if parentnode.func == :xsd

        typeglobaldefs[node.args[:name]] = node

      end

    end

    process_xsd_complextype_nodes(node, typeglobaldefs)

  end

end

function process_xsd_compositor_nodes(parentnode::ASTNode)

  # TODO handling of choice and group - does choice apply within a referenced group (I think not)?

  for node in parentnode.children

    if node.func in [:sequence, :choice, :all]

      # TODO: not all/choice/sequence when within group (but defaulting of occurs is sufficient)
      # TODO: any

      if node.func==:all

        warn(node, "treating all node as a sequence")
        node.func = :sequence

      end

    end

    process_xsd_compositor_nodes(node)

  end

end

function process_xsd_substitution_groups(ast::ASTNode, substitutiongroups, substitutiongroupdefs)

  for substitutiongroup in substitutiongroups

    substitutiongroupname = substitutiongroup[1]
    elementnames = substitutiongroup[2]

    if !((length(elementnames)==1) && (substitutiongroupname==elementnames[1]))

      substitutiongroupnode = ASTNode(:substitutionGroup)
      substitutiongroupnode.args[:name] = substitutiongroupname

      for elementname in elementnames

        elementrefnode = ASTNode(:elementRef)
        elementrefnode.args[:ref] = elementname
        if elementname==substitutiongroup[1]
          elementrefnode.args[:allowsubstitutiongroup] = false
        end
        push!(substitutiongroupnode.children, elementrefnode)

      end

      if isempty(elementnames)
        warn(node, "substitution group $(substitutiongroupname) has no concrete elements")
      end

      push!(ast.children, substitutiongroupnode)

      substitutiongroupdefs[substitutiongroupname] = substitutiongroupnode

    end

  end

end

function process_xsd_element_refs(parentnode::ASTNode, elementglobaldefs, substitutiongroupdefs)

  for node in parentnode.children

    if node.func in [:elementRef, :elementSub]

      ref = node.args[:ref]
      ref = process_xsd_handle_namespace(ref)

      if (node.func==:elementRef) && haskey(substitutiongroupdefs, ref) && !(haskey(node.args, :allowsubstitutiongroup) && (node.args[:allowsubstitutiongroup]==false))
        push!(node.refs, substitutiongroupdefs[ref])
      elseif haskey(elementglobaldefs, ref)
        push!(node.refs, elementglobaldefs[ref])
      else
        error(node, "reference to an element, $(ref), that is not defined")
      end

    end

    process_xsd_element_refs(node, elementglobaldefs, substitutiongroupdefs)

  end

end

function process_xsd_attribute_refs(parentnode::ASTNode, attributeglobaldefs)

  for node in parentnode.children

    if node.func == :attributeRef

        ref = node.args[:ref]
        ref = process_xsd_handle_namespace(ref)
        if haskey(attributeglobaldefs, ref)
          push!(node.refs, attributeglobaldefs[ref])
        else
          error(node, "reference to an attribute, $(ref), that is not defined")
        end

    end

    process_xsd_attribute_refs(node, attributeglobaldefs)

  end

end

function process_xsd_group_refs(parentnode::ASTNode, groupglobaldefs)

  for node in parentnode.children

    if node.func == :groupRef

        ref = node.args[:ref]
        ref = process_xsd_handle_namespace(ref)
        if haskey(groupglobaldefs, ref)
          push!(node.refs, groupglobaldefs[ref])
        else
          error(node, "reference to a group, $(ref), that is not defined")
        end

    end

    process_xsd_group_refs(node, groupglobaldefs)

  end

end

function process_xsd_attributegroup_refs(parentnode::ASTNode, attributegroupglobaldefs)

  for node in parentnode.children

    if node.func == :attributeGroupRef

        ref = node.args[:ref]
        ref = process_xsd_handle_namespace(ref)
        if haskey(attributegroupglobaldefs, ref)
          push!(node.refs, attributegroupglobaldefs[ref])
        else
          error(node, "reference to an attribute group, $(ref), that is not defined")
        end

    end

    process_xsd_attributegroup_refs(node, attributegroupglobaldefs)

  end

end

function process_xsd_type_refs(ast::ASTNode, parentnode::ASTNode, typeglobaldefs)

  # println("Processing node $(parentnode.func) $(haskey(parentnode.args,:name) ? parentnode.args[:name] : '-')")

  for node in parentnode.children

    typeref = ""
    if node.func in [:typeRef]
      typeref = node.args[:ref]
    elseif node.func in [:list]
      if haskey(node.args, :itemType)
        typeref = node.args[:itemType]
      end
    end

    if !isempty(typeref)

      typeref = process_xsd_handle_namespace(typeref)

      if haskey(typeglobaldefs, typeref)

          push!(node.refs, typeglobaldefs[typeref])

      elseif process_xsd_is_builtin_type(typeref)

          push!(node.refs, process_xsd_define_builtin_type(ast, typeglobaldefs, typeref))

      else

          error(node, "A $(node.func) node references a type, $(typeref), that is not defined")

      end

    end

    process_xsd_type_refs(ast, node, typeglobaldefs)

  end

end

function process_xsd_extension_nodes(parentnode::ASTNode)

  for node in parentnode.children

    if node.func==:extension

      if (parentnode.func==:complexContent) && haskey(parentnode.args, :mixed)
        warn(node, "In complex content derived by extension, the mixed attribute of the complex content or complex type element may not be honoured")
      end

      if haskey(node.args, :base)
        typerefnode = ASTNode(:typeRef)
        typerefnode.args[:ref] = node.args[:base]
        insert!(node.children, 1, typerefnode)
      end

    end

    process_xsd_extension_nodes(node)

  end

end

function process_xsd_restriction_nodes(parentnode::ASTNode)

  for node in parentnode.children

    if node.func==:restriction

      if parentnode.func == :complexContent

        warn(node, "Cannot process restriction to produce complex content: base type will be used unrestricted, and the mixed attribute is ignored")

        typerefnode = ASTNode(:typeRef)
        typerefnode.args[:ref] = node.args[:base]

        for child in node.children
          append!(node.source, child.source)
        end

        node.children = [typerefnode]

      else

        facetfuncs = [:enumeration, :pattern, :fractionDigits, :length, :maxExclusive, :maxInclusive, :maxLength, :minExclusive, :minInclusive, :minLength, :totalDigits, :whiteSpace]
        typefuncs = [:simpleType]

        facetchildren = filter(child->child.func in facetfuncs, node.children)
        typechildren = filter(child->child.func in typefuncs, node.children)
        otherchildren = filter(child->!(child.func in facetfuncs || child.func in typefuncs), node.children)

        facets = map(child->child.func, facetchildren)

        if isempty(facets)

          if haskey(node.args, :base)

            typerefnode = ASTNode(:typeRef)
            typerefnode.args[:ref] = node.args[:base]

            node.children = [typerefnode; otherchildren]

          end

        elseif all(facet->facet in [:enumeration, :pattern], facets)

          choicenode = ASTNode(:choice)
          choicenode.children = facetchildren

          for typechild in typechildren
            append!(node.source, typechild.source)
          end

          node.children = [choicenode; otherchildren]

          for facetchild in facetchildren
            if facetchild.func == :pattern
              if !haskey(node.args, :base) || (node.args[:base] != "xs:string")
                warn(facetchild, "restriction pattern applied to type other than xs:string - will ignore base type")
              end
            end
          end

        else

          if !isempty(typechildren)

            warn(node, "Cannot process restriction of local defined type with facets $(join(facets,",")) - locally defined type will be used unrestricted")

            for facetchild in facetchildren
              append!(node.source, facetchild.source)
            end

            node.children = [typechildren; otherchildren]

          else

            @assert haskey(node.args, :base)

            if process_xsd_is_builtin_type(node.args[:base])

              pattern = process_xsd_derive_builtin_pattern(node.args[:base], node, facetchildren)

              patternnode = ASTNode(:pattern)
              patternnode.args[:value] = pattern
              patternnode.children = (ASTNode)[]

              for facetchild in facetchildren
                append!(node.source, facetchild.source)
              end

              node.children = [patternnode; otherchildren]

            else

              warn(node, "Cannot process restriction of type $(node.args[:base]) with facets $(join(facets,",")) - base type will be used unrestricted")

              typerefnode = ASTNode(:typeRef)
              typerefnode.args[:ref] = node.args[:base]

              for facetchild in facetchildren
                append!(node.source, facetchild.source)
              end

              node.children = [typerefnode; otherchildren]

            end

          end

        end

      end

    end

    process_xsd_restriction_nodes(node)

  end

end

function process_xsd_union_nodes(parentnode::ASTNode)

  for node in parentnode.children

    if node.func==:union

      if haskey(node.args, :memberTypes)

        membertypes = split(node.args[:memberTypes], " ")
        for membertype in membertypes

          typerefnode = ASTNode(:typeRef)
          typerefnode.args[:ref] = membertype

          push!(node.children, typerefnode)

        end

      end

    end

    process_xsd_union_nodes(node)

  end

end

function process_xsd_list_nodes(parentnode::ASTNode)

  for node in parentnode.children

    if node.func==:list

      if haskey(node.args, :itemType)
        typerefnode = ASTNode(:typeRef)
        typerefnode.args[:ref] = node.args[:itemType]
        push!(node.children, typerefnode)
      end

    end

    process_xsd_list_nodes(node)

  end

end

function process_xsd_any_nodes(parentnode::ASTNode)

  for node in parentnode.children

    if node.func==:any

      warn(node, "Cannot currently process any node - nothing will be output")

    end

    process_xsd_any_nodes(node)

  end

end

function process_xsd_anyattribute_nodes(parentnode::ASTNode)

  for node in parentnode.children

    if node.func==:anyAttribute

      warn(node, "Cannot currently process anyattribute node - nothing will be output")

    end

    process_xsd_anyattribute_nodes(node)

  end

end

function process_xsd_derive_builtin_pattern(typename, node::ASTNode, facetnodes=(ASTNode)[])

  pattern = ""

	# TODO: *could* add whitespace around all datatypes since this *should* be stripped off by processing

	# For the string-derived classes, we try to use character classes where possible to permit easy transition to Unicode tokenspaces from the current ASCII
	# Note:
	# \c is XSD extension to regex character classes - name characters: letters, digits, plus hyphen, period, colon
	# \i is XSD extension to regex character classes - initial name characters: letters, plus hyphen

  if isempty(facetnodes)
    if typename == "xs:string"
      pattern = ""		# interpreted by GT as a string of any valid characters, including CR/LF
    elseif typename == "xs:normalizedString"
      pattern = "" 		# serialisation space is the same as xs:string - the difference is the whitespace processing by the schema processor, not in the document itself
    elseif typename == "xs:token"
      pattern = "" 		# serialisation space is the same as xs:string - the difference is the whitespace processing by the schema processor, not in the document itself
    elseif typename == "xs:NMTOKEN"
      pattern = "\\c+"
    elseif typename == "xs:NMTOKENS"
      pattern = "\\c+(, \\c+)*"
    # elseif typename == "xs:Name"
    #   warn(node, "check pattern for built-in type $(typename)")
    #   pattern = "([A-Z]|[a-z]|[-:])+([A-Z]|[a-z]|[0-9]|[-:.])*" # TODO what other punctuation is allowed?
    elseif typename == "xs:NCName"
      pattern = "\\i\\c*"
    elseif typename == "xs:ID"
      warn(node, "uniqueness of xs:ID not enforced")
      pattern = "\\i\\c*"
    elseif typename == "xs:IDREF"
      warn(node, "existence of xs:IDREF not enforced")
      pattern = "\\i\\c*"
    elseif typename == "xs:decimal"
      pattern = "[+-]?([0-9]+|[0-9]*\\.[0-9]+)" # TODO permits -0
    elseif typename == "xs:integer"
      pattern = "[+-]?[0-9]+"  # TODO permits -0
    elseif typename == "xs:positiveInteger"
      pattern = "\\+?0+[1-9][0-9]+"
    elseif typename == "xs:negativeInteger"
      pattern = "-?0+[1-9][0-9]+"
    elseif typename == "xs:nonPositiveInteger"
      pattern = "0|-0+[1-9][0-9]+"
    elseif typename == "xs:nonNegativeInteger"
      pattern = "\\+?[0-9]+"
    elseif typename == "xs:float"
      pattern = "[+-]?([0-9]+|[0-9]*\\.[0-9]+)(E[+-]?[0-9]+)?|INF|-INF|NaN" # TODO not restricted to 32-bit float? permits E-0
    elseif typename == "xs:double"
      pattern = "[+-]?([0-9]+|[0-9]*\\.[0-9]+)(E[+-]?[0-9]+)?|INF|-INF|NaN" # TODO not restricted to 64-bit float? permits E-0
    elseif typename == "xs:boolean"
      pattern = "true|false"
    elseif typename in ["xs:date","xs:time","xs:dateTime"]
      datepattern = "[+-]?([0-9]{4}-((01|03|05|07|08|10|12)-(0[1-9]|[1-2][0-9]|3[0-1])|(04|06|09|11)-(0[1-9]|[1-2][0-9]|30)|02-(0[1-9]|1[0-9]|2[0-8]))|[0-9]{2}((1|3|5|7|9)(2|6)|(2|4|6|8)(0|4|8)|0(4|8))-02-29|((0|2|4|8)(0|4|8)|(1|3|5|7)(2|6))00-02-29)"
			timezonepattern = "(Z|[+-](0[0-9]|1[0-4]):[0-5][0-9])?"
			timepattern = "([0-1][0-9]|2[0-3]):[0-5][0-9]:[0-5][0-9](\\.[0-9]+)?"
			# TODO - this is very expensive in terms of choice points: consider using choose() from an appropriate datatype for this (and for numeric types) rather
			# than as a regex
			if typename == "xs:dateTime"
				pattern = datepattern * "T" * timepattern * timezonepattern
			elseif typename == "xs:date"
				pattern = datepattern * timezonepattern
			elseif typename == "xs:time"
				pattern = timepattern * timezonepattern
			end
		else
		  warn(node, "Do not know how to generate built-in type $(typename) - reverting to xs:string")
			pattern = process_xsd_derive_builtin_pattern("xs:string", node)
    end
	else
		# TODO - handle some facets
    warn(node, "Cannot process restriction of type $(typename) with facets $(join(map(child->child.func, facetchildren),",")) - base type will be used unrestricted")
		pattern = process_xsd_derive_builtin_pattern(typename, node)
  end

  pattern

end

function process_xsd_define_builtin_type(ast::ASTNode, typeglobaldefs, typename)

  simpletypenode = ASTNode(:simpleType)
  simpletypenode.args[:name] = typename

  pattern = process_xsd_derive_builtin_pattern(typename, simpletypenode)

  patternnode = ASTNode(:pattern)
  patternnode.args[:value] = pattern
  push!(simpletypenode.children, patternnode)

  push!(ast.children, simpletypenode)

  typeglobaldefs[typename] = simpletypenode

  simpletypenode

end

function process_xsd_occurs_args(parentnode::ASTNode)

  newchildren = ASTNode[]

  for node in parentnode.children

    newchild = node

    if (node.func in [:element, :elementRef, :all, :choice, :sequence, :typeRef, :any, :group]) && (parentnode.func != :xsd)

      lowerbound = 1
      if haskey(node.args, :minOccurs)
        lowerbound = parse(Int, node.args[:minOccurs])
      end

      upperbound = 1
      if haskey(node.args, :maxOccurs)
        if node.args[:maxOccurs] == "unbounded"
          upperbound = Inf
        else
          upperbound = parse(Int, node.args[:maxOccurs])
        end
      end

      if !((lowerbound==1) && (upperbound==1))
        quantifiernode = ASTNode(:quantifier)
        quantifiernode.children = [node]
        quantifiernode.args[:lowerbound] = lowerbound
        quantifiernode.args[:upperbound] = upperbound
        quantifiernode.source = node.source
        newchild = quantifiernode
      end

    end

    push!(newchildren, newchild)

    process_xsd_occurs_args(node)

  end

  parentnode.children = newchildren

end

function process_xsd_use_arg(parentnode::ASTNode)

  newchildren = ASTNode[]

  for node in parentnode.children

    newchild = node

    if (node.func in [:attribute, :attributeRef]) && (parentnode.func != :xsd)

			optional = true
      if haskey(node.args, :use)
        if node.args[:use] == "prohibited"
          node.func = :nothing
					optional = false
        elseif node.args[:use] == "required"
          optional = false
        end
      end

      if optional
        optionalnode = ASTNode(:optional)
        optionalnode.children = [node]
        optionalnode.source = node.source
        newchild = optionalnode
      end

    end

    push!(newchildren, newchild)

    process_xsd_use_arg(node)

  end

  parentnode.children = newchildren

end

function process_xsd_mixed_arg(parentnode::ASTNode)
  for node in parentnode.children
		if node.func in [:complexType]
			add_optional_text_nodes(node)
		end
    process_xsd_mixed_arg(node)
  end
end

function add_optional_text_nodes(node::ASTNode, parentmixedcontent=false, parentafterfirstelement=false)

	mixedcontent = parentmixedcontent
	afterfirstelement = parentafterfirstelement

	if (node.func in [:complexContent, :complexType]) && haskey(node.args, :mixed)
		mixedcontent = (node.args[:mixed] == "true")
	end

	# process children first so that we process elements (and groups) in physical order so that
	# we correctly detect the 'first' element in this mixed contentType subtree
  for child in node.children
    if child.func in [:sequence, :choice, :all, :complexType, :restriction, :extension]
		  afterfirstelement = add_optional_text_nodes(child, mixedcontent, afterfirstelement)
		end
	end

	if mixedcontent

    mixedchildren = (ASTNode)[]

    for child in node.children

      if child.func in [:element, :elementRef, :groupRef]
				# if within this mixed contentType subtree we have not yet reached first element
				# then allow text before this first element as well as after
				if !afterfirstelement
          push!(mixedchildren, create_optional_text_node())
					afterfirstelement = true
				end
        push!(mixedchildren, child)
        push!(mixedchildren, create_optional_text_node())
			else
        push!(mixedchildren, child)
      end

    end

    node.children = mixedchildren

  end

	afterfirstelement

end

function create_optional_text_node()
	textnode = ASTNode(:typeRef)
	textnode.args[:ref] = "xs:string"
	optionalnode = ASTNode(:optional)
	optionalnode.children = [textnode]
	optionalnode
end


text:
const _ARFF_SPACE       = UInt8(' ')
const _ARFF_COMMENT     = UInt8('%')
const _ARFF_AT          = UInt8('@')
const _ARFF_SEP         = UInt8(',')
const _ARFF_NEWLINE     = UInt8('\n')
const _ARFF_NOMSTART    = UInt8('{')
const _ARFF_NOMEND      = UInt8('}')
const _ARFF_ESC         = UInt8('\\')
const _ARFF_MISSING     = UInt8('?')
const _ARFF_RELMARK     = UInt8('\'')

function readARFF(p::String)
    open(p, "r") do io
        lines = readlines(io)
        for i in 1:length(lines)
            line = lines[i]
            # If not empty line or comment
            if !isempty(line) && UInt8(line[1]) != _ARFF_COMMENT
                sline = split(line, " ")
                # println(sline[1][1])
                # If the first symbol is @
                if UInt8(sline[1][1]) == _ARFF_AT
                    # If @relation
                    if sline[1][2:end] == "relation"
                        println(i)
                        println("Relation: " * sline[2])
                    end
                end
            end
            #error(1)
        end
    end
end

# for (i, x) in enumerate([1,2,3,4])
#     println(i, " ", x)
#     i=3
# end

# @grammar begin
#     S   = UnivGlobal(F) | ExistGlobal(F)
#     F   = UnivL(implies(F,F))
#     F   = ExistL(and(F,F))
#     F   = implies(F,F)
#     F   = and(F,F)
#     F   = f(A)
# end
#
# @grammar begin
#     S   = ExistentialIntervalRelation{:L}(and(F,F))
#     F   = f(A)
# end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module J2CArray

using Compat

export to_j2c_array, from_j2c_array, j2c_array_delete, from_ascii_string
import CompilerTools.Helper.isArrayType
import ..getPackageRoot

#eval(x) = Core.eval(J2CArray, x)

function getLib()
    package_root = getPackageRoot()
    if Compat.is_windows()
        return string(package_root, "/deps/libj2carray.dll")
    else
        return string(package_root, "/deps/libj2carray.so.1.0")
    end
end

if !isfile(getLib())
  error("libj2carray not compiled, build with: julia -e 'Pkg.build(\"ParallelAccelerator\")'")
end

function __init__()
  dyn_lib = getLib()

  @eval begin
    # Create a new j2c array object with element size in bytes and given dimension.
    # It will share the data pointer of the given inp array, and if inp is nothing,
    # the j2c array will allocate fresh memory to hold data.
    # NOTE: when elem_bytes is 0, it means the elements must be j2c array type
#    function j2c_array_new(elem_bytes::Int, inp::Union{Array, Void}, ndim::Int, dims::Tuple)
#      # note that C interface mandates Int64 for dimension data
#      _dims = Int64[ convert(Int64, x) for x in dims ]
#      _inp = is(inp, nothing) ? C_NULL : convert(Ptr{Void}, pointer(inp))
#
#      #ccall((:j2c_array_new, $dyn_lib), Ptr{Void}, (Cint, Ptr{Void}, Cuint, Ptr{UInt64}),
#      #      convert(Cint, elem_bytes), _inp, convert(Cuint, ndim), pointer(_dims))
#      ccall((:j2c_array_new, cgen_lib), Ptr{Void}, (Cint, Ptr{Void}, Cuint, Ptr{UInt64}),
#            convert(Cint, elem_bytes), _inp, convert(Cuint, ndim), pointer(_dims))
#    end

    # Array size in the given dimension.
    function j2c_array_size(arr::Ptr{Void}, dim::Int)
      l = ccall((:j2c_array_size, $dyn_lib), Cuint, (Ptr{Void}, Cuint),
                arr, convert(Cuint, dim))
      return convert(Int, l)
    end

    # Retrieve j2c array data pointer, and parameter "own" means caller will
    # handle the memory deallocation of this pointer.
    function j2c_array_to_pointer(arr::Ptr{Void}, own::Bool)
      ccall((:j2c_array_to_pointer,$dyn_lib), Ptr{Void}, (Ptr{Void}, Bool), arr, own)
    end

    # Read the j2c array element of given type at the given (linear) index.
    # If T is Ptr{Void}, treat the element type as j2c array, and the
    # returned array is merely a pointer, not a new object.
    function j2c_array_get(arr::Ptr{Void}, idx::Int, T::Type)
      nbytes = is(T, Ptr{Void}) ? 0 : sizeof(T)
      _value = Array(T, 1)
      ccall((:j2c_array_get,$dyn_lib), Void, (Cint, Ptr{Void}, Cuint, Ptr{Void}),
            convert(Cint, nbytes), arr, convert(Cuint, idx), convert(Ptr{Void}, pointer(_value)))
      return _value[1]
    end

    # Set the j2c array element at the given (linear) index to the given value.
    # If T is Ptr{Void}, treat value as a pointer to j2c array.
    function j2c_array_set{T}(arr::Ptr{Void}, idx::Int, value::T)
      nbytes = is(T, Ptr{Void}) ? 0 : sizeof(T)
      _value = nbytes == 0 ? value : convert(Ptr{Void}, pointer(T[ value ]))
      ccall((:j2c_array_set, $dyn_lib), Void, (Cint, Ptr{Void}, Cuint, Ptr{Void}),
            convert(Cint, nbytes), arr, convert(Cuint, idx), _value)
    end

    # Delete an j2c array object.
    # Note that this only works for scalar array or arrays of
    # objects whose derefrence will definite not trigger nested
    # deletion (either data pointer is NULL, or refcount > 1).
    # Currently there is no way to cleanly delete an nested j2c
    # array without first converting back to a julia array.
    function j2c_array_delete(arr::Ptr{Void})
      ccall((:j2c_array_delete,$dyn_lib), Void, (Ptr{Void},), arr)
    end

    # Dereference a j2c array data pointer.
    # Require that the j2c array data pointer points to either
    # a scalar array, or an array of already dereferenced array
    # (whose data pointers are NULL).
    function j2c_array_deref(arr::Ptr{Void})
      ccall((:j2c_array_deref,$dyn_lib), Void, (Ptr{Void},), arr)
    end

    # convert ASCIIString to J2C_Array and to C's ASCIIString
    function to_j2c_array(inp::AbstractString, ptr_array_dict, mapAtypeKey, j2c_array_new)
      arr = to_j2c_array(inp.data, ptr_array_dict, mapAtypeKey, j2c_array_new)
      ccall((:new_ascii_string, $dyn_lib), Ptr{Void}, (Ptr{Void}, ), arr)
    end

    # convert C's ASCIIString to Julia's Array{UInt8,1}
    function from_ascii_string(str::Ptr{Void}, ptr_array_dict)
      data = ccall((:from_ascii_string, $dyn_lib), Ptr{Void}, (Ptr{Void}, ), str)
      arr = _from_j2c_array(data, UInt8, 1, ptr_array_dict)
      ccall((:delete_ascii_string, $dyn_lib), Void, (Ptr{Void}, ), str)
      return arr
    end
  end
end

# Convert Julia array to J2C array object.
# Note that Julia array data are not copied but shared by the J2C array
# The caller needs to make sure these arrays stay alive so long as the
# returned j2c array is alive.
function to_j2c_array{T, N}(inp :: Array{T, N}, ptr_array_dict :: Dict{Ptr{Void}, Array}, mapAtypeKey :: Dict{DataType,Int64}, j2c_array_new)
  allocation_key = mapAtypeKey[typeof(inp)]
  dims = size(inp)
  _isbits = isbits(T)
  nbytes = _isbits ? sizeof(T) : 0
  _inp = _isbits ? inp : nothing
  #arr = j2c_array_new(nbytes, _inp, N, dims)
  arr = j2c_array_new(allocation_key, _inp, N, dims)
  ptr_array_dict[convert(Ptr{Void}, pointer(inp))] = inp  # establish a mapping between pointer and the original array
  if !(_isbits)
    for i = 1:length(inp)
      obj = to_j2c_array(inp[i], ptr_array_dict, mapAtypeKey, j2c_array_new) # obj is a new j2c array
      j2c_array_set(arr, i, obj) # obj is duplicated during this set
      j2c_array_delete(obj)      # safe to delete obj without triggering free
    end
  end
  return arr
end

# Convert J2C array object to Julia array.
# Note that:
# 1. We assume the input j2c array object contains no data pointer aliases.
# 2. The returned Julia array will share the pointer to J2C array data at leaf level.
function _from_j2c_array(inp::Ptr{Void}, elem_typ::DataType, N::Int, ptr_array_dict :: Dict{Ptr{Void}, Array})
  dims = Array(Int, N)
  len  = 1
  for i = 1:N
    dims[i] = j2c_array_size(inp, i)
    len = len * dims[i]
  end
  if isbits(elem_typ)
    array_ptr = convert(Ptr{Void}, j2c_array_to_pointer(inp, true))
    if haskey(ptr_array_dict, array_ptr)
      arr = ptr_array_dict[array_ptr]
    else
if VERSION > v"0.5.0-dev+3260"
      arr = unsafe_wrap(Array,convert(Ptr{elem_typ}, array_ptr), tuple(dims...), true)
else
      arr = pointer_to_array(convert(Ptr{elem_typ}, array_ptr), tuple(dims...), true)
end
    end
  elseif isArrayType(elem_typ)
    arr = Array(elem_typ, dims...)
    sub_type = elem_typ.parameters[1]
    sub_dim  = elem_typ.parameters[2]
    for i = 1:len
      ptr = j2c_array_get(inp, i, Ptr{Void})
      arr[i] = _from_j2c_array(ptr, sub_type, sub_dim, ptr_array_dict)
    end
  end
  return arr
end

function from_j2c_array(inp::Ptr{Void}, elem_typ::DataType, N::Int, ptr_array_dict :: Dict{Ptr{Void}, Array})
   arr = _from_j2c_array(inp, elem_typ, N, ptr_array_dict)
   j2c_array_delete(inp)
   return arr
end

end


text:
using DataFrames, CSV

ifn = "csv_data_manipulation_in.dat"
ofn = "csv_data_manipulation_out.dat"

df = CSV.read(ifn)
df.SUM = sum.(eachrow(df))
CSV.write(ofn, df)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# --------------- Base type for model:
module Models

using Flux, LinearAlgebra

include("base.jl")

export AbstractFittedModel, AbstractDifferentiableModel,
    FluxModel, LogisticModel, BayesianLogisticModel,
    RTorchModel, PyTorchModel,
    probs, logits

end


text:
@with_kw immutable DiscreteVDPTagMDP <: MDP{Int, Int}
    cmdp::VDPTagMDP     = VDPTagMDP()
    n_bins::Int         = 60
    grid_lim::Float64   = 3.0
    n_angles::Int       = 10
end

@with_kw immutable DiscreteVDPTagPOMDP <: POMDP{Int, Int, Int}
    cpomdp::VDPTagPOMDP = VDPTagPOMDP()
    n_bins::Int         = 60
    grid_lim::Float64   = 3.0
    n_angles::Int       = 10
    n_obs_angles::Int   = 10
end

@with_kw immutable AODiscreteVDPTagPOMDP <: POMDP{TagState, Int, Int}
    cpomdp::VDPTagPOMDP = VDPTagPOMDP()
    n_angles::Int       = 10
    n_obs_angles::Int   = 10
end

@with_kw immutable ADiscreteVDPTagPOMDP <: POMDP{TagState, Int, Float64}
    cpomdp::VDPTagPOMDP = VDPTagPOMDP()
    n_angles::Int       = 10
end

const DiscreteVDPTagProblem = Union{DiscreteVDPTagMDP, DiscreteVDPTagPOMDP, AODiscreteVDPTagPOMDP, ADiscreteVDPTagPOMDP}

mdp(p::DiscreteVDPTagMDP) = p
mdp(p::DiscreteVDPTagPOMDP) = DiscreteVDPTagMDP(p.cpomdp.mdp, p.n_bins, p.grid_lim, p.n_angles)
cproblem(p::DiscreteVDPTagMDP) = p.cmdp
cproblem(p::DiscreteVDPTagPOMDP) = p.cpomdp
cproblem(p::AODiscreteVDPTagPOMDP) = p.cpomdp
cproblem(p::ADiscreteVDPTagPOMDP) = p.cpomdp

convert_s{T}(::Type{T}, x::T, p) = x
convert_a{T}(::Type{T}, x::T, p) = x
convert_o{T}(::Type{T}, x::T, p) = x

# state
function convert_s(::Type{Int}, s::TagState, p::DiscreteVDPTagProblem)
    n = p.n_bins
    factor = n/(2*p.grid_lim)
    ai = clamp(ceil(Int, (s.agent[1]+p.grid_lim)*factor), 1, n)
    aj = clamp(ceil(Int, (s.agent[2]+p.grid_lim)*factor), 1, n)
    ti = clamp(ceil(Int, (s.target[1]+p.grid_lim)*factor), 1, n)
    tj = clamp(ceil(Int, (s.target[2]+p.grid_lim)*factor), 1, n)
    return sub2ind((n,n,n,n), ai, aj, ti, tj)
end
function convert_s(::Type{TagState}, s::Int, p::DiscreteVDPTagProblem)
    n = p.n_bins
    factor = 2*p.grid_lim/n
    ai, aj, ti, tj = ind2sub((n,n,n,n), s)
    return TagState((Vec2(ai, aj)-0.5)*factor-p.grid_lim, (Vec2(ti, tj)-0.5)*factor-p.grid_lim)
end

# action
function convert_a(::Type{Int}, a::Float64, p::DiscreteVDPTagProblem)
    i = ceil(Int, a*p.n_angles/(2*pi))
    while i > p.n_angles
        i -= p.n_angles
    end
    while i < 1
        i += p.n_angles
    end
    return i
end
convert_a(::Type{Float64}, a::Int, p::DiscreteVDPTagProblem) = (a-0.5)*2*pi/p.n_angles

function convert_a(T::Type{Int}, a::TagAction, p::DiscreteVDPTagProblem)
    i = convert_a(T, a.angle, p)
    if a.look
        return i + p.n_angles
    else
        return i
    end
end
function convert_a(::Type{TagAction}, a::Int, p::DiscreteVDPTagProblem)
    return TagAction(a > p.n_angles, convert_a(Float64, a % p.n_angles, p))
end

# observation
function convert_o(::Type{Int}, o::Float64, p::DiscreteVDPTagProblem)
    i = ceil(Int, o*p.n_obs_angles/(2*pi))
    while i > p.n_obs_angles
        i -= p.n_obs_angles
    end
    while i < 1
        i += p.n_obs_angles
    end
    return i
end
convert_o(::Type{Float64}, o::Int, p::DiscreteVDPTagProblem) = (o-0.5)*2*pi/p.n_obs_angles

n_states(p::DiscreteVDPTagProblem) = mdp(p).n_bins^4
n_states(p::AODiscreteVDPTagPOMDP) = Inf
n_actions(p::DiscreteVDPTagMDP) = p.n_angles
n_actions(p::DiscreteVDPTagProblem) = 2*p.n_angles
n_observations(p::DiscreteVDPTagProblem) = p.n_obs_angles
discount(p::DiscreteVDPTagProblem) = discount(cproblem(p))
isterminal(p::DiscreteVDPTagProblem, s) = isterminal(cproblem(p), convert_s(TagState, s, p))
observations(p::DiscreteVDPTagProblem) = 1:p.n_angles

function generate_s(p::DiscreteVDPTagProblem, s::Int, a::Int, rng::AbstractRNG)
    cs = convert_s(TagState, s, p)
    ca = convert_a(action_type(cproblem(p)), a, p)
    csp = generate_s(cproblem(p), cs, ca, rng)
    return convert_s(Int, csp, p)
end

function generate_sr(p::DiscreteVDPTagProblem, s::Int, a::Int, rng::AbstractRNG)
    cs = convert_s(TagState, s, p)
    ca = convert_a(action_type(cproblem(p)), a, p)
    csp = generate_s(cproblem(p), cs, ca, rng)
    r = reward(cproblem(p), cs, ca, csp)
    return (convert_s(Int, csp, p), r)
end

function generate_sor(p::DiscreteVDPTagPOMDP, s::Int, a::Int, rng::AbstractRNG)
    cs = convert_s(TagState, s, p)
    ca = convert_a(action_type(cproblem(p)), a, p)
    csor = generate_sor(cproblem(p), cs, ca, rng)
    return (convert_s(Int, csor[1], p), convert_o(Int, csor[2], p), csor[3])
end

actions(p::DiscreteVDPTagProblem) = 1:n_actions(p)

immutable DiscreteVDPInitDist
    p::DiscreteVDPTagProblem
end
sampletype(::Type{DiscreteVDPInitDist}) = Int
function rand(rng::AbstractRNG, d::DiscreteVDPInitDist)
    cs = rand(rng, VDPInitDist())
    return convert_s(Int, cs, d.p)
end
initial_state_distribution(p::DiscreteVDPTagProblem) = DiscreteVDPInitDist(p)

function generate_s(p::DiscreteVDPTagProblem, s::TagState, a::Int, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    return generate_s(cproblem(p), s, ca, rng)
end

function generate_sr(p::DiscreteVDPTagProblem, s::TagState, a::Int, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    sp = generate_s(cproblem(p), s, ca, rng)
    r = reward(cproblem(p), s, ca, sp)
    return (sp, r)
end

function generate_sor(p::ADiscreteVDPTagPOMDP, s::TagState, a::Int, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    return generate_sor(cproblem(p), s, ca, rng)
end

function generate_o(p::ADiscreteVDPTagPOMDP, s::TagState, a::Int, sp::TagState, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    return generate_o(cproblem(p), s, ca, sp, rng)
end

function generate_sor(p::AODiscreteVDPTagPOMDP, s::TagState, a::Int, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    csor = generate_sor(cproblem(p), s, ca, rng)
    return (csor[1], convert_o(Int, csor[2], p), csor[3])
end

function generate_o(p::AODiscreteVDPTagPOMDP, s::TagState, a::Int, sp::TagState, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    co = generate_o(cproblem(p), s, ca, sp, rng)
    return convert_o(Int, co, p)
end

initial_state_distribution(p::AODiscreteVDPTagPOMDP) = VDPInitDist()
initial_state_distribution(p::ADiscreteVDPTagPOMDP) = VDPInitDist()

gauss_cdf(mean, std, x) = 0.5*(1.0+erf((x-mean)/(std*sqrt(2))))

function obs_weight(p::AODiscreteVDPTagPOMDP, a::Int, sp::TagState, o::Int)
    cp = cproblem(p)
    @assert cp.bearing_std <= 2*pi/6.0 "obs_weight assumes σ <= $(2*pi/6.0)"
    ca = convert_a(action_type(cp), a, p)
    co = convert_o(obs_type(cp), o, p) # float between 0 and 2pi
    upper = co + 0.5*2*pi/p.n_angles
    lower = co - 0.5*2*pi/p.n_angles
    if ca.look
        diff = sp.target - sp.agent
        bearing = atan2(diff[2], diff[1])
        # three cases: o is in bin, below, or above
        if bearing <= upper && bearing > lower
            cdf_up = gauss_cdf(bearing, cp.bearing_std, upper)
            cdf_low = gauss_cdf(bearing, cp.bearing_std, lower)
            prob = cdf_up - cdf_low
        elseif bearing <= lower
            cdf_up = gauss_cdf(bearing, cp.bearing_std, upper)
            cdf_low = gauss_cdf(bearing, cp.bearing_std, lower)
            below_cdf_up = gauss_cdf(bearing, cp.bearing_std, upper-2*pi)
            below_cdf_low = gauss_cdf(bearing, cp.bearing_std, lower-2*pi)
            prob = cdf_up - cdf_low + below_cdf_up - below_cdf_low
        else # bearing > upper
            cdf_up = gauss_cdf(bearing, cp.bearing_std, upper)
            cdf_low = gauss_cdf(bearing, cp.bearing_std, lower)
            above_cdf_up = gauss_cdf(bearing, cp.bearing_std, upper+2*pi)
            above_cdf_low = gauss_cdf(bearing, cp.bearing_std, lower+2*pi)
            prob = cdf_up - cdf_low + above_cdf_up - above_cdf_low
        end
        return prob
    else
        return 1.0
    end
end

function obs_weight(p::ADiscreteVDPTagPOMDP, a::Int, sp::TagState, o::Float64)
    ca = convert_a(TagAction, a, p)
    return obs_weight(cproblem(p), ca, sp, o)
end


text:
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "output_file"
        help = "Output location to save the new data bundle."
        arg_type = ASCIIString
        required = true
    "input_config"
        help = "File listing which external data files to integrate into the network. These can be"*
               "either BAM files or BED files."*
               "For each line trailing fields may be omitted and follow the tab separated format: "*
               "FILE_NAME SHORT_TITLE LONG_TITLE CELL_TYPE LAB EXPERIMENT_ID ANTIBODY_ID TREATMENTS ORGANISM LIFE_STAGE LINK"
        arg_type = ASCIIString
        required = true
    "--assembly"
        help = "Which assembly genome coordinates are in (currently GRCh38, GRCh38_alt, hg19, mm10, mm9)."
        arg_type = ASCIIString
        default = "GRCh38_alt"
    "--quiet", "-q"
        help = "Don't print anything"
        action = :store_true
end
args = parse_args(s)

using ChromNet
using JSON
using GZip
using JLD
using SamIO
using ProgressMeter

# parse the input arguments
inputConfig = args["input_config"]
outputFile = args["output_file"]
if !ismatch(r"\.ChromNet\.jld$", outputFile)
    outputFile = "$outputFile.ChromNet.jld"
end
quiet = args["quiet"]

# load the metadata
inStream = STDIN
dataFileRoot = "."
if inputConfig != "-"
    inStream = open(inputConfig)
    dataFileRoot = dirname(inputConfig)
end
metadata = parse_config(inStream, dataFileRoot)
inputConfig != "-" && close(inStream)

# extract all the data files and custom dataset ids
fileTypes = ASCIIString[]
fileNames = ASCIIString[]
ids = ASCIIString[]
for (k,v) in metadata
    if v["fileType"] == "unknown"
        quiet || println(STDERR, "Ignoring file of unknown type: ", v["dataFile"])
        continue
    end
    push!(fileTypes, v["fileType"])
    push!(fileNames, v["dataFile"])
    push!(ids, k)
end

referenceContigs = SamIO.assembly[args["assembly"]]

# build the custom data matrix
numBins = ceil(Int64, sum(referenceContigs.sizes) / 1000)
data = zeros(Float32, length(ids), numBins)
@showprogress "Binning BAM/BED files: " for i in 1:length(ids)
    if fileTypes[i] == "bed" || fileTypes[i] == "narrowPeak"
        reader = BinningMap(BedReader(fileNames[i], referenceContigs), 1000)
        while !eof(reader)
            data[i,position(reader)] = value(reader)
            advance!(reader)
        end
        close(reader)

    elseif fileTypes[i] == "bed.gz" || fileTypes[i] == "narrowPeak.gz"
        f = GZip.open(fileNames[i])
        reader = BinningMap(BedReader(f, referenceContigs), 1000)
        while !eof(reader)
            data[i,position(reader)] = value(reader)
            advance!(reader)
        end
        close(reader)

    elseif fileTypes[i] == "bam"
        reader = BinningMap(BamReader(fileNames[i], :any, referenceContigs), 1000)
        while !eof(reader)
            data[i,position(reader)] = value(reader)
            advance!(reader)
        end
        close(reader)
    end
end

# build summary statistics
!quiet && println(STDERR, "Computing means and norms...")
means = vec(mean(data, 2))
norms = Float64[norm(data[i,:] - means[i]) for i in 1:size(data)[1]]

# write the JLD file
!quiet && println(STDERR, "Writing results to $outputFile...")
jldopen(outputFile, "w") do f
    f["data"] = data
    f["means"] = vec(mean(data, 2))
    f["norms"] = Float64[norm(data[i,:] - means[i]) for i in 1:size(data)[1]]
    f["ids"] = ids
    f["P"] = length(ids)
    f["N"] = size(data)[2]
    f["metadata"] = metadata
    f["assembly"] = args["assembly"]
end