text:
"""
    abstract type AbstractMetricParams{T} end

Abstract type used to dispatch different geodesic problems.
"""
abstract type AbstractMetricParams{T} end

# contains the full metric components (this type needed for DiffGeoSymbolics)
abstract type AbstractMetric{T} <: AbstractMatrix{T} end

metric_params(m::AbstractMetric{T}) where {T} =
    error("Not implemented for metric $(typeof(m))")

"""
    geodesic_eq(m::AbstractMetricParams{T}, u, v)
    geodesic_eq!(m::AbstractMetricParams{T}, u, v)

Calculate the acceleration components of the geodesic equation given a position `u`, a velocity `v`, and a metric `m`.
"""
geodesic_eq(m::AbstractMetricParams{T}, u, v) where {T} =
    error("Not implemented for metric parameters $(typeof(m))")
geodesic_eq!(m::AbstractMetricParams{T}, u, v) where {T} =
    error("Not implemented for metric parameters $(typeof(m))")

"""
    constrain(m::AbstractMetricParams{T}, u, v; μ::T=0.0f0)

Give time component which would constrain a velocity vector `v` at position `x` to be a
geodesic with mass `μ`.
"""
constrain(m::AbstractMetricParams{T}, u, v; μ::T = 0.0) where {T} =
    error("Not implemented for metric parameters $(typeof(m))")

"""
    on_chart(m::AbstractMetricParams{T}, u)

Check if point `u` is a valid point for the metric described by `m`.

Returns false is `u` is a singularity.
"""
on_chart(m::AbstractMetricParams{T}, u) where {T} = !(sum(u) ≈ 0)

"""
    inner_radius(m::AbstractMetricParams{T})

Return the innermost valid coordinate relative to the origin, for use in geodesic tracing.

This usually represents some property of the metric, e.g. event horizon radius in Kerr/Schwarzschild metrics, or
throat diameter in worm hole metrics.
"""
inner_radius(m::AbstractMetricParams{T}) where {T} = convert(T, 0.0)

"""
    metric_type(m::AbstractMetricParams{T})

Return the [`AbstractMetric`](@ref) type associated with the metric parameters `m`.
"""
metric_type(m::AbstractMetricParams{T}) where {T} =
    error("Not implemented for metric parameters $(typeof(m))")

"""
    metric(m::AbstractMetricParams{T}, u)

Numerically evaluate the corresponding metric for [`AbstractMetricParams`](@ref), given parameter values `m`
and some point `u`.
"""
metric(m::AbstractMetricParams{T}, u) where {T} =
    error("Not implemented for metric $(typeof(m))")

# do we actually want to support this?
# since if it's a symbolic matrix, you can subs other ways better?
#"""
#    metric(m::AbstractMetric{T}, u)
#
#Evaluate the metric at a point `u`.
#"""
#metric(m::AbstractMetric{T}, u) where {T} = error("Not implemented for metric $(typeof(m))")

export AbstractMetricParams,
    geodesic_eq, geodesic_eq!, constrain, on_chart, inner_radius, metric_type


text:
f = open("jawiki-country.txt", "r")
dict = Dict()

for line in readlines(f)
  if ismatch(r"^\|.+\s=\s",line)
    m = match(r"^\|(.+)\s=\s(.+)",line)
    dict[m.captures[1]] = m.captures[2]
  end
end

println("key\tvalue")
for k in keys(dict)
  println(k, "\t", dict[k])
end

ext: jl


text:
module Dataset

struct CacheConfig
    l1d_size::String
    l1i_size::String
    l1d_assoc::Int
    l1i_assoc::Int
    cacheline_size::Int
end

const base_cacheconfig = CacheConfig("8kB","8kB",4,4,32)

struct Gem5Build
    build_root::String
    gem5opt::String
    gem5fast::String
    se_script::String
end

struct MiBenchBuild
    build_root::String
end

function clone_gem5(build_root)
    cwd = pwd()
    cd(build_root)
    if isdir("gem5")
        println("gem5 directory exists, skipping")
    else
        run(`git clone https://github.com/emcores/gem5.git`)
    end
    cd(cwd)
end

function build_gem5_x86(build_root)
    gem5b = get_gem5build(build_root)
    if verify_gem5build(gem5b)
        println("gem5 build exists, returning")
    else
        cwd = pwd()
        cd(build_root)
        cd("gem5")
        run(`git checkout bcf6983bc605b884fba52ec74634bddfd395cd5e`)
        # patch gem5 SConstruct to allow building gem5.fast
        filename = "SConstruct"
        run(`git checkout $filename`)
        all_text = read(filename,String)
        open(filename, "w") do f
            write(f, replace(all_text,
                "'-Wno-unused-parameter'" =>
                "'-Wno-unused-parameter','-Wno-unused-variable'"))
         end
        # build gem5 using scons
        @time run(`scons build/X86/gem5.opt build/X86/gem5.fast -j$(Sys.CPU_THREADS)`)
        cd(cwd)
    end
    return gem5b
end

@inline function get_gem5build(build_root)
    gem5opt = joinpath(build_root,"gem5","build","X86","gem5.opt")
    gem5fast = joinpath(build_root,"gem5","build","X86","gem5.fast")
    gem5scriptse = joinpath(build_root,"gem5","configs","example","se.py")
    return Gem5Build(build_root, gem5opt, gem5fast, gem5scriptse)
end

@inline function verify_gem5build(gem5b::Gem5Build)
    if ispath(gem5b.gem5opt) && ispath(gem5b.gem5fast) && ispath(gem5b.se_script)
        return true
    else
        return false
    end
end

function clone_mibench(build_root)
    cwd = pwd()
    cd(build_root)
    if isdir("mibench")
        println("mibench directory exists, skipping")
    else
        run(`git clone https://github.com/emcores/MiBench.git`)
    end
    cd(cwd)
end

function build_mibench_x86(build_root)
    cwd = pwd()
    cd(build_root)

    cd("MiBench")
    # build MiBench using defaults
    println(pwd())
    juliaexe = Base.julia_exename()
    buildscript = "build.jl"
    cmd = `$juliaexe $buildscript`
    println(cmd)
    @time run(cmd)
    mibenchb = MiBenchBuild(
        build_root
    )
    cd(cwd)
    return mibenchb
end

function exec_gem5opt_command(gem5b::Gem5Build,
    cacheconfig::CacheConfig, outdir::String, exec::String, options::String, maxinsts::Int64=0)

    c = cacheconfig
    command = `$(gem5b.gem5opt) --outdir=$outdir
        $(gem5b.se_script) -c $exec -o "$options"
        --l1d_size $(c.l1d_size) --l1i_size $(c.l1i_size)
        --l1d_assoc $(c.l1d_assoc) --l1i_assoc $(c.l1i_assoc)
        --cacheline_size $(c.cacheline_size)
        --caches --cpu-type TimingSimpleCPU`
    if maxinsts > 0
        command = `$command --maxinsts $maxinsts`
    end
    run(command)
end

end

max_stars_repo_name: emcores/CacheConfigurations

max_issues_repo_name: emcores/CacheConfigurations

max_forks_repo_name: emcores/CacheConfigurations


text:
import GR

macro codes(n)
  return :( [GR.PATH_MOVETO; repeat([GR.PATH_CURVE4], 3 * $n)] )
end

j = [ 0.1186, 0.2254,
      0.1186, 0.2740, 0.1186, 0.3226, 0.1185, 0.3712,
      0.1185, 0.3737, 0.1190, 0.3750, 0.1218, 0.3757,
      0.1469, 0.3822, 0.1719, 0.3889, 0.1970, 0.3955,
      0.2013, 0.3966, 0.2013, 0.3966, 0.2013, 0.3924,
      0.2013, 0.3134, 0.2013, 0.2344, 0.2012, 0.1554,
      0.2012, 0.1385, 0.2020, 0.1215, 0.2004, 0.1046,
      0.1984, 0.0841, 0.1935, 0.0645, 0.1798, 0.0479,
      0.1678, 0.0332, 0.1515, 0.0256, 0.1328, 0.0223,
      0.1118, 0.0186, 0.0906, 0.0176, 0.0693, 0.0194,
      0.0546, 0.0207, 0.0405, 0.0240, 0.0284, 0.0328,
      0.0193, 0.0395, 0.0142, 0.0479, 0.0165, 0.0592,
      0.0188, 0.0709, 0.0271, 0.0775, 0.0383, 0.0811,
      0.0455, 0.0833, 0.0531, 0.0837, 0.0606, 0.0827,
      0.0664, 0.0819, 0.0715, 0.0795, 0.0758, 0.0757,
      0.0802, 0.0717, 0.0842, 0.0675, 0.0879, 0.0629,
      0.0906, 0.0593, 0.0936, 0.0559, 0.0974, 0.0534,
      0.1039, 0.0490, 0.1101, 0.0504, 0.1140, 0.0571,
      0.1166, 0.0616, 0.1172, 0.0667, 0.1177, 0.0716,
      0.1185, 0.0788, 0.1185, 0.0861, 0.1185, 0.0933,
      0.1186, 0.1373, 0.1186, 0.1814, 0.1186, 0.2254 ]

c1 = [ 0.2134, 0.4701,
       0.2134, 0.4421, 0.1899, 0.4195, 0.1609, 0.4195,
       0.1318, 0.4195, 0.1083, 0.4421, 0.1083, 0.4701,
       0.1083, 0.4980, 0.1318, 0.5207, 0.1609, 0.5207,
       0.1899, 0.5207, 0.2134, 0.4980, 0.2134, 0.4701 ]

u = [ 0.3999, 0.1517,
      0.3906, 0.1443, 0.3815, 0.1384, 0.3716, 0.1335,
      0.3386, 0.1172, 0.3063, 0.1197, 0.2750, 0.1374,
      0.2536, 0.1495, 0.2399, 0.1676, 0.2351, 0.1914,
      0.2344, 0.1951, 0.2340, 0.1989, 0.2340, 0.2027,
      0.2340, 0.2634, 0.2340, 0.3241, 0.2339, 0.3849,
      0.2339, 0.3882, 0.2348, 0.3890, 0.2382, 0.3890,
      0.2631, 0.3889, 0.2879, 0.3889, 0.3128, 0.3890,
      0.3158, 0.3890, 0.3165, 0.3882, 0.3165, 0.3854,
      0.3164, 0.3260, 0.3164, 0.2666, 0.3165, 0.2073,
      0.3165, 0.1878, 0.3279, 0.1737, 0.3473, 0.1692,
      0.3541, 0.1676, 0.3606, 0.1679, 0.3669, 0.1710,
      0.3805, 0.1776, 0.3907, 0.1877, 0.3992, 0.1995,
      0.4002, 0.2010, 0.3999, 0.2026, 0.3999, 0.2042,
      0.3999, 0.2644, 0.3999, 0.3246, 0.3998, 0.3848,
      0.3998, 0.3882, 0.4006, 0.3890, 0.4041, 0.3890,
      0.4291, 0.3889, 0.4541, 0.3889, 0.4790, 0.3890,
      0.4819, 0.3890, 0.4825, 0.3883, 0.4825, 0.3856,
      0.4825, 0.3007, 0.4825, 0.2158, 0.4825, 0.1309,
      0.4825, 0.1280, 0.4817, 0.1275, 0.4789, 0.1276,
      0.4537, 0.1276, 0.4285, 0.1276, 0.4032, 0.1276,
      0.4006, 0.1275, 0.3997, 0.1281, 0.3998, 0.1307,
      0.4000, 0.1374, 0.3999, 0.1441, 0.3999, 0.1517 ]

l = [ 0.5970, 0.3216,
      0.5970, 0.2583, 0.5970, 0.1950, 0.5971, 0.1317,
      0.5971, 0.1283, 0.5961, 0.1275, 0.5927, 0.1275,
      0.5679, 0.1277, 0.5430, 0.1277, 0.5182, 0.1276,
      0.5153, 0.1275, 0.5146, 0.1282, 0.5146, 0.1310,
      0.5147, 0.2508, 0.5147, 0.3707, 0.5146, 0.4905,
      0.5146, 0.4930, 0.5152, 0.4941, 0.5178, 0.4948,
      0.5429, 0.5014, 0.5679, 0.5081, 0.5930, 0.5149,
      0.5963, 0.5158, 0.5971, 0.5154, 0.5971, 0.5118,
      0.5970, 0.4484, 0.5970, 0.3850, 0.5970, 0.3216 ]

i = [ 0.7118, 0.2621,
      0.7118, 0.2187, 0.7118, 0.1752, 0.7119, 0.1318,
      0.7119, 0.1287, 0.7114, 0.1275, 0.7077, 0.1275,
      0.6827, 0.1277, 0.6577, 0.1277, 0.6327, 0.1275,
      0.6299, 0.1275, 0.6292, 0.1281, 0.6292, 0.1309,
      0.6293, 0.2111, 0.6293, 0.2914, 0.6292, 0.3716,
      0.6292, 0.3739, 0.6298, 0.3749, 0.6322, 0.3756,
      0.6576, 0.3822, 0.6831, 0.3889, 0.7084, 0.3958,
      0.7118, 0.3967, 0.7119, 0.3957, 0.7119, 0.3931,
      0.7118, 0.3494, 0.7118, 0.3058, 0.7118, 0.2621 ]

a = [ 0.8846, 0.1482,
      0.8821, 0.1462, 0.8801, 0.1448, 0.8782, 0.1432,
      0.8643, 0.1317, 0.8482, 0.1251, 0.8299, 0.1235,
      0.8115, 0.1219, 0.7934, 0.1228, 0.7760, 0.1298,
      0.7375, 0.1453, 0.7264, 0.1825, 0.7380, 0.2140,
      0.7437, 0.2296, 0.7554, 0.2412, 0.7693, 0.2506,
      0.7874, 0.2629, 0.8075, 0.2712, 0.8282, 0.2784,
      0.8460, 0.2846, 0.8639, 0.2905, 0.8822, 0.2948,
      0.8836, 0.2951, 0.8847, 0.2956, 0.8846, 0.2974,
      0.8842, 0.3109, 0.8855, 0.3244, 0.8831, 0.3378,
      0.8792, 0.3599, 0.8620, 0.3685, 0.8420, 0.3665,
      0.8404, 0.3663, 0.8388, 0.3661, 0.8371, 0.3658,
      0.8216, 0.3624, 0.8139, 0.3522, 0.8137, 0.3347,
      0.8136, 0.3294, 0.8129, 0.3242, 0.8115, 0.3191,
      0.8086, 0.3092, 0.8020, 0.3027, 0.7918, 0.2999,
      0.7809, 0.2969, 0.7698, 0.2968, 0.7592, 0.3008,
      0.7384, 0.3086, 0.7309, 0.3320, 0.7429, 0.3515,
      0.7501, 0.3631, 0.7608, 0.3711, 0.7730, 0.3772,
      0.7929, 0.3872, 0.8142, 0.3925, 0.8364, 0.3942,
      0.8605, 0.3961, 0.8845, 0.3953, 0.9081, 0.3896,
      0.9217, 0.3863, 0.9345, 0.3811, 0.9455, 0.3723,
      0.9597, 0.3608, 0.9671, 0.3457, 0.9672, 0.3282,
      0.9676, 0.2623, 0.9675, 0.1964, 0.9676, 0.1306,
      0.9676, 0.1281, 0.9668, 0.1275, 0.9644, 0.1276,
      0.9388, 0.1276, 0.9132, 0.1276, 0.8876, 0.1276,
      0.8853, 0.1275, 0.8845, 0.1280, 0.8846, 0.1304,
      0.8847, 0.1361, 0.8846, 0.1418, 0.8846, 0.1482 ]

ai = [ 0.8846, 0.2623,
       0.8693, 0.2559, 0.8551, 0.2489, 0.8423, 0.2396,
       0.8354, 0.2346, 0.8294, 0.2289, 0.8240, 0.2225,
       0.8105, 0.2065, 0.8177, 0.1816, 0.8294, 0.1716,
       0.8345, 0.1673, 0.8401, 0.1657, 0.8467, 0.1669,
       0.8518, 0.1679, 0.8567, 0.1696, 0.8612, 0.1722,
       0.8660, 0.1752, 0.8706, 0.1786, 0.8756, 0.1811,
       0.8830, 0.1847, 0.8851, 0.1899, 0.8848, 0.1978,
       0.8842, 0.2191, 0.8846, 0.2405, 0.8846, 0.2623 ]

c2 = [ 0.7293, 0.4652,
       0.7293, 0.4373, 0.7058, 0.4146, 0.6768, 0.4146,
       0.6478, 0.4146, 0.6242, 0.4373, 0.6242, 0.4652,
       0.6242, 0.4932, 0.6478, 0.5158, 0.6768, 0.5158,
       0.7058, 0.5158, 0.7293, 0.4932, 0.7293, 0.4652 ]

c3 = [ 0.8663, 0.4660,
       0.8663, 0.4380, 0.8428, 0.4154, 0.8138, 0.4154,
       0.7847, 0.4154, 0.7612, 0.4380, 0.7612, 0.4660,
       0.7612, 0.4939, 0.7847, 0.5166, 0.8138, 0.5166,
       0.8428, 0.5166, 0.8663, 0.4939, 0.8663, 0.4660 ]

c4 = [ 0.7974, 0.5768,
       0.7974, 0.5488, 0.7739, 0.5262, 0.7449, 0.5262,
       0.7158, 0.5262, 0.6923, 0.5488, 0.6923, 0.5768,
       0.6923, 0.6047, 0.7158, 0.6274, 0.7449, 0.6274,
       0.7739, 0.6274, 0.7974, 0.6047, 0.7974, 0.5768 ]

GR.setcolorrep(1, 0.14, 0.14, 0.14)
GR.setcolorrep(2, 0.251, 0.388, 0.847) # darker blue
GR.setcolorrep(3, 0.796, 0.235, 0.2)   # darker red
GR.setcolorrep(4, 0.584, 0.345, 0.698) # darker purple
GR.setcolorrep(5, 0.22, 0.596, 0.149)  # darker green
GR.setcolorrep(6, 0.4, 0.51, 0.878)    # lighter blue
GR.setcolorrep(7, 0.835, 0.388, 0.361) # lighter red
GR.setcolorrep(8, 0.667, 0.475, 0.757) # lighter purple
GR.setcolorrep(9, 0.376, 0.678, 0.318) # lighter green

GR.setviewport(0, 1, 0, 1)
GR.setwindow(0, 1, 0, 1)
GR.setfillintstyle(GR.INTSTYLE_SOLID)
GR.updatews()

for s in 0.1:0.01:0.5
  GR.clearws()
  GR.setviewport(0.5 - s, 0.5 + s, 0.5 - s, 0.5 + s)
  GR.setlinewidth(s * 8)

  GR.setfillcolorind(1)
  GR.drawpath(j,  @codes(19), 1)
  GR.setfillcolorind(6)
  GR.drawpath(c1, @codes(4), 1)
  GR.setlinecolorind(2)
  GR.drawpath(c1, @codes(4), 0)
  GR.setfillcolorind(1)

  GR.drawpath(u,  @codes(21), 1)

  GR.drawpath(l,  @codes(9), 1)

  GR.drawpath(i,  @codes(9), 1)

  GR.drawpath(a,  @codes(26), 1)
  GR.setfillcolorind(0)
  GR.drawpath(ai, @codes(8), 1)

  GR.setfillcolorind(7)
  GR.drawpath(c2, @codes(4), 1)
  GR.setlinecolorind(3)
  GR.drawpath(c2, @codes(4), 0)
  GR.setfillcolorind(8)
  GR.drawpath(c3, @codes(4), 1)
  GR.setlinecolorind(4)
  GR.drawpath(c3, @codes(4), 0)
  GR.setfillcolorind(9)
  GR.drawpath(c4, @codes(4), 1)
  GR.setlinecolorind(5)
  GR.drawpath(c4, @codes(4), 0)

  GR.updatews()

  sleep(0.01)
  end


text:
############################################################################################
"""
$(TYPEDEF)
Abstract super type for parameter constraints.
"""
abstract type AbstractConstraint end

#!NOTE: Needs to be overloaded for new AbstractConstraint
"Check if (constraint, val) combination is valid."
function _check end

"Construct a transform and inverse transform function for given (constraint, val) "
function construct_transform end

############################################################################################
"""
$(TYPEDEF)
Abstract super type for parameter transformations.
"""
abstract type AbstractTransformer end

"Constrain `val` with given `constraint`"
function constrain end

"Unconstrain `val` with given `constraint`"
function unconstrain end

"Compute log(abs(determinant(jacobian(`x`)))) for given transformer."
function log_abs_det_jac end

############################################################################################
"""
$(TYPEDEF)

Utility struct to define Parameter in a way ModelWrappers.jl can handle. Will be separated in ModelWrapper struct for type stability.

# Fields
$(TYPEDFIELDS)
"""
struct Param{A,B}
    val::A
    constraint::B
    function Param(_rng::Random.AbstractRNG, val::A, constraint::B) where {A,B}
        ArgCheck.@argcheck _check(_rng, constraint, val) string(
            "Val and constraint do not match for ", val, " and ", constraint, "."
        )
        return new{A,B}(val, constraint)
    end
end
Param(val::A, constraint::B) where {A,B} = Param(Random.GLOBAL_RNG, val, constraint)

############################################################################################
#!NOTE: Needs to be overloaded for new AbstractTransformer

############################################################################################
# Export
export
    Param,
    AbstractConstraint,
    AbstractTransformer,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
module Tulip

export Model, PrimalDualPoint, readmps

# Cholesky module
    include("Cholesky/Cholesky.jl")

# package code goes here
    include("model.jl")
    include("ipm.jl")
    include("readmps.jl")

end # module


text:
abstract type AbstractMesh1D end
abstract type AbstractPeriodicMesh1D <: AbstractMesh1D end

isperiodic(mesh::AbstractMesh1D) = false
isperiodic(mesh::AbstractPeriodicMesh1D) = true

function Base.show(io::IO, mesh::AbstractMesh1D)
  print(io,
    typeof(mesh), " with ", numcells(mesh), " cells in ", bounds(mesh))
end

"""
    UniformPeriodicMesh1D(xmin::Real, xmax::Real, Nx::Integer)
    UniformPeriodicMesh1D(; xmin::Real, xmax::Real, Nx::Integer)

A uniform periodic mesh in one space dimension of `Nx` cells between
`xmin` and `xmax`.
"""
struct UniformPeriodicMesh1D{T<:Real} <: AbstractPeriodicMesh1D
  xmin::T
  xmax::T
  Nx::Int

  Δx::T

  function UniformPeriodicMesh1D{T}(xmin::T, xmax::T, Nx::Int) where T
    Nx > 0 || throw(ArgumentError("The number of elements `Nx` must be positive [`Nx == $Nx`]."))
    xmin < xmax || throw(ArgumentError("`xmin` must be smaller than `xmax` [`xmin == $xmin, xmax == $xmax`]."))

    new(xmin, xmax, Nx, (xmax-xmin)/Nx)
  end
end

function UniformPeriodicMesh1D(_xmin::Real, _xmax::Real, Nx::Integer)
  xmin, xmax = promote(_xmin, _xmax)
  UniformPeriodicMesh1D{typeof(xmin)}(xmin, xmax, Int(Nx))
end

function UniformPeriodicMesh1D(; xmin::Real, xmax::Real, Nx::Integer)
  UniformPeriodicMesh1D(xmin, xmax, Nx)
end

"""
    UniformMesh1D(xmin::Real, xmax::Real, Nx::Integer)
    UniformMesh1D(; xmin::Real, xmax::Real, Nx::Integer)

A uniform mesh in one space dimension of `Nx` cells between `xmin` and `xmax`.
"""
struct UniformMesh1D{T<:Real} <: AbstractMesh1D
  xmin::T
  xmax::T
  Nx::Int

  Δx::T

  function UniformMesh1D{T}(xmin::T, xmax::T, Nx::Int) where T
    Nx > 0 || throw(ArgumentError("The number of elements `Nx` must be positive [`Nx == $Nx`]."))
    xmin < xmax || throw(ArgumentError("`xmin` must be smaller than `xmax` [`xmin == $xmin, xmax == $xmax`]."))

    new(xmin, xmax, Nx, (xmax-xmin)/Nx)
  end
end

function UniformMesh1D(_xmin::Real, _xmax::Real, Nx::Integer)
  xmin, xmax = promote(_xmin, _xmax)
  UniformMesh1D{typeof(xmin)}(xmin, xmax, Int(Nx))
end

function UniformMesh1D(; xmin::Real, xmax::Real, Nx::Integer)
  UniformMesh1D(xmin, xmax, Nx)
end

@inline numcells(mesh::Union{UniformMesh1D,UniformPeriodicMesh1D}) = mesh.Nx

function bounds(cell::Int, mesh::Union{UniformMesh1D,UniformPeriodicMesh1D})
  @unpack Δx = mesh

  xmin = mesh.xmin + (cell-1)*Δx
  xmax = mesh.xmin + cell*Δx

  xmin, xmax
end

function bounds(cell::Int, mesh::Union{UniformMesh1D{T},UniformPeriodicMesh1D{T}}) where {T<:AbstractFloat}
  @unpack Δx = mesh

  # fix for functions with possible discontinuities on the edges
  xmin = mesh.xmin + (cell-1)*Δx
  xmin = nextfloat(xmin)

  xmax = mesh.xmin + cell*Δx
  xmax = prevfloat(xmax)

  xmin, xmax
end

@inline bounds(mesh::Union{UniformMesh1D,UniformPeriodicMesh1D}) = mesh.xmin, mesh.xmax

function left_cell(cell::Int, mesh::AbstractMesh1D)
  N = numcells(mesh)
  @boundscheck begin
    @assert(1 <= cell <= N)
  end
  isperiodic(mesh) && cell == 1 ? N : cell-1
end

function right_cell(cell::Int, mesh::AbstractMesh1D)
  N = numcells(mesh)
  @boundscheck begin
    @assert(1 <= cell <= N)
  end
  isperiodic(mesh) && cell == N ? 1 : cell+1
end

struct UniformMeshGrid1D{T, Mesh<:Union{UniformMesh1D,UniformPeriodicMesh1D}, Grid<:AbstractVector{T}} <: AbstractArray{T,1}
  mesh::Mesh
  grid::Grid
  continuous::Bool
end

isperiodic(meshgrid::UniformMeshGrid1D) = isperiodic(meshgrid.mesh)
numcells(meshgrid::UniformMeshGrid1D) = numcells(meshgrid.mesh)
bounds(meshgrid::UniformMeshGrid1D) = bounds(meshgrid.mesh)
bounds(i::Int, meshgrid::UniformMeshGrid1D) = bounds(i, meshgrid.mesh)
iscontinuous(meshgrid::UniformMeshGrid1D) = meshgrid.continuous

function Base.show(io::IO, meshgrid::UniformMeshGrid1D)
  print(io,
    "UniformMeshGrid1D with ", numcells(meshgrid), " cells with ", length(meshgrid.grid), " nodes in ", bounds(meshgrid))
end

function Base.length(meshgrid::UniformMeshGrid1D)
  if iscontinuous(meshgrid)
    if isperiodic(meshgrid)
      numcells(meshgrid) * (length(meshgrid.grid) - 1)
    else
      numcells(meshgrid) * (length(meshgrid.grid) - 1) + 1
    end
  else
    numcells(meshgrid) * length(meshgrid.grid)
  end
end

function Base.getindex(meshgrid::UniformMeshGrid1D, i::Int)
  N = length(meshgrid)
  @boundscheck begin
    @argcheck i > 0
    @argcheck i <= N
  end
  if iscontinuous(meshgrid)
    num_nodes_per_cell = length(meshgrid.grid) - 1
    if i <= num_nodes_per_cell + 1
      cell = 1
      node = i
    else
      cell = (i - 2) ÷ num_nodes_per_cell + 1
      node = i - (cell-1) * num_nodes_per_cell
    end
  else
    num_nodes_per_cell = length(meshgrid.grid)
    cell = (i-1) ÷ num_nodes_per_cell + 1
    node = i - (cell-1) * num_nodes_per_cell
  end
  xmin, xmax = bounds(cell, meshgrid)
  ymin, ymax = first(meshgrid.grid), last(meshgrid.grid)
  (meshgrid.grid[node] - ymin) * (xmax - xmin) / (ymax - ymin) + xmin
end

Base.size(meshgrid::UniformMeshGrid1D) = (length(meshgrid),)

struct UniformNonperiodicCoupledOperator{T, Dtype<:AbstractNonperiodicDerivativeOperator{T}, Mesh<:UniformMesh1D{T}, MeshGrid<:UniformMeshGrid1D{T, Mesh}, Coupling<:Union{Val{:continuous}, Val{:plus}, Val{:central}, Val{:minus}}, DerOrder} <: AbstractNonperiodicDerivativeOperator{T}
  D::Dtype
  meshgrid::MeshGrid
  coupling::Coupling
  der_order::DerOrder

  function UniformNonperiodicCoupledOperator(D::Dtype, mesh::Mesh, coupling::Coupling) where {T, Dtype<:AbstractNonperiodicDerivativeOperator{T}, Mesh<:UniformMesh1D{T}, Coupling<:Union{Val{:continuous}, Val{:plus}, Val{:central}, Val{:minus}}}
    meshgrid = UniformMeshGrid1D(mesh, grid(D), coupling===Val(:continuous))
    if !(derivative_order(D) == 1 || (derivative_order(D) == 2 && coupling===Val(:continuous)))
      throw(ArgumentError("Not implemented yet"))
    end
    der_order = Val(derivative_order(D))
    new{T, Dtype, Mesh, typeof(meshgrid), Coupling, typeof(der_order)}(D, meshgrid, coupling, der_order)
  end
end

struct UniformPeriodicCoupledOperator{T, Dtype<:AbstractNonperiodicDerivativeOperator{T}, Mesh<:UniformPeriodicMesh1D{T}, MeshGrid<:UniformMeshGrid1D{T, Mesh}, Coupling<:Union{Val{:continuous}, Val{:plus}, Val{:central}, Val{:minus}}, DerOrder} <: AbstractPeriodicDerivativeOperator{T}
  D::Dtype
  meshgrid::MeshGrid
  coupling::Coupling
  der_order::DerOrder

  function UniformPeriodicCoupledOperator(D::Dtype, mesh::Mesh, coupling::Coupling) where {T, Dtype<:AbstractNonperiodicDerivativeOperator{T}, Mesh<:UniformPeriodicMesh1D{T}, Coupling<:Union{Val{:continuous}, Val{:plus}, Val{:central}, Val{:minus}}}
    meshgrid = UniformMeshGrid1D(mesh, grid(D), coupling===Val(:continuous))
    if !(derivative_order(D) == 1 || (derivative_order(D) == 2 && coupling===Val(:continuous)))
      throw(ArgumentError("Not implemented yet"))
    end
    der_order = Val(derivative_order(D))
    new{T, Dtype, Mesh, typeof(meshgrid), Coupling, typeof(der_order)}(D, meshgrid, coupling, der_order)
  end
end

const UniformCoupledOperator = Union{UniformNonperiodicCoupledOperator, UniformPeriodicCoupledOperator}

iscontinuous(cD::UniformCoupledOperator) = iscontinuous(cD.meshgrid)

function Base.show(io::IO, cD::UniformCoupledOperator)

  print(io, cD.D)
  if get(io, :compact, false) == false
    if iscontinuous(cD)
      print(io, "\ncoupled continuously")
    else
      print(io, "\ncoupled discontinuously (upwind: ", cD.coupling, ")")
    end
  end
  print(io, " on ")
  print(io, cD.meshgrid.mesh)
end

Base.eltype(cD::UniformCoupledOperator) = eltype(cD.D)
derivative_order(cD::UniformCoupledOperator) = derivative_order(cD.D)
accuracy_order(cD::UniformCoupledOperator) = accuracy_order(cD.D)
source_of_coefficients(cD::UniformCoupledOperator) = source_of_coefficients(cD.D)

function left_boundary_weight(cD::UniformCoupledOperator)
  @unpack D, meshgrid = cD
  @unpack mesh, grid = meshgrid
  ymin, ymax = first(grid), last(grid)
  cell = 1
  xmin, xmax = bounds(cell, mesh)
  factor = (xmax - xmin) / (ymax - ymin)
  factor * left_boundary_weight(D)
end

function right_boundary_weight(cD::UniformCoupledOperator)
  @unpack D, meshgrid = cD
  @unpack mesh, grid = meshgrid
  ymin, ymax = first(grid), last(grid)
  cell = numcells(mesh)
  xmin, xmax = bounds(cell, mesh)
  factor = (xmax - xmin) / (ymax - ymin)
  factor * right_boundary_weight(D)
end

function couple_continuously(D::AbstractNonperiodicDerivativeOperator, mesh::UniformMesh1D)
  UniformNonperiodicCoupledOperator(D, mesh, Val(:continuous))
end

function couple_continuously(D::AbstractNonperiodicDerivativeOperator, mesh::UniformPeriodicMesh1D)
  UniformPeriodicCoupledOperator(D, mesh, Val(:continuous))
end

function couple_discontinuously(D::AbstractNonperiodicDerivativeOperator, mesh::UniformMesh1D, coupling::Union{Val{:plus}, Val{:central}, Val{:minus}}=Val(:central))
  UniformNonperiodicCoupledOperator(D, mesh, coupling)
end

function couple_discontinuously(D::AbstractNonperiodicDerivativeOperator, mesh::UniformPeriodicMesh1D, coupling::Union{Val{:plus}, Val{:central}, Val{:minus}}=Val(:central))
  UniformPeriodicCoupledOperator(D, mesh, coupling)
end

grid(cD::UniformCoupledOperator) = cD.meshgrid

function mass_matrix(cD::UniformCoupledOperator)
  m = ones(eltype(cD), size(cD, 1))
  scale_by_mass_matrix!(m, cD)
  Diagonal(vec(m))
end

for (fname, op) in ((:scale_by_mass_matrix!, Base.:*), (:scale_by_inverse_mass_matrix!, Base.:/))
  @eval function $fname(u::AbstractVector, cD::UniformCoupledOperator, factor=true)
    @unpack D, meshgrid = cD
    @unpack mesh, grid = meshgrid
    ymin, ymax = first(grid), last(grid)

    if iscontinuous(meshgrid)
      num_nodes_per_cell = length(grid) - 1

      cell = 1
      xmin, xmax = bounds(cell, mesh)
      jac = (xmax - xmin) / (ymax - ymin)
      if isperiodic(mesh)
        next_cell = left_cell(cell, mesh)
        xmin, xmax = bounds(next_cell, mesh)
        next_jac = (xmax - xmin) / (ymax - ymin)
        u[1] *= $op(factor, jac * get_weight(D, 1) + next_jac * get_weight(D, num_nodes_per_cell+1))
      else
        u[1] *= $op(factor, jac * get_weight(D, 1))
      end
      for i in 2:num_nodes_per_cell
        u[i] *= $op(factor, jac * get_weight(D, i))
      end
      if !isperiodic(mesh) && numcells(mesh) == 1
        u[num_nodes_per_cell+1] *= $op(factor, jac * get_weight(D, num_nodes_per_cell+1))
      end
      if numcells(mesh) == 1
        return u
      end

      for cell in 2:numcells(mesh)-1
        xmin, xmax = bounds(cell, mesh)
        jac = (xmax - xmin) / (ymax - ymin)
        next_cell = left_cell(cell, mesh)
        xmin, xmax = bounds(next_cell, mesh)
        next_jac = (xmax - xmin) / (ymax - ymin)
        u[(cell-1)*num_nodes_per_cell+1] *= $op(factor, jac * get_weight(D, 1) + next_jac * get_weight(D, num_nodes_per_cell+1))
        for i in 2:num_nodes_per_cell
          u[(cell-1)*num_nodes_per_cell+i] *= $op(factor, jac * get_weight(D, i))
        end
      end

      cell = numcells(mesh)
      xmin, xmax = bounds(cell, mesh)
      jac = (xmax - xmin) / (ymax - ymin)
      next_cell = left_cell(cell, mesh)
      xmin, xmax = bounds(next_cell, mesh)
      next_jac = (xmax - xmin) / (ymax - ymin)
      u[(cell-1)*num_nodes_per_cell+1] *= $op(factor, jac * get_weight(D, 1) + next_jac * get_weight(D, num_nodes_per_cell+1))
      for i in 2:num_nodes_per_cell
        u[(cell-1)*num_nodes_per_cell+i] *= $op(factor, jac * get_weight(D, i))
      end
      if !isperiodic(mesh)
        u[(cell-1)*num_nodes_per_cell+num_nodes_per_cell+1] *= $op(factor, jac * get_weight(D, num_nodes_per_cell+1))
      end
    else # discontinuous coupling
      num_nodes_per_cell = length(grid)
      for cell in 1:numcells(mesh)
        xmin, xmax = bounds(cell, mesh)
        jac = (xmax - xmin) / (ymax - ymin)
        idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell
        $fname(view(u, idx), D, $op(factor, jac))
      end
    end

    u
  end
end

function integrate(func, u::AbstractVector, cD::UniformCoupledOperator)
  @unpack D, meshgrid = cD
  @unpack mesh, grid = meshgrid
  ymin, ymax = first(grid), last(grid)

  @inbounds if iscontinuous(meshgrid)
    num_nodes_per_cell = length(grid) - 1

    cell = 1
    xmin, xmax = bounds(cell, mesh)
    jac = (xmax - xmin) / (ymax - ymin)
    if isperiodic(mesh)
      next_cell = left_cell(cell, mesh)
      xmin, xmax = bounds(next_cell, mesh)
      next_jac = (xmax - xmin) / (ymax - ymin)
      res = func(u[1]) * (jac * get_weight(D, 1) + next_jac * get_weight(D, num_nodes_per_cell+1))
    else
      res = func(u[1]) * (jac * get_weight(D, 1))
    end
    @simd for i in 2:num_nodes_per_cell
      res += func(u[i]) * (jac * get_weight(D, i))
    end
    if !isperiodic(mesh) && numcells(mesh) == 1
      res += func(u[num_nodes_per_cell+1]) * (jac * get_weight(D, num_nodes_per_cell+1))
    end
    if numcells(mesh) == 1
      return res
    end

    for cell in 2:numcells(mesh)-1
      xmin, xmax = bounds(cell, mesh)
      jac = (xmax - xmin) / (ymax - ymin)
      next_cell = left_cell(cell, mesh)
      xmin, xmax = bounds(next_cell, mesh)
      next_jac = (xmax - xmin) / (ymax - ymin)
      res += func(u[(cell-1)*num_nodes_per_cell+1]) * (jac * get_weight(D, 1) + next_jac * get_weight(D, num_nodes_per_cell+1))
      @simd for i in 2:num_nodes_per_cell
        res += func(u[(cell-1)*num_nodes_per_cell+i]) * (jac * get_weight(D, i))
      end
    end

    cell = numcells(mesh)
    xmin, xmax = bounds(cell, mesh)
    jac = (xmax - xmin) / (ymax - ymin)
    next_cell = left_cell(cell, mesh)
    xmin, xmax = bounds(next_cell, mesh)
    next_jac = (xmax - xmin) / (ymax - ymin)
    res += func(u[(cell-1)*num_nodes_per_cell+1]) * (jac * get_weight(D, 1) + next_jac * get_weight(D, num_nodes_per_cell+1))
    @simd for i in 2:num_nodes_per_cell
      res += func(u[(cell-1)*num_nodes_per_cell+i]) * (jac * get_weight(D, i))
    end
    if !isperiodic(mesh)
      res += func(u[(cell-1)*num_nodes_per_cell+num_nodes_per_cell+1]) * (jac * get_weight(D, num_nodes_per_cell+1))
    end
  else # discontinuous coupling
    num_nodes_per_cell = length(grid)
    cell = 1
    xmin, xmax = bounds(cell, mesh)
    jac = (xmax - xmin) / (ymax - ymin)
    idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell
    res = jac * integrate(func, view(u, idx), D)
    for cell in 2:numcells(mesh)
      xmin, xmax = bounds(cell, mesh)
      jac = (xmax - xmin) / (ymax - ymin)
      idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell
      res += jac * integrate(func, view(u, idx), D)
    end
  end

  res
end

function mul!(dest::AbstractVector, cD::UniformCoupledOperator, u::AbstractVector, α=true)
  N, _ = size(cD)
  @boundscheck begin
    @argcheck N == length(u)
    @argcheck N == length(dest)
  end

  @unpack D, meshgrid, coupling, der_order = cD
  if coupling === Val(:continuous)
    mul!(dest, D, meshgrid, coupling, der_order, u, α)
    scale_by_inverse_mass_matrix!(dest, cD)
  else
    mul!(dest, D, meshgrid, coupling, der_order, u, α)
  end
  dest
end

function mul!(_dest::AbstractVector, D::AbstractNonperiodicDerivativeOperator, meshgrid::UniformMeshGrid1D, coupling::Union{Val{:plus}, Val{:central}, Val{:minus}}, der_order::Val{1}, _u::AbstractVector, α=true)
  @unpack mesh, grid = meshgrid
  dest = reshape(_dest, length(grid), numcells(mesh))
  u    = reshape(_u,    length(grid), numcells(mesh))
  ymin, ymax = first(grid), last(grid)
  half = one(eltype(D)) / 2

  cell = 1
  xmin, xmax = bounds(cell, mesh)
  factor = (ymax - ymin) / (xmax - xmin)
  mul!(view(dest, :, cell), D, view(u, :, cell), α*factor)
  if numcells(mesh) == 1 && !isperiodic(mesh)
    return _dest
  end
  if coupling === Val(:central)
    dest[end, cell] += half * (u[1, right_cell(cell, mesh)] - u[end, cell]) * α*factor / right_boundary_weight(D)
  elseif coupling === Val(:plus)
    dest[end, cell] +=        (u[1, right_cell(cell, mesh)] - u[end, cell]) * α*factor / right_boundary_weight(D)
  end
  if isperiodic(mesh)
    if coupling === Val(:central)
      dest[1  , cell] += half * (u[1, cell] - u[end,  left_cell(cell, mesh)]) * α*factor / left_boundary_weight(D)
    elseif coupling === Val(:minus)
      dest[1  , cell] +=        (u[1, cell] - u[end,  left_cell(cell, mesh)]) * α*factor / left_boundary_weight(D)
    end
  end
  if numcells(mesh) == 1
    return _dest
  end

  for cell in 2:numcells(mesh)-1
    xmin, xmax = bounds(cell, mesh)
    factor = (ymax - ymin) / (xmax - xmin)
    mul!(view(dest, :, cell), D, view(u, :, cell), α*factor)
    if coupling === Val(:central)
      dest[1  , cell] += half * (u[1, cell] - u[end,  left_cell(cell, mesh)]) * α*factor / left_boundary_weight(D)
      dest[end, cell] += half * (u[1, right_cell(cell, mesh)] - u[end, cell]) * α*factor / right_boundary_weight(D)
    elseif coupling === Val(:plus)
      dest[end, cell] +=        (u[1, right_cell(cell, mesh)] - u[end, cell]) * α*factor / right_boundary_weight(D)
    elseif coupling === Val(:minus)
      dest[1  , cell] +=        (u[1, cell] - u[end,  left_cell(cell, mesh)]) * α*factor / left_boundary_weight(D)
    end
  end

  cell = numcells(mesh)
  xmin, xmax = bounds(cell, mesh)
  factor = (ymax - ymin) / (xmax - xmin)
  mul!(view(dest, :, cell), D, view(u, :, cell), α*factor)
  if coupling === Val(:central)
    dest[1  , cell] += half * (u[1, cell] - u[end,  left_cell(cell, mesh)]) * α*factor / left_boundary_weight(D)
  elseif coupling === Val(:minus)
    dest[1  , cell] +=        (u[1, cell] - u[end,  left_cell(cell, mesh)]) * α*factor / left_boundary_weight(D)
  end
  if isperiodic(mesh)
    if coupling === Val(:central)
      dest[end, cell] += half * (u[1, right_cell(cell, mesh)] - u[end, cell]) * α*factor / right_boundary_weight(D)
    elseif coupling === Val(:plus)
      dest[end, cell] +=        (u[1, right_cell(cell, mesh)] - u[end, cell]) * α*factor / right_boundary_weight(D)
    end
  end

  _dest
end

function mul!(dest::AbstractVector, D::AbstractNonperiodicDerivativeOperator, meshgrid::UniformMeshGrid1D, coupling::Val{:continuous}, der_order::Val{1}, u::AbstractVector, α=true)
  @unpack mesh, grid = meshgrid
  ymin, ymax = first(grid), last(grid)
  num_nodes_per_cell = length(grid) - 1
  # TODO: remove these allocations?
  utmp = similar(u, length(grid))
  desttmp = similar(dest, length(grid))

  cell = 1
  xmin, xmax = bounds(cell, mesh)
  factor = (ymax - ymin) / (xmax - xmin)
  if isperiodic(mesh) && numcells(mesh) == 1
    utmp[1:end-1] .= u
    utmp[end]      = u[1]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor))
    @. dest = desttmp[1:end-1]
    dest[1] += desttmp[end]
  else
    idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell+1
    utmp .= u[idx]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor))
    @. dest[idx] = desttmp
  end
  if numcells(mesh) == 1
    return dest
  end

  for cell in 2:numcells(mesh)-1
    xmin, xmax = bounds(cell, mesh)
    factor = (ymax - ymin) / (xmax - xmin)
    idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell+1
    utmp .= u[idx]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor))
    dest[idx[1]] += desttmp[1]
    @. dest[idx[2:end]] = desttmp[2:end]
  end

  cell = numcells(mesh)
  xmin, xmax = bounds(cell, mesh)
  factor = (ymax - ymin) / (xmax - xmin)
  if isperiodic(mesh)
    idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell
    utmp[1:end-1] .= u[idx]
    utmp[end]      = u[1]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor))
    dest[idx[1]] += desttmp[1]
    @. dest[idx[2:end]] = desttmp[2:end-1]
    dest[1] += desttmp[end]
  else
    idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell+1
    utmp .= u[idx]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor))
    dest[idx[1]] += desttmp[1]
    @. dest[idx[2:end]] = desttmp[2:end]
  end

  dest
end

function mul!(dest::AbstractVector, D::AbstractNonperiodicDerivativeOperator, meshgrid::UniformMeshGrid1D, coupling::Val{:continuous}, der_order::Val{2}, u::AbstractVector, α=true)
  @unpack mesh, grid = meshgrid
  ymin, ymax = first(grid), last(grid)
  num_nodes_per_cell = length(grid) - 1
  # TODO: remove these allocations?
  utmp = similar(u, length(grid))
  desttmp = similar(dest, length(grid))

  cell = 1
  xmin, xmax = bounds(cell, mesh)
  factor1 = (ymax - ymin) / (xmax - xmin)
  factor = factor1^2
  if isperiodic(mesh) && numcells(mesh) == 1
    utmp[1:end-1] .= u
    utmp[end]      = u[1]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor1))
    desttmp[1]   +=  α * factor1 * derivative_left(D, utmp, Val(1))
    desttmp[end] -=  α * factor1 * derivative_right(D, utmp, Val(1))
    @. dest = desttmp[1:end-1]
    dest[1] += desttmp[end]
  else
    idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell+1
    utmp .= u[idx]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor1))
    if isperiodic(mesh)
      desttmp[1]   +=  α * factor1 * derivative_left(D, utmp, Val(1))
    end
    if numcells(mesh) > 1
      desttmp[end] -=  α * factor1 * derivative_right(D, utmp, Val(1))
    end
    @. dest[idx] = desttmp
  end
  if numcells(mesh) == 1
    return dest
  end

  for cell in 2:numcells(mesh)-1
    xmin, xmax = bounds(cell, mesh)
    factor1 = (ymax - ymin) / (xmax - xmin)
    factor = factor1^2
    idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell+1
    utmp .= u[idx]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor1))
    desttmp[1]   +=  α * factor1 * derivative_left(D, utmp, Val(1))
    desttmp[end] -=  α * factor1 * derivative_right(D, utmp, Val(1))
    dest[idx[1]] += desttmp[1]
    @. dest[idx[2:end]] = desttmp[2:end]
  end

  cell = numcells(mesh)
  xmin, xmax = bounds(cell, mesh)
  factor1 = (ymax - ymin) / (xmax - xmin)
  factor = factor1^2
  if isperiodic(mesh)
    idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell
    utmp[1:end-1] .= u[idx]
    utmp[end]      = u[1]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor1))
    desttmp[1]   +=  α * factor1 * derivative_left(D, utmp, Val(1))
    desttmp[end] -=  α * factor1 * derivative_right(D, utmp, Val(1))
    dest[idx[1]] += desttmp[1]
    @. dest[idx[2:end]] = desttmp[2:end-1]
    dest[1] += desttmp[end]
  else
    idx = (cell-1)*num_nodes_per_cell+1:cell*num_nodes_per_cell+1
    utmp .= u[idx]
    mul!(desttmp, D, utmp, α*factor)
    scale_by_mass_matrix!(desttmp, D, inv(factor1))
    desttmp[1]   +=  α * factor1 * derivative_left(D, utmp, Val(1))
    dest[idx[1]] += desttmp[1]
    @. dest[idx[2:end]] = desttmp[2:end]
  end

  dest
end

function lower_bandwidth(cD::UniformNonperiodicCoupledOperator)
  @unpack D = cD
  lower_bandwidth(D)
end

function upper_bandwidth(cD::UniformNonperiodicCoupledOperator)
  @unpack D = cD
  upper_bandwidth(D)
end


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

@testset "Basic constraint tests" begin
    model = MOIU.Model{Float64}()
    MOIT.basic_constraint_tests(model, MOIT.TestConfig())
end

# We need to test this in a module at the top level because it can't be defined
# in a testset. If it runs without error, then we're okay.
module TestExternalModel
    using MathOptInterface
    struct NewSet <: MathOptInterface.AbstractScalarSet end
    struct NewFunction <: MathOptInterface.AbstractScalarFunction end
    MathOptInterface.Utilities.canonicalize!(f::NewFunction) = f
    Base.copy(::NewFunction) = NewFunction()
    Base.copy(::NewSet) = NewSet()
    MathOptInterface.Utilities.@model(ExternalModel,
        (MathOptInterface.ZeroOne, NewSet,),
        (),
        (),
        (),
        (NewFunction,),
        (),
        (),
        ()
    )
    MathOptInterface.Utilities.@model(ExternalOptimizer,
        (MathOptInterface.ZeroOne, NewSet,),
        (),
        (),
        (),
        (NewFunction,),
        (),
        (),
        (),
        true
    )
end

@testset "Super-types" begin
    model = TestExternalModel.ExternalModel{Float64}()
    optimizer = TestExternalModel.ExternalOptimizer{Float64}()
    @test isa(model, MOIU.AbstractModelLike{Float64})
    @test !isa(model, MOIU.AbstractOptimizer{Float64})
    @test !isa(model, MOI.AbstractOptimizer)
    @test !isa(optimizer, MOIU.AbstractModelLike{Float64})
    @test isa(optimizer, MOIU.AbstractOptimizer{Float64})
    @test isa(optimizer, MOI.AbstractOptimizer)
end

@testset "External @model" begin
    model = TestExternalModel.ExternalModel{Float64}()
    c = MOI.add_constraint(
        model, TestExternalModel.NewFunction(), TestExternalModel.NewSet())
    @test typeof(c) == MOI.ConstraintIndex{TestExternalModel.NewFunction,
        TestExternalModel.NewSet}
    c2 = MOI.add_constraint(
        model, TestExternalModel.NewFunction(), MOI.ZeroOne())
    @test typeof(c2) ==
        MOI.ConstraintIndex{TestExternalModel.NewFunction, MOI.ZeroOne}
end

struct DummyFunction <: MOI.AbstractScalarFunction end
struct DummySet <: MOI.AbstractScalarSet end

@testset "Supports with $T" for T in [Float64, Int]
    model = MOIU.Model{T}()
    @testset "ObjectiveFunction" begin
        @test !MOI.supports(model, MOI.ObjectiveFunction{DummyFunction}())
        for F in [MOI.SingleVariable,
                  MOI.ScalarAffineFunction{T},
                  MOI.ScalarQuadraticFunction{T}]
            @test MOI.supports(model, MOI.ObjectiveFunction{F}())
        end
        U = Float32
        for F in [MOI.ScalarAffineFunction{U}, MOI.ScalarQuadraticFunction{U}]
            @test !MOI.supports(model, MOI.ObjectiveFunction{F}())
        end
    end
    @testset "`SingleVariable`-in-`S`" begin
        @test !MOI.supports_constraint(model, DummyFunction, DummySet)
        @test !MOI.supports_constraint(model, MOI.SingleVariable, DummySet)
        @test !MOI.supports_add_constrained_variable(model, DummySet)
        for S in [MOI.EqualTo{T}, MOI.GreaterThan{T}, MOI.LessThan{T},
                  MOI.Interval{T}, MOI.Integer, MOI.ZeroOne,
                  MOI.Semicontinuous{T}, MOI.Semiinteger{T}]
            @test MOI.supports_constraint(model, MOI.SingleVariable, S)
            @test MOI.supports_add_constrained_variable(model, S)
        end
        U = Float32
        for S in [MOI.EqualTo{U}, MOI.GreaterThan{U}, MOI.LessThan{U},
                  MOI.Interval{U}, MOI.Semicontinuous{U}, MOI.Semiinteger{U}]
            @test !MOI.supports_constraint(model, MOI.SingleVariable, S)
            @test !MOI.supports_add_constrained_variable(model, S)
        end
    end
end

@testset "Setting lower/upper bound twice" begin
    @testset "flag_to_set_type" begin
        err = ErrorException("Invalid flag `0x11`.")
        T = Int
        @test_throws err MOIU.flag_to_set_type(0x11, T)
        @test MOIU.flag_to_set_type(0x10, T) == MOI.Integer
        @test MOIU.flag_to_set_type(0x20, T) == MOI.ZeroOne
    end
    @testset "$T" for T in [Int, Float64]
        model = MOIU.Model{T}()
        MOIT.set_lower_bound_twice(model, T)
        MOIT.set_upper_bound_twice(model, T)
    end
end

@testset "Name test" begin
    MOIT.nametest(MOIU.Model{Float64}())
end

@testset "Valid test" begin
    MOIT.validtest(MOIU.Model{Float64}())
end

@testset "Delete test" begin
    MOIT.delete_test(MOIU.Model{Float64}())
end

@testset "Empty test" begin
    MOIT.emptytest(MOIU.Model{Float64}())
end

@testset "supports_constraint test" begin
    MOIT.supports_constrainttest(MOIU.Model{Float64}(), Float64, Int)
    MOIT.supports_constrainttest(MOIU.Model{Int}(), Int, Float64)
end

@testset "OrderedIndices" begin
    MOIT.orderedindicestest(MOIU.Model{Float64}())
end

@testset "Continuous Linear tests" begin
    config = MOIT.TestConfig(solve=false)
    exclude = ["partial_start"] # Model doesn't support VariablePrimalStart.
    MOIT.contlineartest(MOIU.Model{Float64}(), config, exclude)
end

@testset "Continuous Conic tests" begin
    config = MOIT.TestConfig(solve=false)
    MOIT.contconictest(MOIU.Model{Float64}(), config)
end

@testset "Quadratic functions" begin

    model = MOIU.Model{Int}()

    x, y = MOI.add_variables(model, 2)
    @test 2 == @inferred MOI.get(model, MOI.NumberOfVariables())

    f1 = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.([3], [x]), MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]), 7)
    c1 = MOI.add_constraint(model, f1, MOI.Interval(-1, 1))

    @test 1 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Int},MOI.Interval{Int}}())
    @test (@inferred MOI.get(model, MOI.ListOfConstraintIndices{MOI.ScalarQuadraticFunction{Int},MOI.Interval{Int}}())) == [c1]

    f2 = MOI.VectorQuadraticFunction(MOI.VectorAffineTerm.([1, 2, 2], MOI.ScalarAffineTerm.([3, 1, 2], [x, x, y])), MOI.VectorQuadraticTerm.([1, 1, 2], MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y])), [7, 3, 4])
    c2 = MOI.add_constraint(model, f2, MOI.PositiveSemidefiniteConeTriangle(3))

    @test 1 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.VectorQuadraticFunction{Int},MOI.PositiveSemidefiniteConeTriangle}())
    @test (@inferred MOI.get(model, MOI.ListOfConstraintIndices{MOI.VectorQuadraticFunction{Int},MOI.PositiveSemidefiniteConeTriangle}())) == [c2]

    loc = MOI.get(model, MOI.ListOfConstraints())
    @test length(loc) == 2
    @test (MOI.VectorQuadraticFunction{Int},MOI.PositiveSemidefiniteConeTriangle) in loc
    @test (MOI.VectorQuadraticFunction{Int},MOI.PositiveSemidefiniteConeTriangle) in loc

    c3 = MOI.add_constraint(model, f1, MOI.Interval(-1, 1))

    change_1 = MOI.ScalarConstantChange(9)
    f3 = MOIU.modify_function(f1, change_1)
    change_2 = MOI.ScalarCoefficientChange(y, 2)
    f3 = MOIU.modify_function(f3, change_2)

    @test !(MOI.get(model, MOI.ConstraintFunction(), c1) ≈ f3)
    @test !(MOI.get(model, MOI.ConstraintFunction(), c3) ≈ f3)
    MOI.set(model, MOI.ConstraintFunction(), c1, f3)
    MOI.modify(model, c3, change_1)
    MOI.modify(model, c3, change_2)
    F1 = MOI.get(model, MOI.CanonicalConstraintFunction(), c1)
    @test F1 ≈ f3
    @test F1 !== MOI.get(model, MOI.ConstraintFunction(), c1)
    @test MOI.Utilities.is_canonical(F1)
    F3 = MOI.get(model, MOI.CanonicalConstraintFunction(), c3)
    @test F3 ≈ f3
    @test F3 === MOI.get(model, MOI.ConstraintFunction(), c3)
    @test MOI.Utilities.is_canonical(F3)
    @test MOI.get(model, MOI.CanonicalConstraintFunction(), c1) ≈ f3
    @test MOI.Utilities.is_canonical(MOI.get(model, MOI.CanonicalConstraintFunction(), c1))

    f4 = MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 1, 2], MOI.ScalarAffineTerm.([2, 4, 3], [x, y, y])), [5, 7])
    c4 = MOI.add_constraint(model, f4, MOI.SecondOrderCone(2))
    @test 1 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Int},MOI.SecondOrderCone}())

    f5 = MOI.VectorOfVariables([x])
    c5 = MOI.add_constraint(model, f5, MOI.RotatedSecondOrderCone(1))
    @test 1 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.RotatedSecondOrderCone}())

    f6 = MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([2, 9], [x, y])), [6, 8])
    c6 = MOI.add_constraint(model, f6, MOI.SecondOrderCone(2))
    @test 2 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Int},MOI.SecondOrderCone}())

    f7 = MOI.VectorOfVariables([x, y])
    c7 = MOI.add_constraint(model, f7, MOI.Nonpositives(2))
    @test 1 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.Nonpositives}())

    f8 = MOI.VectorOfVariables([x, y])
    c8 = MOI.add_constraint(model, f8, MOI.SecondOrderCone(2))
    @test 1 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.SecondOrderCone}())

    loc1 = MOI.get(model, MOI.ListOfConstraints())
    loc2 = Vector{Tuple{DataType, DataType}}()
    function _pushloc(constrs::Vector{MOIU.ConstraintEntry{F, S}}) where {F, S}
        if !isempty(constrs)
            push!(loc2, (F, S))
        end
    end
    MOIU.broadcastcall(_pushloc, model)
    for loc in (loc1, loc2)
        @test length(loc) == 6
        @test (MOI.VectorQuadraticFunction{Int},MOI.PositiveSemidefiniteConeTriangle) in loc
        @test (MOI.VectorQuadraticFunction{Int},MOI.PositiveSemidefiniteConeTriangle) in loc
        @test (MOI.VectorOfVariables,MOI.RotatedSecondOrderCone) in loc
        @test (MOI.VectorAffineFunction{Int},MOI.SecondOrderCone) in loc
        @test (MOI.VectorOfVariables,MOI.Nonpositives) in loc
    end

    @test MOI.is_valid(model, c4)
    MOI.delete(model, c4)
    @test !MOI.is_valid(model, c4)

    @test 1 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Int},MOI.SecondOrderCone}())
    @test MOI.get(model, MOI.ConstraintFunction(), c6).constants == f6.constants

    fx = MOI.SingleVariable(x)
    fy = MOI.SingleVariable(y)
    obj = 1fx + 2fy
    MOI.set(model, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    message = string("Cannot delete variable as it is constrained with other",
                     " variables in a `MOI.VectorOfVariables`.")
    err = MOI.DeleteNotAllowed(y, message)
    @test_throws err MOI.delete(model, y)
    @test MOI.get(model, MOI.ObjectiveFunction{typeof(obj)}()) ≈ 1fx + 2fy

    @test MOI.is_valid(model, c8)
    MOI.delete(model, c8)
    @test !MOI.is_valid(model, c8)
    @test 0 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.SecondOrderCone}())

    MOI.delete(model, y)
    @test MOI.get(model, MOI.ObjectiveFunction{typeof(obj)}()) ≈ 1fx

    f = MOI.get(model, MOI.ConstraintFunction(), c2)
    @test f.affine_terms == MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([3, 1], [x, x]))
    @test f.quadratic_terms == MOI.VectorQuadraticTerm.([1], MOI.ScalarQuadraticTerm.([1], [x], [x]))
    @test f.constants == [7, 3, 4]

    f = MOI.get(model, MOI.ConstraintFunction(), c6)
    @test f.terms == MOI.VectorAffineTerm.([1], MOI.ScalarAffineTerm.([2], [x]))
    @test f.constants == [6, 8]

    @test 1 == @inferred MOI.get(model, MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.Nonpositives}())
    @test [c7] == @inferred MOI.get(model, MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonpositives}())

    f =  MOI.get(model, MOI.ConstraintFunction(), c7)
    @test f.variables == [x]

    s =  MOI.get(model, MOI.ConstraintSet(), c7)
    @test MOI.dimension(s) == 1
end

# We create a new function and set to test catching errors if users create their
# own sets and functions
struct FunctionNotSupportedBySolvers <: MOI.AbstractScalarFunction
    x::MOI.VariableIndex
end
struct SetNotSupportedBySolvers <: MOI.AbstractSet end
@testset "Default fallbacks" begin
    @testset "set" begin
        model = MOIU.Model{Float64}()
        x = MOI.add_variable(model)
        func = convert(MOI.ScalarAffineFunction{Float64}, MOI.SingleVariable(x))
        c = MOI.add_constraint(model, func, MOI.LessThan(0.0))
        @test !MOI.supports_constraint(model, FunctionNotSupportedBySolvers, SetNotSupportedBySolvers)

        # set of different type
        @test_throws Exception MOI.set(model, MOI.ConstraintSet(), c, MOI.GreaterThan(0.0))
        # set not implemented
        @test_throws Exception MOI.set(model, MOI.ConstraintSet(), c, SetNotSupportedBySolvers())

        # function of different type
        @test_throws Exception MOI.set(model, MOI.ConstraintFunction(), c, MOI.VectorOfVariables([x]))
        # function not implemented
        @test_throws Exception MOI.set(model, MOI.ConstraintFunction(), c, FunctionNotSupportedBySolvers(x))
    end
end

@testset "ListOfSupportedConstraints" begin
    @testset "$set" for set in (
        MOI.EqualTo(1.0), MOI.GreaterThan(1.0), MOI.LessThan(1.0),
        MOI.Interval(1.0, 2.0), MOI.Semicontinuous(1.0, 2.0),
        MOI.Semiinteger(1.0, 2.0), MOI.Integer(), MOI.ZeroOne()
    )
        model = MOIU.Model{Float64}()
        x = MOI.add_variable(model)
        MOI.add_constraint(model, MOI.SingleVariable(x), set)
        @test MOI.get(model, MOI.ListOfConstraints()) ==
            [(MOI.SingleVariable, typeof(set))]
    end
end

@testset "Extension dictionary" begin
    model = MOIU.Model{Float64}()
    model.ext[:my_store] = 1
    @test model.ext[:my_store] == 1
    MOI.empty!(model)
    @test model.ext[:my_store] == 1
    model.ext[:my_store] = 2
    dest = MOIU.Model{Float64}()
    MOI.copy_to(dest, model)
    @test !haskey(dest.ext, :my_store)
    @test model.ext[:my_store] == 2
end


max_stars_repo_name: imweixin/Ostriv

max_issues_repo_name: imweixin/Ostriv

max_forks_repo_name: imweixin/Ostriv


text:
module Xct

export sctenifoldxct

const NCOMP1,NCOMP2=3,5
const NLAYERS,NCELLS=10,500

using Statistics, LinearAlgebra, Distributions, Random
using MultipleTesting, Random, SparseArrays, TSVD, TensorToolbox
# import KrylovKit

include("support.jl")

function sctenifoldxct(X::AbstractMatrix{T}, Y::AbstractMatrix{T}; donorm::Bool=false) where T<:Real

end

end

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
facts("Images") do
  context("type creation") do
    pic = image("this"; volume="/a", volume="/b")
    @fact pic.base.name => "ubuntu:14.04"
    @fact pic.volumes => ["/a", "/b"]

    for name in names(RudeOil.BuildImage)
      if name ∉ [:name, :base, :volumes, :ppas]
        @fact getfield(pic, name) => RudeOil.DEFAULT_IMAGE[name]
      end
    end
  end

  context("dockerfile creation") do
    context("empty") do
      pic = RudeOil.BuildImage("that")
      @fact dockerfile(pic) => ismatch("FROM\\s+$(pic.base.name)")
      @fact dockerfile(pic) => ismatch("RUN\\s+pip") |> not
      @fact dockerfile(pic) => ismatch("RUN\\s+apt-get install") |> not
      @fact dockerfile(pic) => ismatch("RUN\\s+mkdir") |> not
    end
    context("packages") do
      pic = image("that", packages=["pack", "age"])
      match = "RUN\\s+apt-get\\s+install\\s+-y\\s+pack\\s+age"
      @fact dockerfile(pic) => ismatch(match)
    end
    context("pips") do
      pic = image("that", pips=["pack", "age"])
      @fact dockerfile(pic) => ismatch("RUN pip install pack age")
    end
    context("volumes") do
      pic = image("that", volume="/a", volume="/b")
      @fact dockerfile(pic) => ismatch("RUN mkdir -p /a /b")
      @fact dockerfile(pic) => ismatch("VOLUME /a /b")
    end
  end

  context("docker image creation") do
    first = image("mynewimage"; workdir="/myvol", volume="/myvol", packages=["julia"])
    context("single image") do
      activate(machine) do vm
        vm |> first |> run
        images = split(readchomp(RudeOil.command(vm, "images", first.name)), '\n')[2:end]
        images = [split(u, ' ')[1] for u in images]
        @fact images => contains(first.name)
      end
    end
    context("chained image") do
      second = image("mysecondimage"; packages=["python"])
      activate(machine) do vm
        vm |> first |> second |> run
        images = split(readchomp(RudeOil.command(vm, "images", "")), '\n')[2:end]
        images = [split(u, ' ')[1] for u in images]
        @fact images => contains(first.name)
        @fact images => contains(second.name)
      end
    end
  end

  context("Upstream docker image") do
    pic = image("ubuntu:14.04")
    @fact names(pic) => [:name]
  end
end


text:
using Printf

function print_header(data::Data)
    if data.verbosity == 1
        @printf("Iter \t  Objective  \t Inf Linear \t Radius \t Gradient res \t D[end, end] \t AC\n")
    elseif data.verbosity == 2
        @printf("Iter \t  Objective  \t Inf Linear \t Radius \t Gradient res \t D[end, end] \t AC\t\t ||LDL'-Z'PZ||   ||Q*Q'-I|| \t ||Q1*R1-A1||\n")
    end
end

function print_info(data::Data)
    if data.verbosity == 1
        @printf("%d \t  %.5e \t %.5e \t %.5e \t %s \t %.5e \t %d/%d   \n",
            data.iteration,
            0.5*data.x'*data.F.P*data.x + dot(data.q, data.x),
            max(maximum(data.A*data.x - data.b), 0),
            norm(data.x),
            isnan(data.residual) ? string(@sprintf("%.5e", data.residual), "   ") : @sprintf("%.5e", data.residual),
            length(data.F.D) > 0 ? data.F.D[end] : NaN,
            length(data.working_set),
            data.F.artificial_constraints
        )
    elseif data.verbosity == 2
        @printf("%d \t  %.5e \t %.5e \t %.5e \t %s \t %.5e \t %d/%d    \t %.5e \t %.5e \t %.5e\n",
            data.iteration,
            0.5*data.x'*data.F.P*data.x + dot(data.q, data.x),
            max(maximum(data.A*data.x - data.b), 0),
            norm(data.x),
            isnan(data.residual) ? string(@sprintf("%.5e", data.residual), "   ") : @sprintf("%.5e", data.residual),
            length(data.F.D) > 0 ? data.F.D[end] : NaN,
            length(data.working_set),
            data.F.artificial_constraints,
            norm(data.F.U'*data.F.D*data.F.U - data.F.Z'*data.F.P*data.F.Z),
            norm(data.F.QR.Q'*data.F.QR.Q - I),
            norm(data.F.QR.Q1*data.F.QR.R1 - data.A[data.working_set, :]'),
        )
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ f39ca040-710d-11eb-102e-c5521f152a14
begin
	using DataStructures
	using Logging
	const MAXLOG=100
end

# ╔═╡ cf5232c8-710c-11eb-00bb-7d9bcefd8660
md"""
# Uninformed search
A problem consists of five parts: the **initial state**, a set of **actions**, a **transition model** describing the results of those actions, a **goal test function**, and a **path cost function**. The environment of the problem is represented by a state space. A path through the state space from the initial state to a goal state is a solution

Determine the problem formulation for the following:

1. You start with the sequence EACEAABAAB, or in general any sequence made from A, B, C, and E. You can transform this sequence using the following equalities: AC = E, AB = BC, BB = E, and Ex = x for any x. For example, ABBC can be transformed into AEC, and then AC, and then E. Your goal is to produce the sequence E.
    * initial state: start sequence e.g. "ABABAECCEC"
    * actions: transform a pair
    * transition model: AC = E, AB = BC, BB = E & Ex = x for any x
    * goal test function: sequence == "E"
    * path cost function: number of iterations
2. There are six glass boxes in a row, each with a lock. Each of the first five boxes holds a key unlocking the next box in line; the last box holds a banana. You have the key to the first box, and you want the banana.
    * initial state: closed boxes
    * actions: open box, obtain key
    * transition model: next box can be opened with newly obtained key
    * goal test function: obtain banana
    * path cost function: number of actions

For the first problem, implement the breadth-first search algorithm and test it for different start sequences.

*Tip:* given the fact that breadth-first is based on a FIFO-queue, you might want to take a look at `DataStructures.jl`
"""

# ╔═╡ 5e9aaef6-75e0-11eb-209d-fb7e9c6f00a5
# to get you started:
let
	comb = Dict(  "AC"=>"E",
                    "AB"=>"BC",
                    "BB"=>"E",
                    "EA"=>"A",
                    "EB"=>"B",
                    "EC"=>"C",
                    "EE"=>"E")

    struct Node end;
    struct Problem end;

	function goaltest_sequence() end;
	function findactions_sequence() end;
	function solution_sequence() end;
	function treesearch(p::Problem; goaltest::Function, findactions::Function, solution::Function)
		rootnode = p.rootnode();
		fringe = DataStructures.Queue{p.Node}()       # create FIFO Queue
		enqueue!(fringe, rootnode)                  # place root node in the Queue
		while !isempty(fringe)
			node = dequeue!(fringe)
			for action in findactions(node)
				child = Node()
				if child ∉ fringe
					if goaltest(p, child)
						return solution()
					end
					enqueue!(fringe, child)
				end
			end
		end
		@warn "failed to find a solution"
	end;

end

# ╔═╡ Cell order:
# ╠═f39ca040-710d-11eb-102e-c5521f152a14
# ╟─cf5232c8-710c-11eb-00bb-7d9bcefd8660
# ╠═5e9aaef6-75e0-11eb-209d-fb7e9c6f00a5

max_stars_repo_path: PS/PS02 - uninformed search .jl

max_issues_repo_path: PS/PS02 - uninformed search .jl

max_forks_repo_path: PS/PS02 - uninformed search .jl


text:
module ArbFloats

import Base: hash, convert, promote_rule, isa,
    string, show, showcompact, showall, parse,
    finalizer, decompose, precision, setprecision,
    copy, deepcopy,
    zero, one, isinteger, ldexp, frexp, eps,
    isequal, isless, (==),(!=),(<),(<=),(>=),(>), contains,
    min, max, minmax,
    isnan, isinf, isfinite, issubnormal,
    signbit, sign, flipsign, copysign, abs, inv,
    (+),(-),(*),(/),(\),(%),(^), sqrt, hypot,
    trunc, round, ceil, floor,
    fld, cld, div, mod, rem, divrem, fldmod,
    muladd, fma,
    exp, expm1, log, log1p, log2, log10,
    sin, cos, tan, csc, sec, cot, asin, acos, atan, atan2,
    sinh, cosh, tanh, csch, sech, coth, asinh, acosh, atanh,
    sinc, gamma, lgamma, digamma, zeta, factorial,
    BigInt, BigFloat,
    Cint

export ArbFloat,      # co-matched decimal rounding, n | round(hi,n,10) == round(lo,n,10)
       ArbFloat512, ArbFloat256, ArbFloat128, ArbFloat64, ArbFloat32, ArbFloat16,
       @ArbFloat,     # converts string form of argument, precision is optional first arg in two arg form
       midpoint, radius, upperbound, lowerbound, bounds,
       stringCompact, stringAll, stringAllCompact,
       smartvalue, smartstring, showsmart, showallcompact,
       two, three, four, copymidpoint, copyradius, deepcopyradius,
       epsilon, trim, decompose, isexact, notexact,
       iszero, notzero, nonzero, isone, notone, notinteger,
       ispositive, notpositive, isnegative, notnegative,
       includesAnInteger, excludesIntegers, includesZero, excludesZero,
       includesPositive, excludesPositive, includesNegative, excludesNegative,
       includesNonpositive,  includesNonnegative,
       notequal, approxeq, ≊,
       overlap, donotoverlap, iscontainedby, doesnotcontain, isnotcontainedby,
       invsqrt, pow, root, tanpi, cotpi, logbase, sincos, sincospi, sinhcosh,
       doublefactorial, risingfactorial, rgamma, agm, polylog,
       relativeError, relativeAccuracy, midpointPrecision, trimmedAccuracy,
       PI,SQRTPI,LOG2,LOG10,EXP1,EULER,CATALAN,KHINCHIN,GLAISHER,APERY, # constants
       MullerKahanChallenge

# ensure the requisite libraries are available

isdir(Pkg.dir("Nemo")) || throw(ErrorException("Nemo not found"))

libDir = Pkg.dir("Nemo/local/lib");
libFiles = readdir(libDir);
libarb   = joinpath(libDir,libFiles[findfirst([startswith(x,"libarb") for x in libFiles])])
libflint = joinpath(libDir,libFiles[findfirst([startswith(x,"libflint") for x in libFiles])])
isfile(libarb)   || throw(ErrorException("libarb not found"))
isfile(libflint) || throw(ErrorException("libflint not found"))

@static if is_linux() || is_bsd() || is_unix()
    libarb = String(split(libarb,".so")[1])
    libflint = String(split(libflint,".so")[1])
end
@static if is_apple()
    libarb = String(split(libarb,".dynlib")[1])
    libflint = String(split(libflint,".dynlib")[1])
end
@static if is_windows()
    libarb = String(split(libarb,".dll")[1])
    libflint = String(split(libflint,".dll")[1])
end

macro libarb(sym)
    (:($sym), libarb)
end

macro libflint(sym)
    (:($sym), libflint)
end

NotImplemented(info::AbstractString="") = error(string("this is not implemented\n\t",info,"\n"))

include("type/ArfFloat.jl")
include("type/ArbFloat.jl")

include("type/primitive.jl")
include("type/predicates.jl")
include("type/string.jl")
include("type/convert.jl")
include("type/compare.jl")
include("type/io.jl")

include("math/arith.jl")
include("math/round.jl")
include("math/elementary.jl")
include("math/constants.jl")
include("math/special.jl")

end # ArbFloats


text:
import NLPModels: Counters

"""
NLPAtX contains the important information concerning a non linear problem at
the iteration x. Basic information is:
 - x the current candidate for solution to our original problem
 - f(x) which is the funciton evaluation at x
 - g(x) which is the gradient evaluation at x
 - Hx which is the hessian representation at x

 - mu : Lagrange multiplier of the bounds constraints

 - cx : evaluation of the constraint function at x
 - Jx : jacobian matrix of the constraint function at x
 - lambda : Lagrange multiplier of the constraints

 - start_time : Default is a NaN, can be updated to fit the start of the algorithm.
 - evals : number of evaluations of the function (import the type NLPModels.Counters)

 All these information (except for x) are optionnal and need to be update when
 required. The update is done trhough the update! function.
"""
mutable struct 	NLPAtX <: AbstractState

#Unconstrained State
    x            :: Iterate     # current point
    fx           :: FloatVoid   # objective function
    gx           :: Iterate     # gradient
    Hx           :: MatrixType  # Accurate?

#Bounds State
    mu           :: Iterate     # Lagrange multipliers with bounds

#Constrained State
    cx           :: Iterate     # vector of constraints lc <= c(x) <= uc
    Jx           :: MatrixType  # jacobian matrix
    lambda       :: Iterate     # Lagrange multipliers

 #Resources State
    start_time   :: FloatVoid
    evals        :: Counters

 function NLPAtX(x          :: Iterate,
                 lambda     :: Iterate;
                 fx         :: FloatVoid    = NaN,
                 gx         :: Iterate      = NaN * fill(1.0, size(x)),
                 Hx         :: MatrixType   = zeros(0,0),
                 mu         :: Iterate      = NaN * fill(1.0, size(x)),
                 cx         :: Iterate      = NaN * fill(1, size(lambda)),
                 Jx         :: MatrixType   = zeros(length(lambda),length(x)),
                 start_time :: FloatVoid    = NaN,
                 evals      :: Counters     = Counters())

  return new(x, fx, gx, Hx, mu, cx, Jx, lambda, start_time, evals)
 end
end

"""
An additional constructor for unconstrained problems
"""
function NLPAtX(x          :: Iterate;
                fx         :: FloatVoid    = NaN,
                gx         :: Iterate      = NaN * fill(1.0, size(x)),
                Hx         :: MatrixType   = zeros(0,0),
                start_time :: FloatVoid    = NaN,
                evals      :: Counters     = Counters())

	return NLPAtX(x, zeros(0), fx = fx, gx = gx,
                  Hx = Hx, start_time = start_time, evals = evals)
end

"""
Updates the (desired) values of an object of type NLPAtX.
Inputs:
 - An NLPAtX object
 - Any keywords that needs to be updated.
"""
function update!(nlpatx :: NLPAtX;
                 x      :: Iterate    = nothing,
                 fx     :: FloatVoid  = nothing,
                 gx     :: Iterate    = nothing,
                 Hx     :: MatrixType = nothing,
                 mu     :: Iterate    = nothing,
                 cx     :: Iterate    = nothing,
                 Jx     :: MatrixType = nothing,
                 lambda :: Iterate    = nothing,
                 tmps   :: FloatVoid  = nothing,
                 evals  :: Union{Counters, Nothing}  = nothing)

    nlpatx.x   = x  == nothing  ? nlpatx.x   : x
    nlpatx.fx  = fx == nothing  ? nlpatx.fx  : fx
    nlpatx.gx  = gx == nothing  ? nlpatx.gx  : gx
    nlpatx.Hx  = Hx == nothing  ? nlpatx.Hx  : Hx
    nlpatx.mu  = mu == nothing  ? nlpatx.mu  : mu
    nlpatx.cx  = cx == nothing  ? nlpatx.cx  : cx
    nlpatx.Jx  = Jx == nothing  ? nlpatx.Jx  : Jx

    nlpatx.lambda     = lambda == nothing  ? nlpatx.lambda    : lambda

    nlpatx.start_time = tmps   == nothing ? nlpatx.start_time : tmps
    nlpatx.evals      = evals  == nothing ? nlpatx.evals      : evals

    return nlpatx
end

# function convert_nlp(T,  nlpatx :: NLPAtX)
#
#     nlpatxT         = NLPAtX(zeros(T, length(nlpatx.x)))
#     nlpatxT.x       = typeof(nlpatx.x)      != Nothing ? convert.(T, nlpatx.x)      : nlpatx.x
#     nlpatxT.fx      = typeof(nlpatx.fx)     != Nothing ? convert.(T, nlpatx.fx)     : nlpatx.fx
#     nlpatxT.gx      = typeof(nlpatx.gx)     != Nothing ? convert.(T, nlpatx.gx)     : nlpatx.gx
#     nlpatxT.Hx      = typeof(nlpatx.Hx)     != Nothing ? convert.(T, nlpatx.Hx)     : nlpatx.Hx
#     nlpatxT.mu      = typeof(nlpatx.mu)     != Nothing ? convert.(T, nlpatx.mu)     : nlpatx.mu
#     nlpatxT.cx      = typeof(nlpatx.cx)     != Nothing ? convert.(T, nlpatx.cx)     : nlpatx.cx
#     nlpatxT.Jx      = typeof(nlpatx.Jx)     != Nothing ? convert.(T, nlpatx.Jx)     : nlpatx.Jx
#     nlpatxT.lambda  = typeof(nlpatx.lambda) != Nothing ? convert.(T, nlpatx.lambda) : nlpatx.lambda
#
#     return nlpatxT
# end


text:
using DataFrames, CSV
data = CSV.read("data/purchaseData.csv", DataFrame, copycols = true)

println(size(data),"\n")
println(names(data),"\n")
println(first(data, 6),"\n")
println(describe(data),"\n")

max_stars_repo_name: newptcai/StatsWithJuliaBook

max_issues_repo_name: newptcai/StatsWithJuliaBook

max_forks_repo_name: newptcai/StatsWithJuliaBook


text:
abstract type AbstractResult end

abstract type SpectralResult <: AbstractResult end

const Projection{T <: Real} = AbstractMatrix{T}


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module VirtualNetworkClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_DefaultApi.jl")
include("api_SubnetsApi.jl")
include("api_VirtualNetworkPeeringsApi.jl")
include("api_VirtualNetworksApi.jl")

# export models
export convert, AddressSpace
export convert, DhcpOptions
export convert, IPAddressAvailabilityResult
export convert, Resource
export convert, ResourceNavigationLinkFormat
export convert, SubResource
export convert, SubnetListResult
export convert, SubnetPropertiesFormat
export convert, VirtualNetworkListResult
export convert, VirtualNetworkListUsageResult
export convert, VirtualNetworkPeeringListResult
export convert, VirtualNetworkPeeringPropertiesFormat
export convert, VirtualNetworkPropertiesFormat
export convert, VirtualNetworkUsage
export convert, VirtualNetworkUsageName
export convert, IPConfiguration
export convert, NetworkSecurityGroup
export convert, ResourceNavigationLink
export convert, RouteTable
export convert, Subnet
export convert, VirtualNetwork
export convert, VirtualNetworkPeering

# export operations
export convert, DefaultApi, SubnetsApi, VirtualNetworkPeeringsApi, VirtualNetworksApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
using TomoForward
using Images
using Plots
using XfromProjections
using StaticArrays
using Logging
using ImageTransformations
include("../../phantoms/simple_phantoms.jl")

replace_nan(v) = map(x -> isnan(x) ? zero(x) : x, v)

#Translation
function transformation(ϕ)
    v1(x) = x
    v2(x) = x - ϕ*10

    t(x) = SArray{Tuple{2},Float64,1,2}(v1(x[1]),v2(x[2]))
    return t
end

H, W = 50,50#128, 128
img = zeros(H,W)
img = disc_phantom(img,15,25,8)
#img = imresize(img, H, W)
#p_img = plot(Gray.(img), aspect_ratio=:equal, framestyle=:none, title="Original")

function radon_operator(img)
    nangles = 2
    detcount = Int(floor(size(img,1)*1.4))
    angles = rand(0.0:0.001:π, nangles)
    proj_geom = ProjGeom(1.0, detcount, angles)
    A = fp_op_parallel2d_line(proj_geom, size(img, 1), size(img, 2))
    return A
end

frames = zeros(H,W,20)
map(t -> frames[:,:,t+1]=replace_nan(warp(img, transformation(t*0.1), axes(img))), 0:size(frames)[3]-1)
As = map(t -> radon_operator(frames[:,:,t]),1:size(frames)[3])
bs = zeros(size(As[1])[1],size(frames)[3])
map(t -> bs[:,t] = As[t]*vec(frames[:,:,t]), 1:size(frames)[3])
niter=500
u0s = zeros(H,W,size(frames)[3])

mask = disc_phantom(u0s[:,:,1], 25.5, 25.5, 23)
mask = repeat(mask, outer=[1,1,20])

w_tv = 0.1
w_flow  = 0.5

@info "Reconstructing using joint motion estimation and reconstruction"
c=1.1
us_flow = recon2d_tv_primaldual_flow(As, bs, u0s, 20, niter, w_tv, w_flow, c, mask, 0.012)

@info "Reconstruction using tv regularization frame by frame"
us_tv = zeros(H,W,size(frames)[3])
for t = 1:size(frames)[3]
    A = As[t]
    p = bs[:,t]
    us_tv[:,:,t] .= recon2d_tv_primaldual!(us_tv[:,:,t], A, p, niter, w_tv, c)
end

@info "Preparing results in human readable format"
anim = @animate for t=1:size(frames)[3]
    l = @layout [a b c]
    p1 = plot(Gray.(frames[:,:,t]), aspect_ratio=:equal, framestyle=:none, title="Ground truth")
    p2 = plot(Gray.(us_flow[:,:,t]), aspect_ratio=:equal, framestyle=:none, title="Flow")
    p3 = plot(Gray.(us_tv[:,:,t]), aspect_ratio=:equal, framestyle=:none, title="TV")
    plot(p1, p2, p3, layout = l)
end

cwd = @__DIR__
path = normpath(joinpath(@__DIR__, "result"))
!isdir(path) && mkdir(path)
cd(path)
gif(anim, "reconstruction_flow.gif", fps = 1)
cd(cwd)

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*-----------------------------------------------------------------------
#* Author: Tiago Os򱨯rio
#* $Id$
#*---------------------------------------------------------------------
type Ratio
	Ratio()=begin
		new(
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"input signal 1"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"input signal 2"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"output signal"
			)),
			[
				:(input2*output=input1),
				:(output=input1*1e5),
				:(output=-input1*1e5),
			],
			[
				"Calculate output","Calculate output","Calculate output",
			],
			[:input1,:input2,:output,]
		)
	end
	input1::DanaReal
	input2::DanaReal
	output::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Ratio
function setEquationFlow(in::Ratio)
	if abs(input2)>1e-5
		addEquation(1)
	else
		if input2>0
			addEquation(2)
		else
			addEquation(3)
		end
	end
end
function atributes(in::Ratio,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/Ratio"
	fields[:Brief]="Model Ratio."
	fields[:Info]="== Inputs ==
* Two input signals.

== Outputs ==
* One output signal.
"
	drive!(fields,_)
	return fields
end
Ratio(_::Dict{Symbol,Any})=begin
	newModel=Ratio()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
## Тест проверяет корректность решения прямой задачи. Алгоритм описан в /docs/src/direct/direct_check.md
## Возвращает решение, аналитическое решение, сетку по X, сетку по T.
using NonLinearReactionAdvectionDiffusionWithFrontData;
using NonLinearReactionAdvectionDiffusionWithFrontData: phidetermination, Φ;
using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;
using ForwardDiff;

## Зададим параметры для прямой задачи
u_l(t) = 0;                     # ГУ удовлетворяющие модельной функции
u_r(t) = 0;                     # ГУ удовлетворяющие модельной функции
qf(x) = 4*sin(3 * π * x);       # Коэффициент линейного усиления
ε = 0.2;                        # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0, 1;                   # Область по T
N, M = 50, 80;                  # Кол-во разбиений по X, T
h = (b-a)/N;                    # шаг по X
τ = (T-t₀)/M;                   # шаг по T
Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
qₙ =      qf.(Xₙ);              # Сеточные значения коэффициента линейного усиления
ulₘ= u_l.(Tₘ);                  # Сеточные значения левого  ГУ
urₘ= u_r.(Tₘ);                  # Сеточные значения правого ГУ

## Зададим модельную функцию и невязку, получаемую после подстановки `g` в исходное уравнение
function g(x, m)
    t = Tₘ[m]
    (1 - 2t)*sin(π*x)
end
function g_d(x::Real, m::Int)
    t = Tₘ[m];
    - ε * π^2 * (1 - 2t) * sin(π * x) + π * (1 - 2t)^2 * sin(π * x) * cos(π * x) - qf(x) * (1 -2t) * sin(π * x) + 2 * sin(π * x)
end

y₀ = g.(Xₙ, 1);               # Начальные условия
## Модельное решение найденное с помощью известного аналитического решения
u_model = [ g(x, m) for x in Xₙ, m in 1:M+1];

## Создадим функцию, которая будет вычислять вектор правой части с добавлением невязки
function RP(y, m, Xₙ, N, ε, ulₘ, urₘ, qₙ)
    d = [ g_d(x, m) for x in Xₙ[2:N] ]
    NonLinearReactionAdvectionDiffusionWithFrontData.directRP(y, m, Xₙ, N, ε, ulₘ, urₘ, qₙ) - d
end
## Хоть мы и конструируем якобиан с помощью автоматического дифференцирования,
## примите во внимание, что Якобиан ``f_y`` при добавлении `g_d` останется
## без изменений, т.к. `g_d` зависит только от ``x,t``.
## То, что он не зависит от добавления `g_d` мы убедимся через 7 строк.
j(y, m, Xₙ, N, ε, ulₘ, urₘ, qₙ) = ForwardDiff.jacobian( z -> RP(z, m, Xₙ, N, ε, ulₘ, urₘ, qₙ), y)

## С использованием автоматического дифференцирования
u, XX, TP = solve(y₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ, RP, j);
@test all(isapprox.(u_model, u, atol = 0.01))

## С использованием трехдиагонального якобиана без изменений.
u, XX, TP = solve(y₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ, RP, NonLinearReactionAdvectionDiffusionWithFrontData.∂DRP_∂y);
@test all(isapprox.(u_model, u, atol = 0.01))

(u, u_model, Xₙ, Tₘ)
nothing #hide


text:
@io struct TKey32
    fNbytes::Int32
    fVersion::Int16
    fObjlen::Int32
    fDatime::UInt32
    fKeylen::Int16
    fCycle::Int16
    fSeekKey::Int32
    fSeekPdir::Int32
    fClassName::String
    fName::String
    fTitle::String
end

@io struct TKey64
    fNbytes::Int32
    fVersion::Int16
    fObjlen::Int32
    fDatime::UInt32
    fKeylen::Int16
    fCycle::Int16
    fSeekKey::Int64
    fSeekPdir::Int64
    fClassName::String
    fName::String
    fTitle::String
end

const TKey = Union{TKey32, TKey64}

function unpack(io, ::Type{TKey})
    start = position(io)
    skip(io, 4)
    fVersion = readtype(io, Int16)
    seek(io, start)
    if fVersion <= 1000
        return unpack(io, TKey32)
    end
    unpack(io, TKey64)
end

@with_kw struct TBasketKey
    fNbytes::Int32
    fVersion::Int16
    fObjlen::Int32
    fDatime::UInt32
    fKeylen::Int16
    fCycle::Int16
    fSeekKey::Int64
    fSeekPdir::Int64
    fClassName::String
    fName::String
    fTitle::String
    fBufferSize::Int32
    fNevBufSize::Int32
    fNevBuf::Int32
    fLast::Int32
end

function unpack(io, T::Type{TBasketKey})
    start = position(io)
    fields = Dict{Symbol, Union{Integer, String}}()
    fields[:fNbytes] = readtype(io, Int32)
    fields[:fVersion] = readtype(io, Int16)  # FIXME if "complete" it's UInt16 (acc. uproot)

    inttype = fields[:fVersion] <= 1000 ? Int32 : Int64

    fields[:fObjlen] = readtype(io, Int32)
    fields[:fDatime] = readtype(io, UInt32)
    fields[:fKeylen] = readtype(io, Int16)
    fields[:fCycle] = readtype(io, Int16)
    fields[:fSeekKey] = readtype(io, inttype)
    fields[:fSeekPdir] = readtype(io, inttype)
    fields[:fClassName] = readtype(io, String)
    fields[:fName] = readtype(io, String)
    fields[:fTitle] = readtype(io, String)

    # if complete (which is true for compressed, it seems?)
    seek(io, start + fields[:fKeylen] - 18 - 1)
    fields[:fVersion] = readtype(io, Int16)  # FIXME if "complete" it's UInt16 (acc. uproot)
    fields[:fBufferSize] = readtype(io, Int32)
    fields[:fNevBufSize] = readtype(io, Int32)
    fields[:fNevBuf] = readtype(io, Int32)
    fields[:fLast] = readtype(io, Int32)

    T(; fields...)
end

iscompressed(t::T) where T<:Union{TKey, TBasketKey} = t.fObjlen != t.fNbytes - t.fKeylen
origin(t::T) where T<:Union{TKey, TBasketKey} = iscompressed(t) ? -t.fKeylen : t.fSeekKey
seekstart(io, t::T) where T<:Union{TKey, TBasketKey} = seek(io, t.fSeekKey + t.fKeylen)

datastream(io, tkey::TKey) = IOBuffer(decompress_datastreambytes(compressed_datastream(io, tkey), tkey))

"""
    compressed_datastream(io, tkey)

Extract all [compressionheader][rawbytes] from a `TKey`. This is an isolated function
because we want to compartmentalize disk I/O as much as possible.

See also: [`decompress_datastreambytes`](@ref)
"""
function compressed_datastream(io, tkey)
    if !iscompressed(tkey)
        @debug ("Uncompressed datastream of $(tkey.fObjlen) bytes " *
                "at $start (TKey '$(tkey.fName)' ($(tkey.fClassName)))")
        skip(io, 1)   # ???
        return read(io, tkey.fObjlen)
    end
    seekstart(io, tkey)
    return read(io, tkey.fNbytes - tkey.fKeylen)
end

"""
    decompress_datastreambytes(compbytes, tkey)

Process the compressed bytes `compbytes` which was read out by `compressed_datastream` and
pointed to from `tkey`. This function simply return uncompressed bytes according to
the compression algorithm detected (or the lack of).
"""
function decompress_datastreambytes(compbytes, tkey)
    # not compressed
    iscompressed(tkey) || return compbytes

    # compressed
    io = IOBuffer(compbytes)
    fufilled = 0
    uncomp_data = Vector{UInt8}(undef, tkey.fObjlen)
    while fufilled < tkey.fObjlen # careful with 0/1-based index when thinking about offsets
        compression_header = unpack(io, CompressionHeader)
        cname, _, compbytes, uncompbytes = unpack(compression_header)
        rawbytes = read(io, compbytes)

        # indexing `0+1 to 0+2` are two bytes, no need to +1 in the second term
        @view(uncomp_data[fufilled+1:fufilled+uncompbytes]) .= if cname == "ZL"
            transcode(ZlibDecompressor, rawbytes)
        elseif cname == "XZ"
            transcode(XzDecompressor, rawbytes)
        elseif cname == "ZS"
            transcode(ZstdDecompressor, rawbytes)
        elseif cname == "L4"
            # skip checksum which is 8 bytes
            lz4_decompress(rawbytes[9:end], uncompbytes)
        else
            error("Unsupported compression type '$(String(compression_header.algo))'")
        end
        fufilled += uncompbytes
    end
    return uncomp_data
end
@io struct FilePreamble
    identifier::SVector{4, UInt8}  # Root file identifier ("root")
    fVersion::Int32                # File format version
end

# https://root.cern/doc/v624/RMiniFile_8cxx_source.html#l00239
@io struct FileHeader32
    fBEGIN::Int32                  # Pointer to first data record
    fEND::UInt32                   # Pointer to first free word at the EOF
    fSeekFree::UInt32              # Pointer to FREE data record
    fNbytesFree::Int32             # Number of bytes in FREE data record
    nfree::Int32                   # Number of free data records
    fNbytesName::Int32             # Number of bytes in TNamed at creation time
    fUnits::UInt8                  # Number of bytes for file pointers
    fCompress::Int32               # Compression level and algorithm
    fSeekInfo::UInt32              # Pointer to TStreamerInfo record
    fNbytesInfo::Int32             # Number of bytes in TStreamerInfo record
    fUUID::SVector{18, UInt8}      # Universal Unique ID
end

@io struct FileHeader64
    fBEGIN::Int32                  # Pointer to first data record
    fEND::UInt64                   # Pointer to first free word at the EOF
    fSeekFree::UInt64              # Pointer to FREE data record
    fNbytesFree::Int32             # Number of bytes in FREE data record
    nfree::Int32                   # Number of free data records
    fNbytesName::Int32             # Number of bytes in TNamed at creation time
    fUnits::UInt8                  # Number of bytes for file pointers
    fCompress::Int32               # Compression level and algorithm
    fSeekInfo::UInt64              # Pointer to TStreamerInfo record
    fNbytesInfo::Int32             # Number of bytes in TStreamerInfo record
    fUUID::SVector{18, UInt8}      # Universal Unique ID
end

const FileHeader = Union{FileHeader32, FileHeader64}

@io struct ROOTDirectoryHeader32
    fVersion::Int16
    fDatimeC::UInt32
    fDatimeM::UInt32
    fNbytesKeys::Int32
    fNbytesName::Int32
    fSeekDir::Int32
    fSeekParent::Int32
    fSeekKeys::Int32
end

@io struct ROOTDirectoryHeader64
    fVersion::Int16
    fDatimeC::UInt32
    fDatimeM::UInt32
    fNbytesKeys::Int32
    fNbytesName::Int32
    fSeekDir::Int64
    fSeekParent::Int64
    fSeekKeys::Int64
end

const ROOTDirectoryHeader = Union{ROOTDirectoryHeader32, ROOTDirectoryHeader64}

function unpack(io::IOStream, ::Type{ROOTDirectoryHeader})
    fVersion = readtype(io, Int16)
    skip(io, -2)

    if fVersion <= 1000
        return unpack(io, ROOTDirectoryHeader32)
    else
        return unpack(io, ROOTDirectoryHeader64)
    end

end

@io struct CompressionHeader
    algo::SVector{2, UInt8}
    method::UInt8
    c1::UInt8
    c2::UInt8
    c3::UInt8
    u1::UInt8
    u2::UInt8
    u3::UInt8
end

# Built-in types
# function THashList end
function TRef end
function TArray end
function TArrayC end
function TArrayS end
function TArrayL end
function TArrayL64 end
function TRefArray end

function aliasfor(classname)
    if classname == "ROOT::TIOFeatures"
        return ROOT_3a3a_TIOFeatures
    else
        nothing
    end
end


text:
import PolyLog

n = 10_000_000
z_min = -5
z_max = 5
cmpl_data = (z_max - z_min)*rand(ComplexF64, n) + z_min*(1.0 + 1.0im)*ones(n)
real_data = map(real, cmpl_data)

println("Benchmarking li2::Float64")

map(PolyLog.li2, real_data)       # trigger compilation
@time map(PolyLog.li2, real_data)
@time map(PolyLog.li2, real_data)
@time map(PolyLog.li2, real_data)

println("Benchmarking li2::ComplexF64")

map(PolyLog.li2, cmpl_data)       # trigger compilation
@time map(PolyLog.li2, cmpl_data)
@time map(PolyLog.li2, cmpl_data)
@time map(PolyLog.li2, cmpl_data)

println("Benchmarking li3::ComplexF64")

map(PolyLog.li3, cmpl_data)       # trigger compilation
@time map(PolyLog.li3, cmpl_data)
@time map(PolyLog.li3, cmpl_data)
@time map(PolyLog.li3, cmpl_data)

println("Benchmarking li4::ComplexF64")

map(PolyLog.li4, cmpl_data)       # trigger compilation
@time map(PolyLog.li4, cmpl_data)
@time map(PolyLog.li4, cmpl_data)
@time map(PolyLog.li4, cmpl_data)

println("Benchmarking li5::ComplexF64")

map(PolyLog.li5, cmpl_data)       # trigger compilation
@time map(PolyLog.li5, cmpl_data)
@time map(PolyLog.li5, cmpl_data)
@time map(PolyLog.li5, cmpl_data)

println("Benchmarking li6::ComplexF64")

map(PolyLog.li6, cmpl_data)       # trigger compilation
@time map(PolyLog.li6, cmpl_data)
@time map(PolyLog.li6, cmpl_data)
@time map(PolyLog.li6, cmpl_data)


text:
module Textures

using Images

const WalPalette = (
	(0x00, 0x00, 0x00), (0x1f, 0x1f, 0x1f),
	(0x3f, 0x3f, 0x3f), (0x5b, 0x5b, 0x5b),
	(0x7b, 0x7b, 0x7b), (0x9b, 0x9b, 0x9b),
	(0xbb, 0xbb, 0xbb), (0xdb, 0xdb, 0xdb),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xcf, 0x97, 0x4b), (0xa7, 0x7b, 0x3b),
	(0xa7, 0x7b, 0x3b), (0xa7, 0x7b, 0x3b),
	(0x8b, 0x67, 0x2f), (0x8b, 0x67, 0x2f),
	(0x6f, 0x53, 0x27), (0x63, 0x4b, 0x23),
	(0x63, 0x4b, 0x23), (0x53, 0x3f, 0x1f),
	(0x4f, 0x3b, 0x1b), (0x43, 0x2b, 0x17),
	(0x33, 0x27, 0x13), (0x2b, 0x1f, 0x13),
	(0x27, 0x1b, 0x0f), (0x1f, 0x17, 0x0f),
	(0xb3, 0xc7, 0xd3), (0xb3, 0xc7, 0xd3),
	(0xbb, 0xbb, 0xbb), (0xab, 0xab, 0xab),
	(0x9b, 0x9b, 0x9b), (0x9b, 0x9b, 0x9b),
	(0x8b, 0x8b, 0x8b), (0x7b, 0x7b, 0x7b),
	(0x6b, 0x6b, 0x6b), (0x5b, 0x5b, 0x5b),
	(0x5b, 0x5b, 0x5b), (0x4b, 0x4b, 0x4b),
	(0x47, 0x3f, 0x43), (0x3b, 0x37, 0x37),
	(0x2f, 0x2f, 0x2f), (0x27, 0x27, 0x27),
	(0xff, 0xff, 0xa7), (0xeb, 0x97, 0x7f),
	(0xeb, 0x97, 0x7f), (0xcf, 0x97, 0x4b),
	(0xff, 0xff, 0xa7), (0xff, 0xff, 0x7f),
	(0xff, 0xff, 0x53), (0xcf, 0x97, 0x4b),
	(0xff, 0xff, 0x53), (0xff, 0xff, 0x53),
	(0xff, 0xff, 0x53), (0xff, 0xd7, 0x17),
	(0xeb, 0x9f, 0x27), (0xaf, 0x77, 0x1f),
	(0x77, 0x4f, 0x17), (0x43, 0x2b, 0x17),
	(0xeb, 0x97, 0x7f), (0xff, 0x93, 0x00),
	(0xef, 0x7f, 0x00), (0xe3, 0x6b, 0x00),
	(0xd3, 0x57, 0x00), (0xc7, 0x47, 0x00),
	(0xc7, 0x47, 0x00), (0xab, 0x2b, 0x00),
	(0x9b, 0x1f, 0x00), (0x8f, 0x17, 0x00),
	(0x73, 0x17, 0x0b), (0x67, 0x17, 0x07),
	(0x57, 0x13, 0x00), (0x43, 0x0f, 0x00),
	(0x33, 0x0b, 0x00), (0x23, 0x0b, 0x00),
	(0xd7, 0xbb, 0xb7), (0xeb, 0x97, 0x7f),
	(0xeb, 0x97, 0x7f), (0xcb, 0x9b, 0x93),
	(0xbf, 0x7b, 0x6f), (0xa7, 0x8b, 0x77),
	(0x8f, 0x77, 0x53), (0x8f, 0x77, 0x53),
	(0x87, 0x6b, 0x57), (0x7b, 0x5f, 0x4b),
	(0x67, 0x4f, 0x3b), (0x5f, 0x47, 0x37),
	(0x4b, 0x37, 0x2b), (0x3f, 0x2f, 0x23),
	(0x2b, 0x1f, 0x13), (0x1f, 0x17, 0x0f),
	(0xcb, 0x8b, 0x23), (0xaf, 0x77, 0x1f),
	(0x9f, 0x57, 0x33), (0x8b, 0x67, 0x2f),
	(0x63, 0x4b, 0x23), (0x4f, 0x3b, 0x1b),
	(0x33, 0x27, 0x13), (0x1f, 0x17, 0x0f),
	(0xff, 0xff, 0xa7), (0xff, 0xff, 0xd3),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xcb, 0xd7, 0xdf), (0x9f, 0xb7, 0xc3),
	(0x77, 0x7b, 0xcf), (0x5b, 0x87, 0x9b),
	(0x47, 0x77, 0x8b), (0x2f, 0x67, 0x7f),
	(0x2f, 0x67, 0x7f), (0x17, 0x53, 0x6f),
	(0x13, 0x4b, 0x67), (0x0b, 0x3f, 0x53),
	(0x07, 0x2f, 0x3f), (0x00, 0x1f, 0x2b),
	(0x00, 0x0f, 0x13), (0x00, 0x00, 0x00),
	(0xeb, 0xd3, 0xc7), (0xeb, 0x97, 0x7f),
	(0xeb, 0x97, 0x7f), (0xeb, 0x97, 0x7f),
	(0xbf, 0x7b, 0x6f), (0xc3, 0x73, 0x53),
	(0xb3, 0x5b, 0x4f), (0xb3, 0x5b, 0x4f),
	(0x9f, 0x4b, 0x3f), (0x7b, 0x47, 0x47),
	(0x63, 0x33, 0x33), (0x57, 0x2b, 0x2b),
	(0x3f, 0x1f, 0x1f), (0x27, 0x1b, 0x13),
	(0x17, 0x0f, 0x0b), (0x00, 0x00, 0x00),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xd3),
	(0xff, 0xff, 0xd3), (0xff, 0xff, 0xd3),
	(0xff, 0xff, 0xd3), (0xeb, 0xd3, 0xc7),
	(0xd7, 0xbb, 0xb7), (0xc7, 0xab, 0x9b),
	(0xc7, 0xab, 0x9b), (0x97, 0x9f, 0x7b),
	(0x87, 0x8b, 0x6b), (0x73, 0x73, 0x57),
	(0x5b, 0x5b, 0x43), (0x43, 0x43, 0x33),
	(0x2f, 0x2f, 0x23), (0x1b, 0x1b, 0x17),
	(0xeb, 0x97, 0x7f), (0xeb, 0x97, 0x7f),
	(0xeb, 0x97, 0x7f), (0xc3, 0x73, 0x53),
	(0xc3, 0x73, 0x53), (0xb3, 0x5b, 0x4f),
	(0xa7, 0x3b, 0x2b), (0xa7, 0x3b, 0x2b),
	(0x9f, 0x2f, 0x23), (0x8b, 0x27, 0x13),
	(0x6b, 0x2b, 0x1b), (0x57, 0x1f, 0x13),
	(0x43, 0x17, 0x0b), (0x2b, 0x0b, 0x00),
	(0x1b, 0x00, 0x00), (0x00, 0x00, 0x00),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xeb, 0xeb, 0xeb),
	(0xcb, 0xd7, 0xdf), (0xb3, 0xc7, 0xd3),
	(0x9f, 0xb7, 0xc3), (0x77, 0x7b, 0xcf),
	(0x77, 0x7b, 0xcf), (0x67, 0x6b, 0xb7),
	(0x5b, 0x5b, 0x9b), (0x4b, 0x4f, 0x7f),
	(0x3f, 0x3f, 0x67), (0x2f, 0x2f, 0x4b),
	(0x23, 0x1f, 0x2f), (0x17, 0x0f, 0x0b),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xd3),
	(0xff, 0xff, 0xd3), (0xff, 0xff, 0xa7),
	(0xff, 0xff, 0xa7), (0xff, 0xff, 0x7f),
	(0x9b, 0xab, 0x7b), (0x9b, 0xab, 0x7b),
	(0x87, 0x97, 0x63), (0x5f, 0xa7, 0x2f),
	(0x5f, 0x8f, 0x33), (0x5f, 0x7b, 0x33),
	(0x3f, 0x4f, 0x1b), (0x2f, 0x3b, 0x0b),
	(0x23, 0x2f, 0x07), (0x1b, 0x23, 0x00),
	(0x00, 0xff, 0x00), (0x00, 0xff, 0x00),
	(0xff, 0xff, 0x27), (0xff, 0xff, 0x53),
	(0xff, 0xff, 0x53), (0xff, 0xff, 0x53),
	(0xff, 0xff, 0x53), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xa7),
	(0xff, 0xff, 0x53), (0xff, 0xff, 0x53),
	(0xff, 0xff, 0x27), (0xff, 0xff, 0x27),
	(0xff, 0xff, 0x27), (0xff, 0xff, 0x27),
	(0xff, 0xeb, 0x1f), (0xff, 0xd7, 0x17),
	(0xff, 0xab, 0x07), (0xff, 0x93, 0x00),
	(0xff, 0x93, 0x00), (0xff, 0x93, 0x00),
	(0xff, 0x00, 0x00), (0xff, 0x00, 0x00),
	(0xff, 0x00, 0x00), (0xef, 0x00, 0x00),
	(0x9b, 0x1f, 0x00), (0x7f, 0x0f, 0x00),
	(0x5f, 0x00, 0x00), (0x2f, 0x00, 0x00),
	(0xff, 0x00, 0x00), (0x37, 0x37, 0xff),
	(0xff, 0x00, 0x00), (0x00, 0x00, 0xff),
	(0x5b, 0x5b, 0x43), (0x37, 0x37, 0x2b),
	(0x23, 0x23, 0x1b), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xa7), (0xeb, 0x97, 0x7f),
	(0xeb, 0x9f, 0x27), (0xff, 0xff, 0xff),
	(0xff, 0xff, 0xff), (0xff, 0xff, 0xff),
	(0xeb, 0xd3, 0xc7), (0x9f, 0x5b, 0x53),
)

type WAL <: Images.ImageFileType end

# FIXME: https://github.com/timholy/Images.jl/issues/10
# add_image_file_format(".wal", b"", WAL)

import Images.imread
function imread{S<:IO}(io::S, ::Type{WAL})
	# name of texture
	seek(io, 0)
	name = lowercase(rstrip(bytestring(read(io, Uint8, 32)),"\0"))

	# size
	width = read(io, Uint32)
	height = read(io, Uint32)

	# mipmap byte offsets
	offset = read(io, Int32, 4)

	# name of next texture (if animated)
	next_name = lowercase(rstrip(bytestring(read(io, Uint8, 32)),"\0"))

	flags = read(io, Uint32)
	contents = read(io, Uint32)
	value = read(io, Uint32)

	# decode image data
	seek(io, offset[1])
	pixels = Array(Uint8, 3, width, height)
	i = 1
	while i < length(pixels)
		c = read(io, Uint8)+1
		pixels[i]   = WalPalette[c][1]
		pixels[i+1] = WalPalette[c][2]
		pixels[i+2] = WalPalette[c][3]
		i += 3
	end

	prop = {"colorspace" => "RGB"}
	Image(pixels, prop)
end

end