text:
################################################################
#                  Triangular surface facet
################################################################

"""
    struct Face{T1, T2, T3, T4, T5, T6}

Triangular surface facet of a polyhedral shape model.

Note that the mesh normal indicates outward the polyhedron.

# Fields
- `A::T1` : Position of 1st vertex
- `B::T1` : Position of 2nd vertex
- `C::T1` : Position of 3rd vertex

- `center::T1` : Position of mesh center
- `normal::T1` : Normal vector to mesh
- `area  ::T2`   : Area of mesh

- `visiblefacets::T3` : 1-D array of `VisibleFacet`
- `flux         ::T4` : Energy flux from surrounding facets
- `Tz           ::T5` : Temperature profile in depth direction
- `force        ::T6` : Photon recoil force
"""
struct Facet{T1, T2, T3, T4, T5, T6}
    A::T1
    B::T1
    C::T1

    center::T1
    normal::T1
    area  ::T2

    visiblefacets::T3
    flux         ::T4
    Tz           ::T5
    force        ::T6
end

Facet(A, B, C) = Facet([A, B, C])
Facet(vs) = Facet(
    vs[1], vs[2], vs[3],
    getcenter(vs), getnormal(vs), getarea(vs),
    StructArray(VisibleFacet[]), Flux(), Float64[], zeros(3)
)

function Base.show(io::IO, facet::Facet)
    println(io, "Surface facet")
    println("-------------")

    println("Vertices")
    println("    A : ", facet.A)
    println("    B : ", facet.B)
    println("    C : ", facet.C)

    println("Center : ", facet.center)
    println("Normal : ", facet.normal)
    println("Area   : ", facet.area)

    if length(facet.visiblefacets) == 0
        println("No visible facets.")
    else
        @unpack visiblefacets = facet
        println(length(visiblefacets), " facets are visible:")
        df = DataFrame(id=visiblefacets.id, f=visiblefacets.f, d=visiblefacets.d)
        println(df)
    end
end

"""
Array of `Facet`, converted from arrays of nodes and faces of a shape model
"""
getfacets(nodes, faces) = StructArray([Facet(nodes[face]) for face in faces])

# ################################################################
# #                      Facet properties
# ################################################################

getcenter(vs) = getcenter(vs...)
getcenter(v1, v2, v3) = (v1 + v2 + v3) / 3

getnormal(vs) = getnormal(vs...)
getnormal(v1, v2, v3) = normalize((v2 - v1) × (v3 - v2))

getarea(vs) = getarea(vs...)
getarea(v1, v2, v3) = norm((v2 - v1) × (v3 - v2)) * 0.5

################################################################
#                 Face-to-face interactions
################################################################

"""
    struct VisibleFacet{T1, T2, T3}

Index of an interfacing facet and its view factor

# Fields
- `id` : Index of the interfacing facet
- `f`  : View factor from facet i to j
- `d`  : Distance from facet i to j
- `d̂`  : Normal vector from facet i to j
"""
struct VisibleFacet
    id::Int64
    f ::Float64
    d ::Float64
    d̂ ::SVector{3, Float64}
end

VisibleFacet(i::Facet, j::Facet, id) = VisibleFacet(id, view_factor(i, j)...)

"""
    view_factor(i::Facet, j::Facet) -> fᵢⱼ, dᵢⱼ, d̂ᵢⱼ

View factor from facet i to j, assuming Lambertian emission
"""
function view_factor(i::Facet, j::Facet)
    d⃗ᵢⱼ = j.center - i.center  # vector from facet i to j
    dᵢⱼ = norm(d⃗ᵢⱼ)
    d̂ᵢⱼ = normalize(d⃗ᵢⱼ)

    cosθᵢ = i.normal ⋅ d̂ᵢⱼ
    cosθⱼ = j.normal ⋅ (-d̂ᵢⱼ)

    fᵢⱼ = view_factor(cosθᵢ, cosθⱼ, dᵢⱼ, j.area)
    fᵢⱼ, dᵢⱼ, d̂ᵢⱼ
end

view_factor(cosθᵢ, cosθⱼ, dᵢⱼ, aⱼ) = cosθᵢ * cosθⱼ / (π * dᵢⱼ^2) * aⱼ

"""
    mutable struct Flux{T}

Energy flux from/to a facet

# Fields
- `sun ::T`  : Flux of solar radiation,   F_sun
- `scat::T` : Flux of scattered sunlight, F_scat
- `rad ::T`  : Flux of thermal radiation, F_rad
"""
mutable struct Flux{T}
    sun ::T
    scat::T
    rad ::T
end

Flux() = Flux(0., 0., 0.)

################################################################
#                           Orinet3D
################################################################

"""
    isFace(obs::Facet, tar::Facet) -> true/false

Determine if the two facets are facing each other
"""
isFace(obs::Facet, tar::Facet) = (tar.center - obs.center) ⋅ tar.normal < 0 ? true : false

"""
    isAbove(A, B, C, D)             -> Bool
    isAbove(facet::Facet, D)        -> Bool
    isAbove(obs::Facet, tar::Facet) -> Bool

Determine if point D is above triangle facet ABC.
"""
function isAbove(A, B, C, D)
    G = SA_F64[
        A[1]-D[1] A[2]-D[2] A[3]-D[3]
        B[1]-D[1] B[2]-D[2] B[3]-D[3]
        C[1]-D[1] C[2]-D[2] C[3]-D[3]
    ]

    det(G) < 0 ? true : false
end

"""
    isBelow(A, B, C, D)             -> Bool
    isBelow(facet::Facet, D)        -> Bool
    isBelow(obs::Facet, tar::Facet) -> Bool

Determine if point D is below triangle facet ABC.
"""
function isBelow(A, B, C, D)
    G = SA_F64[
        A[1]-D[1] A[2]-D[2] A[3]-D[3]
        B[1]-D[1] B[2]-D[2] B[3]-D[3]
        C[1]-D[1] C[2]-D[2] C[3]-D[3]
    ]

    det(G) > 0 ? true : false
end

isAbove(facet::Facet, D) = isAbove(facet.A, facet.B, facet.C, D)
isBelow(facet::Facet, D) = isBelow(facet.A, facet.B, facet.C, D)

isAbove(obs::Facet, tar::Facet) = isAbove(obs.A, obs.B, obs.C, tar.center)
isBelow(obs::Facet, tar::Facet) = isBelow(obs.A, obs.B, obs.C, tar.center)

################################################################
#                           Raycast
################################################################

"""
    raycast(A, B, C, R)                           -> Bool
    raycast(facet::Facet, R)                      -> Bool
    raycast(facet::Facet, R, obs::AbstractVector) -> Bool
    raycast(facet::Facet, R, obs::Facet)          -> Bool

Ray-triangle intersection detection.
The ray R starts from the origin.
"""
function raycast(A, B, C, R)
    E1 = B - A
    E2 = C - A
    T  = - A

    P = R × E2
    Q = T × E1

    P_dot_E1 = P ⋅ E1

    u = (P ⋅ T)  / P_dot_E1
    v = (Q ⋅ R)  / P_dot_E1
    t = (Q ⋅ E2) / P_dot_E1  # Distance to triangle (if t < 0, the ray hits from the back)

    0 ≤ u ≤ 1 && 0 ≤ v ≤ 1 && 0 ≤ u + v ≤ 1 && t > 0 ? true : false
end

raycast(facet::Facet, R) = raycast(facet.A, facet.B, facet.C, R)
raycast(facet::Facet, R, obs::AbstractVector) = raycast(facet.A - obs, facet.B - obs, facet.C - obs, R)
raycast(facet::Facet, R, obs::Facet) = raycast(facet, R, obs.center)

"""
    findVisibleFacets!(obs::Facet, facets)

Find facets that is visible from the facet where the observer is located.

# Parameters
- `obs`    : Facet where the observer stands
- `facets` : Array of `Facet`
"""
function findVisibleFacets!(obs::Facet, facets)
    ids = Int64[]
    for (id, tar) in enumerate(facets)
        isAbove(obs, tar) && isFace(obs, tar) && push!(ids, id)
    end

    ii = copy(ids)

    for i in ii
        tar_i = facets[i]
        Rᵢ = tar_i.center - obs.center
        dᵢ = norm(Rᵢ)      # distance to facet i
        for j in ii
            i == j && continue

            tar_j = facets[j]
            Rⱼ = tar_j.center - obs.center
            dⱼ = norm(Rⱼ)  # distance to facet j

            raycast(tar_j, Rᵢ, obs) && (dᵢ < dⱼ ? filter!(x->x≠j, ids) : filter!(x->x≠i, ids))
        end
    end

    for id in ids
        tar = facets[id]
        push!(obs.visiblefacets, VisibleFacet(obs, tar, id))
    end
end

"""
    findVisibleFacets!(facets)

Find facets that is visible from each facet

# Parameters
- `facets` : Array of `Facet`
"""
function findVisibleFacets!(facets)
    for obs in facets
        findVisibleFacets!(obs, facets)
    end
end

"""
    isIlluminated(obs::AbstractVector, r̂☉, facets)
    isIlluminated(obs::Facet, r̂☉, facets) -> Bool

Return if the observation point/facet is illuminated by the direct sunlight or not
"""
function isIlluminated(obs::AbstractVector, r̂☉, facets)
    for facet in facets
        raycast(facet, r̂☉, obs) && return false
    end
    return true
end

function isIlluminated(obs::Facet, r̂☉, facets)
    obs.normal ⋅ r̂☉ < 0 && return false
    for id in obs.visiblefacets.id
        raycast(facets[id], r̂☉, obs) && return false
    end
    return true
end

"""
    isAboveHorizon(facet::Facet) -> Bool

Return if the facet is above its local horizon or not
"""
isAboveHorizon(facet::Facet) = length(facet.visiblefaces) == 0

"""
    getIlluminatedFacets(r̂☉, facets; ray_trace=true) -> ids

# Arguments
- `r̂☉`        : Normal vector to indicate solar direction
- `facets`    : Array of `Facet`
- `ray_trace` : Option to turn on/off of ray-trace

# Return
- `ids` : Indices of illuminated facets
"""
function getIlluminatedFacets(r̂☉, facets; ray_trace=true)
    ids = Int64[]
    for (id, obs) in enumerate(facets)
        ray_trace == true  && isIlluminated(obs, r̂☉, facets) && push!(ids, id)
        ray_trace == false && obs.normal ⋅ r̂☉ > 0            && push!(ids, id)  # Pseudo-convex
    end
    ids
end

################################################################
#                    Solid angle of facet
################################################################

"""
    getSolidAngle(facet::Facet, obs::AbstractVector) -> Ω

Solid angle of a triangular facet,
equal to area of the corresponding spherical triangle
"""
function getSolidAngle(facet::Facet, obs::AbstractVector)

    ## Vectors from observer to facet vertices
    A = facet.A - obs
    B = facet.B - obs
    C = facet.C - obs

    AOB = getangle(A, B)
    BOC = getangle(B, C)
    COA = getangle(C, A)

    Ω = getSphericalExcess(AOB, BOC, COA)
end

"""
    getSphericalExcess(a, b, c) -> E

Area of a spherical triangle
c.f. L'Huilier's Theorem

# Arguments
- `a`, `b`, `c` : side lengths of a spherical traiangle
"""
function getSphericalExcess(a, b, c)
    s = (a + b + c) * 0.5  # semiperimeter

    E = tan(s*0.5)
    E *= tan((s - a)*0.5)
    E *= tan((s - b)*0.5)
    E *= tan((s - c)*0.5)
    E = sqrt(E)
    E = 4 * atan(E)
end

getangle(v1, v2) = acos(normalize(v1) ⋅ normalize(v2))

max_stars_repo_name: MasanoriKanamaru/Astroshaper

max_issues_repo_name: MasanoriKanamaru/Astroshaper

max_forks_repo_name: MasanoriKanamaru/Astroshaper


text:
using LeastSquaresSVM
using Test
using MLJBase, MLJ, MLJModels

@testset "LeastSquaresSVM.jl" begin
    include("typetests.jl")
    include("trainingtests.jl")
    include("integrationtests.jl")
end

ext: jl

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
using Test, BioSequences

@testset "Encoding" begin
    seq = LongDNA{4}("ACGTA")
    @test wgregseq.utils.onehot_encoder(seq) == [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1; 1 0 0 0]
    seq = LongDNA{4}("ACGTN")
    @test_throws ArgumentError wgregseq.utils.onehot_encoder(seq)
end

@testset "Energy Matrix" begin
    seq1 = LongDNA{4}("ACGTAA")
    seq2 = LongDNA{4}("AAAAAA")
    emat = zeros(6, 4)
    emat[:, 1] .= 1

    @test wgregseq.utils.eval_emat(seq1, emat) == 3
    @test wgregseq.utils.eval_emat(seq2, emat) == 6

    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, zeros(5, 4))
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, zeros(6, 3))

    emat = zeros(4, 4)
    emat[:, 1] .= 1
    @test wgregseq.utils.eval_emat(seq2, emat, start_ind_seq=1, stop_ind_seq=4) == 4

    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_seq=1)
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_seq=-1, stop_ind_seq=4)
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_seq=-1, stop_ind_seq=4)
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_seq=3, stop_ind_seq=4)
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_seq=3, stop_ind_seq=4)

    emat = zeros(6, 4)
    emat[:, 1] .= 1
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_emat=-1, stop_ind_emat=4)
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_emat=-1, stop_ind_emat=4)
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_emat=3, stop_ind_emat=4)
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_emat=3, stop_ind_emat=4)

    @test wgregseq.utils.eval_emat(seq2, emat, start_ind_seq=3, stop_ind_seq=4, start_ind_emat=3, stop_ind_emat=4) == 2
    @test_throws ArgumentError wgregseq.utils.eval_emat(seq2, emat, start_ind_seq=3, stop_ind_seq=4, start_ind_emat=2, stop_ind_emat=4)

    seq_string = "ACGTAC"
    @test wgregseq.utils.eval_emat(seq_string, emat) == 2
    @test typeof(seq_string) == String

end


text:
export read_H5SeisCon

function read_H5SeisCon(con::H5SeisCon, blocks::Array{Int,1})
  if isnothing(con.seis)
    @error "Seis is Nothing"
    return
  end

  if isnothing(con.pkey) || isnothing(con.xkey) || isnothing(con.ykey) || isnothing(con.zkey)
    @error "pkey/xkey/ykey/zkey is Nothing"
    return
  end

  if isnothing(con.pkeyvals)
    @error "pkeyvals is Nothing"
    return
  end

  ind = Array{Integer,1}(undef, 0)
  # data = Array{Array{Float32,2},1}(undef, length(blocks))
  data = Array{Float32,2}(undef, con.seis.getNSamp(), 0)  # without preallocating rows I'm unable to `hcat`
  keylist = [con.pkey]
  for block in blocks
    minlist = [con.pkeyvals[block]]
    maxlist = [con.pkeyvals[block]]
    TRACE, _, block_ind = con.seis.getSortedData(keylist, minlist, maxlist)
    # data[block], _, block_ind = con.seis.getSortedData(keylist, minlist, maxlist)
    ind = vcat(ind, block_ind)
    data = hcat(data, TRACE)
  end

  return H5Seis(con.seis, ind, data)
end

# function read_H5SeisCon(con::H5SeisCon, blocks::Array{Int,1})
#   if isnothing(con.seis)
#     @error "Seis is Nothing"
#     return
#   end

#   if isnothing(con.pkey) || isnothing(con.xkey) || isnothing(con.ykey) || isnothing(con.zkey)
#     @error "pkey/xkey/ykey/zkey is Nothing"
#     return
#   end

#   if isnothing(con.pkeyvals)
#     @error "pkeyvals is Nothing"
#     return
#   end

#   data = Array{Array{Float32,2},1}(undef, length(blocks))
#   keylist = [con.pkey]
#   for block in blocks
#     minlist = [con.pkeyvals[block]]
#     maxlist = [con.pkeyvals[block]]
#     data[block], _, _ = con.seis.getSortedData(keylist, minlist, maxlist)
#   end

#   return data
# end

# RANGES & INT
function read_H5SeisCon(con::H5SeisCon, blocks::TR) where {TR<:AbstractRange}
  read_H5SeisCon(con, Array(blocks))
end
function read_H5SeisCon(con::H5SeisCon, blocks::Integer)
  read_H5SeisCon(con, [blocks])
end


text:
Base.Broadcast.broadcastable(l::SupervisedLoss) = Ref(l)

# --------------------------------------------------------------
# aggmode support

for (FUN, DESC, EXAMPLE) in (
        (:value,
         "result of the loss function",
         """
         3-element Array{Float64,1}:
           1.0
           9.0
          25.0
         """
        ),
        (:deriv,
         "derivative of the loss function in respect to the output",
         """
         3-element Array{Float64,1}:
            2.0
            6.0
          -10.0
         """
        ),
        (:deriv2,
         "second derivative of the loss function in respect to the output",
         """
         3-element Array{Float64,1}:
          2.0
          2.0
          2.0
         """
        ))
    @eval begin

        # ------------------------------------------------------
        # FALLBACK

        # By default compute the element-wise result
        @doc """
            $($(string(FUN)))(loss, targets::AbstractArray, outputs::AbstractArray)

        Compute the $($(DESC)) for each index-pair in `targets`
        and `outputs` individually and return the result as an
        array of the appropriate size.

        In the case that the two parameters are arrays with a
        different number of dimensions, broadcast will be
        performed. Note that the given parameters are expected to
        have the same size in the dimensions they share.

        Note: This function should always be type-stable. If it
        isn't, you likely found a bug.

        # Arguments

        - `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

        - `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

        - `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

        # Examples

        ```jldoctest
        julia> $($(string(FUN)))(L2DistLoss(), [1.0, 2.0, 3.0], [2, 5, -2])
        $($(EXAMPLE))
        ```
        """
        @inline function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray,
                output::AbstractArray)
            ($FUN)(loss, target, output, AggMode.None())
        end

        # (mutating) By default compute the element-wise result
        @doc """
            $($(string(FUN)))!(buffer::AbstractArray, loss, targets::AbstractArray, outputs::AbstractArray) -> buffer

        Compute the $($(DESC)) for each index-pair in `targets`
        and `outputs` individually, and store them in the
        preallocated `buffer`. Note that `buffer` has to be of
        the appropriate size.

        In the case that the two parameters, `targets` and
        `outputs`, are arrays with a different number of
        dimensions, broadcast will be performed. Note that the
        given parameters are expected to have the same size in
        the dimensions they share.

        Note: This function should always be type-stable. If it
        isn't, you likely found a bug.

        # Arguments

        - `buffer::AbstractArray`: Array to store the computed
          values in. Old values will be overwritten and lost.

        - `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

        - `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

        - `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

        # Examples

        ```jldoctest
        julia> buffer = zeros(3); # preallocate a buffer

        julia> $($(string(FUN)))!(buffer, L2DistLoss(), [1.0, 2.0, 3.0], [2, 5, -2])
        $($(EXAMPLE))
        ```
        """
        @inline function ($(Symbol(FUN,:!)))(
                buffer::AbstractArray,
                loss::SupervisedLoss,
                target::AbstractArray,
                output::AbstractArray)
            ($(Symbol(FUN,:!)))(buffer, loss, target, output, AggMode.None())
        end

        # Translate ObsDim.Last to the correct ObsDim.Constant (for code reduction)
        @inline function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray,
                output::AbstractArray{T,N},
                avg::AggregateMode,
                ::ObsDim.Last = ObsDim.Last()) where {T,N}
            ($FUN)(loss, target, output, avg, ObsDim.Constant{N}())
        end

        # (mutating) Translate ObsDim.Last to the correct ObsDim.Constant (for code reduction)
        @inline function ($(Symbol(FUN,:!)))(
                buffer::AbstractArray,
                loss::SupervisedLoss,
                target::AbstractArray,
                output::AbstractArray{T,N},
                avg::AggregateMode,
                ::ObsDim.Last = ObsDim.Last()) where {T,N}
            ($(Symbol(FUN,:!)))(buffer, loss, target, output, avg, ObsDim.Constant{N}())
        end

        # ------------------------------------------------------
        # ELEMENT-WISE BROADCAST

        # Compute element-wise (returns an array)
        @generated function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,M},
                output::AbstractArray{T,N},
                ::AggMode.None) where {Q,M,T,N}
            quote
                $(Expr(:meta, :inline))
                S = typeof(($($FUN))(loss, one(Q), one(T)))
                ($($FUN)).(Ref(loss), target, output)
            end
        end

        # (mutating) Compute element-wise (returns an array)
        function ($(Symbol(FUN,:!)))(
                buffer::AbstractArray,
                loss::SupervisedLoss,
                target::AbstractArray{Q,M},
                output::AbstractArray{T,N},
                ::AggMode.None) where {Q,M,T,N}
            buffer .= ($FUN).(Ref(loss), target, output)
            buffer
        end

        # ------------------------------------------------------
        # REDUCE TO NUMBER - BROADCAST

        # Compute the mean (returns a Number)
        @generated function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,M},
                output::AbstractArray{T,N},
                ::AggMode.Mean) where {Q,M,T,N}
            bigger = M > N ? :target : :output
            S, B = min(M,N), max(M,N)
            P = promote_type(Q,T)
            quote
                @nexprs $S (n)->@dimcheck(size(target, n) == size(output, n))
                nrm = 1 / $P(length($bigger))
                out = zero(($($FUN))(loss, one(Q), one(T)) * nrm)
                @inbounds @simd for I in CartesianIndices(size($bigger))
                    @nexprs $B n->(i_n = I[n])
                    out += ($($FUN))(loss, @nref($M,target,i), @nref($N,output,i)) * nrm
                end
                out
            end
        end

        # Compute the sum (returns a Number)
        @generated function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,M},
                output::AbstractArray{T,N},
                ::AggMode.Sum) where {Q,M,T,N}
            bigger = M > N ? :target : :output
            S, B = min(M,N), max(M,N)
            quote
                @nexprs $S (n)->@dimcheck(size(target, n) == size(output, n))
                out = zero(($($FUN))(loss, one(Q), one(T)))
                @inbounds @simd for I in CartesianIndices(size($bigger))
                    @nexprs $B n->(i_n = I[n])
                    out += ($($FUN))(loss, @nref($M,target,i), @nref($N,output,i))
                end
                out
            end
        end

        # ------------------------------------------------------
        # REDUCE TO NUMBER - SAME SIZE

        # Compute the total weighted mean (returns a Number)
        function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                avg::AggMode.WeightedMean,
                ::ObsDim.Constant{O}) where {Q,T,N,O}
            O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
            @dimcheck size(target) == size(output)
            @dimcheck size(output, O) == length(avg.weights)
            k = prod(n != O ? size(output,n) : 1 for n in 1:N)::Int
            nrm = avg.normalize ? inv(k * sum(avg.weights)) : inv(k * one(sum(avg.weights)))
            out = zero(($FUN)(loss, one(Q), one(T)) * (avg.weights[1] * nrm))
            @inbounds @simd for I in CartesianIndices(size(output))
                out += ($FUN)(loss, target[I], output[I]) * (avg.weights[I[O]] * nrm)
            end
            out
        end

        # Compute the total weighted sum (returns a Number)
        function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                avg::AggMode.WeightedSum,
                ::ObsDim.Constant{O}) where {Q,T,N,O}
            O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
            @dimcheck size(target) == size(output)
            @dimcheck size(output, O) == length(avg.weights)
            nrm = avg.normalize ? inv(sum(avg.weights)) : inv(one(sum(avg.weights)))
            out = zero(($FUN)(loss, one(Q), one(T)) * (avg.weights[1] * nrm))
            @inbounds @simd for I in CartesianIndices(size(output))
                out += ($FUN)(loss, target[I], output[I]) * (avg.weights[I[O]] * nrm)
            end
            out
        end

        # ------------------------------------------------------
        # PER OBSERVATION - SAME SIZE - TO VECTOR

        # Compute the mean per observation (returns a Vector)
        function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                avg::AggMode.Mean,
                obsdim::ObsDim.Constant{O}) where {Q,T,N,O}
            S = typeof(($FUN)(loss, one(Q), one(T)) / one(Int))
            buffer = zeros(S, size(output, O))
            ($(Symbol(FUN,:!)))(buffer, loss, target, output, avg, obsdim)
        end

        # (mutating) Compute the mean per observation (returns a Vector)
        function ($(Symbol(FUN,:!)))(
                buffer::AbstractVector{B},
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                ::AggMode.Mean,
                ::ObsDim.Constant{O}) where {B,Q,T,N,O}
            N == 1 && throw(ArgumentError("Mean per observation non sensible for two Vectors. Try omitting the obsdim"))
            O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
            @dimcheck size(target) == size(output)
            @dimcheck length(buffer) == size(output, O)
            fill!(buffer, zero(B))
            P = promote_type(Q,T)
            k = P(prod(size(output,n) for n in 1:N if n != O))
            nrm = 1 / k
            @inbounds @simd for I in CartesianIndices(size(output))
                buffer[I[O]] += ($FUN)(loss, target[I], output[I]) * nrm
            end
            buffer
        end

        # Compute the sum per observation (returns a Vector)
        function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                avg::AggMode.Sum,
                obsdim::ObsDim.Constant{O}) where {Q,T,N,O}
            S = typeof(($FUN)(loss, one(Q), one(T)))
            buffer = zeros(S, size(output, O))
            ($(Symbol(FUN,:!)))(buffer, loss, target, output, avg, obsdim)
        end

        # (mutating) Compute the sum per observation (returns a Vector)
        function ($(Symbol(FUN,:!)))(
                buffer::AbstractVector{B},
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                ::AggMode.Sum,
                ::ObsDim.Constant{O}) where {B,Q,T,N,O}
            N == 1 && throw(ArgumentError("Sum per observation non sensible for two Vectors. Try omitting the obsdim"))
            O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
            @dimcheck size(target) == size(output)
            @dimcheck length(buffer) == size(output, O)
            fill!(buffer, zero(B))
            @inbounds @simd for I in CartesianIndices(size(output))
                buffer[I[O]] += ($FUN)(loss, target[I], output[I])
            end
            buffer
        end
    end

    for KIND in (:MarginLoss, :DistanceLoss)
        @eval begin

            # --------------------------------------------------
            # FALLBACK

            # By default compute the element-wise result
            @inline function ($FUN)(loss::$KIND, numbers::AbstractArray)
                ($FUN)(loss, numbers, AggMode.None())
            end

            # (mutating) By default compute the element-wise result
            @inline function ($(Symbol(FUN,:!)))(
                    buffer::AbstractArray,
                    loss::$KIND,
                    numbers::AbstractArray)
                ($(Symbol(FUN,:!)))(buffer, loss, numbers, AggMode.None())
            end

            # Translate ObsDim.Last to the correct ObsDim.Constant (for code reduction)
            @inline function ($FUN)(
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    avg::AggregateMode,
                    ::ObsDim.Last = ObsDim.Last()) where {T,N}
                ($FUN)(loss, numbers, avg, ObsDim.Constant{N}())
            end

            # (mutating) Translate ObsDim.Last to the correct ObsDim.Constant (for code reduction)
            @inline function ($(Symbol(FUN,:!)))(
                    buffer::AbstractArray,
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    avg::AggregateMode,
                    ::ObsDim.Last = ObsDim.Last()) where {T,N}
                ($(Symbol(FUN,:!)))(buffer, loss, numbers, avg, ObsDim.Constant{N}())
            end

            # --------------------------------------------------
            # ELEMENT-WISE BROADCAST

            # Compute element-wise (returns an array)
            @inline function ($FUN)(
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    ::AggMode.None) where {T,N}
                S = typeof(($FUN)(loss, one(T)))
                ($FUN).(Ref(loss), numbers)::Array{S,N}
            end

            # (mutating) Compute element-wise (returns an array)
            function ($(Symbol(FUN,:!)))(
                    buffer::AbstractArray,
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    ::AggMode.None) where {T,N}
                buffer .= ($FUN).(Ref(loss), numbers)
                buffer
            end

            # --------------------------------------------------
            # REDUCE TO NUMBER

            # Compute the mean (returns a Number)
            function ($FUN)(
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    ::AggMode.Mean) where {T,N}
                nrm = 1 / length(numbers)
                mapreduce(x -> ($FUN)(loss, x) * nrm, +, numbers)
            end

            # Compute the sum (returns a Number)
            function ($FUN)(
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    ::AggMode.Sum) where {T,N}
                sum(x -> ($FUN)(loss, x), numbers)
            end

            # Compute the total weighted mean (returns a Number)
            function ($FUN)(
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    avg::AggMode.WeightedMean,
                    ::ObsDim.Constant{O}) where {T,N,O}
                @dimcheck size(numbers, O) == length(avg.weights)
                k = prod(n != O ? size(numbers,n) : 1 for n in 1:N)::Int
                nrm = avg.normalize ? inv(k * sum(avg.weights)) : inv(k * one(sum(avg.weights)))
                out = zero(($FUN)(loss, one(T)) * (avg.weights[1] * nrm))
                @inbounds @simd for I in CartesianIndices(size(numbers))
                    out += ($FUN)(loss, numbers[I]) * (avg.weights[I[O]] * nrm)
                end
                out
            end

            # Compute the total weighted sum (returns a Number)
            function ($FUN)(
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    avg::AggMode.WeightedSum,
                    ::ObsDim.Constant{O}) where {T,N,O}
                @dimcheck size(numbers, O) == length(avg.weights)
                nrm = avg.normalize ? inv(sum(avg.weights)) : inv(one(sum(avg.weights)))
                out = zero(($FUN)(loss, one(T)) * (avg.weights[1] * nrm))
                @inbounds @simd for I in CartesianIndices(size(numbers))
                    out += ($FUN)(loss, numbers[I]) * (avg.weights[I[O]] * nrm)
                end
                out
            end

            # --------------------------------------------------
            # PER OBSERVATION - SAME SIZE - TO VECTOR

            # Compute the mean per observation (returns a Vector)
            function ($FUN)(
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    avg::AggMode.Mean,
                    obsdim::ObsDim.Constant{O}) where {T,N,O}
                S = typeof(($FUN)(loss, one(T)) / one(Int))
                buffer = zeros(S, size(numbers, O))
                ($(Symbol(FUN,:!)))(buffer, loss, numbers, avg, obsdim)
            end

            # (mutating) Compute the mean per observation (returns a Vector)
            function ($(Symbol(FUN,:!)))(
                    buffer::AbstractVector{B},
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    ::AggMode.Mean,
                    ::ObsDim.Constant{O}) where {B,T,N,O}
                N == 1 && throw(ArgumentError("Mean per observation non sensible for two Vectors. Try omitting the obsdim"))
                O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
                @dimcheck length(buffer) == size(numbers, O)
                fill!(buffer, zero(B))
                k = prod(size(numbers,n) for n in 1:N if n != O)::Int
                nrm = 1 / k
                @inbounds @simd for I in CartesianIndices(size(numbers))
                    buffer[I[O]] += ($FUN)(loss, numbers[I]) * nrm
                end
                buffer
            end

            # Compute the sum per observation (returns a Vector)
            function ($FUN)(
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    avg::AggMode.Sum,
                    obsdim::ObsDim.Constant{O}) where {T,N,O}
                S = typeof(($FUN)(loss, one(T)))
                buffer = zeros(S, size(numbers, O))
                ($(Symbol(FUN,:!)))(buffer, loss, numbers, avg, obsdim)
            end

            # (mutating) Compute the sum per observation (returns a Vector)
            function ($(Symbol(FUN,:!)))(
                    buffer::AbstractVector{B},
                    loss::$KIND,
                    numbers::AbstractArray{T,N},
                    ::AggMode.Sum,
                    ::ObsDim.Constant{O}) where {B,T,N,O}
                N == 1 && throw(ArgumentError("Sum per observation non sensible for two Vectors. Try omitting the obsdim"))
                O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
                @dimcheck length(buffer) == size(numbers, O)
                fill!(buffer, zero(B))
                @inbounds @simd for I in CartesianIndices(size(numbers))
                    buffer[I[O]] += ($FUN)(loss, numbers[I])
                end
                buffer
            end
        end
    end
end

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

# abstract MarginLoss <: SupervisedLoss

value(loss::MarginLoss, target::Number, output::Number) = value(loss, target * output)
deriv(loss::MarginLoss, target::Number, output::Number) = target * deriv(loss, target * output)
deriv2(loss::MarginLoss, target::Number, output::Number) = deriv2(loss, target * output)
function value_deriv(loss::MarginLoss, target::Number, output::Number)
    v, d = value_deriv(loss, target * output)
    (v, target*d)
end

# Fallback for losses that don't want to take advantage of this
value_deriv(loss::MarginLoss, agreement::Number) = (value(loss, agreement), deriv(loss, agreement))

isunivfishercons(::MarginLoss) = false
isfishercons(loss::MarginLoss) = isunivfishercons(loss)
isnemitski(::MarginLoss) = true
ismarginbased(::MarginLoss) = true

# For every convex margin based loss L(a) the following statements
# are equivalent:
#   - L is classification calibrated
#   - L is differentiable at 0 and L'(0) < 0
# We use this here as fallback implementation.
# The other direction of the implication is not implemented however.
isclasscalibrated(loss::MarginLoss) = isconvex(loss) && isdifferentiable(loss, 0) && deriv(loss, 0) < 0

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

# abstract DistanceLoss <: SupervisedLoss

value(loss::DistanceLoss, target::Number, output::Number) = value(loss, output - target)
deriv(loss::DistanceLoss, target::Number, output::Number) = deriv(loss, output - target)
deriv2(loss::DistanceLoss, target::Number, output::Number) = deriv2(loss, output - target)
value_deriv(loss::DistanceLoss, target::Number, output::Number) = value_deriv(loss, output - target)

# Fallback for losses that don't want to take advantage of this
value_deriv(loss::DistanceLoss, difference::Number) = (value(loss, difference), deriv(loss, difference))

isdistancebased(::DistanceLoss) = true
issymmetric(::DistanceLoss) = false
isclipable(::DistanceLoss) = true

# --------------------------------------------------------------
# Fallback implementations and docstrings

@doc doc"""
    value(loss, target::Number, output::Number) -> Number

Compute the (non-negative) numeric result for the loss-function
denoted by the parameter `loss` and return it. Note that `target`
and `output` can be of different numeric type, in which case
promotion is performed in the manner appropriate for the given
loss.

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

```math
L : Y \times \mathbb{R} \rightarrow [0,\infty)
```

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we want to
  compute the value with.

- `target::Number`: The ground truth ``y \in Y`` of the observation.

- `output::Number`: The predicted output ``\hat{y} \in \mathbb{R}``
  for the observation.

# Examples

```jldoctest
#               loss        y    ŷ
julia> value(L1DistLoss(), 1.0, 2.0)
1.0

julia> value(L1DistLoss(), 1, 2)
1

julia> value(L1HingeLoss(), -1, 2)
3

julia> value(L1HingeLoss(), -1f0, 2f0)
3.0f0
```
"""
value(l::SupervisedLoss, target::Number, output::Number) =
    MethodError(value, (l, target, output))

"""
    value(loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode) -> Number

Compute the weighted or unweighted sum or mean (depending on
`avgmode`) of the individual values of the loss function for each
pair in `targets` and `outputs`. This method will not allocate a
temporary array.

In the case that the two parameters are arrays with a different
number of dimensions, broadcast will be performed. Note that the
given parameters are expected to have the same size in the
dimensions they share.

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

- `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

- `avgmode::AggregateMode`: Must be one of the following: [`AggMode.Sum()`](@ref),
  [`AggMode.Mean()`](@ref), [`AggMode.WeightedSum`](@ref), or
  [`AggMode.WeightedMean`](@ref).

# Examples

```jldoctest
julia> value(L1DistLoss(), [1,2,3], [2,5,-2], AggMode.Sum())
9

julia> value(L1DistLoss(), [1.,2,3], [2,5,-2], AggMode.Sum())
9.0

julia> value(L1DistLoss(), [1,2,3], [2,5,-2], AggMode.Mean())
3.0

julia> value(L1DistLoss(), Float32[1,2,3], Float32[2,5,-2], AggMode.Mean())
3.0f0
```
"""
value(l::Loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode) =
    MethodError(value, (l, target, output, avgmode))

"""
    value(loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::ObsDimension) -> AbstractVector

Compute the values of the loss function for each pair in
`targets` and `outputs` individually, and return either the
weighted or unweighted sum or mean for each observation (depending on
`avgmode`). This method will not allocate a temporary array, but
it will allocate the resulting vector.

Both arrays have to be of the same shape and size. Furthermore
they have to have at least two array dimensions (i.e. they must
not be vectors).

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

- `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

- `avgmode::AggregateMode`: Must be one of the following: [`AggMode.Sum()`](@ref),
  [`AggMode.Mean()`](@ref), [`AggMode.WeightedSum`](@ref), or
  [`AggMode.WeightedMean`](@ref).

- `obsdim::ObsDimension`: Specifies which of the array dimensions
  denotes the observations. see `?ObsDim` for more information.
"""
value(l::Loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::LearnBase.ObsDimension) =
    MethodError(value, (l, target, output, avgmode, obsdim))

"""
    value!(buffer::AbstractArray, loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::ObsDimension) -> buffer

Compute the values of the loss function for each pair in
`targets` and `outputs` individually, and return either the
weighted or unweighted sum or mean for each observation, depending on
`avgmode`. The results are stored into the given vector `buffer`.
This method will not allocate a temporary array.

Both arrays have to be of the same shape and size. Furthermore
they have to have at least two array dimensions (i.e. so they
must not be vectors).

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `buffer::AbstractArray`: Array to store the computed values in.
  Old values will be overwritten and lost.

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

- `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

- `avgmode::AggregateMode`: Must be one of the following: [`AggMode.Sum()`](@ref),
  [`AggMode.Mean()`](@ref), [`AggMode.WeightedSum`](@ref), or
  [`AggMode.WeightedMean`](@ref).

- `obsdim::ObsDimension`: Specifies which of the array dimensions
  denotes the observations. see `?ObsDim` for more information.

# Examples

```jldoctest
julia> targets = reshape(1:8, (2, 4)) ./ 8;

julia> outputs = reshape(1:2:16, (2, 4)) ./ 8;

julia> buffer = zeros(2);

julia> value!(buffer, L1DistLoss(), targets, outputs, AggMode.Sum(), ObsDim.First())
2-element Array{Float64,1}:
 1.5
 2.0

julia> buffer = zeros(4);

julia> value!(buffer, L1DistLoss(), targets, outputs, AggMode.Sum(), ObsDim.Last())
4-element Array{Float64,1}:
 0.125
 0.625
 1.125
 1.625
```
"""
value!(buffer::AbstractArray, l::Loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::LearnBase.ObsDimension) =
    MethodError(value!, (buffer, l, target, output, avgmode, obsdim))

@doc doc"""
    deriv(loss, target::Number, output::Number) -> Number

Compute the derivative for the loss-function (denoted by the
parameter `loss`) in respect to the `output`. Note that
`target` and `output` can be of different numeric type,
in which case promotion is performed in the manner
appropriate for the given loss.

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we want to
  compute the derivative with.

- `target::Number`: The ground truth ``y \in Y`` of the observation.

- `output::Number`: The predicted output ``\hat{y} \in \mathbb{R}``
  for the observation.

# Examples

```jldoctest
#               loss        y    ŷ
julia> deriv(L2DistLoss(), 1.0, 2.0)
2.0

julia> deriv(L2DistLoss(), 1, 2)
2

julia> deriv(L2HingeLoss(), -1, 2)
6

julia> deriv(L2HingeLoss(), -1f0, 2f0)
6.0f0
```
"""
deriv(l::SupervisedLoss, target::Number, output::Number) =
    MethodError(deriv, (l, target, output))

"""
    deriv(loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode) -> Number

Compute the weighted or unweighted sum or mean (depending on `avgmode`)
of the individual derivatives of the loss function for each pair
in `targets` and `outputs`. This method will not allocate a
temporary array.

In the case that the two parameters are arrays with a different
number of dimensions, broadcast will be performed. Note that the
given parameters are expected to have the same size in the
dimensions they share.

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

- `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

- `avgmode::AggregateMode`: Must be one of the following: [`AggMode.Sum()`](@ref),
  [`AggMode.Mean()`](@ref), [`AggMode.WeightedSum`](@ref), or
  [`AggMode.WeightedMean`](@ref).

# Examples

```jldoctest
julia> deriv(L2DistLoss(), [1,2,3], [2,5,-2], AggMode.Sum())
-2

julia> deriv(L2DistLoss(), [1,2,3], [2,5,-2], AggMode.Mean())
-0.6666666666666666
```
"""
deriv(l::Loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode) =
    MethodError(deriv, (l, target, output, avgmode))

"""
    deriv(loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::ObsDimension) -> AbstractVector

Compute the derivative of the loss function for each pair in
`targets` and `outputs` individually, and return either the
weighted or unweighted sum or mean for each observation (depending on
`avgmode`). This method will not allocate a temporary array, but
it will allocate the resulting vector.

Both arrays have to be of the same shape and size. Furthermore
they have to have at least two array dimensions (i.e. they must
not be vectors).

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

- `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

- `avgmode::AggregateMode`: Must be one of the following: [`AggMode.Sum()`](@ref),
  [`AggMode.Mean()`](@ref), [`AggMode.WeightedSum`](@ref), or
  [`AggMode.WeightedMean`](@ref).

- `obsdim::ObsDimension`: Specifies which of the array dimensions
  denotes the observations. see `?ObsDim` for more information.
"""
deriv(l::Loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::LearnBase.ObsDimension) =
    MethodError(deriv, (l, target, output, avgmode, obsdim))

"""
    deriv!(buffer::AbstractArray, loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::ObsDimension) -> buffer

Compute the derivative of the loss function for each pair in
`targets` and `outputs` individually, and return either the
weighted or unweighted sum or mean for each observation, depending on
`avgmode`. The results are stored into the given vector `buffer`.
This method will not allocate a temporary array.

Both arrays have to be of the same shape and size. Furthermore
they have to have at least two array dimensions (i.e. so they
must not be vectors).

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `buffer::AbstractArray`: Array to store the computed values in.
  Old values will be overwritten and lost.

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

- `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

- `avgmode::AggregateMode`: Must be one of the following: [`AggMode.Sum()`](@ref),
  [`AggMode.Mean()`](@ref), [`AggMode.WeightedSum`](@ref), or
  [`AggMode.WeightedMean`](@ref).

- `obsdim::ObsDimension`: Specifies which of the array dimensions
  denotes the observations. see `?ObsDim` for more information.

# Examples

```jldoctest
julia> targets = reshape(1:8, (2, 4)) ./ 8;

julia> outputs = reshape(1:2:16, (2, 4)) ./ 8;

julia> buffer = zeros(2);

julia> deriv!(buffer, L1DistLoss(), targets, outputs, AggMode.Sum(), ObsDim.First())
2-element Array{Float64,1}:
 3.0
 4.0

julia> buffer = zeros(4);

julia> deriv!(buffer, L1DistLoss(), targets, outputs, AggMode.Sum(), ObsDim.Last())
4-element Array{Float64,1}:
 1.0
 2.0
 2.0
 2.0
```
"""
deriv!(buffer::AbstractArray, l::Loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::LearnBase.ObsDimension) =
    MethodError(deriv!, (buffer, l, target, output, avgmode, obsdim))

@doc doc"""
    deriv2(loss, target::Number, output::Number) -> Number

Compute the second derivative for the loss-function (denoted by
the parameter `loss`) in respect to the `output`. Note that
`target` and `output` can be of different numeric type, in which
case promotion is performed in the manner appropriate for the
given loss.

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we want to
  compute the second derivative with.

- `target::Number`: The ground truth ``y \in Y`` of the observation.

- `output::Number`: The predicted output ``\hat{y} \in \mathbb{R}``
  for the observation.

# Examples

```jldoctest
#               loss             y    ŷ
julia> deriv2(LogitDistLoss(), -0.5, 0.3)
0.42781939304058886

julia> deriv2(LogitMarginLoss(), -1f0, 2f0)
0.104993574f0
```
"""
deriv2(l::SupervisedLoss, target::Number, output::Number) =
    MethodError(deriv2, (l, target, output))

"""
    deriv2(loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode) -> Number

Compute the weighted or unweighted sum or mean (depending on `avgmode`)
of the individual second derivatives of the loss function for
each pair in `targets` and `outputs`. This method will not
allocate a temporary array.

In the case that the two parameters are arrays with a different
number of dimensions, broadcast will be performed. Note that the
given parameters are expected to have the same size in the
dimensions they share.

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

- `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

- `avgmode::AggregateMode`: Must be one of the following: [`AggMode.Sum()`](@ref),
  [`AggMode.Mean()`](@ref), [`AggMode.WeightedSum`](@ref), or
  [`AggMode.WeightedMean`](@ref).

# Examples

```jldoctest
julia> deriv2(LogitDistLoss(), [1.,2,3], [2,5,-2], AggMode.Sum())
0.49687329928636825

julia> deriv2(LogitDistLoss(), [1.,2,3], [2,5,-2], AggMode.Mean())
0.1656244330954561
```
"""
deriv2(l::Loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode) =
    MethodError(deriv2, (l, target, output, avgmode))

"""
    deriv2(loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::ObsDimension) -> AbstractVector

Compute the second derivative of the loss function for each pair in
`targets` and `outputs` individually, and return either the
weighted or unweighted sum or mean for each observation (depending on
`avgmode`). This method will not allocate a temporary array, but
it will allocate the resulting vector.

Both arrays have to be of the same shape and size. Furthermore
they have to have at least two array dimensions (i.e. they must
not be vectors).

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

- `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

- `avgmode::AggregateMode`: Must be one of the following: [`AggMode.Sum()`](@ref),
  [`AggMode.Mean()`](@ref), [`AggMode.WeightedSum`](@ref), or
  [`AggMode.WeightedMean`](@ref).

- `obsdim::ObsDimension`: Specifies which of the array dimensions
  denotes the observations. see `?ObsDim` for more information.
"""
deriv2(l::Loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::LearnBase.ObsDimension) =
    MethodError(deriv2, (l, target, output, avgmode, obsdim))

"""
    deriv2!(buffer::AbstractArray, loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::ObsDimension) -> buffer

Compute the second derivative of the loss function for each pair in
`targets` and `outputs` individually, and return either the
weighted or unweighted sum or mean for each observation, depending on
`avgmode`. The results are stored into the given vector `buffer`.
This method will not allocate a temporary array.

Both arrays have to be of the same shape and size. Furthermore
they have to have at least two array dimensions (i.e. so they
must not be vectors).

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `buffer::AbstractArray`: Array to store the computed values in.
  Old values will be overwritten and lost.

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `targets::AbstractArray`: The array of ground truths ``\\mathbf{y}``.

- `outputs::AbstractArray`: The array of predicted outputs ``\\mathbf{\\hat{y}}``.

- `avgmode::AggregateMode`: Must be one of the following: [`AggMode.Sum()`](@ref),
  [`AggMode.Mean()`](@ref), [`AggMode.WeightedSum`](@ref), or
  [`AggMode.WeightedMean`](@ref).

- `obsdim::ObsDimension`: Specifies which of the array dimensions
  denotes the observations. see `?ObsDim` for more information.

# Examples

```jldoctest
julia> targets = reshape(1:8, (2, 4)) ./ 8;

julia> outputs = reshape(1:2:16, (2, 4)) ./ 8;

julia> buffer = zeros(2);

julia> deriv2!(buffer, L2DistLoss(), targets, outputs, AggMode.Sum(), ObsDim.First())
2-element Array{Float64,1}:
 8.0
 8.0

julia> buffer = zeros(4);

julia> deriv2!(buffer, L2DistLoss(), targets, outputs, AggMode.Sum(), ObsDim.Last())
4-element Array{Float64,1}:
 4.0
 4.0
 4.0
 4.0
```
"""
deriv2!(buffer::AbstractArray, l::Loss, target::AbstractArray, output::AbstractArray, avgmode::AggregateMode, obsdim::LearnBase.ObsDimension) =
    MethodError(deriv2!, (buffer, l, target, output, avgmode, obsdim))

@doc doc"""
    value_deriv(loss, target::Number, output::Number) -> Tuple

Return the results of [`value`](@ref) and [`deriv`](@ref) as a
tuple, in which the first element is the value and the second
element the derivative.

In some cases this function can yield better performance, because
the losses can make use of shared variables when computing
the results. Note that `target` and `output` can be of
different numeric type, in which case promotion is performed
in the manner appropriate for the given loss.

Note: This function should always be type-stable. If it isn't,
you likely found a bug.

# Arguments

- `loss::SupervisedLoss`: The loss-function ``L`` we are working with.

- `target::Number`: The ground truth ``y \in Y`` of the observation.

- `output::Number`: The predicted output ``\hat{y} \in \mathbb{R}``

# Examples

```jldoctest
#                     loss         y    ŷ
julia> value_deriv(L2DistLoss(), -1.0, 3.0)
(16.0, 8.0)
```
"""
function value_deriv(l::SupervisedLoss, target::Number, output::Number)
    value(l, target, output), deriv(l, target, output)
end

@doc doc"""
    isconvex(loss::SupervisedLoss) -> Bool

Return `true` if the given `loss` denotes a convex function.
A function ``f : \mathbb{R}^n \rightarrow \mathbb{R}`` is convex if
its domain is a convex set and if for all ``x, y`` in that
domain, with ``\theta`` such that for ``0 \leq \theta \leq 1``,
we have

```math
f(\theta x + (1 - \theta) y) \leq \theta f(x) + (1 - \theta) f(y)
```

# Examples

```jldoctest
julia> isconvex(LPDistLoss(0.5))
false

julia> isconvex(ZeroOneLoss())
false

julia> isconvex(L1DistLoss())
true

julia> isconvex(L2DistLoss())
true
```
"""
isconvex

@doc doc"""
    isstrictlyconvex(loss::SupervisedLoss) -> Bool

Return `true` if the given `loss` denotes a strictly convex function.
A function ``f : \mathbb{R}^n \rightarrow \mathbb{R}`` is
strictly convex if its domain is a convex set and if for all
``x, y`` in that domain where ``x \neq y``, with
``\theta`` such that for ``0 < \theta < 1``, we have

```math
f(\theta x + (1 - \theta) y) < \theta f(x) + (1 - \theta) f(y)
```

# Examples

```jldoctest
julia> isstrictlyconvex(L1DistLoss())
false

julia> isstrictlyconvex(LogitDistLoss())
true

julia> isstrictlyconvex(L2DistLoss())
true
```
"""
isstrictlyconvex

@doc doc"""
    isstronglyconvex(loss::SupervisedLoss) -> Bool

Return `true` if the given `loss` denotes a strongly convex function.
A function ``f : \mathbb{R}^n \rightarrow \mathbb{R}`` is
``m``-strongly convex if its domain is a convex set and if
``\forall x,y \in`` **dom** ``f`` where ``x \neq y``,
and ``\theta`` such that for ``0 \le \theta \le 1`` , we have

```math
f(\theta x + (1 - \theta)y) < \theta f(x) + (1 - \theta) f(y) - 0.5 m \cdot \theta (1 - \theta) {\| x - y \|}_2^2
```

In a more familiar setting, if the loss function is
differentiable we have

```math
\left( \nabla f(x) - \nabla f(y) \right)^\top (x - y) \ge m {\| x - y\|}_2^2
```

# Examples

```jldoctest
julia> isstronglyconvex(L1DistLoss())
false

julia> isstronglyconvex(LogitDistLoss())
false

julia> isstronglyconvex(L2DistLoss())
true
```
"""
isstronglyconvex(::SupervisedLoss) = false

@doc doc"""
    isdifferentiable(loss::SupervisedLoss, [x::Number]) -> Bool

Return `true` if the given `loss` is differentiable
(optionally limited to the given point `x` if specified).

A function ``f : \mathbb{R}^{n} \rightarrow \mathbb{R}^{m}`` is
differentiable at a point ``x \in`` **int dom** ``f``, if there
exists a matrix ``Df(x) \in \mathbb{R}^{m \times n}`` such that
it satisfies:

```math
\lim_{z \neq x, z \to x} \frac{{\|f(z) - f(x) - Df(x)(z-x)\|}_2}{{\|z - x\|}_2} = 0
```

A function is differentiable if its domain is open and it is
differentiable at every point ``x``.

# Examples

```jldoctest
julia> isdifferentiable(L1DistLoss())
false

julia> isdifferentiable(L1DistLoss(), 1)
true

julia> isdifferentiable(L2DistLoss())
true
```
"""
isdifferentiable(l::SupervisedLoss) = istwicedifferentiable(l)
isdifferentiable(l::SupervisedLoss, at) = isdifferentiable(l)

@doc doc"""
    istwicedifferentiable(loss::SupervisedLoss, [x::Number]) -> Bool

Return `true` if the given `loss` is differentiable
(optionally limited to the given point `x` if specified).

A function ``f : \mathbb{R}^{n} \rightarrow \mathbb{R}`` is
said to be twice differentiable at a point ``x \in`` **int
dom** ``f``, if the function derivative for ``\nabla f``
exists at ``x``.

```math
\nabla^2 f(x) = D \nabla f(x)
```

A function is twice differentiable if its domain is open and it
is twice differentiable at every point ``x``.

# Examples

```jldoctest
julia> isdifferentiable(L1DistLoss())
false

julia> isdifferentiable(L1DistLoss(), 1)
true

julia> isdifferentiable(L2DistLoss())
true
```
"""
istwicedifferentiable(::SupervisedLoss) = false
istwicedifferentiable(l::SupervisedLoss, at) = istwicedifferentiable(l)

@doc doc"""
    islocallylipschitzcont(loss::SupervisedLoss) -> Bool

Return `true` if the given `loss` function is locally-Lipschitz
continous.

A supervised loss ``L : Y \times \mathbb{R} \rightarrow [0, \infty)``
is called locally Lipschitz continuous if ``\forall a \ge 0``
there exists a constant :math:`c_a \ge 0`, such that

```math
\sup_{y \in Y} \left| L(y,t) − L(y,t′) \right| \le c_a |t − t′|,  \qquad  t,t′ \in [−a,a]
```

Every convex function is locally lipschitz continuous

# Examples

```jldoctest
julia> islocallylipschitzcont(ExpLoss())
true

julia> islocallylipschitzcont(SigmoidLoss())
true
```
"""
islocallylipschitzcont(l::SupervisedLoss) = isconvex(l) || islipschitzcont(l)

@doc doc"""
    islipschitzcont(loss::SupervisedLoss) -> Bool

Return `true` if the given `loss` function is Lipschitz continuous.

A supervised loss function ``L : Y \times \mathbb{R} \rightarrow
[0, \infty)`` is Lipschitz continous, if there exists a finite
constant ``M < \infty`` such that

```math
|L(y, t) - L(y, t′)| \le M |t - t′|,  \qquad  \forall (y, t) \in Y \times \mathbb{R}
```

# Examples

```jldoctest
julia> islipschitzcont(SigmoidLoss())
true

julia> islipschitzcont(ExpLoss())
false
```
"""
islipschitzcont(::SupervisedLoss) = false

@doc doc"""
    isnemitski(loss::SupervisedLoss) -> Bool

Return `true` if the given `loss` denotes a Nemitski loss function.

We call a supervised loss function ``L : Y \times \mathbb{R}
\rightarrow [0,\infty)`` a Nemitski loss if there exist a
measurable function ``b : Y \rightarrow [0, \infty)`` and an
increasing function ``h : [0, \infty) \rightarrow [0, \infty)``
such that

```math
L(y,\hat{y}) \le b(y) + h(|\hat{y}|),  \qquad  (y, \hat{y}) \in Y \times \mathbb{R}.
```

If a loss if locally lipsschitz continuous then it is a Nemitski loss
"""
isnemitski(l::SupervisedLoss) = islocallylipschitzcont(l)

@doc doc"""
    isclipable(loss::SupervisedLoss) -> Bool

Return `true` if the given `loss` function is clipable. A
supervised loss ``L : Y \times \mathbb{R} \rightarrow [0,
\infty)`` can be clipped at ``M > 0`` if, for all ``(y,t)
\in Y \times \mathbb{R}``,

```math
L(y, \hat{t}) \le L(y, t)
```

where ``\hat{t}`` denotes the clipped value of ``t`` at
``\pm M``. That is

```math
\hat{t} = \begin{cases} -M & \quad \text{if } t < -M \\ t & \quad \text{if } t \in [-M, M] \\ M & \quad \text{if } t > M \end{cases}
```

# Examples

```jldoctest
julia> isclipable(ExpLoss())
false

julia> isclipable(L2DistLoss())
true
```
"""
isclipable(::SupervisedLoss) = false

@doc doc"""
    ismarginbased(loss::SupervisedLoss) -> Bool

Return `true` if the given `loss` is a margin-based loss.

A supervised loss function ``L : Y \times \mathbb{R} \rightarrow
[0, \infty)`` is said to be **margin-based**, if there exists a
representing function ``\psi : \mathbb{R} \rightarrow [0, \infty)``
satisfying

```math
L(y, \hat{y}) = \psi (y \cdot \hat{y}),  \qquad  (y, \hat{y}) \in Y \times \mathbb{R}
```

# Examples

```jldoctest
julia> ismarginbased(HuberLoss(2))
false

julia> ismarginbased(L2MarginLoss())
true
```
"""
ismarginbased(::SupervisedLoss) = false

"""
    isclasscalibrated(loss::SupervisedLoss) -> Bool
"""
isclasscalibrated(::SupervisedLoss) = false

@doc doc"""
    isdistancebased(loss::SupervisedLoss) -> Bool

Return `true` ifthe given `loss` is a distance-based loss.

A supervised loss function ``L : Y \times \mathbb{R} \rightarrow
[0, \infty)`` is said to be **distance-based**, if there exists a
representing function ``\psi : \mathbb{R} \rightarrow [0, \infty)``
satisfying ``\psi (0) = 0`` and

```math
L(y, \hat{y}) = \psi (\hat{y} - y),  \qquad  (y, \hat{y}) \in Y \times \mathbb{R}
```

# Examples

```jldoctest
julia> isdistancebased(HuberLoss(2))
true

julia> isdistancebased(L2MarginLoss())
false
```
"""
isdistancebased(::SupervisedLoss) = false

@doc doc"""
    issymmetric(loss::SupervisedLoss) -> Bool

Return `true` if the given loss is a symmetric loss.

A function ``f : \mathbb{R} \rightarrow [0,\infty)`` is said
to be symmetric about origin if we have

```math
f(x) = f(-x), \qquad  \forall x \in \mathbb{R}
```

A distance-based loss is said to be symmetric if its representing
function is symmetric.

# Examples

```jldoctest
julia> issymmetric(QuantileLoss(0.2))
false

julia> issymmetric(LPDistLoss(2))
true
```
"""
issymmetric(::SupervisedLoss) = false

isminimizable(l::SupervisedLoss) = isconvex(l)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "HOHQMesh"
version = v"1.3.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/trixi-framework/HOHQMesh/releases/download/v$(version)/HOHQMesh-v$(version).tar.gz", "da5c1420c7ba312769d0abffbec6b1b7271dcf9a8f7a3f94465b2c1b1d17c418"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/HOHQMesh-*

# Build HOHQMesh
make -j${nproc} F90=$FC

# Copy executable to binary directory
mkdir -p ${bindir}
if [[ "${bb_full_target}" == *-w64-mingw32-libgfortran[34]* ]]; then
  # For some reason, on i686 with libgfortran3 or libgfortran4,
  # the executable is created without the `.exe` extension
  cp HOHQMesh ${bindir}/HOHQMesh${exeext}
else
  cp HOHQMesh${exeext} ${bindir}
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc = "glibc"),
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("aarch64", "linux"; libc = "glibc"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "glibc"),
    Platform("powerpc64le", "linux"; libc = "glibc"),
    Platform("i686", "linux"; libc = "musl"),
    Platform("x86_64", "linux"; libc = "musl"),
    Platform("x86_64", "macos"; ),
    Platform("x86_64", "freebsd"; ),
    Platform("i686", "windows"; ),
    Platform("x86_64", "windows"; ),
    Platform("aarch64", "macos"),
    Platform("armv6l", "linux"; call_abi = "eabihf", libc = "glibc"),
]

platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("HOHQMesh", :HOHQMesh)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9.1.0", julia_compat="1.6")

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
using PkgBenchmark

function ask_version()
  print("What version would you like to compare to? ")
  String(chomp(strip(readline(STDIN))))
end

results = judge("ImageQuilting", ask_version())

showall(PkgBenchmark.benchmarkgroup(results))

println()


text:
using Random
function sgd!(net::Network,train::Vector{Tuple{Vector{Float64},Vector{Float64}}})::Nothing
    for e=1:net.epochs
        train = shuffle(train)
        batches = [train[(1+(i*net.batch)):((i+1)*net.batch)] for i=0:(Int(length(train)/net.batch)-1)]

        for b in batches
            errors = []
            activations = []

            for x in b
                push!(errors,backprop(net,x[1],x[2])[1])
                push!(activations,backprop(net,x[1],x[2])[2])
            end

            # change w & b using errors = []
            for w1=1:length(net.weights)
                net.weights[wl] = net.weights[wl].-(net.η/net.batch).*sum([errors[e][wl]*transpose(activations[e][wl-1]) for e=1:length(errors)])

            end
            for bl=length(net.biases)
                net.biases[bl] = net.biases[bl].-(net.η/net.batch).*sum([errors[e][bl] for e=1:length(errors)])

            end
        end
    end

    return nothing
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
## Token interface to a container.  A token is the address
## of an item in a container.  The token has two parts: the
## container and the item's address.  The address is of type
## AbstractSemiToken.

module Tokens

import Base.isless
import Base.isequal
import Base.colon

abstract AbstractSemiToken

immutable IntSemiToken <: AbstractSemiToken
    address::Int
end

immutable Token{T, S <: AbstractSemiToken}
    container::T
    semitoken::S
end

## The following two operations extract the two parts of a token.
semi(i::Token) = i.semitoken
container(i::Token) = i.container

## This operation puts a token together from its parts.
assemble(m, s::AbstractSemiToken) = Token(m,s)

## The remaining functions operate on the container
## data via the token.  They are undefined (i.e., return
## nothing) at this abstract level.  It is up to
## specific implementations to define them.

## Return a key indexed by a token.
deref_key(i::Token) = nothing

## Return a value indexed by a token.
deref_value(i::Token) = nothing

## Return a key-value pair indexed by a token.
deref(i::Token) = nothing

## Return the status of a token (is it valid)
status(i::Token) = nothing

## Delete the item addressed by a token.
delete!(i::Token) = nothing

## Advance the token (return a new token) to the
## next item in the container.
advance(i::Token) = nothing

## Regress the token (return a new token) to the
## previous item in the container.
regress(i::Token) = nothing

## Compare tokens
isless(i::Token, j::Token) = nothing
isequal(i::Token, j::Token) = nothing

## Index a range of a container using a start/end
## token pair.  The colon operator is inclusive;
## the excludelast operation excludes the end token.
colon(i::Token, j::Token) = nothing
excludelast(i::Token, j::Token) = nothing

export semi, container, assemble, deref_key, deref_value
export deref, status, delete!, advance, regress
export excludelast

end

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
# Alternative Krylov phiv methods using error estimates of Saad to automatically
# terminate Arnoldi/Lanczos iterations.
# Currently only expv for Lanczos is implemented.

########################################
# Cache types
abstract type SubspaceCache{T} end
abstract type HermitianSubspaceCache{T} <: SubspaceCache{T} end

mutable struct StegrCache{T,R<:Real} <: HermitianSubspaceCache{T}
    v::Vector{T} # Subspace-propagated vector
    w::Vector{T}
    sw::Stegr.StegrWork{R}
    StegrCache(::Type{T}, n::Integer) where T = new{T,real(T)}(
        Vector{T}(undef, n), Vector{T}(undef, n),
        Stegr.StegrWork(real(T), n))
end

"""
    expT!(α, β, t, cache)

Calculate the subspace exponential `exp(t*T)` for a tridiagonal
subspace matrix `T` with `α` on the diagonal and `β` on the
super-/subdiagonal, diagonalizing via `stegr!`.
"""
function expT!(α::AbstractVector{R}, β::AbstractVector{R}, t::Number,
               cache::StegrCache{T,R}) where {T,R<:Real}
    LAPACK.stegr!(α, β, cache.sw)
    sel = 1:length(α)
    @inbounds for i = sel
        cache.w[i] = exp(t*cache.sw.w[i])*cache.sw.Z[1,i]
    end
    mul!(@view(cache.v[sel]), @view(cache.sw.Z[sel,sel]), @view(cache.w[sel]))
end

get_subspace_cache(Ks::KrylovSubspace{T,U}) where {T,U<:Complex} =
    error("Subspace exponential caches not yet available for non-Hermitian matrices.")
get_subspace_cache(Ks::KrylovSubspace{T,U}) where {T,U<:Real} =
    StegrCache(T, Ks.maxiter)

########################################
# Phiv with error estimate as termination condition
"""
    expv!(w, t, A, b, Ks, cache)

Alternative interface for calculating the action of `exp(t*A)` on the
vector `b`, storing the result in `w`. The Krylov iteration is
terminated when an error estimate for the matrix exponential in the
generated subspace is below the requested tolerance. `Ks` is a
`KrylovSubspace` and `typeof(cache)<:HermitianSubspaceCache`, the
exact type decides which algorithm is used to compute the subspace
exponential.
"""
function expv!(w::AbstractVector{T}, t::Number, A, b::AbstractVector{T},
               Ks::KrylovSubspace{T, B, B}, cache::HSC;
               atol::B=1.0e-8, rtol::B=1.0e-4,
               m=min(Ks.maxiter, size(A,1)),
               verbose::Bool=false) where {B, T <: Number, HSC <: HermitianSubspaceCache}
    if m > Ks.maxiter
        resize!(Ks, m)
    else
        Ks.m = m # might change if error estimate is below requested tolerance
    end

    V, H = getV(Ks), getH(Ks)
    Ks.beta = norm(b)
    @. V[:, 1] = b / Ks.beta

    ε = atol + rtol * Ks.beta
    verbose && @printf("Initial norm: β₀ %e, stopping threshold: %e\n", Ks.beta, ε)

    α = @diagview(H)
    β = @diagview(H,-1)
    n = size(V, 1)

    for j ∈ 1:m
        lanczos_step!(j, A, V, α, β)
        expT!(@view(α[1:j]), @view(β[1:j]), t, cache)

        # This is practical error estimate Er₂ from
        #
        #   Saad, Y. (1992). Analysis of some Krylov subspace
        #   approximations. SIAM Journal on Numerical Analysis.
        σ = β[j]*Ks.beta*abs(cache.v[j])
        verbose && @printf("iter %d, α[%d] %e, β[%d] %e, σ %e\n",j, j, α[j], j, β[j], σ)
        if σ < ε
            Ks.m = j
            break
        end
    end
    verbose && println("Krylov subspace size: ", Ks.m)

    lmul!(Ks.beta, mul!(w, @view(Ks.V[:,1:Ks.m]), @view(cache.v[1:Ks.m])))
end


text:
using Oceananigans.Models: AbstractModel
using Oceananigans.Utils: prettykeys

mutable struct NaNChecker{F}
    fields :: F
    erroring :: Bool
end

NaNChecker(fields) = NaNChecker(fields, false) # default

function Base.summary(nc::NaNChecker)
    fieldnames = prettykeys(nc.fields)
    if nc.erroring
        return "Erroring NaNChecker for $fieldnames"
    else
        return "NaNChecker for $fieldnames"
    end
end

Base.show(io, nc::NaNChecker) = print(io, summary(nc))

"""
    NaNChecker(; fields, erroring=false)

Return a `NaNChecker`, which sets `sim.running=false` if a `NaN` is detected
in any member of `fields` when `NaNChecker(sim)` is called. `fields` should be
a container with key-value pairs like a dictionary or `NamedTuple`.

If `erroring=true`, the `NaNChecker` will throw an error on NaN detection.
"""
NaNChecker(; fields, erroring=false) = NaNChecker(fields, erroring)

hasnan(field) = any(isnan, parent(field))
hasnan(model::AbstractModel) = hasnan(first(fields(model)))

function (nc::NaNChecker)(simulation)
    for (name, field) in pairs(nc.fields)
        if hasnan(field)
            simulation.running = false

            if nc.erroring
                clock = simulation.model.clock
                error("time = $(clock.time), iteration = $(clock.iteration): NaN found in field $name. Aborting simulation.")
            else
                @info "NaN found in field $name. Stopping simulation."
            end
        end
    end
    return nothing
end

"""
    erroring_NaNChecker!(simulation)

Toggles `simulation`'s `NaNChecker` to throw an error when a `NaN` is detected.
"""
function erroring_NaNChecker!(simulation)
    simulation.callbacks[:nan_checker].func.erroring = true
    return nothing
end


text:
zlsb250zhudong = PrimTransmission(18,73,"曲轴")
zlsb250zhuzhou = DrivingTransmission(11,14,16,18,20,26,"主轴")
zlsb250fuzhou = DrivenTransmission(32,25,22,21,18,21,"副轴")
zlsb250flywheel = Flywheel(14,41,"大小飞")

zlsb250chuan = Transmission(zlsb250zhudong,zlsb250zhuzhou,zlsb250fuzhou,zlsb250flywheel)

zlsb250 = vehiclespeed("ZLsb250",zlsb250chuan)

x_sl250 = 1:6
y_sl250 = zlsb250.second_ratio

max_stars_repo_name: powerzf/VehicleCalculate

max_issues_repo_name: powerzf/VehicleCalculate

max_forks_repo_name: powerzf/VehicleCalculate


text:
using Documenter
using fastPHASE

makedocs(
    sitename = "fastPHASE",
    format = Documenter.HTML(),
    modules = [fastPHASE],
    pages = [
        "Home" => "index.md",
        "API" => "man/api.md",
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo   = "github.com/biona001/fastPHASE.jl.git",
    target = "build"
)


text:
using DynamicPolynomials
using NamedArrays

struct gensys_sys
    Gamma0::AbstractArray
    Gamma1::AbstractArray
    Psi::AbstractArray
    Pi::AbstractArray
    backward::AbstractArray
    forward::AbstractArray
end

function find_expectations(variables)
    vars_string = [repr(vars) for vars in variables]

    #check which variables are expectations and get their length

    expec_vars_test = match.(r"_e",vars_string)
    expec_vars_test = expec_vars_test .!= nothing
    length_expec = sum(expec_vars_test)
    println("Found ",length_expec, " expectations variables")
    expec_vars = variables[expec_vars_test]
    return expec_vars
end

function find_lags(variables)
    vars_string = [repr(vars) for vars in variables]

    #check which variables are expectations and get their length

    lag_vars_test = match.(r"_l",vars_string)
    lag_vars_test = lag_vars_test .!= nothing
    length_lag = sum(lag_vars_test)
    println("Found ",length_lag, " lags variables")
    lag_vars = variables[lag_vars_test]
    return lag_vars
end

function get_current(variables)
    vars_string = [repr(vars) for vars in variables]

    #check which variables are expectations and get their length

    cur_vars_test = match.(r"_",vars_string)
    cur_vars_test = cur_vars_test .=== nothing
    length_cur = sum(cur_vars_test)
    println("Found ",length_cur, " current variables")
    cur_vars = variables[cur_vars_test]
    return cur_vars
end

function get_current_from_lag(variables,lags)
    lags_string = [repr(vars) for vars in lags]
    vars_string = [repr(vars) for vars in variables]

    stripped_vars = [rstrip(vv,['_','l']) for vv in lags_string]

    current_vars_index = zeros(Int64,length(lags))

    for j in 1:length(stripped_vars)
        test = stripped_vars[j] .== vars_string
        if sum(test) == 0
            @error "Lag but no current"
        else
            current_vars_index[j] = findfirst(test)
        end
    end
    return variables[current_vars_index]
end

function get_current_from_foward(variables,fowards)
    fowards_string = [repr(vars) for vars in fowards]
    vars_string = [repr(vars) for vars in variables]

    stripped_vars = [rstrip(vv,['_','f']) for vv in fowards_string]

    current_vars_index = zeros(Int64,length(fowards))

    for j in 1:length(stripped_vars)
        test = stripped_vars[j] .== vars_string
        if sum(test) == 0
            @error "Expectation but no current"
        else
            current_vars_index[j] = findfirst(test)
        end
    end
    return variables[current_vars_index]
end

function declare_aux_variables(expec_vars)
    expec_vars_string = [repr(vars) for vars in expec_vars]
     #transform variables in strings
    stripped_vars = [rstrip(vv,['_','e']) for vv in expec_vars_string]
    foward_vars = string.(stripped_vars,"_f")
    expectation_errors = string.("nu_",stripped_vars)
    expec_vars_aux = string.(expec_vars_string,"_aux")
    length_expec = length(expec_vars)

    all_vars = [reshape(foward_vars,1,length_expec) reshape(expectation_errors,1,length_expec)]

    def = "@polyvar"

    for vars in all_vars
        def = join([def,vars]," ")
    end
    eval(Meta.parse(def))

    def_aux = string.(expec_vars_aux,"=",foward_vars,"+",expectation_errors)
    eval.(Meta.parse.(def_aux))
    return foward_vars,expectation_errors,expec_vars_aux
end

function model2gensys(model,shocks)
    keys_model = collect(keys(model))
    varis = Dict()

    for i in keys(model)
        merge!(varis,Dict(i=>variables(model[i])))
    end

    vari_all = copy(varis[keys_model[1]])

    for i in keys(varis)
        union!(vari_all,variables(model[i]))
    end

    expecs = find_expectations(vari_all)
    lags = find_lags(vari_all)
    f_vars, expec_errors, expecs_aux = declare_aux_variables(expecs)

    expec_errors_names = expec_errors

    expecs_aux = eval.(Meta.parse.(expecs_aux))
    f_vars = eval.(Meta.parse.(f_vars))
    expec_errors = eval.(Meta.parse.(expec_errors))

    currents = get_current(vari_all)
    cur_w_lags = get_current_from_lag(vari_all,lags)
    cur_w_f = get_current_from_foward(vari_all,f_vars)
    only_cur = setdiff(currents,[cur_w_lags;cur_w_f;shocks])

    #print(only_cur)

    #SUbstituting the expectations variable by their counterpart true value + expectation error

    for j in 1:length(expecs)
        for i in keys(model)
            model[i] = subs(model[i],expecs[j] => expecs_aux[j])
        end
    end

    foward_variables = [f_vars;setdiff(currents,[cur_w_f; shocks; only_cur])]
    backward_variables = [setdiff(currents,[cur_w_lags; shocks;only_cur]); lags; only_cur]

    if length(foward_variables) >= length(backward_variables)
        var_names = [repr(vars) for vars in foward_variables]
    else
        var_names = [repr(vars) for vars in backward_variables]
    end

    shock_names = [repr(vars) for vars in shocks]

    #println(backward_variables)
    #println(foward_variables)

    #Building each matrix

    Gamma1 = NamedArray(zeros(length(model),length(model)))
    setnames!(Gamma1,var_names,1)
    setnames!(Gamma1,var_names,2)
    Gamma0 = NamedArray(zeros(length(model),length(model)))
    setnames!(Gamma0,var_names,1)
    setnames!(Gamma0,var_names,1)
    Psi = NamedArray(zeros(length(model),length(shocks)))
    setnames!(Psi,var_names,1)
    setnames!(Psi,shock_names,2)
    Pi =  NamedArray(zeros(length(model),length(expec_errors)))
    setnames!(Pi,var_names,1)
    setnames!(Pi,expec_errors_names,2)

    #Gamma 0

    for j in 1:length(foward_variables)
        for i in keys(model)
            index = findfirst(keys_model .== i)
            Gamma0[index,j] = -coefficient(model[i],foward_variables[j])
        end
    end

    #Gamma 1

    for j in 1:length(backward_variables)
        for i in keys(model)
            index = findfirst(keys_model .== i)
            Gamma1[index,j] = coefficient(model[i],backward_variables[j])
        end
    end

    #Psi

    for j in 1:length(shocks)
        for i in keys(model)
            index = findfirst(keys_model .== i)
            Psi[index,j] = coefficient(model[i],shocks[j])
        end
    end

    #Pi

    for j in 1:length(expec_errors)
        for i in keys(model)
            index = findfirst(keys_model .== i)
            Pi[index,j] = coefficient(model[i],expec_errors[j])
        end
    end
    return gensys_sys(Gamma0,Gamma1,Psi,Pi,backward_variables,foward_variables)
end

function gensys(gensys_sys::gensys_sys)
    var_names = names(gensys_sys.Gamma0,1)
    ans = gensys(Array(gensys_sys.Gamma0),Array(gensys_sys.Gamma1),Array(gensys_sys.Psi),Array(gensys_sys.Pi))
    T1 = NamedArray(ans.Theta1)
    setnames!(T1,var_names,1)
    T2 = NamedArray(ans.Theta2)
    setnames!(T2,var_names,1)
    T3 = NamedArray(ans.Theta3)
    setnames!(T3,var_names,1)
    return Sims(T1,T2,T3,ans.eu)
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
# speed of light in a vacuum in cm s^-1
const c0 = 2.99792458e10

# Gravity constant in g^-1 cm^3 s^-2
const G = 6.6738399e-8

# atomic mass unit in g
const m_H = 1.660539040e-24

# Boltzmann constant in erg K^-1
const k_B = 1.38064852e-16

# Stefan-Boltzmann constant in erg cm^-2 s^-1 K^-4
const sigma_sb = 5.6704e-5

# radiation density constant in erg cm^-3 K^-4
const a_rad = 4 * sigma_sb / c0

# solar mass in g
const Msun = 1.98e33

# solar radius in cm
const Rsun = 6.95508e10

# solar luminosity in erg s^-1
const Lsun = 3.845999e33


text:
module OnlineAITests

import OnlineAI
import OnlineStats
import FactCheck
FactCheck.clear_results()

# sev = OnlineStats.log_severity()
# OnlineStats.log_severity!(OnlineStats.ErrorSeverity)  # turn off most logging

include("nnet_tests.jl")
include("ensembles_test.jl")

# # put logging back the way it was
# OnlineStats.log_severity!(sev)

FactCheck.exitstatus()

end # module


text:
module Operations

using TOML
using Base.Random: UUID
using Base: LibGit2
using Base: Pkg
using Base.Pkg.Types: VersionSet
using Pkg3.Types

function parse_toml(path::String...; fakeit::Bool=false)
    p = joinpath(path...)
    !fakeit || isfile(p) ? TOML.parsefile(p) : Dict{Any,Any}()
end

user_depot() = abspath(homedir(), ".julia")
depots() = Base.Loading.DEPOTS

function registries(depot::String)
    d = joinpath(depot, "registries")
    regs = filter!(readdir(d)) do r
        isfile(joinpath(d, r, "registry.toml"))
    end
    return map(reg->joinpath(depot, "registries", reg), regs)
end
registries() = [r for d in depots() for r in registries(d)]

const line_re = r"""
    ^ \s*
    ([0-9a-f]{8}-[0-9a-f]{4}-[0-9a-f]{4}-[0-9a-f]{4}-[0-9a-f]{12})
    \s* = \s* \{
    \s* name \s* = \s* "([^"]*)" \s*,
    \s* path \s* = \s* "([^"]*)" \s*,?
    \s* \} \s* $
"""x

function find_registered(names::Vector{String})
    # name --> uuid --> paths
    where = Dict{String,Dict{UUID,Vector{String}}}()
    isempty(names) && return where
    names_re = let p = "\\bname\\s*=\\s*\"(?:$(names[1])"
        for i in 2:length(names)
            p *= "|$(names[i])"
        end
        p *= ")\""
        Regex(p)
    end
    for registry in registries()
        packages_file = joinpath(registry, "registry.toml")
        open(packages_file) do io
            for line in eachline(io)
                ismatch(r"^\s*\[\s*packages\s*\]\s*$", line) && break
            end
            for line in eachline(io)
                ismatch(names_re, line) || continue
                m = match(line_re, line)
                m == nothing && error("misformated packages.toml file: $line")
                uuid = UUID(m.captures[1])
                name = Base.unescape_string(m.captures[2])
                path = Base.unescape_string(m.captures[3])
                if !haskey(where, name)
                    where[name] = Dict{UUID,Vector{String}}()
                end
                if !haskey(where[name], uuid)
                    where[name][uuid] = String[]
                end
                push!(where[name][uuid], abspath(registry, path))
            end
        end
    end
    return where
end
find_registered(name::String) = find_registered([name])[name]

function find_installed(uuid::UUID, sha1::SHA1)
    for depot in depots()
        path = abspath(depot, "packages", string(uuid), string(sha1))
        ispath(path) && return path
    end
    return abspath(user_depot(), "packages", string(uuid), string(sha1))
end

include("libgit2_discover.jl")

default_env() = get(ENV, "JULIA_ENV", nothing)
find_config() = find_config(default_env())

const config_names = ["JuliaConfig.toml", "Config.toml"]
const manifest_names = ["JuliaManifest.toml", "Manifest.toml"]
const default_envs = [
    "v$(VERSION.major).$(VERSION.minor).$(VERSION.patch)",
    "v$(VERSION.major).$(VERSION.minor)",
    "v$(VERSION.major)",
    "default",
]

function find_project_env(start_path::String = pwd())
    path = LibGit2.discover(start_path, ceiling = homedir())
    repo = LibGit2.GitRepo(path)
    work = LibGit2.workdir(repo)
    for name in config_names
        path = abspath(work, name)
        isfile(path) && return path
    end
    return abspath(work, config_names[end])
end

function find_default_env()
    for depot in depots(), env in default_envs, name in config_names
        path = joinpath(depot, "environments", env, name)
        isfile(path) && return path
    end
    env = VERSION.major == 0 ? default_envs[2] : default_envs[3]
    return joinpath(user_depot(), "environments", env, config_names[end])
end

function find_config(env::String)
    if isempty(env)
        error("invalid environment name: \"\"")
    elseif env == "/"
        return find_default_env()
    elseif env == "."
        return find_project_env()
    elseif startswith(env, "/") || startswith(env, "./")
        # path to config file or project directory
        splitext(env)[2] == ".toml" && return abspath(env)
        for name in config_names
            path = abspath(env, name)
            isfile(path) && return path
        end
        return abspath(env, config_names[end])
    else # named environment
        for depot in depots()
            path = joinpath(depot, "environments", env, "Config.toml")
            isfile(path) && return path
        end
        return joinpath(user_depot(), "environments", env, "Config.toml")
    end
end

function find_config(::Void)
    try
        return find_project_env()
    catch err
        err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
    end
    return find_default_env()
end

function find_manifest(env::Union{Void,String}=default_env())
    config_file = find_config(env)
    if isfile(config_file)
        config = parse_toml(config_file)
        haskey(config, "manifest") &&
        return abspath(config["manifest"])
    end
    names = ["JuliaManifest.toml", "Manifest.toml"]
    dir = dirname(config_file)
    for name in names
        path = joinpath(dir, name)
        isfile(path) && return path
    end
    return joinpath(dir, names[end])
end

function load_manifest(env::Union{Void,String}=default_env())
    manifest = find_manifest(env)
    T = Dict{String,Dict{String,String}}
    return isfile(manifest) ? convert(T, parse_toml(manifest)) : T()
end

inspec(v::VersionNumber, s::VersionSpec) = v in s
inspec(v::VersionNumber, s::VersionNumber) = v == s

function parse_version_set(s::String)
    parts = split(s, '-')
    length(parts) == 1 && return VersionSet(VersionSpec(parts[1]))
    length(parts) != 2 && error("invalid version spec: ", repr(s))
    lower = VersionSet(VersionSpec(parts[1])).intervals[1].lower
    upper = VersionSet(VersionSpec(parts[2])).intervals[1].upper
    return VersionSet(lower, upper)
end

function add(names...; kwargs...)
    pkgs = Dict{String,Union{VersionNumber,VersionSpec}}()
    for name in names
        pkgs[string(name)] = vs"*"
    end
    for (key, val) in kwargs
        pkgs[string(key)] = val
    end
    return add(pkgs)
end

function add(pkgs::Dict{String})
    orig_pkgs = copy(pkgs)
    names = sort!(collect(keys(pkgs)))
    regs = find_registered(names)
    manifest = load_manifest()
    uuids = Dict{String,UUID}(name => info["uuid"] for (name, info) in manifest)
    # disambiguate package names
    info("Resolving package UUIDs")
    let ambig = false
        for name in names
            if haskey(uuids, name)
                uuid = uuids[name]
                if haskey(regs, name)
                    for uuid′ in collect(keys(regs[name]))
                        uuid′ == uuid || delete!(regs[name], uuid)
                    end
                end
                haskey(regs, name) && !isempty(regs[name]) && continue
                error("""
                $name/$uuid found in manifest but not in registries;
                To install a different $name package `pkg rm $name` and then do `pkg add $name` again.
                """)
            elseif length(regs[name]) == 1
                uuids[name] = first(first(regs[name]))
            else
                msg = "$name is ambiguous, it could refer to:\n"
                for (i, (uuid, paths)) in enumerate(sort!(collect(regs[name]), by=first))
                    msg *= " [$i] $uuid"
                    for path in paths
                        info = parse_toml(path, "package.toml")
                        msg *= " – $(info["repo"])"
                        break
                    end
                    msg *= "\n"
                end
                info(msg)
                ambig = true
            end
        end
        ambig && error("interactive package choice not yet implemented")
    end
    merge!(regs, find_registered(setdiff(keys(uuids), names)))
    where = isempty(regs) ? Dict{String,Vector{String}}() :
        Dict(name => paths for (name, info) in regs for (uuid, paths) in info)

    # copy manifest versions to pkgs
    # if already in specified version set, leave as installed
    for (name, info) in manifest
        ver = VersionNumber(info["version"])
        if !haskey(pkgs, name) || inspec(ver, pkgs[name])
            pkgs[name] = ver
        end
    end

    # compute reqs & deps for Pkg.Resolve.resolve
    #  - reqs: what we need to choose versions for
    #  - deps: relevant portion of dependency graph

    # reqs :: String --> VersionSet
    reqs = convert(Pkg.Types.Requires, pkgs)

    info("Resolving package versions")
    # deps :: String --> VersionNumber --> (SHA1, String --> VersionSet)
    deps = Dict{String,Dict{VersionNumber,Pkg.Types.Available}}()
    names = sort!(collect(keys(uuids)))
    for name in names
        spec = get(pkgs, name, VersionSpec())
        uuid, paths = uuids[name], where[name]
        deps[name] = Dict{VersionNumber,Pkg.Types.Available}()
        for path in paths
            versions = parse_toml(path, "versions.toml")
            requires = parse_toml(path, "requirements.toml", fakeit=true)
            for (v, d) in versions
                ver = VersionNumber(v)
                inspec(ver, spec) || continue
                r = haskey(requires, v) ?
                    Dict(dep => parse_version_set(r) for (dep, r) in requires[v]) :
                    Pkg.Types.Requires()
                if haskey(r, "julia")
                    VERSION in r["julia"] || continue
                    delete!(r, "julia")
                end
                x = deps[name][ver] = Pkg.Types.Available(SHA1(d["hash-sha1"]), r)
                for dep in keys(x.requires)
                    (dep in names || dep == "julia") && continue
                    found = find_registered(dep)
                    @assert length(found) == 1 # TODO: use UUIDs to disambiguate
                    uuids[dep] = first(found)[1]
                    where[dep] = first(found)[2]
                    push!(names, dep)
                end
            end
        end
    end
    deps = Pkg.Query.prune_dependencies(reqs, deps)
    vers = Pkg.Resolve.resolve(reqs, deps)

    # find repos and hashes for each package & version
    repos = Dict{String,Vector{String}}()
    hashes = Dict{String,SHA1}()
    for (name, ver) in vers
        repos[name] = String[]
        v = string(ver)
        for path in where[name]
            package = parse_toml(path, "package.toml")
            repo = package["repo"]
            repo in repos[name] || push!(repos[name], repo)
            versions = parse_toml(path, "versions.toml")
            if haskey(versions, v)
                sha1 = SHA1(versions[v]["hash-sha1"])
                if haskey(hashes, name)
                    sha1 == hashes[name] || warn("$name: hash mismatch for version $v!")
                else
                    hashes[name] = sha1
                end
            end
        end
        @assert haskey(hashes, name)
    end
    foreach(sort!, values(repos))

    # clone or update repos and find or create source trees
    for (name, sha1) in hashes
        uuid = uuids[name]
        version_path = find_installed(uuid, sha1)
        ispath(version_path) && continue
        repo_path = joinpath(user_depot(), "upstream", string(uuid))
        urls = copy(repos[name])
        git_hash = LibGit2.GitHash(sha1.bytes)
        repo = ispath(repo_path) ? LibGit2.GitRepo(repo_path) : begin
            info("Cloning [$uuid] $name")
            LibGit2.clone(shift!(urls), repo_path, isbare=true)
        end
        refspecs = ["+refs/*:refs/remotes/cache/*"]
        while !isempty(urls)
            try
                LibGit2.GitObject(repo, git_hash)
                break
            catch err
                err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
            end
            info("Updating $name from $(urls[1])")
            LibGit2.fetch(repo, remoteurl=shift!(urls), refspecs=refspecs)
        end
        obj = try LibGit2.GitObject(repo, git_hash)
        catch err
            err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
            error("$name: git object $(string(sha1)) could not be found")
        end
        tree = LibGit2.peel(LibGit2.GitTree, obj)
        sha1′ = SHA1(string(LibGit2.GitHash(tree)))
        if sha1 != sha1′
            # SHA1 was a commit instead of tree
            sha1 = hashes[name] = sha1′
            version_path = find_installed(uuid, sha1)
            ispath(version_path) && continue
        end
        mkpath(version_path)
        opts = LibGit2.CheckoutOptions(
            checkout_strategy = LibGit2.Consts.CHECKOUT_FORCE,
            target_directory = Base.unsafe_convert(Cstring, version_path)
        )
        info("Installing $name at $(string(sha1))")
        LibGit2.checkout_tree(repo, tree, options=opts)
    end

    # update config and manifest files
    config_file = find_config()
    config = parse_toml(config_file, fakeit=true)
    if !haskey(config, "deps")
        config["deps"] = Dict{String,String}()
    end
    for (name, sha1) in hashes
        uuid = string(uuids[name])
        version = string(vers[name])
        config["deps"][name] = uuid
        manifest[name] = Dict{String,String}(
            "uuid"      => uuid,
            "hash-sha1" => sha1,
            "version"   => version,
        )
    end
    isempty(config["deps"]) && delete!(config, "deps")

    if !isempty(config) || ispath(config_file)
        mkpath(dirname(config_file))
        info("Updating config file $config_file")
        open(config_file, "w") do io
            TOML.print(io, config, sorted=true)
        end
    end
    manifest_file = find_manifest()
    if !isempty(manifest) || ispath(manifest_file)
        info("Updating manifest file $manifest_file")
        open(manifest_file, "w") do io
            TOML.print(io, manifest, sorted=true)
        end
    end
end

end # module