text:
baremodule FoldsChainRules

module Internal

import Folds
using ChainRulesCore: #
    @thunk,
    ChainRulesCore,
    HasReverseMode,
    InplaceableThunk,
    NoTangent,
    ProjectTo,
    RuleConfig,
    add!!,
    rrule_via_ad,
    unthunk
using Core: Typeof
using Referenceables: referenceable
using Transducers: #
    BottomRF,
    Completing,
    Executor,
    R_,
    Reduced,
    Transducer,
    Transducers,
    __foldl__,
    complete,
    foldl_nocomplete,
    inner,
    next,
    reduced,
    unreduced,
    xform

include("utils.jl")
include("foldl.jl")
include("sum.jl")

end  # module Internal

end  # baremodule FoldsChainRules


text:
##############################################################################
##
## Cluster Covariance
##
##############################################################################

struct ClusterCovariance <: CovarianceEstimator
    clusternames::Tuple{Symbol, Vararg{Symbol}} # names are not allowed to be empty
    clusters::Union{Tuple{Vararg{GroupedArray}}, Nothing}
end

"""
    cluster(names::Symbol...)

Estimate variance-covariance matrix with a cluster-robust estimator.

# Arguments
- `names::Symbol...`: column names of variables that represent the clusters.
"""
cluster(ns::Symbol...) = ClusterCovariance((ns...,), nothing)
cluster(ns::NTuple{N, Symbol}, gs::NTuple{N, GroupedArray}) where N =
    ClusterCovariance(ns, gs)

"""
    names(vce::ClusterCovariance)

Return column names of variables used to form clusters for `vce`.
"""
names(v::ClusterCovariance) = v.clusternames

length(v::ClusterCovariance) = length(names(v))

show(io::IO, ::ClusterCovariance) = print(io, "Cluster-robust covariance estimator")

function show(io::IO, ::MIME"text/plain", v::ClusterCovariance)
    print(io, length(v), "-way cluster-robust covariance estimator:")
    for n in names(v)
        print(io, "\n  ", n)
    end
end

function Vcov.completecases(table, v::ClusterCovariance)
    Tables.istable(table) || throw(ArgumentError("completecases requires a table input"))
    N = length(Tables.rows(table))
    out = trues(N)
    # See https://github.com/JuliaData/DataFrames.jl/pull/2726
    aux = BitVector(undef, N)
    columns = Tables.columns(table)
    for name in names(v)
        col = Tables.getcolumn(columns, name)
        if Missing <: eltype(col)
            aux .= .!ismissing.(col)
            out .&= aux
        end
    end
    return out
end

function materialize(table, v::ClusterCovariance)
    Tables.istable(table) || throw(ArgumentError("materialize requires a table input"))
    ns = names(v)
    return cluster(ns, map(n -> GroupedArray(Tables.getcolumn(table, n)), ns))
end

"""
    nclusters(vce::ClusterCovariance)

Return the number of clusters for each dimension/way of clustering.
"""
function nclusters(v::ClusterCovariance)
    NamedTuple{names(v)}(map(x -> x.ngroups, v.clusters))
end

function df_FStat(x::RegressionModel, v::ClusterCovariance, ::Bool)
    minimum(nclusters(v)) - 1
end

function S_hat(x::RegressionModel, v::ClusterCovariance)
    # Cameron, Gelbach, & Miller (2011): section 2.3
    dim = size(modelmatrix(x), 2) * size(residuals(x), 2)
    S = zeros(dim, dim)
    for c in combinations(1:length(v))
        if length(c) == 1
            g = v.clusters[c[1]]
        else
            g = GroupedArray((v.clusters[i] for i in c)...)
        end
        S += (-1)^(length(c) - 1) * helper_cluster(modelmatrix(x), residuals(x), g)
    end
    # scale total vcov estimate by ((N-1)/(N-K)) * (G/(G-1))
    # another option would be to adjust each matrix given by helper_cluster by number of its categories
    # both methods are presented in Cameron, Gelbach and Miller (2011), section 2.3
    # I choose the first option following reghdfe
    G = minimum(nclusters(v))
    rmul!(S, (size(modelmatrix(x), 1) - 1) / dof_residual(x) * G / (G - 1))
end

# res is a Vector in OLS, Matrix in IV
function helper_cluster(X::Matrix, res::Union{Vector, Matrix}, g::GroupedArray)
    X2 = zeros(eltype(X), g.ngroups, size(X, 2) * size(res, 2))
    idx = 0
    for k in 1:size(res, 2)
        for j in 1:size(X, 2)
            idx += 1
            @inbounds @simd for i in 1:size(X, 1)
                X2[g.groups[i], idx] += X[i, j] * res[i, k]
            end
        end
    end
    return Symmetric(X2' * X2)
end

function StatsBase.vcov(x::RegressionModel, v::ClusterCovariance)
    xtx = inv(crossmodelmatrix(x))
    pinvertible(Symmetric(xtx * S_hat(x, v) * xtx))
end

ext: jl


text:
using SparsityOperators
using Test
using Random
using LinearAlgebra
using FFTW

include("testNormalOp.jl")
include("testOperators.jl")


text:
module InfiniteOpt

# Import and export JuMP
import Reexport
Reexport.@reexport using JuMP

# Import the necessary packages.
import Distributions
import DataStructures
import FastGaussQuadrature
import AbstractTrees
import LeftChildRightSiblingTrees
import LinearAlgebra
import MutableArithmetics
using Base.Meta

# Import and export SpecialFunctions for the NLP interface
Reexport.@reexport using SpecialFunctions

# Make useful aliases (note we get MOI and MOIU from JuMP)
const JuMPC = JuMP.Containers
const MOIUC = MOIU.CleverDicts
const _LCRST = LeftChildRightSiblingTrees
const _MA = MutableArithmetics
export JuMPC, MOIUC # this makes these accessible to the submodules

# Import the Collections module
include("Collections/Collections.jl")

# Import core methods
include("datatypes.jl")
include("infinite_domains.jl")
include("scalar_parameters.jl")
include("array_parameters.jl")
include("variable_basics.jl")
include("infinite_variables.jl")
include("semi_infinite_variables.jl")
include("point_variables.jl")
include("finite_variables.jl")
include("nlp.jl")
include("expressions.jl")
include("measures.jl")

# Import and export MeasureToolbox
include("MeasureToolbox/MeasureToolbox.jl")
Reexport.@reexport using .MeasureToolbox

# import more core methods
include("derivatives.jl")
include("constraints.jl")
include("macros.jl")
include("objective.jl")
include("measure_expansions.jl")
include("derivative_evaluations.jl")
include("optimize.jl")
include("results.jl")
include("show.jl")
include("utilities.jl")
include("general_variables.jl")

# Import and export TranscriptionOpt
include("TranscriptionOpt/TranscriptionOpt.jl")
Reexport.@reexport using .TranscriptionOpt

# Define additional stuff that should not be exported
const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

# Following JuMP, export everything that doesn't start with a _
for sym in names(@__MODULE__, all = true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_") || startswith(sym_string, "@_")
        continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") && Base.isidentifier(sym_string[2:end])))
        continue
    end
    @eval export $sym
end

end # module


text:
# u'=f(t,u,v,du) with 0 = g(t,u,v)
type ConstrainedODEProblem{uType,duType,tType,isinplace,F,G,C,MM} <: AbstractDAEProblem{uType,duType,tType,isinplace}
  f::F
  g::G
  u0::uType
  v0::duType
  tspan::Tuple{tType,tType}
  callback::C
  mass_matrix::MM
end

function ConstrainedODEProblem(f,g,u0,v0,tspan; iip = isinplace(f,4), callback = nothing,mass_matrix=I)
  ConstrainedODEProblem{typeof(u0),typeof(v0),eltype(tspan),iip,typeof(f),typeof(g),typeof(callback),typeof(mass_matrix)}(f,g,u0,v0,tspan,callback,mass_matrix)
end

# u' = f[1](t,u,v,du) + f[2](t,u,v,du) + ... and a single 0 = g(t,u,v)
type SplitConstrainedODEProblem{uType,duType,tType,isinplace,F,G,C,MM} <: AbstractDAEProblem{uType,duType,tType,isinplace}
  f::F
  g::G
  u0::uType
  v0::duType
  tspan::Tuple{tType,tType}
  callback::C
  mass_matrix::MM
end

function SplitConstrainedODEProblem(f,g,u0,v0,tspan; iip = isinplace(f,4), callback = nothing,mass_matrix=I)
  @assert typeof(f) <: Tuple
  @assert typeof(g) <: Tuple
  SplitConstrainedODEProblem{typeof(u0),typeof(v0),eltype(tspan),iip,typeof(f),typeof(g),typeof(callback),typeof(_mm)}(f,g,u0,v0,tspan,callback,mass_matrix)
end

# Tuple of u[i]'=f[1] and a single 0 = g(t,u,v)
type PartitionedConstrainedODEProblem{uType,duType,tType,isinplace,F,G,C,MM} <: AbstractDAEProblem{uType,duType,tType,isinplace}
  f::F
  g::G
  u0::uType
  v0::duType
  tspan::Tuple{tType,tType}
  callback::C
  mass_matrix::MM
end

function PartitionedConstrainedODEProblem(f,g,u0,v0,tspan; iip = isinplace(f,4), callback = nothing,mass_matrix=nothing)
  @assert typeof(f) <: Tuple
  @assert typeof(g) <: Tuple
  if mass_matrix == nothing
    _mm = ((I for i in 1:length(f))...)
  end
  PartitionedConstrainedODEProblem{typeof(u0),typeof(v0),eltype(tspan),iip,typeof(f),typeof(g),typeof(callback),typeof(_mm)}(f,g,u0,v0,tspan,callback,_mm)
end


text:
MP.promote_rule_constant(::Type{T}, ::Type{<:DMonomialLike{C}}) where {C, T} = Term{C, promote_type(T, Int)}
MP.promote_rule_constant(::Type{S}, ::Type{Term{C, T}}) where {S, C, T} = Term{C, promote_type(S, T)}
MP.promote_rule_constant(::Type{S}, ::Type{<:TermPoly{C, T}}) where {S, C, T} = Polynomial{C, promote_type(S, T)}
MP.promote_rule_constant(::Type, ::Type{Term{C}}) where {C} = Any
MP.promote_rule_constant(::Type, ::Type{Polynomial{C}}) where {C} = Any
Base.promote_rule(::Type{Term{C}}, ::Type{Term{C, T}}) where {C, T} = Term{C}
Base.promote_rule(::Type{Term{C, T}}, ::Type{Term{C}}) where {C, T} = Term{C}
Base.promote_rule(::Type{Term{C}}, ::Type{<:DMonomialLike{C}}) where {C} = Term{C}
Base.promote_rule(::Type{<:DMonomialLike{C}}, ::Type{Term{C}}) where {C} = Term{C}

Base.convert(::Type{Term{C}}, m::DMonomialLike{C}) where {C} = convert(Term{C,Int}, m)
Base.convert(::Type{Polynomial{C}}, t::Union{TermPoly{C},DMonomialLike{C}}) where {C} = MP.polynomial(t)


text:
immutable BrownianBridge
end

function cumsim(::BrownianBridge, t, k)
    b = cumsim(BrownianMotion(), t, k)
    t1 = t[1]
    dt = t[end] - t1
    result = zeros(b)
    for i in eachindex(t)
        result[i,:] .= b[i,:] .- b[end,:] * (t[i] - t1) / dt
    end
    result
end


text:
module RSDSGE

import SymPy
import SymPy: Sym, jacobian # Used a lot, and doesn't lead to conflict
import Base.copy
import NLsolve

include("constructor/modeltypes.jl")
include("constructor/utils.jl")
include("constructor/constructor.jl")
include("constructor/parsefile.jl")
include("constructor/perturbation.jl")
include("modelutils.jl")
include("solve/solve.jl")

# Functions to be exported
export RSDSGEModel, RSDSGESolution, solve, copy, updateparameters!, findSS!

end


text:
using Revise
using ADCME
using NNFEM
using JLD2
using PyCall
using LinearAlgebra
using Random, Distributions

reset_default_graph()

#ρ = 1153.4 kg/m^3 = 1153.4 kg/m^3 * 0.000076 m = 0.0876584 kg/m^2
#E = 944836855.0kg/m/s^2 * 0.000076m = 71807.60098kg/s^2 = 0.07180760098kg/ms^2
prop = Dict("name"=> "PlaneStress", "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)

stress_scale = 1.0
strain_scale = 1.0

mat = PlaneStress(prop)
H0 = mat.H

nntype = "ae_scaled"
idx = 1
H_function = "None"
include("nnutil.jl")

n = 50
X_1d = LinRange(-2.0, 2.0, n)
X = zeros(Float64,  n^3, 3)
for i1 = 1:n
    for i2 = 1:n
        for i3 = 1:n
            X[(i1-1)*n^2 + (i2-1)*n + i3, :] .= X_1d[i1], X_1d[i2], X_1d[i3]
        end
    end
end
ε_μ = 0.0; ε_σ = 0.02
noise = rand(Normal(ε_μ, ε_σ), (n^3,3))
Y = X * H0 + noise

loss = constant(0.0)
x = constant(X)
dummy = x
y = nn(x[:,1:3], dummy, dummy)
loss += mean((y-Y)^2)

sess = Session(); init(sess)
@show run(sess, loss)
for i = 1:10
    BFGS!(sess, loss, 1000)
    ADCME.save(sess, "Data/NNLearn_$(nntype)_ite$(i).mat")
end


text:
# ------------------------------------------
# Adjoint julia script
# ------------------------------------------

t = ccall(:clock, Int32, ())
println(t)

prim = zeros(3)
w = [1.0, 0.0, 1.0]
γ = 1.4
ccall(
    #(:hello, :libkit), # name of C function and library
    (:hello, "./c/libkit.so"), # name of C function and library
    Nothing,              # output type
    (),                        # tuple of input types
                                  # name of Julia variable to pass in
)
ccall(
    (:conserve_prim_1d, "./c/libkit.so"), # name of C function and library
    #(:conserve_prim_1d, "./c/kinetic.so"), # name of C function and library
    Nothing,              # output type
    (Ref{Cdouble}, Ref{Cdouble}, Cdouble,),                        # tuple of input types
    prim, w, γ                                  # name of Julia variable to pass in
)
println(prim)

#=
using Kinetic, BenchmarkTools

u = collect(-5.0:0.1:5.0)
nu = length(u)
weights = ones(nu) .* 0.5

fw = zeros(3)
fh = zeros(nu)
fb = zeros(nu)

inK = 2
γ = 5.0 / 3.0
primL = [1., 0., 1.]
wL = prim_conserve(primL, γ)
hL = maxwellian(u, primL) |> Array;
bL = hL .* 2 ./ (2.)
shL = zeros(nu)
sbL = zeros(nu)
lenL = 0.1

primR = [0.5, 0., 1.]
wR = prim_conserve(primR, γ)
hR = maxwellian(u, primR) |> Array;
bR = hR .* 2 ./ (2.)
shR = zeros(nu)
sbR = zeros(nu)
lenR = 0.1

muref = 0.001
omega = 0.72
prandtl = 1.0
dt = 1e-4

output_ptr = ccall(
        (:gsl_permutation_alloc, :libgsl), # name of C function and library
        Ptr{gsl_permutation},              # output type
        (Csize_t,),                        # tuple of input types
        n                                  # name of Julia variable to pass in
    )

#--- kfvs ---#
@btime ccall(
    (:__kinetic_MOD_flux_kfvs_2f1v, "./fortran/kitmod.so"),
    Nothing,
    (Ref{Float64}, Ref{Float64}, Ref{Float64},
    Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
    Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
    Ref{Int}, Ref{Float64}, Ref{Float64}, Ref{Float64}),
    fw,
    fh,
    fb,
    hL,
    bL,
    shL,
    sbL,
    lenL,
    hR,
    bR,
    shR,
    sbR,
    lenR,
    nu,
    u,
    weights,
    dt,
)

@btime flux_kfvs!(fw, fh, fb, hL, bL, hR, bR, u, weights, dt, shL, sbL, shR, sbR)

#--- ugks ---#
@btime ccall(
    (:__kinetic_MOD_flux_ugks_2f1v, "./fortran/kitmod.so"),
    Nothing,
    (Ref{Float64}, Ref{Float64}, Ref{Float64},
    Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
    Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
    Ref{Int}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}),
    fw,
    fh,
    fb,
    wL,
    hL,
    bL,
    shL,
    sbL,
    lenL,
    wR,
    hR,
    bR,
    shR,
    sbR,
    lenR,
    nu,
    u,
    weights,
    inK,
    γ,
    muref,
    omega,
    prandtl,
    dt,
)

@btime flux_ugks!(fw, fh, fb, wL, hL, bL, wR, hR, bR, u, weights, inK, γ, muref, omega, prandtl, dt, lenL, lenR, shL, sbL, shR, sbR)

"""
Results on my NUC with i7-8559U
@101 velocity points

Kinetic.jl
- kfvs ~ 6.747 μs (13 allocations: 11.38 KiB)
- ugks ~ 13.344 μs (123 allocations: 20.94 KiB) on my NUC

Kifort
- kfvs ~ 5.421 μs (37 allocations: 800 bytes)
- ugks ~ 11.413 μs (55 allocations: 1.09 KiB) on my NUC

"""
=#

max_stars_repo_name: vavrines/Kic

max_issues_repo_name: vavrines/Kic

max_forks_repo_name: vavrines/Kic


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct MonitoredLine <: ACBranch
        name::String
        available::Bool
        active_power_flow::Float64
        reactive_power_flow::Float64
        arc::Arc
        r::Float64
        x::Float64
        b::NamedTuple{(:from, :to), Tuple{Float64, Float64}}
        flow_limits::NamedTuple{(:from_to, :to_from), Tuple{Float64, Float64}}
        rate::Float64
        angle_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        services::Vector{Service}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `active_power_flow::Float64`
- `reactive_power_flow::Float64`
- `arc::Arc`
- `r::Float64`: System per-unit value, validation range: `(0, 4)`, action if invalid: `warn`
- `x::Float64`: System per-unit value, validation range: `(0, 4)`, action if invalid: `warn`
- `b::NamedTuple{(:from, :to), Tuple{Float64, Float64}}`: System per-unit value, validation range: `(0, 2)`, action if invalid: `warn`
- `flow_limits::NamedTuple{(:from_to, :to_from), Tuple{Float64, Float64}}`: TODO: throw warning above max SIL
- `rate::Float64`: TODO: compare to SIL (warn) (theoretical limit)
- `angle_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`, validation range: `(-1.571, 1.571)`, action if invalid: `error`
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct MonitoredLine <: ACBranch
    name::String
    available::Bool
    active_power_flow::Float64
    reactive_power_flow::Float64
    arc::Arc
    "System per-unit value"
    r::Float64
    "System per-unit value"
    x::Float64
    "System per-unit value"
    b::NamedTuple{(:from, :to), Tuple{Float64, Float64}}
    "TODO: throw warning above max SIL"
    flow_limits::NamedTuple{(:from_to, :to_from), Tuple{Float64, Float64}}
    "TODO: compare to SIL (warn) (theoretical limit)"
    rate::Float64
    angle_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function MonitoredLine(name, available, active_power_flow, reactive_power_flow, arc, r, x, b, flow_limits, rate, angle_limits, services=Device[], ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), )
    MonitoredLine(name, available, active_power_flow, reactive_power_flow, arc, r, x, b, flow_limits, rate, angle_limits, services, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function MonitoredLine(; name, available, active_power_flow, reactive_power_flow, arc, r, x, b, flow_limits, rate, angle_limits, services=Device[], ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), )
    MonitoredLine(name, available, active_power_flow, reactive_power_flow, arc, r, x, b, flow_limits, rate, angle_limits, services, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function MonitoredLine(::Nothing)
    MonitoredLine(;
        name="init",
        available=false,
        active_power_flow=0.0,
        reactive_power_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        r=0.0,
        x=0.0,
        b=(from=0.0, to=0.0),
        flow_limits=(from_to=0.0, to_from=0.0),
        rate=0.0,
        angle_limits=(min=-1.571, max=1.571),
        services=Device[],
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

"""Get [`MonitoredLine`](@ref) `name`."""
get_name(value::MonitoredLine) = value.name
"""Get [`MonitoredLine`](@ref) `available`."""
get_available(value::MonitoredLine) = value.available
"""Get [`MonitoredLine`](@ref) `active_power_flow`."""
get_active_power_flow(value::MonitoredLine) = get_value(value, value.active_power_flow)
"""Get [`MonitoredLine`](@ref) `reactive_power_flow`."""
get_reactive_power_flow(value::MonitoredLine) = get_value(value, value.reactive_power_flow)
"""Get [`MonitoredLine`](@ref) `arc`."""
get_arc(value::MonitoredLine) = value.arc
"""Get [`MonitoredLine`](@ref) `r`."""
get_r(value::MonitoredLine) = value.r
"""Get [`MonitoredLine`](@ref) `x`."""
get_x(value::MonitoredLine) = value.x
"""Get [`MonitoredLine`](@ref) `b`."""
get_b(value::MonitoredLine) = value.b
"""Get [`MonitoredLine`](@ref) `flow_limits`."""
get_flow_limits(value::MonitoredLine) = get_value(value, value.flow_limits)
"""Get [`MonitoredLine`](@ref) `rate`."""
get_rate(value::MonitoredLine) = get_value(value, value.rate)
"""Get [`MonitoredLine`](@ref) `angle_limits`."""
get_angle_limits(value::MonitoredLine) = value.angle_limits
"""Get [`MonitoredLine`](@ref) `services`."""
get_services(value::MonitoredLine) = value.services
"""Get [`MonitoredLine`](@ref) `ext`."""
get_ext(value::MonitoredLine) = value.ext
"""Get [`MonitoredLine`](@ref) `time_series_container`."""
get_time_series_container(value::MonitoredLine) = value.time_series_container
"""Get [`MonitoredLine`](@ref) `internal`."""
get_internal(value::MonitoredLine) = value.internal

"""Set [`MonitoredLine`](@ref) `name`."""
set_name!(value::MonitoredLine, val) = value.name = val
"""Set [`MonitoredLine`](@ref) `available`."""
set_available!(value::MonitoredLine, val) = value.available = val
"""Set [`MonitoredLine`](@ref) `active_power_flow`."""
set_active_power_flow!(value::MonitoredLine, val) = value.active_power_flow = set_value(value, val)
"""Set [`MonitoredLine`](@ref) `reactive_power_flow`."""
set_reactive_power_flow!(value::MonitoredLine, val) = value.reactive_power_flow = set_value(value, val)
"""Set [`MonitoredLine`](@ref) `arc`."""
set_arc!(value::MonitoredLine, val) = value.arc = val
"""Set [`MonitoredLine`](@ref) `r`."""
set_r!(value::MonitoredLine, val) = value.r = val
"""Set [`MonitoredLine`](@ref) `x`."""
set_x!(value::MonitoredLine, val) = value.x = val
"""Set [`MonitoredLine`](@ref) `b`."""
set_b!(value::MonitoredLine, val) = value.b = val
"""Set [`MonitoredLine`](@ref) `flow_limits`."""
set_flow_limits!(value::MonitoredLine, val) = value.flow_limits = set_value(value, val)
"""Set [`MonitoredLine`](@ref) `rate`."""
set_rate!(value::MonitoredLine, val) = value.rate = set_value(value, val)
"""Set [`MonitoredLine`](@ref) `angle_limits`."""
set_angle_limits!(value::MonitoredLine, val) = value.angle_limits = val
"""Set [`MonitoredLine`](@ref) `services`."""
set_services!(value::MonitoredLine, val) = value.services = val
"""Set [`MonitoredLine`](@ref) `ext`."""
set_ext!(value::MonitoredLine, val) = value.ext = val
"""Set [`MonitoredLine`](@ref) `time_series_container`."""
set_time_series_container!(value::MonitoredLine, val) = value.time_series_container = val


text:
#  ======================================================
#  Structures
#  ======================================================

"""
    EdgeValGraph{V <: Integer, E_VALS <: AbstractTuple, E_VALS_C}

A type representing an undirected simple graph with edge values.

The element type `V` specifies the type of the vertex indices and `E_VAL`
specifies the type of the edge values. User should usually not specify `E_VAL_C`
by themself but rather let a constructor do that.
"""
mutable struct EdgeValGraph{  V <: Integer,
                            E_VALS <: AbstractTuple,
                            E_VALS_C
                         } <: AbstractEdgeValGraph{V, E_VALS}

    ne::Int
    fadjlist::Adjlist{V}
    edgevals::E_VALS_C
end

"""
    EdgeValOutDiGraph{V <: Integer, E_VAL <: AbstractTuple, E_VAL_C}

A type representing a directed simple graph with edge values.

The element type `V` specifies the type of the vertex indices and `E_VAL`
specifies the type of the edge values. User should usually not specify `E_VAL_C`
by themself but rather let a constructor do that.
"""
mutable struct EdgeValOutDiGraph{ V <: Integer,
                                E_VALS <: AbstractTuple,
                                E_VALS_C
                              } <: AbstractEdgeValGraph{V, E_VALS}

    ne::Int
    fadjlist::Adjlist{V}
    edgevals::E_VALS_C
end

"""
    EdgeValDiGraph{V <: Integer, E_VAL <: AbstractTuple, E_VAL_C}

A type representing a directed simple graph with edge values.

The element type `V` specifies the type of the vertex indices and `E_VAL`
specifies the type of the edge values. User should usually not specify `E_VAL_C`
but rather let a constructor do that.
"""
mutable struct EdgeValDiGraph{    V <: Integer,
                                E_VALS <: AbstractTuple,
                                E_VALS_C
                           } <: AbstractEdgeValGraph{V, E_VALS}

    ne::Int
    fadjlist::Adjlist{V}
    badjlist::Adjlist{V}
    edgevals::E_VALS_C
    redgevals::E_VALS_C
end

#  ======================================================
#  Constructors
#  ======================================================

# TODO this helpers maybe belong somewhere else and need a better name
construct_E_VAL(edgeval_types::Tuple) = Tuple{ (T for T in edgeval_types)... }
construct_E_VAL(edgeval_types::NamedTuple) =
    NamedTuple{ Tuple(typeof(edgeval_types).names), Tuple{ (T for T in edgeval_types)... }}

const default_edgeval_types = (weight=Float64,)

function create_edgevals(n, E_VAL::Type{<:Tuple})
    return Tuple( Adjlist{T}(n) for T in E_VAL.types )
end

function create_edgevals(n, E_VAL::Type{<:NamedTuple})
    return NamedTuple{Tuple(E_VAL.names)}(Tuple( Adjlist{T}(n) for T in E_VAL.types ))
end

#  ------------------------------------------------------
#  Constructors for empty graphs
#  ------------------------------------------------------

"""
    EdgeValGraph{V}(n, edgeval_types=$(default_edgeval_types))
    EdgeValGraph{V, E_VAL}

Construct a `EdgeValGraph` with `n` vertices and 0 edges with of types
`edgeval_types`.

If omitted, the element type `V` is the type of `n`.
"""
function EdgeValGraph{V, E_VAL}(n::Integer) where {V <: Integer, E_VAL <: AbstractTuple}

    fadjlist = Adjlist{V}(n)
    edgevals = create_edgevals(n, E_VAL)
    E_VAL_C = typeof(edgevals)

    return EdgeValGraph{V, E_VAL, E_VAL_C}(0, fadjlist, edgevals)
end

"""
    EdgeValOutDiGraph{V}(n, edgeval_types=$(default_edgeval_types))

Construct a `EdgeValOutDiGraph` with `n` vertices and 0 edges of types
`edgeval_types`.
If omitted, the element type `V` is the type of `n`.

"""
function EdgeValOutDiGraph{V, E_VAL}(n::Integer) where {V<:Integer, E_VAL}

    fadjlist = Adjlist{V}(n)
    edgevals = create_edgevals(n, E_VAL)
    E_VAL_C = typeof(edgevals)

    return EdgeValOutDiGraph{V, E_VAL, E_VAL_C}(0, fadjlist, edgevals)
end

"""
    EdgeValDiGraph{V}(n, edgeval_types=$(default_edgeval_types))

Construct a `EdgeValDiGraph` with `n` vertices and 0 edges with value-types
`edgeval_types`.

If omitted, the element type `V` is the type of `n`.
"""
function EdgeValDiGraph{V, E_VAL}(n::Integer) where {V<:Integer, E_VAL}
    fadjlist = Adjlist{V}(n)
    badjlist = Adjlist{V}(n)
    edgevals = create_edgevals(n, E_VAL)
    redgevals = create_edgevals(n, E_VAL)
    E_VAL_C = typeof(edgevals)

    return EdgeValDiGraph{V, E_VAL, E_VAL_C}(
                0, fadjlist, badjlist, edgevals, redgevals)
end

for G in (:EdgeValGraph, :EdgeValOutDiGraph, :EdgeValDiGraph)
    @eval function $G(n::V, edgeval_types::AbstractTupleOfTypes=default_edgeval_types) where {V <: Integer}
        E_VAL = construct_E_VAL(edgeval_types)
        return $G{V, E_VAL}(n)
    end

    @eval function $G{V}(n, edgeval_types::AbstractTupleOfTypes=default_edgeval_types) where {V <: Integer}
        E_VAL = construct_E_VAL(edgeval_types)
        return $G{V, E_VAL}(V(n))
    end
end

#  ------------------------------------------------------
#  Constructor from other graphs
#  ------------------------------------------------------

"""
    EdgeValGraph(undef, g::SimpleGraph, edgeval_types=$(default_edgeval_types))
    EdgeValGraph([edgeval_initializer,] g::SimpleGraph, edgeval_types=$(default_edgeval_types))

Construct a `EdgeValGraph` with the same structure as `g` with uninitialized edge
values of types `edgeval_types`.
Construct a `EdgeValGraph` with the same structure as `g`.
The optional argument `edgeval_initializer` takes a
function that assigns to each edge (s, d) an edge value. If it is not given,
then each edge gets the value `default_edgeval(edgeval_types)`.
"""
function EdgeValGraph end

function EdgeValGraph(
                ::UndefInitializer,
                g::SimpleGraph{V},
                edgeval_types::AbstractTupleOfTypes=default_edgeval_types) where {V}

    n = nv(g)
    fadjlist = deepcopy_adjlist(g.fadjlist)

    E_VAL = construct_E_VAL(edgeval_types)
    E_VAL_C = edgevals_container_type(Val(E_VAL)) # TODO ?
    edgevals = E_VAL_C( Adjlist{T}(undef, n) for T in E_VAL.types )
    for s in OneTo(n)
        for (i, T) in enumerate(E_VAL.types)
            edgevals[i][s] = Vector{T}(undef, length(fadjlist[s]))
        end
    end
    return EdgeValGraph{V, E_VAL, E_VAL_C}(ne(g), fadjlist, edgevals)
end

"""
    EdgeValOutDiGraph(undef, g::SimpleDiGraph, edgeval_types=$(default_edgeval_types))
    EdgeValOutDiGraph([edgeval_initializer,] g::SimpleDiGraph, edgeval_types=$(default_edgeval_types))

Construct a `EdgeValOutDiGraph` with the same structure as `g` with uninitialized edge values of types `edgeval_types`.
Construct a `EdgeValOutDiGraph` with the same structure as `g`.
The optional argument `edgeval_initializer` takes a
function that assigns to each edge (s, d) an edge value. If it is not given,
then each edge gets the value `default_edgeval(edgeval_types)`.
"""
function EdgeValOutDiGraph end

function EdgeValOutDiGraph(
                     ::UndefInitializer,
                    g::SimpleDiGraph{V},
        edgeval_types::AbstractTupleOfTypes=default_edgeval_types) where {V}

    n = nv(g)
    fadjlist = deepcopy_adjlist(g.fadjlist)

    E_VAL = construct_E_VAL(edgeval_types)
    E_VAL_C = edgevals_container_type(Val(E_VAL)) # TODO ?
    edgevals = E_VAL_C( Adjlist{T}(undef, n) for T in E_VAL.types )
    for s in OneTo(n)
        for (i, T) in enumerate(E_VAL.types)
            edgevals[i][s] = Vector{T}(undef, length(fadjlist[s]))
        end
    end
    return EdgeValOutDiGraph{V, E_VAL, E_VAL_C}(ne(g), fadjlist, edgevals)
end

"""
    EdgeValDiGraph(undef, g::SimpleDiGraph, edgeval_types=$(default_edgeval_types))
    EdgeValDiGraph([edgeval_initializer,] g::SimpleDiGraph, edgeval_types=$(default_edgeval_types))

Construct a `EdgeValDiGraph` with the same structure as `g` with uninitialized edge values of types `edgeval_types`.
Construct a `EdgeValDiGraph` with the same structure as `g`.
The optional argument `edgeval_initializer` takes a
function that assigns to each edge (s, d) an edge value. If it is not given,
then each edge gets the value `default_edgeval(edgeval_types)`.
"""
function EdgeValDiGraph end

function EdgeValDiGraph(
                     ::UndefInitializer,
        g            ::SimpleDiGraph{V},
        edgeval_types::AbstractTupleOfTypes=default_edgeval_types
        ) where {V}

    n = nv(g)
    fadjlist = deepcopy_adjlist(g.fadjlist)
    badjlist = deepcopy_adjlist(g.badjlist)

    E_VAL = construct_E_VAL(edgeval_types)
    E_VAL_C = edgevals_container_type(Val(E_VAL)) # TODO ?
    edgevals = E_VAL_C( Adjlist{T}(undef, n) for T in E_VAL.types )
    redgevals = E_VAL_C( Adjlist{T}(undef, n) for T in E_VAL.types )
    for s in OneTo(n)
        for (i, T) in enumerate(E_VAL.types)
            edgevals[i][s] = Vector{T}(undef, length(fadjlist[s]))
            redgevals[i][s] = Vector{T}(undef, length(badjlist[s]))
        end
    end
    return EdgeValDiGraph{V, E_VAL, E_VAL_C}(ne(g), fadjlist, badjlist, edgevals, redgevals)
end

function EdgeValGraph(
        edgeval_initializer::Base.Callable,
        g                  ::SimpleGraph,
        edgeval_types      ::AbstractTupleOfTypes   = default_edgeval_types
        )

    gv = EdgeValGraph(undef, g, edgeval_types)

    n = nv(g)
    # TODO there is a more efficient method for this
    for e in edges(g)
        s, d = Tuple(e)
        set_edgevals!(gv, s, d, edgeval_initializer(s, d))
    end

    return gv
end

for G in (:EdgeValOutDiGraph, :EdgeValDiGraph)
    @eval function $G(edgeval_initializer::Base.Callable,
                      g::SimpleDiGraph,
                      edgeval_types::AbstractTupleOfTypes=default_edgeval_types)
        gv = $G(undef, g, edgeval_types)

        n = nv(g)
        # TODO there is a more efficient method for this
        for e in edges(g)
            s, d = Tuple(e)
            set_edgevals!(gv, s, d, edgeval_initializer(s, d))
        end

        return gv
    end
end

# TODO these constructors do not work
#=
EdgeValGraph(g::SimpleGraph, edgeval_types::AbstractTupleOfTypes=default_edgeval_types) =
    EdgeValGraph(default_edgevals, g, edgeval_types)

EdgeValOutDiGraph(g::SimpleDiGraph, edgeval_types::AbstractTupleOfTypes=default_edgeval_types) =
    EdgeValOutDiGraph(default_edgevals, g, edgeval_types)

EdgeValDiGraph(g::SimpleDiGraph, edgeval_types::AbstractTupleOfTypes=default_edgeval_types) =
    EdgeValDiGraph(default_edgevals, g::SimpleGraph, edgeval_types)
=#

# =========================================================
# Interface
# =========================================================

#  ------------------------------------------------------
#  add_edge!
#  ------------------------------------------------------

"""
    add_edge!(g::AbstractEdgeValGraph{V, E_VALS}, s, d, value=default_edgeval(E_VALS))

Add an edge `e = (s, d, [value])` to a graph `g` and set the edge value.

Return `true` if the edge was added successfully, otherwise return `false`.
If the edge already exists, return `false` but still change the edge value.
"""
function add_edge! end

function add_edge!(g::EdgeValGraph{V, E_VALS},
                   s::Integer,
                   d::Integer,
                   values::E_VALS=default_edgeval(E_VALS)) where {V, E_VALS}

    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds

    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        # edge already there, replace value, but return false
        set_values_for_index!(edgevals, s, index, values)
        s == d && return false # selfloop
        index = searchsortedfirst(g.fadjlist[d], s)
        set_values_for_index!(edgevals, d, index, values)
        return false
    end

    insert!(list, index, d)
    insert_values_for_index!(edgevals, s, index, values)
    g.ne += 1

    s == d && return true # selfloop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    insert_values_for_index!(edgevals, d, index, values)
    return true # edge successfully added
end

function add_edge!(g::EdgeValOutDiGraph{V, E_VALS},
                   s::Integer,
                   d::Integer,
                   value::E_VALS=default_edgeval(E_VALS)) where {V, E_VALS}

    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        # edge already there, replace value, but return false
        set_values_for_index!(edgevals, s, index, value)
        return false
    end

    insert!(list, index, d)
    insert_values_for_index!(edgevals, s, index, value)
    g.ne += 1
    return true # edge successfully added
end

function add_edge!(g::EdgeValDiGraph{V, E_VALS},
                   s::Integer,
                   d::Integer,
                   value::E_VALS=default_edgeval(E_VALS)) where {V, E_VALS}
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        # edge already there, replace value, but return false
        set_values_for_index!(g.edgevals, s, index, value)
        return false
    end

    insert!(list, index, d)
    insert_values_for_index!(g.edgevals, s, index, value)
    g.ne += 1

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    set_values_for_index!(g.redgevals, d, index, value)

    return true # edge successfully added
end

#  ------------------------------------------------------
#  rem_edge!
#  ------------------------------------------------------

function rem_edge!(g::EdgeValGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false
    deleteat!(list, index)
    delete_values_for_index!(edgevals, s, index)

    g.ne -= 1
    s == d && return true # self-loop

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)
    delete_values_for_index!(edgevals, d, index)

    return true
end

function rem_edge!(g::EdgeValOutDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    edgevals = g.edgevals
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false
    deleteat!(list, index)
    delete_values_for_index!(edgevals, s, index)

    g.ne -= 1
    return true
end

function rem_edge!(g::EdgeValDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false
    deleteat!(list, index)
    delete_values_for_index!(g.edgevals, s, index)

    g.ne -= 1

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)
    delete_values_for_index!(g.redgevals, d, index)

    return true
end

#  ------------------------------------------------------
#  has_edge
#  ------------------------------------------------------

"""
    has_edge!(g::AbstractEdgeValGraph, s, d)

Return `true` if `g` has an edge from node `s` to `d`.
"""
function has_edge end

function has_edge(g::EdgeValGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        d = s
        list_s = list_d
    end
    return LightGraphs.insorted(d, list_s)
end

function has_edge(g::EdgeValOutDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list_s = g.fadjlist[s]

    return LightGraphs.insorted(d, list_s)
end

function has_edge(g::EdgeValDiGraph, s::Integer, d::Integer)
    verts = vertices(g)
    (s in verts && d in verts) || return false # edge out of bounds
    @inbounds list_fadj = g.fadjlist[s]
    @inbounds list_badj = g.badjlist[d]
    if length(list_fadj) > length(list_badj)
        d = s
        list_fadj = list_badj
    end
    return LightGraphs.insorted(d, list_fadj)
end

#  ------------------------------------------------------
#  get_edgevals
#  ------------------------------------------------------

# TODO more specific error

function get_edgeval(g::EdgeValGraph, s::Integer, d::Integer; key=nokey)

    if g isa OneEdgeValGraph && key === nokey
        key = 1
    end

    validkey_or_throw(g, key)

    verts = vertices(g)

    (s in verts && d in verts) || error("No such edge")
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end

    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return g.edgevals[key][s][index]
    end
    error("No such edge")
end

function get_edgeval(g::EdgeValOutDiGraph, s::Integer, d::Integer; key=nokey)

    if g isa OneEdgeValGraph && key === nokey
        key = 1
    end

    validkey_or_throw(g, key)

    verts = vertices(g)
    (s in verts && d in verts) || error("No such edge")
    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return g.edgevals[key][s][index]
    end
    error("No such edge")
end

function get_edgeval(g::EdgeValDiGraph, s::Integer, d::Integer; key=nokey)

    if g isa OneEdgeValGraph && key === nokey
        key = 1
    end

    validkey_or_throw(g, key)

    verts = vertices(g)
    (s in verts && d in verts) || error("No such edge")
    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return g.edgevals[key][s][index]
    end
    error("No such edge")
end

function get_edgeval_or(g::EdgeValGraph, s::Integer, d::Integer, alternative; key=nokey)

    if g isa OneEdgeValGraph && key === nokey
        key = 1
    end

    validkey_or_throw(g, key)

    verts = vertices(g)

    (s in verts && d in verts) || return alternative
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end

    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return g.edgevals[key][s][index]
    end
    return alternative
end

function get_edgeval_or(g::EdgeValOutDiGraph, s::Integer, d::Integer, alternative; key=nokey)

    if g isa OneEdgeValGraph && key === nokey
        key = 1
    end

    validkey_or_throw(g, key)

    verts = vertices(g)
    (s in verts && d in verts) || return alternative
    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return g.edgevals[key][s][index]
    end
    return alternative
end

function get_edgeval_or(g::EdgeValDiGraph, s::Integer, d::Integer, alternative; key=nokey)

    if g isa OneEdgeValGraph && key === nokey
        key = 1
    end

    validkey_or_throw(g, key)

    verts = vertices(g)
    (s in verts && d in verts) || return alternative
    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return g.edgevals[key][s][index]
    end
    return alternative
end

"""
    get_edgevals(g::AbstractEdgeValGraph, s, d, default=nothing)
    get_edgevals(g::AbstractEdgeValGraph, e, default=nothing)
Return the edge value for the edge `e: s => d` and `default` if that edge does not exist.
"""
function get_edgevals(g::EdgeValGraph{V, E_VAL}, s::Integer, d::Integer) where {V, E_VAL}

    verts = vertices(g)
    (s in verts && d in verts) || error("Values not found")
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return values_for_index(g.edgevals, E_VAL, s, index)
    end
    (s in verts && d in verts) || error("Values not found")
end

function get_edgevals(g::EdgeValOutDiGraph{V, E_VAL}, s::Integer, d::Integer) where {V, E_VAL}

    verts = vertices(g)
    (s in verts && d in verts) ||  error("Values not found")

    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return values_for_index(g.edgevals, E_VAL, s, index)
    end
    error("Values not found")
end

# TODO could probably be made faster by checking the shorter list
function get_edgevals(g::EdgeValDiGraph{V, E_VAL}, s::Integer, d::Integer) where
{V, E_VAL}

    verts = vertices(g)
    (s in verts && d in verts) || error("Values not found")

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        return values_for_index(g.edgevals, E_VAL, s, index)
    end
    error("Values not found")
end

function get_edgevals_or(g::EdgeValGraph{V, E_VAL}, s::Integer, d::Integer, default) where {V, E_VAL}

    verts = vertices(g)
    (s in verts && d in verts) || return default
    @inbounds list_s = g.fadjlist[s]
    @inbounds list_d = g.fadjlist[d]
    if length(list_s) > length(list_d)
        s, d = d, s
        list_s = list_d
    end
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return values_for_index(g.edgevals, E_VAL, s, index)
    end
    return default
end

function get_edgevals_or(g::EdgeValOutDiGraph{V, E_VAL}, s::Integer, d::Integer, default) where {V, E_VAL}

    verts = vertices(g)
    (s in verts && d in verts) ||  return default

    @inbounds list_s = g.fadjlist[s]
    index = searchsortedfirst(list_s, d)
    @inbounds if index <= length(list_s) && list_s[index] == d
        return values_for_index(g.edgevals, E_VAL, s, index)
    end
    return default
end

# TODO could probably be made faster by checking the shorter list
function get_edgevals_or(g::EdgeValDiGraph{V, E_VAL}, s::Integer, d::Integer, default) where
{V, E_VAL}

    verts = vertices(g)
    (s in verts && d in verts) || return default

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds if index <= length(list) && list[index] == d
        return values_for_index(g.edgevals, E_VAL, s, index)
    end
    return default
end

#  ------------------------------------------------------
#  set_edgevals!
#  ------------------------------------------------------

function set_edgeval!(g::EdgeValGraph, s::Integer, d::Integer, value; key=nokey)

    if g isa OneEdgeValGraph && key === nokey
        key = 1
    end

    validkey_or_throw(g, key)

    verts = vertices(g)
    edgevals = g.edgevals[key]

    (s in verts && d in verts) || return false

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    @inbounds edgevals[s][index] = value

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    @inbounds edgevals[d][index] = value

    return true
end

function set_edgeval!(g::EdgeValOutDiGraph, s::Integer, d::Integer, value; key=nokey)

    if g isa OneEdgeValGraph && key === nokey
        key = 1
    end

    validkey_or_throw(g, key)

    verts = vertices(g)
    edgevals = g.edgevals[key]

    (s in verts && d in verts) || return false

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    @inbounds edgevals[s][index] = value

    return true

end

function set_edgeval!(g::EdgeValDiGraph, s::Integer, d::Integer, value; key=nokey)

    if g isa OneEdgeValGraph && key === nokey
        key = 1
    end

    validkey_or_throw(g, key)

    verts = vertices(g)
    edgevals = g.edgevals[key]

    (s in verts && d in verts) || return false

    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    @inbounds edgevals[s][index] = value

    redgevals = g.redgevals[key]

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    @inbounds redgevals[d][index] = value

    return true
end

"""
    set_edgevals!(g::AbstractEdgeValGraph, s, d, value)
    set_edgevals!(g::AbstractEdgeValGraph, e, value)
Set the values of the edge `e: s -> d` to `value`. Return `true` if such an edge exists and
`false` otherwise.
"""
function set_edgevals!(g::EdgeValGraph, s::Integer, d::Integer, values)
    verts = vertices(g)
    (s in verts && d in verts) || return false
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    edgevals = g.edgevals
    set_values_for_index!(edgevals, s, index, values)

    @inbounds list = g.fadjlist[d]
    index = searchsortedfirst(list, s)
    set_values_for_index!(edgevals, d, index, values)

    return true
end

function set_edgevals!(g::EdgeValOutDiGraph, s::Integer, d::Integer, values)
    verts = vertices(g)
    (s in verts && d in verts) || return false
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    edgevals = g.edgevals
    set_values_for_index!(edgevals, s, index, values)

    return true
end

function set_edgevals!(g::EdgeValDiGraph, s::Integer, d::Integer, values)
    verts = vertices(g)
    (s in verts && d in verts) || return false
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds index <= length(list) && list[index] == d || return false
    edgevals = g.edgevals
    set_values_for_index!(g.edgevals, s, index, values)

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    set_values_for_index!(g.redgevals, d, index, values)

    return true
end

#=
set_edgevals!(g::AbstractEdgeValGraph, e::SimpleEdge, u) = set_edgevals!(g, src(e), dst(e), u)
=#

#  ------------------------------------------------------
#  is_directed
#  ------------------------------------------------------

is_directed(::Type{<:EdgeValGraph}) = false
is_directed(::Type{<:EdgeValOutDiGraph}) = true
is_directed(::Type{<:EdgeValDiGraph}) = true

is_directed(g::EdgeValGraph) = false
is_directed(g::EdgeValOutDiGraph) = true
is_directed(g::EdgeValDiGraph) = true

#  ------------------------------------------------------
#  outneighbors
#  ------------------------------------------------------

outneighbors(g::AbstractEdgeValGraph, v::Integer) = g.fadjlist[v]

#  ------------------------------------------------------
#  inneighbors
#  ------------------------------------------------------

inneighbors(g::EdgeValGraph, v::Integer) = outneighbors(g, v)
inneighbors(g::EdgeValDiGraph, v::Integer) = g.badjlist[v]

#  ------------------------------------------------------
#  outedgevals
#  ------------------------------------------------------

"""
    outedgevals(g::AbstractEdgeValGraph, v)
Return an iterator the edge values of outgoing edges from `v` to its neighbors.
The order of the neighbors is the same as for `outneighbors(g, v)`.
"""
outedgevals(g::AbstractEdgeValGraph{V, E_VAL},
            v::Integer
            ) where {V, E_VAL, E_VAL_C} =
    EdgevalsIterator(Int(v), outdegree(g, v), g.edgevals)

#  ------------------------------------------------------
#  intedgevals
#  ------------------------------------------------------

"""
    inedgevals(g::EdgeValGraph, v)
Return an iterator the edge values of incoming edges from `v` to its neighbors.
The order of the neighbors is the same as for `inneighbors(g, v)`.
"""
inedgevals(g::EdgeValGraph, v::Integer) = outedgevals(g, v)

inedgevals(g::EdgeValDiGraph{V, E_VAL, E_VAL_C}, v::Integer) where {V, E_VAL, E_VAL_C} =
    EdgevalsIterator{E_VAL, E_VAL_C}(Int(v), indegree(g, v), g.redgevals)

# ====================================================================
# Iterators
# ====================================================================

@inline function iterate(eit::ValEdgeIter{<:EdgeValGraph{V, E_VAL}}, state=(one(V), 1) ) where {V, E_VAL}
    g = eit.g
    fadjlist = g.fadjlist
    edgevals = g.edgevals
    n = V(nv(g))
    u, i = state

    @inbounds while u < n
        list_u = fadjlist[u]
        if i > length(list_u)
            u += one(u)
            i = searchsortedfirst(fadjlist[u], u)
            continue
        end
        e = ValEdge(u, list_u[i], values_for_index(edgevals, E_VAL, u, i))
        state = (u, i + 1)
        return e, state
    end

    # i > length(fadjlist[end]) || fadlist[end][i] == n

    @inbounds (n == 0 || i > length(fadjlist[n])) && return nothing

    e = ValEdge(n, n, values_for_index(edgevals, E_VAL, u, i))
    state = (u, i + 1)
    return e, state
end

function iterate(
            iter::ValEdgeIter{<:Union{EdgeValOutDiGraph{V, E_VAL}, EdgeValDiGraph{V, E_VAL}}},
            state=(one(V), 1)
    ) where {V, E_VAL}

    g = iter.g
    fadjlist = g.fadjlist
    edgevals = g.edgevals
    n = V(nv(g))
    u, i = state

    @inbounds while u <= n
        if i > length(fadjlist[u])
            u == n && return nothing

            u += one(u)
            i = 1
            continue
        end
        e = ValDiEdge(u, fadjlist[u][i], values_for_index(edgevals, E_VAL, u, i))
        return e, (u, i + 1)
    end

    return nothing
end


text:
import brainflow

# specify logging library to use
brainflow.enable_dev_brainflow_logger(Integer(brainflow.BOARD_CONTROLLER))

params = brainflow.BrainFlowInputParams()
board_shim = brainflow.BoardShim(Integer(brainflow.SYNTHETIC_BOARD), params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(256, board_shim)
brainflow.release_session(board_shim)

max_stars_repo_name: Cognixion-inc/brainflow

max_issues_repo_name: Cognixion-inc/brainflow

max_forks_repo_name: Cognixion-inc/brainflow


text:
_merger(f) = f
_merger(f::OnlineStat) = merge
_merger(f::Tup) = map(_merger, f)

function reduce(f, t::DDataset; dims = nothing, cache = false, kw...)
    if dims !== nothing
        if !isa(t, DNDSparse)
            error("`dims` only supported for NDSparse")
        end
        if !isempty(kw)
            error("unsupported keyword arguments $kw")
        end
        return _reducedim(f, t, dims, cache)
    end
    if haskey(kw, :init)
        return reduce_inited(f, t, kw.init, get(kw, :select, nothing))
    end
    select = get(kw, :select, valuenames(t))
    xs = delayedmap(t.chunks) do x
        f = isa(f, OnlineStat) ? copy(f) : f # required for > 1 chunks on the same proc
        reduce(f, x; select=select)
    end
    if f isa Tup
        g, _ = IndexedTables.init_funcs(f, false)
    else
        g = f
    end
    h = (a,b)->IndexedTables._apply(_merger(g), a,b)
    collect(get_context(), treereduce(delayed(h), xs))
end

function reduce_inited(f, t::DDataset, v0, select)
    xs = delayedmap(t.chunks) do x
        f = isa(f, OnlineStat) ? copy(f) : f
        reduce(f, x; select=select === nothing ? rows(x) : select)
    end
    g = isa(f, OnlineStat) ? merge : f
    merge(collect(get_context(), treereduce(delayed(g), xs)), v0)
end

function groupreduce(f, t::DDataset, by=pkeynames(t); kwargs...)
    @noinline function groupchunk(x)
        groupreduce(f, x, by; kwargs...)
    end
    if f isa Tup || t isa DIndexedTable
        g, _ = IndexedTables.init_funcs(f, false)
    else
        g = f
    end
    h = _merger(g)
    if (f isa ApplyColwise) && f.functions isa Union{Function, Type}
        mergef = (x,y) -> map(f.functions, x,y)
    else
        mergef = (x,y) -> IndexedTables._apply(h, x,y)
    end
    @noinline function mergechunk(x, y)
        # use NDSparse's merge
        if x isa IndexedTable
            z = merge(_convert(NDSparse, x), _convert(NDSparse, y), agg=mergef)
            _convert(IndexedTable, z)
        else
            merge(x,y, agg=mergef)
        end
    end

    t1 = fromchunks(delayedmap(groupchunk, t.chunks))
    with_overlaps(t1, true) do cs
        treereduce(delayed(mergechunk), cs)
    end
end

function groupby(f, t::DDataset, by=pkeynames(t);
                 select=(t isa DNDSparse ? valuenames(t) : Not(by)),
                 kwargs...)
    by = lowerselection(t, by)
    select = lowerselection(t, select)
    if (by isa Tup) && isempty(by)
        @noinline function _groupby(x)
            groupby(f, x, by;  select=select, kwargs...)
        end
        collect(get_context(), delayed(_groupby)(
            treereduce(delayed(_merge), t.chunks)
            )
       )
   elseif by != lowerselection(t, Keys()) || has_overlaps(t.domains, closed=true)
       subsel = Tuple(setdiff(select, by))
       # translate the new selection
       t = rechunk(t, by, subsel)
       newselect = map(select) do x
           if x in by
               findall(in(x), by)[1]
           else
               x1 = findall(in(x), subsel)
               x1[1] + length(by)
           end
       end
       return groupby(f, t; select=newselect, kwargs...)
    else
        @noinline function groupchunk(x)
            groupby(f, x, by; select=select, kwargs...)
        end

        return fromchunks(delayedmap(groupchunk, t.chunks))
    end
end

function _reducedim(f, x::DNDSparse, dims, cache)
    keep = setdiff([1:ndims(x);], dims) # TODO: Allow symbols
    if isempty(keep)
        throw(ArgumentError("to remove all dimensions, use `reduce(f, A)`"))
    end
    groupreduce(f, x, (keep...,), select=valuenames(x), cache=cache)
end

_reducedim(f, x::DNDSparse, dims::Symbol, cache) = _reducedim(f, x, [dims], cache)

"""
    reducedim_vec(f::Function, t::DNDSparse, dims)

Like `reducedim`, except uses a function mapping a vector of values to a scalar instead of a 2-argument scalar function.

See also [`reducedim`](@ref).
"""
function reducedim_vec(f, x::DNDSparse, dims)
    keep = setdiff([1:ndims(x);], dims)
    if isempty(keep)
        throw(ArgumentError("to remove all dimensions, use `reduce(f, A)`"))
    end

    t = selectkeys(x, (keep...,); agg=nothing)
    cache_thunks(groupby(f, t))
end

reducedim_vec(f, x::DNDSparse, dims::Symbol) = reducedim_vec(f, x, [dims])
Base.@deprecate aggregate_stats(s, t; by=pkeynames(t), with=valuenames(t)) groupreduce(s, t, by; select=with)


text:
module YTDL

using HTTP, LazyJSON, PrettyTables, DataFrames
export ytdl, available_codecs

abstract type AV end

struct Audio <: AV
    url::String
    length::Int64
    ext::String
    temp_ext::String
    function Audio(nt, ext)
        # Int32 reach its ceiling at 2^31-1 (≒ 256 MiB)
        # Int64 accepts up to 2^63-1 (≒ 1 EiB)
        len = haskey(nt, :length) ? parse(Int64, nt.length) : 0
        new(nt.url, len, ext.a, ext.atemp)
    end
end

struct Video <: AV
    url::String
    length::Int64
    ext::String
    temp_ext::String
    function Video(nt, ext)
        len = haskey(nt, :length) ? parse(Int64, nt.length) : 0
        new(nt.url, len, ext.v, ext.vtemp)
    end
    Video(t::Nothing, ext) = nothing
end

include("codec.jl")
include("deps.jl")
include("meta.jl")
include("downloader.jl")
include("presets.jl")
include("utils.jl")
include("youtube-dl.jl")

const USER_JL = joinpath(homedir(), ".julia", "config", "ytdl.jl")
if isfile(USER_JL)
    include(USER_JL)
end

# global setting
function ytdl(url::AbstractString; kwargs...)
    d = get(ENV, "YTDL_DIR", homedir())
    opts = haskey(kwargs, :dir) ? kwargs.data : merge((;dir=d), kwargs)
    isdir(opts.dir) || return error("$(opts.dir) does not exist.")
    get(opts, :useragent, nothing) |> HTTP.setuseragent!
    v_id = queries(url)["v"]
    ytdl_body(v_id, opts)
    return nothing
end

# download setting each video
function ytdl_body(v_id, opts)
    info = get_metadata(v_id)
    if isnothing(info)
        if get(opts, :enable_backend, false)
            return youtube_dl(v_id, opts)
        else
            println("""
                This video type (cipherd video source URL) is not yet available,
                but you can download it using youtube-dl backend-mode.
                Add option `enable_backend=true` and call `ytdl` again!
            """)
            return nothing
        end
    end
    path = joinpath(opts.dir, safename(info.title, opts))
    @async get_thumbnail(v_id, path)
    c = sellect_format(info.codecs, opts)
    ext = sellect_ext(c.a, c.v)
    downloader(
        Audio(info.codecs[c.a], ext),
        Video(get(info.codecs, c.v, nothing), ext),
        path,
        info2cmd(info, v_id, opts),
        opts
    )
end

end # module


text:
abstract type AbstractServiceFormulation end

mutable struct ServiceModel{D<:PSY.Service,
                            B<:AbstractServiceFormulation}
    service::Type{D}
    formulation::Type{B}
end

function construct_service!(canonical::Canonical,
                           service_model::ServiceModel,
                           system_formulation::Type{S},
                           sys::PSY.System;
                           kwargs...) where {S<:PM.AbstractPowerModel}

    construct_service!(canonical,
                      service_model.service,
                      service_model.formulation,
                      system_formulation,
                      sys;
                      kwargs...)

    return

end


text:
using StanSample, DiffEqBayesStan, MCMCChains
using OrdinaryDiffEq
using ParameterizedFunctions
using RecursiveArrayTools
using Distributions
using StatsPlots

function SIR!(du,u,pars,t)
    β, γ = pars
    S, I, R = u

    dS = -β*S*I
    dI = β*S*I - γ*I
    dR = γ*I
    du .= [dS, dI, dR]
end

t_obs = collect(1.0:1.0:218)
prob = ODEProblem(SIR!, [0.99, 0.01, 0.0], (0.0, t_obs[end]), [0.20, 0.15])
sol = solve(prob, Tsit5(), reltol = 1e-6, abstol = 1e-6, saveat = t_obs)

begin
    y = DataFrame(sol',[:S,:I,:R])
    plot(t_obs, y."S", lab="S")
    plot!(t_obs, y."I", lab="I")
    plot!(t_obs, y."R", lab="R")
    savefig(joinpath(@__DIR__, "sir.png"))
end

init_params = [1.0, 1.0]
data = convert(Array, sol)
priors = [Normal(), Normal()]

bayesian_result = stan_inference(prob, t_obs, data, priors; output_format=:dataframe)

if success(bayesian_result.return_code)
    df = bayesian_result.chains
    df |> display
end


text:
@testset "serialization" begin

    xx = "hey"; some_js = js"var"; x = [1f0, 2f0]

    js_str = js"console.log($xx); $x; $((2, 4)); $(some_js) hello = 1;"

    expect = "console.log('hey'); JSServe.deserialize_js({\"__javascript_type__\":\"TypedVector\",\"payload\":[1.0,2.0]}); JSServe.deserialize_js([2,4]); var hello = 1;"
    @test string(js_str) == expect

    asset = JSServe.Asset("file.dun_exist"; check_isfile=false)
    test_throw() = JSServe.include_asset(JSServe.Asset("file.dun_exist"))
    Test.@test_throws ErrorException("Unrecognized asset media type: dun_exist") test_throw()

    function test_handler(session, request)
        obs1 = Observable(Float16(2.0))
        obs2 = Observable(DOM.div("Data: ", obs1, dataTestId="hey"))
        return DOM.div(obs2)
    end
    testsession(test_handler, port=8555) do app
        hey = query_testid("hey")
        @test evaljs(app, js"$(hey).innerText") == "Data: 2.0"
        float16_obs = children(children(app.dom)[1][])[2]
        float16_obs[] = Float16(77)
        @test evaljs(app, js"$(hey).innerText") == "Data: 77"
    end
end

@testset "http" begin
    @test_throws ErrorException("Invalid sessionid: lol") JSServe.request_to_sessionid((target="lol",))
    @test JSServe.request_to_sessionid((target="lol",), throw=false) === nothing
end

@testset "hyperscript" begin
    function handler(session, request)
        the_script = DOM.script("window.testglobal = 42")
        s1 = Hyperscript.Style(css("p", fontWeight="bold"), css("span", color="red"))
        the_style = DOM.style(Hyperscript.styles(s1))
        return DOM.div(:hello, the_style, the_script, dataTestId="hello")
    end

    testsession(handler) do app
        @test evaljs(app, js"window.testglobal")  == 42
        hello_div = query_testid("hello")
        @test evaljs(app, js"$(hello_div).innerText")  == "hello"
        @test evaljs(app, js"$(hello_div).children.length") == 3
        @test evaljs(app, js"$(hello_div).children[0].tagName") == "P"
        @test evaljs(app, js"$(hello_div).children[1].tagName") == "STYLE"
        @test evaljs(app, js"$(hello_div).children[2].tagName") == "SCRIPT"
    end
end

@testset "async messages" begin
    obs = Observable(0); counter = Observable(0)
    testing_started = Ref(false)
    function handler(session, request)
        # Dont start this!
        testing_started[] && return DOM.div()
        onjs(session, obs, js"""function (v) {
            var t = JSServe.update_obs($(counter), JSServe.get_observable($(counter)) + 1);
        }""")

        for i in 1:2
            obs[] += 1
        end
        @async begin
            yield()
            for i in 1:2
                obs[] += 1
                yield()
            end
        end
        @async begin
            yield()
            for i in 1:2
                obs[] += 1
                yield()
            end
        end
        return DOM.div(obs, counter)
    end
    # Ugh, ElectronTests loads the handler multiple times to make sure it works
    # and doesn't get stuck, so we need to do this manually
    @isdefined(app) && close(app)
    app = JSServe.Server(handler, "0.0.0.0", 8558)
    try
        eapp = Electron.Application()
        window = Electron.Window(eapp)
        try
            @test obs[] == 0
            @test counter[] == 0
            testing_started[] = true
            Electron.load(window, URI(string("http://localhost:", 8558)))
            @wait_for counter[] == obs[]
        finally
            close(eapp)
        end
    finally
        close(app)
    end
end

@testset "Dependencies" begin
    jss = js"""function (v) {
        console.log($(JSTest));
    }"""
    div = DOM.div(onclick=jss)
    s = JSServe.Session()
    JSServe.register_resource!(s, div)
    @test JSTest.assets[1] in s.dependencies
end


text:
Servers.Server(j::Job) = Server(j.server)

function Jobs.Job(dir::AbstractString, s = "localhost"; version::Int = -1)
    server = Servers.maybe_start_server(s)
    if !ispath(server, dir)
        dir = request_job_dir(dir, server)
        dir === nothing && return
    end
    if version == last_version(dir, server)
        dir = Jobs.main_job_dir(dir)
    elseif !occursin(Jobs.VERSION_DIR_NAME, dir) && version != -1
        dir = Jobs.version_dir(Jobs.main_job_dir(dir), version)
    end
    resp = HTTP.get(server, "/jobs/" * dir)
    # Supplied dir was not a valid path, so we ask
    # previously registered jobs on the server that
    # contain dir.
    job = JSON3.read(resp.body, Job)
    job.server = server.name
    # if haskey(job.metadata, :timestamp)
    #     job.metadata[:timestamp] = DateTime(job.metadata[:timestamp])
    # end
    return job
end

function request_job_dir(dir::String, server::Server)
    resp = HTTP.get(server, "/registered_jobs/" * dir)
    matching_jobs = reverse(JSON3.read(resp.body, Vector{Tuple{String,DateTime}}))
    choices = ["$j -- $t" for (j, t) in matching_jobs]
    if length(matching_jobs) == 1
        return matching_jobs[1][1]
    elseif length(matching_jobs) == 0
        error("No jobs found matching $dir")
    elseif isdefined(Base, :active_repl)
        menu = RadioMenu(choices)
        choice = request("Multiple matching jobs were found, choose one:", menu)
        if choice != -1
            return matching_jobs[choice][1]
        else
            return nothing
        end
    else
        err_msg = "No concrete job for $dir found, closest matches are:"
        for m in choices
            err_msg = join([err_msg, m], "\n")
        end
        error(err_msg)
    end
end

function save(job::Job, server::Server = Servers.maybe_start_server(job))
    # First we check whether the job is trying to be saved in a archived directory, absolutely not allowed
    @assert !isrunning(job, server) "Can't save a job in a directory where another is running."

    @assert !Jobs.isarchived(job)
    "Not allowed to save a job in a archived directory, please specify a different directory."

    if isempty(job.name)
        @warn "Job had no name, changed it to: noname"
        job.name = "noname"
    end
    Structures.sanitize!(job.structure)
    Calculations.sanitize_flags!(job.calculations, job.structure, job.name,
                                 "./"*Jobs.TEMP_CALC_DIR)

    Jobs.sanitize_cutoffs!(job)

    curver = job.version

    tmpdir = mkpath(tempname())
    apath = abspath(job)

    if job.environment != ""
        environment = get_environment(job.environment, job.server)
        @assert environment !== nothing "Environment with name $(job.environment) not found!"
    else
        environment = nothing
    end
    job.dir = tmpdir
    FileIO.write_job_files(job, environment)
    files_to_send = Dict([f => read(joinpath(tmpdir, f), String) for f in readdir(tmpdir)])

    job.dir = apath
    rm(tmpdir, recursive=true)

    resp_job_version = JSON3.read(HTTP.post(server, "/jobs/" * apath, [], JSON3.write(files_to_send)).body,
                          Int)
    @info "Job version: $(curver) => $(resp_job_version)."
    job.version = resp_job_version
    Calculations.rm_tmp_flags!.(job.calculations)
    return job
end

function submit(job::Job)
    server = Servers.maybe_start_server(job)
    verify_execs(job, server)
    save(job)
    return HTTP.put(server, "/jobs/" * abspath(job))
end

function submit(jobs::Vector{Job}, run = true)
    # To verify all execs
    server_names = unique(map(x -> x.server, jobs))
    buckets = [jobs[findall(x -> x.server == s, jobs)] for s in server_names]
    outbuckets = [Vector{Job}(undef, length(b)) for b in buckets]
    Threads.@threads for i in 1:length(server_names)
        server = Servers.maybe_start_server(server_names[i])
        bucket = buckets[i]
        outbucket = outbuckets[i]
        execs = unique(vcat([map(x->x.exec, j.calculations) for j in bucket]...))

        replacements = verify_execs(execs, server)
        Threads.@threads for ij in 1:length(bucket)
            job = bucket[ij]
            for c in job.calculations
                for (e, rep) in replacements
                    if c.exec == e
                        c.exec.dir = rep.dir
                        c.exec.modules = rep.modules
                    end
                end
            end
            outbucket[ij] = save(job, server)
            run && HTTP.put(server, "/jobs/" * abspath(job))
        end
    end
    return outbuckets
end

"""
    state(job::Job, [server::Server])

Returns the state of a job.
"""
function state(job::Job, server::Server = Servers.maybe_start_server(job))
    return JSON3.read(HTTP.get(server, "/job_state/" * abspath(job)).body, Jobs.JobState)
end

"""
    isrunning(job::Job, [server::Server])

Returns whether a job is running or not. If the job was
submitted using `slurm`, a `QUEUED` status also counts as
running.
"""
function isrunning(job::Job, server::Server = Servers.maybe_start_server(job))
    s = state(job, server)
    return s == Jobs.Running || s == Jobs.Pending
end

"""
    versions(job::Job)

Returs the valid versions of `job`.
"""
versions(job::Job) = versions(abspath(job), job.server)
function versions(dir, s = "localhost")
    server=Servers.maybe_start_server(s)
    return JSON3.read(HTTP.get(server, "/job_versions/" * Jobs.main_job_dir(dir)).body, Vector{Int})
end

"""
    last_version(job::Job)

Returns the last version number of `job`.
"""
last_version(job::Job) = last_version(abspath(job), job.server)
function last_version(dir, s="localhost")
    t = versions(dir, s)
    return isempty(t) ? 0 : t[end]
end

function submission_time(job::Job)
    server = Servers.maybe_start_server(job)
    resp = HTTP.get(server, "/job_submission_time/" * abspath(job))
    return JSON3.read(resp.body, Float64)
end

"""
    last_running_calculation(job::Job)

Returns the last `Calculation` for which an output file was created.
"""
function last_running_calculation(job::Job)
    server = Servers.maybe_start_server(job)
    resp = HTTP.get(server, "/last_running_calculation/" * abspath(job))
    if resp.status == 204
        return nothing
    else
        return job[JSON3.read(resp.body, Int)]
    end
end

outputdata(job::Job, calcs::Vector{String}=String[]; kwargs...) =
    outputdata(job.dir, job.server, calcs; kwargs...)

function outputdata(jobdir::String, s::String = "localhost", calcs::Vector{String}=String[]; extra_parse_funcs = nothing)
    server = Servers.maybe_start_server(s)
    jobdir = isabspath(jobdir) ? jobdir : joinpath(server, jobdir)
    resp = HTTP.get(server, "/outputdata/" * jobdir, [], JSON3.write(calcs))
    if resp.status == 204
        error("No outputdata found yet. Is the job running?")
    end
    tmp_path = JSON3.read(resp.body,
                          String)

    local_temp = tempname() * ".jld2"
    Servers.pull(server, tmp_path, local_temp)
    dat = JLD2.load(local_temp)
    rm(local_temp)
    out = dat["outputdata"]
    if extra_parse_funcs !== nothing
        for k in keys(out)
            if k ∈ calcs
                n = c.name
                try
                    f = joinpath(jobdir, c.outfile)
                    local_f = tempname()
                    Servers.pull(server, f, local_f)
                    FileIO.parse_file(local_f, extra_parse_funcs, out = out[n])
                    rm(local_f)
                catch
                    nothing
                end
            end
        end
    end
    return out
end

function known_execs(e::String, dir::String = ""; server = Server("localhost"))
    s = Servers.maybe_start_server(server)
    return JSON3.read(HTTP.get(s, "/known_execs/", [], JSON3.write(Dict("exec" => e, "dir" => dir))).body, Vector{Calculations.Exec})
end
known_execs(e::Calculations.Exec; kwargs...) = known_execs(e.exec,e.dir; kwargs...)

function get_exec(name::String; server=Server("localhost"))
    s = Servers.maybe_start_server(server)
    return JSON3.read(HTTP.get(s, "/exec/$name").body, Calculations.Exec)
end
function save(e::Exec; server=Server("localhost"))
    s = Servers.maybe_start_server(server)
    return JSON3.read(HTTP.post(s, "/exec/", [], JSON3.write(e)).body, Calculations.Exec)
end

function verify_execs(job::Job, server::Server)
    replacements = verify_execs(unique(map(x->x.exec, job.calculations)), server)
    for (e, rep) in replacements
        for c in job.calculations
            if c.exec == e
                c.exec.dir = rep.dir
                c.exec.modules = rep.modules
            end
        end
    end
end

function verify_execs(execs::Vector{Exec}, server::Server)
    replacements = Dict{Exec, Exec}()
    for e in execs
        if !JSON3.read(HTTP.get(server, "/verify_exec/", [], JSON3.write(e)).body, Bool)
            possibilities = known_execs(e, server=server)
            curn = 0
            while length(possibilities) > 1
                possibilities = filter(x -> all(splitpath(x.dir)[end-curn:end] .== splitpath(e.dir)[end-curn:end]), possibilities)
                curn += 1
            end
            if !isempty(possibilities)
                replacement = possibilities[1]
                @warn "Executable ($(e.exec)) in dir ($(e.dir)) not runnable,\n but found a matching replacement executable in dir ($(replacement.dir)).\nUsing that one..."
                replacements[e] = replacement
            else
                error("$e is not a valid executable on server $(server.name).\nReplace it.")
            end
        end
    end
    return replacements
end

#TODO: work this
"""
    abort(job::Job)

Will try to remove the job from the scheduler's queue.
If the last running calculation happened to be a `Calculation{QE}`, the correct abort file will be written.
For other codes the process is not smooth, and restarting is not guaranteed.
"""
function abort(job::Job)
    lastrunning = job.calculations[last_running_calculation(job)]
    if lastrunning == nothing
        error("Is this job running?")
    end
    if eltype(lastrunning) == QE
        length(filter(x -> eltype(x) == QE, job.calculations)) > 1 &&
            @warn "It's absolutely impossible to guarantee a graceful abort of a multi job script with QE."

        abortpath = writeabortfile(job, lastrunning)
        while ispath(abortpath)
            continue
        end
        qdel(job)
    else
        qdel(job)
    end
end

function registered_jobs(fuzzy::String, server::Server)
    server = Servers.maybe_start_server(server)
    resp = HTTP.get(server, "/registered_jobs/" * fuzzy)
    return reverse(JSON3.read(resp.body, Vector{Tuple{String,DateTime}}))
end

registered_jobs(fuzzy::String, server::String) =
    registered_jobs(fuzzy, Servers.Server(server))

function registered_jobs(fuzzy::String = "")
    all_jobs = Dict{String,Vector{Tuple{String, DateTime}}}()
    for s in Servers.known_servers()
        jobs = registered_jobs(fuzzy, s)
        if !isempty(jobs)
            all_jobs[s.name] = jobs
        end
    end
    return all_jobs
end

function running_jobs(fuzzy="", s="localhost")
    server = Servers.maybe_start_server(s)
    resp = HTTP.get(server, "/running_jobs/" * fuzzy)
    return reverse(JSON3.read(resp.body, Vector{String}))
end

function switch_version!(job::Job, version::Int)
    allvers = versions(job)
    if !(version in allvers)
        error("Version $version does not exist.")
    end
    tj = Job(Jobs.main_job_dir(job), job.server, version=version)
    for f in fieldnames(Job)
        setfield!(job, f, getfield(tj, f))
    end
    return job
end

function rm_version!(job::Job, version::Int)
    server = Servers.maybe_start_server(job.server)
    allvers = versions(job)
    if !(version in allvers)
        error("Version $version does not exist.")
    end

    HTTP.put(server, "/rm_version/$version", [], JSON3.write(job))
    if version == job.version
        @warn "Job version is the same as the one to be removed, switching to last known version."
        lv = last_version(job)
        if lv == version
            lv = versions(job)[end-1]
        end
        switch_version!(job, lv)
    end
end

function environment_from_jobscript(scriptpath, s="localhost")
    server = Servers.maybe_start_server(s)
    tmp = tempname()
    Servers.pull(server, scriptpath, tmp)
    return Jobs.environment_from_jobscript(tmp)
end

function get_environment(name, s="localhost")
    server = Servers.maybe_start_server(s)
    return JSON3.read(HTTP.get(server, "/environment/$name").body, Environment)
end

function add_environment(env::Environment, name::String, s="localhost")
    server = Servers.maybe_start_server(s)
    return HTTP.post(server, "/environment/$name", [], JSON3.write(env))
end

function rm_environment!(name::String, s="localhost")
    server = Servers.maybe_start_server(s)
    return HTTP.put(server, "/environment/$name")
end


text:
"""
The pinched hemisphere benchmark. Half over fully closed sphere.

"""
module hemisphere_examples

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json

function _execute(n = 2, visualize = true)
    E = 6.825e7;
    nu = 0.3;
    thickness  =  0.04;
    # analytical solution for the vertical deflection under the load
    analyt_sol = 0.0924;
    R = 10.0;
    formul = FEMMShellT3FFModule

    tolerance = R/n/1000
    fens, fes = Q4spheren(R, n)
    fens, fes = Q4toT3(fens, fes)

    spherical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
        r = vec(XYZ);
        csmatout[:, 3] .= vec(r)/norm(vec(r))
        csmatout[:, 2] .= (0.0, 0.0, 1.0)
        cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
        csmatout[:, 1] .= vec(view(csmatout, :, 1))/norm(vec(view(csmatout, :, 1)))
        cross3!(view(csmatout, :, 2), view(csmatout, :, 3), view(csmatout, :, 1))
        return csmatout
    end
    ocsys = CSys(3, 3, spherical!)

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    # Report
    @info "Mesh: $n elements per side"

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    femm.drilling_stiffness_scale = 0.1
    femm.threshold_angle = 45.0
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,5,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    # top
    l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    nl = selectnode(fens; box = Float64[0 0 R R 0 0], tolerance = tolerance)
    loadbdry = FESetP1(reshape(nl, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
    fi = ForceIntensity(FFlt[0, -1, 0, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);
    nl = selectnode(fens; box = Float64[R R 0 0 0 0], tolerance = tolerance)
    loadbdry = FESetP1(reshape(nl, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
    fi = ForceIntensity(FFlt[1, 0, 0, 0, 0, 0]);
    F += distribloads(lfemm, geom0, dchi, fi, 3);

    # @infiltrate
    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    resultpercent =  dchi.values[nl, 1][1]*100
    @info "Solution: $(round(resultpercent/analyt_sol, digits = 4))%"

    # formul._resultant_check(femm, geom0, u0, Rfield0, dchi)

    # Visualization
    if !visualize
        return true
    end
    scattersysvec!(dchi, (R/4)/maximum(abs.(U)).*U)
    update_rotation_field!(Rfield0, dchi)
    plots = cat(plot_space_box([[0 0 -R]; [R R R]]),
        #plot_nodes(fens),
        plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    dims = 1)
    pl = render(plots)
    return true
end

function test_convergence()
    @info "Hemisphere benchmark"
    # for n in [2, ]
    for n in [2, 4, 8, 16, 32]
        _execute(n, false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
include(joinpath(@__DIR__, "pedagogy_graph.jl"))
using Random

dims = [5, 5]
g = LG.grid(dims)
β = 0.5

pp = PRS.HardCoreGraph(g, β)
rng = Random.MersenneTwister(123)

hard_core_graph_pedagogy(pp, dims, rng)


text:
module ContinuousTimePolicyGradients
    using LinearAlgebra, Statistics
    using OrdinaryDiffEq, DiffEqFlux, GalacticOptim
    using UnPack, Plots

    export CTPG_train, view_result

    include("construct_CTPG.jl")
end


text:
# Magnetic pendulum basins of attraction
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

α = 0.2; ω = 1.0; d = 0.3
gx = gy = range(-5, 5; length = 1500)
config = @strdict(α, ω, d, gx, gy)

function magnetic_basins(config)
    @unpack α, ω, d, gx, gy = config
    ma = Systems.magnetic_pendulum(; α, ω, d)
    if haskey(config, "attractors")
        attractors = config["attractors"]
        @time basins, attractors = basins_of_attraction((gx, gy), ma;
        mx_chk_att = 1, mx_chk_fnd_att = 10, attractors)
    else
        @time basins, attractors = basins_of_attraction((gx, gy), ma;
        mx_chk_att = 1, mx_chk_fnd_att = 10, mx_chk_lost = 10)
    end
    return @strdict(gx, gy, basins, attractors)
end

# produce high quality
file, _ = produce_or_load(datadir(), config, magnetic_basins; prefix = "magnetic_basins", tag = false)
attractors = file["attractors"]

# produce zoomed version
gx = range(1.80, 1.95; length = 1000)
gy = range(0, 0.12; length = 1000)
config = @strdict(α, ω, d, gx, gy, attractors)
produce_or_load(datadir(), config, magnetic_basins; prefix = "magnetic_basins_zoomed", tag = false)

# %% Plot this
@unpack gx, gy, basins, attractors = wload(datadir(savename("magnetic_basins", config, "jld2")))
ma = Systems.magnetic_pendulum(; α, ω, d)

LC =  matplotlib.colors.ListedColormap
cmap = LC([matplotlib.colors.to_rgb("C$k") for k in 0:2])

fig = figure(figsize=(0.33figx, figy))
ax = gca()
ax.pcolormesh(gx, gy, basins'; cmap, shading = "gouraud")
ax.set_aspect("auto")
xticks([-5, 5])
yticks([-5, 5])
xlabel("\$x\$", labelpad=-32)
ylabel("\$y\$", labelpad=-32)
for m in values(attractors)
    m1, m2 = columns(m)
    scatter(m1, m2; color = "white", edgecolor = "black", zorder = 99, s = 100)
end
# wsave(plotsdir("magneticpendulum"), fig)

# %% Plot zoomed version
@unpack gx, gy, basins, attractors = wload(datadir(savename("magnetic_basins_zoomed", config, "jld2")))
# fig = figure(figsize=(figx/2, figx/2))
x0 = 0.55; wx = 0.42; y0 = 0.65; wy = 0.32;
axin = ax.inset_axes([x0, y0, wx, wy])

axin.pcolormesh(gx, gy, basins'; cmap, shading = "gouraud")
# axin.set_aspect("equal")

# xticks([1.8, 1.95], size = 20)
# yticks([0, 0.12], size = 20)
zbox = ((gx[1], gy[1]), (gx[end], gy[end]))
axis_zoomin!(axin, ax, zbox, zbox, "C3"; dir = :top, lw = 4)

axin.axis("off")

fig.tight_layout(pad=0.3)

wsave(plotsdir("5", "magneticpendulum.png"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
######
#
# 2-dimensional stuff.
#

function regulargrid2d(box, res)
    xmin, ymin, xmax, ymax = box
    rx, ry = res
    dx = (xmax-xmin)/(rx-1)
    dy = (ymax-ymin)/(ry-1)

    vs_cnt = rx*ry # vertices count
    es_cnt = (rx-1)*ry + rx*(ry-1)+ (rx-1)*(rx-1) # Horizontal + vertical + diagonal
    fs_cnt = 2*(rx-1)*(ry-1) # down + up

    ps = zeros(vs_cnt, 3) # points
    vs = zeros(Int,vs_cnt,1) # vertices
    es = zeros(Int,es_cnt,2) # edges
    fs = zeros(Int,fs_cnt,3) # faces

    # rows, columns and values for edge-vertex adjacency
    e_rows = zeros(Int, 2*es_cnt)
    e_cols = zeros(Int, 2*es_cnt)
    e_vals = zeros(Int, 2*es_cnt)

    # rows, columns and values for face-edge adjacency
    f_rows = zeros(Int, 3*fs_cnt)
    f_cols = zeros(Int, 3*fs_cnt)
    f_vals = zeros(Int, 3*fs_cnt)

    esc = 1
    fsc = 1
    for r in 0:(ry-1)
        for c in 0:(rx-1)
            # vertices indices
            v00 = r*rx+ c
            v01 = v00 + 1
            v10 = v00 + rx
            v11 = v01 + rx
            # horizontal edges
            he00 = r*(rx-1) + c
            he01 = he00 + 1
            he10 = he00 + (rx-1)
            he11 = he01 + (rx-1)
            # vertical edges
            ve00 = (rx-1)*ry + r*rx + c
            ve01 = ve00 + 1
            ve10 = ve00 + rx
            ve11 = ve01 + rx
            # diagonal edges
            de00 = (rx-1)*ry + rx*(ry-1) + r*(rx-1) + c
            de01 = de00 + 1
            de10 = de00 + (rx-1)
            de11 = de01 + (rx-1)
            # down faces
            df00 = r*(rx-1) + c
            df01 = df00 + 1
            df10 = df00 + (rx-1)
            df11 = df01 + (rx-1)
            # up faces
            uf00 = (rx-1)*(ry-1) + r*(rx-1) + c
            uf01 = uf00 + 1
            uf10 = uf00 + (rx-1)
            uf11 = uf01 + (rx-1)

            # setting points
            ps[v00+1,:] = [xmin + c*dx, ymin + r*dy, 0.0]
            # setting vertices
            vs[v00+1] = v00
            # Setting edges and faces
            if c < (rx-1) # horizontal edges
                es[he00+1,:] = [v00,v01]
                e_rows[esc] = he00; e_cols[esc] = v00; e_vals[esc] = -1; esc = esc + 1
                e_rows[esc] = he00; e_cols[esc] = v01; e_vals[esc] = 1; esc = esc + 1
            end
            if r < (ry-1) # vertical edges
                es[ve00+1,:] = [v00,v10]
                e_rows[esc] = ve00; e_cols[esc] = v00; e_vals[esc] = -1; esc = esc + 1
                e_rows[esc] = ve00; e_cols[esc] = v10; e_vals[esc] = 1; esc = esc + 1
            end
            if r <(ry-1) && c < (rx-1)
                fs[df00+1,:] = [v00,v01,v11] # down faces
                f_rows[fsc] = df00; f_cols[fsc] = he00; f_vals[fsc] = 1; fsc = fsc + 1
                f_rows[fsc] = df00; f_cols[fsc] = ve01; f_vals[fsc] = 1; fsc = fsc + 1
                f_rows[fsc] = df00; f_cols[fsc] = de00; f_vals[fsc] = -1; fsc = fsc + 1

                fs[uf00+1,:] = [v11,v10,v00] # up faces
                f_rows[fsc] = uf00; f_cols[fsc] = de00; f_vals[fsc] = 1; fsc = fsc + 1
                f_rows[fsc] = uf00; f_cols[fsc] = he10; f_vals[fsc] = -1; fsc = fsc + 1
                f_rows[fsc] = uf00; f_cols[fsc] = ve00; f_vals[fsc] = -1; fsc = fsc + 1

                es[de00+1,:] = [v00, v11]    # diagonal edges
                e_rows[esc] = de00; e_cols[esc] = v00; e_vals[esc] = -1; esc = esc + 1
                e_rows[esc] = de00; e_cols[esc] = v11; e_vals[esc] = 1; esc = esc + 1
            end
        end
    end
    ps, vs.+1, es.+1, fs.+1
end

######
#
# 3-dimensional stuff.
#
function regulargrid3d(box, res)
    xmin, ymin, zmin, xmax, ymax, zmax = box
    resx, resy, resz = res
    resxy = resx*resy
    resxyz = resxy*resz
    dx = (xmax-xmin)/(resx-1)
    dy = (ymax-ymin)/(resy-1)
    dz = (zmax-zmin)/(resz-1)
    ps = zeros(resxyz, 3)
    vs = zeros(Int,resxyz,1)
    ts = zeros(Int,6*(resx-1)*(resy-1)*(resz-1),4)
    tc = 1
    for k in 0:(resz-1)
        for j in 0:(resy-1)
            for i in 0:(resx-1)
                id = resxy*k+ resx*j + i
                vs[id+1] = id
                ps[id+1,:] = [xmin+i*dx,ymin+j*dy,zmin+k*dz]
                v0 = resxy*k + j*resx + i
                v1 = resxy*k + j*resx + i + 1
                v2 = resxy*k + (j+1)*resx + i+1
                v3 = resxy*k + (j+1)*resx + i

                v4 = resxy*(k+1) + j*resx + i
                v5 = resxy*(k+1) + j*resx + i + 1
                v6 = resxy*(k+1) + (j+1)*resx + i + 1
                v7 = resxy*(k+1) + (j+1)*resx + i
                if i <(resx-1) && j < (resy-1) && k < (resz-1)
                    ts[tc+0,:] = [v2,v0,v3,v7]
                    ts[tc+1,:] = [v0,v2,v6,v7]
                    ts[tc+2,:] = [v4,v0,v6,v7]
                    ts[tc+3,:] = [v6,v0,v1,v2]
                    ts[tc+4,:] = [v0,v6,v1,v4]
                    ts[tc+5,:] = [v6,v5,v1,v4]
                    tc = tc + 6
                end
            end
        end
    end
    ps, ts.+1
end

function compute_subfaces(faces::Array{Int,2})
    fs_cnt, k = size(faces)
    mask = ones(Bool,k)
    sfs = Array{Array{Int64,2}}(undef,k*fs_cnt)
    for i in 1:fs_cnt
        for j in 1:k
            mask[j] = false
            sfs_index = (i-1)*k + j
            sfs[sfs_index] = sort(faces[[i],mask],dims=2)
            mask[j] = true
        end
    end
    subfaces = vcat(unique(sfs)...)
    subfaces
end


text:
# continuous 1-d; different scalings should lead to different sizes, different boundary condtions shouldn't
@testset "Construction Types" begin
    xSizes = (33, 67, 128)
    boundaries = (DEFAULT_BOUNDARY, padded, NaivePer)
    sVals = [1, 2, 3.5, 8, 16]
    βs = [1, 1.5,4.0]
    waveTypes =(ContinuousWavelets.morl, dog0, paul4, cDb2, ContinuousWavelets.morse)
    @testset "xSz=$xSize, b=$boundary, s=$s, β=$β, wfc=$(wave)" for xSize in xSizes, boundary in boundaries, s in sVals, β in βs, wave in waveTypes
        wfc = ContinuousWavelets.wavelet(wave, s=s, boundary=boundary,β=β)
        xr = rand(Float64, xSize);
        x32 = rand(Float32, xSize);
        xc = rand(ComplexF64, xSize);
        # the sizes are of course broken at this size, so no warnings needed
        yc = 3; yr = 3; y32 = 3
        with_logger(ConsoleLogger(stderr, Logging.Error)) do
            yr = ContinuousWavelets.cwt(xr, wfc);
            y32 = ContinuousWavelets.cwt(x32, wfc);
            yc = ContinuousWavelets.cwt(xc, wfc);
        end
        if isAnalytic(wfc.waveType)
            @test eltype(yr) == ComplexF64
            @test eltype(y32) == ComplexF32
            @test eltype(yc) == ComplexF64
        else
            @test eltype(yr) == Float64
            @test eltype(y32) == Float32
            @test eltype(yc) == ComplexF64
        end
        nOctaves, totalWavelets, sRanges, sWidths = ContinuousWavelets.getNWavelets(xSize, wfc);
        @test totalWavelets == size(yc, 2)
    end
end
# xSize = xSizes[3]; boundary = boundaries[3]; s = sVals[end]; β = βs[3]; wave = waveTypes[4]
# xSize = 33; boundary = SymBoundary(); s = 1; β = 1; wave = dog2