text:
module hyperelasticity

using FenicsPy

mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

# Define Dirichlet boundary (x = 0 or x = 1)
c = Constant((0.0, 0.0, 0.0))
r = Expression(("scale*0.0",
                "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
                "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
                scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3, degree=2)

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = FeFunction(V)                 # Displacement from previous iteration
B  = Constant((0.0, -0.5, 0.0))  # Body force per unit volume
T  = Constant((0.1,  0.0, 0.0))  # Traction force on the boundary

# Kinematics
d = len(u)
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))^2

# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Solve variational problem
solve(F == 0, u, bcs, J=J)

# Save solution in VTK format
file = File("hyperelasticity/displacement.pvd")
file << u

# Plot solution
#plot(u)

end # module hyperelasticity


text:
"""
    UniformBlockDiagonal{T}

Homogeneous block diagonal matrices.  `k` diagonal blocks each of size `m×m`
"""
struct UniformBlockDiagonal{T} <: AbstractMatrix{T}
    data::Array{T,3}
    facevec::Vector{SubArray{T,2,Array{T,3}}}
end

function UniformBlockDiagonal(dat::Array{T,3}) where {T}
    UniformBlockDiagonal(
        dat,
        SubArray{T,2,Array{T,3}}[view(dat, :, :, i) for i = 1:size(dat, 3)],
    )
end

function Base.copyto!(dest::UniformBlockDiagonal{T}, src::UniformBlockDiagonal{T}) where {T}
    sdat = src.data
    ddat = dest.data
    size(ddat) == size(sdat) || throw(DimensionMismatch(""))
    copyto!(ddat, sdat)
    dest
end

function Base.copyto!(dest::Matrix{T}, src::UniformBlockDiagonal{T}) where {T}
    size(dest) == size(src) || throw(DimensionMismatch(""))
    fill!(dest, zero(T))
    sdat = src.data
    m, n, l = size(sdat)
    for k = 1:l
        ioffset = (k - 1) * m
        joffset = (k - 1) * n
        for j = 1:n
            jind = joffset + j
            for i = 1:m
                dest[ioffset+i, jind] = sdat[i, j, k]
            end
        end
    end
    dest
end

function Base.getindex(A::UniformBlockDiagonal{T}, i::Int, j::Int) where {T}
    Ad = A.data
    m, n, l = size(Ad)
    (0 < i ≤ l * m && 0 < j ≤ l * n) ||
    throw(IndexError("attempt to access $(l*m) × $(l*n) array at index [$i, $j]"))
    iblk, ioffset = divrem(i - 1, m)
    jblk, joffset = divrem(j - 1, n)
    iblk == jblk ? Ad[ioffset+1, joffset+1, iblk+1] : zero(T)
end

function LinearAlgebra.Matrix(A::UniformBlockDiagonal{T}) where {T}
    Ad = A.data
    m, n, l = size(Ad)
    mat = zeros(T, (m * l, n * l))
    @inbounds for k = 0:(l-1)
        kp1 = k + 1
        km = k * m
        kn = k * n
        for j = 1:n
            knpj = kn + j
            for i = 1:m
                mat[km+i, knpj] = Ad[i, j, kp1]
            end
        end
    end
    mat
end

function Base.size(A::UniformBlockDiagonal)
    m, n, l = size(A.data)
    (l * m, l * n)
end

"""
    BlockedSparse{Tv,S,P}

A `SparseMatrixCSC` whose nonzeros form blocks of rows or columns or both.

# Members
* `cscmat`: `SparseMatrixCSC{Tv, Int32}` representation for general calculations
* `nzasmat`: nonzeros of `cscmat` as a dense matrix
* `colblkptr`: pattern of blocks of columns

The only time these are created are as products of `ReMat`s.
"""
mutable struct BlockedSparse{T,S,P} <: AbstractMatrix{T}
    cscmat::SparseMatrixCSC{T,Int32}
    nzsasmat::Matrix{T}
    colblkptr::Vector{Int32}
end

Base.size(A::BlockedSparse) = size(A.cscmat)

Base.size(A::BlockedSparse, d) = size(A.cscmat, d)

Base.getindex(A::BlockedSparse, i::Integer, j::Integer) = getindex(A.cscmat, i, j)

LinearAlgebra.Matrix(A::BlockedSparse) = Matrix(A.cscmat)

SparseArrays.sparse(A::BlockedSparse) = A.cscmat

SparseArrays.nnz(A::BlockedSparse) = nnz(A.cscmat)

function Base.copyto!(L::BlockedSparse{T}, A::SparseMatrixCSC{T}) where {T}
    size(L) == size(A) && nnz(L) == nnz(A) ||
    throw(DimensionMismatch("size(L) ≠ size(A) or nnz(L) ≠ nnz(A"))
    copyto!(nonzeros(L.cscmat), nonzeros(A))
    L
end

ext: jl


text:
# TODO: combine into one function to avoid code repetition
# pass integration 'method' parameter, enum?

function left_riemann(f::Function, a::Real, b::Real, n::Int)

    delta = (b - a) / n # Calculate delta between samples
    X = a .+(0:n-1) * delta # Get range of samples
    fX = map(f, X) # Calculate f(x) for each sample in range
    return sum(fX) * delta # Sum and divide by sample width

end

function right_riemann(f::Function, a::Real, b::Real, n::Int)

    delta = (b - a) / n
    X = a .+(1:n) * delta
    fX = map(f, X)
    return sum(fX) * delta

end

function trapezoid(f::Function, a::Real, b::Real, n::Int)

    delta = (b - a) / n
    X = a .+(0:n) * delta
    fX = map(f, X)

    # must calculate mean height of left and right
    return sum(fX[2:end] + diff(fX) * 0.5) * delta

end


text:
# functions for PETSc TS (time stepping) algorithms
export TS, set_ic, set_times, solve!, set_rhs_function, set_rhs_jac
export set_lhs_function, set_lhs_jac
export ComputeRHSFunctionLinear, ComputeRHSJacobianConstant

mutable struct TS{T}
  p::C.TS{T}
  data::Array{Any, 1}  # hold the various ctx tuples
  function TS{T}(p::C.TS{T}, data=nothing; first_instance=false) where {T}
    data_arr = Any[]
    ts = new{T}(p, data_arr)
    push!(ts.data, data)

    if first_instance
      if (!deactivate_finalizers)
        finalizer(PetscDestroy, ts)
      end
    end

    return ts
  end
end

function PetscDestroy(ts::TS{T}) where {T}
  if !PetscFinalized(T) && !isfinalized(ts)
    ts_ref = Ref(ts.p)
    chk(C.TSDestroy(ts_ref))
    ts.p = C.TS{T}(C_NULL)
  end
end

function isfinalized(ts::TS)
  return isfinalized(ts.p)
end

function isfinalized(ts::C.TS)
  return ts.pobj == C_NULL
end

#get the internal KSP object for this TS
function KSP(ts::TS{T}) where {T}
  ksp_c = Ref{C.KSP{T}}()
  chk(C.TSGetKSP(ts.p, ksp_c))
  return KSP{T}(ksp_c[])
end

"""
  Most preferred constructor: take ProblemType, method from options
  database
"""
function TS(::Type{T} ;comm=MPI.COMM_WORLD) where {T<:Scalar}
  ts = Ref{C.TS{T}}()
  chk(C.TSCreate(comm, ts))
  return TS{T}(ts[])
end

"""
  Preferred constructor: set problem type explicitly, get method from
  options database
"""
function TS(::Type{T}, tsptype::C.TSProblemType; comm=MPI.COMM_WORLD) where {T<:Scalar}

  ts = TS(T, comm=comm)
  chk(C.TSSetProblemType(ts.p, tsptype))

  return ts
end

"""
  More explicit constructor: set problem type, method directly

`tsptype` sets the problem type and can be one of the following:
* `TS_LINEAR` - a linear set of ODEs
* `TS_NONLINEAR` - a nonlinear set of ODEs or DEAs

`tstype` sets the method used to solve the problem. More information
about the possible methods is available at the official PETSc [docs](http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSType.html).
"""
function TS(::Type{T}, tsptype::C.TSProblemType, tstype::C.TSType;
                       comm=MPI.COMM_WORLD) where {T<:Scalar}

  ts = TS(T, tsptype, comm=comm)
  chk(C.TSSetType(ts.p, tstype))

  return ts
end

function set_ic(ts::TS{T}, u::Vec{T}) where {T<:Scalar}
  chk(C.TSSetSolution(ts.p, u.p))
end

"""
  Set the times related quantities:
    t0 : initial time value
    dt0: initial time step
    nsteps: maximum number of steps
    tmax: maximum time value
"""
function set_times(ts::TS{T}, t0, dt0,  nsteps::Integer, tmax) where {T<:Scalar}
  TR = real(T)  # PetscReal

  chk(C.TSSetInitialTimeStep(ts.p, TR(t0), TR(dt0)))
  chk(C.TSSetDuration(ts.p, nsteps, TR(tmax)))
end

function petscview(ts::TS{T}) where {T<:Scalar}
  viewer = C.PetscViewer{T}(C_NULL)
  chk(C.TSView(ts.p, viewer))
end

"""
  Solve the system using the initial condition provided in vec
"""
function solve!(ts::TS{T}, vec::Vec{T}) where {T}

  chk(C.TSSolve(ts.p, vec.p))
end

"""
  Solve the system using the intitial condition proived by set_ic
"""
function solve!(ts::TS{T}) where {T}

  vecp = C.Vec{T}(C_NULL)
  chk(C.TSSolve(ts.p, vecp))
end

###############################################################################
# right hand side function

"""
  Wrapper for the right hand side function.  This function is always passed
  to PETSc as the right hand side function, and calls the user supplied
  function internally.  The user supplied function must be the first
  component of the ctx tuple

"""
function rhs_wrapper(ts::C.TS{T}, t, u::C.Vec{T}, F::C.Vec{T}, ctx_ptr::Ptr{Cvoid}) where {T}

  Treal = real(T)
  # transform into high level objects
  bigts = TS{T}(ts, first_instance=false)
  bigu = Vec{T}(u, first_instance=false)
  bigF = Vec{T}(F, first_instance=false)

  ctx = unsafe_pointer_to_objref(ctx_ptr)
  func = ctx[1]
  ctx_inner = ctx[2]  # the user provided ctx

  ret_status = func(bigts, Treal(t), bigu, bigF, ctx_inner)
  return PetscErrorCode(ret_status)
end

"""
  Sets the function that evalutes u_t = g(u, t) for an ODE.
  The function must have the signature:

  f(TS, t, U, F, ctx)

  where TS is a TS object,
  t is the current time
  u is the current state vector
  F is the vector to be populated with u_t
  ctx is the user supplied context tuple (empty tuple if not provided)
"""
function set_rhs_function(ts::TS{T}, r::Vec{T}, f::Function, ctx=()) where {T}

  ctx_outer = (f, ctx)
  push!(ts.data, ctx_outer)
  #ctx_ptr = pointer_from_objref(ctx_outer)
  Treal = real(T)
  # this this pre-compilable?
  #fptr = @cfunction(rhs_wrapper, PetscErrorCode, (C.TS{T}, Treal, C.Vec{T}, C.Vec{T}, Ptr{Cvoid}))
  fptr = @cfunction(rhs_wrapper, PetscErrorCode, (C.TS{T}, T, C.Vec{T}, C.Vec{T}, Ptr{Cvoid}))

  chk(C.TSSetRHSFunction(ts.p, r.p, fptr, ctx_ptr))
end

# a PETSc provided rhs function for the linear, time invarient coefficient
# matrix case

function ComputeRHSFunctionLinear(ts::TS, t, u::Vec, F::Vec, ctx::Tuple)

  # this is amusing, wrapping things just to unwrap them again
  C.TSComputeRHSFunctionLinear(ts.p, t, u.p, F.p, C_NULL)
end

###############################################################################
# right hand side jacobian
function rhs_jac_wrapper(ts::C.TS{T}, t, u::C.Vec{T}, A::C.Mat{T}, B::C.Mat{T}, ctx_ptr::Ptr{Cvoid}) where {T}

  Treal = real(T)
  bigts = TS{T}(ts, first_instance=false)
  bigu = Vec{T}(u, first_instance=false)
  bigA = Mat{T}(A, first_instance=false)
  bigB = Mat{T}(B, first_instance=false)

  ctx = unsafe_pointer_to_objref(ctx_ptr)
  func = ctx[1]
  ctx_inner = ctx[2]

  ret_status = func(bigts, Treal(t), bigu, bigA, bigB, ctx_inner)

  return PetscErrorCode(ret_status)
end

function set_rhs_jac(ts::TS{T}, A::Mat{T}, B::Mat{T}, f::Function, ctx=()) where {T}

  ctx_outer = (f, ctx)
  push!(ts.data, ctx_outer)
  ctx_ptr = pointer_from_objref(ctx_outer)
  Treal = real(T)
  # this this pre-compilable?
  #fptr = @cfunction(rhs_jac_wrapper, PetscErrorCode, (C.TS{T}, Treal, C.Vec{T}, C.Mat{T}, C.Mat{T}, Ptr{Cvoid}))
  fptr = @cfunction(rhs_jac_wrapper, PetscErrorCode, (C.TS{T}, T, C.Vec{T}, C.Mat{T}, C.Mat{T}, Ptr{Cvoid}))

  chk(C.TSSetRHSJacobian(ts.p, A.p, B.p, fptr, ctx_ptr))
end

# a PETSc provided rhs jacobian function for time invarient jacobians
function ComputeRHSJacobianConstant(ts::TS, t, u::Vec, A::Mat, B::Mat, ctx::Tuple)

  chk(C.TSComputeRHSJacobianConstant(ts.p, t, u.p, A.p, B.p, C_NULL))
end

###############################################################################
# left hand side function

function set_lhs_function(ts::TS{T}, res::Vec{T}, f::Function, ctx::Tuple=()) where {T}

  ctx_outer = (f, ctx)
  push!(ts.data, ctx_outer)
  ctx_ptr = pointer_from_objref(ctx_outer)
  Treal = real(T)

  #fptr = @cfunction(lhs_wrapper, PetscErrorCode, (C.TS{T}, Treal, C.Vec{T}, C.Vec{T}, C.Vec{T}, Ptr{Cvoid}))
  fptr = @cfunction(rhs_jac_wrapper, PetscErrorCode, (C.TS{T}, T, C.Vec{T}, C.Mat{T}, C.Mat{T}, Ptr{Cvoid}))

  chk(C.TSSetIFunction(ts.p, res.p, fptr, ctx_ptr))
  return nothing
end

function lhs_wrapper(ts::C.TS{T}, t, u::C.Vec{T}, ut::C.Vec{T}, F::C.Vec{T}, ctx_ptr::Ptr{Cvoid}) where {T}

  Treal = real(T)
  # transform into high level objects
  bigts = TS{T}(ts, first_instance=false)
  bigu = Vec{T}(u, first_instance=false)
  bigut = Vec{T}(ut, first_instance=false)
  bigF = Vec{T}(F, first_instance=false)

  ctx = unsafe_pointer_to_objref(ctx_ptr)
  func = ctx[1]
  ctx_inner = ctx[2]  # the user provided ctx

  ret_status = func(bigts, Treal(t), bigu, bigut, bigF, ctx_inner)
  return PetscErrorCode(ret_status)
end

###############################################################################
# left hand side jacobian

function lhs_jac_wrapper(ts::C.TS{T}, t, u::C.Vec{T}, ut::C.Vec{T}, a, A::C.Mat{T}, B::C.Mat{T}, ctx_ptr::Ptr{Cvoid}) where {T}

  Treal = real(T)
  bigts = TS{T}(ts, first_instance=false)
  bigu = Vec{T}(u, first_instance=false)
  bigut = Vec{T}(ut, first_instance=false)
  bigA = Mat{T}(A, first_instance=false)
  bigB = Mat{T}(B, first_instance=false)

  ctx = unsafe_pointer_to_objref(ctx_ptr)
  func = ctx[1]
  ctx_inner = ctx[2]  # the user provided ctx

  ret_status = func(bigts, bigu, bigut, Treal(a), bigA, bigB, ctx_inner)
  return PetscErrorCode(ret_status)
end

function set_lhs_jac(ts::TS{T}, A::Mat{T}, B::Mat{T}, f::Function, ctx::Tuple=()) where {T}

  ctx_outer = (f, ctx)
  push!(ts.data, ctx_outer)
  ctx_ptr = pointer_from_objref(ctx_outer)
  Treal = real(T)

  #fptr = @cfunction(lhs_jac_wrapper, PetscErrorCode, (C.TS{T}, Treal, C.Vec{T}, C.Vec{T}, Treal, C.Mat{T}, C.Mat{T}, Ptr{Cvoid}))
  fptr = @cfunction(lhs_jac_wrapper, PetscErrorCode, (C.TS{T}, T, C.Vec{T}, C.Vec{T}, T, C.Mat{T}, C.Mat{T}, Ptr{Cvoid}))

  chk(C.TSSetIJacobian(ts.p, A.p, B.p, fptr, ctx_ptr))
end


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

# Import
import Base: convert, size, getindex, setindex!, zero, similar, zero, one, fill, \

# Export types
export RasterBases, NodalCoeff, SolutionCoeff, RasterCoeff, IntermediateCoeff

export gradient

@enum MapTransforms _nb _ng _ns _bn _bg _bs _sn _sb _sg _gn _gs _gb

const Reltypes = Union{Float64, Complex{Float64}}

#
# Coefficent types
#

# Nodal coefficients represent functions in the mesh basis
"""
    NodalCoeff(mesh, coeff)

Return a nodal coefficient array for a basis defined on the `mesh` initialised
with the `coeff` vector of length equal to the number of nodes in the basis.
The element type of `coeff` may be `Float64`, or `Complex{Float64}`.
"""
mutable struct NodalCoeff{T} <: AbstractArray{T, 1}
  mesh::Mesh
  data::Vector{T}
  function NodalCoeff(mesh::Mesh, data::Vector{T}) where {T<:Reltypes}
    @assert length(data) == nodecount(mesh)
    return new{T}(mesh,data)
  end
end

"""
    NodalCoeff(mesh)

Return an uninitialised nodal coefficient array for a basis defined on the `mesh`.
"""
NodalCoeff(mesh::Mesh) = NodalCoeff(mesh, Float64)
NodalCoeff(mesh::Mesh, ::Type{T}) where T = NodalCoeff(mesh, Vector{T}(undef, nodecount(mesh)))

@compat Base.IndexStyle(::Type{<:NodalCoeff}) = IndexLinear()
Base.similar(coeff::NodalCoeff{T}, ::Type{Te}, dims::Dims) where {T,Te} = NodalCoeff(coeff.mesh, T)
size(coeff::NodalCoeff) = (length(coeff.data),)
getindex(coeff::NodalCoeff, i::Int) = coeff.data[i]
setindex!(coeff::NodalCoeff, v, i::Int) = (coeff.data[i] = v)

\(A::SparseMatrixCSC, q::NodalCoeff) = A\q.data

"""
    one(NodalCoeff, mesh)

Return the multiplicative identity element for a basis defined on the `mesh`.
"""
one(::Type{NodalCoeff}, mesh::Mesh) = NodalCoeff(mesh, ones(nodecount(mesh)))

"""
    zero(NodalCoeff, mesh)

Return the additive identity element for a basis defined on the `mesh`.
"""
zero(::Type{NodalCoeff}, mesh::Mesh) = NodalCoeff(mesh, zeros(nodecount(mesh)))

"""
    fill(NodalCoeff, mesh, value)

Return a constant function of `value` for a basis defined on the `mesh`.
"""
function fill(::Type{NodalCoeff}, mesh::Mesh, value::T) where {T}
  NodalCoeff(mesh, fill(value, nodecount(mesh)))
end

# Raster coefficients represent functions in one of the raster bases
@compat abstract type RasterCoeffTypes{T} <: AbstractArray{T, 1} end
@compat Base.IndexStyle(::Type{<:RasterCoeffTypes}) = IndexLinear()
getindex(coeff::T, i::Int) where {T<:RasterCoeffTypes} = coeff.data[i]
setindex!(coeff::T, v, i::Int) where {T<:RasterCoeffTypes} = (coeff.data[i] = v)

"""
    one(RasterCoeffTypes, Raster)

Return the multiplicative identity element for a raster basis defined by `raster`.
"""
one(::Type{T}, rast::Raster) where {T<:RasterCoeffTypes} = T(rast, ones(size(T, rast)...))

"""
    zero(RasterCoeffTypes, mesh)

Return the additive identity element for a raster basis defined by `raster`.
"""
zero(::Type{T}, rast::Raster) where {T<:RasterCoeffTypes} = T(rast, zeros(size(T, rast)...))

"""
    fill(RasterCoeffTypes, mesh, value)

Return a constant function of `value` for a basis defined on the `mesh`.
"""
function one(::Type{T}, rast::Raster, value::Float64) where {T<:RasterCoeffTypes}
  T(rast, fill(size(T, rast)...))
end

# Solution coefficients represent a function expressed in the solution basis,
# which does not include raster points which fall outside of the support of the
# associated mesh.
"""
    SolutionCoeff(raster, coeff)

Return a coefficient array for the solution basis defined by the `raster`,
initialised with the `coeff` vector of length equal to slen(raster). The element
type of `coeff` may be `Float64`, or `Complex{Float64}`.
"""
mutable struct SolutionCoeff{T} <: RasterCoeffTypes{T}
  rast::Raster
  data::Vector{T}
  function SolutionCoeff(rast::Raster, data::Vector{T}) where {T<:Reltypes}
    @assert length(data) == slen(rast)
    return new{T}(rast,data)
  end
end

Base.similar(coeff::SolutionCoeff{T}, ::Type{Te}, dims::Dims) where {T, Te} = SolutionCoeff(coeff.rast, T)
size(coeff::SolutionCoeff) = (length(coeff.data),)
size(::Type{SolutionCoeff}, raster::Raster) = (slen(raster),)

"""
    SolutionCoeff(raster)

Return an uninitialised coefficient array for solution basis defined by the `raster`.
"""
function SolutionCoeff(rast::Raster, ci::NodalCoeff{T}) where {T}
  co = SolutionCoeff(rast, T)
  map!(co,ci)
  return co
end

SolutionCoeff(rast::R) where {R<:Raster} = SolutionCoeff(rast, Float64)
SolutionCoeff(rast::R, ::Type{T}) where {T, R<:Raster} = SolutionCoeff(rast, Vector{T}(undef, slen(rast)))

# Raster coefficients represent a function expressed in the rasterised basis,
# which is defined over a square or cuboid redion, and may include superfluous
# elements which are outside of the support of the mesh.
"""
    RasterCoeff(raster, coeff)

Return a coefficient array for the raster basis defined by the `raster`,
initialised with the `coeff` vector of length equal to blen(raster). The element
type of `coeff` may be `Float64`, or `Complex{Float64}`.
"""
mutable struct RasterCoeff{T} <: RasterCoeffTypes{T}
  rast::Raster
  data::Vector{Float64}
  function RasterCoeff(rast::Raster, data::Vector{T}) where {T<:Reltypes}
    @assert length(data) == blen(rast)
    return new{T}(rast,data)
  end
end

Base.similar(coeff::RasterCoeff{T}, ::Type{Te}, dims::Dims) where {T,Te} = RasterCoeff(coeff.rast, T)
size(coeff::RasterCoeff) = (length(coeff.data),)
size(::Type{RasterCoeff}, raster::Raster) = (blen(raster),)

"""
    RasterCoeff(raster)

Return an uninitialised coefficient array for raster basis defined by the `raster`.
"""
function RasterCoeff(rast::Raster, ci::NodalCoeff{T}) where {T}
  co = RasterCoeff(rast, T)
  map!(co,ci)
  return co
end

RasterCoeff(rast::R) where {R<:Raster} = RasterCoeff(rast, Float64)
RasterCoeff(rast::R, ::Type{T}) where {T, R<:Raster} = RasterCoeff(rast, Vector{T}(undef, blen(rast)))

# Intermediate coefficients represent a function expressed in a higher resolution
# version of the raster basis, and may include superfluous elements which are
# outside of the support of the mesh.
mutable struct IntermediateCoeff{T} <: RasterCoeffTypes{T}
  rast::Raster
  data::Vector{T}
  function IntermediateCoeff(rast::Raster, data::Vector{T}) where {T<:Reltypes}
    @assert length(data) == glen(rast)
    return new{T}(rast,data)
  end
end

Base.similar(coeff::IntermediateCoeff{T}, ::Type{Te}, dims::Dims) where {T,Te} = RasterCoeff(coeff.rast, T)
size(coeff::IntermediateCoeff) = (length(coeff.data),)
size(::Type{IntermediateCoeff}, raster::Raster) = (glen(raster),)

"""
    IntermediateCoeff(raster, coeff)

Return a coefficient array for the intermediate basis defined by the `raster`,
initialised with the `coeff` vector of length equal to glen(raster). The element
type of `coeff` may be `Float64`, or `Complex{Float64}`.
"""
function IntermediateCoeff(rast::Raster, ci::NodalCoeff{T}) where {T}
  co = IntermediateCoeff(rast, T)
  map!(co,ci)
  return co
end

"""
    IntermediateCoeff(raster)

Return an uninitialised coefficient array for intermediate basis defined by the `raster`.
"""
IntermediateCoeff(rast::R) where {R<:Raster} = IntermediateCoeff(rast, Float64)
IntermediateCoeff(rast::R, ::Type{T}) where {T, R<:Raster} = IntermediateCoeff(rast, Vector{T}(undef, glen(rast)))

#
# Coefficient mapping (and construction)
#

# Convert everything to nodal coefficients
function convert(::Type{NodalCoeff}, ci::T) where {T<:RasterCoeffTypes}
  co = NodalCoeff(ci.rast.mesh, eltype(ci.data))
  map!(co, ci)
  return co
end

"""
  map!(out::NodalCoeff, in::RasterCoeffTypes)

Map the function defined by the `input` coefficients defined on a raster basis
to a nodal basis defined on the mesh, overwriting `out` in place.
"""
map!(co::NodalCoeff, ci::SolutionCoeff)  = _map!(co.data, ci.data, ci.rast, _sn)
map!(co::NodalCoeff, ci::RasterCoeff) = _map!(co.data, ci.data, ci.rast, _bn)
map!(co::NodalCoeff, ci::IntermediateCoeff) = _map!(co.data, ci.data, ci.rast, _gn)

# Convert everything to raster coefficients
function convert(::Type{RasterCoeff}, ci::T) where {T<:RasterCoeffTypes}
  co = RasterCoeff(ci.rast, eltype(ci.data))
  map!(co,ci)
  return co
end

map!(co::RasterCoeff, ci::NodalCoeff) = _map!(co.data, ci.data, co.rast, _nb)
map!(co::RasterCoeff, ci::SolutionCoeff)  = _map!(co.data, ci.data, ci.rast, _sb)
map!(co::RasterCoeff, ci::IntermediateCoeff) = _map!(co.data, ci.data, ci.rast, _gb)

# Convert everything to solution coefficients
function convert(::Type{SolutionCoeff}, ci::T) where {T<:RasterCoeffTypes}
  co = SolutionCoeff(ci.rast, eltype(ci.data))
  map!(co,ci)
  return co
end

map!(co::SolutionCoeff, ci::NodalCoeff) = _map!(co.data, ci.data, co.rast, _ns)
map!(co::SolutionCoeff, ci::RasterCoeff)  = _map!(co.data, ci.data, ci.rast, _bs)
map!(co::SolutionCoeff, ci::IntermediateCoeff) = _map!(co.data, ci.data, ci.rast, _gs)

# Convert everything to intermediate coefficients
function convert(::Type{IntermediateCoeff}, ci::T) where {T<:RasterCoeffTypes}
  co = IntermediateCoeff(ci.rast, eltype(ci.data))
  map!(co,ci)
  return co
end

map!(co::IntermediateCoeff, ci::NodalCoeff) = _map!(co.data, ci.data, co.rast, _ng)
map!(co::IntermediateCoeff, ci::RasterCoeff)  = _map!(co.data, ci.data, ci.rast, _bg)
map!(co::IntermediateCoeff, ci::SolutionCoeff) = _map!(co.data, ci.data, ci.rast, _sg)

# Low level mapping function
function _map!(ovec::Vector{Float64},
               ivec::Vector{Float64},
               rast::Raster,
               mode::MapTransforms)

  modeint = Int(mode)
  ovecptr = pointer(ovec)
  ivecptr = pointer(ivec)
  rastptr = rast.ptr

  ilen = length(ivec)
  olen = length(ovec)

  icxx"""

      RVector iprm($(ilen), $(ivecptr), SHALLOW_COPY);
      RVector oprm($(olen), $(ovecptr), SHALLOW_COPY);

      switch($(modeint))
      {
          case 0:
              $(rastptr)->Map_MeshToBasis(iprm, oprm);
              break;
          case 1:
              $(rastptr)->Map_MeshToGrid(iprm, oprm);
              break;
          case 2:
              $(rastptr)->Map_MeshToSol(iprm, oprm);
              break;

          case 3:
              $(rastptr)->Map_BasisToMesh(iprm, oprm);
              break;
          case 4:
              $(rastptr)->Map_BasisToGrid(iprm, oprm);
              break;
          case 5:
              $(rastptr)->Map_BasisToSol(iprm, oprm);
              break;

          case 6:
              $(rastptr)->Map_SolToMesh(iprm, oprm);
              break;
          case 7:
              $(rastptr)->Map_SolToBasis(iprm, oprm);
              break;
          case 8:
              $(rastptr)->Map_SolToGrid(iprm, oprm);
              break;

          case 9:
              $(rastptr)->Map_GridToMesh(iprm, oprm);
              break;
          case 10:
              $(rastptr)->Map_GridToSol(iprm, oprm);
              break;
          case 11:
              $(rastptr)->Map_GridToBasis(iprm, oprm);
              break;
      }
  """

end

# Low level mapping function
function _map!(ovec::Vector{Complex{Float64}},
               ivec::Vector{Complex{Float64}},
               rast::Raster,
               mode::MapTransforms)

  modeint = Int(mode)
  rastptr = rast.ptr

  ilen = length(ivec)

  rivec = real(ivec)
  iivec = imag(ivec)
  rivecptr = pointer(rivec)
  iivecptr = pointer(iivec)

  olen = length(ovec)

  rovec = real(ovec)
  iovec = imag(ovec)
  rovecptr = pointer(rovec)
  iovecptr = pointer(iovec)

  icxx"""

      CVector iprm($(ilen));
      CVector oprm($(olen));

      std::complex<double> *val;

      val = iprm.data_buffer();

      for (int i = 0; i < $(ilen); i++)
          val[i] = std::complex<double> ($(rivecptr)[i], $(iivecptr)[i]);

      switch($(modeint))
      {
          case 0:
              $(rastptr)->Map_MeshToBasis(iprm, oprm);
              break;
          case 1:
              $(rastptr)->Map_MeshToGrid(iprm, oprm);
              break;
          case 2:
              $(rastptr)->Map_MeshToSol(iprm, oprm);
              break;

          case 3:
              $(rastptr)->Map_BasisToMesh(iprm, oprm);
              break;
          case 4:
              $(rastptr)->Map_BasisToGrid(iprm, oprm);
              break;
          case 5:
              $(rastptr)->Map_BasisToSol(iprm, oprm);
              break;

          case 6:
              $(rastptr)->Map_SolToMesh(iprm, oprm);
              break;
          case 7:
              $(rastptr)->Map_SolToBasis(iprm, oprm);
              break;
          case 8:
              $(rastptr)->Map_SolToGrid(iprm, oprm);
              break;

          case 9:
              $(rastptr)->Map_GridToMesh(iprm, oprm);
              break;
          case 10:
              $(rastptr)->Map_GridToSol(iprm, oprm);
              break;
          case 11:
              $(rastptr)->Map_GridToBasis(iprm, oprm);
              break;
      }

      val = oprm.data_buffer();

      for (int i = 0; i < $(olen); i++)
      {
        $(rovecptr)[i] = real(val[i]);
        $(iovecptr)[i] = imag(val[i]);
      }
  """

  ovec .= rovec .+ im.*iovec

  return ovec

end

"""
    gradient(coeff)

Compute the spatial gradient of `coeff` which must be expressed in an
`IntermediateCoeff` basis.
"""
function gradient(coeff::IntermediateCoeff)

  # TODO: Native implementaiton from ADMM.jl
  ndim = dimensions(coeff.rast.mesh)
  len = glen(coeff.rast)
  rptr = coeff.rast.ptr
  coptr = pointer(coeff.data)

  ∇coeff = Array{Float64}(undef, len, ndim)
  gcoptr = pointer(∇coeff)

  icxx"""
    const IVector &gdim = $(rptr)->GDim();
    const RVector &gsize = $(rptr)->GSize();

    RVector img($(len), $(coptr), SHALLOW_COPY);
    RVector *imgrad = new RVector[$(ndim)];
    ImageGradient(gdim, gsize, img, imgrad, $(rptr)->Elref());

    for(int i = 0; i < $(len); i++)
      $(gcoptr)[i] = imgrad[0][i];

    for(int i=0; i < $(len); i++)
      $(gcoptr)[i+$(len)] = imgrad[1][i];

    if($ndim > 2)
    {
      for(int i=0; i < $(len); i++)
        $(gcoptr)[i+(2*$(len))] = imgrad[2][i];
    }

    delete []imgrad;
  """

  return ∇coeff

end

function gradient(coeff::IntermediateCoeff{Complex{Float64}})

  # TODO: Native implementaiton from ADMM.jl
  ndim = dimensions(coeff.rast.mesh)
  len = glen(coeff.rast)
  rptr = coeff.rast.ptr

  rcoeff = real(coeff.data)
  icoeff = imag(coeff.data)
  rcoptr = pointer(rcoeff)
  icoptr = pointer(icoeff)

  ∇coeff = Array{Float64}(undef, len*2, ndim)
  gcoptr = pointer(∇coeff)

  icxx"""
    const IVector &gdim = $(rptr)->GDim();
    const RVector &gsize = $(rptr)->GSize();

    CVector img($(len));

    std::complex<double> *val;

    val = img.data_buffer();

    for (int i = 0; i < $(len); i++)
        val[i] = std::complex<double> ($(rcoptr)[i], $(icoptr)[i]);

    CVector *imgrad = new CVector[$(ndim)];
    ImageGradient(gdim, gsize, img, imgrad, $(rptr)->Elref());

    for(int i = 0; i < $(len); i++)
    {
      $(gcoptr)[2*i]   = real(imgrad[0][i]);
      $(gcoptr)[2*i+1] = imag(imgrad[0][i]);
    }

    for(int i=0; i < $(len); i++)
    {
      $(gcoptr)[2*i+(2*$(len))]   = real(imgrad[1][i]);
      $(gcoptr)[2*i+(2*$(len))+1] = imag(imgrad[1][i]);
    }

    if($ndim > 2)
    {
      for(int i=0; i < $(len); i++)
      {
        $(gcoptr)[2*i+(4*$(len))]   = real(imgrad[2][i]);
        $(gcoptr)[2*i+(4*$(len))+1] = real(imgrad[2][i]);
      }
    }

    delete []imgrad;
  """

  ∇coeffout = ∇coeff[1:2:end, :] + im*∇coeff[2:2:end, :]
  return ∇coeffout

end


text:
@testset "node2vec" begin
    clusters = [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

    g = smallgraph(:karate)
    fg = FeaturedGraph(g)
    vectors = node2vec(fg; walks_per_node=10, len=80, p=1.0, q=1.0)
    R = kmeans(vectors, 2)

    learned_clusters = copy(assignments(R))
    # ensure that the cluster containing node 1 is cluster 1
    if assignments(R)[1] != 1
        learned_clusters = [i == 1 ? 2 : 1 for i in assignments(R)]
    end

    incorrect = sum(learned_clusters .!= clusters)
    @test incorrect < 4
end


text:
function update_weights(data, e, closest_same, closest_other, weights, weights_mult, closest_same_weights, closest_other_weights, m, k, max_f_vals, min_f_vals)

	# Penalty term
	penalty = sum(closest_same_weights .* (abs.(e .- closest_same)./((max_f_vals .- min_f_vals) .+ eps(Float64))), dims=1)

	# Reward term
	reward = sum(weights_mult .* (closest_other_weights .* (abs.(e .- closest_other)./((max_f_vals .- min_f_vals) .+ eps(Float64)))), dims=1)

	# Weights update
	weights = weights .- penalty./m .+ reward./m

	# Return updated weights.
	return vec(weights)
end

max_stars_repo_name: Architecton/diplomska-naloga-koda

max_issues_repo_name: Architecton/diplomska-naloga-koda

max_forks_repo_name: Architecton/diplomska-naloga-koda


text:
using Projections, CSV, DataFrames, Plots, Distributions, LinearAlgebra, Random

path     = @__DIR__
savepath = path

include(joinpath(path, "experiments.jl"))

## initialization
save = false

l1 = unique(vcat(2:0.1:3, 3:0.05:6));
l2 = unique(vcat(2:0.1:3, 3:0.05:4));
N1 = @. ceil(Int64, 10^l1)
N2 = @. ceil(Int64, 10^l2)

## benchmarks
comparison_solver(Our(), N1; save = false, maxevals = 1);
@time table1 = comparison_solver(Our(), N1; save = save,  savepath = savepath, maxevals = 100);

comparison_solver(General(), N2; save = false, maxevals = 1);
@time table2 = comparison_solver(General(), N2; save = save,  savepath = savepath, maxevals = 100);

comparison_philpott(N2; save = false, maxevals = 1);
@time table3 = comparison_philpott(N2; save = save, savepath = savepath, maxevals = 100);

## tables
savetable(vcat(table1, table2, table3), save = save, savepath = savepath, savename = "comparison")

## Kullback-Leibler divergence
m = model_DRO(10, KullbackLeibler(), 0.1)
λ = 1:0.01:3
h(λ,m; save = save, savepath = savepath)

## Burg entropy
m = model_DRO(10, Burg(), 0.1)
λ = 1.9:0.01:3
h(λ, m; save = save, savepath = savepath)

## Hellinger distance
m = model_DRO(10, Hellinger(), 0.1)
λ = 1.9:0.01:5
h(λ, m; save = save, savepath = savepath)

## χ²-distance distance
m = model_DRO(10, ChiSquare(), 0.1)
λ = 1.9:0.01:3
h(λ, m; save = save, savepath = savepath)

## Modified χ²-distance distance
m = model_DRO(10, ModifiedChiSquare(), 0.1)
λ = -1.8:0.01:3
h(λ, m; save = save, savepath = savepath)

## l-2 norm
m = model_DRO(10, Ltwo(), 0.1)
λ = 0:0.01:30
h(λ, m; save = save, savepath = savepath)

max_stars_repo_name: VaclavMacha/Projections

max_issues_repo_name: VaclavMacha/Projections

max_forks_repo_name: VaclavMacha/Projections


text:
function pprint(io::GarishIO, ::MIME"text/plain", t::Tuple)
    pprint_list_like(io, t, "(", ")")
end

function pprint(io::GarishIO, ::MIME"text/plain", s::Pair)
    pprint(io, s.first)
    print_operator(io, "=>")
    if !io.compact
        buf = IOBuffer()
        pprint(GarishIO(IOContext(buf, :color=>false), io), s.first)
        io.state.offset = length(String(take!(buf))) + 4
    end

    within_nextlevel(io) do
        upperlevel_type = io.state.type
        upperlevel_noindent_in_first_line = io.state.noindent_in_first_line
        io.state.type = StructField
        io.state.noindent_in_first_line = true
        pprint(io, s.second)
        io.state.noindent_in_first_line = upperlevel_noindent_in_first_line
        io.state.type = upperlevel_type
    end
    return
end

function pprint(io::GarishIO, ::MIME"text/plain", bool::Bool)
    print_token(show, io, :constant, bool)
end

function pprint(io::GarishIO, ::MIME"text/plain", s::AbstractString)
    print_token(show, io, :string, s)
end

function pprint(io::GarishIO, ::MIME"text/plain", ::UndefInitializer)
    print_token(io, :undef, "undef")
end

function pprint(io::GarishIO, ::MIME"text/plain", ::Nothing)
    print_token(io, :constant, "nothing")
end

function pprint(io::GarishIO, ::MIME"text/plain", ::Missing)
    print_token(io, :constant, "missing")
end

"""
    print_indent(io::GarishIO)

Print an indentation. This should be only used under `MIME"text/plain"` or equivalent.
"""
function print_indent(io::GarishIO)
    io.compact && return
    io.state.level > 0 || return

    io.show_indent || return print(io, " "^(io.indent * io.state.level))
    for _ in 1:io.state.level
        print_token(io, :comment, "│")
        print(io, " "^(io.indent - 1))
    end
end

"""
    print_operator(io::GarishIO, op)

Print an operator, such as `=`, `+`, `=>` etc. This should be only used under `MIME"text/plain"` or equivalent.
"""
function print_operator(io::GarishIO, op)
    io.compact || print(io, " ")
    print_token(io, :operator, op)
    io.compact || print(io, " ")
end


text:
@testset "operations" begin

    @testset "offsetcombine" begin
        r = rand(10)
        x = rand(10)
        y = rand(10)

        @test offsetcombine!(r, +, x, y) ≈ x .+ y
        @test offsetcombine!(r, *, x, y) ≈ x .* y

        @test offsetcombine!(r, +, x, y, 2) ≈ x .+ circshift(y, 2)
        @test offsetcombine!(r, *, x, y, 2) ≈ x .* circshift(y, 2)

        @test offsetcombine(+, x, y) ≈ x .+ y
        @test offsetcombine(*, x, y) ≈ x .* y

        @test offsetcombine(+, x, y, 2) ≈ x .+ circshift(y, 2)
        @test offsetcombine(*, x, y, 2) ≈ x .* circshift(y, 2)
    end

    @testset "Bipolar" begin
        vectors = [[-1, 1, 1, -1],
                   [1, -1, -1, 1],
                   [-1, 1, 1, -1]]

        hdvs = BipolarHDV.(vectors)

        @test (hdvs[1] + hdvs[2] .== [0, 0, 0, 0]) |> all
        @test aggregate(hdvs) == last(hdvs)

        @test (hdvs[1] * hdvs[2] .== [-1, -1, -1, -1]) |> all
        @test bind(hdvs) == BipolarHDV([1, -1, -1, 1])

        hdv = first(hdvs)
        v = first(vectors) |> copy

        @test Π(hdv, 2) == circshift(v, 2)
        @test Π!(hdv, 2) == circshift(v, 2)
        @test resetoffset!(hdv) == v
        @test hdv.v == v
    end

    @testset "Binary" begin
        vectors = [[true, false, false, true],
                   [true, false, false, true],
                   [false, true, true, false]]

        hdvs = BinaryHDV.(vectors)

        @test aggregate(hdvs) == BinaryHDV([true, false, false, true])

        @test (hdvs[1] * hdvs[2] .== [false, false, false, false]) |> all
        @test bind(hdvs[2:3]) == BinaryHDV([true, true, true, true])

        hdv = first(hdvs)
        v = first(vectors) |> copy

        @test Π(hdv, 2) == circshift(v, 2)
        @test Π!(hdv, 2) == circshift(v, 2)
        @test resetoffset!(hdv) == v
        @test hdv.v == v
    end

    @testset "GradedBipolarHDV" begin
        vectors = [[0, 0, 0, 0.0],
                   [0.8, 0.1, -0.3, 0.2],
                   [-1.0, 1, 1, -1]]

        hdvs = GradedBipolarHDV.(vectors)

        @test hdvs[1] + hdvs[2] ≈ hdvs[2]
        @test aggregate(hdvs) ≈ last(hdvs)

        @test hdvs[1] * hdvs[2] ≈ hdvs[1]
        @test bind(hdvs) ≈ BipolarHDV([0.0, 0.0, 0.0, 0.0])

        hdv = first(hdvs)
        v = first(vectors) |> copy

        @test Π(hdv, 2) == circshift(v, 2)
        @test Π!(hdv, 2) == circshift(v, 2)
        @test resetoffset!(hdv) == v
        @test hdv.v == v
    end

    @testset "GradedHDV" begin
        vectors = [[0.5, 0.5, 0.5, 0.5],
                   [0.8, 0.1, 0.3, 0.2],
                   [1.0, 0, 1, 0]]

        hdvs = GradedHDV.(vectors)

        @test hdvs[1] + hdvs[2] ≈ hdvs[2]
        @test aggregate(hdvs) ≈ last(hdvs)

        @test hdvs[1] * hdvs[2] ≈ hdvs[1]
        @test bind(hdvs) ≈ hdvs[1]

        hdv = first(hdvs)
        v = first(vectors) |> copy

        @test Π(hdv, 2) == circshift(v, 2)
        @test Π!(hdv, 2) == circshift(v, 2)
        @test resetoffset!(hdv) == v
        @test hdv.v == v
    end

    @testset "RealHDV" begin
        vectors = [[0, 0, 0, 0.0],
                   [0.8, 0.1, -0.3, 0.2],
                   [-1.0, 1, 1, -1]]

        hdvs = RealHDV.(vectors)

        @test hdvs[1] + hdvs[2] ≈ hdvs[2] / sqrt(2)
        @test aggregate(hdvs) ≈ RealHDV(sum(vectors) / √(3))

        @test hdvs[1] * hdvs[2] ≈ hdvs[1]
        @test bind(hdvs) ≈ hdvs[1] .* hdvs[2] .* hdvs[3]

        hdv = first(hdvs)
        v = first(vectors) |> copy

        @test Π(hdv, 2) == circshift(v, 2)
        @test Π!(hdv, 2) == circshift(v, 2)
        @test resetoffset!(hdv) == v
        @test hdv.v == v
    end
end


text:
"""
The gene-reaction rules (grr) are written as
[[gene1 and gene2...] or [gene3 and...] ...]
so for a reaction to be available it is sufficient that one group
is available, but inside a group all of the genes need to be available
"""

@testset "knockout_single_gene" begin
    m = StandardModel()
    add_metabolite!(m, Metabolite("A"))
    add_metabolite!(m, Metabolite("B"))
    add_gene!(m, Gene("g1"))
    add_gene!(m, Gene("g2"))
    add_reaction!(
        m,
        Reaction("v1", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"]]),
    )
    add_reaction!(
        m,
        Reaction("v2", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1", "g2"]]),
    )
    add_reaction!(
        m,
        Reaction("v3", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"], ["g2"]]),
    )
    add_reaction!(
        m,
        Reaction(
            "v4",
            metabolites = Dict("A" => -1.0, "B" => 1.0),
            grr = [["g1", "g2"], ["g2"]],
        ),
    )

    remove_gene!(m, "g1", knockout_reactions = true)

    @test length(m.reactions) == 2
    @test !haskey(m.reactions, "v1")
    @test !haskey(m.reactions, "v2")
end

@testset "knockout_multiple_genes" begin
    m = StandardModel()
    add_metabolite!(m, Metabolite("A"))
    add_metabolite!(m, Metabolite("B"))
    add_gene!(m, Gene("g1"))
    add_gene!(m, Gene("g2"))
    add_gene!(m, Gene("g3"))
    add_reaction!(
        m,
        Reaction(
            "v1",
            metabolites = Dict("A" => -1.0, "B" => 1.0),
            grr = [["g1"], ["g2"], ["g3"]],
        ),
    )
    add_reaction!(
        m,
        Reaction("v2", metabolites = Dict("A" => 1.0, "B" => -1.0), grr = [["g1"], ["g3"]]),
    )

    remove_genes!(m, ["g1", "g3"], knockout_reactions = true)

    @test haskey(m.reactions, "v1")
    @test !haskey(m.reactions, "v2")
end


text:
# vim: set ft=julia sw=4 ts=4:
#
# Fibonacci.jl
#
# Calculate Fibonacci numbers in Julia

module Fibonacci

export fib_naive, fib

"Naive implementation of Fibonacci number calculation"
function fib_naive(n)
	if n < 2
		n
	else
		fib_naive(n-2) + fib_naive(n-1)
	end
end

"Tail-call optimized implementation"
function fib(n)
	function aux(n, a, b)
		if n < 1
			a
		else
			aux(n-1, b, a+b)
		end
	end
	aux(n, BigInt(0), BigInt(1))
end

end # module

max_stars_repo_name: kkirstein/proglang-playground

max_issues_repo_name: kkirstein/proglang-playground

max_forks_repo_name: kkirstein/proglang-playground


text:
_scitypes(X) = NamedTuple{schema(X).names}(schema(X).scitypes)

abstract type MLJTask <: MLJType end # `Task` already has meaning in Base

mutable struct UnsupervisedTask <: MLJTask
    X
    input_scitypes::NamedTuple
    input_scitype
end

function UnsupervisedTask(; data=nothing, ignore=Symbol[], verbosity=1)

    data != nothing || error("You must specify data=... ")

    if ignore isa Symbol
        ignore = [ignore, ]
    end

    names = schema(data).names |> collect

    issubset(ignore, names) || error("ignore=$ignore contains a column name "*
                                     "not encountered in supplied data.")

    features = filter(names) do ftr
        !(ftr in ignore)
    end

    if length(features) == 1
        input_is_multivariate = false
        X = selectcols(data, features[1]) # an abstract vector
    else
        input_is_multivariate = true
        X = selectcols(data, features) # a table
    end

    input_scitype = scitype(X)

    if input_is_multivariate
        input_scitypes = _scitypes(X)
    else
        input_scitypes = (input=scitype_union(X), )
    end

    if Unknown <: input_scitype
        @warn "An input feature with unknown scitype has been encountered. "*
        "Check the input_scitypes field of your task. "

    end

    if verbosity > 0
        @info "input_scitype = $input_scitype"
    end

    return UnsupervisedTask(X, input_scitypes, input_scitype)
end

# X and y can be different views of a common object.
mutable struct SupervisedTask <: MLJTask
    X
    y
    is_probabilistic
    input_scitypes
    target
    input_scitype
    target_scitype
end

"""
    contains_unknown(X)

In the case that X is a table, determine if one of the columns has
elements of `Unknown` scitype. In the case of an `AbstractVector` of
tuples, determine if one of the tuples has an element of `Unknown`
scitype. In the case of any other `AbstractVector`, determine if an
element has `Unknown` scitype. In all other cases, return `false`.

"""
function contains_unknown(X)
    target_scitype = scitype(X)
    if target_scitype <: Table
        types = values(_scitypes(X))
        unknown_detected = !all([!(Unknown <: t) for t in types])
    elseif target_scitype <: AbstractArray{<:Tuple}
        types = scitype_union(X).types
        unknown_detected = !all([!(Unknown <: t) for t in types])
    elseif target_scitype <: AbstractArray
        unknown_detected = Unknown <: scitype_union(X)
    else
        unknown_detected = false
    end
    return unknown_detected
end

function SupervisedTask(X, y::AbstractVector;
                        is_probabilistic=false,
                        target=:target,
                        verbosity=1)

    # is_probabilistic != nothing ||
    #     error("You must specify is_probabilistic=true or is_probabilistic=false. ")

    target_scitype = scitype(y)
    input_scitype = scitype(X)

    if input_scitype <: Table
        input_scitypes = _scitypes(X)
    else
        input_scitypes = (input=scitype_union(X),)
    end

    unknown_detected = contains_unknown(input_scitype)||
    contains_unknown(target_scitype)

    if unknown_detected
        @warn "An Unknown scitype has been encountered. "*
        "Check the input_scitypes and target_scitype fields of "*
        "your task. "
    end

    if verbosity > 0
        @info "\nis_probabilistic = $is_probabilistic\ninput_scitype = "*
        "$input_scitype \ntarget_scitype = $target_scitype"
    end

    return SupervisedTask(X, y,
                          is_probabilistic,
                          input_scitypes,
                          target,
                          input_scitype, target_scitype)
end

function SupervisedTask(; data=nothing, is_probabilistic=false, target=nothing, ignore=Symbol[], verbosity=1)

    target_as_input = target

    data != nothing ||
        error("You must specify data=... or use SupervisedTask(X, y, ...).")

    is_probabilistic != nothing ||
        error("You must specify is_probabilistic=true or is_proabilistic=false. ")

    if ignore isa Symbol
        ignore = [ignore, ]
    end

    names = schema(data).names |> collect

    issubset(ignore, names) ||
        error("ignore=$ignore contains a column name not encountered "*
              "in supplied data.")

    target != nothing ||
        error("You must specify target=... (use Symbol or Vector{Symbol}) ")

    if target isa Vector
        issubset(target, names) ||
            error("One or more specified target columns missing from "*
                  "supplied data. ")
        yrowtable = Tables.rowtable(selectcols(data, target))
        y = values.(yrowtable)
    else
        target in schema(data).names ||
            error("Specificed target not found in data. ")
        y = selectcols(data, target)
        target = [target, ]
    end

    features = filter(names) do ftr
        !(ftr in target) && !(ftr in ignore)
    end

    if length(features) == 1
        X = selectcols(data, features[1])
        input_is_multivariate = false
    else
        X = selectcols(data, features)
        input_is_multivariate = true
    end

    return SupervisedTask(X, y;
                          is_probabilistic=is_probabilistic,
                          target=target_as_input,
                          verbosity=verbosity)
end

## RUDIMENTARY TASK OPERATIONS

nrows(task::MLJTask) = nrows(task.X)
Base.eachindex(task::MLJTask) = Base.OneTo(nrows(task))
nfeatures(task::MLJTask) = length(schema(task.X).names)

X_(task::MLJTask) = task.X
y_(task::SupervisedTask) = task.y

X_and_y(task::SupervisedTask) = (X_(task), y_(task))

# make tasks callable:
(task::UnsupervisedTask)() = X_(task)
(task::SupervisedTask)() = X_and_y(task)


text:
# dirpath = "C:/data/Performance_All/"
# largest_file = "Performance_2000Q4.txt"
# outpath = "c:/data"
#

# download("https://packages.revolutionanalytics.com/datasets/AirOnTime87to12/AirOnTimeCSV.zip")
# ;unzip c:/data/AirOnTimeCSV AirOnTimeCSV.zip

# Uncomment for debugging
dirpath = "c:/data/AirOnTimeCSV/"
largest_file = "airOT199302.csv"
outpath = "c:/data/jdf-bench/airOT199302.csv"
data_label = "Air On Time 199302"
delim = ','
header = true

include("C:/Users/RTX2080/git/JDF/benchmarks/benchmarks.jl")
@time res = gen_benchmark("c:/data/AirOnTimeCSV/", "airOT199302.csv", "c:/data/jdf-bench/airOT199302.csv", data_label)
# gen_benchmark("C:/data/Performance_All/", "Performance_2000Q4.txt", "c:/data/jdf-bench/Performance_2000Q4.txt", "Fannie Mae Performance 2000Q4")

sizedf = sum_file_size(outpath)

using StatsPlots

p = plot(
    sizedf.pkg,
    sizedf.fs/1024^3,
    linetype = :bar,
    ylab = "Size (GB)",
    legend = false,
    title = "On-disk file Size for various formats\n $data_label data")
savefig(p, joinpath(outpath, largest_file*"_filesize.png"))


text:
"""
    add_constraints!(ata_model::AbstractModel; constraints_file = "constraints.csv", constraints_delim = ";")

# Description

Add categorical and sum constraints to the `ATA.AbstractModel` as specified in the `constraints_file`.
Alternatively, the constraint dataframe can be passed using the argument `constraints`.
It requires the [`start_ata`](#ATA.start_ata) step.

# Arguments

- **`ata_model::AbstractModel`** : Required. The model built with `start_ata()` and with settings loaded by [`start_ata`](#ATA.start_ata) function.
- **`constraints::DataFrames.DataFrame = DataFrames.DataFrame()**`. Optional. Default: DataFrames.DataFrame(). The dataframe containing the categorical and quantitative constraints.
- **`constraints_file`** : Optional. Default: "constraints.csv". The path of the file containing the categorical and sum constraints in the form of custom-separated values.
- **`constraints_delim`** : Optional. Default: ";". The custom-separator for the `constraints_file`.

"""
function add_constraints!(
    ata_model::AbstractModel;
    constraints::DataFrames.DataFrame = DataFrames.DataFrame(),
    constraints_file = "constraints.csv",
    constraints_delim = ";",
)
    if size(constraints, 1) > 0
        categorical_constraints = constraints
        success!(ata_model, "Constraints dataframe loaded.")
    elseif !isfile(constraints_file)
        error!(
            ata_model,
            string(
                "Constraints file with name ",
                constraints_file,
                " does not exist. Provide a valid constraints dataframe or a name of an existing file.",
            ),
        )
        return nothing
    else
        try
            categorical_constraints =
                CSV.read(constraints_file, DataFrames.DataFrame, delim = constraints_delim)
        catch e
            error!(ata_model, "Error in reading constraints file.")
            return nothing
        end
        success!(ata_model, "Constraints file loaded.")
    end
    try
        n_items = ata_model.settings.n_items
        bank = ata_model.settings.bank
        A = Vector{Matrix{Float64}}(undef, ata_model.settings.T)
        b = Vector{Vector{Float64}}(undef, ata_model.settings.T)
        if size(ata_model.constraints[1].constr_b, 1) > 0
            for t = 1:ata_model.settings.T
                A[t] = copy(ata_model.constraints[t].constr_A)
                b[t] = copy(ata_model.constraints[t].constr_b)
            end
        else
            for t = 1:ata_model.settings.T
                A[t] = zeros(Float64, 0, ata_model.settings.n_items)
                b[t] = Vector{Float64}(undef, 0)
            end
        end
        x_forced0 = ata_model.settings.forced0
        if size(categorical_constraints, 1) == 0
            error!(
                ata_model,
                string("No lines in file ", constraints_file, " nothing added."),
            )
        else
            categorical_constraints.var = Symbol.(categorical_constraints.var)
            CatCons = Vector{Vector{Float64}}(undef, ata_model.settings.T)
            for g = 1:ata_model.settings.n_groups
                if g > 1
                    tests = collect(
                        (sum(
                            ata_model.settings.Tg[1:(g-1)],
                        )+1):(sum(
                            ata_model.settings.Tg[1:(g-1)],
                        )+ata_model.settings.Tg[g]),
                    )
                else
                    tests = collect(1:ata_model.settings.Tg[1])
                end
                GroupCatCons = findall((categorical_constraints.group .== g))
                for con in GroupCatCons
                    for t in tests
                        if !ismissing(categorical_constraints[con, :min])
                            indices =
                                string.(
                                    bank[!, Symbol.(categorical_constraints[con, :var])],
                                ) .== string(categorical_constraints[con, :value])
                            indices[ismissing.(indices)] .= false
                            indices = findall(indices .== true)
                            A[t] = vcat(A[t], zeros(Float64, 1, n_items))
                            if ismissing(categorical_constraints[con, :weight])
                                A[t][end, indices] .= -1.0
                                push!(b[t], -categorical_constraints[con, :min])
                            else
                                A[t][end, indices] .=
                                    -(categorical_constraints[con, :weight])
                                push!(
                                    b[t],
                                    -categorical_constraints[con, :weight] *
                                    categorical_constraints[con, :min],
                                )
                            end
                        end
                        if !ismissing(categorical_constraints[con, :max])
                            indices =
                                string.(
                                    ata_model.settings.bank[
                                        !,
                                        Symbol.(categorical_constraints[con, :var]),
                                    ],
                                ) .== string(categorical_constraints[con, :value])
                            indices[ismissing.(indices)] .= false
                            indices = findall(indices .== true)
                            if categorical_constraints[con, :max] == 0
                                x_forced0[t][indices] .= false
                            else
                                A[t] = vcat(A[t], zeros(Float64, 1, n_items))
                                if ismissing(categorical_constraints[con, :weight])
                                    A[t][end, indices] .= 1.0
                                    push!(b[t], categorical_constraints[con, :max])
                                else
                                    A[t][end, indices] .=
                                        (categorical_constraints[con, :weight])
                                    push!(
                                        b[t],
                                        categorical_constraints[con, :weight] *
                                        categorical_constraints[con, :max],
                                    )
                                end
                            end
                        end
                    end
                end
            end
            success!(ata_model, "Linear constraints (categorical and quantitative) added.")
        end

        for t = 1:ata_model.settings.T
            DelimitedFiles.writedlm("opt/A_$t.csv", A[t])
            DelimitedFiles.writedlm("opt/b_$t.csv", b[t])
        end
        DelimitedFiles.writedlm("opt/x_forced0.txt", x_forced0)
        #update model
        ata_model.settings.forced0 = x_forced0
        for t = 1:ata_model.settings.T
            ata_model.constraints[t].constr_b = b[t]
            ata_model.constraints[t].constr_A = A[t]
        end
        JLD2.@save "opt/ata_model.jld2" ata_model

    catch e
        error!(ata_model, string("", sprint(showerror, e)))
    end
    return nothing
end


text:
"""
Functor to update JS part when an observable changes.
We make this a Functor, so we can clearly identify it and don't sent
any updates, if the JS side requires to update an Observable
(so we don't get an endless update cycle)
"""
struct JSUpdateObservable
    session::Session
    id::String
end

function (x::JSUpdateObservable)(value)
    # Sent an update event
    send(x.session, payload = value, id = x.id, type = UpdateObservable)
end

"""
Update the value of an observable, without sending changes to the JS frontend.
This will be used to update updates from the forntend.
"""
function update_nocycle!(obs::Observable, value)
    setindex!(
        obs, value,
        notify = (f-> !(f isa JSUpdateObservable))
    )
end

"""
Register all Observables of a session with the connected Javascript frontend.
It updates all the values on the JS side accordingly,
and enables two way communication
"""
function register_obs!(session::Session)
    for (id, (registered, observable)) in session.observables
        if !registered
            # Register on the JS side by sending the current value
            send(
                session,
                type = UpdateObservable,
                id = id,
                payload = observable[]
            )
            updater = JSUpdateObservable(session, id)
            # Make sure we update the Javascript values!
            on(updater, observable)
            session.observables[id] = (true, observable)
        end
    end
    return
end

function jsrender(session::Session, obs::Observable)
    html = map(obs) do value
        repr_richest(value)
    end
    dom = DOM.m_unesc("span", html[])
    onjs(session, html, js"""
        function (html){
            var dom = $(dom);
            if(dom){
                dom.innerHTML = html;
                return true;
            }else{
                //deregister the callback if the observable dom is gone
                return false;
            }
        }
    """)
    return dom
end


text:
@testset "kernelproduct" begin
    rng = MersenneTwister(123456)
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)

    k1 = LinearKernel()
    k2 = SqExponentialKernel()
    k3 = RationalQuadraticKernel()

    k = KernelProduct([k1, k2])
    @test length(k) == 2
    @test k(v1, v2) == (k1 * k2)(v1, v2)
    @test (k * k)(v1, v2) ≈ k(v1, v2)^2
    @test (k * k3)(v1, v2) ≈ (k3 * k)(v1, v2)

    @testset "kernelmatrix" begin
        rng = MersenneTwister(123456)

        Nx = 5
        Ny = 4
        D = 3

        w1 = rand(rng) + 1e-3
        w2 = rand(rng) + 1e-3
        k1 = SqExponentialKernel()
        k2 = LinearKernel()
        k = k1 * k2

        @testset "$(typeof(x))" for (x, y) in [
            (ColVecs(randn(rng, D, Nx)), ColVecs(randn(rng, D, Ny))),
            (RowVecs(randn(rng, Nx, D)), RowVecs(randn(rng, Ny, D))),
        ]
            @test kernelmatrix(k, x, y) ≈ kernelmatrix(k1, x, y) .* kernelmatrix(k2, x, y)

            @test kernelmatrix(k, x) ≈ kernelmatrix(k1, x) .* kernelmatrix(k2, x)

            K_diag_manual = kerneldiagmatrix(k1, x) .* kerneldiagmatrix(k2, x)
            @test kerneldiagmatrix(k, x) ≈ K_diag_manual

            tmp = Matrix{Float64}(undef, length(x), length(y))
            @test kernelmatrix!(tmp, k, x, y) ≈ kernelmatrix(k, x, y)

            tmp_square = Matrix{Float64}(undef, length(x), length(x))
            @test kernelmatrix!(tmp_square, k, x) ≈ kernelmatrix(k, x)

            tmp_diag = Vector{Float64}(undef, length(x))
            @test kerneldiagmatrix!(tmp_diag, k, x) ≈ kerneldiagmatrix(k, x)
        end
    end
    test_ADs(x->SqExponentialKernel() * LinearKernel(c= x[1]), rand(1), ADs = [:ForwardDiff, :ReverseDiff])
    @test_broken "Zygote issue"
end


text:
if basename(pwd()) == "aoc"
    cd("2018/23")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
"""
$(TYPEDEF)

XYZ Euclidean manifold variable node softtype.

Example
-------
```julia
p3 = Point3()
```
"""
struct Point3 <: IncrementalInference.InferenceVariable
  dims::Int
  labels::Vector{String}
  manifolds::Tuple{Symbol, Symbol, Symbol}
  Point3() = new(3, String[], (:Euclid,:Euclid,:Euclid))
end


text:
function check_yml_params(d)
  species_keys = ["migration threshold", "number of genes", "number of phenotypes", "abiotic phenotypes", "biotic phenotypes", "migration phenotype", "migration threshold", "ploidy", "epistasis matrix", "pleiotropy matrix", "growth rate", "expression array", "selection coefficient", "mutation probabilities", "mutation magnitudes", "N", "vision radius", "check fraction", "environmental noise", "optimal phenotype values", "optimal phenotypes", "age", "recombination", "initial energy", "bottleneck function", "bottleneck times", "reproduction start age", "reproduction end age"]
  model_keys = ["generations", "space", "metric", "periodic", "resources", "interactions", "food sources", "seed"]
  for (k, sp) in d["species"]
    @assert haskey(sp, "name") "Species name field missing."
    for kk in species_keys
      @assert haskey(sp, kk) "Species $(sp["name"]) does not have field $(kk)."
    end
  end
  for kk in model_keys
    @assert haskey(d["model"], kk) "Model parameter $kk is missing."
  end
end

function check_param_shapes(d)
  nspecies = length(d["species"])

  for species in 1:nspecies
    spacesize = prod(d["model"]["space"])
    dd = d["species"][species]
    # Size of optimal phenotypes
    for opt in dd["optimal phenotype values"]
      @assert length(opt) == length(dd["abiotic phenotypes"]) * spacesize "Not as many optimal phenotypes as abiotic phenotypes. Species: $species"
      # for mat in opt
      #   @assert length(mat) == spacesize "Optimal phenotpyes per trait is not the same size as the number of sites"
      # end
    end
    # Biotic and abiotic phenotypes are Array
    @assert typeof(dd["abiotic phenotypes"]) <: AbstractArray "Abiotic phenotypes should be array"
    @assert typeof(dd["biotic phenotypes"]) <: AbstractArray "Biotic phenotypes should be array"
    # Ploidy
    @assert dd["ploidy"] < 3  "Ploidy more than 2 is not implemented"
    # epistasis matrix
    # @assert size(dd["epistasis matrix"], 2) % dd["ploidy"] == 0 "Number of columns in epistasisMat are not correct. They should a factor of ploidy"
    @assert length(dd["epistasis matrix"])  == (dd["ploidy"] * dd["number of genes"])^2 "epistasis matrix does not have correct number of elements in species $species."
    # name is string
    @assert typeof(dd["name"]) <: AbstractString "name should be a string"
    # age is integer
    @assert typeof(dd["age"]) <: Int "Age of species $species should be integer."
    @assert typeof(dd["reproduction start age"]) <: Int "reproduction start age of species $species should be integer."
    @assert typeof(dd["reproduction end age"]) <: Int "reproduction end age of species $species should be integer."
    # Check fraction is float
    @assert typeof(dd["check fraction"]) <: Real "Check fraction of species $species should be number."
    # recombination is Bool
    @assert typeof(dd["recombination"]) <: Real "recombination of species $species should be type numeric"
    @assert typeof(dd["initial energy"]) <: Real "Initial energy should be a number"
    # bottleneck should be nothing or array/string
    @assert typeof(dd["bottleneck function"]) <: AbstractString || isnothing(dd["bottleneck function"])
    if !isnothing(dd["bottleneck function"])
      @assert isfile(dd["bottleneck function"]) "bottleneck function points to a nonexistant file"
    end
    @assert typeof(dd["bottleneck times"]) <: AbstractArray{Int} || isnothing(dd["bottleneck times"])
  end
  # Space should be 2D
  if !isnothing(d["model"]["space"]) && d["model"]["space"] != "nothing"
    @assert length(d["model"]["space"]) == 2 "Space should be only 2D"
  end
  # Resources
  @assert typeof(d["model"]["resources"]) <: AbstractArray{Int} "Resources should be array of Integers"
  @assert length(d["model"]["resources"]) == prod(d["model"]["space"]) "Resources should be as long as number of sites"
  # # Area
  # @assert typeof(d["model"]["area"]) <: AbstractArray "Area should be array"
  # @assert eltype(d["model"]["area"]) <: Real "Area elements should be numbers"
  # @assert length(d["model"]["area"]) == prod(d["model"]["space"]) "Area should be as long as number of sites"
  # Interactions
  @assert typeof(d["model"]["interactions"]) <: AbstractArray "Interactions should be an array"
  @assert eltype(d["model"]["interactions"]) <: AbstractFloat "Elements of interactions should be floating numbers"
  @assert length(d["model"]["interactions"]) == nspecies*nspecies "Interactions should be as many as number of species times number of species"
  # Food sources
  @assert typeof(d["model"]["food sources"]) <: AbstractArray "Food sources should be array"
  @assert eltype(d["model"]["food sources"]) <: AbstractFloat "Elements of food sources should be floating numbers"
  @assert length(d["model"]["food sources"]) == nspecies*nspecies "Food sources should be as many as number of species times number of species"
end

function reformat_params!(d)
  nspecies = length(d["species"])

  # Fix formats and shapes
  for species in 1:nspecies
    # 1. Reshape and convert pleiotropy matrix
    pmat = d["species"][species]["pleiotropy matrix"]
    ntraits = d["species"][species]["number of phenotypes"]
    ngenes = d["species"][species]["number of genes"] * d["species"][species]["ploidy"]
    d["species"][species]["pleiotropy matrix"] = Bool.(reshape(pmat, ntraits, ngenes))
    # 2. Reshape epistasis matrix
    d["species"][species]["epistasis matrix"] = reshape(d["species"][species]["epistasis matrix"], ngenes, ngenes)
    # 3. add key ngenes
    d["species"][species]["ngenes"] = ngenes
    # 4. reshape optimal phenotype values
    spacesize = prod(d["model"]["space"])
    optphens = d["species"][species]["optimal phenotype values"]
    output = [Array{Float64, 2}[] for i in  1:length(optphens)]
    for optind in 1:length(optphens)
      start = -spacesize + 1
      for phen in 1:length(d["species"][species]["abiotic phenotypes"])
        first = start + (phen * spacesize)
        last = first + spacesize - 1
        push!(output[optind], reshape(optphens[optind][first:last], Tuple(d["model"]["space"])...))
      end
    end
    d["species"][species]["optimal phenotype values"] = output
    # Check fraction to float
    d["species"][species]["check fraction"] = Float64(d["species"][species]["check fraction"])
    # include bottleneck function
    if !isnothing(d["species"][species]["bottleneck function"])
      include(d["species"][species]["bottleneck function"])
      # reshape bottleneck times
      d["species"][species]["bottleneck times"] = Bool.(reshape(d["species"][species]["bottleneck times"], prod(d["model"]["space"]), d["model"]["generations"]))
    else
      d["species"][species]["bottleneck times"] = Bool.(zeros(Int, prod(d["model"]["space"]), d["model"]["generations"]))
    end
  end

  # Metric to Symbol
  d["model"]["metric"] = Symbol(d["model"]["metric"])

  # "nothing" to nothing for seed and space
  if d["model"]["seed"] == "nothing"
    d["model"]["seed"] = nothing
  end
  if d["model"]["space"] == "nothing"
    d["model"]["space"] = nothing
  end

  # space to tuple
  if !isnothing(d["model"]["space"])
    d["model"]["space"] = Tuple(d["model"]["space"])
  end

  # reshape and convert resources
  d["model"]["resources"] = reshape(d["model"]["resources"], d["model"]["space"]...)
  # # reshape and convert area
  # d["model"]["area"] = reshape(d["model"]["area"], d["model"]["space"]...)
  # reshape and convert interactions
  d["model"]["interactions"] = reshape(d["model"]["interactions"], nspecies, nspecies)
  # reshape and convert food sources
  d["model"]["food sources"] = reshape(d["model"]["food sources"], nspecies, nspecies)
end


text:
struct PlanetaryConstants{T<:AbstractFloat} <: ConstantsCollection
    Ω::T  # Rotation rate of the planet [rad/s].
    f::T  # Nominal value for the Coriolis frequency [rad/s].
    g::T  # Standard acceleration due to gravity [m/s²].
end

PlanetaryConstants(; Ω=1.0, f=0.0, g=1.0) = PlanetaryConstants(Ω, f, g)

function choose_f(Ω, f, lat)
end

function Earth(; f=nothing, lat=nothing)
    Ω = 7.2921150e-5
    g = 9.80665

    if isnothing(f) && isnothing(lat)
        f′ = 1e-4  # Corresponds to a latitude of 43.29°N.
    elseif isnothing(lat)
        f′ = f
    elseif isnothing(f)
        f′ = 2*Ω*sind(lat)
    else
        throw(ArgumentError("Cannot specify both f and lat!"))
    end

    abs(f′) <= 2Ω || throw(ArgumentError("Coriolis parameter |f| cannot be larger than 2Ω!"))

    PlanetaryConstants(Ω, f′, g)
end

function EarthStationary()
    Ω = 0
    f = 0
    g = 9.80665
    PlanetaryConstants(Ω, f, g)
end

function Europa(; f=nothing, lat=nothing)
    # Values taken from Soderlund (Table 1, 2019) [arXiv:1901.04093].
    Ω = 2.1e-5
    g = 1.3

    if isnothing(f) && isnothing(lat)
        throw(ArgumentError("Must specify Coriolis parameter f or latitude!"))
    elseif isnothing(lat)
        f′ = f
    elseif isnothing(f)
        f′ = 2*Ω*sind(lat)
    else
        throw(ArgumentError("Cannot specify both f and lat!"))
    end

    abs(f′) <= 2Ω || throw(ArgumentError("Coriolis parameter |f| cannot be larger than 2Ω!"))

    PlanetaryConstants(Ω, f, g)
end

function Enceladus(; f=nothing, lat=nothing)
    # Values taken from Soderlund (Table 1, 2019) [arXiv:1901.04093].
    Ω = 5.3e-5
    g = 0.1

    if isnothing(f) && isnothing(lat)
        throw(ArgumentError("Must specify Coriolis parameter f or latitude!"))
    elseif isnothing(lat)
        f′ = f
    elseif isnothing(f)
        f′ = 2*Ω*sind(lat)
    else
        throw(ArgumentError("Cannot specify both f and lat!"))
    end

    abs(f′) <= 2Ω || throw(ArgumentError("Coriolis parameter |f| cannot be larger than 2Ω!"))

    PlanetaryConstants(Ω, f, g)
end


text:
using BenchmarkLite
using Iterators
using Compat

type Collect{Itr} <: Proc end

Base.string{Itr}(::Collect{Itr}) = string(Itr)
Base.length(::Collect, n::Int) = n
Base.isvalid(::Collect, n::Int) = n > 0
Base.done(::Collect, n, s) = nothing

abstract Take

type Take1 <: Take end
type Take2 <: Take end
type Take3 <: Take end
type Take4 <: Take end
type Take5 <: Take end
type Take6 <: Take end

collect(::Take, v, n) = Base.collect(take(v, n))

function Base.run{Itr<:Take}(::Collect{Itr}, n::Int, v)
    op = Itr()
    collect(op, v, n)
end

Base.start(::Collect{Take1}, n::Int) = count()
Base.start(::Collect{Take2}, n::Int) = repeated(1)
Base.start(::Collect{Take3}, n::Int) = cycle(1:3)
Base.start(::Collect{Take4}, n::Int) = cycle(1:500)
Base.start(::Collect{Take5}, n::Int) = chain(filter(iseven, 1:n), count(2,2))
Base.start(::Collect{Take6}, n::Int) = partition(cycle(1:3), 2, 1)

############

abstract Partition

type Partition1 <: Partition end
type Partition2 <: Partition end
type Partition3 <: Partition end

collect(::Partition, v, m, n) = Base.collect(partition(v, m, n))

function Base.run{Itr<:Partition}(p::Collect{Itr}, n::Int, state)
    v,x,y = state
    op = Itr()
    collect(op, v, x, y)
end

Base.start(::Collect{Partition1}, n::Int) = (1:n+1, 2, 1)
Base.start(::Collect{Partition2}, n::Int) = (1:3n+3, 2, 3)
Base.start(::Collect{Partition3}, n::Int) = (1:10n+10, 2, 10)

############

abstract GroupBy

type GroupBy1 <: GroupBy end
type GroupBy2 <: GroupBy end

collect(::GroupBy, f, v) = Base.collect(groupby(f, v))

function Base.run{Itr<:GroupBy}(p::Collect{Itr}, n::Int, state)
    f, v = state
    op = Itr()
    collect(op, f, v)
end

Base.start(::Collect{GroupBy1}, n::Int) = (x->x[1], ["abc"[[rand(1:3), rand(1:3)]] for i = 1:n])
Base.start(::Collect{GroupBy2}, n::Int) = (iseven, [1:n])

############

type IMap end

collect(::IMap, op, vs) = Base.collect(imap(op, vs...))

function Base.run(p::Collect{IMap}, n::Int, state)
    iop, vs = state
    op = IMap()
    collect(op, iop, vs)
end

Base.start(::Collect{IMap}, n::Int) = (+, (1:n, 2n:-2:0))

############

type Distinct end
collect(::Distinct, v) = Base.collect(distinct(v))

type Subsets end
collect(::Subsets, v) = Base.collect(subsets(v))

@compat function Base.run{Itr<:Union{Subsets, Distinct}}(p::Collect{Itr}, n::Int, v)
    op = Itr()
    collect(op, v)
end

Base.start(::Collect{Distinct}, n::Int) = rand(1:n, n)
Base.start(::Collect{Subsets}, n::Int) = 1:n

#############

type Chain end
collect(::Chain, vs) = Base.collect(chain(vs...))

type Product end
collect(::Product, vs) = Base.collect(product(vs...))

@compat function Base.run{Itr<:Union{Chain, Product}}(p::Collect{Itr}, n::Int, vs)
    op = Itr()
    collect(op, vs)
end

Base.start(::Collect{Chain}, n::Int) = (1:3:n, 2:3:n+1, 3:3:n+2)
Base.start(::Collect{Product}, n::Int) = (1:n>>1, 1:n>>2, 1:n>>3)

################

procs1 = Proc[
              Collect{Take1}(),
              Collect{Take2}(),
              Collect{Take3}(),
              Collect{Take4}(),
              Collect{Take5}(),
              Collect{Take6}(),
              Collect{Partition1}(),
              Collect{Partition2}(),
              Collect{Partition3}(),
              Collect{GroupBy1}(),
              Collect{GroupBy2}(),
              Collect{Distinct}(),
              Collect{Chain}(),
             ]

cfgs1 = 4 .^ (6:10)

procs2 = Proc[
              Collect{Subsets}(),
              ]

cfgs2 = [4:2:16]

procs3 = Proc[
              Collect{Product}(),
              ]

cfgs3 = 2 .^ (5:7)

###########

rtable1 = run(procs1, cfgs1)
rtable2 = run(procs2, cfgs2)
rtable3 = run(procs3, cfgs3)

show(rtable1; unit=:msec)
println()
show(rtable2; unit=:msec)
println()
show(rtable3; unit=:msec)
println()


text:
function traverse(posx::Float64,posy::Float64,velx::Float64,vely::Float64,
                    x0::Float64,x1::Float64,y0::Float64,y1::Float64)
    # Given a particle in a cell with position (posx,posy)
    # and velocity (velx,vely) where the cell is defined by
    # [x0,x1] x [y0,y1],
    # this function returns nexti, nextj = the direction of the cell
    # that the particle will, as well as nextx,nexty = the intersection point
    # and dist = the travelled distance.
    # move into after having passed the current cell.
    #
    # The layout is as follows:
    #
    #                     j-1            j            j+1
    #     \  y
    #      \
    #     x \
    #
    #   i+1
    #
    #                 ---------------x1---------------
    #                 |                              |
    #                 |                              |
    #     i           y0                             y1
    #                 |                              |
    #                 |                              |
    #                 ---------------x0---------------
    #
    #   i-1
    #
    # The coordinate system is
    #
    #  ^  x
    #  |
    #  |
    #  |
    #    --------> y
    #

    # assert that we are in the domain
    #eps = 1e-14
    #@assert posx>= x0-eps
    #@assert posx<= x1+eps
    #@assert posy>= y0-eps
    #@assert posy<= y1+eps

    # assert a non zero velocity
    normv = sqrt(vely^2+velx^2)
    #@assert normv>0

    # are we moving north east, south east, south west or north west
    vertical   = ifelse(velx>=0,+1,-1)
    horizontal = ifelse(vely>=0,+1,-1)

    # which face (left right top bottom) will we hit
    # i.e. to which y or x do we need to calculate the distance
    x = ifelse(vertical==+1,x1,x0)
    y = ifelse(horizontal==+1,y1,y0)

    # calculate distance to that point considering the velocity
    distx = abs((x-posx)/velx)
    disty = abs((y-posy)/vely)

    # check which face we hit first and save the corresponding i/j shift
    # as well as the distance travelled
    if disty<distx
        nexti = 0
        nextj = horizontal
        dist = disty
    elseif disty>distx
        nexti = vertical
        nextj = 0
        dist = distx
    else
        nexti = vertical
        nextj = horizontal
        dist = distx # == disty
    end

    # at what point to we cross the interface
    nextx = posx + dist*velx
    nexty = posy + dist*vely
    return nexti,nextj, dist*normv, nextx,nexty,velx,vely
end


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

### Standard map display "layers." ###
const LAYER_STANDARD = Compat.@Dict(
    1 => Style(0x7BB6EF, 4), # Soft blue
    2 => Style(0x66C266, 3), # Soft green
    3 => Style(0xE68080, 3), # Soft red
    4 => Style(0xFF9900, 3), # Soft orange
    5 => Style(0xDADA47, 3), # Dark yellow
    6 => Style(0x999999, 2), # Dark gray
    7 => Style(0xE0E0E0, 2), # Light gray
    8 => Style(0x999999, 1)) # Dark gray

const LAYER_CYCLE = Compat.@Dict(
    1 => Style(0x006600, 3), # Green
    2 => Style(0x5C85FF, 3), # Blue
    3 => Style(0x5C85FF, 2), # Blue
    4 => Style(0x999999, 2)) # Dark gray

const LAYER_PED = Compat.@Dict(
    1 => Style(0x999999, 3), # Dark gray
    2 => Style(0x999999, 3), # Dark gray
    3 => Style(0x999999, 2), # Dark gray
    4 => Style(0xE0E0E0, 2)) # Light gray

const LAYER_FEATURES = Compat.@Dict(
    1 => Style(0x9966FF, 1.5, "."), # Lavender
    2 => Style(0xFF0000, 1.5, "."), # Red
    3 => Style(0x000000, 1.5, "."), # Black
    4 => Style(0xFF66FF, 1.5, "."), # Pink
    5 => Style(0x996633, 1.5, "."), # Brown
    6 => Style(0xFF9900, 2.0, "."), # Orange
    7 => Style(0xCC00CC, 1.5, ".")) # Brown

const LAYER_BUILDINGS = Compat.@Dict(
    1 => Style(0xE1E1EB, 1, "-"), # Lighter gray
    2 => Style(0xB8DBFF, 1, "-"), # Light blue
    3 => Style(0xB5B5CE, 1, "-"), # Light gray
    4 => Style(0xFFFF99, 1, "-"), # Pale yellow
    5 => Style(0x006600, 1, "-")) # Green


text:
using Ipopt
using SparseArrays
function solve_ipopt(P::AbstractMatrix, q, A::AbstractMatrix, b,
    r_max=Inf, r_min=-Inf, x_init=zeros(0);
    max_iter=5000, print_level=0)

    n = size(P, 1)
    A_ = [A; ones(n)']
    m = size(A_, 1)

    eval_f(x) = dot(x, Symmetric(P)*x)/2 + dot(q, x)
    eval_grad_f(x, grad_f) = grad_f[:] = P*x + q

    function eval_g(x, g)
        g[1:end-1] = A*x - b
        g[end] = dot(x, x)
    end

    function eval_jac_g(x, mode, rows, cols, values)
        if mode == :Structure
            k = 1
            for j = 1:n, i=1:m rows[k] = i
                cols[k] = j
                k += 1
            end
        else
            A_[end, :] .= 2*x
            values[:] .= reshape(A_, length(A_))
        end
    end

    function eval_h(x, mode, rows, cols, obj_factor, lambda, values)
        if mode == :Structure
            k = 1
            @inbounds for j = 1:n, i=j:n
                rows[k] = i; cols[k] = j
                k += 1
            end
        else
            k = 1
            @inbounds for j=1:n, i=j:n
                values[k] = obj_factor*P[j, i]
                if i == j
                    values[k] += 2*lambda[end]
                end
                k +=1
            end
        end
    end

    x_L = -Inf*ones(n)
    x_U = Inf*ones(n)

    g_L = -Inf*ones(m)
    g_L[end] = r_min^2
    g_U = zeros(m)
    g_U[end] = r_max^2

    nz_P = Int(n*(n + 1)/2)
    nz_A_ = m*n

    prob = createProblem(n, x_L, x_U, m, g_L, g_U, nz_A_, nz_P,
                eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
    addOption(prob, "max_iter", max_iter)
    addOption(prob, "print_level", print_level)
    if length(x_init) == n
        prob.x .= x_init
    end

    timing = @elapsed status = solveProblem(prob)
    # println(Ipopt.ApplicationReturnStatus[status])
    # println(prob.obj_val)
    prob.mult_g[end] *= 2
    return prob.x, prob.mult_g, timing
end

function solve_ipopt(P::SparseMatrixCSC, q, A::SparseMatrixCSC, b,
    r_max=Inf, r_min=-Inf, x_init=zeros(0);
    max_iter=5000, print_level=0)

    # Ensure diagonal is present
    n = size(P, 1)
    P = P + 1e-60*I
    P_lower = SparseMatrixCSC(LowerTriangular(P))
    A_ = [A; ones(1, n)]
    m, n = size(A_)

    eval_f(x) = dot(x, P*x)/2 + dot(q, x)
    eval_grad_f(x, grad_f) = grad_f[:] = P*x + q

    function eval_g(x, g)
        # Bad: g    = zeros(2)  # Allocates new array
        # OK:  g[:] = zeros(2)  # Modifies 'in place'
        g[1:end-1] = A*x - b
        g[end] = dot(x, x)
    end

    function eval_jac_g(x, mode, rows, cols, values)
        if mode == :Structure
            rows[:], cols[:] = findnz(A_)
        else
            m, n = size(A_)
            @inbounds for i = 1:n
                A_[end, i] = 2*x[i]
            end
            values[:] .= A_.nzval
        end
    end

    function eval_h(x, mode, rows, cols, obj_factor, lambda, values)
        if mode == :Structure
            rows[:], cols[:] = findnz(P_lower)
        else
            P_ = obj_factor*P_lower + 2*lambda[end]*I
            values[:] .= P_.nzval
        end
    end

    x_L = -Inf*ones(n)
    x_U = Inf*ones(n)

    g_L = -Inf*ones(m)
    g_L[end] = r_min^2
    g_U = zeros(m)
    g_U[end] = r_max^2

    nz_P = nnz(P_lower)
    nz_A_ = nnz(A_)

    prob = createProblem(n, x_L, x_U, m, g_L, g_U, nz_A_, nz_P,
                eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
    addOption(prob, "max_iter", max_iter)
    addOption(prob, "print_level", print_level)
    if length(x_init) == n
        prob.x .= x_init
    end
    timing = @elapsed status = solveProblem(prob)

    # println(Ipopt.ApplicationReturnStatus[status])
    # println(prob.obj_val)
    # println(prob.mult_g)
    prob.mult_g[end] *= 2
    return prob.x, prob.mult_g, timing
end


text:
#!/usr/bin/env julia

using DataFrames
import JLD

import Celeste.AccuracyBenchmark
import Celeste.ArgumentParse
import Celeste.DeterministicVI
import Celeste.Log
import Celeste.Model
import Celeste.ParallelRun
import Celeste.SDSSIO
import Celeste.Transform

const TEST_DATA_DIR = joinpath(Pkg.dir("Celeste"), "test", "data")
const FIELD_EXTENTS_FITS = joinpath(Pkg.dir("Celeste", "test", "data", "field_extents.fits"))
const CELESTE_CATALOG_DIR = "/global/cscratch1/sd/jregier/celeste_sdds_catalog_jld"
const CATALOG_FILE_REGEXP = r"^celeste-(-?[0-9\.]+)-(-?[0-9\.]+)-(-?[0-9\.]+)-(-?[0-9\.]+)\.jld$"

parser = ArgumentParse.ArgumentParser()
ArgumentParse.add_argument(
    parser,
    "--catalog-dir",
    help="Path to Celeste catalog files (celeste-<ramin>-<ramax>-<decmin>-<decmax>.jld)",
    default=CELESTE_CATALOG_DIR,
)
ArgumentParse.add_argument(
    parser,
    "--subtract-sdss",
    action=:store_true,
    help="Write Celeste expectation - SDSS imagery instead of just expectation",
)
ArgumentParse.add_argument(
    parser,
    "run",
    help="SDSS run #",
    arg_type=Int,
)
ArgumentParse.add_argument(
    parser,
    "camcol",
    help="SDSS camcol #",
    arg_type=Int,
)
ArgumentParse.add_argument(
    parser,
    "field",
    help="SDSS field #",
    arg_type=Int,
)
parsed_args = ArgumentParse.parse_args(parser, ARGS)

function find_relevant_catalog_files(
    ra_min_deg::Float64, ra_max_deg::Float64, dec_min_deg::Float64, dec_max_deg::Float64
)
    @assert ra_max_deg >= ra_min_deg
    @assert dec_max_deg >= dec_min_deg

    relevant_paths = String[]
    for filename in readdir(parsed_args["catalog-dir"])
        match_values = match(CATALOG_FILE_REGEXP, filename)
        if match_values === nothing
            continue
        else
            ra_start_deg = parse(Float64, match_values.captures[1])
            ra_end_deg = parse(Float64, match_values.captures[2])
            dec_start_deg = parse(Float64, match_values.captures[3])
            dec_end_deg = parse(Float64, match_values.captures[4])
            if ra_start_deg > ra_max_deg || ra_end_deg < ra_min_deg
                continue
            elseif dec_start_deg > dec_max_deg || dec_end_deg < dec_min_deg
                continue
            else
                push!(relevant_paths, joinpath(parsed_args["catalog-dir"], filename))
            end
        end
    end
    relevant_paths
end

function load_sources(
    ra_min_deg::Float64, ra_max_deg::Float64, dec_min_deg::Float64, dec_max_deg::Float64
)
    file_paths = find_relevant_catalog_files(ra_min_deg, ra_max_deg, dec_min_deg, dec_max_deg)
    Log.info("Found $(length(file_paths)) relevant catalog files")
    all_sources = ParallelRun.OptimizedSource[]
    for file_path in file_paths
        Log.info("Loading $file_path")
        sources::Vector{ParallelRun.OptimizedSource} = JLD.load(file_path)["results"]
        Log.info("  Found $(length(sources)) sources in $file_path")
        for source in sources
            source_ra_deg = source.vs[Model.ids.pos[1]]
            source_dec_deg = source.vs[Model.ids.pos[2]]
            if source_ra_deg < ra_min_deg || source_ra_deg > ra_max_deg
                continue
            elseif source_dec_deg < dec_min_deg || source_dec_deg > dec_max_deg
                continue
            else
                push!(all_sources, source)
            end
        end
    end
    all_sources
end

function catalog_entry_from_variational_params(variational_parameters)
    data = AccuracyBenchmark.variational_parameters_to_data_frame_row("0", variational_parameters)
    AccuracyBenchmark.make_catalog_entry(first(eachrow(data)))
end

# based on code in DeterministicVI.elbo_likelihood()
function fill_celeste_expectation!(
    images::Vector{Model.Image}, sources::Vector{ParallelRun.OptimizedSource}
)
    catalog = [catalog_entry_from_variational_params(source.vs) for source in sources]
    patches = Model.get_sky_patches(images, catalog, radius_override_pix=25.0)
    variational_params = Vector{Float64}[source.vs for source in sources]

    active_sources = [1]
    elbo_args = DeterministicVI.ElboArgs(images, patches, active_sources)
    elbo_vars = DeterministicVI.ElboIntermediateVariables(Float64, elbo_args.Sa, false, false)
    bvn_bundle = Model.BvnBundle{Float64}(elbo_args.psf_K, elbo_args.S)
    source_brightnesses = DeterministicVI.load_source_brightnesses(
        elbo_args,
        variational_params,
        calculate_gradient=false,
        calculate_hessian=false,
    )
    for image_index in 1:elbo_args.N
        Model.load_bvn_mixtures!(
            bvn_bundle.star_mcs,
            bvn_bundle.gal_mcs,
            elbo_args.S,
            elbo_args.patches,
            variational_params,
            elbo_args.active_sources,
            elbo_args.psf_K,
            image_index,
            elbo_vars.elbo.has_gradient,
            elbo_vars.elbo.has_hessian,
        )
        image = images[image_index]
        for h in 1:image.H, w in 1:image.W
            # fills elbo_vars.E_G in-place
            DeterministicVI.add_pixel_term!(
                elbo_args,
                variational_params,
                image_index, h, w,
                bvn_bundle,
                source_brightnesses,
                elbo_vars,
            )
            image.pixels[h, w] += elbo_vars.E_G.v[] - image.sky[h, w]
        end
    end
end

function subtract_sdss_values!(
    celeste_images::Vector{Model.Image}, sdss_images::Vector{Model.Image}
)
    for band in 1:5, h in 1:celeste_images[band].H, w in 1:celeste_images[band].W
        celeste_images[band].pixels[h, w] -= sdss_images[band].pixels[h, w]
    end
end

function main()
    rcf = SDSSIO.RunCamcolField(
        parsed_args["run"],
        parsed_args["camcol"],
        parsed_args["field"],
    )

    working_directory = pwd()
    cd(TEST_DATA_DIR)
    command = `make RUN=$(rcf.run) CAMCOL=$(rcf.camcol) FIELD=$(rcf.field)`
    Log.info("Downloading SDSS frame files ($command in $TEST_DATA_DIR)...")
    run(command)
    cd(working_directory)

    Log.info("Opening field extents: $FIELD_EXTENTS_FITS")
    field_extents_data = AccuracyBenchmark.load_fits_table_as_data_frame(
        FIELD_EXTENTS_FITS,
        2,
        [:run, :camcol, :field, :ramin, :ramax, :decmin, :decmax],
    )
    row_selector = (
        (field_extents_data[:run] .== rcf.run)
        .& (field_extents_data[:camcol] .== rcf.camcol)
        .& (field_extents_data[:field] .== rcf.field)
    )
    field_row = field_extents_data[row_selector, :]
    Log.info(field_row)

    full_sources = load_sources(
        field_row[1, :ramin],
        field_row[1, :ramax],
        field_row[1, :decmin],
        field_row[1, :decmax],
    )
    Log.info("Found $(length(full_sources)) sources total")

    strategy = SDSSIO.PlainFITSStrategy(AccuracyBenchmark.SDSS_DATA_DIR)
    sdss_images = SDSSIO.load_field_images(strategy, [rcf])
    images = deepcopy(sdss_images)

    for band in 1:5, h in 1:images[band].H, w in 1:images[band].W
        images[band].pixels[h, w] = 0
    end

    fill_celeste_expectation!(images, full_sources)

    label = "expectations"
    if parsed_args["subtract-sdss"]
        subtract_sdss_values!(images, sdss_images)
        label = "errors"
    end

    output_filename = "celeste_$(label)_$(rcf.run)_$(rcf.camcol)_$(rcf.field).fits"
    AccuracyBenchmark.save_images_to_fits(output_filename, images)
end

main()


text:
module FaultTree

include("_faulttree.jl")
include("_operations.jl")
include("_dot.jl")
include("_tobdd.jl")
include("_comp.jl")

include("_mcs.jl")

end


text:
const Optional{T} = Union{Nothing, T}

struct AxisAspect
    aspect::Float32
end

struct DataAspect end

"""
LinearTicks with ideally a number of `n_ideal` tick marks.
"""
struct LinearTicks
    n_ideal::Int

    function LinearTicks(n_ideal)
        if n_ideal <= 0
            error("Ideal number of ticks can't be smaller than 0, but is $n_ideal")
        end
        new(n_ideal)
    end
end

struct WilkinsonTicks
    k_ideal::Int
    k_min::Int
    k_max::Int
    Q::Vector{Tuple{Float64, Float64}}
    granularity_weight::Float64
    simplicity_weight::Float64
    coverage_weight::Float64
    niceness_weight::Float64
    min_px_dist::Float64
end

"""
Like LinearTicks but for multiples of `multiple`.
Example where approximately 5 numbers should be found
that are multiples of pi, printed like "1π", "2π", etc.:

```
MultiplesTicks(5, pi, "π")
```
"""
struct MultiplesTicks
    n_ideal::Int
    multiple::Float64
    suffix::String
end

mutable struct LineAxis
    parent::Scene
    protrusion::Node{Float32}
    attributes::Attributes
    elements::Dict{Symbol, Any}
    tickpositions::Node{Vector{Point2f0}}
    tickvalues::Node{Vector{Float32}}
    ticklabels::Node{Vector{String}}
end

struct LimitReset end

mutable struct RectangleZoom
    active::Bool
    restrict_x::Bool
    restrict_y::Bool
    from::Union{Nothing, Point2f0}
    to::Union{Nothing, Point2f0}
    rectnode::Observable{FRect2D}
    plots::Vector{AbstractPlot}
end

struct ScrollZoom
    speed::Float32
    reset_timer::Ref{Any}
    prev_xticklabelspace::Ref{Any}
    prev_yticklabelspace::Ref{Any}
    reset_delay::Float32
end

struct DragPan
    reset_timer::Ref{Any}
    prev_xticklabelspace::Ref{Any}
    prev_yticklabelspace::Ref{Any}
    reset_delay::Float32
end

struct ScrollEvent
    x::Float32
    y::Float32
end

struct KeysEvent
    keys::Set{AbstractPlotting.Keyboard.Button}
end

@Layoutable Axis begin
    scene::Scene
    xaxislinks::Vector{Axis}
    yaxislinks::Vector{Axis}
    limits::Node{FRect2D}
    block_limit_linking::Node{Bool}
    mouseeventhandle::MouseEventHandle
    scrollevents::Observable{ScrollEvent}
    keysevents::Observable{KeysEvent}
    interactions::Dict{Symbol, Tuple{Bool, Any}}
end

@Layoutable Colorbar

@Layoutable Label

@Layoutable Box

@Layoutable Slider

@Layoutable Button

@Layoutable Toggle

@Layoutable Menu

abstract type LegendElement end

struct LineElement <: LegendElement
    attributes::Attributes
end

struct MarkerElement <: LegendElement
    attributes::Attributes
end

struct PolyElement <: LegendElement
    attributes::Attributes
end

struct LegendEntry
    elements::Vector{LegendElement}
    attributes::Attributes
end

const EntryGroup = Tuple{Optional{String}, Vector{LegendEntry}}

@Layoutable Legend begin
    entrygroups::Node{Vector{EntryGroup}}
end

@Layoutable LScene begin
    scene::Scene
end

@Layoutable Textbox begin
    cursorindex::Node{Int}
    cursoranimtask
end


text:
#function for fitting a linear regression. Takes as input
#  tds::Array{ResponseData}      Training data set
#  fm::FeatureMap                Feature map  as defined in ML_models.jl
#  options                       Contains the regularization parameters and other optional arguments
# NOTE: You can not use a kernal combined with an L1 regularizor since we can no longer use the "kernel trick". But we can still use  feature maps.
function findfolder(strname)
  strs= ["/../../$strname","/../$strname","/$strname" ]
  strfolder = filter(isdir, map(str -> string(pwd(),str),strs))
  if length(strfolder) == 0
    return ""
  elseif length(strfolder) > 1
    warn("found more than on $strname folder")
  end
  return strfolder[1]
end

function findfiles(strdata)
  strbase = basename(strdata)
  strdir = dirname(strdata)
  strdir =  (strdir=="") ? "" : string("/",strdir)
  # will look for data in three folders:
  folders =[string(pwd(),"/..",strdir),string(pwd(),strdir)]
  if strdir != "/data" && strdir != "/Data"
    push!(folders, findfolder("data"))
    push!(folders, findfolder("Data"))
  end
  filter!(f -> isdir(f), folders)
  files = vcat(map(str -> map(s -> string(str,"/",s), readdir(str)), folders)...)
  filter!(f -> contains(f, strbase), files)
end

function jlsTojld(strdata)
  if ((str =split(strdata,'.')[end]) != "jls")
    error("can only convert jls, not type $str, to jld")
  end
  rds = load_rds(strdata)
  rd_pairs = [
    ("(volfrac,a)=($(rd.volfrac),$(rd.particles[1].r))",rd)
  for rd in rds]
  strfile=string(split(strdata,'.')[1],".jld")
  save(strfile, Dict(rd_pairs))
  return Dict(rd_pairs)
end

function load_moments(strfile)
  if !isfile(strfile)
    strfile =  filter(s -> split(s,'.')[length(split(s,'.'))-1]=="mnts" ,findfiles(strfile))[1]
  end
  println("loading:", strfile)
  strs = split(strfile,'.')
  if strs[end] == "jls"
    moments = open(strfile, "r") do f deserialize(f) end
  elseif strs[end] == "jld"
    moments = Array{StatisticalMoments{Float64}}(collect(values(load(strfile)))[1])
  else error("file type not recognised")
  end
  return moments
end