text:
export default_type, default_library, similar_library, library
export getlibrary, getlibraryfor

"""
    default_type(::FullDim{N}, ::Type{T}) where {N, T}

Returns the default polyhedron type for `N`-dimensional polyhedron of coefficient type `T`.
"""
function default_type end

default_type(::FullDim{N}, ::Type{T}) where {N, T} = SimplePolyhedron{N, T, Intersection{N, T, Vector{T}}, Hull{N, T, Vector{T}}}
default_type(::FullDim{1}, ::Type{T}) where T = Interval{T, SVector{1, T}}

"""
    default_library(::FullDim{N}, ::Type{T}) where {N, T}

Returns the default polyhedral library for `N`-dimensional polyhedron of coefficient type `T`.
"""
function default_library end

_default_type(::Type{T}) where T = T
# See https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/35
_default_type(::Type{AbstractFloat}) = Float64
default_library(::FullDim, ::Type{T}) where T = SimplePolyhedraLibrary{_default_type(T)}()
default_library(::FullDim{1}, ::Type{T}) where T = IntervalLibrary{_default_type(T)}()

"""
    similar_library(lib::PolyhedraLibrary, d::FullDim{N}, ::Type{T}) where {N, T}

Returns a library that supports polyhedra of full dimension `T` with coefficient type `T`. If `lib` does not support it, this commonly calls `default_library(d, T)`.
"""
function similar_library end

# Shortcuts
similar_library(p::Polyhedron{N}, d::FullDim, ::Type{T}) where {N, T} = similar_library(library(p), d, T)
similar_library(p::Polyhedron{N}, ::Type{T}) where {N, T} = similar_library(p, FullDim{N}(), T)
similar_library(p::Polyhedron{N, T}, d::FullDim) where {N, T} = similar_library(p, d, T)

"""
    library(p::Polyhedron)

Returns the library used by `p`.
"""
function library end

function getlibrary(args...)
    Base.depwarn("getlibrary is deprecated, use library instead", :getlibrary)
    library(args...)
end

function getlibraryfor(args...)
    Base.depwarn("getlibraryfor is deprecated, use similar_library instead. Note that the dimension `N` now needs to be given as `FullDim{N}()`.", :getlibraryfor)
    similar_library(args...)
end

"""
    default_solver(p::Rep)

Returns a default linear programming solver for the polyhedron `p` (e.g. CDD has an internal solver which is used by default).
"""
default_solver(p::Rep) = nothing
function default_solver(p::Rep, ps::Rep...)
    s = default_solver(p)
    if s === nothing
        default_solver(ps...)
    else
        s
    end
end

"""
    solver(p::Rep, solver::MathProgBase.AbstractMathProgSolver=default_solver(p))

If the V-representation of `p` has been computed, returns `VRepSolver()`, otherwise, returns `solver`.
"""
function solver(p::Rep, solver::Union{Nothing, MPB.AbstractMathProgSolver}=default_solver(p))
    if vrepiscomputed(p)
        VRepSolver()
    else
        solver
    end
end

solver(v::VRepresentation, solver::Union{Nothing, MPB.AbstractMathProgSolver}=default_solver(v)) = VRepSolver()
solver(h::HRepresentation, solver::Union{Nothing, MPB.AbstractMathProgSolver}=default_solver(h)) = solver

_promote_reptype(P1::Type{<:HRep}, ::Type{<:HRep}) = P1
_promote_reptype(P1::Type{<:VRep}, ::Type{<:VRep}) = P1
function promote_reptype(P1::Type{<:Rep}, P2::Type{<:Rep})
    _promote_reptype(P1, P2)
end
promote_reptype(P1::Type{<:Rep}, P2::Type{<:Rep}, P::Type{<:Rep}...) = promote_reptype(promote_reptype(P1, P2), P...)
promote_reptype(P::Type{<:Rep}) = P

# whether Rep has the constructor Rep(::It...; solver=...)
supportssolver(::Type{<:Rep}) = false

function constructpolyhedron(RepT::Type{<:Rep{N, T}}, p::Tuple{Vararg{Rep}}, it::It{N, T}...) where {N, T}
    if supportssolver(RepT)
        solver = default_solver(p...)
        if solver !== nothing
            return RepT(it..., solver=solver)
        end
    end
    RepT(it...)::RepT # FIXME without this type annotation even convexhull(::PointsHull{2,Int64,Array{Int64,1}}, ::PointsHull{2,Int64,Array{Int64,1}}) is not type stable, why ?
end

function default_similar(p::Tuple{Vararg{Rep}}, d::FullDim{N}, ::Type{T}, it::It{N, T}...) where {N, T}
    # Some types in p may not support `d` or `T` so we call `similar_type` after `promote_reptype`
    RepT = similar_type(promote_reptype(typeof.(p)...), d, T)
    constructpolyhedron(RepT, p, it...)
end

"""
    similar(p::Tuple{Vararg{Polyhedra.Rep}}, ::Polyhedra.FullDim{N}, ::Type{T}, it::Polyhedra.It{N, T}...)

Creates a representation with a type similar to `p` of a polyhedron of full dimension `N`, element type `T` and initialize it with the iterators `it`.
The type of the result will be chosen closer to the type of `p[1]`.
"""
Base.similar(p::Tuple{Vararg{Rep}}, d::FullDim{N}, ::Type{T}, it::It{N, T}...) where {N, T} = default_similar(p, d, T, it...)
function promote_coefficienttype(p::Tuple{Vararg{Rep}})
    promote_type(MultivariatePolynomials.coefficienttype.(p)...)
end
function Base.similar(p::Tuple{Vararg{Rep{N}}}, d::FullDim{N}, it::It{N}...) where N
    T = promote_coefficienttype(p)
    similar(p, d, T, it...)
end
function Base.similar(p::Tuple{Vararg{Rep{N}}}, it::It{N}...) where N
    similar(p, FullDim{N}(), it...)
end
Base.similar(p::Rep, args...) = similar((p,), args...)


text:
export funCompR2

"""
GenSPGL (this function is used for parallel_weighted spgl1, it is a traditional spgl2 function)

Use:    f,g1,g2 = funCompositeR(A, r, funForward, funPenalty, params)
"""
function funCompR2(A::TA,
                       x::AbstractArray,
                       r::AbstractArray,
                       funForward::Function, funPenalty::Function,
                       nProdAt::Int64,
                       params::Dict{String,Any}) where
                         {TA<:Union{joAbstractLinearOperator,AbstractArray}}

    nProdAt += one(Int64)
    f,v = funPenalty(r, params)

    if ~(params["proxy"])
        g1 = funForward(A, x, -v, params)
        g2 = g1[:]
    else
        g1 = funForward(A, x, -v, params)
	g2 = g1[:]
    end
    return f,g1,g2
end

ext: jl


text:
export @class

function allocclass(name, super)
  ptr = ccall(:objc_allocateClassPair, Ptr{Void}, (Ptr{Void}, Ptr{Cchar}, Csize_t),
              super, name, 0)
  ptr == C_NULL && error("Couldn't allocate class $name")
  return Class(ptr)
end

function register(class::Class)
  ccall(:objc_registerClassPair, Void, (Ptr{Void},),
        class)
  return class
end

createclass(name, super) = allocclass(name, super) |> register

getmethod(class::Class, sel::Selector) =
  ccall(:class_getInstanceMethod, Ptr{Void}, (Ptr{Void}, Ptr{Void}),
        class, sel)

methodtypeenc(method::Ptr) =
  ccall(:method_getTypeEncoding, Ptr{Cchar}, (Ptr{Void},),
        method) |> bytestring

methodtypeenc(class::Class, sel::Selector) = methodtypeenc(getmethod(class, sel))

methodtype(args...) = methodtypeenc(args...) |> parseencoding

replacemethod(class::Class, sel::Selector, imp::Ptr{Void}, types::String) =
  ccall(:class_replaceMethod, Bool, (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Cchar}),
        class, sel, imp, types)

function setmethod(class::Class, sel::Selector, imp::Ptr{Void}, types::String)
  meth = getmethod(class, sel)
  meth ≠ C_NULL && methodtype(meth) != parseencoding(types) &&
    error("New method $(name(sel)) of $class must match $(methodtype(meth))")
  replacemethod(class, sel, imp, types)
end

# Syntax

function createdef(f, args, Ts, body, ret)
  def = :(function $f($(map((x,T)->:($x::$(ctype(T))), args, Ts)...)) end)
  for (arg, T) in zip(args, Ts)
    T in (Object, Selector) && push!(def.args[2].args, :($arg = $T($arg)))
  end
  push!(def.args[2].args, ret == Void ? # Ugly workaround for 0.3
                            :($body; nothing) :
                            :(convert($(ctype(ret)), $body)))
  return def
end

function createmethod(class, ex)
  instance = ex.args[1] == symbol("@-")
  ret = eval(current_module(), ex.args[2])
  params = ex.args[3:end-1]
  sel = @_ params map(x->x.args[1], _) join(_, ":") string(_, ":") Selector
  Ts = @>> params map(x->x.args[2]) map(x->isa(x, Symbol) ? Object : x.args[2]) map(T->eval(current_module(), T))
  args = @>> params map(x->x.args[2]) map(x->isa(x, Symbol) ? x : x.args[3])
  body = ex.args[end]

  Ts = [Object, Selector, Ts...]
  args = [:self, :_cmd, args...]
  typ = encodetype(ret, Ts...)

  f = "$(class)_$(name(sel))" |> symbol
  return quote
    $(esc(createdef(f, args, Ts, body, ret)))
    setmethod($(instance ? esc(class) : :(class($(esc(class))))), $sel,
              cfunction($(esc(f)), $(ctype(ret)), $(Expr(:tuple, map(ctype, Ts)...))),
              $typ)
  end
end

macro class (def)
  name = namify(def.args[2])
  sym = Expr(:quote, name)
  super = isexpr(def.args[2], :(<:)) ? def.args[2].args[2] : :NSObject
  expr = quote
    classexists($sym) || createclass($sym, Class($(Expr(:quote, super))))
    isdefined($sym) || (const $(esc(name)) = Class($sym))
  end
  for ex in def.args[3].args
    if isexpr(ex, :macrocall) && ex.args[1] in (symbol("@+"), symbol("@-"))
      ex = createmethod(name, ex)
    end
    push!(expr.args, ex)
  end
  :(@objc $expr; $(esc(name)))
end

# quote
#   @class type Foo
#     @- (Cdouble) multiply:(Cdouble)x by:(Cdouble)y begin
#       x*y
#     end
#   end
# end |> macroexpand

# @class type Foo
#   @- (Cdouble) multiply:(Cdouble)x by:(Cdouble)y begin
#     x*y
#   end
# end

# @objc begin
#   foo = [Foo new]
#   x = [foo multiply:28 by:26]
#   [foo release]
#   x
# end


text:
# time series

# Makie
using CairoMakie, TimeSeries, Dates
let
    # dummy data
    dates = Date(2018, 1, 1):Day(1):Date(2018, 12, 31)
    ta = TimeArray(dates, rand(length(dates)))
    vals = 0.5*values(ta)
    tempo = string.(timestamp(ta))
    lentime = length(tempo)
    slice_dates = range(1,lentime, step = lentime÷8)

    fig = Figure(resolution = (600, 400), backgroundcolor = :honeydew, font =:sans)
    ax = Axis(fig[1,1], xlabel = "Date", ylabel = "value")
    line1 = lines!(ax, 1:lentime, vals; color = :black, linewidth = 0.85)
    ax.xticks = (slice_dates, tempo[slice_dates])
    ax.xticklabelrotation = π/4
    ax.xticklabelalign = (:right, :center)
    display(fig)
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# Used "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy as
# a reference.  Note that there were some typos in that document so the code
# here may not correspond exactly.

immutable NormalWishart <: Distribution
    dim::Int
    zeromean::Bool
    mu::Vector{Float64}
    kappa::Float64
    Tchol::Cholesky{Float64}  # Precision matrix (well, sqrt of one)
    nu::Float64

    function NormalWishart(mu::Vector{Float64}, kappa::Real,
                                  Tchol::Cholesky{Float64}, nu::Real)
        # Probably should put some error checking in here
        d = length(mu)
        zmean::Bool = true
        for i = 1:d
            if mu[i] != 0.
                zmean = false
                break
            end
        end
        new(d, zmean, mu, float64(kappa), Tchol, float64(nu))
    end
end

function NormalWishart(mu::Vector{Float64}, kappa::Real,
                       T::Matrix{Float64}, nu::Real)
    NormalWishart(mu, kappa, cholfact(T), nu)

end

function insupport(::Type{NormalWishart}, x::Vector{Float64}, Lam::Matrix{Float64})
    return (all(isfinite(x)) &&
           size(Lam, 1) == size(Lam, 2) &&
           isApproxSymmmetric(Lam) &&
           size(Lam, 1) == length(x) &&
           hasCholesky(Lam))
end

pdf(nw::NormalWishart, x::Vector{Float64}, Lam::Matrix{Float64}) =
        exp(logpdf(nw, x, Lam))

function logpdf(nw::NormalWishart, x::Vector{Float64}, Lam::Matrix{Float64})
    if !insupport(NormalWishart, x, Lam)
        return -Inf
    else
        p = length(x)

        nu = nw.nu
        kappa = nw.kappa
        mu = nw.mu
        Tchol = nw.Tchol
        hnu = 0.5 * nu
        hp = 0.5 * p

        # Normalization
        logp::Float64 = hp*(log(kappa) - float64(log2π))
        logp -= hnu * logdet(Tchol)
        logp -= hnu * p * log(2.)
        logp -= lpgamma(p, hnu)

        # Wishart (MvNormal contributes 0.5 as well)
        logp += (hnu - hp) * logdet(Lam)
        logp -= 0.5 * trace(Tchol \ Lam)

        # Normal
        z = nw.zeromean ? x : x - mu
        logp -= 0.5 * kappa * dot(z, Lam * z)

        return logp

    end
end

function rand(nw::NormalWishart)
    Lam = rand(Wishart(nw.nu, nw.Tchol))
    mu = rand(MvNormal(nw.mu, inv(Lam) ./ nw.kappa))
    return (mu, Lam)
end


text:
module VisualOptim
include("update.jl")
include("lossfunction.jl")
include("optimizer.jl")
include("utils.jl")
end


text:
# Source: http://salvi.chaosnet.org/snippets/i-patch.html
"""
write_volume(voxels, filename)

Writes the 3D matrix `voxels` into `filename.raw`
with the metaimage header `filename.mhd`.
"""
function write_volume(voxels, filename)
    n = size(voxels)
    open("$(filename).mhd", "w") do f
        println(f, "NDims = 3")
        println(f, "DimSize = $(n[1]) $(n[2]) $(n[3])")
        println(f, "ElementSize = 1.0 1.0 1.0")
        println(f, "ElementSpacing = 1.0 1.0 1.0")
        println(f, "ElementType = MET_DOUBLE")
        println(f, "ElementByteOrderMSB = False") # LittleEndian
        println(f, "ElementDataFile = $(filename).raw")
    end
    open("$(filename).raw", "w") do f
        for val in voxels
            write(f, val)
        end
    end
end

function writeparaviewformat(surface, fname, scaletuple)
    @unpack mincorner, maxcorner, edgelength = scaletuple
    voxels = Array{Float64}(undef, edgelength, edgelength, edgelength)
    r = range(mincorner, stop=maxcorner, length=edgelength)
    for i in 1:edgelength, j in 1:edgelength, k in 1:edgelength
        voxels[i,j,k] = value(evaluate(surface, SVector{3}([r[i], r[j], r[k]])))
    end
    write_volume(voxels, fname)
end

function readobj(fname, scaletuple)
    @unpack mincorner, maxcorner, edgelength = scaletuple
    m = load(fname)
    vs = vertices(m)./(edgelength-1)
    vs = vs.*(abs(mincorner)+abs(maxcorner))
    diffv = fill(mincorner, 3)
    vs = [v+diffv for v in vs]
    ns = -1 .* normalize.(normals(m))
    return vs, ns
end

struct NiceTree{A<:AbstractArray}
    data::String
    children::A
end

# write to nicer format
AbstractTrees.children(tree::NiceTree) = tree.children
AbstractTrees.printnode(io::IO, tree::NiceTree) = print(io, tree.data)

function toNiceTree(node::CachedCSGNode)
    if isempty(node.children)
        return NiceTree(_name(node.data), [])
    else
        op = node.data
        cs = node.children
        if op === :complement
            return NiceTree(string(op), [toNiceTree(cs[1])])
        else
            return NiceTree(string(op), [toNiceTree(cs[1]), toNiceTree(cs[2])])
        end
    end
end

# export to JSON

function toJSON(fname, tree)
    try
        open(fname, "w") do io
            println(io, JSON.json(tree))
        end
    catch e
        @show e
    end
    return nothing
end

rmse(n, ps) = sqrt((evaldistance(n, ps)/length(ps)))


text:
"""
PAGE_EQREFS

Dictionary to keep track of equations that are labelled on a page to allow
references within the page.
"""
const PAGE_EQREFS = LittleDict{String, Int}()

"""
PAGE_EQREFS_COUNTER

Counter to keep track of equation numbers as they appear along the page, this
helps with equation referencing. (The `_XC0q` is just a random string to avoid
clashes).
"""
const PAGE_EQREFS_COUNTER = "COUNTER_XC0q"

"""
$(SIGNATURES)

Reset the PAGE_EQREFS dictionary.
"""
 function def_PAGE_EQREFS!()
    empty!(PAGE_EQREFS)
    PAGE_EQREFS[PAGE_EQREFS_COUNTER] = 0
    return nothing
end

"""
PAGE_BIBREFS

Dictionary to keep track of bibliographical references on a page to allow
citation within the page.
"""
const PAGE_BIBREFS = LittleDict{String, String}()

"""
$(SIGNATURES)

Reset the PAGE_BIBREFS dictionary.
"""
def_PAGE_BIBREFS!() = (empty!(PAGE_BIBREFS); nothing)

"""
$(SIGNATURES)

Given a latex command such as `\\eqref{abc}`, hash the reference (here `abc`),
check if the given dictionary `d` has an entry corresponding to that hash and
return the appropriate HTML for it.
"""
function form_href(lxc::LxCom, dname::String;
                   parens="("=>")", class="href")::String
    cont  = content(lxc.braces[1])   # "r1, r2, r3"
    refs  = strip.(split(cont, ",")) # ["r1", "r2", "r3"]
    names = refstring.(refs)
    nkeys = length(names)
    # construct the partial link with appropriate parens, it will be
    # resolved at the second pass (HTML pass) whence the introduction of {{..}}
    # inner will be "{{href $dn $hr1}}, {{href $dn $hr2}}, {{href $dn $hr3}}"
    # where $hr1 is the hash of r1 etc.
    inner = prod("{{href $dname $name}}$(ifelse(i < nkeys, ", ", ""))"
                    for (i, name) ∈ enumerate(names))
    # encapsulate in a span for potential css-styling
    return html_span(class, parens.first * inner * parens.second)
end

lx_eqref(lxc::LxCom, _) = form_href(lxc, "EQR";  class="eqref")
lx_cite(lxc::LxCom, _)  = form_href(lxc, "BIBR"; parens=""=>"", class="bibref")
lx_citet(lxc::LxCom, _) = form_href(lxc, "BIBR"; parens=""=>"", class="bibref")
lx_citep(lxc::LxCom, _) = form_href(lxc, "BIBR"; class="bibref")

function lx_label(lxc::LxCom, _)
    refs = content(lxc.braces[1]) |> strip |> refstring
    return ""
end

function lx_biblabel(lxc::LxCom, _)::String
    name = refstring(strip(content(lxc.braces[1])))
    PAGE_BIBREFS[name] = content(lxc.braces[2])
    return ""
end

function lx_toc(::LxCom, _)
    minlevel = LOCAL_VARS["mintoclevel"].first
    maxlevel = LOCAL_VARS["maxtoclevel"].first
    return "{{toc $minlevel $maxlevel}}"
end


text:
module GraphViz
    using Compat

    if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
        include("../deps/deps.jl")
    else
        error("GraphViz not properly installed. Please run Pkg.build(\"LLRBVisualize\").")
    end

    # Plugin Struct

    immutable gvdevice_engine_t
        initialize::Ptr{Void}
        format::Ptr{Void}
        finalize::Ptr{Void}
    end

    immutable Pointf
        x::Float64
        y::Float64
    end

    immutable gvdevice_features_t
        flags::Cint
        default_margin_x::Float64
        default_margin_y::Float64
        default_pagesize_x::Float64
        default_pagesize_y::Float64
        default_dpi_x::Float64
        default_dpi_y::Float64
    end

    immutable gvplugin_installed_t
        id::Cint
        ctype::Ptr{UInt8}
        quality::Cint
        engine::Ptr{gvdevice_engine_t}
        features::Ptr{gvdevice_features_t}
    end

    const API_render        = @compat Int32(0)
    const API_layout        = @compat Int32(1)
    const API_textlayout    = @compat Int32(2)
    const API_device        = @compat Int32(3)
    const API_loadimage     = @compat Int32(4)

    const EMIT_SORTED                   = (1<<0)
    const EMIT_COLORS                   = (1<<1)
    const EMIT_CLUSTERS_LAST            = (1<<2)
    const EMIT_PREORDER                 = (1<<3)
    const EMIT_EDGE_SORTED              = (1<<4)
    const GVDEVICE_DOES_PAGES           = (1<<5)
    const GVDEVICE_DOES_LAYERS          = (1<<6)
    const GVDEVICE_EVENTS               = (1<<7)
    const GVDEVICE_DOES_TRUECOLOR       = (1<<8)
    const GVDEVICE_BINARY_FORMAT        = (1<<9)
    const GVDEVICE_COMPRESSED_FORMAT    = (1<<10)
    const GVDEVICE_NO_WRITER            = (1<<11)
    const GVRENDER_Y_GOES_DOWN          = (1<<12)
    const GVRENDER_DOES_TRANSFORM       = (1<<13)
    const GVRENDER_DOES_ARROWS          = (1<<14)
    const GVRENDER_DOES_LABELS          = (1<<15)
    const GVRENDER_DOES_MAPS            = (1<<16)
    const GVRENDER_DOES_MAP_RECTANGLE   = (1<<17)
    const GVRENDER_DOES_MAP_CIRCLE      = (1<<18)
    const GVRENDER_DOES_MAP_POLYGON     = (1<<19)
    const GVRENDER_DOES_MAP_ELLIPSE     = (1<<20)
    const GVRENDER_DOES_MAP_BSPLINE     = (1<<21)
    const GVRENDER_DOES_TOOLTIPS        = (1<<22)
    const GVRENDER_DOES_TARGETS         = (1<<23)
    const GVRENDER_DOES_Z               = (1<<24)
    const GVRENDER_NO_WHITE_BG          = (1<<25)
    const LAYOUT_NOT_REQUIRED           = (1<<26)
    const OUTPUT_NOT_REQUIRED           = (1<<27)

    immutable gvplugin_api_t
        api::Cint
        types::Ptr{gvplugin_installed_t}
    end

    immutable gvplugin_library_t
        name::Ptr{UInt8}
        apis::Ptr{gvplugin_api_t}
    end

    # Job Struct

    immutable gvplugin_active_device_t
        engine::Ptr{gvdevice_engine_t}
        id::Cint
        features::Ptr{gvdevice_features_t}
        ctype::Ptr{UInt8}
    end

    immutable gvplugin_active_render_t
        engine::Ptr{Void}
        id::Cint
        features::Ptr{Void}
        ctype::Ptr{UInt8}
    end

    immutable gvplugin_active_loadimage_t
        engine::Ptr{Void}
        id::Cint
        ctype::Ptr{UInt8}
    end

    immutable gv_argvlist_t
        argv::Ptr{Ptr{UInt8}}
        argc::Cint;
        alloc::Cint;
    end

    immutable Point{T}
        x::T
        y::T
    end

    immutable Box{T}
        topleft::Point{T}
        bottomright::Point{T}
    end

    immutable gvdevice_callback_t
        refresh::Ptr{Void}          # void (*refresh) (GVJ_t * job);
        button_press::Ptr{Void}     # void (*button_press) (GVJ_t * job, int button, pointf pointer);
        button_release::Ptr{Void}   # void (*button_release) (GVJ_t * job, int button, pointf pointer);
        motion::Ptr{Void}           # void (*motion) (GVJ_t * job, pointf pointer);
        modify::Ptr{Void}           # void (*modify) (GVJ_t * job, const char *name, const char *value);
        del::Ptr{Void}              # void (*del) (GVJ_t * job);  /* can't use "delete" 'cos C++ stole it */
        read::Ptr{Void}             # void (*read) (GVJ_t * job, const char *filename, const char *layout);
        layout::Ptr{Void}           # void (*layout) (GVJ_t * job, const char *layout);
        render::Ptr{Void}           # void (*render) (GVJ_t * job, const char *format, const char *filename);
    end

    # TODO: These are probably wrong
    type GVCOMMON_s
        info::Ptr{Ptr{UInt8}}
        cmdname::Ptr{UInt8}
        verbose::Cint
        config::UInt8
        auto_outfile_names::UInt8
        errorfn::Ptr{Void}
        show_boxes::Ptr{Ptr{Void}}
        lib::Ptr{Ptr{Void}}
        viewNum::Cint
        builtins::Ptr{Void}
        demand_loading::Cint
    end

    type GVC_s
        common::GVCOMMON_s

        config_path::Ptr{UInt8}
        config_found::UInt8

        input_filenames::Ptr{Ptr{UInt8}}

        gvgs::Ptr{Void}
        gvg::Ptr{Void}

        # Hack until tuples are properly inlined into types
        apis0::Ptr{Void}
        apis1::Ptr{Void}
        apis2::Ptr{Void}
        apis3::Ptr{Void}
        apis4::Ptr{Void}
        api0::Ptr{Void}
        api1::Ptr{Void}
        api2::Ptr{Void}
        api3::Ptr{Void}
        api4::Ptr{Void}
        packages::Ptr{Void}

        #  size_t (*write_fn) (GVJ_t *job, const char *s, size_t len);
        write_fn::Ptr{Void}

        # More stuff I don't need right now
    end

    type GVJ_s
        gvc::Ptr{GVC_s}
        next::Ptr{GVJ_s}
        next_active::Ptr{GVJ_s}

        common::Ptr{Void}

        obj_state::Ptr{Void}
        input_filename::Ptr{UInt8}
        graph_index::Cint

        layout_type::Ptr{UInt8}

        output_filename::Ptr{UInt8}
        output_file::Ptr{Void}
        output_data::Ptr{UInt8}
        output_data_allocated::Cuint
        output_data_position::Cuint

        output_langname::Ptr{UInt8}
        output_lang::Ptr{Cint}

        render::gvplugin_active_render_t
        device::gvplugin_active_device_t
        loadimage::gvplugin_active_loadimage_t

        callbacks::Ptr{gvdevice_callback_t}

        device_dpi::Pointf
        device_sets_dpi::UInt8

        displat::Ptr{Void}
        screen::Cint

        context::Ptr{Void}
        external_context::UInt8

        imagedata::Ptr{UInt8}

        flags::Cint

        numLayers::Cint
        layerNum::Cint

        pagesArraySize::Point{Int32}
        pagesArrayFirst::Point{Int32}
        pagesArrayMajor::Point{Int32}
        pagesArrayMinor::Point{Int32}
        pagesArrayElem::Point{Int32}
        numPages::Cint

        bb::Box{Float64}
        pad::Point{Float64}
        clip::Box{Float64}
        pageBox::Box{Float64}
        pageSize::Point{Float64}
        focus::Point{Float64}

        zoom::Float64
        rotation::Cint

        view::Point{Float64}
        canvasBox::Box{Float64}
        margin::Point{Float64}

        dpi::Point{Float64}

        width::UInt32
        height::UInt32

        pageBoundingBox::Box{Int32}
        boundingBox::Box{Int32}

        scale::Point{Float64}
        translation::Point{Float64}
        devscale::Point{Float64}

        fit_mode::UInt8
        needs_refresh::UInt8
        click::UInt8
        has_grown::UInt8
        has_been_rendered::UInt8

        button::UInt8
        pointer::Point{Float64}
        oldpointer::Point{Float64}

        current_obj::Ptr{Void}
        selected_obj::Ptr{Void}

        active_tooltip::Ptr{UInt8}
        selected_href::Ptr{UInt8}

        selected_obj_type_name::gv_argvlist_t
        selected_obj_attributes::gv_argvlist_t

        window::Ptr{Void}

        keybindings::Ptr{Void}
        numkeys::Cint
        keycodes::Ptr{Void}
    end

    # Disciplines

    immutable Agmemdisc_s
        #void *(*open) (Agdisc_t*);  /* independent of other resources */
        open::Ptr{Void}
        #void *(*alloc) (void *state, size_t req);
        alloc::Ptr{Void}
        #void *(*resize) (void *state, void *ptr, size_t old, size_t req);
        resize::Ptr{Void}
        # void (*free) (void *state, void *ptr);
        free::Ptr{Void}
        # void (*close) (void *state);
        close::Ptr{Void}
    end

    immutable Agiddisc_s
        # void *(*open) (Agraph_t * g, Agdisc_t*);    /* associated with a graph */
        open::Ptr{Void}
        # long (*map) (void *state, int objtype, char *str, unsigned long *id, int createflag);
        map::Ptr{Void}
        # long (*alloc) (void *state, int objtype, unsigned long id);
        alloc::Ptr{Void}
        #void (*free) (void *state, int objtype, unsigned long id);
        free::Ptr{Void}
        #char *(*print) (void *state, int objtype, unsigned long id);
        print::Ptr{Void}
        #void (*close) (void *state);
        close::Ptr{Void}
        #void (*idregister) (void *state, int objtype, void *obj);
        idregister::Ptr{Void}
    end

    immutable Agiodisc_s
        # int (*afread) (void *chan, char *buf, int bufsize);
        afread::Ptr{Void}
        # int (*putstr) (void *chan, const char *str);
        putstr::Ptr{Void}
        # int (*flush) (void *chan);  /* sync */
        flush::Ptr{Void}
    end

    immutable Agdisc_s
        mem::Ptr{Agmemdisc_s}
        id::Ptr{Agiddisc_s}
        io::Ptr{Agiodisc_s}
    end

    function jl_afread(io::Ptr{Void}, buf::Ptr{UInt8}, bufsize::Cint)
        #@show (io,buf,bufsize)
        ret = readbytes!(unsafe_pointer_to_objref(io)::IO,pointer_to_array(buf,@compat Int(bufsize)))
        #@show ret
        convert(Cint,ret)
    end

    function jl_putstr(io::Ptr{Void}, str::Ptr{UInt8})
        #@show (io,str)
        convert(Cint,write(unsafe_pointer_to_objref(io)::IO,pointer_to_array(str,@compat Int(ccall(:strlen,Csize_t,(Ptr{UInt8},),str)))))::Cint
    end

    jl_flush(io::Ptr{Void}) = convert(Cint,0)

    const JuliaIODisc = [Agiodisc_s(
        cfunction(jl_afread,Cint,(Ptr{Void},Ptr{UInt8},Cint)),
        cfunction(jl_putstr,Cint,(Ptr{Void},Ptr{UInt8})),
        cfunction(jl_flush,Cint,(Ptr{Void},))
    )]

    null(::Type{gvplugin_installed_t}) = gvplugin_installed_t((@compat Int32(0)),convert(Ptr{UInt8},0),
        (@compat Int32(0)),convert(Ptr{gvdevice_engine_t},0),convert(Ptr{gvdevice_features_t},0))
    null(::Type{gvplugin_api_t}) = gvplugin_api_t((@compat Int32(0)),convert(Ptr{gvplugin_installed_t},0))

    # Memory interface

    # I/O interface

    # API - Context

    type Context
        handle::Ptr{Void}
        function Context()
            this = new(ccall((:gvContext,gvc),Ptr{Void},()))
            finalizer(this,free)
            this
        end
    end

    function free(t::Context)
        if t.handle != C_NULL
            ccall((:gvFreeContext,gvc), Void, (Ptr{Void},), t.handle)
        end
        t.handle = C_NULL
    end

    # API - Graph

    export Graph

    type Graph
        handle::Ptr{Void}
        didlayout::Bool
        function Graph(p::Ptr{Void})
            this = new(p,false)
            finalizer(this,free)
            this
        end
    end

    function free(g::Graph)
        if g.handle != C_NULL
            ccall((:agclose,cgraph), Cint, (Ptr{Void},), g.handle)
        end
        g.handle = C_NULL
    end

    Graph(graph::IO) = Graph(ccall((:agread,cgraph),Ptr{Void},(Any,Ptr{Void}),graph,[Agdisc_s(
        cglobal((:AgMemDisc,cgraph)),
        cglobal((:AgIdDisc,cgraph)),
        pointer(JuliaIODisc)
        )]))
    Graph(graph::Vector{UInt8}) = Graph(IOBuffer(graph))
    Graph(graph::AbstractString) = Graph(bytestring(graph).data)

    function layout!(g::Graph;engine="dot", context = default_context)
        @assert g.handle != C_NULL
        ccall((:gvLayout,gvc),Cint,(Ptr{Void},Ptr{Void},Ptr{UInt8}),context.handle,g.handle,engine)
        g.didlayout = true
    end

    render_x11(c::Context,g::Graph) = ccall((:gvRender,gvc),Cint,(Ptr{Void},Ptr{Void},Ptr{UInt8},Ptr{Void}),c.handle,g.handle,"x11",C_NULL)
    render_jobs(c::Context,g::Graph) = ccall((:gvRenderJobs,gvc),Cint,(Ptr{Void},Ptr{Void}),c.handle,g.handle)

    # Render

    # IO device

    type IODeviceState
        io::IO
        oldwritefn::Ptr{Void}
    end

    const active_devices = ObjectIdDict()

    function jlio_write(job::Ptr{Void},s::Ptr{UInt8},len::Csize_t)
        job = unsafe_load(convert(Ptr{GVJ_s},job))
        ioc = unsafe_pointer_to_objref(job.context)::IODeviceState
        write(ioc.io,pointer_to_array(s,@compat Int(len)))
        len #Julia doesn't do half things :)
    end

    # determined by counting bytes ;)
    const WRITEFN_OFFSET = 200

    function julia_io_initialize(firstjob::Ptr{Void})
        #@show firstjob
        firstjob = convert(Ptr{GVJ_s},firstjob)
        job = unsafe_load(firstjob)
        # Temporarily put in our custom write function
        ioc = unsafe_pointer_to_objref(job.context)::IODeviceState
        writefnptr = convert(Ptr{Ptr{Void}},job.gvc+WRITEFN_OFFSET)
        ioc.oldwritefn = unsafe_load(writefnptr)
        unsafe_store!(writefnptr,cfunction(jlio_write,Csize_t,(Ptr{Void},Ptr{UInt8},Csize_t)))
        # This function has void return
        nothing
    end
    function julia_io_finalize(firstjob::Ptr{Void})
        # Reset the write pointer we changed in julia_io_initialize
        firstjob = convert(Ptr{GVJ_s},firstjob)
        job = unsafe_load(firstjob)
        ioc = unsafe_pointer_to_objref(job.context)::IODeviceState
        writefnptr = convert(Ptr{Ptr{Void}},job.gvc+WRITEFN_OFFSET)
        unsafe_store!(writefnptr,ioc.oldwritefn)
        # Also remove it from the gc preserve dict
        haskey(active_devices,ioc) && pop!(active_devices,ioc)
        nothing
    end

    const default_context = GraphViz.Context()

    const julia_io_engine = [ gvdevice_engine_t(cfunction(julia_io_initialize,Void,(Ptr{Void},)),C_NULL,cfunction(julia_io_finalize,Void,(Ptr{Void},))) ]
    const julia_io_features = [ gvdevice_features_t((@compat Int32(GVDEVICE_DOES_TRUECOLOR|GVDEVICE_DOES_LAYERS)),0.,0.,0.,0.,72.,72.) ]
    const julia_io_name = "julia_io:svg".data
    const julia_io_libname = "julia_io".data
    const julia_io_device =
    [
      gvplugin_installed_t((@compat Int32(0)),pointer(julia_io_name), (@compat Int32(0)), pointer(julia_io_engine), pointer(julia_io_features));
      null(gvplugin_installed_t)
    ]
    const julia_io_api =
    [
        gvplugin_api_t(API_device, pointer(julia_io_device))
        null(gvplugin_api_t)
    ]

    add_julia_io!(c::Context) = ccall((:gvAddLibrary,gvc),Void,(Ptr{Void},Ptr{gvplugin_library_t}),c.handle,[gvplugin_library_t(pointer(julia_io_libname),pointer(julia_io_api))])

    function render(io::IO,g::GraphViz.Graph; context = default_context, format="julia_io:svg")
        GraphViz.add_julia_io!(context)
        if !g.didlayout
            error("Must call layout before calling render!")
        end
        state = IODeviceState(io,C_NULL)
        active_devices[state] = state
        ccall((:gvRenderContext,GraphViz.gvc),Cint,(Ptr{Void},Ptr{Void},Ptr{UInt8},Any),context.handle,g.handle,format,state)
    end

    import Base: writemime

    function writemime(io::IO, ::MIME"image/svg+xml", x::Graph)
        if !x.didlayout
            layout!(x,engine="dot")
        end
        render(io,x)
    end

    # Cairo device

    if isdir(Pkg.dir("Cairo"))
        using Cairo

        function cairo_initialize(firstjob::Ptr{Void})
            firstjob = convert(Ptr{GVJ_s},firstjob)
            job = unsafe_load(firstjob)
            if job.context != C_NULL
                c = unsafe_pointer_to_objref(job.context)::CairoContext
                job.context = c.ptr
                job.external_context    = 0x1
                job.width = width(c.surface)
                job.height = height(c.surface)
                unsafe_store!(firstjob,job)
            else
                job.external_context = 1
                unsafe_store!(firstjob,job)
                global last_surface = firstjob
            end
            nothing
        end

        global last_surface = nothing

        function cairo_finalize(firstjob::Ptr{Void})
            #=firstjob = convert(Ptr{GVJ_s},firstjob)
            job = unsafe_load(firstjob)
            if last_surface == firstjob
                surface = ccall((:cairo_get_target,Cairo._jl_libcairo),Ptr{Void},(Ptr{Void},),job.context)
                last_surface = CairoSurface(surface, job.width, job.height)
            end=#
            nothing
        end

        function cairo_format(firstjob::Ptr{Void})
            global last_surface
            firstjob = convert(Ptr{GVJ_s},firstjob)
            job = unsafe_load(firstjob)
            if last_surface == firstjob
                surface = ccall((:cairo_get_target,Cairo._jl_libcairo),Ptr{Void},(Ptr{Void},),job.context)
                last_surface = CairoSurface(surface, job.width, job.height)
            end
            nothing
        end

        const generic_cairo_engine = [ gvdevice_engine_t(cfunction(cairo_initialize,Void,(Ptr{Void},)),cfunction(cairo_format,Void,(Ptr{Void},)),cfunction(cairo_finalize,Void,(Ptr{Void},))) ]
        const generic_cairo_features = [ gvdevice_features_t((@compat Int32(0)),0.,0.,0.,0.,96.,96.) ]
        const generic_cairo_features_interactive = [ gvdevice_features_t((@compat Int32(0)),0.,0.,0.,0.,96.,96.) ]
        const generic_cairo_name = "julia:cairo".data
        const generic_cairo_libname = "julia:cairo".data
        const generic_cairo_device =
        [
          gvplugin_installed_t((@compat Int32(0)),pointer(generic_cairo_name), (@compat Int32(0)), pointer(generic_cairo_engine), pointer(generic_cairo_features));
          null(gvplugin_installed_t)
        ]
        const generic_cairo_api =
        [
            gvplugin_api_t(API_device, pointer(generic_cairo_device))
            null(gvplugin_api_t)
        ]

        add_julia_cairo!(c::Context) = ccall((:gvAddLibrary,gvc),Void,(Ptr{Void},Ptr{gvplugin_library_t}),c.handle,[gvplugin_library_t(pointer(generic_cairo_libname),pointer(generic_cairo_api))])

        function render(c::CairoContext,g::GraphViz.Graph; context = default_context, format="julia:cairo")
            GraphViz.add_julia_cairo!(context)
            if !g.didlayout
                error("Must call layout before calling render!")
            end
            ccall((:gvRenderContext,GraphViz.gvc),Cint,(Ptr{Void},Ptr{Void},Ptr{UInt8},Any),context.handle,g.handle,format,c)
        end

        function cairo_render(g::GraphViz.Graph; context = default_context, format="julia:cairo")
            global last_surface
            GraphViz.add_julia_cairo!(context)
            if !g.didlayout
                error("Must call layout before calling render!")
            end
            ccall((:gvRenderContext,GraphViz.gvc),Cint,(Ptr{Void},Ptr{Void},Ptr{UInt8},Ptr{Void}),context.handle,g.handle,format,C_NULL)
            surface = last_surface
            last_surface = nothing
            return surface
        end

        function writemime(io::IO, m::MIME"image/png", x::Graph)
            if !x.didlayout
                layout!(x,engine="dot")
            end
            writemime(io, m, cairo_render(x))
        end
        #=
        if isdir(Pkg.dir("Gtk"))
            using Gtk

            function gtk_initialize(firstjob::Ptr{Void})
                firstjob = convert(Ptr{GVJ_s},firstjob)
                job = unsafe_load(firstjob)
                c = unsafe_pointer_to_objref(job.context)::Gtk.Canvas
                c.data = firstjob
                job.context = getgc(c).ptr
                job.window = pointer_from_objref(c)
                unsafe_store!(firstjob,job)
                nothing
            end

            function gtk_finalize(firstjob::Ptr{Void})
                firstjob = convert(Ptr{GVJ_s},firstjob)
                job = unsafe_load(firstjob)
                c = unsafe_pointer_to_objref(job.window)::Gtk.Canvas
                draw(gtk_update,c)
                wait() #TODO: Make conditional on closing the widget
                nothing
            end

            function gtk_update(c::Canvas)
                jobp = c.data::Ptr{GVJ_s}
                job = unsafe_load(jobp)
                job.height              = height(c)
                job.width               = width(c)
                job.context             = getgc(c).ptr
                job.external_context    = 0x1
                unsafe_store!(jobp,job)
                println(job.callbacks)
                println(unsafe_load(job.callbacks).refresh)
                ccall(unsafe_load(job.callbacks).refresh,Void,(Ptr{Void},),jobp)
            end
            const gtk_engine = [ gvdevice_engine_t(cfunction(gtk_initialize,Void,(Ptr{Void},)),C_NULL,cfunction(gtk_finalize,Void,(Ptr{Void},))) ]
            const gtk_features = [ gvdevice_features_t(Int32(GVDEVICE_EVENTS),0.,0.,0.,0.,96.,96.) ]
            const gtk_name = "julia_gtk:cairo".data
            const gtk_libname = "julia_gtk:cairo".data
            const gtk_device =
            [
              gvplugin_installed_t(Int32(0),pointer(gtk_name), Int32(0), pointer(gtk_engine), pointer(gtk_features));
              null(gvplugin_installed_t)
            ]
            const gtk_api =
            [
                gvplugin_api_t(API_device, pointer(gtk_device))
                null(gvplugin_api_t)
            ]

            add_julia_gtk!(c::Context) = ccall((:gvAddLibrary,gvc),Void,(Ptr{Void},Ptr{gvplugin_library_t}),c.handle,[gvplugin_library_t(pointer(gtk_libname),pointer(gtk_api))])

            function render(c::Gtk.Canvas,cg::Context,g::Graph)
                @async begin
                    add_julia_gtk!(cg)
                    ccall((:gvRenderContext,gvc),Cint,(Ptr{Void},Ptr{Void},Ptr{UInt8},Any),cg.handle,g.handle,"julia_gtk",c)
                end
                nothing
            end
        end
        =#
    end

    graph_plugins(c::Context) = Graph(ccall((:gvPluginsGraph,gvc),Ptr{Void},(Ptr{Void},),c.handle))

    function listPlugins(c,kind)
        s = Array(Cint,1)
        r = ccall((:gvPluginList,gvc),Ptr{Ptr{UInt8}},(Ptr{Void},Ptr{UInt8},Ptr{Cint},Ptr{UInt8}),c.handle,kind,s,C_NULL)
        if r == C_NULL
            error("No Plugins available")
        end
        ret = Array(ByteString,s[1])
        for i = 1:s[1]
            ret[i] = bytestring(unsafe_load(r,i))
            c_free(unsafe_load(r,i))
        end
        c_free(r)
        ret
    end
end


text:
module GalewskyShallowWaterIMEX

using Gridap
using GridapGeosciences

include("GalewskyInitialConditions.jl")

# Solves the Galewsky test case for the shallow water equations on a sphere
# of physical radius 6371220m. Involves a shear flow instability of a zonal
# jet triggered by an initial gravity wave.
# reference:
#   Galewsky, Scott and Polvani (2004) Tellus, 56A 429-440

order  = 1
degree = 4

n      = 48
dt     = 60.0
nstep  = Int(24*60^2*20/dt) # 20 days
model = CubedSphereDiscreteModel(n; radius=rₑ)

hf, uf = shallow_water_imex_time_stepper(model, order, degree,
                                    h₀, u₀, f, g,
                                    dt, 0.5*dt, nstep;
                                    write_solution=true,
                                    write_solution_freq=240,
                                    write_diagnostics=true,
                                    write_diagnostics_freq=1,
                                    dump_diagnostics_on_screen=true)

end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
# ---
# title: 959. Regions Cut By Slashes
# id: problem959
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Union Find, Graph
# link: <https://leetcode.com/problems/regions-cut-by-slashes/description/>
# hidden: true
# ---
#
# In a N x N `grid` composed of 1 x 1 squares, each 1 x 1 square consists of a
# `/`, `\`, or blank space.  These characters divide the square into contiguous
# regions.
#
# (Note that backslash characters are escaped, so a `\` is represented as
# `"\\"`.)
#
# Return the number of regions.
#
#
#
# **Example 1:**
#
#
#
#     Input: [
#       " /",
#       "/ "
#     ]
#     Output: 2
#     Explanation: The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/1.png)
#
#
# **Example 2:**
#
#
#
#     Input: [
#       " /",
#       "  "
#     ]
#     Output: 1
#     Explanation: The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/2.png)
#
#
# **Example 3:**
#
#
#
#     Input: [
#       "\\/",
#       "/\\"
#     ]
#     Output: 4
#     Explanation: (Recall that because \ characters are escaped, "\\/" refers to \/, and "/\\" refers to /\.)
#     The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/3.png)
#
#
# **Example 4:**
#
#
#
#     Input: [
#       "/\\",
#       "\\/"
#     ]
#     Output: 5
#     Explanation: (Recall that because \ characters are escaped, "/\\" refers to /\, and "\\/" refers to \/.)
#     The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/4.png)
#
#
# **Example 5:**
#
#
#
#     Input: [
#       "//",
#       "/ "
#     ]
#     Output: 3
#     Explanation: The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/5.png)
#
#
#
#
# **Note:**
#
#   1. `1 <= grid.length == grid[0].length <= 30`
#   2. `grid[i][j]` is either `'/'`, `'\'`, or `' '`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module DEBmicroTrait

using FieldMetadata,
      Unitful,
      Roots,
      LinearAlgebra,
      Statistics,
      CSV,
      DataFrames,
      DataFramesMeta,
      Dierckx

import FieldMetadata: @default, @description, @units, @bounds, @logscaled, @flattenable, @plottable, @selectable, @chain,
                      default, description, units, bounds, logscaled, flattenable, plottable, selectable

# Field metadata columns
@chain columns @units @description

using Unitful: J, kJ, W, L, g, kg, cm, m, s, hr, d, mol, mmol, μmol
Unitful.register(@__MODULE__);
@unit bp "bp" Basepair 1u"m" true;

export AbstractMetabolism, AbstractMetabolismC, MetabolismC
export AbstractAssimilation, AbstractAssimilationC, AssimilationC
export AbstractTurnover, Turnover
export AbstractParams, Params
export AbstractIsolate, IsolateComposition
export AbstractSetup, Setup

include("components/metabolism.jl")
include("components/supeca.jl")
include("components/assimilation.jl")
include("components/turnover.jl")
include("components/thermostoichwizard.jl")
include("components/soil_properties.jl")
include("components/depolymerization.jl")
include("allometry/cell_composition.jl")
include("allometry/cell_assimilation.jl")
include("allometry/cell_metabolism.jl")
include("allometry/cell_turnover.jl")
include("allometry/cell_depolymerization.jl")
include("setup.jl")
include("model.jl")

end

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
# Check that package loads
using StatsKit

using Test

# TODO: run all tests in the same process after 'using StatsKit'
#       to check that they pass when all packages are loaded at the same time
#       This is not possible currently as Pkg.test starts a new environment

# Check that reexported packages do not export different objects under the same name
objects = Dict{Any, Symbol}()
for m in [Bootstrap, CategoricalArrays, Clustering, CSV, DataFrames,
          Distances, Distributions, GLM, HypothesisTests, KernelDensity,
          Loess, MultivariateStats, MixedModels, StatsBase, ShiftedArrays, TimeSeries]
    for name in names(m)
        # Work around undefined exports:
        # JuliaStats/StatsBase#663,
        # JuliaStats/Clustering.j#210,
        # JuliaStats/GLM.jl#405,
        # JuliaStats/StatsModels.jl#218,
        # JuliaStats/MultivariateStats.jl#142,
        # JuliaArrays/ShiftedArrays.jl#42,
        # JuliaStats/TimeSeries.jl#487
        name in (:findat, :kmeans_opts, :Gaussian, :LagTerm, :dropterm,
                 :invsqrtm!, :wmean!, :reduce_vec, :mapreduce_vec, :find) && continue
        push!(objects, getproperty(m, name) => name)
    end
end

exceptions = [:ci, :head, :lead, :predict, :predict, :rename!, :rename, :msd, :evaluate,
              :rmsd, :ci, :rename!, :lag, :lag, :tail, :colwise, :colwise, :head, :meanad,
              :msd, :transform, :tail, :rmsd, :rename, :transform, :evaluate, :lead]
objects2 = filter(((k, v),) -> !(v in exceptions), objects)
if !allunique(values(objects2))
    for (k, v) in objects2
        conflicts = filter(((k2, v2),) -> v == v2, objects2)
        if length(conflicts) > 1
            @show conflicts
        end
    end
    error("Conflicting functions not in known exceptions found")
end


text:
function parseDateTimeN(str)
    date, mmn = split(str, '.')
    date1, time1 = split(date,'T')

    time2 = parse.(Int64, split(time1, ':'))

    mmn1 = mmn * reduce(*, ["0" for i in 1:(9-length(mmn))])

    rd = reverse(digits(parse(Int, mmn1), pad = 9))

    t = reduce(vcat, [
        time2,
        parse(Int, reduce(*, string.(rd[1:3]))),
        parse(Int, reduce(*, string.(rd[4:6]))),
        parse(Int, reduce(*, string.(rd[7:9])))]
        )

    DateTimeN(Date(date1), Time(t...))
end

import Base:show

show(io::IO, dd::DateTimeN) = begin
    print(io, dd.d)
    print(io, dd.t)
end

DateTimeN(str::String) = parseDateTimeN(str)

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
nums = split(readline(), " ")
println(parse(Int, nums[1]) + parse(Int, nums[2]))

max_stars_repo_name: MKolman/advent-of-code

max_issues_repo_name: MKolman/advent-of-code

max_forks_repo_name: MKolman/advent-of-code


text:
"""
    mid_point(M, p, q, x)

Compute the mid point between p and q. If there is more than one mid point
of (not neccessarily minimizing) geodesics (i.e. on the sphere), the one nearest
to x is returned.
"""
mid_point(M::Manifold, p, q, ::Any) = mid_point(M, p, q)
mid_point!(M::Manifold, y, p, q, ::Any) = mid_point!(M, y, p, q)

function mid_point(M::Circle, p, q, x)
    if distance(M, p, q) ≈ π
        X = 0.5 * log(M, p, q)
        Y = log(M, p, x)
        return exp(M, p, (sign(X) == sign(Y) ? 1 : -1) * X)
    end
    return mid_point(M, p, q)
end

function mid_point(M::Sphere, p, q, x)
    if isapprox(M, p, -q)
        X = log(M, p, x) / distance(M, p, x) * π
    else
        X = log(M, p, q)
    end
    return exp(M, p, 0.5 * X)
end
function mid_point!(M::Sphere, y, p, q, x)
    if isapprox(M, p, -q)
        X = log(M, p, x) / distance(M, p, x) * π
    else
        X = log(M, p, q)
    end
    y .= exp(M, p, 0.5 * X)
    return y
end

@doc raw"""
    reflect(M, f, x)

reflect the point `x` from the manifold `M` at the point `f(x)` of the
function $f\colon \mathcal M \to \mathcal M$, i.e.,

````math
    \operatorname{refl}_f(x) = \operatorname{refl}_{f(x)}(x),
````
see also [`reflect`](@ref reflect(M::Manifold, p, x))`(M,p,x)`.
"""

reflect(M::Manifold, pr::Function, x) = reflect(M, pr(x), x)

@doc raw"""
    reflect(M, p, x)

reflect the point `x` from the manifold `M` at point `x`, i.e.

````math
    \operatorname{refl}_p(x) = \exp_p(-\log_p x).
````
where exp and log denote the exponential and logarithmic map on `M`.
"""
reflect(M::Manifold, p, x) = exp(M, p, -log(M, p, x))

@doc raw"""
    sym_rem(x,[T=π])

Compute symmetric remainder of `x` with respect to the interall 2*`T`, i.e.
`(x+T)%2T`, where the default for `T` is $π$
"""
function sym_rem(x::N, T=π) where {N<:Number}
    return (x ≈ T ? convert(N, -T) : rem(x, convert(N, 2 * T), RoundNearest))
end


text:
# ---
# title: 1103. Distribute Candies to People
# id: problem1103
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/distribute-candies-to-people/description/>
# hidden: true
# ---
#
# We distribute some number of `candies`, to a row of **`n = num_people`**
# people in the following way:
#
# We then give 1 candy to the first person, 2 candies to the second person, and
# so on until we give `n` candies to the last person.
#
# Then, we go back to the start of the row, giving `n + 1` candies to the first
# person, `n + 2` candies to the second person, and so on until we give `2 * n`
# candies to the last person.
#
# This process repeats (with us giving one more candy each time, and moving to
# the start of the row after we reach the end) until we run out of candies.  The
# last person will receive all of our remaining candies (not necessarily one
# more than the previous gift).
#
# Return an array (of length `num_people` and sum `candies`) that represents the
# final distribution of candies.
#
#
#
# **Example 1:**
#
#
#
#     Input: candies = 7, num_people = 4
#     Output: [1,2,3,1]
#     Explanation:
#     On the first turn, ans[0] += 1, and the array is [1,0,0,0].
#     On the second turn, ans[1] += 2, and the array is [1,2,0,0].
#     On the third turn, ans[2] += 3, and the array is [1,2,3,0].
#     On the fourth turn, ans[3] += 1 (because there is only one candy left), and the final array is [1,2,3,1].
#
#
# **Example 2:**
#
#
#
#     Input: candies = 10, num_people = 3
#     Output: [5,2,3]
#     Explanation:
#     On the first turn, ans[0] += 1, and the array is [1,0,0].
#     On the second turn, ans[1] += 2, and the array is [1,2,0].
#     On the third turn, ans[2] += 3, and the array is [1,2,3].
#     On the fourth turn, ans[0] += 4, and the final array is [5,2,3].
#
#
#
#
# **Constraints:**
#
#   * 1 <= candies <= 10^9
#   * 1 <= num_people <= 1000
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using CPLEX
# using GLPK
# using Cbc
# using Clp
using JuMP

include("seqata_test_model_util.jl")

# ===========
# Préparartion des arguments
args = Args.parse_commandline()
const TEST_LEVEL = 1
Args.set(:level, TEST_LEVEL) # sera passé temporairement à 0
Args.set(:outdir, "/tmp") # eqivalent à  -d/tmp
Args.set(:itermax, 300)

# ===========
inst = instance_build_mini10()
@test inst.nb_planes == 10

# Création d'une struct encapsulant la spécifiation d'un test LP
specs = [
   LpSolverSpec(:lp,  :cplex, LpTimingSolver,  CPLEX.Optimizer),
#    LpSolverSpec(:lp2, :cplex, Lp2TimingSolver, CPLEX.Optimizer),
#    LpSolverSpec(:lp3, :cplex, Lp3TimingSolver, CPLEX.Optimizer),
#    LpSolverSpec(:lp4, :cplex, Lp4TimingSolver, CPLEX.Optimizer),
]

for spec in specs
    test_name = "Test descente LP $(spec.algo):$(spec.external_mip_solver)"
    print("$test_name ... ")
    @testset "$test_name" begin
        Args.set(:level, 0)
        test_one_lp_descent(spec)
        Args.set(:level, TEST_LEVEL)
    end
    println(" => fait.")
end


text:
"""
    AbstractIndex

An `AbstractVector` subtype optimized for indexing.
"""
abstract type AbstractIndex{name,K,V<:Integer,Ks<:AbstractVector{K},Vs<:AbstractUnitRange{V}} <: AbstractUnitRange{V} end

Base.valtype(::Type{<:AbstractIndex{name,K,V,Ks,Vs}}) where {name,K,V,Ks,Vs} = V

"""
    values_type(::AbstractIndex)
"""
values_type(::T) where {T} = values_type(T)
values_type(::Type{<:AbstractIndex{name,K,V,Ks,Vs}}) where {name,K,V,Ks,Vs} = Vs

Base.keytype(::Type{<:AbstractIndex{name,K}}) where {name,K} = K

"""
    keys_type(::AbstractIndex)
"""
keys_type(::T) where {T} = keys_type(T)
keys_type(::Type{<:AbstractIndex{name,K,V,Ks,Vs}}) where {name,K,V,Ks,Vs} = Ks

Base.size(a::AbstractIndex) = (length(a),)

Base.first(a::AbstractIndex) = first(values(a))

Base.last(a::AbstractIndex) = last(values(a))

Base.step(a::AbstractIndex) = step(values(a))

Base.firstindex(a::AbstractIndex) = first(keys(a))

Base.lastindex(a::AbstractIndex) = last(keys(a))

Base.haskey(a::AbstractIndex{name,K}, key::K) where {name,K} = key in keys(a)

Base.allunique(a::AbstractIndex) = true

Base.length(a::AbstractIndex) = length(values(a))

Base.isempty(a::AbstractIndex) = isempty(values(a))

Base.in(a, itr::AbstractIndex) = in(x, values(a))

Base.eachindex(a::AbstractIndex) = keys(a)

Base.pairs(a::AbstractIndex) = Base.Iterators.Pairs(a, keys(a))

function StaticRanges.is_dynamic(::Type{T}) where {T<:AbstractIndex}
    return is_dynamic(values_type(T)) & is_dynamic(keys_type(T))
end

function StaticRanges.is_fixed(::Type{T}) where {T<:AbstractIndex}
    return is_fixed(values_type(T)) & is_fixed(keys_type(T))
end

function StaticRanges.is_static(::Type{T}) where {T<:AbstractIndex}
    return is_static(values_type(T)) & is_static(keys_type(T))
end

function StaticRanges.can_set_first(::Type{T}) where {T<:AbstractIndex}
    return can_set_first(values_type(T)) & can_set_first(keys_type(T))
end

function StaticRanges.can_set_last(::Type{T}) where {T<:AbstractIndex}
    return can_set_last(values_type(T)) & can_set_last(keys_type(T))
end

function StaticRanges.can_set_length(::Type{T}) where {T<:AbstractIndex}
    return can_set_length(values_type(T)) & can_set_length(keys_type(T))
end

# FIXME
#Base.Slice(a::AbstractIndex) = x

for (f) in (:(==), :isequal)
    @eval begin
        Base.$(f)(a::AbstractIndex, b::AbstractIndex) = $(f)(values(a), values(b))
        Base.$(f)(a::AbstractIndex, b::AbstractVector) = $(f)(values(a), b)
        Base.$(f)(a::AbstractVector, b::AbstractIndex) = $(f)(a, values(b))
        Base.$(f)(a::AbstractIndex, b::OrdinalRange) = $(f)(values(a), b)
        Base.$(f)(a::OrdinalRange, b::AbstractIndex) = $(f)(a, values(b))
    end
end

for f in (:+, :-)
    @eval begin
        function Base.$(f)(x::AbstractIndex, y::AbstractIndex)
            if same_type(x, y)
                return similar_type(x)(combine_keys(x, y), +(values(x), values(y)))
            else
                return $(f)(promote(x, y)...)
            end
        end

        Base.$(f)(x::AbstractIndex, y::AbstractVector) = $(f)(promote(x, y)...)

        Base.$(f)(x::AbstractVector, y::AbstractIndex) = $(f)(promote(x, y)...)

        Base.$(f)(x::AbstractIndex, y::AbstractUnitRange) = $(f)(promote(x, y)...)

        Base.$(f)(x::AbstractUnitRange, y::AbstractIndex) = $(f)(promote(x, y)...)
    end
end

function StaticRanges.set_length!(a::AbstractIndex, len::Int)
    can_set_length(a) || error("Cannot use set_length! for instances of typeof $(typeof(x)).")
    set_length!(keys(a), len)
    set_length!(values(a), len)
    return a
end

StaticRanges.Size(::Type{T}) where {T<:AbstractIndex} = Size(values_type(T))


text:
# We calculate the expected total number of collisions,
# from legitimate submissions: n - D + D((D-1)/D)^n
# n = # values, D = size of the output space
# source: http://matt.might.net/articles/counting-hash-collisions/
# To run this script install julia and run "julia collisions.jl".

bits = big"128"  # length of hash
ppl = 8000000000 # world population estimate
cens = 4*24*14   # CENs per person for two weeks, assuming rotation frequency of 15 minutes
n = ppl*cens     # total number of CENs for two weeks
H = 2^bits       # size of hash output space, uniform distribution is assumed

print(n - H + H*((H-1)/H)^n)


text:
function PascalDistribution(arg0::jint, arg1::jdouble)
    return PascalDistribution((jint, jdouble), arg0, arg1)
end

function cumulative_probability(obj::PascalDistribution, arg0::jint)
    return jcall(obj, "cumulativeProbability", jdouble, (jint,), arg0)
end

function get_number_of_successes(obj::PascalDistribution)
    return jcall(obj, "getNumberOfSuccesses", jint, ())
end

function get_numerical_mean(obj::PascalDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::PascalDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_probability_of_success(obj::PascalDistribution)
    return jcall(obj, "getProbabilityOfSuccess", jdouble, ())
end

function get_support_lower_bound(obj::PascalDistribution)
    return jcall(obj, "getSupportLowerBound", jint, ())
end

function get_support_upper_bound(obj::PascalDistribution)
    return jcall(obj, "getSupportUpperBound", jint, ())
end

function is_support_connected(obj::PascalDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_probability(obj::PascalDistribution, arg0::jint)
    return jcall(obj, "logProbability", jdouble, (jint,), arg0)
end

function probability(obj::PascalDistribution, arg0::jint)
    return jcall(obj, "probability", jdouble, (jint,), arg0)
end


text:
"""
    struct SobolSampler <: AbstractSamplingAlgorithm

*Experimental feature, not part of stable public API.*

Sample from Sobol sequence. Also see [Sobol.jl](https://github.com/stevengj/Sobol.jl).

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct SobolSampler{TR<:AbstractDensityTransformTarget} <: AbstractSamplingAlgorithm
    trafo::TR = PriorToUniform()
    nsamples::Int = 10^5
end
export SobolSampler

"""
    struct GridSampler <: AbstractSamplingAlgorithm

*Experimental feature, not part of stable public API.*

Sample from equidistantly distributed points in each dimension.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct GridSampler{TR<:AbstractDensityTransformTarget} <: AbstractSamplingAlgorithm
    trafo::TR = PriorToUniform()
    ppa::Int = 100
end
export GridSampler

function bat_sample_impl(
    rng::AbstractRNG,
    target::AnyDensityLike,
    algorithm::Union{SobolSampler, GridSampler}
)
    density_notrafo = convert(AbstractDensity, target)
    shaped_density, trafo = bat_transform(algorithm.trafo, density_notrafo)
    shape = varshape(shaped_density)
    density = unshaped(shaped_density)

    samples = _gen_samples(density, algorithm)

    logvals = map(logdensityof(density), samples)
    weights = exp.(logvals)

    vol = exp(BigFloat(log_volume(spatialvolume(var_bounds(density)))))
    est_integral = mean(weights) * vol
    # ToDo: Add integral error estimate

    samples_trafo = shape.(DensitySampleVector(samples, logvals, weight = weights))
    samples_notrafo = inverse(trafo).(samples_trafo)

    return (result = samples_notrafo, result_trafo = samples_trafo, trafo = trafo, integral = est_integral)
end

function _gen_samples(density::AbstractDensity, algorithm::SobolSampler)
    bounds = var_bounds(density)
    isinf(bounds) && throw(ArgumentError("SobolSampler doesn't support densities with infinite support"))
    sobol = Sobol.SobolSeq(bounds.vol.lo, bounds.vol.hi)
    p = vcat([[Sobol.next!(sobol)] for i in 1:algorithm.nsamples]...)
    return p
end

function _gen_samples(density::AbstractDensity, algorithm::GridSampler)
    bounds = var_bounds(density)
    isinf(bounds) && throw(ArgumentError("SobolSampler doesn't support densities with infinite support"))
    dim = totalndof(density)
    ppa = algorithm.ppa
    ranges = [range(bounds.vol.lo[i], bounds.vol.hi[i], length = trunc(Int, ppa)) for i in 1:dim]
    p = vec(collect(Iterators.product(ranges...)))
    return [collect(p[i]) for i in 1:length(p)]
end

"""
    struct PriorImportanceSampler <: AbstractSamplingAlgorithm

*Experimental feature, not part of stable public API.*

Importance sampler using IID samples from the prior.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct PriorImportanceSampler <: AbstractSamplingAlgorithm
    nsamples::Int = 10^5
end
export PriorImportanceSampler

function bat_sample_impl(
    rng::AbstractRNG,
    posterior::AbstractPosteriorDensity,
    algorithm::PriorImportanceSampler
)
    shape = varshape(posterior)

    prior = getprior(posterior)
    prior_samples = bat_sample(prior, IIDSampling(nsamples = algorithm.nsamples)).result
    unshaped_prior_samples = unshaped.(prior_samples)

    v = unshaped_prior_samples.v
    prior_weight = unshaped_prior_samples.weight
    posterior_logd = map(logdensityof(unshaped(posterior)), v)
    weight = exp.(posterior_logd - unshaped_prior_samples.logd) .* prior_weight

    est_integral = mean(weight)
    # ToDo: Add integral error estimate

    posterior_samples = shape.(DensitySampleVector(v, posterior_logd, weight = weight))

    return (result = posterior_samples, prior_samples = prior_samples, integral = est_integral)
end


text:
# This file was generated, do not modify it. # hide
capacity = select(data, [:country, :primary_fuel, :capacity_mw]);
first(capacity, 5)


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis

path = "examples/examples_files/atlas_simple.urdf"
mech = Mechanism(path, floating=false, g = -.5)

storage = simulate!(mech, 10., record = true)
visualize(mech, storage)


text:
# Utilities for linear programming

require("linprog.jl")

module SymbolicLP
import GLPK
using LinProgGLPK  # for lpsolve

import Base.(+), Base.(-), Base.(*), Base.(/), Base.split

export LPBlock, LPData, LPParsed, addconstraints, addobjective, lpeval, lpparse, lpsolve

## Problem specification syntax
# The total LP problem can be split up into blocks. The LPBlock type
#   stores data for a specific block.
# Variables within a block are referenced by a symbol, e.g., :height, :width
# Variables can be referenced across blocks as an expression :($obj[height]),
#   where obj is an LPBlock object
# Constraints can be expressed as :(height > 5) or :(height == 2*width) or
#   :($xlabel1[baseline] == $xlabel2[baseline]). Comparisons can be chained,
#   e.g., :(left < middle < right).
# The objective function to be minimized can be written :(6*width - height)
# There is automatic support for simple "soft" constraints,
#   :(soft(height >= 5, 7))
#   This allows height values less than 5, but any such values are penalized.
#   This is represented in the following way:
#     A new variable, which here will be called :delta
#     New "hard" (strict) constraints: :(height + delta >= 5), :(delta >= 0)
#     New term added to objective function: :(7*delta)
#   Soft constraints can be expressed as equalities,
#   :(soft(height == 2*width, 12))
#   This introduces two new variables, :deltapos and :deltaneg, four new
#     constraints, and an additional sum of two terms to the objective.
#   Soft constraints cannot be chained.
#   There cannot be any cross-references to auxillary variables for
#     soft constraints
# The LP is "compiled" (parsed) from LPBlocks using
#   lpparse(T, block1, block2, ...), merging all mutually-referencing blocks into
#   a single LP. This ultimately produces f, A, b, Aeq, ..., with the possible
#   exception of function calls (see next point). A parsed LP is stored in
#   the LPParsed type.
# Constraint expressions can include function calls, e.g., :(x <= sin(pi/2)).
#   These function calls are not evaluated when the LP problem is compiled; that
#   is done on a separate lpeval(parsedlp) step, resulting in an object of
#   type LPData. This has the actual, f, A, b, Aeq, ... that one feeds to the
#   LP solver.

# TODO? efficient enabling/disabling particular constraints. Do this by adding names for individual constrains in LPBlock.

type LPBlock
    syms::Vector{Symbol}
    constraints::Vector{Expr}
    objective::Vector{Expr}
end
LPBlock(syms::Vector{Symbol}, constraints::Vector{Expr}, obj::Expr) = LPBlock(syms, constraints, [obj])
LPBlock(syms::Vector{Symbol}, constraints::Vector{Expr}) = LPBlock(syms, constraints, Expr[])
LPBlock(syms::Vector{Symbol}) = LPBlock(syms, Expr[], Expr[])

function addconstraints(lpb::LPBlock, exa::Expr...)
    for ex in exa
        push!(lpb.constraints, ex)
    end
    lpb
end

function addobjective(lpb::LPBlock, ex::Expr)
    push!(lpb.objective, ex)
    lpb
end

type LPData{T}
    f::Vector{T}
    A::AbstractSparseMatrix{T}
    b::Vector{T}
    Aeq::AbstractSparseMatrix{T}
    beq::Vector{T}
    lb::Vector{T}
    ub::Vector{T}
end

type LPParsed{T}
    # Storage for explicit values
    vals::LPData{T}
    # Storage for function calls. Only one of vals.b[i] and funs.b[i] can
    # be non-zero.
    funs::LPData{Function}
    findex::Vector{Int}
    rowineqindex::Vector{Int}
    roweqindex::Vector{Int}
end

function lpeval{T}(::Type{T}, M::SparseMatrixCSC{Function})
    n = length(M.nzval)
    nzval = Array(T, n)
    nzvalfun = M.nzval
    for i = 1:n
        nzval[i] = nzvalfun[i]()
    end
    SparseMatrixCSC(M.m, M.n, copy(M.colptr), copy(M.rowval), nzval)
end
lpeval{T}(::Type{T}, fa::Vector{Function}) = T[convert(T, f()) for f in fa]
function stuff{T}(v::Vector{T}, indx::Vector{Int}, n::Int)
    a = zeros(T, n)
    a[indx] = v
    a
end

lpeval{T}(lpp::LPParsed{T}) = LPData{T}(
    isempty(lpp.findex) ? lpp.vals.f : lpp.vals.f + stuff(lpeval(T, lpp.funs.f), lpp.findex, length(lpp.vals.f)),
    lpp.vals.A + lpeval(T, lpp.funs.A),
    isempty(lpp.rowineqindex) ? lpp.vals.b : lpp.vals.b + stuff(lpeval(T, lpp.funs.b), lpp.rowineqindex, length(lpp.vals.b)),
    lpp.vals.Aeq + lpeval(T, lpp.funs.Aeq),
    isempty(lpp.roweqindex) ? lpp.vals.beq : lpp.vals.beq + stuff(lpeval(T, lpp.funs.beq), lpp.roweqindex, length(lpp.vals.beq)),
    lpp.vals.lb,
    lpp.vals.ub)

function lpsolve(lpd::LPData, pvtuple...)
    pv = {pvtuple...}
    intindx = find(pv[1:2:end] .== "Int") # test for integer programming
    useint = !isempty(intindx)
    local param
    if useint
        param = GLPK.IntoptParam()
        colkind = fill(GLPK.CV, length(lpd.f))
        colkind[pv[intindx[end]+1]] = GLPK.IV
        keep = trues(length(pv))
        keep[intindx] = false
        keep[intindx+1] = false
        pv = pv[keep]
    else
        param = GLPK.SimplexParam()
    end
    param["msg_lev"] = GLPK.MSG_ERR
    param["presolve"] = GLPK.OFF
    for i = 1:2:length(pv)
        param[pv[i]] = pv[i+1]
    end
    if useint
        return mixintprog(lpd.f, lpd.A, lpd.b, lpd.Aeq, lpd.beq, lpd.lb, lpd.ub, colkind, param)
    else
        return linprog_simplex(lpd.f, lpd.A, lpd.b, lpd.Aeq, lpd.beq, lpd.lb, lpd.ub, param)
    end
end

function lpsolve(lpb::LPBlock, pv...)
    lpp, chunks = lpparse(Float64, lpb)
    lpd = lpeval(lpp)
    z, x, flag = lpsolve(lpd, pv...)
    x = x[chunks[1]]
    z, x, flag
end

# Cannonical forms of comparison operators.
# For LP, <= is the same as <, >= is the same as >.
function getop(s::Symbol)
    s == :(<=) ? :(<) :
    s == :(>=) ? :(>) : s
end

# Parsing the LP problem
function parsesoft(blk::LPBlock)
    out = LPBlock(copy(blk.syms), Expr[], copy(blk.objective))
    for ex in blk.constraints
        if ex.head == :call && ex.args[1] == :soft
            exc::Expr = ex.args[2]
            pval = ex.args[3]
            if exc.head != :comparison || length(exc.args) != 3
                error("soft constraints must be simple comparisons (with <, >, <=, >=, or ==), but this constraint does not qualify: ", exc)
            end
            args = exc.args
            op = getop(args[2])
            if op == :(==)
                deltaneg = gensym()
                deltapos = gensym()
                push!(out.syms, deltaneg)
                push!(out.syms, deltapos)
                push!(out.constraints, :($deltaneg < 0))
                push!(out.constraints, :($deltapos > 0))
                push!(out.constraints, :($(args[1]) - $deltaneg > $(args[3])))
                push!(out.constraints, :($(args[1]) - $deltapos < $(args[3])))
                push!(out.objective, :($pval*($deltapos-$deltaneg)))
            else
                delta = gensym()
                push!(out.syms, delta)
                push!(out.constraints, expr(:comparison, Any[delta, op, 0])) # :($delta $op 0))
                push!(out.constraints, expr(:comparison, Any[:($(args[1]) + $delta), op, args[3]]))
                if op == :(<)
                    push!(out.objective, :(-$pval*$delta))
                else
                    push!(out.objective, :($pval*$delta))
                end
            end
        else
            # This is not a soft constraint, just copy it
            push!(out.constraints, ex)
        end
    end
    out
end

# Type for algebra on numbers and thunks
typealias NumOrThunk Union(Number, Function)

type LinearIndexExpr{T}
    indx::Vector{Int}
    coef::Vector{T}
    rhs::T
end
LinearIndexExpr{T}(::Type{T}, indx::Vector{Int}, coef::Vector, rhs) = LinearIndexExpr{T}(indx, convert(Array{T}, coef), convert(T, rhs))
indx(l::LinearIndexExpr) = l.indx
indx(l::Nothing) = Array(Int, 0)
coef{T}(::Type{T}, l::LinearIndexExpr{T}) = l.coef
coef{T}(::Type{T}, l::Nothing) = Array(T, 0)

function lpparse{T<:FloatingPoint}(::Type{T}, blocks::LPBlock...)
    nblocks = length(blocks)
    # Parse soft constraints
    sblocks = LPBlock[parsesoft(b) for b in blocks]
    # Build dictionaries to resolve references, including cross-block references
    len = Int[length(b.syms) for b in sblocks]
    offset = [0, cumsum(len)]
    blockindex = Dict{LPBlock, Int}(blocks, 1:nblocks)
    symindexes = Array(Dict{Symbol, Int}, nblocks)
    chunks = Array(Range1{Int}, nblocks)
    for i = 1:nblocks
        symindexes[i] = Dict(sblocks[i].syms, offset[i]+1:offset[i+1])
        chunks[i] = offset[i]+1:offset[i]+length(blocks[i].syms)
    end
    n = offset[end]  # total number of variables in LP problem

    eqrows = Array(LinearIndexExpr{NumOrThunk}, 0)
    ineqrows = Array(LinearIndexExpr{NumOrThunk}, 0)
    lb = fill(-inf(T), n)
    ub = fill(inf(T), n)
    objective = nothing
    for iblock in 1:nblocks
        block = sblocks[iblock]
        # Parse the constraints
        for ex in block.constraints
            if ex.head == :(=)
                error("Use ==, not =, in comparisons: ", ex)
            end
            if ex.head != :comparison
                error("Expression is not an equality or inequality: ex.head = ", ex.head, "\n  ex = ", ex)
            end
            # Parse the individual expressions
            pargs = Array(Any, length(ex.args))
            for i = 1:length(ex.args)
                if iseven(i)
                    pargs[i] = getop(ex.args[i])  # copy the comparison operator
                else
                    pargs[i] = resolve(NumOrThunk, ex.args[i], iblock, blockindex, symindexes)
                end
            end
            # Create the individual constraints from chained comparisons
            for i = 2:2:length(ex.args)
                op = pargs[i]
                darg = pargs[i-1] - pargs[i+1]
                if op == :(==)
                    # equality
                    push!(eqrows, darg)
                else
                    # inequality
                    gflag = op == :(>)
                    if length(darg.indx) == 0
                        println("Warning: expression that uses no variables: ", ex.args[i-1:i+1])
                    elseif length(darg.indx) == 1 && isa(darg.coef[1], Number) && isa(darg.rhs, Number)
                        # upper/lower bound
                        if darg.coef[1] < 0
                            gflag = !gflag
                        end
                        val = -darg.rhs/darg.coef[1]
                        j = darg.indx[1]
                        if gflag
                            lb[j] = max(val, lb[j])
                        else
                            ub[j] = min(val, ub[j])
                        end
                    else
                        # inequality
                        if gflag
                            darg.coef = _neg(darg.coef)
                            darg.rhs = _neg(darg.rhs)
                        end
                        push!(ineqrows, darg)
                    end
                end
            end
        end
        # Parse the objective
        for term in block.objective
            objective += resolve(NumOrThunk, term, iblock, blockindex, symindexes)
        end
    end
    # Split the thunks from the values
    # We maintain thunks as a smaller array, so we don't have so many ()->0 to evaluate
    valineqrows = Array(LinearIndexExpr{T}, length(ineqrows))
    thunkineqrows = Array(LinearIndexExpr{Function}, 0)
    thunkineqindx = Array(Int, 0)  # stores the row#s that actually have thunks
    valeqrows = Array(LinearIndexExpr{T}, length(eqrows))
    thunkeqrows = Array(LinearIndexExpr{Function}, 0)
    thunkeqindx = Array(Int, 0)
    for i = 1:length(ineqrows)
        valineqrows[i], tmp = split(T, ineqrows[i])
        if !(tmp === nothing)
            push!(thunkineqrows, tmp)
            push!(thunkineqindx, i)
        end
    end
    for i = 1:length(eqrows)
        valeqrows[i], tmp = split(T, eqrows[i])
        if !(tmp === nothing)
            push!(thunkeqrows, tmp)
            push!(thunkeqindx, i)
        end
    end
    valf, thunkf = split(T, objective)
    # Create the matrix representation
    rows, cols, val, valb = sparsedata(valineqrows)
    valA = sparse(rows, cols, val, length(ineqrows), n)
    rows, cols, val, valbeq = sparsedata(valeqrows)
    valAeq = sparse(rows, cols, val, length(eqrows), n)
    lpval = LPData(tovector(T, valf, n), valA, valb, valAeq, valbeq, lb, ub)
    rows, cols, val, thunkb = sparsedata(thunkineqrows)
    thunkA = sparse(thunkineqindx[rows], cols, val, length(ineqrows), n)
    rows, cols, val, thunkbeq = sparsedata(thunkeqrows)
    thunkAeq = sparse(thunkeqindx[rows], cols, val, length(eqrows), n)
    lpthunk = LPData(coef(Function, thunkf), thunkA, thunkb, thunkAeq, thunkbeq, Function[], Function[])
    LPParsed(lpval, lpthunk, indx(thunkf), thunkineqindx, thunkeqindx), chunks
end

# Resolve references and encode linear operations
resolve{T}(::Type{T}, arg::Number, ind::Int, blockindex::Dict{LPBlock, Int}, symindexes::Vector{Dict{Symbol, Int}}) = convert(T, arg)
resolve{T, Tv<:Number}(::Type{T}, arg::Vector{Tv}, ind::Int, blockindex::Dict{LPBlock, Int}, symindexes::Vector{Dict{Symbol, Int}}) = convert(Array{T}, arg)
resolve{T}(::Type{T}, arg::Function, ind::Int, blockindex::Dict{LPBlock, Int}, symindexes::Vector{Dict{Symbol, Int}}) = get(opmap, arg, arg)
function resolve{T}(::Type{T}, arg::Symbol, ind::Int, blockindex::Dict{LPBlock, Int}, symindexes::Vector{Dict{Symbol, Int}})
    i = get(symindexes[ind], arg, -1)
    if i == -1
        return eval(arg)  # try to evaluate it, perhaps it's a constant
    end
    LinearIndexExpr(T, [i], [one(T)], zero(T))
end
function resolve{T}(::Type{T}, ex::Expr, ind::Int, blockindex::Dict{LPBlock, Int}, symindexes::Vector{Dict{Symbol, Int}})
    local ret
    if ex.head == :ref
        i = get(blockindex, ex.args[1], -1)
        if i == -1
            error("Referenced object not resolvable. Is the list of blocks complete?")
        end
        ret = resolve(T, ex.args[2], i, blockindex, symindexes)
    elseif ex.head == :call
        if contains((:(+), :(-), :(*), :(/)), ex.args[1])
            n = length(ex.args)
            args = Array(Any, n)
            for i = 1:n
                args[i] = resolve(T, ex.args[i], ind, blockindex, symindexes)
            end
            ret = eval(expr(ex.head, args))
        else
            ret = () -> eval(ex)
        end
    end
    ret
end

# Types and operations for intermediate storage during parsing of constraints, references, and function calls
# Algebra with numbers and thunks. We don't want these exported, so give them private names.
_add(n1::Number, n2::Number) = n1+n2
_add(f::Function, n::Number) = () -> f()+n
_add(n::Number, f::Function) = () -> n+f()
_add(f1::Function, f2::Function) = () -> f1()+f2()
_sub(n1::Number, n2::Number) = n1-n2
_sub(f::Function, n::Number) = () -> f()-n
_sub(n::Number, f::Function) = () -> n-f()
_sub(f1::Function, f2::Function) = () -> f1()-f2()
_mul(n1::Number, n2::Number) = n1*n2
_mul(f::Function, n::Number) = () -> f()*n
_mul(n::Number, f::Function) = () -> n*f()
_mul(f1::Function, f2::Function) = () -> f1()*f2()
_div(n1::Number, n2::Number) = n1/n2
_div(f::Function, n::Number) = () -> f()/n
_div(n::Number, f::Function) = () -> n/f()
_div(f1::Function, f2::Function) = () -> f1()/f2()
_neg(n::Number) = -n
_neg(f::Function) = () -> -f()
_mul{T}(n::NumOrThunk, v::Vector{T}) = T[_mul(n, item) for item in v]
_neg{T}(v::Vector{T}) = T[_neg(item) for item in v]
_div{T}(v::Vector{T}, n::NumOrThunk) = T[_div(item, n) for item in v]
opmap = Dict(((+),(-),(*),(/)), (_add, _sub, _mul, _div))

(+){T}(l::LinearIndexExpr{T}, n::NumOrThunk) = LinearIndexExpr(T, l.indx, l.coef, _add(l.rhs, n))
(+)(n::NumOrThunk, l::LinearIndexExpr) = l+n
(+)(l::LinearIndexExpr, n::Nothing) = l
(+)(n::Nothing, l::LinearIndexExpr) = l
(*){T}(l::LinearIndexExpr{T}, n::NumOrThunk) = LinearIndexExpr(T, l.indx, _mul(n, l.coef), _mul(n, l.rhs))
(*)(n::NumOrThunk, l::LinearIndexExpr) = l*n
function (+){T}(l1::LinearIndexExpr{T}, l2::LinearIndexExpr{T})
    indx = Array(Int, 0)
    coef = Array(T, 0)
    ind1 = l1.indx
    ind2 = l2.indx
    i1 = 1
    i2 = 1
    ii1 = ind1[i1]
    ii2 = ind2[i2]
    while i1 <= length(ind1) || i2 <= length(ind2)
        if ii1 < ii2
            push!(indx, ii1)
            push!(coef, l1.coef[i1])
            i1 += 1
            ii1 = i1 > length(ind1) ? typemax(Int) : ind1[i1]
        elseif ii1 > ii2
            push!(indx, ii2)
            push!(coef, l2.coef[i2])
            i2 += 1
            ii2 = i2 > length(ind2) ? typemax(Int) : ind2[i2]
        else
            push!(indx, ii1)
            push!(coef, _add(l1.coef[i1], l2.coef[i2]))
            i1 += 1
            i2 += 1
            ii1 = i1 > length(ind1) ? typemax(Int) : ind1[i1]
            ii2 = i2 > length(ind2) ? typemax(Int) : ind2[i2]
        end
    end
    LinearIndexExpr(T, indx, coef, _add(l1.rhs,l2.rhs))
end
(-){T}(l::LinearIndexExpr{T}) = LinearIndexExpr(T, l.indx, _neg(l.coef), _neg(l.rhs))
(-)(l1::LinearIndexExpr, l2::LinearIndexExpr) = l1 + (-l2)
(-){T}(l::LinearIndexExpr{T}, n::NumOrThunk) = LinearIndexExpr(T, l.indx, l.coef, _sub(l.rhs,n))
(-){T}(n::NumOrThunk, l::LinearIndexExpr) = LinearIndexExpr(T, l.indx, _neg(l.coef), _sub(n,l.rhs))
(/){T}(l::LinearIndexExpr{T}, n::NumOrThunk) = LinearIndexExpr(T, l.indx, _div(l.coef, n), _div(l.rhs, n))

# Separate thunks from explicit values
getnum(x::Number) = x
getnum(x::Function) = 0
getthunk(x::Number) = () -> 0
getthunk(x::Function) = x
function split{T}(::Type{T}, l::LinearIndexExpr{NumOrThunk})
    retnum = LinearIndexExpr(T, l.indx, [getnum(c) for c in l.coef], getnum(l.rhs))
    isthunkcoef = [isa(c, Function) for c in l.coef]
    local retthunk
    if any(isthunkcoef) || isa(l.rhs, Function)
        coefthunk = [getthunk(c) for c in l.coef[isthunkcoef]]
        retthunk = LinearIndexExpr(Function, l.indx[isthunkcoef], coefthunk, getthunk(l.rhs))
    else
        retthunk = nothing
    end
    retnum, retthunk
end
# split{T, Tl<:Number}(::Type{T}, l::LinearIndexExpr{Tl}) = l, nothing
split{T}(::Type{T}, l::Nothing) = nothing, nothing

# Convert a list of LinearIndexExpr to format needed to create the sparse matrix
function sparsedata{T}(la::Vector{LinearIndexExpr{T}})
    len = [length(l.indx) for l in la]
    n = sum(len)
    rows = Array(Int, n)
    cols = Array(Int, n)
    val = Array(T, n)
    m = length(la)
    valb = Array(T, m)
    offset = 0
    for i = 1:m
        l = la[i]
        k = length(l.indx)
        r = offset+(1:k)
        rows[r] = i
        cols[r] = l.indx
        val[r] = l.coef
        valb[i] = _neg(l.rhs)
        offset += k
    end
    rows, cols, val, valb
end

function tovector(::Type{Function}, l::LinearIndexExpr{Function}, n::Int)
    f = fill(()->0, n)
    f[l.indx] = l.coef
    f
end
function tovector{T}(::Type{T}, l::LinearIndexExpr{T}, n::Int)
    f = zeros(T, n)
    f[l.indx] = l.coef
    f
end
tovector(::Type{Function}, l::Nothing, n::Int) = fill(()->0, n)
tovector{T}(::Type{T}, l::Nothing, n::Int) = zeros(T, n)

end # module


text:
#=
# 5. Viscous flow about a moving body
In this notebook we will demonstrate the simulation of a moving body. It is straightforward
to set up a moving body. The main caveat is that the simulation is slower,
because the integrator must update the operators continuously throughout the simulation.

We will demonstrate this on an oscillating flat plate.
=#

using ViscousFlow
#-
using Plots

#=
### Problem specification and discretization
For simplicity, we will not create a free stream in this problem. Everything
here is the usual.
=#
my_params = Dict()
my_params["Re"] = 200
#-
xlim = (-1.0,1.0)
ylim = (-1.0,1.0)
my_params["grid Re"] = 4.0
g = setup_grid(xlim,ylim,my_params)

Δs = surface_point_spacing(g,my_params)

#=
### Set up body
Set up the plate and place it at the origin. (We don't actually have
to move it, since it defaults to the origin, but it's helpful to put
this here in case we wish to initialize it differently.)
=#
body = Plate(1.0,Δs)
T = RigidTransform((0,0),0)
T(body)

#=
### Set the body motion
Now we specify the body motion. We will use oscillatory pitch-heave kinematics for this:
=#
a = 0.25 # location of pitch axis, a = 0.5 is leading edge
ϕp = -π/2  # phase lag of pitch
ϕh = 0.0  # phase lag of heave
A = 0.25  # amplitude/chord
fstar = 1/π # fc/U
α₀ = 0 # mean angle of attack
Δα = 10π/180 # amplitude of pitching
U₀ = 0.0 # translational motion (set to zero in place of free stream)
K = π*fstar # reduced frequency, K = πfc/U

oscil1 = RigidBodyTools.PitchHeave(U₀,a,K,ϕp,α₀,Δα,A,ϕh)
motion = RigidBodyMotion(oscil1)

# We can inspect the kinematics in this `motion` by plotting them:
plot(motion)

#=
### Define the boundary condition functions
Instead of using the default boundary condition functions, we define
special ones here that provide the instantaneous surface velocity (i.e. the velocity
of every surface point) from the prescribed
motion. Every surface has an "exterior" and "interior" side. For
a flat plate, these two sides are the upper and lower sides, and both sides
are next to the fluid, so both of them are assigned the prescribed velocity
of the plate. (For closed bodies, we would assign this velocity to only
one of the sides, and zero to the other side. We will see an example of this in a later case.)
We pack these into a special dictionary and
pass these to the system construction.
=#
function my_vsplus(t,base_cache,phys_params,motions)
  vsplus = zeros_surface(base_cache)
  surface_velocity!(vsplus,base_cache,motions,t)
  return vsplus
end

function my_vsminus(t,base_cache,phys_params,motions)
  vsminus = zeros_surface(base_cache)
  surface_velocity!(vsminus,base_cache,motions,t)
  return vsminus
end

bcdict = Dict("exterior" => my_vsplus, "interior" => my_vsminus)

#=
### Construct the system structure
Here, we supply the motion and boundary condition functions as additional arguments.
=#
sys = viscousflow_system(g,body,phys_params=my_params,motions=motion,bc=bcdict);

#=
Before we solve the problem, it is useful to note that the Reynolds number
we specified earlier may not be the most physically-meaningful Reynolds number.
More relevant in this problem is the Reynolds number based on the maximum
body speed.
=#
Umax, imax, tmax, bmax = maxlistvelocity(sys)
Re_eff = my_params["Re"]*Umax

#-
u0 = init_sol(sys)
tspan = (0.0,10.0)
integrator = init(u0,tspan,sys)

#=
### Solve
This takes longer per time step than it does for stationary bodies. Here, we only
run it for 1.5 time units just to demonstrate it.
=#
@time step!(integrator,1.5)

#=
### Examine the solution
Let's look at a few snapshots of the vorticity field. Note that the
plotting here requires us to explicitly call the [`surfaces`](@ref)
function to generate the instantaneous configuration of the plate.
=#
sol = integrator.sol
plt = plot(layout = (1,3), size = (800, 300), legend=:false)
tsnap = 0.5:0.5:1.5
for (i, t) in enumerate(tsnap)
    plot!(plt[i],vorticity(sol,sys,t),sys,layers=false,title="t = $(round(t,digits=2))",clim=(-5,5),levels=range(-5,5,length=30),color = :RdBu)
    plot!(plt[i],surfaces(sol,sys,t))
end
plt
# and the forces
sol = integrator.sol
fx, fy = force(sol,sys,1);
#-
plot(
plot(sol.t,2*fx,xlim=(0,Inf),ylim=(-3,3),xlabel="Convective time",ylabel="\$C_D\$",legend=:false),
plot(sol.t,2*fy,xlim=(0,Inf),ylim=(-6,6),xlabel="Convective time",ylabel="\$C_L\$",legend=:false),
    size=(800,350)
)


text:
function mean_sqrt_of_positive_diagonals(x)
    vals = diag(x)
    return mean(sqrt.(vals[vals .> 0]))^2
end

"""
Note that this is not exported with a different name to avoid confusion with
eigenvalue_clean(mat::Hermitian, eigenvalue_threshold::Real).
"""
function _eigenvalue_clean(mat::Hermitian, obs::Real)
    if sum(isnan.(mat)) > 0 return mat end # If someone inputs a matrix involving a NaN
    N = size(mat)[1]
    eigenvalues, eigenvectors = eigen(mat)
    sigma2 = 1 - maximum(eigenvalues)/N
    q = obs/N
    eigenvalue_threshold = (sigma2 * (1 + 1/q + 2*sqrt(1/q))) * mean_sqrt_of_positive_diagonals(mat)
    return eigenvalue_clean(eigenvalues, eigenvectors, eigenvalue_threshold)
end

"""
    eigenvalue_clean(eigenvalues::Vector{<:Real}, eigenvectors::Matrix{<:Real},
                     eigenvalue_threshold::Real)

This takes the small eigenvalues (with values below eigenvalue_threshold). It sets them to the
greater of their average or eigenvalue_threshold/(4*number_of_small_eigens). Then the matrix is reconstructed and returned (as a `Hermitian`)
### Inputs
* `eigenvalues` - The eigenvalues of a matrix.
* `eigenvectors` - The eigenvectors  of a matrix.
* `eigenvalue_threshold` - The threshold for a eigenvalue to be altered.
### Returns
* A `Hermitian`.

    eigenvalue_clean(mat::Hermitian, eigenvalue_threshold::Real)

This splits a matrix into its eigenvalues and eigenvectors. Then takes the small eigenvalues (with values below `eigenvalue_threshold`). It sets them to the
greater of their average or `eigenvalue_threshold/(4*number_of_small_eigens)`. Then the matrix is reconstructed and returned (as a `Hermitian`)
### Inputs
* `mat` - A matrix that you want to regularise with eigenvalue regularisation.
* `eigenvalue_threshold` - The threshold for a eigenvalue to be altered.
### Returns
* A `Hermitian`.

    eigenvalue_clean(mat::Hermitian, ts::SortedDataFrame)

Similarly to the above two methods these functions regularise a matrix by setting small eigenvalues to near zero.
The method of Laloux, Cizeau, Bouchaud & Potters 2000 is used to choose a threshold.
### Inputs
* `mat` - A matrix that you want to regularise with eigenvalue regularisation.
* `ts` - The tick data.
### Returns
* A `Hermitian`.

    eigenvalue_clean(covariance_matrix::CovarianceMatrix, ts::SortedDataFrame;
                     apply_to_covariance::Bool = true)
### Inputs
* `mat` - A matrix that you want to regularise with eigenvalue regularisation.
* `ts` - The tick data.
* `apply_to_covariance` Should regularisation be applied to the covariance matrix or the correlation matrix.
### Returns
* A `CovarianceMatrix`.

Note that if the input matrices include any NaN terms then regularisation is not possible. The matrix will be silently returned (as these NaNs will generally be
from upstream problems so it is useful to return the matrix rather than throw at this point).As a result outputs should be checked.

## References
Laloux, L., Cizeau, P., Bouchaud J. , Potters, M. 2000. "Random matrix theory and financial correlations" International Journal of Theoretical Applied FInance, 3, 391-397.
"""
function eigenvalue_clean(eigenvalues::Vector{<:Real},
                          eigenvectors::Matrix{<:Real}, eigenvalue_threshold::Real)
    number_of_small_eigens = sum(eigenvalues .< eigenvalue_threshold)
    av_small_eigens = max(eigenvalue_threshold/(4*number_of_small_eigens) , mean(eigenvalues[eigenvalues .< eigenvalue_threshold]) )
    eigenvalues[eigenvalues .< eigenvalue_threshold] .= av_small_eigens
    regularised_mat = construct_matrix_from_eigen(eigenvalues, eigenvectors)
    return regularised_mat
end
function eigenvalue_clean(mat::Hermitian, eigenvalue_threshold::Real)
    if sum(isnan.(mat)) > 0 return mat end # If someone inputs a matrix involving a NaN
    eigenvalues, eigenvectors = eigen(mat)
    return eigenvalue_clean(eigenvalues, eigenvectors, eigenvalue_threshold)
end
function eigenvalue_clean(mat::Hermitian, ts::SortedDataFrame)
    dims = size(mat)[1]
    obs = nrow(ts.df)/dims
    regularised_mat = _eigenvalue_clean(mat, obs)
    return Hermitian(regularised_mat)
end
function eigenvalue_clean(covariance_matrix::CovarianceMatrix, ts::SortedDataFrame;
                          apply_to_covariance::Bool = true)
    if apply_to_covariance
        regularised_covariance = eigenvalue_clean(covariance(covariance_matrix), ts)
        corr, vols = cov2cor_and_vol(regularised_covariance, covariance_matrix.time_period_per_unit, covariance_matrix.time_period_per_unit)
        return CovarianceMatrix(corr, vols, covariance_matrix.labels, covariance_matrix.time_period_per_unit)
    else
        return CovarianceMatrix(Hermitian(eigenvalue_clean(covariance_matrix.correlation, ts)), covariance_matrix.volatility, covariance_matrix.labels, covariance_matrix.time_period_per_unit)
    end
end


text:
using ImageMagick  # maybe 16.04 specific requirement
using Arena.Amphitheatre
using RoME
using GraffSDK

#Point the giraffe in the right direction
config = loadGraffConfig()
config.userId = ""
config.robotId = ""
config.sessionId = ""

#Set up amphis
visdatasets = AbstractAmphitheatre[]

# visualize factor graph
push!(visdatasets, BasicGraffPose(config,
                                    meanmax=:max,
                                    poseProp = plDrawProp(0.15, 0.05, RGBA(1,1,0,0.5)),
                                    landmarkProp = plDrawProp(0.2, 0.1, RGBA(0,1,0,0.5))))

# visualize april tags on poses
push!(visdatasets, GraffTagOnPose(config) )

# visualize point cloud on poses
push!(visdatasets, Amphitheatre.GraffCloudOnPose(config, CameraModel(640, 480, 387.205, [322.042, 238.544])))

# run visualzer
vis, vistask = visualize(visdatasets, quat=Amphitheatre.Quat(0.0,1.0,0.0,0.0))
vistask

## stop visualizer
stopAmphiVis!()


text:
# TODO replace all @adjoints in matrices by rrules once Tangent gradients become available
# https://github.com/FluxML/Zygote.jl/issues/603

function _check_mul(A::AbstractMatrix, b::AbstractVector)
    if size(A, 2) != length(b)
        DimensionMismatch(
            "Number of columns of A ($(size(A, 2))) must correspond with length of b ($(length(b)))"
        ) |> throw
    end
end

function _check_mul(A::AbstractMatrix, B::AbstractMatrix)
    if size(A, 2) != size(B, 1)
        DimensionMismatch(
              "Number of columns of A ($(size(A, 2))) must correspond with number of rows of B ($(size(B, 1)))"
        ) |> throw
    end
end

include("maybe_hot_vector.jl")
include("maybe_hot_matrix.jl")
include("ngram_matrix.jl")
include("pre_imputing_matrix.jl")
include("post_imputing_matrix.jl")

const MaybeHotArray{T} = Union{MaybeHotVector{T}, MaybeHotMatrix{T}}

const ImputingMatrix{T, U} = Union{PreImputingMatrix{T, U}, PostImputingMatrix{T, U}}
const PreImputingDense = Dense{T, <: PreImputingMatrix} where T
const PostImputingDense = Dense{T, <: PostImputingMatrix} where T

Base.zero(X::T) where T <: ImputingMatrix = T(zero(X.W), zero(X.ψ))
Base.similar(X::T) where T <: ImputingMatrix = T(similar(X.W), similar(X.ψ))

function _split_bc(bc::Base.Broadcast.Broadcasted{Broadcast.ArrayStyle{T}}) where T <: ImputingMatrix
    bc1, bc2 = _split_bc(bc.args)
    Base.broadcasted(bc.f, bc1...), Base.broadcasted(bc.f, bc2...)
end
_split_bc(A::ImputingMatrix) = (A.W, A.ψ)
_split_bc(x::Tuple) = tuple(zip(_split_bc.(x)...)...)
_split_bc(x) = (x, x)

Base.BroadcastStyle(::Type{T}) where T <: ImputingMatrix = Broadcast.ArrayStyle{T}()

function Base.copy(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{T}}) where T <: ImputingMatrix
    T(map(Base.copy, _split_bc(bc))...)
end
function Base.copyto!(A::T, bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{T}}) where T <: ImputingMatrix
    bc1, bc2 = _split_bc(bc)
    Base.copyto!(A.W, bc1)
    Base.copyto!(A.ψ, bc2)
    return A
end

_print_params(io::IO, A::PreImputingMatrix) = Base.print_array(io, A.ψ |> permutedims)
_print_params(io::IO, A::PostImputingMatrix) = Base.print_array(io, A.ψ)

function Base.print_array(io::IO, A::ImputingMatrix)
    println(io, "W:")
    Base.print_array(io, A.W)
    println(io, "\n\nψ:")
    _print_params(io, A)
end

Base.print_array(io::IO, A::NGramMatrix) = Base.print_array(io, A.S)

function Flux.update!(opt, x::ImputingMatrix, x̄)
    if !isnothing(x̄.W)
        x.W .-= Flux.Optimise.apply!(opt, x.W, x̄.W)
    end
    if !isnothing(x̄.ψ)
        x.ψ .-= Flux.Optimise.apply!(opt, x.ψ, x̄.ψ)
    end
end
function Flux.update!(x::ImputingMatrix, x̄)
    if !isnothing(x̄.W)
        x.W .-= x̄.W
    end
    if !isnothing(x̄.ψ)
        x.ψ .-= x̄.ψ
    end
end

function Base.show(io::IO, X::T) where T <: Union{ImputingMatrix, MaybeHotArray, NGramMatrix}
    if get(io, :compact, false)
        if ndims(X) == 1
            print(io, length(X), "-element ", nameof(T))
        else
            print(io, join(size(X), "×"), " ", nameof(T))
        end
    else
        _show_fields(io, X)
    end
end

"""
    preimputing_dense(in, out, σ)

Like `Flux.Dense`, but use a [`PreImputingMatrix`](@ref) instead of a standard matrix.

# Examples
```jldoctest
julia> d = preimputing_dense(2, 3)
[pre_imputing]Dense(2, 3)

julia> typeof(d.weight)
PreImputingMatrix{Float32, Matrix{Float32}, Vector{Float32}}

julia> typeof(d.bias)
Vector{Float32} (alias for Array{Float32, 1})
```

See also: [`PreImputingMatrix`](@ref), [`postimputing_dense`](@ref), [`PostImputingMatrix`](@ref).
"""
preimputing_dense(d::Dense) = Dense(PreImputingMatrix(d.weight), d.bias, d.σ)
preimputing_dense(args...) = preimputing_dense(Dense(args...))

"""
    postimputing_dense(d_in, d_out, σ)

Like `Flux.Dense`, but use a [`PostImputingMatrix`](@ref) instead of a standard matrix.

# Examples
```jldoctest
julia> d = postimputing_dense(2, 3)
[post_imputing]Dense(2, 3)

julia> typeof(d.weight)
PostImputingMatrix{Float32, Matrix{Float32}, Vector{Float32}}

julia> typeof(d.bias)
Vector{Float32} (alias for Array{Float32, 1})
```

See also: [`PostImputingMatrix`](@ref), [`preimputing_dense`](@ref), [`PreImputingMatrix`](@ref).
"""
postimputing_dense(d::Dense) = Dense(PostImputingMatrix(d.weight), d.bias, d.σ)
postimputing_dense(args...) = postimputing_dense(Dense(args...))

_name(::PreImputingMatrix) = "[pre_imputing]"
_name(::PostImputingMatrix) = "[post_imputing]"
function Base.show(io::IO, l::Dense{F, <:ImputingMatrix}) where F
    print(io, "$(_name(l.weight))Dense(", size(l.weight, 2), ", ", size(l.weight, 1))
    l.σ == identity || print(io, ", ", l.σ)
    print(io, ")")
end