text:
"""
    Calculators

This module exists to bridge the gap between the `Models` and the `Dynamics`.

Here, we provide functions and types for evaluating and storing quantities obtained from the
`Models`.
In addition any further manipulation of those quantities, such as computing eigenvalues,
is included here.

This module is largely needed to facilitate integration of both ring polymer and
classical dynamics to allow using the same models and functions for both.
Specific ring polymer types are provided that have the extra fields and methods needed
to evaluate the quantities for each bead.
"""
module Calculators

using LinearAlgebra: LinearAlgebra, Hermitian, I, Eigen, tr
using StaticArrays: SMatrix, SVector
using RingPolymerArrays: get_centroid!

using NQCModels: NQCModels, Model, nstates
using NQCModels.AdiabaticModels: AdiabaticModel
using NQCModels.DiabaticModels: DiabaticModel, DiabaticFrictionModel
using NQCModels.FrictionModels: AdiabaticFrictionModel

using NQCDynamics: ndofs

"""
    AbstractCalculator{M<:Model}

Top-level type for all calculators.

Each concrete calculator contains the `Model` and the fields to store the quantities
obtained from the model.
"""
abstract type AbstractCalculator{T,M<:Model} end
abstract type AbstractAdiabaticCalculator{T,M<:AdiabaticModel} <: AbstractCalculator{T,M} end
abstract type AbstractDiabaticCalculator{T,M<:Union{DiabaticFrictionModel,DiabaticModel}} <: AbstractCalculator{T,M} end
abstract type AbstractStaticDiabaticCalculator{T,M} <: AbstractDiabaticCalculator{T,M} end
abstract type AbstractFrictionCalculator{T,M<:AdiabaticFrictionModel} <: AbstractCalculator{T,M} end

NQCModels.nstates(calc::AbstractCalculator) = NQCModels.nstates(calc.model)
NQCModels.eachstate(calc::AbstractCalculator) = NQCModels.eachstate(calc.model)
NQCModels.nelectrons(calc::AbstractCalculator) = NQCModels.nelectrons(calc.model)
NQCModels.eachelectron(calc::AbstractCalculator) = NQCModels.eachelectron(calc.model)
NQCModels.mobileatoms(calc::AbstractCalculator) = NQCModels.mobileatoms(calc.model, size(calc.derivative, 2))
NQCModels.dofs(calc::AbstractCalculator) = NQCModels.dofs(calc.model)

Base.eltype(::AbstractCalculator{T}) where {T} = T

mutable struct DependentField{V,R}
    value::V
    position::R
end
needsupdate(field::DependentField, r) = field.position != r

function Base.getproperty(calc::AbstractCalculator, name::Symbol)
    field = getfield(calc, name)
    if field isa DependentField
        return field.value
    else
        return field
    end
end

function Base.setproperty!(calc::AbstractCalculator, name::Symbol, x)
    field = getfield(calc, name)
    if field isa DependentField
        setfield!(field, :value, x)
    else
        setfield!(calc, name, x)
    end
end

"""
Each of the quantities specified here has functions:
`get_quantity(calculator, r)`
`evaluate_quantity!(calculator, r)!`

The user should access only the former.
This will ensure quantities are correctly evaluated and cached accordingly.

The latter is called by the former and is where the details required to calculate the quantity are found.
"""
const quantities = [
    :potential,
    :derivative,
    :eigen,
    :adiabatic_derivative,
    :nonadiabatic_coupling,

    :traceless_potential,
    :V̄,
    :traceless_derivative,
    :D̄,
    :traceless_adiabatic_derivative,

    :centroid,
    :centroid_potential,
    :centroid_derivative,
    :centroid_eigen,
    :centroid_adiabatic_derivative,
    :centroid_nonadiabatic_coupling,

    :friction,
]

for quantity in quantities
    get_quantity = Symbol(:get_, quantity)
    evaluate_quantity! = Symbol(:evaluate_, quantity, :!)
    field = Expr(:call, :getfield, :calculator, QuoteNode(quantity))

    @eval function $(get_quantity)(calculator, r)

        if needsupdate($field, r)
            copyto!($field.position, r)
            $(evaluate_quantity!)(calculator, r)
        end
        return $field.value
    end
end

struct AdiabaticCalculator{T,M} <: AbstractAdiabaticCalculator{T,M}
    model::M
    potential::DependentField{T,Matrix{T}}
    derivative::DependentField{Matrix{T},Matrix{T}}
    stats::Dict{Symbol,Int}
    function AdiabaticCalculator{T}(model::M, atoms::Integer) where {T,M<:Model}
        potential = zero(T)
        derivative = zeros(T, ndofs(model), atoms)
        position = fill(NaN, ndofs(model), atoms)

        stats = Dict{Symbol,Int}(
            :potential=>0,
            :derivative=>0,
        )

        new{T,M}(
            model,
            DependentField(potential, copy(position)),
            DependentField(derivative, copy(position)),
            stats
        )
    end
end

struct RingPolymerAdiabaticCalculator{T,M} <: AbstractAdiabaticCalculator{T,M}
    model::M
    potential::DependentField{Vector{T},Array{T,3}}
    derivative::DependentField{Array{T,3},Array{T,3}}
    stats::Dict{Symbol,Int}
    function RingPolymerAdiabaticCalculator{T}(model::M, atoms::Integer, beads::Integer) where {T,M<:Model}
        potential = zeros(T, beads)
        derivative = zeros(T, ndofs(model), atoms, beads)
        position = fill(NaN, ndofs(model), atoms, beads)

        stats = Dict{Symbol,Int}(
            :potential=>0,
            :derivative=>0,
        )

        new{T,M}(
            model,
            DependentField(potential, copy(position)),
            DependentField(derivative, copy(position)),
            stats
        )
    end
end

struct DiabaticCalculator{T,M,S,L} <: AbstractDiabaticCalculator{T,M}
    model::M
    potential::DependentField{Hermitian{T,SMatrix{S,S,T,L}},Matrix{T}}
    derivative::DependentField{Matrix{Hermitian{T,SMatrix{S,S,T,L}}},Matrix{T}}
    eigen::DependentField{LinearAlgebra.Eigen{T,T,SMatrix{S,S,T,L},SVector{S,T}},Matrix{T}}
    adiabatic_derivative::DependentField{Matrix{SMatrix{S,S,T,L}},Matrix{T}}
    nonadiabatic_coupling::DependentField{Matrix{SMatrix{S,S,T,L}},Matrix{T}}
    stats::Dict{Symbol,Int}
    function DiabaticCalculator{T}(model::M, atoms::Integer) where {T,M<:Model}
        n = nstates(model)
        mat = NQCModels.DiabaticModels.matrix_template(model, T)
        vec = NQCModels.DiabaticModels.vector_template(model, T)

        potential = Hermitian(zero(mat))
        derivative = [Hermitian(zero(mat)) for _=1:ndofs(model), _=1:atoms]
        eigen = Eigen(zero(vec), zero(mat) + I)
        adiabatic_derivative = [zero(mat) for _ in CartesianIndices(derivative)]
        nonadiabatic_coupling = [zero(mat) for _ in CartesianIndices(derivative)]
        position = fill(NaN, ndofs(model), atoms)

        stats = Dict{Symbol,Int}(
            :potential=>0,
            :derivative=>0,
            :eigen=>0,
            :adiabatic_derivative=>0,
            :nonadiabatic_coupling=>0
        )

        new{T,M,n,n^2}(
            model,
            DependentField(potential, copy(position)),
            DependentField(derivative, copy(position)),
            DependentField(eigen, copy(position)),
            DependentField(adiabatic_derivative, copy(position)),
            DependentField(nonadiabatic_coupling, copy(position)),
            stats
        )
    end
end

struct RingPolymerDiabaticCalculator{T,M,S,L} <: AbstractDiabaticCalculator{T,M}
    model::M
    potential::DependentField{Vector{Hermitian{T,SMatrix{S,S,T,L}}},Array{T,3}}
    derivative::DependentField{Array{Hermitian{T,SMatrix{S,S,T,L}},3},Array{T,3}}
    eigen::DependentField{Vector{LinearAlgebra.Eigen{T,T,SMatrix{S,S,T,L},SVector{S,T}}},Array{T,3}}
    adiabatic_derivative::DependentField{Array{SMatrix{S,S,T,L},3},Array{T,3}}
    nonadiabatic_coupling::DependentField{Array{SMatrix{S,S,T,L},3},Array{T,3}}

    traceless_potential::DependentField{Vector{Hermitian{T,SMatrix{S,S,T,L}}},Array{T,3}}
    V̄::DependentField{Vector{T},Array{T,3}}
    traceless_derivative::DependentField{Array{Hermitian{T,SMatrix{S,S,T,L}},3},Array{T,3}}
    D̄::DependentField{Array{T,3},Array{T,3}}
    traceless_adiabatic_derivative::DependentField{Array{SMatrix{S,S,T,L},3},Array{T,3}}

    centroid::DependentField{Matrix{T},Array{T,3}}
    centroid_potential::DependentField{Hermitian{T,SMatrix{S,S,T,L}},Array{T,3}}
    centroid_derivative::DependentField{Matrix{Hermitian{T,SMatrix{S,S,T,L}}},Array{T,3}}
    centroid_eigen::DependentField{LinearAlgebra.Eigen{T,T,SMatrix{S,S,T,L},SVector{S,T}},Array{T,3}}
    centroid_adiabatic_derivative::DependentField{Matrix{SMatrix{S,S,T,L}},Array{T,3}}
    centroid_nonadiabatic_coupling::DependentField{Matrix{SMatrix{S,S,T,L}},Array{T,3}}

    stats::Dict{Symbol,Int}
    function RingPolymerDiabaticCalculator{T}(model::M, atoms::Integer, beads::Integer) where {T,M<:Model}
        n = nstates(model)
        mat = NQCModels.DiabaticModels.matrix_template(model, T)
        vec = NQCModels.DiabaticModels.vector_template(model, T)

        potential = [Hermitian(zero(mat)) for _=1:beads]
        derivative = [Hermitian(zero(mat)) for _=1:ndofs(model), _=1:atoms, _=1:beads]
        adiabatic_derivative = [zero(mat) for _=1:ndofs(model), _=1:atoms, _=1:beads]
        eigen = [Eigen(zero(vec), zero(mat) + I) for _=1:beads]
        nonadiabatic_coupling = [zero(mat) for _=1:ndofs(model), _=1:atoms, _=1:beads]

        traceless_potential = [Hermitian(zero(mat)) for _=1:beads]
        V̄ = zeros(beads)
        traceless_derivative = [Hermitian(zero(mat)) for _=1:ndofs(model), _=1:atoms, _=1:beads]
        D̄ = zeros(ndofs(model), atoms, beads)
        traceless_adiabatic_derivative = [zero(mat) for _=1:ndofs(model), _=1:atoms, _=1:beads]

        centroid = zeros(T, ndofs(model), atoms)
        centroid_potential = Hermitian(zero(mat))
        centroid_derivative = [Hermitian(zero(mat)) for _=1:ndofs(model), _=1:atoms]
        centroid_eigen = Eigen(zero(vec), zero(mat) + I)
        centroid_adiabatic_derivative = [zero(mat) for _ in CartesianIndices(centroid_derivative)]
        centroid_nonadiabatic_coupling = [zero(mat) for _ in CartesianIndices(centroid_derivative)]

        position = fill(NaN, ndofs(model), atoms, beads)

        stats = Dict{Symbol,Int}(
            :potential=>0,
            :derivative=>0,
            :eigen=>0,
            :adiabatic_derivative=>0,
            :nonadiabatic_coupling=>0,
            :traceless_potential=>0,
            :V̄=>0,
            :traceless_derivative=>0,
            :D̄=>0,
            :traceless_adiabatic_derivative=>0,
            :centroid=>0,
            :centroid_potential=>0,
            :centroid_derivative=>0,
            :centroid_eigen=>0,
            :centroid_adiabatic_derivative=>0,
            :centroid_nonadiabatic_coupling=>0,
        )

        new{T,M,n,n^2}(
            model,
            DependentField(potential, copy(position)),
            DependentField(derivative, copy(position)),
            DependentField(eigen, copy(position)),
            DependentField(adiabatic_derivative, copy(position)),
            DependentField(nonadiabatic_coupling, copy(position)),

            DependentField(traceless_potential, copy(position)),
            DependentField(V̄, copy(position)),
            DependentField(traceless_derivative, copy(position)),
            DependentField(D̄, copy(position)),
            DependentField(traceless_adiabatic_derivative, copy(position)),

            DependentField(centroid, copy(position)),
            DependentField(centroid_potential, copy(position)),
            DependentField(centroid_derivative, copy(position)),
            DependentField(centroid_eigen, copy(position)),
            DependentField(centroid_adiabatic_derivative, copy(position)),
            DependentField(centroid_nonadiabatic_coupling, copy(position)),
            stats
        )
    end
end

function Calculator(model::DiabaticModel, atoms::Integer, t::Type{T}) where {T}
    DiabaticCalculator{t}(model, atoms)
end
function Calculator(model::AdiabaticModel, atoms::Integer, t::Type{T}) where {T}
    AdiabaticCalculator{t}(model, atoms)
end
function Calculator(model::DiabaticModel, atoms::Integer, beads::Integer, t::Type{T}) where {T}
    RingPolymerDiabaticCalculator{t}(model, atoms, beads)
end
function Calculator(model::AdiabaticModel, atoms::Integer, beads::Integer, t::Type{T}) where {T}
    RingPolymerAdiabaticCalculator{t}(model, atoms, beads)
end

function evaluate_potential!(calc::AbstractCalculator, R)
    calc.stats[:potential] += 1
    calc.potential = NQCModels.potential(calc.model, R)
    return nothing
end

function evaluate_potential!(calc::AbstractCalculator, R::AbstractArray{T,3}) where {T}
    calc.stats[:potential] += 1
    @views for i in axes(R, 3)
        calc.potential[i] = NQCModels.potential(calc.model, R[:,:,i])
    end
    return nothing
end

function evaluate_V̄!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:V̄] += 1
    potential = get_potential(calc, r)
    for i in 1:length(calc.V̄)
        calc.V̄[i] = tr(potential[i]) / nstates(calc.model)
    end
end

function evaluate_traceless_potential!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:traceless_potential] += 1
    n = nstates(calc.model)
    potential = get_potential(calc, r)
    V̄ = get_V̄(calc, r)
    for I in eachindex(potential)
        calc.traceless_potential[I] = Hermitian(SMatrix{n,n}(
            i != j ? potential[I][j,i] : potential[I][j,i] - V̄[I] for j=1:n, i=1:n
        ))
    end
end

function evaluate_centroid!(calc::AbstractCalculator, r::AbstractArray{T,3}) where {T}
    calc.stats[:centroid] += 1
    get_centroid!(calc.centroid, r)
end

function evaluate_centroid_potential!(calc::AbstractCalculator, r::AbstractArray{T,3}) where {T}
    calc.stats[:centroid_potential] += 1
    centroid = get_centroid(calc, r)
    calc.centroid_potential = NQCModels.potential(calc.model, centroid)
end

function evaluate_derivative!(calc::AbstractCalculator, R)
    calc.stats[:derivative] += 1
    NQCModels.derivative!(calc.model, calc.derivative, R)
end

function evaluate_derivative!(calc::AbstractCalculator, R::AbstractArray{T,3}) where {T}
    calc.stats[:derivative] += 1
    @views for i in axes(R, 3)
        NQCModels.derivative!(calc.model, calc.derivative[:,:,i], R[:,:,i])
    end
end

function evaluate_D̄!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:D̄] += 1
    derivative = get_derivative(calc, r)
    for I in eachindex(derivative)
        calc.D̄[I] = tr(derivative[I]) / nstates(calc.model)
    end
end

function evaluate_traceless_derivative!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:traceless_derivative] += 1
    n = nstates(calc.model)
    derivative = get_derivative(calc, r)
    D̄ = get_D̄(calc, r)
    for I in eachindex(derivative)
        calc.traceless_derivative[I] = Hermitian(SMatrix{n,n}(
            i != j ? derivative[I][j,i] : derivative[I][j,i] - D̄[I] for j=1:n, i=1:n
        ))
    end
end

function evaluate_traceless_adiabatic_derivative!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:traceless_adiabatic_derivative] += 1
    n = nstates(calc.model)
    adiabatic_derivative = get_adiabatic_derivative(calc, r)
    D̄ = get_D̄(calc, r)
    for I in eachindex(D̄)
        calc.traceless_adiabatic_derivative[I] = SMatrix{n,n}(
            i != j ? adiabatic_derivative[I][j,i] : adiabatic_derivative[I][j,i] - D̄[I] for j=1:n, i=1:n
        )
    end
end

function evaluate_centroid_derivative!(calc::AbstractCalculator, r::AbstractArray{T,3}) where {T}
    calc.stats[:centroid_derivative] += 1
    centroid = get_centroid(calc, r)
    NQCModels.derivative!(calc.model, calc.centroid_derivative, centroid)
end

function evaluate_eigen!(calc::AbstractDiabaticCalculator, r)
    calc.stats[:eigen] += 1
    potential = get_potential(calc, r)
    eig = LinearAlgebra.eigen(potential)
    corrected_vectors = correct_phase(eig.vectors, calc.eigen.vectors)
    calc.eigen = Eigen(eig.values, corrected_vectors)
    return nothing
end

function evaluate_centroid_eigen!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:centroid_eigen] += 1
    potential = get_centroid_potential(calc, r)
    eig = LinearAlgebra.eigen(potential)
    corrected_vectors = correct_phase(eig.vectors, calc.centroid_eigen.vectors)
    calc.centroid_eigen = Eigen(eig.values, corrected_vectors)
    return nothing
end

function evaluate_eigen!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:eigen] += 1
    potential = get_potential(calc, r)
    for i=1:length(potential)
        eig = LinearAlgebra.eigen(potential[i])
        corrected_vectors = correct_phase(eig.vectors, calc.eigen[i].vectors)
        calc.eigen[i] = Eigen(eig.values, corrected_vectors)
    end
    return nothing
end

function correct_phase(new_vectors::SMatrix, old_vectors::SMatrix)
    n = size(new_vectors, 1)
    vect = SVector{n}(sign(LinearAlgebra.dot(new_vectors[:,i], old_vectors[:,i])) for i=1:n)
    return new_vectors .* vect'
end

function evaluate_adiabatic_derivative!(calc::AbstractDiabaticCalculator, r)
    calc.stats[:adiabatic_derivative] += 1
    U = get_eigen(calc, r).vectors
    diabatic_derivative = get_derivative(calc, r)
    for I in eachindex(diabatic_derivative)
        calc.adiabatic_derivative[I] = U' * diabatic_derivative[I] * U
    end
end

function evaluate_centroid_adiabatic_derivative!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:centroid_adiabatic_derivative] += 1
    centroid_derivative = get_centroid_derivative(calc, r)
    centroid_eigen = get_centroid_eigen(calc, r)
    for I in eachindex(centroid_derivative)
        calc.centroid_adiabatic_derivative[I] = centroid_eigen.vectors' * centroid_derivative[I] * centroid_eigen.vectors
    end
end

function evaluate_adiabatic_derivative!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:adiabatic_derivative] += 1
    derivative = get_derivative(calc, r)
    eigen = get_eigen(calc, r)
    for i in axes(derivative, 3) # Beads
        for j in axes(derivative, 2) # Atoms
            for k in axes(derivative, 1) # DoFs
                calc.adiabatic_derivative[k,j,i] = eigen[i].vectors' * derivative[k,j,i] * eigen[i].vectors
            end
        end
    end
end

function evaluate_nonadiabatic_coupling!(calc::AbstractDiabaticCalculator, r)
    calc.stats[:nonadiabatic_coupling] += 1
    adiabatic_derivative = get_adiabatic_derivative(calc, r)
    eigen = get_eigen(calc, r)
    for I in eachindex(calc.adiabatic_derivative)
        calc.nonadiabatic_coupling[I] = evaluate_nonadiabatic_coupling(adiabatic_derivative[I], eigen.values)
    end
end

function evaluate_centroid_nonadiabatic_coupling!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:centroid_nonadiabatic_coupling] += 1
    centroid_adiabatic_derivative = get_centroid_adiabatic_derivative(calc, r)
    centroid_eigen = get_centroid_eigen(calc, r)
    for I in eachindex(centroid_adiabatic_derivative)
        calc.centroid_nonadiabatic_coupling[I] = evaluate_nonadiabatic_coupling(centroid_adiabatic_derivative[I], centroid_eigen.values)
    end
end

function evaluate_nonadiabatic_coupling!(calc::RingPolymerDiabaticCalculator, r)
    calc.stats[:nonadiabatic_coupling] += 1
    adiabatic_derivative = get_adiabatic_derivative(calc, r)
    eigen = get_eigen(calc, r)
    for i in eachindex(eigen) # Beads
        for I in CartesianIndices(size(adiabatic_derivative)[1:2])
            calc.nonadiabatic_coupling[I,i] = evaluate_nonadiabatic_coupling(adiabatic_derivative[I,i], eigen[i].values)
        end
    end
end

"""
# References

- HammesSchifferTully_JChemPhys_101_4657_1994 Eq. (32)
- SubotnikBellonzi_AnnuRevPhyschem_67_387_2016, section 2.3
"""
function evaluate_nonadiabatic_coupling(adiabatic_derivative::SMatrix, eigenvalues::SVector)
    n = length(eigenvalues)
    SMatrix{n,n}(
        (i != j ? adiabatic_derivative[j,i] / (eigenvalues[i] - eigenvalues[j]) : 0
        for j=1:n, i=1:n))
end

include("friction.jl")
include("large_diabatic.jl")

"""
Evaluates all electronic properties for the current position `r`.
# Properties evaluated:
- Diabatic potential
- Diabatic derivative
- Eigenvalues and eigenvectors
- Adiabatic derivative
- Nonadiabatic coupling

This should no longer be used, instead access the quantities directly with `get_quantity(calc, r)`.
"""
function update_electronics!(calculator::AbstractDiabaticCalculator, r::AbstractArray)
    get_nonadiabatic_coupling(calculator, r)
    return nothing
end

function update_electronics!(calculator::RingPolymerDiabaticCalculator, r::AbstractArray{T,3}) where {T}
    get_nonadiabatic_coupling(calculator, r)
    update_centroid_electronics!(calculator, r)
    return nothing
end

function update_centroid_electronics!(calculator::RingPolymerDiabaticCalculator, r::AbstractArray{T,3}) where {T}
    get_centroid_nonadiabatic_coupling(calculator, r)
    return nothing
end

end # module


text:
function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
    a = zeros(3, 3)
    a[1, 1] = inradius
    a[1, 2] = -inradius
    a[1, 3] = -1.0 / inradius
    a[2, 1] = 2 * (lc + mc)
    a[2, 2] = -2 * (ls + ms)
    a[2, 3] = 2ms / inradius^2
    a[3, 2] = 2(ls + ms)
    a[3, 3] = -2ms / outradius^2
    return a
end

function analytical_coefficient_rhs(ls, ms, theta0)
    r = zeros(3)
    Ks = bulk_modulus(ls, ms)
    r[2] = -Ks * theta0
    r[3] = Ks * theta0
    return r
end

struct AnalyticalSolution
    inradius::Any
    outradius::Any
    center::Any
    A1c::Any
    A1s::Any
    A2s::Any
    ls::Any
    ms::Any
    lc::Any
    mc::Any
    theta0::Any
    function AnalyticalSolution(
        inradius,
        outradius,
        center,
        ls,
        ms,
        lc,
        mc,
        theta0,
    )
        a = analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
        r = analytical_coefficient_rhs(ls, ms, theta0)
        coeffs = a \ r
        new(
            inradius,
            outradius,
            center,
            coeffs[1],
            coeffs[2],
            coeffs[3],
            ls,
            ms,
            lc,
            mc,
            theta0,
        )
    end
end

function radial_displacement(A::AnalyticalSolution, r)
    if r <= A.inradius
        return A.A1c * r
    else
        return A.A1s * r + A.A2s / r
    end
end

function (A::AnalyticalSolution)(x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    ur = radial_displacement(A, r)
    if ur ≈ 0.0
        [0.0, 0.0]
    else
        costheta = (x[1] - A.center[1]) / r
        sintheta = (x[2] - A.center[2]) / r
        u1 = ur * costheta
        u2 = ur * sintheta
        return [u1, u2]
    end
end

function shell_radial_stress(ls, ms, theta0, A1, A2, r)
    return (ls + 2ms) * (A1 - A2 / r^2) + ls * (A1 + A2 / r^2) -
           (ls + 2ms / 3) * theta0
end

function shell_circumferential_stress(ls, ms, theta0, A1, A2, r)
    return ls * (A1 - A2 / r^2) + (ls + 2ms) * (A1 + A2 / r^2) -
           (ls + 2ms / 3) * theta0
end

function shell_out_of_plane_stress(ls, ms, A1, theta0)
    return 2 * ls * A1 - (ls + 2ms / 3) * theta0
end

function core_in_plane_stress(lc, mc, A1)
    return (lc + 2mc) * A1 + lc * A1
end

function core_out_of_plane_stress(lc, A1)
    return 2 * lc * A1
end

function rotation_matrix(x, r)
    costheta = x[1] / r
    sintheta = x[2] / r
    Q = [
        costheta -sintheta
        sintheta costheta
    ]
    return Q
end

function shell_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    Q = rotation_matrix(relpos, r)

    srr = shell_radial_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)
    stt = shell_circumferential_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)

    cylstress = [
        srr 0.0
        0.0 stt
    ]

    cartstress = Q * cylstress * Q'
    s11 = cartstress[1, 1]
    s22 = cartstress[2, 2]
    s12 = cartstress[1, 2]
    s33 = shell_out_of_plane_stress(A.ls, A.ms, A.A1s, A.theta0)

    return [s11, s22, s12, s33]
end

function core_stress(A::AnalyticalSolution)
    s11 = core_in_plane_stress(A.lc, A.mc, A.A1c)
    s33 = core_out_of_plane_stress(A.lc, A.A1c)
    return [s11, s11, 0.0, s33]
end

function exact_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    if r < A.inradius
        return core_stress(A)
    else
        return shell_stress(A, x)
    end
end

ext: jl


text:
#__precompile__(true)

module IntervalConstraintProgramming

using ValidatedNumerics
using MacroTools
using Compat
using FixedSizeArrays: setindex

import Base:
    show, ∩, ∪, !, ⊆, setdiff

export
    @contractor,
    Separator, separator, @separator, @constraint,
    @function,
    SubPaving, Paving,
    pave, refine!,
    Vol,
    show_code

include("reverse_mode.jl")
include("contractor.jl")
include("separator.jl")
include("paving.jl")
include("setinversion.jl")
include("volume.jl")
include("functions.jl")
include("ast.jl")

end # module


text:
# Command: Element Send Keys
"""
    clear!(element::Element)::Nothing

The Element Send Keys command scrolls into view the form control element and then sends the provided keys to the element.
In case the element is not keyboard-interactable, an element not interactable error is returned.

!!! note

    Keyboard keys are not implemented yet.
"""
function element_keys!(element::Element, value::AbstractString)::Nothing
    @unpack addr, id = element.session
    element_id = element.id
    response = HTTP.post(
        "$addr/session/$id/element/$element_id/value",
        [("Content-Type" => "application/json")],
        JSON3.write("value" => split(value, "")),
    )
    @assert response.status == 200
end


text:
println("いつも2525")
println("AtCoderくん")

max_stars_repo_name: sugitanishi/competitive-programming

max_issues_repo_name: sugitanishi/competitive-programming

max_forks_repo_name: sugitanishi/competitive-programming


text:
type Tree
    op
    value
    subtrees :: Array{Tree, 1}
end

function Base.show(io :: IO, m :: Tree)
    if length(m.subtrees) > 0
        print(io, "[")
    end
    print(io, string(m.op))
    if length(m.subtrees) > 0
        print(io, ", ")
    end
    subs = join([string(subm) for subm in m.subtrees], ", ")
    print(io, subs)
    if length(m.subtrees) > 0
        print(io, "]")
    end
end

function beautify(m :: Tree, indent :: Int64)
    indentation = join(["    " for i in 1:indent])
    string(indentation, m.op, "\n", join([beautify(sub, indent + 1) for sub in m.subtrees]))
end

beautify(m :: Tree) = beautify(m, 0)

pprint(m) = print(beautify(m))

function bottom_up(func :: Function, tree :: Tree)
    func(tree, [bottom_up(func, t) for t in tree.subtrees])
end

function toTree(op)
    Tree(op, Dict(), [])
end

function toTree(skeleton :: Array)
    op = skeleton[1]
    subs = skeleton[2:end]
    Tree(op, Dict(), [toTree(s) for s in subs])
end

function add!(dict :: Dict, dict1 :: Dict)
    for key in keys(dict1)
        dict[key] = vcat(get!(dict, key, []), dict1[key])
    end
    dict
end

function _index(tree :: Tree, inds :: Array)
    ind = Dict{Any, Any}(tree.op => [tree])
    for ind1 in inds
        add!(ind, ind1)
    end
    ind
end

index(tree :: Tree) = bottom_up(_index, tree)

max_stars_repo_name: Non-Contradiction/NextGeneration

max_issues_repo_name: Non-Contradiction/NextGeneration

max_forks_repo_name: Non-Contradiction/NextGeneration


text:
using LightGraphs

function sample_graphs(n, p, k)
    gs = Vector{SimpleGraph{Int}}(undef, k)
    for i in 1:k
        gs[i] = erdos_renyi(n, p, seed=i)
    end
    gs
end

function sample_graphs_t(n, p, k)
    gs = Vector{SimpleGraph{Int}}(undef, k)
    Threads.@threads for i in 1:k
        gs[i] = erdos_renyi(n, p, seed=i)
    end
    gs
end

println("threads: ", Threads.nthreads())
@time sample_graphs(10^5, 10^-4, 32)
@time sample_graphs(10^5, 10^-4, 32)
@time sample_graphs_t(10^5, 10^-4, 32)
@time sample_graphs_t(10^5, 10^-4, 32)


text:
@testset "Linear Stretching" begin

    @testset "constructors" begin
        @test LinearStretching() === LinearStretching(nothing, nothing, 0.0f0, 1.0f0, nothing, nothing, false)
        @test LinearStretching(src_minval=0.1f0, src_maxval=0.9f0, dst_minval=0.05f0, dst_maxval=0.95f0, no_clamp=true) ===
              LinearStretching(0.1f0, 0.9f0, 0.05f0, 0.95f0, nothing, nothing, true)

        @test LinearStretching((0.1f0, 0.9f0)=>(0.2f0, 0.8f0)) === LinearStretching(0.1f0, 0.9f0, 0.2f0, 0.8f0)
        @test LinearStretching((0.1f0, 0.9f0)=>(0.2f0, 0.8f0), no_clamp=true) ===
              LinearStretching(0.1f0, 0.9f0, 0.2f0, 0.8f0, nothing, nothing, true)

        @test LinearStretching(nothing=>(0.2f0, 0.8f0)) === LinearStretching((nothing, nothing)=>(0.2f0, 0.8f0))
        @test LinearStretching(nothing=>(0.2f0, 0.8f0), no_clamp=true) ===
              LinearStretching((nothing, nothing)=>(0.2f0, 0.8f0), no_clamp=true)

        @test LinearStretching((0.1f0, 0.9f0)=>nothing, no_clamp=true) ===
              LinearStretching(0.1f0, 0.9f0, 0.0f0, 1.0f0, nothing, nothing, true)

        @test_throws MethodError LinearStretching(0.1f0, 0.9f0)
        @test_throws MethodError LinearStretching((0.1f0, 0.9f0), (0.0f0, 1.0f0))
        @test_throws ArgumentError LinearStretching((0.9f0, 0.1f0)=>nothing)
        @test_throws ArgumentError LinearStretching(nothing=>(0.9f0, 0.1f0))
        @test_throws ArgumentError LinearStretching((0.9f0, 0.1f0)=>(0.9f0, 0.1f0))
    end

    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        #=
        Stretching an image consisting of a linear ramp should not change the image
        if the specified minimum and maximum values match the minimum and maximum of
        the image.
        =#
        img = T.(collect(reshape(1/100:1/100:1, 10, 10)))
        minval = minimum(img)
        maxval = maximum(img)
        ret =  adjust_histogram(img, LinearStretching(nothing=>(minval, maxval)))
        if T <: Gray{Float32} || T <: Gray{Float64}
            @test all(map((i, r) -> isapprox(i, r), img, ret))
        else
            @test ret == img
        end

        # Verify that NaN is also handled correctly.
        if T <: Gray{Float32} || T <: Gray{Float64}
            img[10] = NaN
            ret = adjust_histogram(img, LinearStretching(nothing=>(minval, maxval)))
            @test isapprox(first(img), first(ret))
            @test isapprox(last(img), last(ret))
            @test isnan(ret[10])
        end

        # Verify that the smallest and largest values match the specified minval and maxval.
        img = T.(collect(reshape(1/100:1/100:1, 10, 10)))
        minval = minimum(img)
        maxval = maximum(img)
        ret =  adjust_histogram(img, LinearStretching(nothing=>(0, 1)))
        @test isapprox(0, first(ret))
        @test isapprox(1, last(ret))
        @test isapprox(0, minimum(ret[.!isnan.(ret)]))
        @test isapprox(1, maximum(ret[.!isnan.(ret)]))

        # Verify that the return type matches the input type.
        img = T.(testimage("mandril_gray"))
        ret = adjust_histogram(img, LinearStretching(nothing=>(0, 1)))
        @test eltype(ret) == eltype(img)
        @test isapprox(0, minimum(ret))
        @test isapprox(1, maximum(ret))

        ret = adjust_histogram(img, LinearStretching(nothing=>(0.2, 0.8)))
        @test eltype(ret) == eltype(img)
        @test isapprox(0.2, minimum(ret))
        @test isapprox(0.8, maximum(ret))

        # Verify that results are correctly clamped to [0.2, 0.9] if it exceeds the range
        ret = adjust_histogram(img, LinearStretching((0.1, 0.8)=>(0.2, 0.9)))
        @test eltype(ret) == eltype(img)
        @test isapprox(T(0.2), minimum(ret))
        @test isapprox(T(0.9), maximum(ret), atol=1e-2)
    end

    # Verify that no_clamp option handles different input types correctly without ArgumentError
    img = Float32.(testimage("mandril_gray"))
    ret_clamp = adjust_histogram(img, LinearStretching((0.3, 0.8)=>(0.1, 1.1), no_clamp=true))
    @test eltype(ret_clamp) == eltype(img)
    @test isapprox(-0.5f0, minimum(ret_clamp))
    @test isapprox(1.272549f0, maximum(ret_clamp))
    ret_noclamp = adjust_histogram(img, LinearStretching((0.3, 0.8)=>(0.1, 1.1), no_clamp=false))
    @test eltype(ret_noclamp) == eltype(img)
    @test isapprox(0.1f0, minimum(ret_noclamp))
    @test isapprox(1.1f0, maximum(ret_noclamp))
    # when no_clamp==true, the output is still clamped by (typemin(T), typemax(T))
    img = N0f8.(testimage("mandril_gray"))
    ret_clamp = adjust_histogram(img, LinearStretching((0.3, 0.8)=>(0.1, 1.1), no_clamp=true))
    @test eltype(ret_clamp) == eltype(img)
    @test isapprox(0.0N0f8, minimum(ret_clamp))
    @test isapprox(1.0N0f8, maximum(ret_clamp))
    ret_noclamp = adjust_histogram(img, LinearStretching((0.3, 0.8)=>(0.1, 1.1), no_clamp=false))
    @test eltype(ret_noclamp) == eltype(img)
    @test isapprox(0.1N0f8, minimum(ret_noclamp))
    @test isapprox(1.0N0f8, maximum(ret_noclamp))
    @test ret_clamp != ret_noclamp

    for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
        #=
        Create a color image that spans a narrow graylevel range.  Then
        quantize the 256 bins down to 32 and determine how many bins have
        non-zero counts.
        =#
        imgg = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10])
        img = colorview(RGB,imgg,imgg,imgg)
        img = T.(img)
        _, counts_before = build_histogram(img,32, minval = 0, maxval = 1)
        nonzero_before = sum(counts_before .!= 0)

        #=
        Stretch the histogram. Then quantize the 256 bins down to 32 and
        verify that all 32 bins have non-zero counts. This will confirm that
        the dynamic range of the original image has been increased.
        =#
        ret = adjust_histogram(img, LinearStretching(nothing=>(0, 1)))
        edges, counts_after = build_histogram(ret, 32, minval = 0, maxval = 1)
        nonzero_after = sum(counts_after .!= 0)
        @test nonzero_before < nonzero_after
        @test nonzero_after == 32
        @test eltype(img) == eltype(ret)
    end

    @testset "deprecations" begin
        @info "four depwarns are expected"
        @test LinearStretching(minval = 0.1) === LinearStretching(dst_minval = 0.1)
        @test LinearStretching(maxval = 0.9) === LinearStretching(dst_maxval = 0.9)
        @test LinearStretching(minval = 0.1, maxval = 0.9) === LinearStretching(dst_minval = 0.1, dst_maxval = 0.9)
    end
end


text:
# A test problem coming from LSQR.

using LinearOperators

function lstp(nrow :: Int, ncol :: Int, ndupl :: Int, npower :: Int, λ :: Real, x :: Array)

  # LSTP  generates a sparse least-squares test problem of the form
  #
  #   minimize ‖ [A ] x - [b] ‖
  #            ‖ [λI]     [0] ‖
  #
  # having a specified solution x.  The matrix A is constructed
  # in the form A = HY*D*HZ, where D is an nrow by ncol diagonal matrix,
  # and HY and HZ are Householder transformations.

  @assert(nrow >= ncol);

  # Construct two unit vectors for the Householder transformations.
  # fourpi = 4π;
  fourpi = 4 * 3.141592;    # This is the approximation used in Mike's original subroutine.
  α = fourpi / nrow;        # 4π / nrow;
  β = fourpi / ncol;        # 4π / ncol;
  hy = map(sin, [1:nrow;] * α)
  hz = map(cos, [1:ncol;] * β)

  α = norm(hy); hy /= α; HY = opHouseholder(hy);  # HY is nrow x nrow.
  β = norm(hz); hz /= β; HZ = opHouseholder(hz);  # HZ is ncol x ncol.

  # Set the diagonal matrix D containing the singular values of A.
  d = (div.(([0:ncol-1;] .+ ndupl), ndupl) * ndupl / ncol).^npower;  # Integer div!
  D = opDiagonal(nrow, ncol, d);
  A = HY * D * HZ;

  Acond = abs(d[ncol] / d[1]);

  # Compute residual vector.
  r = zeros(nrow);
  r[1:ncol]  = HZ * x ./ d;
  t = 1.0;
  for i = ncol + 1 : nrow
    j = i - ncol;
    r[i] = t * j / nrow;
    t = -t;
  end
  r = HY * r;

  # Compute right-hand side b = r + Ax.
  rnorm = norm(r);
  b = r + A * x;

  return (b, A, D, HY, HZ, Acond, rnorm)
end

function test(nrow, ncol, ndupl, npower, damp)

  x = ncol .- [1:ncol;];  # Desired solution.
  return lstp(nrow, ncol, ndupl, npower, damp, x);
end

function testall()
  damp = 0;  # Must be zero for this problem to be consistent.
  test(40, 40, 4, 1, damp);
  test(40, 40, 4, 2, damp);
  test(40, 40, 4, 3, damp);
  test(40, 40, 4, 4, damp);
end


text:
abstract type Tensor <: SymbolicObject end
struct TensorIndex <: SymbolicObject
    __pyobject__::PyCall.PyObject
end
struct TensorIndexType <: SymbolicObject
    __pyobject__::PyCall.PyObject
end
struct TensorSymmetry <: SymbolicObject
    __pyobject__::PyCall.PyObject
end
struct TensorHead <: SymbolicObject
    __pyobject__::PyCall.PyObject
end

#struct TensScalar <: Tensor
#    var::Sym;
#    expr::Sym;
#    __pyobject__::PyCall.PyObject;
#end

struct TensMul <: Tensor
    __pyobject__::PyCall.PyObject
end
struct TensAdd <: Tensor
    __pyobject__::PyCall.PyObject
end
struct IndexedTensor <: Tensor
    __pyobject__::PyCall.PyObject
end

function SymPy._convert(::Val{:MutableDenseNDimArray}, x)
    sh = x.shape
    reshape(x.tolist(), sh)
end

function sympy_type_convert(pyexp)
    if typeof(pyexp) in  [Sym,PyObject]
        cname = pyexp.__class__.__name__
        if cname == "TensMul"
            return convert(TensMul,pyexp)
        elseif cname == "TensAdd"
            return convert(TensAdd,pyexp)
        elseif cname == "Tensor"
            return convert(IndexedTensor,pyexp)
        elseif cname == "TensorSymmetry"
            return convert(TensorSymmetry,pyexp)
        elseif cname == "MutableDenseMatrix"
            return convert(Array{Sym},pyexp)
        elseif cname == "MutableDenseNDimArray"
            return convert(Array{Sym},pyexp)
        elseif cname == "ImmutableDenseNDimArray"
            if typeof(pyexp) == Sym
                return pyexp.__pyobject__.tolist()
            else
                return pyexp.tolist()
            end
        elseif cname == "Zero"
            return 0
        elseif cname == "Half"
            return 1//2
        elseif cname == "One"
            return 1
        elseif cname == "Rational"
            return convert(Rational,pyexp)
        elseif cname == "Integer"
            return convert(Int,pyexp)
        elseif cname == "Float"
            return convert(Float64,pyexp)
        else
            return pyexp
        end
    else
        return pyexp
    end
end


text:
if VERSION < v"1.0.1"
    # Workaround for JuliaLang/julia/pull/28625
    if Base.HOME_PROJECT[] !== nothing
        Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
    end
end

using Documenter, MaxPlus

makedocs(
    modules = [MaxPlus],
    sitename = "MaxPlus",
    format = Documenter.HTML(),
    pages = [
        "Introduction" => "index.md",
    ]
)

deploydocs(
    repo   = "github.com/Lecrapouille/MaxPlus.jl.git",
)


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
  Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter, GeoStats

makedocs(
  format = :html,
  sitename = "GeoStats.jl",
  authors = "Júlio Hoffimann Mendes",
  assets = ["assets/style.css"],
  pages = [
    "Home" => "index.md",
    "Manual" => [
      "Problems and solvers" => "problems_and_solvers.md",
      "Spatial data" => "spatialdata.md",
      "Domains" => "domains.md",
      "Variography" => [
        "empirical_variograms.md",
        "theoretical_variograms.md",
        "fitting_variograms.md"
      ],
      "Kriging estimators" => "estimators.md",
      "Solver comparisons" => "comparisons.md",
      "Plotting" => "plotting.md"
    ],
    "Examples" => "examples.md",
    "Contributing" => "contributing.md",
    "About" => [
      "Community" => "about/community.md",
      "License" => "about/license.md",
      "Citing" => "about/citing.md"
    ]
  ]
)

deploydocs(
  repo  = "github.com/juliohm/GeoStats.jl.git",
  target = "build",
  deps = nothing,
  make = nothing,
  julia = "1.0"
)


text:
# One of the Julia's key features is based on its type system.
# You can create and commpute with custom types without any overhead.

# How ever you need to understand the basic types of Julia to composite it.

# ## Types

# We can think about CRUD for each types.

# ### Array

# #### create an Array
# how to create an array? you can just create one. with `[]`
# one thing you need to remember is one demension array is just a vector
array = []
typeof(array)
array = [1, 2, 3]

# #### read
# read is just like the python.

# #### update
# you can modify it with index.
array[1] = 12
array
# you also asign value with same size of vector.
array[1:2] = [23, 22]
array

# But one thing you need to remember,
# unlikely list in python, one it has its type, and then it should get certain type of input.
array[2] = "test"

# ### Dictionary
# Dictionary is just like the one in python. But you need to use Dict() not {}.

dict_A = Dict()
typeof(dict_A)

dict_A["test"] = 1
dict_A
dict_A["test"]
# ### Vector

# ### Matrix

# ### DataFrame?

# ### UnitRange

# ## QnA

# ### How can I find the mathods of the type?

# You can use `methodswith` function to find methods.

using InteractiveUtils # hide

methodswith(Matrix)

# ### How can I see the type hierachy?

# You can use `subtypes` function to see subtypes.

subtypes(AbstractMatrix)

# You can define your custom function to display the whole type structure.

# Also, there are some useful packages to produce type structures.

# One of the packages is `AbstractTrees`

using AbstractTrees

AbstractTrees.children(x::Type) = subtypes(x)
print_tree(AbstractMatrix)

# or you can use `TypeTree` package.

using TypeTree

print(join(tt(AbstractMatrix), ""))

# test plots

using Plots

plot(1:10, 1:10)

savefig(joinpath(@OUTPUT, "ISL-lab-2-mpg.svg")) # hide

# \figalt{Simple plot}{ISL-lab-2-mpg.svg}


text:
# Some example runs of the clustering methods and optimization problems
 # TODO:OLD Update this file. This is the old ClustForOpt
 # saves results as jld2 file, which can be imported and analyzed by subsequent functions

using ClustForOpt
using Gurobi
env = Gurobi.Env()

 # default kmeans + centroid run
run_clust("GER","battery";n_init=3)

 #  kmeans + medoid run
run_clust("GER","battery";representation="medoid",n_init=3)

 #  kmedoids + medoid run (partitional)
run_clust("GER","battery";method="kmedoids",representation="medoid",n_init=3)

 # kmedoids + medoid run (exact)
 #QUESTION Shall we force the usage of Gurobi
run_clust("GER","battery";method="kmedoids_exact",representation="medoid",n_init=3,gurobi_env=env)

 #  hierarchical + centroid run
run_clust("GER","battery";method="hierarchical",representation="centroid",n_init=1)

 #  hierarchical + medoid run
run_clust("GER","battery";method="hierarchical",representation="medoid",n_init=1)

 # dbaclust + centroid run
 #iterations=100
 #inner_iterations=30
 #rad_sc_min=0
 #rad_sc_max=5

 #  dbaclust + centroid run (single core, for parallel runs, use parallel version)
run_clust("GER","battery";method="dbaclust",representation="centroid",n_init=3,iterations=50,rad_sc_min=0,rad_sc_max=1,inner_iterations=30)


text:
using Decimals
using Base.Test

# Decimal constructor
d = [
    Decimal(0, 2, -1),
    Decimal(0, 1, -1),
    Decimal(0, 100, -4),
    Decimal(0, 1512, -2),
    Decimal(1, 3, -2),
    Decimal(1, 4, -6),
]

@test isa(d, Array{Decimal,1})


text:
module JuliaEchoServer

using HTTP, HTTP.Messages
using JSON

const HOST = "0.0.0.0"
const PORT = 8081
mutable struct EchoServer
  host::String
  port::Int64
end
EchoServer() = EchoServer(HOST, PORT)
EchoServer(host::String) = EchoServer(host, PORT)
EchoServer(port::Int64) = EchoServer(HOST, port)

listen(self::EchoServer) = begin
  println("start listening on $(self.host):$(self.port)")
  HTTP.serve(self.host, self.port) do req::HTTP.Request
    println("receive request $req")
    headers = [
      "Content-Type" => "application/json"
    ]
    body = Dict(
      "method" => req.method,
      "target" => req.target,
      "headers" => req.headers,
      "body" => String(req.body),
    )
    return HTTP.Response(200, headers, body = JSON.json(body))
  end
end

end # module

max_stars_repo_name: nmatsui/JuliaEchoServer

max_issues_repo_name: nmatsui/JuliaEchoServer

max_forks_repo_name: nmatsui/JuliaEchoServer


text:
function pmap!(func, array)
    buffer = pmap(func, array)
    copy!(array, buffer)
end


text:
include(joinpath("dataAnalysis", "Occurrence_test.jl"))
include(joinpath("dataAnalysis", "Duration_test.jl"))
include(joinpath("dataAnalysis", "testFunctions_test.jl"))

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
function jump_objective()

    model = BilevelModel()

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    ex1 = -4x -3y

    @objective(Upper(model), Min, ex1)

    ex2 = y

    @objective(Lower(model), Min, ex2)

    @constraints(Lower(model), begin
        c1, 2x+y <= 4
        c2, x+2y <= 4
        c3, x >= 0
        c4, y >= 0
    end)

    tp = JuMP.objective_function_type(Lower(model))
    @test JuMP.objective_function(Lower(model), tp) == ex2

    tp = JuMP.objective_function_type(Upper(model))
    @test JuMP.objective_function(Upper(model), tp) == ex1

    @test JuMP.objective_sense(model) == MOI.MIN_SENSE
    @test_throws ErrorException JuMP.relative_gap(model)
    @test_throws ErrorException JuMP.dual_objective_value(model)
    @test_throws ErrorException JuMP.objective_bound(model)
    @test_throws ErrorException JuMP.set_objective(model, MOI.MAX_SENSE, x)

    @test_throws ErrorException JuMP.objective_function_type(model)
    @test_throws ErrorException JuMP.objective_function(model)
    @test_throws ErrorException JuMP.objective_function(model, MOI.SingleVariable)

    @test_throws ErrorException JuMP.optimize!(Upper(model))

end

function jump_constraints()

    model = BilevelModel()

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @constraints(Upper(model), begin
        cup, 2x+y <= 4
    end)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        c1, 2x+y <= 4
        c2, x+2y <= 4
        c3, x >= 0
        c4, y >= 0
    end)

    @test JuMP.normalized_rhs(c1) == 4.0

    @test is_valid(model, x)
    @test is_valid(Upper(model), x)
    @test is_valid(Lower(model), x) # it is in both levels

    @test is_valid(model, c1)
    @test !is_valid(Upper(model), c1)
    @test is_valid(Lower(model), c1)

    # JuMP.constraint_object(c1, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})

    JuMP.set_dual_start_value(c1, 1.2)
    JuMP.dual_start_value(c1) == 1.2

    BilevelJuMP.set_primal_upper_bound_hint(x, 1.7)
    @test BilevelJuMP.get_primal_upper_bound_hint(x) == 1.7
    BilevelJuMP.set_primal_lower_bound_hint(x, 1.8)
    @test BilevelJuMP.get_primal_lower_bound_hint(x) == 1.8

    @variable(Upper(model), 0 <= alpha <= 10, BilevelJuMP.DualOf(c1))
    @test_throws ErrorException @variable(Upper(model), 0 <= alpha <= 10, BilevelJuMP.DualOf(cup))
    @test_throws ErrorException @variable(Upper(model), 0 <= alpha <= 10, BilevelJuMP.DualOf(cup), bad_key_arg = false)
    @test_throws ErrorException @constraint(model, x+2y <= 4)

end

function jump_variables()

    model = BilevelModel()

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @constraints(Upper(model), begin
        cup, 2x+y <= 4
    end)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        c1, 2x+y <= 4
        c2, x+2y <= 4
        c3, x >= 0
        c4, y >= 0
    end)

    @variable(Upper(model), 0 <= alpha <= 10, BilevelJuMP.DualOf(c1))

    JuMP.set_start_value(alpha, 2.2)
    @test JuMP.start_value(alpha) == 2.2

    @test x === copy(x)
    @test JuMP.isequal_canonical(x, copy(x))
    @test !JuMP.isequal_canonical(x, y)

    @test JuMP.variable_type(model) == BilevelJuMP.BilevelVariableRef

    @test JuMP.is_valid(model, x)
end

function jump_objective_solver(optimizer, mode)

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @constraints(Upper(model), begin
        cx, x == 0
    end)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        cy, y == 0
    end)

    optimize!(model)

    @test JuMP.objective_sense(model) == MOI.MIN_SENSE
    # not accepted by cbc
    # @test JuMP.relative_gap(Upper(model)) ≈ 0.0 atol=1e-3
    # @test JuMP.dual_objective_value(model) ≈ 0.0 atol=1e-3
    @test JuMP.objective_value(Upper(model)) ≈ 0.0 atol=1e-3
    JuMP.objective_bound(Upper(model))
    JuMP.set_objective(Upper(model), MOI.MAX_SENSE, x)

end

function jump_display()
    atol = config.atol

    # config.bound_hint = true

    # min -4x -3y
    # s.t.
    # y = argmin_y y
    #      2x + y <= 4
    #       x +2y <= 4
    #       x     >= 0
    #           y >= 0
    #
    # sol: x = 2, y = 0
    # obj_upper = -8
    # obj_lower =  0

    model = BilevelModel()

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        c1, 2x+y <= 4
        c2, x+2y <= 4
        c3, x >= 0
        c4, y >= 0
    end)

    @test JuMP.num_constraints(model) ==
        JuMP.num_constraints(Upper(model)) + JuMP.num_constraints(Lower(model))

    display(x)
    println()
    display(c2)
    println()
    display(model)
    println()
    display(Upper(model))
    println()
    display(Lower(model))
    println()

    xx  = JuMP.variable_by_name(model, "x")
    cc2 = JuMP.constraint_by_name(model, "c2")

    # set_optimizer(model, MOIU.Model{Float64})
    # display(model)
    # println()

end

function invalid_lower_objective(optimizer, mode)
    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @test_throws ErrorException optimize!(model)
    return
end

function invalid_optimizer(optimizer, mode)
    @test_throws ErrorException BilevelModel(optimizer, mode = mode)
    return
end

function jump_display_solver(optimizer, mode)
    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    # config.bound_hint = true

    # min -4x -3y
    # s.t.
    # y = argmin_y y
    #      2x + y <= 4

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        c1, 2x+y <= 4
    end)

    display(model)
    println()
    display(Upper(model))
    println()
    display(Lower(model))
    println()

end

function jump_bounds()
    atol = config.atol

    model = BilevelModel()

    # config.bound_hint = true

    # min -4x -3y
    # s.t.
    # y = argmin_y y
    #      2x + y <= 4

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        c, 2x+y <= 4
    end)

    @variable(Upper(model), l, DualOf(c))

    for var in [x, y, l]
        @test has_lower_bound(var) == false
        set_lower_bound(var, 12)
        @test lower_bound(var) == 12
        @test has_lower_bound(var) == true
        delete_lower_bound(var)
        @test has_lower_bound(var) == false
    end

    for var in [x, y, l]
        @test has_upper_bound(var) == false
        set_upper_bound(var, 12)
        @test upper_bound(var) == 12
        @test has_upper_bound(var) == true
        delete_upper_bound(var)
        @test has_upper_bound(var) == false
    end

    for var in [x, y]
        @test is_fixed(var) == false
        fix(var, 12)
        @test fix_value(var) == 12
        @test is_fixed(var) == true
        unfix(var)
        @test is_fixed(var) == false
    end

    @test is_fixed(l) == false
    @test_throws ErrorException fix(l, 12)
    @test_throws ErrorException unfix(l)

    for var in [x, y]
        set_start_value(var, 12)
        @test start_value(var) == 12
        set_start_value(var, 13)
        @test start_value(var) == 13
    end

    for var in [x, y]
        @test is_binary(var) == false
        set_binary(var)
        @test is_binary(var) == true
        unset_binary(var)
        @test is_binary(var) == false
    end

    @test is_binary(l) == false
    @test_throws ErrorException set_binary(l)
    @test_throws ErrorException unset_binary(l)

    for var in [x, y]
        @test is_integer(var) == false
        set_integer(var)
        @test is_integer(var) == true
        unset_integer(var)
        @test is_integer(var) == false
    end

    @test is_integer(l) == false
    @test_throws ErrorException set_integer(l)
    @test_throws ErrorException unset_integer(l)

end

function jump_attributes()

    model = BilevelModel()

    # min -4x -3y
    # s.t.
    # y = argmin_y y
    #      2x + y <= 4

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        c, 2x+y <= 4
    end)

    BilevelJuMP.set_copy_names(model)
    @test BilevelJuMP.get_copy_names(model)
    BilevelJuMP.unset_copy_names(model)
    @test !BilevelJuMP.get_copy_names(model)

    BilevelJuMP.set_pass_start(model)
    @test BilevelJuMP.get_pass_start(model)
    BilevelJuMP.unset_pass_start(model)
    @test !BilevelJuMP.get_pass_start(model)

    @test isnan(JuMP.solve_time(model))
    @test isnan(BilevelJuMP.build_time(model))

    @test_throws MethodError JuMP.set_optimizer_attributes(model, "weird" => true, "strange" => "yes")
    @test_throws ErrorException JuMP.get_optimizer_attribute(model, "weird")

    return nothing
end

function jump_attributes_solver(optimizer, mode)

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @constraints(Upper(model), begin
        cx, x == 0
    end)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        cy, y == 0
    end)

    optimize!(model)

    @test JuMP.objective_sense(model) == MOI.MIN_SENSE
    # not accepted by cbc
    # @test JuMP.relative_gap(Upper(model)) ≈ 0.0 atol=1e-3
    # @test JuMP.dual_objective_value(model) ≈ 0.0 atol=1e-3
    @test JuMP.objective_value(Upper(model)) ≈ 0.0 atol=1e-3
    JuMP.objective_bound(Upper(model))
    JuMP.set_objective(Upper(model), MOI.MAX_SENSE, x)

    JuMP.solve_time(model)
    BilevelJuMP.build_time(model)

    silent_mode = JuMP.get_optimizer_attribute(model, MOI.Silent())
    JuMP.set_optimizer_attribute(model, MOI.Silent(), true)
    JuMP.unset_silent(model)
    @test !JuMP.get_optimizer_attribute(model, MOI.Silent())
    JuMP.set_silent(model)
    @test JuMP.get_optimizer_attribute(model, MOI.Silent())
    JuMP.set_optimizer_attribute(model, MOI.Silent(), silent_mode)

    JuMP.set_time_limit_sec(model, 3.0)
    @test JuMP.time_limit_sec(model) == 3.0
    JuMP.unset_time_limit_sec(model)

    @test JuMP.result_count(model) == 1
    JuMP.node_count(model)
    # TODO improve this check
    @test JuMP.simplex_iterations(model) >= 0
    @test_throws ArgumentError JuMP.barrier_iterations(model)

    @test_throws MethodError JuMP.set_optimizer_attributes(mode, "weird" => true, "strange" => "yes")

end

function mixed_mode_unit()

    # min -4x -3y
    # s.t.
    # y = argmin_y y
    #      2x + y <= 4
    #       x +2y <= 4
    #       x     >= 0
    #           y >= 0
    #
    # sol: x = 2, y = 0
    # obj_upper = -8
    # obj_lower =  0

    model = BilevelModel()

    @variable(Upper(model), x >= 0)
    @variable(Lower(model), y >= 0)

    @objective(Upper(model), Min, -4x -3y)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        c1, 2x+y <= 4
        c2, x+2y <= 4
    end)

    @test_throws ErrorException BilevelJuMP.set_mode(c1, BilevelJuMP.SOS1Mode())
    @test_throws ErrorException BilevelJuMP.set_mode(c1, BilevelJuMP.IndicatorMode())
    @test_throws ErrorException BilevelJuMP.set_mode(x, BilevelJuMP.SOS1Mode())
    @test_throws ErrorException BilevelJuMP.set_mode(x, BilevelJuMP.IndicatorMode())

    BilevelJuMP.set_mode(model, BilevelJuMP.MixedMode())

    @test_throws ErrorException BilevelJuMP.set_mode(c1, BilevelJuMP.MixedMode())
    @test_throws ErrorException BilevelJuMP.set_mode(c1, BilevelJuMP.StrongDualityMode())
    @test_throws ErrorException BilevelJuMP.set_mode(x, BilevelJuMP.MixedMode())
    @test_throws ErrorException BilevelJuMP.set_mode(x, BilevelJuMP.StrongDualityMode())

    BilevelJuMP.set_mode(x, BilevelJuMP.FortunyAmatMcCarlMode())
    BilevelJuMP.set_mode(y, BilevelJuMP.IndicatorMode())
    BilevelJuMP.set_mode(c1, BilevelJuMP.FortunyAmatMcCarlMode())
    BilevelJuMP.set_mode(c2, BilevelJuMP.IndicatorMode())

    @test typeof(BilevelJuMP.get_mode(y)) <: BilevelJuMP.IndicatorMode
    @test typeof(BilevelJuMP.get_mode(c2)) <: BilevelJuMP.IndicatorMode

    BilevelJuMP.unset_mode(x)
    BilevelJuMP.unset_mode(y)
    BilevelJuMP.unset_mode(c1)
    BilevelJuMP.unset_mode(c2)

    @test BilevelJuMP.get_mode(y) === nothing
    @test BilevelJuMP.get_mode(c2) === nothing

    return nothing
end

function variables_unit()

    model = BilevelModel()

    @variable(Upper(model), w >= 0)
    @variable(Upper(model), x >= 0)
    @variable(Lower(model), y >= 0)
    @variable(Lower(model), z >= 0)

    @test Set(JuMP.all_variables(Upper(model))) == Set([x,y,z,w])
    @test Set(JuMP.all_variables(Lower(model))) == Set([x,y,z,w])
    @test Set(JuMP.all_variables(model)) == Set([x,y,z,w])

    JuMP.delete(model, x)
    JuMP.delete(model, y)

    @test Set(JuMP.all_variables(Upper(model))) == Set([w,z])
    @test Set(JuMP.all_variables(Lower(model))) == Set([w,z])
    @test Set(JuMP.all_variables(model)) == Set([w,z])

    ex = @expression(model, w + z)
    @constraint(Upper(model), ctr, ex >= 0)

    @test 0.0 == JuMP.normalized_rhs(ctr)
    JuMP.set_normalized_rhs(ctr, 4)
    @test 4.0 == JuMP.normalized_rhs(ctr)
    JuMP.add_to_function_constant(ctr, 2)
    @test 2.0 == JuMP.normalized_rhs(ctr)

    @test 1.0 == JuMP.normalized_coefficient(ctr, w)
    JuMP.set_normalized_coefficient(ctr, w, 2.0)
    @test 2.0 == JuMP.normalized_coefficient(ctr, w)

    ex1 = BilevelAffExpr(-1.0)
    add_to_expression!(ex1, 2.0, w)
    add_to_expression!(ex1, 1.0, z)
    @constraint(Lower(model), ex1 >= 0)

    ex2 = w + z + 1
    @constraint(Lower(model), ex2 >= 0)

    ex3 = ex = w^2 + 2 * w * z + z^2 + w + z - 1
    @objective(Lower(model), Min, ex3)

    @test coefficient(ex3, w, z) == 2

    return nothing
end

function jump_no_cb()

    model = BilevelModel()

    @variable(Upper(model), x >= 0)

    @test_throws ErrorException MOI.set(model, MOI.LazyConstraintCallback(), x -> x)
    @test_throws ErrorException MOI.set(model, MOI.UserCutCallback(), x -> x)
    @test_throws ErrorException MOI.set(model, MOI.HeuristicCallback(), x -> x)

    return nothing
end

function constraint_unit()

    model = BilevelModel()

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @constraint(Upper(model), ctru, x == 0)
    @constraint(Lower(model), ctrl, y == 0)

    for (f,s) in JuMP.list_of_constraint_types(Upper(model))
        @test ctru == JuMP.all_constraints(Upper(model), f, s)[]
    end
    for (f,s) in JuMP.list_of_constraint_types(Lower(model))
        @test ctrl == JuMP.all_constraints(Lower(model), f, s)[]
        @test Set([ctrl, ctru]) == Set(JuMP.all_constraints(model, f, s))
    end
    @test JuMP.list_of_constraint_types(Lower(model)) ==
        JuMP.list_of_constraint_types(Upper(model))
    @test JuMP.list_of_constraint_types(Lower(model)) ==
        JuMP.list_of_constraint_types(model)

    for (f,s) in JuMP.list_of_constraint_types(Upper(model))
        @test JuMP.num_constraints(Upper(model), f, s) == 1
    end
    for (f,s) in JuMP.list_of_constraint_types(Lower(model))
        @test JuMP.num_constraints(Lower(model), f, s) == 1
    end
    JuMP.delete(model, ctru)
    @test isempty(JuMP.list_of_constraint_types(Upper(model)))
    @test !isempty(JuMP.list_of_constraint_types(Lower(model)))
    JuMP.delete(model, ctrl)
    @test isempty(JuMP.list_of_constraint_types(Lower(model)))
end


text:
#####
##### Middle Layer
#####

function pooling_forward(
    src::Memory{T,N},
    dims::Dims;
    kind::AbstractKind = Inference(),
    algo = Lib.dnnl_pooling_max,
    opdesc = Ref{Lib.dnnl_pooling_v2_desc_t}()
) where {T,N,M}
    dst_size = output_size(size(src), dims)
    dst_md = memorydesc(T, dst_size, dnnl_format_any())

    @apicall dnnl_pooling_v2_forward_desc_init(
        opdesc,
        kind,
        algo,
        src,
        dst_md,
        dims.strides,
        dims.kernel,
        dims.dilation,
        dims.padding,
        dims.padding,
    )

    return temp_primitive(opdesc, noattributes(), global_engine(), noforward()) do p, pd
        dst_md = query_md(pd, @query(dst))
        dst = similar(src, T, dst_size, dst_md)

        if kind === Inference()
            execute!(p, @dnnl_args src dst)
            return dst
        else
            workspace_md = query_md(pd, @query(workspace))
            workspace = similar(src, UInt8, (getbytes(workspace_md),), workspace_md)
            execute!(p, @dnnl_args src dst workspace)
            return (; dst, workspace, forward = pd)
        end
    end
end

function pooling_backward(
    _diff_dst::Memory{T,N},
    diff_src_md::MemoryDesc,
    dims::Dims;
    algo = Lib.dnnl_pooling_max,
    workspace::Memory,
    forward,
    opdesc = Ref{Lib.dnnl_pooling_v2_desc_t}()
) where {T,N}
    @apicall dnnl_pooling_v2_backward_desc_init(
        opdesc,
        algo,
        diff_src_md,
        toany(_diff_dst),
        dims.strides,
        dims.kernel,
        dims.dilation,
        dims.padding,
        dims.padding,
    )

    return temp_primitive(opdesc, noattributes(), global_engine(), forward) do p, pd
        diff_dst = maybe_reorder(pd, _diff_dst, @query(diff_dst))
        diff_src_md = query_md(pd, @query(diff_src))
        diff_src = similar(diff_dst, T, logicalsize(diff_src_md, Val(N)), diff_src_md)
        execute!(p, @dnnl_args diff_dst diff_src workspace)
        return diff_src
    end
end

#####
##### Higher Level API
#####

const Max = Lib.dnnl_pooling_max
const MeanInclude = Lib.dnnl_pooling_avg_include_padding
const MeanExclude = Lib.dnnl_pooling_avg_exclude_padding

struct Pooling{T,N}
    dims::Dims{N}
    opdesc::Base.RefValue{Lib.dnnl_pooling_v2_desc_t}
end

# Type aliases
const MaxPool{N} = Pooling{Max,N}
const MeanPool{N} = Pooling{MeanExclude,N}
const InclusiveMeanPool{N} = Pooling{MeanInclude,N}

function Pooling{T}(
    kernel::NTuple{N,Integer}; strides = kernel, dilation = 0, padding = 0
) where {T,N}
    strides = expand(Val(N), strides)
    dilation = expand(Val(N), dilation)
    padding = expand(Val(N), padding)
    opdesc = Ref{Lib.dnnl_pooling_v2_desc_t}()
    return Pooling{T,N}(Dims{N}(kernel, strides, dilation, padding), opdesc)
end

function (pool::Pooling{T})(src::AbstractArray; kw...) where {T}
    return pooling_forward(OneDNN.Memory(src), pool.dims; algo = T, pool.opdesc, kw...)
end

function ChainRulesCore.rrule(pool::Pooling{T}, _src) where {T}
    src = Memory(_src)
    src_md_any = toany(src)
    nt = pool(src; kind = Training())
    @unpack workspace, forward = nt

    function pooling_pullback(_diff_dst)
        diff_dst = Memory(_diff_dst)
        diff_src = pooling_backward(
            diff_dst, src_md_any, pool.dims; algo = T, workspace, forward, pool.opdesc
        )
        return (ChainRulesCore.NoTangent(), diff_src)
    end

    return nt.dst, pooling_pullback
end


text:
export ActorProxy, SourceProxy, ActorSourceProxy
export actor_proxy!, source_proxy!
export ProxyObservable, proxy

import Base: show

"""
    ActorProxy

Can be used as a super type for common proxy object. Automatically specifies `ValidActorProxy` trait behavior. Each `ActorProxy` must implement
its own method for `actor_proxy!(::Type, proxy, actor)` function which have to return a valid actor object.

See also: [`proxy`](@ref), [`actor_proxy!`](@ref)
"""
abstract type ActorProxy       end

"""
    SourceProxy

Can be used as a super type for common proxy object. Automatically specifies `ValidSourceProxy` trait behavior. Each `SourceProxy` must implement
its own method for `source_proxy!(::Type, proxy, source)` function which have to return a valid subscribable object.

See also: [`proxy`](@ref), [`source_proxy!`](@ref)
"""
abstract type SourceProxy      end

"""
    ActorSourceProxy

Can be used as a super type for common proxy object. Automatically specifies `ValidActorSourceProxy` trait behavior. Each `ActorSourceProxy` must implement
its own method for `source_proxy!(::Type, proxy, source)` function which have to return a valid subscribable object and also for `actor_proxy!(::Type, proxy, actor)` function which have to return a valid actor object..

See also: [`proxy`](@ref), [`actor_proxy!`](@ref), [`source_proxy!`](@ref)
"""
abstract type ActorSourceProxy end

abstract type ProxyTrait end

struct ValidActorProxy        <: ProxyTrait end
struct ValidSourceProxy       <: ProxyTrait end
struct ValidActorSourceProxy  <: ProxyTrait end
struct InvalidProxy           <: ProxyTrait end

as_proxy(::Type)                        = InvalidProxy()
as_proxy(::Type{ <: ActorProxy })       = ValidActorProxy()
as_proxy(::Type{ <: SourceProxy })      = ValidSourceProxy()
as_proxy(::Type{ <: ActorSourceProxy }) = ValidActorSourceProxy()
as_proxy(::O) where O                   = as_proxy(O)

call_actor_proxy!(::Type{L}, proxy::T,  actor::A)  where { L, T, A } = call_actor_proxy!(L, as_proxy(T), as_actor(A), proxy, actor)
call_source_proxy!(::Type{L}, proxy::T, source::S) where { L, T, S } = call_source_proxy!(L, as_proxy(T), as_subscribable(S), proxy, source)

# Check for invalid actor trait usages
call_actor_proxy!(::Type, ::InvalidProxy, ::InvalidActorTrait,   proxy, actor) = error("Type $(typeof(proxy)) is not a valid proxy type. \nConsider extending your type with one of the ActorProxy, SourceProxy or ActorSourceProxy abstract types or implement as_proxy(::Type{<:$(typeof(proxy))}).")
call_actor_proxy!(::Type, ::InvalidProxy, _,                     proxy, actor) = error("Type $(typeof(proxy)) is not a valid proxy type. \nConsider extending your type with one of the ActorProxy, SourceProxy or ActorSourceProxy abstract types or implement as_proxy(::Type{<:$(typeof(proxy))}).")
call_actor_proxy!(::Type, _,              ::InvalidActorTrait,   proxy, actor) = error("Type $(typeof(actor)) is not a valid actor type. \nConsider extending your actor with one of the abstract actor types <: (Actor{T}, NextActor{T}, ErrorActor{T}, CompletionActor{T}) or implement as_actor(::Type{<:$(typeof(actor))}).")

# Invoke valid proxy and actor
call_actor_proxy!(::Type{L}, ::ValidActorProxy,       ::BaseActorTrait{D},       proxy, actor) where { L, D } = actor_proxy!(L, proxy, actor)
call_actor_proxy!(::Type{L}, ::ValidActorProxy,       ::ErrorActorTrait{D},      proxy, actor) where { L, D } = actor_proxy!(L, proxy, actor)
call_actor_proxy!(::Type{L}, ::ValidActorProxy,       ::NextActorTrait{D},       proxy, actor) where { L, D } = actor_proxy!(L, proxy, actor)
call_actor_proxy!(::Type{L}, ::ValidActorProxy,       ::CompletionActorTrait{D}, proxy, actor) where { L, D } = actor_proxy!(L, proxy, actor)
call_actor_proxy!(::Type{L}, ::ValidSourceProxy,      ::BaseActorTrait{D},       proxy, actor) where { L, D } = actor
call_actor_proxy!(::Type{L}, ::ValidSourceProxy,      ::NextActorTrait{D},       proxy, actor) where { L, D } = actor
call_actor_proxy!(::Type{L}, ::ValidSourceProxy,      ::ErrorActorTrait{D},      proxy, actor) where { L, D } = actor
call_actor_proxy!(::Type{L}, ::ValidSourceProxy,      ::CompletionActorTrait{D}, proxy, actor) where { L, D } = actor
call_actor_proxy!(::Type{L}, ::ValidActorSourceProxy, ::BaseActorTrait{D},       proxy, actor) where { L, D } = actor_proxy!(L, proxy, actor)
call_actor_proxy!(::Type{L}, ::ValidActorSourceProxy, ::NextActorTrait{D},       proxy, actor) where { L, D } = actor_proxy!(L, proxy, actor)
call_actor_proxy!(::Type{L}, ::ValidActorSourceProxy, ::ErrorActorTrait{D},      proxy, actor) where { L, D } = actor_proxy!(L, proxy, actor)
call_actor_proxy!(::Type{L}, ::ValidActorSourceProxy, ::CompletionActorTrait{D}, proxy, actor) where { L, D } = actor_proxy!(L, proxy, actor)

# Check for invalid subscribable trait usages
call_source_proxy!(::Type, ::InvalidProxy, ::InvalidSubscribableTrait, proxy, source) = error("Type $(typeof(proxy)) is not a valid proxy type. \nConsider extending your type with one of the ActorProxy, SourceProxy or ActorSourceProxy abstract types or implement as_proxy(::Type{<:$(typeof(proxy))}).")
call_source_proxy!(::Type, ::InvalidProxy, _,                          proxy, source) = error("Type $(typeof(proxy)) is not a valid proxy type. \nConsider extending your type with one of the ActorProxy, SourceProxy or ActorSourceProxy abstract types or implement as_proxy(::Type{<:$(typeof(proxy))}).")
call_source_proxy!(::Type, _,              ::InvalidSubscribableTrait, proxy, source) = error("Type $(typeof(source)) is not a valid subscribable type. \nConsider extending your subscribable with Subscribable{T} abstract type or implement as_subscribable(::Type{<:$(typeof(source))}).")

# Invoke valid proxy and source
call_source_proxy!(::Type{L}, ::ValidActorProxy,       ::SimpleSubscribableTrait{D},    proxy, source) where { L, D } = source
call_source_proxy!(::Type{L}, ::ValidActorProxy,       ::ScheduledSubscribableTrait{D}, proxy, source) where { L, D } = source
call_source_proxy!(::Type{L}, ::ValidSourceProxy,      ::SimpleSubscribableTrait{D},    proxy, source) where { L, D } = source_proxy!(L, proxy, source)
call_source_proxy!(::Type{L}, ::ValidSourceProxy,      ::ScheduledSubscribableTrait{D}, proxy, source) where { L, D } = source_proxy!(L, proxy, source)
call_source_proxy!(::Type{L}, ::ValidActorSourceProxy, ::SimpleSubscribableTrait{D},    proxy, source) where { L, D } = source_proxy!(L, proxy, source)
call_source_proxy!(::Type{L}, ::ValidActorSourceProxy, ::ScheduledSubscribableTrait{D}, proxy, source) where { L, D } = source_proxy!(L, proxy, source)

"""
    ProxyObservable{L, S, P}(proxied_source::S, proxy::P)

An interface for proxied Observables.

See also: [`proxy`](@ref)
"""
@subscribable struct ProxyObservable{L, S, P} <: Subscribable{L}
    proxied_source :: S
    proxy          :: P
end

function on_subscribe!(observable::ProxyObservable{L}, actor) where L
    return subscribe!(observable.proxied_source, call_actor_proxy!(L, observable.proxy, actor))
end

"""
    actor_proxy!(::Type, proxy, actor)

This is function is used to wrap an actor with its proxied version given a particular proxy object. Must return another actor.
Each valid `ActorProxy` and `ActorSourceProxy` must implement its own method for `actor_proxy!` function. The first argument is the same as
the type of data of the connected proxy observable.

See also: [`proxy`](@ref), [`ActorProxy`](@ref), [`ActorSourceProxy`](@ref)
"""
actor_proxy!(L, proxy, actor) = error("You probably forgot to implement actor_proxy!(::Type, proxy::$(typeof(proxy)), actor::$(typeof(actor)))")

"""
    source_proxy!(::Type, proxy, source)

This is function is used to wrap a source with its proxied version given a particular proxy object. Must return another Observable.
Each valid `SourceProxy` and `ActorSourceProxy` must implement its own method for `source_proxy!` function. The first argument is the same as
the type of data of the connected proxy observable.

See also: [`proxy`](@ref), [`SourceProxy`](@ref), [`ActorSourceProxy`](@ref)
"""
source_proxy!(L, proxy, source) = error("You probably forgot to implement source_proxy!(::Type, proxy::$(typeof(proxy)), source::$(typeof(source)))")

Base.show(io::IO, observable::ProxyObservable{L}) where L = print(io, "ProxyObservable($L, $(observable.proxy))")

# -------------------------------- #
# Helpers                          #
# -------------------------------- #

"""
    proxy(::Type{L}, source, proxy) where L

Creation operator for the `ProxyObservable` with a given source and proxy objects.

# Example

```jldoctest
using Rocket

source = from(1:5)

struct MyCustomProxy <: ActorProxy end

struct MyCustomActor{A} <: Actor{Int}
    actor :: A
end

Rocket.on_next!(actor::MyCustomActor, data::Int) = next!(actor.actor, data ^ 2)
Rocket.on_error!(actor::MyCustomActor, err)      = error!(actor.actor, err)
Rocket.on_complete!(actor::MyCustomActor)        = complete!(actor.actor)

Rocket.actor_proxy!(::Type{Int}, proxy::MyCustomProxy, actor::A) where A = MyCustomActor{A}(actor)

proxied = proxy(Int, source, MyCustomProxy())

subscribe!(proxied, logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 4
[LogActor] Data: 9
[LogActor] Data: 16
[LogActor] Data: 25
[LogActor] Completed
```

See also: [`ProxyObservable`](@ref), [`ActorProxy`](@ref), [`SourceProxy`](@ref), [`ActorSourceProxy`](@ref)
"""
proxy(::Type{L}, source, proxy) where L = as_proxy_observable(L, call_source_proxy!(L, proxy, source), proxy)

as_proxy_observable(::Type{L}, proxied_source::S, proxy) where { L, S } = as_proxy_observable(L, as_subscribable(S), proxied_source, proxy)

as_proxy_observable(::Type{L}, ::InvalidSubscribableTrait,   proxied_source,    _)        where L           = throw(InvalidSubscribableTraitUsageError(proxied_source))
as_proxy_observable(::Type{L}, ::SimpleSubscribableTrait,    proxied_source::S, proxy::P) where { L, S, P } = ProxyObservable{L, S, P}(proxied_source, proxy)
as_proxy_observable(::Type{L}, ::ScheduledSubscribableTrait, proxied_source::S, proxy::P) where { L, S, P } = ProxyObservable{L, S, P}(proxied_source, proxy)


text:
using CompScienceMeshes
using Test

ignd(x,y) = x^2 * cos(x*y)

Iref = 0.082739596439277

p1 = [0.0, 0.0]
p2 = [1.0, 0.0]
p3 = [0.0, 1.0]

N = length(CompScienceMeshes.triangleGaussW)
J = zeros(N)
for _n in 1 : N
    u, w = trgauss(_n)
    for _g in 1 : length(w)
        _p = u[1,_g]*p1 + u[2,_g]*p2 + (1-u[1,_g]-u[2,_g])*p3
        x, y = _p[1], _p[2]
        J[_n] += w[_g] * ignd(x,y)
    end
end

E = [0.3326898477547438,0.00128868370286643,0.0003505830254662883,0.0002990068871080405,0.00026769940415425264,1.4413762136992496e-7,1.8717613253691006e-7,1.8109645370453297e-12,3.4565487801639488e-12,3.0571872387445795e-12,3.3689945069086314e-12,3.3050899511418194e-12,3.3156568461898747e-12]

for _q in 1 : N
    @test abs(J[_q]-Iref)/abs(Iref) <=  E[_q]*2.0
end


text:
using NamedDims
using OffsetArrays
using SparseArrays
using Test

@testset "get the parent array that was wrapped" begin
    for orig in ([1 2; 3 4], spzeros(2, 2))
        @test parent(NamedDimsArray(orig, (:x, :y))) === orig
        @test unname(NamedDimsArray(orig, (:x, :y))) === orig
        @test unname(orig) === orig
    end
end

@testset "function dim(array, dim)" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))
    @test dim(nda, :x) == 1
    @test dim(nda, [:y, :x]) == [2, 1]
    @test dim(nda, 2) == 2
    @test dim(parent(nda), 2) == 2
    @test_throws ArgumentError dim(nda, :nope)
end

@testset "get the named array that was wrapped" begin
    @test dimnames(NamedDimsArray([10 20; 30 40], (:x, :y))) === (:x, :y)
end

@testset "Name-asserting constructor" begin
    orig_full = NamedDimsArray(ones(3, 4, 5), (:a, :b, :c))
    @test dimnames(NamedDimsArray(orig_full, (:a, :b, :c))) == (:a, :b, :c)
    @test dimnames(NamedDimsArray(orig_full, (:a, :b, :_))) == (:a, :b, :c)
    @test_throws DimensionMismatch NamedDimsArray(orig_full, (:a, :b, :wrong))
    @test_throws DimensionMismatch NamedDimsArray(orig_full, (:c, :a, :b))

    orig_partial = NamedDimsArray(ones(3, 4, 5), (:a, :_, :c))
    @test dimnames(NamedDimsArray(orig_partial, (:a, :b, :c))) == (:a, :b, :c)
    @test dimnames(NamedDimsArray(orig_partial, (:a, :_, :c))) == (:a, :_, :c)

    @testset "deprecated: refine_names" begin
        @test dimnames(refine_names(orig_full, (:a, :b, :_))) == (:a, :b, :c)
        @test_throws DimensionMismatch NamedDimsArray(orig_full, (:a, :b, :wrong))
    end
end

@testset "getindex" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test nda[x=1, y=1] == nda[y=1, x=1] == nda[1, 1] == nda[CartesianIndex(1, 1)] == 10
    @test nda[y=end, x=end] == nda[end, end] == 40

    # Unspecified dims become slices
    @test nda[y=1] == nda[y=1, x=:] == nda[:, 1] == [10; 30]

    @test nda[CartesianIndex(1), 1] == nda[1, 1]

    @testset "Views with all dims specified are not scalars" begin
        @test view(nda, x=1, y=1) isa SubArray{Int, 0}
        @test view(nda, 1,1) == view(nda.data, 1,1)
    end

    @testset "No arguments" begin
        ndm = NamedDimsArray([10 20; 30 40], (:x, :y))
        @test_throws BoundsError ndm[]

        ndv = NamedDimsArray([1, 2, 3], (:x,))
        @test_throws BoundsError ndv[]

        nds = NamedDimsArray([4], (:x,))  # 1d 1el array (vector with one element)
        @test nds[] == 4

        nds2 = NamedDimsArray{()}(fill(4)); # 0d 1el array (kind of like a scalar)
        @test nds2[] == 4
    end

    @testset "name preservation" begin
        @test dimnames(nda[y=1]) == (:x, )
        @test dimnames(nda[y=1:1]) == (:x, :y)
    end

    # https://github.com/invenia/NamedDims.jl/issues/8
    @testset "with multiple-wildcards" begin
        nda_mw = NamedDimsArray{(:_, :_, :c)}(ones(10, 20, 30));
        @test nda_mw[c=2] == ones(10, 20)
        @test dimnames(nda_mw[c=2]) == (:_, :_)
    end

    @testset "newaxis" begin
        newaxis = [CartesianIndex()];

        @test dimnames(nda[:, newaxis, :]) == (:x, :_, :y)
        @test size(nda[:, newaxis, :]) == (2, 1, 2)

        @test dimnames(nda[1, newaxis, 1]) == (:_, )
        @test size(nda[1, newaxis, 1]) == (1, )

        @test dimnames(nda[CartesianIndex(1, 1), newaxis]) == (:_, )
        @test size(nda[CartesianIndex(1, 1), newaxis]) == (1, )
        @test nda[CartesianIndex(1, 1), newaxis][1] == nda[1, 1]
    end

    @testset "BitArray / Array{Bool}" begin
        nda = NamedDimsArray(rand(1:9, 10,11,12), (:x, :y, :z))
        ndv = NamedDimsArray(rand(1:9, 10), :x)

        # one of several dimensions:
        @test dimnames(nda[rand(10) .> 0.3, 1, :]) == (:x, :z)          # BitArray
        @test dimnames(nda[collect(rand(10) .> 0.3), 1, :]) == (:x, :z) # Array{Bool}
        @test dimnames(nda[ndv .> 3, 1, :]) == (:x, :z)                 # NamedDimsArray{...}
        @test dimnames(nda[NamedDimsArray(collect(rand(10) .> 0.3), :x), 1, :]) == (:x, :z)

        # only dim of a vector:
        @test dimnames(ndv[rand(10) .> 0.3]) == (:x,)
        @test dimnames(ndv[collect(rand(10) .> 0.3)]) == (:x,)
        @test dimnames(ndv[ndv .> 3]) == (:x,)
        @test dimnames(ndv[NamedDimsArray(collect(ndv .> 3), :x)]) == (:x,)

        # https://github.com/invenia/NamedDims.jl/issues/70
        # flattening of ndims>1 should drop names, as new dim is none of original ones:
        @test nda[rand(10,11,12) .> 0.3] isa Vector
        @test nda[collect(rand(10,11,12) .> 0.3)] isa Vector
        @test nda[nda .> 3] isa Vector
        @test nda[NamedDimsArray(collect(nda .> 3), dimnames(nda))] isa Vector
    end

    @testset "arrays of indices" begin
        nda = NamedDimsArray(rand(1:9, 3,3,3), (:x, :y, :z))
        ndv = NamedDimsArray(rand(1:9, 10), :x)

        # integers
        @test dimnames(nda[:, [1,3], :]) == (:x, :y, :z)
        @test dimnames(nda[:, [1 3; 3 1], :]) == (:x, :_, :_, :z)

        @test dimnames(ndv[[1, 3]]) == (:x,)
        @test dimnames(ndv[[1 3; 3 1]]) == (:_, :_)

        # Vector{CartesianIndex{N}}: for N>1 this makes a new dim, like nda[nda .> 3]
        @test nda[findall(iseven, nda)] isa Vector
        @test dimnames(ndv[findall(iseven, ndv)]) == (:x,)

        # ... as does a vector of integers:
        @test nda[1:end] isa Vector
        @test nda[collect(2:end)] isa Vector
        @test dimnames(ndv[1:end]) == (:x,)

        # ... or equivalent colon:
        @test nda[:] isa Vector
        @test dimnames(ndv[:]) == (:x,)
    end
end

@testset "views" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test @view(nda[y=1]) == @view(nda[y=1, x=:]) == @view(nda[:, 1]) == [10; 30]

    @testset "name preservation" begin
        @test dimnames(nda[y=1]) == (:x, )
        @test dimnames(nda[y=1:1]) == (:x, :y)
    end
end

@testset "setindex!" begin
    @testset "by name" begin
        nda = NamedDimsArray([10 20; 30 40], (:x, :y))

        nda[x=1, y=1] = 100
        @test nda == [100 20; 30 40]

        nda[x=1] .= 1000
        @test nda == [1000 1000; 30 40]
    end

    @testset "no arguments" begin
        nds = NamedDimsArray([4], (:x,))  # 1d 1el array (vector with one element)
        @test (nds[] = 2) == 2
        @test nds[] == 2

        nds2 = NamedDimsArray{()}(fill(4)); # 0d 1el array (kind of like a scalar)
        @test (nds2[] = 2) == 2
        @test nds2[] == 2
    end

    @testset "by position" begin
        nda = NamedDimsArray([10 20; 30 40], (:x, :y))

        nda[1, 1] = 100
        @test nda == [100 20; 30 40]

        nda[1, :] .= 1000
        @test nda == [1000 1000; 30 40]
    end
end

@testset "IndexStyle" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))
    @test IndexStyle(typeof(nda)) == IndexLinear()

    sparse_nda = NamedDimsArray(spzeros(4, 2), (:x, :y))
    @test IndexStyle(typeof(sparse_nda)) == IndexCartesian()
end

@testset "length/size/axes" begin
    nda = NamedDimsArray([10 20; 30 40; 50 60], (:x, :y))

    @test length(nda) == 6

    @test axes(nda) == (1:3, 1:2)
    @test axes(nda, :x) == (1:3) == axes(nda, 1)

    @test size(nda) == (3, 2)
    @test size(nda, :x) == 3 == size(nda, 1)
end

@testset "similar" begin
    nda = NamedDimsArray{(:a, :b, :c, :d)}(ones(10, 20, 30, 40))

    @testset "content" begin
        ndb = similar(nda)
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (10, 20 , 30, 40)
        @test dimnames(ndb) == (:a, :b, :c, :d)
    end
    @testset "eltype" begin
        ndb = similar(nda, Char)
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Char
        @test size(ndb) == (10, 20, 30, 40)
        @test dimnames(ndb) == (:a, :b, :c, :d)
    end
    @testset "size" begin
        ndb = similar(nda, Float64, (15, 25, 35, 45))
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (15, 25, 35, 45)
        @test dimnames(ndb) == (:a, :b, :c, :d)
    end
    @testset "dim names" begin
        ndb = similar(nda, Float64, (:w, :x, :y, :z))
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (10, 20, 30, 40)
        @test dimnames(ndb) == (:w, :x, :y, :z)
    end

    @testset "dimensions" begin
        ndb = similar(nda, Float64, (w=11, x=22))
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (11, 22)
        @test dimnames(ndb) == (:w, :x)
    end
    @testset "parent type" begin
        oa = OffsetArray(ones(10,20,30,40), -5:4, -10:9, -15:14, -20:19)
        ndb = NamedDimsArray(oa, (:a, :b, :c, :d))
        ndc = similar(ndb)
        ndd = similar(ndb, Int)
        @test parent(ndb) !== parent(ndc)
        @test eltype(ndc) == Float64
        @test size(ndc) == (10, 20, 30, 40)
        @test dimnames(ndc) == (:a, :b, :c, :d)
        @test parent(ndc) isa typeof(oa)
        @test parent(ndd) isa OffsetArray
        @test eltype(parent(ndd)) == Int
    end
    @testset "repeated names" begin
        ndr = NamedDimsArray([1 2; 3 4], (:same, :same))
        @test dimnames(similar(ndr)) == (:same, :same)
        @test dimnames(similar(ndr, Float64)) == (:same, :same)
    end
end

@testset "Strided Array Interface" begin
    x = ones(3, 5)
    nda = NamedDimsArray{(:a, :b)}(x)
    @test strides(nda) == (1, 3) == strides(x)
    @test stride(nda, :b) == 3 == stride(nda, 2) == stride(x, 2)
end

@testset "show" begin
    nda = NamedDimsArray([1 2; 3 4], (:x, :y))
    @test occursin("(:x, :y)", string(nda))
    ndv = NamedDimsArray([1, 2, 3], :x)
    @test occursin("], :x)", string(ndv))
    @test eval(Meta.parse(string(ndv))) isa NamedDimsArray{(:x,)}

    str = repr(MIME"text/plain"(), nda) # mime => 3-arg show
    @test occursin("2×2 NamedDimsArray(::", str)
    @test occursin("→ y", str)
end

const cnda = NamedDimsArray([10 20; 30 40], (:x, :y))
@testset "allocations: wrapper" begin
    @test 0 == @ballocated parent(cnda)
    @test 0 == @ballocated dimnames(cnda)

    # These tests use `@allocated` as for some reason `@ballocated` reports 1 alloc
    if VERSION >= v"1.4" # see https://github.com/invenia/NamedDims.jl/issues/115
        @test_broken 0 == @allocated NamedDimsArray(cnda, (:x, :y))
        @test_broken 0 == @allocated NamedDimsArray(cnda, (:x, :_))
    elseif VERSION >= v"1.1"
        @test 0 == @allocated NamedDimsArray(cnda, (:x, :y))
        @test 0 == @allocated NamedDimsArray(cnda, (:x, :_))
    else
        @test 0 == @allocated NamedDimsArray(cnda, (:x, :y))
        @test_broken 0 == @allocated NamedDimsArray(cnda, (:x, :_))
    end

    # indexing
    @test 0 == @ballocated cnda[1,1]
    @test 0 == @ballocated cnda[1,1] = 55

    @test 0 == @ballocated cnda[x=1, y=1]
    @test @ballocated(cnda[x=1]) == @ballocated(cnda[1, :])
    @test 0 == @ballocated cnda[x=1, y=1] = 66
end


text:
using Cbc

include(joinpath(Pkg.dir("MathProgBase"),"test","mixintprog.jl"))

mixintprogtest(CbcSolver(LogLevel=1))


text:
#
# @testset "cholesky_U_deravitive" begin
#     for n = 1:10
#         for k = 1:100
#             A = randn(MMatrix{n,n,Float64,n^2})
#             A = Symmetric(MMatrix(A * A'))  # needs to be PD
#             Ȧ = randn(MMatrix{n,n,Float64,n^2})
#             Ȧ = Symmetric(MMatrix(Ȧ + Ȧ'))  # doesn't need to be PD
#             U, U⁻¹, U̇ = cholesky_U_deravitive(A, Ȧ)
#             Ȧ_formula = U̇' * U + U' * U̇
#             @test sum(abs.(Ȧ_formula - Ȧ)) < (1.0e-10 * n^2)
#         end
#     end
# end


text:
function gamma_cdf(a::Number, x::Number)
    return gamma_inc(a,x,0)[1]
end

function gamma_ccdf(a::Number, x::Number)
    return gamma_inc(a,x,0)[2]
end

#where {T} の意味わかっていない
function gamma_cdf(a::ForwardDiff.Dual{T}, x::Number) where {T}
    y = gamma_inc(ForwardDiff.value(a),x,0)[1]
    y_a = (log(x) - digamma(a))*y-exp(a*log(x)-log(a)-loggamma(a+1))*pFq(SA[a,a], SA[a+1,a+1], -x)
    return ForwardDiff.Dual{T}(y,y_a)
end

function gamma_cdf(a::Number, x::ForwardDiff.Dual{T}) where {T}
    y = gamma_inc(a,ForwardDiff.value(x),0)[1]
    y_x = exp(-x + (a-1)*log(x) - loggamma(a))
    return ForwardDiff.Dual{T}(y,y_x)
end

function gamma_ccdf(a::ForwardDiff.Dual{T}, x::Number) where {T}
    y = gamma_inc(ForwardDiff.value(a),x,0)[2]
    y_a = -((log(x) - digamma(a))*(1-y)-exp(a*log(x)-log(a)-loggamma(a+1))*pFq(SA[a,a], SA[a+1,a+1], -x))
    return ForwardDiff.Dual{T}(y,y_a)
end

function gamma_ccdf(a::Number, x::ForwardDiff.Dual{T}) where {T}
    y = gamma_inc(a,ForwardDiff.value(x),0)[2]
    y_x = -exp(-x + (a-1)*log(x) - loggamma(a))
    return ForwardDiff.Dual{T}(y,y_x)
end

function gamma_cdf(a::ForwardDiff.Dual{T}, x::ForwardDiff.Dual{T}) where {T}
    y = gamma_cdf(ForwardDiff.value(a),ForwardDiff.value(x))
    y_a = (log(x) - digamma(a))*y-exp(a*log(x)-log(a)-loggamma(a+1))*pFq(SA[a,a], SA[a+1,a+1], -x)
    y_x = exp(-x + (a-1)*log(x) - loggamma(a))
    pa = ForwardDiff.partials(a).values
    px = ForwardDiff.partials(x).values
    ga = y_a*pa[1] + y_x*px[1]
    gx = y_a*pa[2] + y_x*px[2]
    return ForwardDiff.Dual{T}(y,ga,gx)
end

function gamma_ccdf(a::ForwardDiff.Dual{T}, x::ForwardDiff.Dual{T}) where {T}
    y = gamma_ccdf(ForwardDiff.value(a),ForwardDiff.value(x))
    y_a = -((log(x) - digamma(a))*(1-y)-exp(a*log(x)-log(a)-loggamma(a+1))*pFq(SA[a,a], SA[a+1,a+1], -x))
    y_x = -exp(-x + (a-1)*log(x) - loggamma(a))
    pa = ForwardDiff.partials(a).values
    px = ForwardDiff.partials(x).values
    ga = y_a*pa[1] + y_x*px[1]
    gx = y_a*pa[2] + y_x*px[2]
    return ForwardDiff.Dual{T}(y,ga,gx)
end

function cdf(d::Gamma, x)
    shp, scl = params(d)
    return gamma_cdf(shp, x/scl)
end

function ccdf(d::Gamma, x)
    shp, scl = params(d)
    return gamma_ccdf(shp, x/scl)
end

#integral erfc((log(x) - a)/b) dx = x erfc((log(x) - a)/b) - e^(a + b^2/4) erf((2 a + b^2 - 2 log(x))/(2 b)) + 定数
function eqcdf(d::LogNormal, x)
    mu, sigma = params(d)
    return 0.5*(1+x*erfc((log(x) - mu)/(sigma*sqrt(2)))*exp(-(mu + (sigma^2)/2)) - erf((mu + sigma^2 - log(x))/(sigma*sqrt(2))))
end

function eqcdf(d::Weibull, x)
    out = zero(x)
    if x >= zero(x)
    shp, scl = params(d)
    out = gamma_cdf(inv(shp),(max(x, 0) / scl) ^ shp)
    end
    return out
end

function eqcdf(d::Gamma, x)
    out = zero(x)
    if x >= zero(x)
        shp, scl = params(d)
        out = gamma_cdf(shp+1, x/scl)+(x/scl)*gamma_ccdf(shp, x/scl)/shp
    end
    return out
end

function eqcdf(d::Exponential, x)
    out = zero(x)
    if x >= zero(x)
        b = rate(d)
        out = b * exp(-b * x)
    end
    return out
end

function  logmean(d::LogNormal)
    mu, sigma = params(d)
    return mu + (sigma^2)/2
end

function  logmean(d::Weibull)
    shp, scl = params(d)
    return log(scl) + loggamma(1+inv(shp))
end

function  logmean(d::Gamma)
    shp, scl= params(d)
    return log(shp) + log(scl)
end

function  logmean(d::Exponential)
    b = scale(d)
    return log(b)
end

function logeqpdf(d::UnivariateDistribution,x)
    log(ccdf2(d,x))-logmean(d)
end

function eqpdf(d::UnivariateDistribution,x)
    ccdf(d,x)/mean(d)
end

max_stars_repo_name: abikoushi/IntervalCensored

max_issues_repo_name: abikoushi/IntervalCensored

max_forks_repo_name: abikoushi/IntervalCensored