text:
# Paths in the representation domain of the BZ, obtained from the [SeeK] publication.
# Methods return **k**-vector coordinates in a primitive reciprocal basis ("reciprocal
# crystallographic primitive cell" in the [SeeK]'s nomenclature).
#
# ^[Seek] http://dx.doi.org/10.1016/j.commatsci.2016.10.015; see also online interface at
#         https://www.materialscloud.org/work/tools/seekpath

module KPaths

# ---------------------------------------------------------------------------------------- #
export irrfbz_path, KPath, points, paths, KPathInterpolant
# ---------------------------------------------------------------------------------------- #
using Bravais:
    Bravais, # so it is possible to call Brillouin.Bravais from the "outside"
    bravaistype,
    boundscheck_sgnum,
    reciprocalbasis,
    primitivize,
    transform,
    ReciprocalBasis,
    DirectBasis,
    primitivebasismatrix
using ..Brillouin:
    AVec,
    BasisLike,
    SHOWDIGITS,
    BasisEnum, CARTESIAN, LATTICE, setting, set_setting!,
    cartesianize, latticize
import ..Brillouin:
    latticize,
    latticize!,
    cartesianize!,
    basis
using LinearAlgebra:
    norm,
    dot,
    ×
using StaticArrays
using DocStringExtensions

import Base: show, summary, getindex, setindex!, IndexStyle, size
# ---------------------------------------------------------------------------------------- #

include("codegen-kpoints.jl")  # defines `get_points_3d`, `pathsd_3d` (& other subfunctions)
include("bravais-branches.jl") # defines `extended_bravais`

# ---------------------------------------------------------------------------------------- #

abstract type AbstractPath{T} <: AbstractVector{T} end

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct KPath{D} <: AbstractPath{Pair{Symbol, SVector{D, Float64}}}
    points  :: Dict{Symbol, SVector{D,Float64}}
    # TODO: Make values of `points` of type `ReciprocalPoint{D}`?
    paths   :: Vector{Vector{Symbol}}
    basis   :: ReciprocalBasis{D}
    setting :: Ref{BasisEnum}
end

"""
    points(kp::KPath{D}) -> Dict{Symbol, SVector{D,Float64}}

Return a dictionary of the **k**-points (values) and associated **k**-labels (keys)
referenced in `kp`.
"""
points(kp::KPath) = kp.points

"""
    paths(kp::KPath) -> Vector{Vector{Symbol}}

Return a vector of vectors, with each vector describing a connected path between between
**k**-points referenced in `kp` (see also [`points(::KPath)`](@ref)).
"""
paths(kp::KPath) = kp.paths

Base.@propagate_inbounds function getindex(kp::KPath, i::Int)
    # index into the `i`th point in the "flattened" `paths(kp)`
    @boundscheck i < 0 && throw(BoundsError(kp, i))
    j = 1
    i′ = i′′ = i
    while (i′′ -= length(paths(kp)[j])) > 0
        i′ = i′′
        j += 1
        @boundscheck j > length(paths(kp)) && throw(BoundsError(kp, i))
    end
    @boundscheck i′ > length(paths(kp)[j]) && throw(BoundsError(kp, i))

    lab = @inbounds paths(kp)[j][i′]
    return lab => points(kp)[lab]
end
IndexStyle(::Type{<:KPath}) = IndexLinear()
size(kp::KPath) = (sum(length, paths(kp)),)

function summary(io::IO, kp::KPath{D}) where D
    print(io, "KPath{", D, "} (", length(points(kp)), " points, ",
              length(paths(kp)), " paths, ", length(kp), " points in paths)")
end
function show(io::IO, ::MIME"text/plain", kp::KPath)
    summary(io, kp)
    println(io, ":")
    # points
    print(io, " points: ")
    foreach(enumerate(points(kp))) do (i,(lab,kv))
        i ≠ 1 && print(io, "         ")
        println(io, ":", lab, " => ", round.(kv, digits=SHOWDIGITS))
    end
    # paths
    print(io, "  paths: ")
    foreach(enumerate(paths(kp))) do (i,path)
        i ≠ 1 && print(io, "         ")
        println(io, path)
    end
    # basis
    print(io, "  basis: ")
    foreach(enumerate(basis(kp))) do (i,v)
        i ≠ 1 && print(io, "         ")
        if i ≠ length(basis(kp))
            println(io, round.(v, digits=SHOWDIGITS))
        else
            print(io, round.(v, digits=SHOWDIGITS))
        end
    end
end
# ---------------------------------------------------------------------------------------- #

"""
    irrfbz_path(sgnum::Integer, Rs, [::Union{Val(D), Integer},]=Val(3))  -->  ::KPath{D}

Returns a **k**-path (`::KPath`) in the (primitive) irreducible Brillouin zone for a space
group with number `sgnum`, (conventional) direct lattice vectors `Rs`, and dimension `D`.
The path includes all distinct high-symmetry lines and points as well as relevant parts of
the Brillouin zone boundary.

The dimension `D` (1, 2, or 3) is specified as the third input argument, preferably as a
static `Val{D}` type parameter (or, type-unstably, as an `<:Integer`). Defaults to `Val(3)`.

`Rs` refers to the direct basis of the conventional unit cell, i.e., not the primitive
direct basis vectors. The setting of `Rs` must agree with the conventional setting choices
in the International Tables of Crystallography, Volume A (the "ITA conventional setting").
If `Rs` is a subtype of a `StaticVector` or `NTuple`, the dimension can be inferred from its
(static) size; in this case, this dimension will take precedence (i.e. override, if
different) over any dimension specified in the third input argument.

## Notes
- The returned **k**-points are given in the basis of the **primitive** reciprocal basis in
  the CDML setting. To obtain the associated transformation matrices between the
  ITA conventional setting and the CDML primitive setting, see `primitivebasismatrix` of
  Bravais.jl](https://thchr.github.io/Crystalline.jl/stable/bravais/) (or, equivalently, the
  relations defined Table 2 of [^1]).
  To transform to a Cartesian basis, see [`cartesianize!`](@ref).
- To interpolate a `KPath`, see [`interpolate(::KPath, ::Integer)`](@ref) and
  [`splice(::KPath, ::Integer)`](@ref).
- All paths currently assume time-reversal symmetry (or, equivalently, inversion symmetry).
  If neither are present, include the "inverted" -**k** paths manually.

## Data and referencing
3D paths are sourced from the SeeK-path publication: please cite the original work
[^2].

## References
[^1] Aroyo et al., [Acta Cryst. A70, 126 (2014)](https://doi.org/10.1107/S205327331303091X).
[^2] Hinuma, Pizzi, Kumagai, Oba, & Tanaka, *Band structure diagram paths based on
    crystallography*, [Comp. Mat. Sci. **128**, 140 (2017)]
    (http://dx.doi.org/10.1016/j.commatsci.2016.10.015).
"""
function irrfbz_path(sgnum::Integer, Rs, Dᵛ::Val{D}=Val(3)) where D
    D′ = length(Rs)
    D′ ≠ D && throw(DimensionMismatch("inconsistent dimensions of `Rs` and `Dᵛ`"))
    any(R -> length(R) ≠ D, Rs) && throw(DimensionMismatch("inconsistent element dimensions in `Rs`"))
    Rs = convert(DirectBasis{D}, Rs)

    # (extended) bravais type
    bt = bravaistype(sgnum, D; normalize=false)
    ext_bt = extended_bravais(sgnum, bt, Rs, Dᵛ)

    # get data about path
    lab2kvs = get_points(ext_bt, Rs, Dᵛ)
    paths = get_paths(ext_bt, Dᵛ)

    # unshuffle HPKOT-standard primitive setting-difference in oA and mC settings
    unshuffle_hpkot_setting!(lab2kvs, bt, D)

    # compute (primitive) reciprocal basis
    cntr = last(bt)
    pGs = primitivize(reciprocalbasis(Rs), cntr)

    return KPath(lab2kvs, paths, pGs, Ref(LATTICE))
end
function irrfbz_path(sgnum::Integer,
            Rs::Union{StaticVector{D, <:AVec{<:Real}}, NTuple{D, <:AVec{<:Real}}}) where D
    return irrfbz_path(sgnum, Rs, Val(D))
end
irrfbz_path(sgnum::Integer, Rs::AVec{<:AVec{<:Real}}, D::Integer) = irrfbz_path(sgnum, Rs, Val(D))

# Note: output k-points are returned in the primitive setting defined in Table 3 of the
# HPKOT paper. This setting differs from the standard crystallographic primitive setting
# for 'oA' and 'mC' Bravais types; "unshuffling" this is done in `unshuffle_hpkot_setting!`
function get_points(ext_bt::Symbol,
                    Rs::Union{Nothing, AVec{<:AVec{<:Real}}},
                    Dᵛ::Val{D}) where D
    if Dᵛ === Val(3)
        return get_points_3d(ext_bt, Rs)
    elseif Dᵛ === Val(2)
        return get_points_2d(ext_bt, Rs)
    elseif Dᵛ === Val(1)
        ext_bt === :lp || throw(DomainError(ext_bt, "invalid extended Bravais type"))
        return Dict(:Γ => SVector(0.0), :X => SVector(0.5))
    else
        throw(DomainError(D, "unsupported input dimension"))
    end
end

function get_paths(ext_bt::Symbol, Dᵛ::Val{D}) where D
    if Dᵛ === Val(3)
        paths = get(pathsd_3d, ext_bt, nothing)
    elseif Dᵛ === Val(2)
        paths = get(pathsd_2d, ext_bt, nothing)
    elseif Dᵛ === Val(1)
        paths = ext_bt === :lp ? [[:Γ, :X]] : nothing
    else
        throw(DomainError(D, "unsupported input dimension"))
    end

    paths === nothing && throw(DomainError(ext_bt, "invalid extended Bravais type"))
    return paths
end

function unshuffle_hpkot_setting!(lab2kvs, bt, D)
    # the points in `labs2kvs` are returned in the reciprocal basis of the primitive HPKOT
    # setting (Table 3, HPKOT paper). That setting differs from the primitive CDML setting
    # that we follow (Table 2, https://doi.org/10.1107/S205327331303091X) for Bravais types
    # 'oA' and 'mC'. Hence, we need to shuffle it back to CDML primitive settings.
    # Specifically, HPKOT has centering matrix `P′` which differs from the CDML centering
    # matrix `P`. To cast the points back to the CDML primitive setting we first transform
    # the points to the conventional setting (by multiplying by `((P′)⁻¹)ᵀ`) and then
    # transform back to the CDML primitive setting (by  multiplying with `Pᵀ`):
    if D == 3 && (bt == "oA" || bt == "mC")
        cntr = last(bt)
        if bt == "oA"
            P′ = @SMatrix [0 0 1; 0.5 0.5 0; -0.5 0.5 0] # HPKOT centering matrix
        else # bt == "mC"
            P′ = @SMatrix [0.5 -0.5 0; 0.5 0.5 0; 0 0 1]
        end
        P = primitivebasismatrix(cntr, Val(3))           # CDML centering matrix
        PᵀP′⁻¹ᵀ = P'*inv(P′)'
        for (lab, kv′) in lab2kvs
            kv = PᵀP′⁻¹ᵀ*kv′  # k-vector coordinates transform w/ transpose of matrix (as
            lab2kvs[lab] = kv # explained in e.g. Bravais.jl's src/transform.jl)
        end
    end
    return lab2kvs
end

# ---------------------------------------------------------------------------------------- #

include("interpolate-paths.jl")
export interpolate, splice, cumdists

# ---------------------------------------------------------------------------------------- #

"""
    cartesianize!(kp::KPath{D})

Transform a **k**-path `kp` in a lattice basis to a Cartesian basis with (primitive)
reciprocal basis vectors `basis(kp)`. Modifies `kp` in-place.
"""
function cartesianize!(kp::KPath{D}) where D
    setting(kp) === CARTESIAN && return kp
    for (lab, kv) in points(kp)
        @inbounds points(kp)[lab] = cartesianize(kv, basis(kp))
    end
    set_setting!(kp, CARTESIAN)
    return kp
end

"""
    latticize!(kp::KPath)

Transform a **k**-path `kp` in a Cartesian basis to a lattice basis with (primitive)
reciprocal lattice vectors `basis(kp)`. Modifies `kp` in-place.
"""
function latticize!(kp::KPath)
    setting(kp) === LATTICE && return kp
    basismatrix = reduce(hcat, basis(kp))
    for (lab, kv) in points(kp)
        @inbounds points(kp)[lab] = latticize(kv, basismatrix)
    end
    set_setting!(kp, LATTICE)
    return kp
end

"""
    latticize(kp::KPath{D}, basis::AbstractVector{<:AbstractVector{<:Real}})

Transform a **k**-path `kp` in a Cartesian basis to a lattice basis with (primitive)
reciprocal lattice vectors `basis`.

If `kp` is not in a Cartesian basis (i.e., if `setting(kp) == LATTICE`), `kp` is returned
as-is.
"""
function latticize(kp::KPath{D}, basis::AVec{<:AVec{<:Real}}) where D
    setting(kp) === LATTICE && return kp
    basismatrix = convert(SMatrix{D,D,Float64,D*D}, reduce(hcat, basis))
    kvs′ = Dict{Symbol, SVector{D,Float64}}()
    for (lab, kv) in points(kp)
        @inbounds kvs′[lab] = latticize(kv, basismatrix)
    end
    return typeof(kp)(kvs′, paths(kp), basis, Ref(LATTICE))
end

# ---------------------------------------------------------------------------------------- #
end # module


text:
"""
    SkipDeepEquilibriumNetwork(model, shortcut, solver; p=nothing, jacobian_regularization::Bool=false,
                               sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)
    SkipDeepEquilibriumNetwork(model, solver; p=nothing, jacobian_regularization::Bool=false,
                               sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)

Skip Deep Equilibrium Network as proposed in [pal2022mixing](@cite)

## Arguments

* `model`: Explicit Neural Network which takes 2 inputs
* `shortcut`: Shortcut for the network (If not given, then we create SkipDEQV2)
* `solver`: Solver for the optimization problem (See: [`ContinuousDEQSolver`](@ref) & [`DiscreteDEQSolver`](@ref))
* `jacobian_regularization`: If true, Jacobian Loss is computed and stored in the [`DeepEquilibriumSolution`](@ref)
* `p`: Optional parameters for the `model`
* `sensealg`: See [`SteadyStateAdjoint`](@ref)
* `kwargs`: Additional Parameters that are directly passed to `solve`

## Example

```julia
# SkipDEQ
model = SkipDeepEquilibriumNetwork(
    Parallel(
        +,
        Dense(2, 2; bias=false),
        Dense(2, 2; bias=false)
    ),
    Dense(2, 2),
    ContinuousDEQSolver(VCABM3(); abstol=0.01f0, reltol=0.01f0)
)

model(rand(Float32, 2, 1))

# SkipDEQV2
model = SkipDeepEquilibriumNetwork(
    Parallel(
        +,
        Dense(2, 2; bias=false),
        Dense(2, 2; bias=false)
    ),
    ContinuousDEQSolver(VCABM3(); abstol=0.01f0, reltol=0.01f0)
)

model(rand(Float32, 2, 1))
```

See also: [`DeepEquilibriumNetwork`](@ref), [`MultiScaleDeepEquilibriumNetwork`](@ref), [`MultiScaleSkipDeepEquilibriumNetwork`](@ref)
"""
struct SkipDeepEquilibriumNetwork{M,S,J,P,RE1,RE2,A,Se,K} <: AbstractDeepEquilibriumNetwork
    jacobian_regularization::Bool
    model::M
    shortcut::S
    p::P
    re1::RE1
    re2::RE2
    split_idx::Int
    solver::A
    kwargs::K
    sensealg::Se
    stats::DEQTrainingStats

    function SkipDeepEquilibriumNetwork(jacobian_regularization, model, shortcut, p, re1, re2,
                                        split_idx, solver, kwargs, sensealg, stats)
        p1, re1 = destructure_parameters(model)
        split_idx = length(p1)
        p2, re2 = shortcut === nothing ? ((eltype(p1))[], nothing) : destructure_parameters(shortcut)

        p = p === nothing ? vcat(p1, p2) : eltype(p1).(p)

        return new{typeof(model),typeof(shortcut),jacobian_regularization,typeof(p),typeof(re1),
                   typeof(re2),typeof(solver),typeof(sensealg),typeof(kwargs)}(jacobian_regularization, model, shortcut, p,
                                                                               re1, re2, split_idx, solver, kwargs,
                                                                               sensealg, stats)
    end
end

Flux.@functor SkipDeepEquilibriumNetwork

function Base.show(io::IO, l::SkipDeepEquilibriumNetwork{M,S,J}) where {M,S,J}
    shortcut_ps = l.split_idx == length(l.p) ? 0 : length(l.p) - l.split_idx
    return print(io, "SkipDeepEquilibriumNetwork(jacobian_regularization = $J, ",
                 "shortcut_parameter_count = $shortcut_ps) ", string(length(l.p)), " Trainable Parameters")
end

function SkipDeepEquilibriumNetwork(model, shortcut, solver; p=nothing, jacobian_regularization::Bool=false,
                                    sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)
    return SkipDeepEquilibriumNetwork(jacobian_regularization, model, shortcut, p, nothing,
                                      nothing, 0, solver, kwargs, sensealg, DEQTrainingStats(0))
end

function SkipDeepEquilibriumNetwork(model, solver; p=nothing, jacobian_regularization::Bool=false,
                                    sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)
    return SkipDeepEquilibriumNetwork(jacobian_regularization, model, nothing, p, nothing,
                                      nothing, 0, solver, kwargs, sensealg, DEQTrainingStats(0))
end

function (deq::SkipDeepEquilibriumNetwork)(x::AbstractArray{T}) where {T}
    p1, p2 = deq.p[1:(deq.split_idx)], deq.p[(deq.split_idx + 1):end]
    z = deq.re2(p2)(x)::typeof(x)

    current_nfe = deq.stats.nfe

    # Dummy call to ensure that mask is generated
    Zygote.@ignore _ = deq.re1(p1)(z, x)

    z_star = solve_steady_state_problem(deq.re1, p1, x, z, deq.sensealg, deq.solver; dudt=nothing,
                                        update_nfe=() -> (deq.stats.nfe += 1), deq.kwargs...)

    jac_loss = (deq.jacobian_regularization ? compute_deq_jacobian_loss(deq.re1, p1, z_star, x) : T(0)) ::T

    residual = Zygote.@ignore z_star .- deq.re1(p1)(z_star, x)

    return z_star, DeepEquilibriumSolution(z_star, z, residual, jac_loss, deq.stats.nfe - current_nfe)
end

function (deq::SkipDeepEquilibriumNetwork{M,Nothing})(x::AbstractArray{T}) where {M,T}
    z = deq.re1(deq.p)(zero(x), x)::typeof(x)

    current_nfe = deq.stats.nfe

    z_star = solve_steady_state_problem(deq.re1, deq.p, x, z, deq.sensealg, deq.solver; dudt=nothing,
                                        update_nfe=() -> (deq.stats.nfe += 1), deq.kwargs...)

    jac_loss = (deq.jacobian_regularization ? compute_deq_jacobian_loss(deq.re1, deq.p, z_star, x) : T(0)) ::T

    residual = Zygote.@ignore z_star .- deq.re1(deq.p)(z_star, x)

    return z_star, DeepEquilibriumSolution(z_star, z, residual, jac_loss, deq.stats.nfe - current_nfe)
end

ext: jl


text:
"""
    Interpolate <: Imputor

Performs linear interpolation between the nearest values in an vector.
"""
struct Interpolate <: Imputor end

"""
    impute!(imp::Interpolate, ctx::Context, data::AbstractVector)

Uses linear interpolation between existing elements of a vector to fill in missing data.

WARNING: Missing values at the head or tail of the array cannot be interpolated if there
are no existing values on both sides. As a result, this method does not guarantee
that all missing values will be imputed.
"""
function impute!(imp::Interpolate, ctx::Context, data::AbstractVector{<:Union{T, Missing}}) where T
    i = findfirst(ctx, data) + 1

    while i < length(data)
        if ismissing(ctx, data[i])
            prev_idx = i - 1
            next_idx = findnext(ctx, data, i + 1)

            if next_idx !== nothing
                gap_sz = (next_idx - prev_idx) - 1

                diff = data[next_idx] - data[prev_idx]
                incr = diff / T(gap_sz + 1)
                start_val = data[prev_idx]
                stop_val = data[next_idx]

                values = Real(start_val):Real(incr):Real(stop_val)

                idx_range = prev_idx:(prev_idx + length(values) - 1)
                # println(collect(idx_range))
                # println(values)

                data[idx_range] = values
                i = next_idx
            else
                break
            end
        end
        i += 1
    end

    return data
end


text:
using FinEtools
using FinEtools.AlgoDeforLinearModule
using FinEtools.MeshExportModule

mu = 756000
nu = 0.49546;
E = 2*(1.0 + nu) * mu;
width = 48.0; height = 44.0; thickness  = 10.0;
free_height  = 16.0;
Mid_edge  = [48.0 52.0 0.0];# Location of tracked  deflection
magn = 1.0/(free_height*thickness);# Density of applied load
nt = 1
n = 16*nt;# number of elements per side
tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance

fens,fes = T10block(width, height, thickness, n, n, nt)

# Reshape into a trapezoidal panel
for i = 1:count(fens)
    fens.xyz[i,2] = fens.xyz[i,2]+(fens.xyz[i,1]/width)*(height -fens.xyz[i,2]/height*(height-free_height));
end

# Clamped edge of the membrane
l1 = selectnode(fens; box=[0.,0.,-Inf, Inf,-Inf, Inf], inflate = tolerance)
ess1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>l1)
ess2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>l1)
ess3 = FDataDict("displacement"=>  0.0, "component"=> 3, "node_list"=>l1)

# Traction on the opposite edge
boundaryfes =  meshboundary(fes);
Toplist  = selectelem(fens, boundaryfes, box= [width, width, -Inf, Inf, -Inf, Inf], inflate=  tolerance);
el1femm = FEMMBase(IntegDomain(subset(boundaryfes, Toplist), SimplexRule(2, 3)))
flux1 = FDataDict("traction_vector"=>[0.0, 0.0, +magn],
    "femm"=>el1femm
    )

# Make the region
MR = DeforModelRed3D
material = MatDeforElastIso(MR,  0.0, E, nu, 0.0)
region1 = FDataDict("femm"=>FEMMDeforLinearMST10(MR,
    IntegDomain(fes, SimplexRule(3, 2)), material))

modeldata = FDataDict("fens"=>fens,
 "regions"=>[region1],
 "essential_bcs"=>[ess1, ess2, ess3],
 "traction_bcs"=>[flux1]
 )

# Call the solver
modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

u = modeldata["u"]
geom = modeldata["geom"]

# Extract the solution
nl = selectnode(fens, box=[Mid_edge[1],Mid_edge[1],Mid_edge[2],Mid_edge[2], -Inf, Inf],
          inflate=tolerance);
theutip = mean(u.values[nl,:])
println("displacement =$(theutip[2]) as compared to converged $convutip")

modeldata["postprocessing"] = FDataDict("file"=>"cook_like-ew",
   "quantity"=>:Cauchy, "component"=>:xy)
modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
fld = modeldata["postprocessing"]["exported"][1]["field"]
println("range of Cauchy_xy = $((minimum(fld.values), maximum(fld.values)))")
File = modeldata["postprocessing"]["exported"][1]["file"]
@async run(`"paraview.exe" $File`)

modeldata["postprocessing"] = FDataDict("file"=>"cook_like-ew-vm",
   "quantity"=>:vm, "component"=>1)
modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
fld = modeldata["postprocessing"]["exported"][1]["field"]
println("range of vm = $((minimum(fld.values), maximum(fld.values)))")
File = modeldata["postprocessing"]["exported"][1]["file"]
@async run(`"paraview.exe" $File`)

modeldata["postprocessing"] = FDataDict("file"=>"cook_like-ew-pressure",
   "quantity"=>:pressure, "component"=>1)
modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
fld = modeldata["postprocessing"]["exported"][1]["field"]
println("range of pressure = $((minimum(fld.values), maximum(fld.values)))")
File = modeldata["postprocessing"]["exported"][1]["file"]
@async run(`"paraview.exe" $File`)

modeldata["postprocessing"] = FDataDict("file"=>"cook_like-ew-princ1",
   "quantity"=>:princCauchy, "component"=>1)
modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
fld = modeldata["postprocessing"]["exported"][1]["field"]
println("range of princCauchy Max = $((minimum(fld.values), maximum(fld.values)))")
File = modeldata["postprocessing"]["exported"][1]["file"]
@async run(`"paraview.exe" $File`)

modeldata["postprocessing"] = FDataDict("file"=>"cook_like-ew-princ3",
   "quantity"=>:princCauchy, "component"=>3)
modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
fld = modeldata["postprocessing"]["exported"][1]["field"]
println("range of princCauchy Min = $((minimum(fld.values), maximum(fld.values)))")
File = modeldata["postprocessing"]["exported"][1]["file"]
@async run(`"paraview.exe" $File`)

true


text:
@testset "when" begin
    @test 2 === @when let 1 = 1
        2
    end
    @test 2 === @when 1 = 1 2

    @test 1 === @when let a = 1
        a
    end

    @test 1 === @when a = 1 a

    @test nothing === @when let (a, b) = 1
        a + b
    end

    @test nothing === @when (a, b) = 1 a + b

    ab = (2, 3)
    @test 5 === @when let (a, b) = ab
        a + b
    end

    @test 5 === @when (a, b) = ab a + b

    @data WhenTest begin
        WhenTest_1(Int)
        WhenTest_2(Int)
    end

    var1 = WhenTest_1(2)
    var2 = WhenTest_2(2)
    @test 200 === @when let WhenTest_1(x) = var1,
                            @inline WhenAction(x) = 100x
        WhenAction(x)
    end

    @test 200 === @when WhenTest_1(x) = var1 begin
        100x
    end

    @test nothing === @when let WhenTest_1(x) = var2,
                                @inline WhenAction(x) = 100x
        WhenAction(x)
    end

end


text:
using DataDeps

register(DataDep(
  "spect-heart",
  "http://archive.ics.uci.edu/ml/datasets/SPECT+Heart",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/spect/SPECT.train",
  "0fd0258bfc8eb986623a45890676465ee92e738a95fdc01a6aef893d9c920115",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:23,
        categoric_indices=2:23,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"""
Base class for tuple-based fitness schemes.

Type parameters:
  * `N` is the number of the objectives
  * `F` is the type of each objective
  * `FA` is the actual type of the multi-objective fitness
  * `MIN` if objectives should be minimized or maximized
  * `AGG` the type of aggregator
"""
abstract type TupleFitnessScheme{N,F<:Number,FA,MIN,AGG} <: FitnessScheme{FA} end

@inline numobjectives(::TupleFitnessScheme{N}) where {N} = N
@inline fitness_eltype(::TupleFitnessScheme{N,F}) where {N,F} = F
@inline is_minimizing(::TupleFitnessScheme{N,F,FA,MIN}) where {N,F,FA,MIN} = MIN

@generated nafitness(::TupleFitnessScheme{N,F,NTuple{N,F}}) where {N,F} = ntuple(_ -> convert(F, NaN), N)
isnafitness(f::NTuple{N,F}, ::TupleFitnessScheme{N,F}) where {N,F} = any(isnan, f)

aggregate(f::NTuple{N,F}, fs::TupleFitnessScheme{N,F}) where {N,F} = fs.aggregator(f)

@inline is_better(f1::NTuple{N,F}, f2::NTuple{N,F}, fs::TupleFitnessScheme{N,F,NTuple{N,F}}) where {N,F} = hat_compare(f1, f2, fs, -1) == -1
@inline is_worse(f1::NTuple{N,F}, f2::NTuple{N,F}, fs::TupleFitnessScheme{N,F,NTuple{N,F}}) where {N,F} = hat_compare(f1, f2, fs, 1) == 1

"""
Pareto dominance for `N`-tuple (`N`≧1) fitnesses.

`aggregator::AGG` is a function mapping tuple fitness to a single numerical value.
Might be used for comparisons (or not, depending on the setup).
Always used when printing fitness vectors though.
"""
struct ParetoFitnessScheme{N,F<:Number,MIN,AGG} <: TupleFitnessScheme{N,F,NTuple{N,F},MIN,AGG}
    aggregator::AGG    # fitness aggregation function

    ParetoFitnessScheme{N,F,MIN,AGG}(; aggregator::AGG=sum) where {N, F<:Number, MIN, AGG} =
        new{N,F,MIN,AGG}(aggregator)

    ParetoFitnessScheme{N,F,MIN}(; aggregator::AGG=sum) where {N, F<:Number, MIN, AGG} =
        new{N,F,MIN,AGG}(aggregator)

    ParetoFitnessScheme{N,F}(; is_minimizing::Bool=true, aggregator::AGG=sum) where {N, F<:Number, AGG} =
        new{N,F,is_minimizing,AGG}(aggregator)

    ParetoFitnessScheme{N}(; fitness_type::Type{F} = Float64,
                            is_minimizing::Bool=true, aggregator::AGG=sum) where {N, F<:Number, AGG} =
        new{N,fitness_type,is_minimizing,AGG}(aggregator)
end

# comparison and for minimizing Pareto scheme
function hat_compare_pareto(u::NTuple{N,F}, v::NTuple{N,F}, expected::Int=0) where {N, F}
    res = 0
    @inbounds for i in 1:N
        delta = u[i] - v[i]
        if delta > 0.0
            if res == 0
                res = 1
                if expected == -1 return res end
            elseif res == -1
                return 0 # non-dominated
            end
        elseif delta < 0.0
            if res == 0
                res = -1
                if expected == 1 return res end
            elseif res == 1
                return 0 # non-dominated
            end
        end
    end
    return res
end

hat_compare(f1::NTuple{N,F}, f2::NTuple{N,F}, fs::ParetoFitnessScheme{N,F,true}, expected::Int=0) where {N,F} =
    hat_compare_pareto(f1, f2, expected)
hat_compare(f1::NTuple{N,F}, f2::NTuple{N,F}, fs::ParetoFitnessScheme{N,F,false}, expected::Int=0) where {N,F} =
    hat_compare_pareto(f2, f1, expected)

"""
ϵ-dominance for `N`-tuple (`N`≧1) fitnesses.

`aggregator::AGG` is a function mapping tuple fitness to a single numerical value.
Might be used for comparisons (or not, depending on the setup).
Always used when printing fitness vectors though.
"""
struct EpsDominanceFitnessScheme{N,F<:Number,MIN,AGG} <: FitnessScheme{NTuple{N,F}}
    ϵ::F              # ɛ-domination threshold
    aggregator::AGG    # fitness aggregation function

    function EpsDominanceFitnessScheme{N,F}(
            ϵ::F; is_minimizing::Bool=true, aggregator::AGG=sum) where {N, F<:Number, AGG}
        ϵ>0.0 || throw(ArgumentError("ϵ must be positive"))
        new{N,F,is_minimizing,AGG}(ϵ, aggregator)
    end
end

EpsDominanceFitnessScheme{N}(ϵ::F; fitness_type::Type{F} = Float64,
        is_minimizing::Bool=true, aggregator::AGG=sum) where {N, F<:Number, AGG} =
    EpsDominanceFitnessScheme{N,fitness_type}(ϵ; is_minimizing=is_minimizing, aggegator=aggregator)

# comparison for minimizing ϵ-dominance scheme
function hat_compare_ϵ(u::NTuple{N,F}, v::NTuple{N,F},
                       ϵ::F, expected::Int=0) where {N,F}
    res = 0 # true if any u[i] < v[i] + ϵ
    @inbounds for i in 1:N
        delta = u[i] - v[i] - ϵ
        if delta > 0.0
            if res == 0
                res = 1
                if expected == -1 return 1 end
            elseif res == -1
                return 0 # non-dominated
            end
        elseif delta < 0.0
            if res == 0
                res = -1
                if expected == 1 return -1 end
            elseif res == 1
                return 0 # non-dominated
            end
        end
    end
    return res
end

hat_compare(f1::NTuple{N,F}, f2::NTuple{N,F}, fs::EpsDominanceFitnessScheme{N,F,true}, expected::Int=0) where {N,F} =
    hat_compare_ϵ(f1, f2, fs.ϵ, expected)
hat_compare(f1::NTuple{N,F}, f2::NTuple{N,F}, fs::EpsDominanceFitnessScheme{N,F,false}, expected::Int=0) where {N,F} =
    hat_compare_ϵ(f2, f1, fs.ϵ, expected)

# ϵ-index of the fitness component for minimizing scheme
@inline function ϵ_index(u::F, ϵ::F, ::Type{Val{true}}) where F
    if isnan(u)
        return (typemax(Int), zero(F))
    else
        u_div_ϵ = clamp(u/ϵ, convert(F, typemin(Int)), convert(F, typemax(Int)))
        ix = floor(Int, u_div_ϵ+10eps(F))
        return (ix, max(zero(F), u_div_ϵ-ix))
    end
end

# ϵ-index of the fitness component for maximizing scheme
@inline function ϵ_index(u::F, ϵ::F, ::Type{Val{false}}) where F
    if isnan(u)
        return (typemin(Int), zero(F))
    else
        u_div_ϵ = clamp(u/ϵ, convert(F, typemin(Int)), convert(F, typemax(Int)))
        ix = ceil(Int, u_div_ϵ)
        return (ix, ix-u_div_ϵ)
    end
end

# vectorized ϵ-index
@generated function ϵ_index(u::NTuple{N,F}, ϵ::Vector{F}, is_minimizing::Type{Val{MIN}}) where {N,F,MIN}
    quote
        pairs = Base.Cartesian.@ntuple $N i -> ϵ_index(u[i], ϵ[i], is_minimizing)
        ix = Base.Cartesian.@ntuple $N i -> pairs[i][1]
        sqrdist = zero(F)
        Base.Cartesian.@nexprs $N i -> sqrdist += pairs[i][2]^2
        return ix, sqrt(sqrdist)
    end
end

"""
ϵ-box indexed representation of the N-tuple fitness.

Used together with `EpsBoxDominanceFitnessScheme`.
"""
struct IndexedTupleFitness{N,F}
    orig::NTuple{N,F}       # original fitness
    agg::F                  # aggregated fitness
    index::NTuple{N,Int}    # ϵ-index vector
    dist::F                 # distance between ϵ-index vector and the original fitness

    function IndexedTupleFitness(u::NTuple{N,F}, agg::F, ϵ::Vector{F}, is_minimizing::Type{Val{MIN}}) where {N, F, MIN}
        ix, dist = ϵ_index(u, ϵ, is_minimizing)
        return new{N,F}(u, agg, ix, dist)
    end
end

IndexedTupleFitness(u::NTuple{N,F}, agg::F, ϵ::F, is_minimizing::Type{Val{MIN}}) where {N, F, MIN} =
    IndexedTupleFitness(u, agg, fill(ϵ, N), is_minimizing)

Base.convert(::Type{NTuple{N,F}}, fitness::IndexedTupleFitness{N,F}) where {N, F} = fitness.orig

@generated function nafitness(::Type{IndexedTupleFitness{N,F}}) where {N, F}
    quote
        IndexedTupleFitness(Base.Cartesian.@ntuple($N, _ -> convert($F, NaN)),
                            NaN, 1.0, Val{true})
    end
end

# comparison for minimizing ϵ-box dominance scheme
"""
Returns a tuple of `u` and `v` comparison:
  * `-1`: u≺v
  * `0`: u and v non-dominating
  * `1`: u≻v
and whether `u` index fully matches `v` index.
"""
function hat_compare_ϵ_box(
        u::IndexedTupleFitness{N,F},
        v::IndexedTupleFitness{N,F},
        is_minimizing::Bool=true, expected::Int=0) where {N,F}
    comp = 0
    @inbounds for (ui, vi) in zip(u.index, v.index)
        if ui > vi
            if comp == 0
                comp = 1
                if expected < 0 return (1, false) end
            elseif comp == -1
                return (0, false)  # non-dominated
            end
        elseif ui < vi
            if comp == 0
                comp = -1
                if expected > 0 return (-1, false) end
            elseif comp == 1
                return (0, false) # non-dominated
            end
        end
    end
    if !is_minimizing
        comp = -comp
    end
    if comp != 0
        return (comp, false)
    else
        uv_diff = u.dist - v.dist
        return (uv_diff < -10eps(F) ? -1 :
                uv_diff > 10eps(F) ? 1 : 0, true)
    end
end

function check_epsbox_ϵ(ϵ::Number, n::Int)
    ϵ>0.0 || throw(ArgumentError("ϵ must be positive"))
    return fill(ϵ, n)
end

function check_epsbox_ϵ(ϵ::AbstractVector{<:Number}, n::Int)
    length(ϵ)==n || throw(ArgumentError("The length of ϵ vector ($(length(ϵ))) does not match the specified fitness dimensions ($n)"))
    all(isposdef, ϵ) || throw(ArgumentError("ϵ must be positive"))
    return ϵ
end

"""
`EpsBoxDominanceFitnessScheme` defines ϵ-box dominance for
`N`-tuple (`N`≧1) fitnesses.
It operates with fitnesses of type `IndexedTupleFitness`.

`aggregator::AGG` is a function mapping tuple fitness to a single numerical value.
Might be used for comparisons (or not, depending on the setup).
Always used when printing fitness vectors though.
"""
struct EpsBoxDominanceFitnessScheme{N,F<:Number,MIN,AGG} <: TupleFitnessScheme{N,F,IndexedTupleFitness{N,F},MIN,AGG}
    ϵ::Vector{F}        # per-objective ɛ-domination thresholds
    aggregator::AGG     # fitness aggregation function

    EpsBoxDominanceFitnessScheme{N,F}(ϵ::Union{F,Vector{F}};
            is_minimizing::Bool=true, aggregator::AGG=sum) where {N,F<:Number,AGG} =
        new{N,F,is_minimizing,AGG}(check_epsbox_ϵ(ϵ, N), aggregator)

    EpsBoxDominanceFitnessScheme{N}(ϵ::Union{F,Vector{F}};
            is_minimizing::Bool=true, aggregator::AGG=sum) where {N,F<:Number,AGG} =
        new{N,F,is_minimizing,AGG}(check_epsbox_ϵ(ϵ, N), aggregator)
end

isnafitness(f::IndexedTupleFitness{N,F},
            fit_scheme::EpsBoxDominanceFitnessScheme{N,F}) where {N,F<:Number} =
    isnafitness(f.orig, fit_scheme)

EpsBoxDominanceFitnessScheme(fs::ParetoFitnessScheme{N,F}, ϵ::F=one(F)) where {N,F} =
    EpsBoxDominanceFitnessScheme{N,F}(ϵ, is_minimizing=is_minimizing(fs), aggregator=fs.aggregator)

EpsBoxDominanceFitnessScheme(fs::ParetoFitnessScheme{N,F}, ϵ::Vector{F}) where {N,F} =
    EpsBoxDominanceFitnessScheme{N,F}(ϵ, is_minimizing=is_minimizing(fs), aggregator=fs.aggregator)

EpsBoxDominanceFitnessScheme(fs::EpsDominanceFitnessScheme{N,F}, ϵ::Union{F,Vector{F}}=fs.ϵ) where {N,F} =
    EpsBoxDominanceFitnessScheme{N,F}(ϵ, is_minimizing=is_minimizing(fs), aggregator=fs.aggregator)

ParetoFitnessScheme(fs::EpsBoxDominanceFitnessScheme{N,F}) where {N,F} =
  ParetoFitnessScheme{N,F}(is_minimizing=is_minimizing(fs), aggregator=fs.aggregator)

Base.convert(::Type{IndexedTupleFitness{N,F}}, fitness::NTuple{N,F},
             fs::EpsBoxDominanceFitnessScheme{N,F}) where {N,F} =
    IndexedTupleFitness(fitness, aggregate(fitness, fs), fs.ϵ, Val{is_minimizing(fs)})

Base.convert(::Type{IndexedTupleFitness}, fitness::NTuple{N,F},
             fs::EpsBoxDominanceFitnessScheme{N,F}) where {N,F} =
    IndexedTupleFitness(fitness, aggregate(fitness, fs), fs.ϵ, Val{is_minimizing(fs)})

Base.convert(::Type{NTuple{N,F}}, fitness::IndexedTupleFitness{N,F},
             fs::EpsBoxDominanceFitnessScheme{N,F}) where {N,F} = fitness.orig

hat_compare(u::IndexedTupleFitness{N,F}, v::IndexedTupleFitness{N,F},
            fs::EpsBoxDominanceFitnessScheme{N,F}, expected::Int=0) where {N,F} =
    hat_compare_ϵ_box(u, v, is_minimizing(fs), expected)

hat_compare(u::NTuple{N,F}, v::IndexedTupleFitness{N,F},
            fs::EpsBoxDominanceFitnessScheme{N,F}, expected::Int=0) where {N,F} =
    hat_compare(convert(IndexedTupleFitness{N,F}, u, fs), v, fs, expected)

hat_compare(u::IndexedTupleFitness{N,F}, v::NTuple{N,F},
            fs::EpsBoxDominanceFitnessScheme{N,F}, expected::Int=0) where {N,F} =
    hat_compare(u, convert(IndexedTupleFitness{N,F}, v, fs), fs, expected)

hat_compare(u::NTuple{N,F}, v::NTuple{N,F},
            fs::EpsBoxDominanceFitnessScheme{N,F}, expected::Int=0) where {N,F} =
    hat_compare(convert(IndexedTupleFitness{N,F}, u, fs),
                convert(IndexedTupleFitness{N,F}, v, fs), fs, expected)

# special overload that strips index equality flag
(hc::HatCompare{FS})(u::IndexedTupleFitness{N,F},
                     v::IndexedTupleFitness{N,F}) where {FS<:EpsBoxDominanceFitnessScheme,N,F} =
    hat_compare(u, v, hc.fs)[1]


text:
export DxGeo, DyGeo, DxExp,DξExp, DxLog, DyLog
export DforwardLogs
@doc doc"""
    DxGeo(M,x,y,t,η)
computes $D_xg(t;x,y)[\eta]$.

# See also
 [`DyGeo`](@ref), [`jacobiField`](@ref)
"""
DxGeo(M::mT,x::P,y::P,t::Number,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = jacobiField(M,x,y,t,η,βDgx)
@doc doc"""
    DyGeo(M,x,y,t,η)
computes $D_yg(t;x,y)[\eta]$.

# See also
 [`DxGeo`](@ref), [`jacobiField`](@ref)
"""
DyGeo(M::mT,x::P,y::P,t::Number,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = jacobiField(M,y,x,1-t,η,βDgx)
@doc doc"""
    DxExp(M,x,ξ,η)
computes $D_x\exp_x\xi[\eta]$.

# See also
 [`DξExp`](@ref), [`jacobiField`](@ref)
"""
DxExp(M::mT,x::P,ξ::T,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector}= jacobiField(M,x,exp(M,x,ξ),1.0,η,βDexpx)
@doc doc"""
    DξExp(M,x,ξ,η)
computes $D_\xi\exp_x\xi[\eta]$.
Note that $\xi\in T_\xi(T_x\mathcal M) = T_x\mathcal M$ is still a tangent vector.

# See also
 [`DxExp`](@ref), [`jacobiField`](@ref)
"""
DξExp(M::mT,x::P,ξ::T,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = jacobiField(M,x,exp(M,x,ξ),1.0,η,βDexpξ)
@doc doc"""
    DxLog(M,x,y,η)
computes $D_xlog_xy[\eta]$.

# See also
 [`DyLog`](@ref), [`jacobiField`](@ref)
"""
DxLog(M::mT,x::P,y::P,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = jacobiField(M,x,y,0.0,η,βDlogx)
@doc doc"""
    DyLog(M,x,y,η)
computes $D_ylog_xy[\eta]$.

# See also
 [`DxLog`](@ref), [`jacobiField`](@ref)
"""
DyLog(M::mT,x::P,y::P,η::T) where {mT <: Manifold, P <: MPoint, T<: TVector} = jacobiField(M,y,x,1.0,η,βDlogy)

@doc doc"""
    ν = DforwardLogs(M,x,ξ)

compute the differenital of [`forwardLogs`](@ref) $F$ on the [`Power`](@ref) manifold
`M` at `x` and direction `ξ` ,
in the power manifold array, the differential of the function

```math
F_i(x) = \sum_{j\in\mathcal I_i} \log_{x_i} x_j$, \quad i \in \mathcal G,
```

where $\mathcal G$ is the set of indices of the [`Power`](@ref) manifold `M`
and $\mathcal I_i$ denotes the forward neighbors of $i$.

# Input
* `M`     – a [`Power`](@ref) manifold
* `x`     – a [`PowPoint`](@ref).
* `ξ`     – a [`PowTVector`](@ref).

# Ouput
* `ν` – resulting tangent vector in $T_x\mathcal N$ representing the differentials of the logs, where
  $\mathcal N$ is thw power manifold with the number of dimensions added to `size(x)`.
"""
function DforwardLogs(M::Power,x::PowPoint,ξ::PowTVector)::PowTVector
  sξ = size(ξ)
  R = CartesianIndices(sξ)
  d = length(sξ)
  maxInd = [last(R).I...]
  d2 = (d>1) ? ones(Int,d+1) + (d-1)*(1:(d+1) .== d+1 ) : d
  N = Power(M.manifold,(sξ...,d))
  ν = zeroTVector(N, repeat(x,inner=d2) )
  for i in R # iterate over all pixel
    for k in 1:d # for all direction combinations
      I = [i.I...] # array of index
      J = I .+ 1 .* (1:d .== k) #i + e_k is j
      if all( J .<= maxInd )
        # this is neighbor in range,
        j = CartesianIndex{d}(J...) # neigbbor index as Cartesian Index
        # collects two, namely in kth direction since xi appears as base and arg
        ν[i,k] = DxLog(M.manifold,x[i],x[j],ξ[i]) + DyLog(M.manifold,x[i],x[j],ξ[j])
      end
    end # directions
  end # i in R
  return ν
end


text:
using Statistics

function fill_refs!(refs::AbstractArray, X::AbstractArray,
                    breaks::AbstractVector, extend::Union{Bool, Missing})
    n = length(breaks)
    lower = first(breaks)
    upper = last(breaks)

    @inbounds for i in eachindex(X)
        x = X[i]

        if ismissing(x)
            refs[i] = 0
        elseif extend === true && x == upper
            refs[i] = n-1
        elseif extend !== true && !(lower <= x < upper)
            extend === missing ||
                throw(ArgumentError("value $x (at index $i) does not fall inside the breaks: " *
                                    "adapt them manually, or pass extend=true or extend=missing"))
            refs[i] = 0
        else
            refs[i] = searchsortedlast(breaks, x)
        end
    end
end

"""
    default_formatter(from, to, i; leftclosed, rightclosed)

Provide the default label format for the `cut(x, breaks)` method.
"""
default_formatter(from, to, i; leftclosed, rightclosed) =
    string(leftclosed ? "[" : "(", from, ", ", to, rightclosed ? "]" : ")")

@doc raw"""
    cut(x::AbstractArray, breaks::AbstractVector;
        labels::Union{AbstractVector{<:AbstractString},Function},
        extend::Union{Bool,Missing}=false, allowempty::Bool=false)

Cut a numeric array into intervals at values `breaks`
and return an ordered `CategoricalArray` indicating
the interval into which each entry falls. Intervals are of the form `[lower, upper)`,
i.e. the lower bound is included and the upper bound is excluded, except for the last
interval which is closed on both ends, i.e. `[lower, upper]`.

If `x` accepts missing values (i.e. `eltype(x) >: Missing`) the returned array will
also accept them.

# Keyword arguments
* `extend::Union{Bool, Missing}=false`: when `false`, an error is raised if some values
  in `x` fall outside of the breaks; when `true`, breaks are automatically added to include
  all values in `x`, and the upper bound is included in the last interval; when `missing`,
  values outside of the breaks generate `missing` entries.
* `labels::Union{AbstractVector,Function}`: a vector of strings giving the names to use for
  the intervals; or a function `f(from, to, i; leftclosed, rightclosed)` that generates
  the labels from the left and right interval boundaries and the group index. Defaults to
  `"[from, to)"` (or `"[from, to]"` for the rightmost interval if `extend == true`).
* `allowempty::Bool=false`: when `false`, an error is raised if some breaks appear
  multiple times, generating empty intervals; when `true`, duplicate breaks are allowed
  and the intervals they generate are kept as unused levels
  (but duplicate labels are not allowed).

# Examples
```jldoctest
julia> using CategoricalArrays

julia> cut(-1:0.5:1, [0, 1], extend=true)
5-element CategoricalArray{String,1,UInt32}:
 "[-1.0, 0.0)"
 "[-1.0, 0.0)"
 "[0.0, 1.0]"
 "[0.0, 1.0]"
 "[0.0, 1.0]"

julia> cut(-1:0.5:1, 2)
5-element CategoricalArray{String,1,UInt32}:
 "Q1: [-1.0, 0.0)"
 "Q1: [-1.0, 0.0)"
 "Q2: [0.0, 1.0]"
 "Q2: [0.0, 1.0]"
 "Q2: [0.0, 1.0]"

julia> cut(-1:0.5:1, 2, labels=["A", "B"])
5-element CategoricalArray{String,1,UInt32}:
 "A"
 "A"
 "B"
 "B"
 "B"

julia> fmt(from, to, i; leftclosed, rightclosed) = "grp $i ($from//$to)"
fmt (generic function with 1 method)

julia> cut(-1:0.5:1, 3, labels=fmt)
5-element CategoricalArray{String,1,UInt32}:
 "grp 1 (-1.0//-0.333333)"
 "grp 1 (-1.0//-0.333333)"
 "grp 2 (-0.333333//0.333333)"
 "grp 3 (0.333333//1.0)"
 "grp 3 (0.333333//1.0)"
```
"""
@inline function cut(x::AbstractArray, breaks::AbstractVector;
                    extend::Union{Bool, Missing}=false,
                    labels::Union{AbstractVector{<:AbstractString},Function}=default_formatter,
                    allowmissing::Union{Bool, Nothing}=nothing,
                    allow_missing::Union{Bool, Nothing}=nothing,
                    allowempty::Bool=false)
    if allow_missing !== nothing
        Base.depwarn("allow_missing argument is deprecated, use extend=missing instead",
                     :cut)
        extend = missing
    end
    if allowmissing !== nothing
        Base.depwarn("allowmissing argument is deprecated, use extend=missing instead",
                     :cut)
         extend = missing
    end
    return _cut(x, breaks, extend, labels, allowempty)
end

# Separate function for inferability (thanks to inlining of cut)
function _cut(x::AbstractArray{T, N}, breaks::AbstractVector,
              extend::Union{Bool, Missing},
              labels::Union{AbstractVector{<:AbstractString},Function},
              allowempty::Bool=false) where {T, N}
    if !allowempty && !allunique(breaks)
        throw(ArgumentError("all breaks must be unique unless `allowempty=true`"))
    end

    if !issorted(breaks)
        breaks = sort(breaks)
    end

    if extend === true
        xnm = T >: Missing ? skipmissing(x) : x
        length(breaks) >= 1 || throw(ArgumentError("at least one break must be provided"))
        local min_x, max_x
        try
            min_x, max_x = extrema(xnm)
        catch err
            if T >: Missing && all(ismissing, xnm)
                if length(breaks) < 2
                    throw(ArgumentError("could not extend breaks as all values are missing: " *
                                        "please specify at least two breaks manually"))
                else
                    min_x, max_x = missing, missing
                end
            else
                rethrow(err)
            end
        end
        if !ismissing(min_x) && breaks[1] > min_x
            breaks = [min_x; breaks]
        end
        if !ismissing(max_x) && breaks[end] < max_x
            breaks = [breaks; max_x]
        end
        length(breaks) > 1 ||
            throw(ArgumentError("could not extend breaks as all values are equal: " *
                                "please specify at least two breaks manually"))
    end

    refs = Array{DefaultRefType, N}(undef, size(x))
    try
        fill_refs!(refs, x, breaks, extend)
    catch err
        # So that the error appears to come from cut() itself,
        # since it refers to its keyword arguments
        if isa(err, ArgumentError)
            throw(err)
        else
            rethrow(err)
        end
    end

    n = length(breaks)
    n >= 2 || throw(ArgumentError("at least two breaks must be provided when extend is not true"))
    if labels isa Function
        from = breaks[1:n-1]
        to = breaks[2:n]
        levs = Vector{String}(undef, n-1)
        for i in 1:n-2
            levs[i] = labels(from[i], to[i], i,
                             leftclosed=breaks[i] != breaks[i+1], rightclosed=false)
        end
        levs[end] = labels(from[end], to[end], n-1,
                           leftclosed=breaks[end-1] != breaks[end],
                           rightclosed=coalesce(extend, false))
    else
        length(labels) == n-1 || throw(ArgumentError("labels must be of length $(n-1), but got length $(length(labels))"))
        # Levels must have element type String for type stability of the result
        levs::Vector{String} = copy(labels)
    end
    if !allunique(levs)
        if labels === default_formatter
            throw(ArgumentError("all labels must be unique, but `breaks` contains duplicates: " *
                                "specify custom `labels` with unique names"))
        else
            throw(ArgumentError("all labels must be unique"))
        end
    end

    pool = CategoricalPool(levs, true)
    S = T >: Missing || extend isa Missing ? Union{String, Missing} : String
    CategoricalArray{S, N}(refs, pool)
end

"""
    quantile_formatter(from, to, i; leftclosed, rightclosed)

Provide the default label format for the `cut(x, ngroups)` method.
"""
quantile_formatter(from, to, i; leftclosed, rightclosed) =
    string("Q", i, ": ", leftclosed ? "[" : "(", from, ", ", to, rightclosed ? "]" : ")")

"""
    cut(x::AbstractArray, ngroups::Integer;
        labels::Union{AbstractVector{<:AbstractString},Function},
        allowempty::Bool=false)

Cut a numeric array into `ngroups` quantiles, determined using `quantile`.

If `x` contains `missing` values, they are automatically skipped when computing
quantiles.

# Keyword arguments
* `labels::Union{AbstractVector,Function}`: a vector of strings giving the names to use for
  the intervals; or a function `f(from, to, i; leftclosed, rightclosed)` that generates
  the labels from the left and right interval boundaries and the group index. Defaults to
  `"Qi: [from, to)"` (or `"Qi: [from, to]"` for the rightmost interval if `extend == true`).
* `allowempty::Bool=false`: when `false`, an error is raised if some quantiles breakpoints
  are equal, generating empty intervals; when `true`, duplicate breaks are allowed
  and the intervals they generate are kept as unused levels
  (but duplicate labels are not allowed).
"""
function cut(x::AbstractArray, ngroups::Integer;
             labels::Union{AbstractVector{<:AbstractString},Function}=quantile_formatter,
             allowempty::Bool=false)
    xnm = eltype(x) >: Missing ? skipmissing(x) : x
    breaks = Statistics.quantile(xnm, (1:ngroups-1)/ngroups)
    if !allowempty && !allunique(breaks)
        n = length(unique(breaks)) - 1
        throw(ArgumentError("cannot compute $ngroups quantiles: `quantile` " *
                            "returned only $n groups due to duplicated values in `x`." *
                            "Pass `allowempty=true` to allow empty quantiles or " *
                            "choose a lower value for `ngroups`."))
    end
    cut(x, breaks; extend=true, labels=labels, allowempty=allowempty)
end


text:
## ----------------------------------------------------------------------------
# beta vs stuff
let
    method = ME_Z_EXPECTED_G_MOVING
    cGLC_plt = plot(;xlabel = "cGLC", ylabel = "beta")
    D_plt = plot(;xlabel = "D", ylabel = "beta")
    for exp in EXPS
        datfile = ME_INDEX[method, :DFILE, exp]
        dat = deserialize(datfile)
        beta = dat[:exp_beta]

        params = (;label = "", color = exp_colors[exp],
            alpha = 0.7, ms = 7
        )
        cGLC = Hd.val("cGLC", exp)
        D = Hd.val("D", exp)
        scatter!(cGLC_plt, [cGLC], [beta]; params...)
        scatter!(D_plt, [D], [beta]; params...)
    end
    mysavefig([cGLC_plt, D_plt], "beta_vs_stuff"; method)
end

## ----------------------------------------------------------------------------
# flux vs beta
let
    method = ME_Z_EXPECTED_G_MOVING
    p = plot(title = iJR.PROJ_IDER, xlabel = "beta", ylabel = "biom")
    for exp in EXPS
        datfile = ME_INDEX[method, :DFILE, exp]
        dat = deserialize(datfile)
        model = dat[:model]
        objidx = ChU.rxnindex(model, iJR.BIOMASS_IDER)
        epouts = dat[:epouts]
        exp_beta = dat[:exp_beta]
        exp_xi = Hd.val("xi", exp)
        scatter!(p, [exp_beta], [Hd.val("D", exp)], ms = 12, color = :white, label = "")

        betas = collect(keys(epouts)) |> sort
        bioms = [ChU.av(model, epouts[beta], objidx) for beta in betas]
        scatter!(p, betas, bioms, label = "", color = :black, alpha = 0.2)

    end
    mysavefig(p, "obj_val_vs_beta"; method)
end


text:
using zad2, Test

@testset "testy " begin

	funk,xL,xU = x->x^2,-400,400
	@test zad22(funk,xL,xU) == (0,0)

	funk,xL,xU = x->sin(x),0.1,pi
	@test zad22(funk,xL,xU) == (0.0,3.14)

	funk,xL,xU = x->x^6+x^4+x^5+x^2+x ,-2, 2
	@test zad22(funk,xL,xU) == (-0.21,-0.34)

end


text:
isdefined(Base, :__precompile__) && __precompile__()

module OptimalTransmissionRouting

# import Compat
import JuMP
import Memento
import Images
import FileIO
import ColorTypes
import Plots

# Create our module level logger (this will get precompiled)
const _LOGGER = Memento.getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via `getlogger(PowerModels)`
# NOTE: If this line is not included then the precompiled `_PM._LOGGER` won't be registered at runtime.
__init__() = Memento.register(_LOGGER)

include("core/define_weights.jl")
include("core/image2weight.jl")
include("core/optimal_routing.jl")
include("core/costs_and_equipment_details.jl")
include("core/a_star.jl")
include("core/impedance_parameters.jl")
include("core/spatial_data_preparation.jl")
include("io/plot_optimal_path.jl")

# Spatial image files
# include("spatial_image_files/clc_agricultural.tif")
# include("spatial_image_files/clc_mountains.tif")
# include("spatial_image_files/clc_natural.tif")
# include("spatial_image_files/clc_urban.tif")
# include("spatial_image_files/clc_graph_Europe.tif")
# include("spatial_image_files/clc_grid_225kv_400kv.tif")
# include("spatial_image_files/natura2000.tif")
# include("spatial_image_files/clc_railroads.tif")
# include("spatial_image_files/clc_roads.tif")
# include("spatial_image_files/sea.tif")

end


text:
@testset "Triangle" begin
    # Make some points between (0,0,0) and (1,1,1) and put in a point cloud
    N = 1000
    pos = rand(SVector{3,Float32}, N)
    triangle = Triangle(SVector(0.0, 0.0), SVector(1.0, 0.0), SVector(0.0, 1.0))

    search_correct = true
    for i = 1:length(pos)
        v = pos[i]
        expected = v[1] + v[2] < 1.0
        actual = v ∈ triangle
        search_correct &= actual == expected
    end
    @test search_correct

    triangle = Triangle(SVector(0.5, 0.5), SVector(0.5, 1.5), SVector(1.5, 1.0))
    @test SVector(0.0, 0.0, 0.0) ∉ triangle
    @test SVector(0.0, 0.0, 1.0) ∉ triangle
    @test SVector(0.0, 0.0, -1.0) ∉ triangle
    @test SVector(1.0, 1.0, 0.0) ∈ triangle
    @test SVector(1.0, 1.0, -1.0) ∈ triangle
    @test SVector(1.0, 1.0, 1.0) ∈ triangle

    @test boundingbox(triangle) == BoundingBox(0.5, 0.5, -Inf, 1.5, 1.5, Inf)

    @test area(triangle) == 0.5
end


text:
const JULIAN_EPOCH = TTEpoch(AstroDates.JULIAN_EPOCH, AstroDates.H12)
const J2000_EPOCH = TTEpoch(AstroDates.J2000_EPOCH, AstroDates.H12)
const MODIFIED_JULIAN_EPOCH = TTEpoch(AstroDates.MODIFIED_JULIAN_EPOCH, AstroDates.H00)
const FIFTIES_EPOCH = TTEpoch(AstroDates.FIFTIES_EPOCH, AstroDates.H00)
const CCSDS_EPOCH = TTEpoch(AstroDates.CCSDS_EPOCH, AstroDates.H00)
const GALILEO_EPOCH = TTEpoch(AstroDates.GALILEO_EPOCH, AstroDates.H00)
const GPS_EPOCH = TTEpoch(AstroDates.GPS_EPOCH, AstroDates.H00)
const UNIX_EPOCH = TAIEpoch(AstroDates.UNIX_EPOCH, Time(0, 0, 10.0))

const PAST_INFINITY = TAIEpoch(UNIX_EPOCH, -Inf)
const FUTURE_INFINITY = TAIEpoch(UNIX_EPOCH, Inf)

const EPOCH_77 = TAIEpoch(1977, 1, 1)

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
using PGFPlotsX, BasisFunctions

f = exp

F1 = Expansion( ChebyshevT(10), rand(ChebyshevT(10)))
F2 = Expansion( Fourier(10), rand(Fourier(10)))

for F in (F1,F2)
    Plot(F)

    Plot(F;n=300)

    Axis(F,f)
    Axis(F,F)
    Axis(f,F)
end

Plot(F2;plot_extension=true)
Axis(Plot(F2;plot_complex=true)...)


text:
module FemtoCleaner

# For interactive development
using Revise

using Base.Distributed
using GitHub
using GitHub: GitHubAPI, GitHubWebAPI, Checks
using HTTP
using Deprecations
using CSTParser
using Deprecations: isexpr
using MbedTLS
using JSON
using AbstractTrees: children
using Base: LibGit2

include("workqueue.jl")

function with_cloned_repo(f, api::GitHubWebAPI, repo, auth)
    creds = LibGit2.UserPasswordCredentials(String(copy(Vector{UInt8}("x-access-token"))), String(copy(Vector{UInt8}(auth.token))))
    repo_url = "https://github.com/$(get(repo.full_name))"
    local_dir = mktempdir()
    try
        enabled = gc_enable(false)
        lrepo = LibGit2.clone(repo_url, local_dir; payload=Nullable(creds))
        gc_enable(enabled)
        f((lrepo, local_dir))
    finally
        rm(local_dir, force=true, recursive=true)
    end
end

function with_pr_branch(f, api, repo, auth)
    with_cloned_repo(api, repo, auth) do x
        LibGit2.branch!(lrepo, "fbot/deps", track=LibGit2.Consts.REMOTE_ORIGIN)
        f(x)
    end
end

if VERSION < v"0.7.0-DEV.695"
    include("blame.jl")
else
    using LibGit2: GitBlame
end

function deprecations_for_repo(lrepo, local_dir, is_julia_itself)
    if is_julia_itself
        ver = readstring(joinpath(local_dir, "VERSION"))
        hunk = GitBlame(lrepo, "VERSION")[1]
        l, r = LibGit2.revcount(lrepo, string(hunk.orig_commit_id), "HEAD")
        vers = Pkg.Reqs.parse(IOBuffer("julia $ver+$(l+r)"))
    else
        vers = Pkg.Reqs.parse(joinpath(local_dir, "REQUIRE"))
    end
    deps = Deprecations.applicable_deprecations(vers)

end

function process_deprecations(lrepo, local_dir; is_julia_itself=false, deps = deprecations_for_repo(lrepo, local_dir, is_julia_itself))
    changed_any = false
    problematic_files = String[]
    all_files = String[]
    all_files = String[]
    for (root, dirs, files) in walkdir(local_dir)
        for file in files
            fpath = joinpath(root, file)
            (endswith(fpath, ".jl") || endswith(fpath, ".md")) || continue
            file == "NEWS.md" && continue
            push!(all_files, fpath)
        end
    end
    max_iterations = 30
    iteration_counter = fill(0, length(all_files))
    # Iterate. Some rewrites may expose others
    while any(x->x != -1, iteration_counter)
        # We need to redo the analysis after every fetmocleaning round, since
        # things may have changes as the result of an applied rewrite.
        analysis = Deprecations.process_all(filter(f->endswith(f, ".jl"), all_files))
        for (i, fpath) in enumerate(all_files)
            iteration_counter[i] == -1 && continue
            problematic_file = false
            file_analysis = endswith(fpath, ".jl") ? (analysis[1], analysis[2][fpath]) : nothing
            try
                if !Deprecations.edit_file(fpath, deps, endswith(fpath, ".jl") ? edit_text : edit_markdown;
                        analysis = file_analysis)
                    # Nothing to change
                    iteration_counter[i] = -1
                elseif iteration_counter[i] > max_iterations
                    warn("Iterations did not converge for file $fpath")
                    problematic_file = true
                else
                    iteration_counter[i] += 1
                    changed_any = true
                end
            catch e
                warn("Exception thrown when fixing file $fpath. Exception was:\n",
                     sprint(showerror, e, catch_backtrace()))
                problematic_file = true
            end
            if problematic_file
                push!(problematic_files, file)
                iteration_counter[i] = -1
            end
            changed_any && !problematic_file && LibGit2.add!(lrepo, relpath(fpath, local_dir))
        end
    end
    changed_any, problematic_files
end

function push_repo(api::GitHubWebAPI, repo, auth; force=true, remote_branch="fbot/deps")
    creds = LibGit2.UserPasswordCredentials(String(copy(Vector{UInt8}("x-access-token"))), String(copy(Vector{UInt8}(auth.token))))
    enabled = gc_enable(false)
    LibGit2.push(repo, refspecs = ["+HEAD:refs/heads/$remote_branch"], force=force,
        payload=Nullable(creds))
    gc_enable(enabled)
end

struct SourceFile
    data::Vector{UInt8}
    offsets::Vector{UInt64}
end
Base.length(file::SourceFile) = length(file.offsets)

function SourceFile(data)
    offsets = UInt64[0]
    buf = IOBuffer(data)
    local line = ""
    while !eof(buf)
        line = readuntil(buf,'\n')
        !eof(buf) && push!(offsets, position(buf))
    end
    if !isempty(line) && line[end] == '\n'
        push!(offsets, position(buf)+1)
    end
    SourceFile(data,offsets)
end

function compute_line(file::SourceFile, offset)
    ind = searchsortedfirst(file.offsets, offset)
    ind <= length(file.offsets) && file.offsets[ind] == offset ? ind : ind - 1
end

function Base.getindex(file::SourceFile, line::Int)
    if line == length(file.offsets)
        return file.data[(file.offsets[end]+1):end]
    else
        # - 1 to skip the '\n'
        return file.data[(file.offsets[line]+1):max(1, file.offsets[line+1]-1)]
    end
end
Base.getindex(file::SourceFile, arr::AbstractArray) = [file[x] for x in arr]

function repl_to_annotation(fpath, file, lrepo, local_dir, repo, repl)
    # Compute blob hash for fpath
    blob_hash = LibGit2.addblob!(lrepo, joinpath(local_dir, fpath))
    # Put together description
    start_line = compute_line(file, first(repl.range))
    message = """
        $(repl.dep === nothing ? "" : repl.dep.description)
        In file $fpath starting at line $(start_line):
            $(strip(String(file[start_line])))
    """
    Checks.Annotation(
        basename(fpath),
        "https://github.com/$(GitHub.name(repo))/blob/$(blob_hash)/$(fpath)",
        compute_line(file, first(repl.range)), compute_line(file, last(repl.range)),
        "notice",
        message,
        string(typeof(repl.dep).name.name)[1:min(end, 40)],
        ""
    )
end

function collect_deprecation_annotations(api::GitHubAPI, lrepo, local_dir, repo, auth; is_julia_itself=false)
    deps = deprecations_for_repo(lrepo, local_dir, is_julia_itself)
    annotations = Checks.Annotation[]
    problematic_files = String[]
    for (root, dirs, files) in walkdir(local_dir)
        for file in files
            fpath = joinpath(root, file)
            (endswith(fpath, ".jl") || endswith(fpath, ".md")) || continue
            file == "NEWS.md" && continue
            contents = readstring(fpath)
            sfile = SourceFile(contents)
            problematic_file = false
            try
                if endswith(fpath, ".md")
                    changed_any, _ = Deprecations.edit_markdown(contents, deps)
                    if changed_any
                        blob_hash = LibGit2.addblob!(lrepo, joinpath(local_dir, fpath))
                        push!(annotations, Checks.Annotation(
                            basename(fpath),
                            "https://github.com/$(GitHub.name(repo))/blob/$(blob_hash)/$(fpath)",
                            1, 1,
                            "notice",
                            """
                            Code changes were found in this Markdown document:
                            $(fpath)
                            """,
                            "MarkdownCode",
                            ""
                        ))
                    end
                else
                    repls = Deprecations.text_replacements(contents, deps)
                    for repl in repls
                        push!(annotations, repl_to_annotation(fpath, sfile, lrepo, local_dir, repo, repl))
                    end
                end
            catch e
                warn("Exception thrown when fixing file $file. Exception was:\n",
                     sprint(showerror, e, catch_backtrace()))
                problematic_file = true
            end
            problematic_file && push!(problematic_files, file)
        end
    end
    annotations, problematic_files
end

function apply_deprecations(api::GitHubAPI, lrepo, local_dir, commit_sig, repo, auth; issue_number = 0)
    is_julia_itself = GitHub.name(repo) == "JuliaLang/julia"
    changed_any, problematic_files = process_deprecations(lrepo, local_dir; is_julia_itself=is_julia_itself)
    if changed_any
        LibGit2.commit(lrepo, "Fix deprecations"; author=commit_sig, committer=commit_sig, parent_ids=[LibGit2.GitHash(lrepo, "HEAD")])
        push_repo(api, lrepo, auth)
    end
    if issue_number != 0
        if changed_any
            create_pull_request(api, repo, auth=auth, params = Dict(
                    :issue => issue_number,
                    :base => get(repo.default_branch),
                    :head => "fbot/deps"
                )
            )
            if !isempty(problematic_files)
                create_comment(api, repo, issue_number, :pr, auth=auth, params = Dict(
                    :body => string("Failed to process the following files: ",
                                    join("`" .* problematic_files .* "`", ", "),
                                    ". :(")
                    )
                )
            end
            println("Created pull request for $(GitHub.name(repo))")
        else
            create_comment(api, repo, issue_number, :issue, params = Dict(
                :body => "No applicable deprecations were found in this repository."
            ), auth=auth)
            println("Processing complete for $(GitHub.name(repo)): no changes made")
        end
    else
        if changed_any
            body = "I fixed a number of deprecations for you"
            if !isempty(problematic_files)
                body *= string(", but I failed to process the following files: ",
                               join("`" .* problematic_files .* "`", ", "),
                               ". :(")
            end
            create_pull_request(api, repo, auth=auth, params = Dict(
                    :title => "Fix deprecations",
                    :body => body,
                    :base => get(repo.default_branch),
                    :head => "fbot/deps"
                )
            )
            println("Created pull request for $(GitHub.name(repo))")
        else
            println("Processing complete for $(GitHub.name(repo)): no changes made")
        end
    end
end

function my_diff_tree(repo::LibGit2.GitRepo, oldtree::LibGit2.GitTree, newtree::LibGit2.GitTree; pathspecs::AbstractString="")
    diff_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @LibGit2.check ccall((:git_diff_tree_to_tree, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{LibGit2.DiffOptionsStruct}),
                   diff_ptr_ptr, repo.ptr, oldtree.ptr, newtree.ptr, isempty(pathspecs) ? C_NULL : pathspecs)
    return LibGit2.GitDiff(repo, diff_ptr_ptr[])
end

function apply_deprecations_if_updated(api::GitHubAPI, lrepo, local_dir, before, after, commit_sig, repo, auth)
    before = LibGit2.GitCommit(lrepo, before)
    after = LibGit2.GitCommit(lrepo, after)
    delta = my_diff_tree(lrepo, LibGit2.peel(before), LibGit2.peel(after); pathspecs="REQUIRE")[1]
    old_blob = LibGit2.GitBlob(lrepo, delta.old_file.id)
    new_blob = LibGit2.GitBlob(lrepo, delta.new_file.id)
    vers_old = Pkg.Reqs.parse(IOBuffer(LibGit2.content(old_blob)))
    vers_new = Pkg.Reqs.parse(IOBuffer(LibGit2.content(new_blob)))
    if vers_new["julia"] != vers_old["julia"]
        apply_deprecations(api, lrepo, local_dir, commit_sig, repo, auth)
    end
end

function cleanrepo(repo_url; show_diff = true, delete_local = true)
    local_dir = mktempdir()
    successful = true
    try
        enabled = gc_enable(false)
        info("Cloning $repo_url to $local_dir...")
        lrepo = LibGit2.clone(repo_url, local_dir)
        gc_enable(enabled)
        info("Processing deprecations...")
        changed_any, problematic_files = process_deprecations(lrepo, local_dir; is_julia_itself=contains(repo_url, "JuliaLang/julia"))
        isempty(problematic_files) || (successful = false)
    catch e
        bt = catch_backtrace()
        Base.display_error(STDERR, e, bt)
        successful = false
    finally
        if show_diff
            cd(local_dir) do
                run(`git status`)
                run(`git diff --cached`)
            end
        end
        if delete_local
            info("Deleting cloned repo from $local_dir...")
            rm(local_dir, force=true, recursive=true)
        end
    end
    return successful
end

include("interactions.jl")
include("autodeployment.jl")

const autodeployment_enabled = haskey(ENV, "FEMTOCLEANER_AUTODEPLOY") ?
    ENV["FEMTOCLEANER_AUTODEPLOY"] == "yes" : false

let app_key = Ref{Any}(nothing)
    global get_auth
    function get_auth(app_id)
        if app_key[] == nothing
            app_key[] = MbedTLS.PKContext()
            MbedTLS.parse_key!(app_key[], haskey(ENV, "FEMTOCLEANER_PRIVKEY") ? ENV["FEMTOCLEANER_PRIVKEY"] : readstring(joinpath(dirname(@__FILE__),"..","privkey.pem")))
        end
        GitHub.JWTAuth(app_id, app_key[])
    end
end

function event_callback(api::GitHubAPI, app_name, app_id, sourcerepo_installation,
                        commit_sig, listener, bug_repository, event)
    # On installation, process every repository we just got installed into
    if event.kind == "installation"
        jwt = get_auth(app_id)
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        for repo in event.payload["repositories"]
            repo = GitHub.repo(api, GitHub.Repo(repo); auth=auth)
            with_cloned_repo(api, repo, auth) do x
                apply_deprecations(api, x..., commit_sig, repo, auth)
            end
        end
    elseif event.kind == "check_run"
        jwt = get_auth(app_id)
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        if event.payload["action"] == "requested_action" && event.payload["requested_action"]["identifier"] == "fix"
            repo = GitHub.Repo(event.payload["repository"])
            pr = PullRequest(event.payload["check_run"]["check_suite"]["pull_requests"][1])
            with_cloned_repo(api, repo, auth) do x
                lrepo, local_dir = x
                LibGit2.checkout!(lrepo, event.payload["check_run"]["check_suite"]["head_sha"])
                is_julia_itself = GitHub.name(repo) == "JuliaLang/julia"
                changed_any, problematic_files = process_deprecations(lrepo, local_dir; is_julia_itself=is_julia_itself)
                if changed_any
                    LibGit2.commit(lrepo, "Fix deprecations"; author=commit_sig, committer=commit_sig, parent_ids=[LibGit2.GitHash(lrepo, "HEAD")])
                    push_repo(api, lrepo, auth; force=false, remote_branch=event.payload["check_run"]["check_suite"]["head_branch"])
                end
            end
        end
    elseif event.kind == "pull_request"
        if !(event.payload["action"] in ("opened", "reopened", "synchronize"))
            return HTTP.Response(200)
        end
        jwt = get_auth(app_id)
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        repo = Repo(event.payload["repository"])
        pr = PullRequest(event.payload["pull_request"])
        local annotations
        with_cloned_repo(api, repo, auth) do x
            lrepo, local_dir = x
            LibGit2.checkout!(lrepo, get(get(pr.head).sha))
            annotations, _ = collect_deprecation_annotations(api, x..., repo, auth)
        end
        actions = Checks.Action[]
        conclusion = "neutral"
        if length(annotations) == 0
            message = """
            No applicable deprecations were detected.
            """
            output = Checks.Output(
                "Femtocleaning",
                message,
                "",
                Checks.Annotation[],
                Checks.Image[]
            )
            conclusion = "success"
        else
            truncated = length(annotations) > 50
            message = """
            Several femtocleaning opportunities were detected
            """
            output = Checks.Output(
                "Femtocleaning",
                message,
                "See below",
                annotations[1:min(50, length(annotations))],
                GitHub.Image[]
            )
            actions = Checks.Action[
                Checks.Action(
                    "Fix it!",
                    "Fixes issues in this PR (adds commit).",
                    "fix"
                )
            ]
        end
        max_annotation = 50
        cr = GitHub.create_check_run(api, repo, auth=auth, params = Dict(
            :name => "femtocleaner",
            :head_branch => get(pr.head).ref,
            :head_sha => get(pr.head).sha,
            :status => "completed",
            :conclusion => conclusion,
            :completed_at => now(),
            :actions => actions,
            :output => output
        ))
        while max_annotation < length(annotations)
            empty!(output.annotations)
            append!(output.annotations, annotations[max_annotation+1:min(max_annotation+50, end)])
            max_annotation += 50
            GitHub.update_check_run(api, repo, get(cr.id), auth=auth, params = Dict(
                :output => output,
                :actions => actions
            ))
        end
    elseif event.kind == "installation_repositories"
        jwt = get_auth(app_id)
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        for repo in event.payload["repositories_added"]
            repo = GitHub.repo(api, GitHub.Repo(repo); auth=auth)
            with_cloned_repo(api, repo, auth) do x
                apply_deprecations(api, x..., commit_sig, repo, auth)
            end
        end
    elseif event.kind == "pull_request_review"
        jwt = get_auth(app_id)
        pr_response(api, event, jwt, commit_sig, app_name, sourcerepo_installation, bug_repository)
    elseif event.kind == "push"
        jwt = get_auth(app_id)
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        repo = Repo(event.payload["repository"])
        # Check if REQUIRE was updated
        for commit in event.payload["commits"]
            if "REQUIRE" in commit["modified"]
                with_cloned_repo(api, repo, auth) do x
                    apply_deprecations_if_updated(api, x...,
                        event.payload["before"], event.payload["after"],
                        commit_sig, repo, auth)
                end
                break
            end
        end
        maybe_autdodeploy(event, listener, jwt, sourcerepo_installation, autodeployment_enabled)
    elseif event.kind == "issues" && event.payload["action"] == "opened"
        jwt = get_auth(app_id)
        iss = Issue(event.payload["issue"])
        repo = Repo(event.payload["repository"])
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        if lowercase(get(iss.title)) == "run femtocleaner"
            with_cloned_repo(api, repo, auth) do x
                apply_deprecations(api, x..., commit_sig, repo, auth; issue_number = get(iss.number))
            end
        end
    end
    return HTTP.Response(200)
end

function run_server()
    app_id = parse(Int, strip(haskey(ENV, "FEMTOCLEANER_APPID") ? ENV["FEMTOCLEANER_APPID"] : readstring(joinpath(dirname(@__FILE__),"..","app_id"))))
    secret = haskey(ENV, "FEMTOCLEANER_SECRET") ? ENV["FEMTOCLEANER_SECRET"] : nothing
    sourcerepo_installation = haskey(ENV, "FEMTOCLEANER_INSTALLATION") ? parse(Int, ENV["FEMTOCLEANER_INSTALLATION"]) : 0
    bug_repository = haskey(ENV, "FEMTOCLEANER_BUGREPO") ? ENV["FEMTOCLEANER_BUGREPO"] : ""
    (secret == nothing) && warn("Webhook secret not set. All events will be accepted. This is an insecure configuration!")
    jwt = get_auth(app_id)
    app_name = get(GitHub.app(; auth=jwt).name)
    commit_sig = LibGit2.Signature("$(app_name)[bot]", "$(app_name)[bot]@users.noreply.github.com")
    api = GitHub.DEFAULT_API
    local listener
    if nprocs() != 1
        #@everywhere filter(x->x != 1, procs()) worker_loop()
        for p in filter(x->x != 1, procs())
            @spawnat p begin
                myid() == 2 && update_existing_repos(api, commit_sig, app_id)
                worker_loop()
            end
        end
    end
    listener = GitHub.EventListener(secret=secret) do event
        queue() do
            revise()
            Base.invokelatest(event_callback, api, app_name, app_id,
                              sourcerepo_installation, commit_sig, listener, bug_repository, event)
        end
        HTTP.Response(200)
    end
    GitHub.run(listener, IPv4(0,0,0,0), 10000+app_id)
    wait()
end

end # module


text:
# print data/sensor_data_*.txt

module PrintSensorData
    using DataFrames, CSV

    function print_data_200()
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_200.txt")
        df_200_mm = CSV.read(data_path, DataFrame,
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')

        println("sensor_data_200.txt")
        println("Data size = ", size(df_200_mm))
        println("Column Names = ", names(df_200_mm))
        println("Data type = ", eltype.(eachcol(df_200_mm)))
        println("Describe = ", describe(df_200_mm))
        println("")

        return df_200_mm
    end

    function print_data_600()
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_600_mm = CSV.read(data_path, DataFrame,
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')

        println("sensor_data_600.txt")
        println("Data size = ", size(df_600_mm))
        println("Column Names = ", names(df_600_mm))
        println("Data type = ", eltype.(eachcol(df_600_mm)))
        println("Describe = ", describe(df_600_mm))
        println("")

        return df_600_mm
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# Test that freeze thaw alone reproduces expected behavior: exponential behavior
# for liquid water content, ice content, and total water conserved

using MPI
using OrderedCollections
using StaticArrays
using Statistics
using Test

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()
using CLIMAParameters.Planet: ρ_cloud_liq
using CLIMAParameters.Planet: ρ_cloud_ice

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Land.SoilHeatParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state

@testset "Freeze thaw alone" begin
    struct tmp_model <: BalanceLaw end
    struct tmp_param_set <: AbstractParameterSet end

    function get_grid_spacing(
        N_poly::Int64,
        nelem_vert::Int64,
        zmax::FT,
        zmin::FT,
    ) where {FT}
        test_config = ClimateMachine.SingleStackConfiguration(
            "TmpModel",
            N_poly,
            nelem_vert,
            zmax,
            tmp_param_set(),
            tmp_model();
            zmin = zmin,
        )

        Δ = min_node_distance(test_config.grid)
        return Δ
    end

    function init_soil_water!(land, state, aux, coordinates, time)
        state.soil.water.ϑ_l = eltype(state)(land.soil.water.initialϑ_l(aux))
        state.soil.water.θ_i = eltype(state)(land.soil.water.initialθ_i(aux))
    end

    FT = Float64
    ClimateMachine.init()

    N_poly = 5
    nelem_vert = 50
    zmax = FT(0)
    zmin = FT(-1)
    t0 = FT(0)
    timeend = FT(30)
    dt = FT(0.05)

    n_outputs = 30
    every_x_simulation_time = ceil(Int, timeend / n_outputs)
    Δ = get_grid_spacing(N_poly, nelem_vert, zmax, zmin)
    cs = FT(3e6)
    κ = FT(1.5)
    τLTE = FT(cs * Δ^FT(2.0) / κ)
    freeze_thaw_source = PhaseChange{FT}(Δt = dt, τLTE = τLTE)

    soil_param_functions =
        SoilParamFunctions{FT}(porosity = 0.75, Ksat = 0.0, S_s = 1e-3)

    bottom_flux = (aux, t) -> eltype(aux)(0.0)
    surface_flux = (aux, t) -> eltype(aux)(0.0)
    surface_state = nothing
    bottom_state = nothing

    bc = GeneralBoundaryConditions(
        Dirichlet(surface_state = surface_state, bottom_state = bottom_state),
        Neumann(surface_flux = surface_flux, bottom_flux = bottom_flux),
    )
    ϑ_l0 = (aux) -> eltype(aux)(1e-10)
    θ_i0 = (aux) -> eltype(aux)(0.33)

    soil_water_model = SoilWaterModel(
        FT;
        initialϑ_l = ϑ_l0,
        initialθ_i = θ_i0,
        boundaries = bc,
    )

    soil_heat_model =
        PrescribedTemperatureModel((aux, t) -> eltype(aux)(276.15))

    m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model)
    sources = (freeze_thaw_source,)
    m = LandModel(
        param_set,
        m_soil;
        source = sources,
        init_state_prognostic = init_soil_water!,
    )

    driver_config = ClimateMachine.SingleStackConfiguration(
        "LandModel",
        N_poly,
        nelem_vert,
        zmax,
        param_set,
        m;
        zmin = zmin,
        numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )
    state_types = (Prognostic(),)
    dons_arr =
        Dict[dict_of_nodal_states(solver_config, state_types; interp = true)]
    time_data = FT[0]
    callback = GenericCallbacks.EveryXSimulationTime(every_x_simulation_time) do
        dons = dict_of_nodal_states(solver_config, state_types; interp = true)
        push!(dons_arr, dons)
        push!(time_data, gettime(solver_config.solver))
        nothing
    end
    ClimateMachine.invoke!(solver_config; user_callbacks = (callback,))
    dons = dict_of_nodal_states(solver_config, state_types; interp = true)
    push!(dons_arr, dons)
    push!(time_data, gettime(solver_config.solver))
    m_liq = [
        ρ_cloud_liq(param_set) * mean(dons_arr[k]["soil.water.ϑ_l"])
        for k in 1:(n_outputs + 1)
    ]
    m_ice = [
        ρ_cloud_ice(param_set) * mean(dons_arr[k]["soil.water.θ_i"])
        for k in 1:(n_outputs + 1)
    ]
    total_water = m_ice + m_liq
    τft = max(dt, τLTE)
    m_ice_of_t = m_ice[1] * exp.(-1.0 .* (time_data .- time_data[1]) ./ τft)
    m_liq_of_t = -m_ice_of_t .+ (m_liq[1] + m_ice[1])
    @test mean(abs.(m_liq .- m_liq_of_t)) < 1e-9
    @test mean(abs.(m_ice .- m_ice_of_t)) < 1e-9
end

max_forks_repo_name: jm-c/CLIMA


text:
const NDIMS = 2

type Body
	mass::Float64
	pos::Vector{Float64}
	vel::Vector{Float64}

	initial_energy

	step
	prev_accel::Array{Float64, 2}

	function Body(mass=0.0, pos=Vector{Float64}(NDIMS), vel=Vector{Float64}(NDIMS))
		accel_hist_len = 4
		accel_hist = zeros(Float64, NDIMS, accel_hist_len)

		new(mass, pos, vel, 0, 0, accel_hist)
	end
end

######## IO ###########

import Base.show
function show(io::IO, a::Body)
	print(io, "mass: ", a.mass, "\n")
	print(io, "pos: ", join(a.pos, ", "), "\n")
	print(io, "vel: ", join(a.vel, ", "), "\n")
end

import Base.print
function print(io::IO, a::Body)
	@printf io "%24.16e\n" a.mass

	map(Array[a.pos, a.vel]) do b
		map(x-> @printf(io,"%24.16e",x), b)
		print(io, "\n")
	end
end

function read_body()
	mass = parse(Float64, readline())
	pos = map(x -> parse(Float64, x), split(readline()))
	vel = map(x -> parse(Float64, x), split(readline()))
	return Body(mass, pos, vel)
end

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

function accel(body::Body)
	r2 = dot(body.pos, body.pos)
	r3 = r2 * sqrt(r2) # we want r^(2/3)
	acc = body.pos * (-body.mass / r3)
	return acc
end

function jerk(body::Body)
	r2 = dot(body.pos, body.pos)
	r3 = r2 * sqrt(r2) # we want r^(2/3)
	j = (body.vel+body.pos*(-3*(dot(body.pos,body.vel))/r2))*(-body.mass/r3)
	return j
end

function prev_accel(body::Body, n::Int)
	return body.prev_accel[:,n]
end

function set_prev_accel!(body::Body, n::Int, prev::Vector)
	body.prev_accel[:,n] = prev
end

####### Energy #########

function kin_energy(body::Body)
	return 0.5 * dot(body.vel, body.vel)
end

function pot_energy(body::Body)
	return -body.mass / sqrt(dot(body.pos, body.pos))
end

function total_energy(body::Body)
	return kin_energy(body) + pot_energy(body)
end

function init_energy!(body::Body)
	body.initial_energy = total_energy(body)
end


text:
export newton

"""
newton(f,df,H,x)

Newton method for solving min_x f(x)

"""
function newton(f::Function,df::Function,H::Function,x::Vector;maxIter=20,atol=1e-8,out::Int=0,storeInterm::Bool=false)

    his = zeros(maxIter,2)
    X = (storeInterm) ? zeros(length(x),maxIter) : []
    i = 1; flag = -1; LL = []
    while i<=maxIter
        fk  = f(x)
        dfk = df(x)
        his[i,:] = [fk norm(dfk)]
        if storeInterm; X[:,i] = x; end;

        if(norm(dfk)<atol)
            flag = 0
            his = his[1:i,:]
            break
        end

        # get search direction
        d2f = H(x)
        try
            LL    = cholfact(d2f)
        catch err
            if isa(err,Base.LinAlg.PosDefException)
                flag = -2
                his = his[1:i,:]
                break
            else
                throw(err)
            end
        end
        pk    = - (LL\dfk)
        if out>0
            @printf "iter=%04d\tf=%1.2e\t|df|=%1.2e\n" i his[i,1] his[i,2]
        end
        # update
        x  += pk
        i+=1
    end
    i = min(maxIter,i)

    if out>=0
        if flag==-1
            println(@sprintf("newton iterated maxIter (=%d) times but reached only atol of %1.2e instead of tol=%1.2e",i,his[i,2],atol))
        elseif flag==-2
            println(@sprintf("newton stopped because the Hessian at iteration %d was not positive definite.",i))
        elseif out>1
            println(@sprintf("newton achieved desired atol of %1.2e at iteration %d.",atol,i))
        end
    end
    if storeInterm; X = X[:,1:i]; end
    return x,flag,his,X
end


text:
# This file is part of TaylorSeries.jl, MIT licensed
#

using TaylorSeries, IntervalArithmetic

using Test
eeuler = Base.MathConstants.e

@testset "Tests Taylor1 and TaylorN expansions over Intervals" begin
    a = 1..2
    b = -1 .. 1
    p4(x, a) = x^4 + 4*a*x^3 + 6*a^2*x^2 + 4*a^3*x + a^4
    p5(x, a) = x^5 + 5*a*x^4 + 10*a^2*x^3 + 10*a^3*x^2 + 5*a^4*x + a^5

    ti = Taylor1(Interval{Float64}, 10)
    x, y = set_variables(Interval{Float64}, "x y")

    @test eltype(ti) == Interval{Float64}
    @test eltype(x) == Interval{Float64}

    @test p4(ti,-a) == (ti-a)^4
    @test p5(ti,-a) == (ti-a)^5
    @test p4(ti,-b) == (ti-b)^4
    @test all((p5(ti,-b)).coeffs .⊆ ((ti-b)^5).coeffs)

    @test p4(x,-y) == (x-y)^4
    @test p5(x,-y) == (x-y)^5
    @test p4(x,-a) == (x-a)^4
    @test p5(x,-a) == (x-a)^5
    @test p4(x,-b) == (x-b)^4
    for ind in eachindex(p5(x,-b))
        @test all((p5(x,-b)[ind]).coeffs .⊆ (((x-b)^5)[ind]).coeffs)
    end

    # Tests `evaluate`
    @test evaluate(p4(x,y), IntervalBox(a,-b)) == p4(a, -b)
    @test (p5(x,y))(IntervalBox(a,b)) == p5(a, b)
    @test (a-b)^4 ⊆ ((x-y)^4)(a × b)
    @test (((x-y)^4)[4])(a × b) == -39 .. 81

    p4n = normalize_taylor(p4(x,y), a × b, true)
    @test (0..16) ⊆ p4n((-1..1)×(-1..1))
    p5n = normalize_taylor(p5(x,y), a × b, true)
    @test (-32 .. 32) ⊆ p5n((-1..1)×(-1..1))

    p4n = normalize_taylor(p4(x,y), a × b, false)
    @test (0..16) ⊆ p4n((0..1)×(0..1))
    p5n = normalize_taylor(p5(x,y), a × b, false)
    @test (0..32) ⊆ p5n((0..1)×(0..1))

    @test evaluate(x*y^3, (-1..1)×(-1..1)) == (-1..1)
    @test evaluate(x*y^2, (-1..1)×(-1..1)) == (-1..1)
    @test evaluate(x^2*y^2, (-1..1)×(-1..1)) == (0..1)

    ii = -1..1
    t = Taylor1(1)
    @test 0..2 ⊆ (1+t)(ii)
    t = Taylor1(2)
    @test 0..4 ⊆ ((1+t)^2)(ii)

    ii = 0..6
    t = Taylor1(4)
    f(x) = 0.1 * x^3 - 0.5*x^2 + 1
    ft = f(t)
    f1 = normalize_taylor(ft, ii, true)
    f2 = normalize_taylor(ft, ii, false)
    @test Interval(-23/27, f(6)) ⊆ f(ii)
    @test Interval(-23/27, f(6)) ⊆ ft(ii)
    @test Interval(-23/27, f(6)) ⊆ f1(-1..1)
    @test Interval(-23/27, f(6)) ⊆ f2(0..1)
    @test f1(-1..1) ⊆ f(ii)
    @test diam(f1(-1..1)) < diam(f2(0..1))

    # An example from Makino's thesis
    ii = 0..1
    t = Taylor1(5)
    g(x) = 1 - x^4 + x^5
    gt = g(t)
    g1 = normalize_taylor(gt, 0..1, true)
    @test Interval(g(4/5),1) ⊆ g(ii)
    @test Interval(g(4/5),1) ⊆ gt(ii)
    @test Interval(g(4/5),1) ⊆ g1(-1..1)
    @test g1(-1..1) ⊂ g(ii)
    @test diam(g1(-1..1)) < diam(gt(ii))
end


text:
############################################################
# Split Layer
############################################################
@defstruct SplitLayer Layer (
  name :: AbstractString = "split",
  no_copy :: Bool = false,
  (tops :: Vector{Symbol} = Symbol[], length(tops) >= 1),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == 1),
)
@characterize_layer(SplitLayer,
  can_do_bp => true
)

type SplitLayerState{N} <: LayerState
  layer      :: SplitLayer
  blobs      :: Vector{Blob}
  blobs_diff :: Vector{Blob}
end

function setup(backend::Backend, layer::SplitLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  N = length(layer.tops)

  if layer.no_copy
    # directly re-use the input blob
    blobs = Blob[inputs[1] for i = 1:N]
  else
    blobs = Array{Blob}(N)
    blobs[1] = inputs[1]
    for i = 2:N
      blobs[i] = make_blob(backend, eltype(inputs[1]), size(inputs[1]))
    end
  end
  blobs_diff = Array{Blob}(N)
  blobs_diff[1] = diffs[1] # re-use the first backward blob
  for i = 2:N
    if isa(diffs[1], NullBlob)
      blobs_diff[i] = NullBlob()
    else
      blobs_diff[i] = make_blob(backend, eltype(inputs[1]), size(inputs[1]))
    end
  end

  return SplitLayerState{N}(layer, blobs, blobs_diff)
end
function shutdown(backend::Backend, state::SplitLayerState)
  # some blobs are shared, but never mind, blob destroy function has
  # a guard that does not cause problems on double destroying
  map(destroy, state.blobs)
  map(destroy, state.blobs_diff)
end

function forward(backend::Backend, state::SplitLayerState, inputs::Vector{Blob})
  if !state.layer.no_copy
    for i = 2:length(state.blobs)
      copy!(state.blobs[i], inputs[1])
    end
  end
end

function backward{N}(backend::CPUBackend, state::SplitLayerState{N}, inputs::Vector{Blob}, diffs::Vector{Blob})
  if !isa(diffs[1], NullBlob)
    diff = diffs[1]
    len = length(diff)
    one = convert(eltype(diff), 1)
    for i = 2:N
      BLAS.axpy!(len, one, state.blobs_diff[i].data, 1, diff.data, 1)
    end
  end
end