text:
# * Generation of one-body Hamiltonian

function potential_matrix(V::Function, R)
    r = axes(R,1)
    R'QuasiDiagonal(V.(r))*R
end

operator(M, R) = applied(*, R, M, R')

"""
    one_body_hamiltonian(::Type{Tuple}, atom, orb)

Return the kinetic and one-body potential energy operators (as a
tuple) for the orbital `orb` of `atom`.
"""
function one_body_hamiltonian(::Type{Tuple}, atom::Atom, orb)
    R = radial_basis(atom)

    D = Derivative(axes(R,1))
    T = R'D'D*R
    T /= -2

    # Add in centrifugal part
    ℓ = getℓ(orb)
    T += potential_matrix(r -> ℓ*(ℓ+1)/(2r^2), R)

    V = potential_matrix(r -> atom.potential(orb, r), R)

    operator(T, R), operator(V, R)
end

"""
    one_body_hamiltonian(::Type{Tuple}, atom, orb)

Return the one-body energy operator for the orbital `orb` of `atom`.
"""
one_body_hamiltonian(atom, orb) = +(one_body_hamiltonian(Tuple, atom, orb)...)

# * One-body Hamiltonian operators

# ** T̂
"""
    KineticEnergyHamiltonian

The kinetic energy part of the one-body Hamiltonian, ubcluding the
centrifugal potential. It is diagonal in spin, i.e. it does not couple
orbitals of opposite spin.
"""
struct KineticEnergyHamiltonian <: OneBodyOperator end
Base.iszero(me::OrbitalMatrixElement{1,A,KineticEnergyHamiltonian,B}) where {A<:SpinOrbital{<:Orbital},B<:SpinOrbital{<:Orbital}} =
    me.a[1].m[2] != me.b[1].m[2]

Base.show(io::IO, ::KineticEnergyHamiltonian) = write(io, "T̂")

# ** V̂
"""
    PotentialEnergyHamiltonian

The potential energy part of the one-body Hamiltonian. It is diagonal
in spin, i.e. it does not couple orbitals of opposite spin.
"""
struct PotentialEnergyHamiltonian <: OneBodyOperator end
Base.iszero(me::OrbitalMatrixElement{1,A,PotentialEnergyHamiltonian,B}) where {A<:SpinOrbital{<:Orbital},B<:SpinOrbital{<:Orbital}} =
    me.a[1].m[2] != me.b[1].m[2]

Base.show(io::IO, ::PotentialEnergyHamiltonian) = write(io, "V̂")

# * AtomicOneBodyHamiltonian types

abstract type AbstractAtomicOneBodyHamiltonian end
Base.iszero(::AbstractAtomicOneBodyHamiltonian) = false

# ** ZeroAtomicOneBodyHamiltonian

struct ZeroAtomicOneBodyHamiltonian{M,N} <: AbstractAtomicOneBodyHamiltonian
    a::Tuple{M,N}
end

Base.axes(ĥ::ZeroAtomicOneBodyHamiltonian) = ĥ.a
Base.axes(ĥ::ZeroAtomicOneBodyHamiltonian,i) = ĥ.a[i]

Base.iszero(::ZeroAtomicOneBodyHamiltonian) = true

Base.show(io::IO, ĥ::ZeroAtomicOneBodyHamiltonian) =
    write(io, "0")

# ** AtomicOneBodyHamiltonian
"""
    AtomicOneBodyHamiltonian(op, orbital)

Structure holding a one-body energy operator `op` acting on its
associated `orbital`.
"""
struct AtomicOneBodyHamiltonian{LT,O} <: AbstractAtomicOneBodyHamiltonian
    op::LT
    orbital::O
end

"""
    AtomicOneBodyHamiltonian(atom, orb)

Create the one-body Hamiltonian corresponding to the orbital `orb` of
`atom`.
"""
AtomicOneBodyHamiltonian(atom::Atom, orb::AbstractOrbital) =
    AtomicOneBodyHamiltonian(one_body_hamiltonian(atom, orb), orb)

Base.axes(ĥ::AtomicOneBodyHamiltonian, args...) =
    axes(ĥ.op, args...)

Base.zero(ĥ::AtomicOneBodyHamiltonian) =
    ZeroAtomicOneBodyHamiltonian(axes(ĥ))

SCF.update!(::AtomicOneBodyHamiltonian) = nothing
SCF.update!(::AtomicOneBodyHamiltonian, ::Atom) = nothing

matrix(ĥ::AtomicOneBodyHamiltonian) = matrix(ĥ.op)

"""
    ĥ ⋆ ϕ

Return the lazy multiplicative action of the
[`AtomicOneBodyHamiltonian`](@ref) `ĥ` on the radial orbital
coefficient vector `ϕ`.
"""
LazyArrays.:(⋆)(ĥ::AtomicOneBodyHamiltonian, ϕ::AbstractVector) =
    ĥ.op.args[2] ⋆ ϕ

"""
    materialize!(::MulAdd{<:Any, <:Any, <:Any, T, <:AtomicOneBodyHamiltonian, Source, Dest})

Materialize the lazy multiplication–addition of the type `y ← α*H*x +
β*y` where `H` is a [`AtomicOneBodyHamiltonian`](@ref) and `x` and `y`
are [`RadialOrbital`](@ref)s.
"""
LazyArrays.materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T, <:AtomicOneBodyHamiltonian, Source, Dest}) where {T,Source,Dest} =
    materialize!(MulAdd(ma.α, ma.A.op.args[2], ma.B.args[2],
                        ma.β, ma.C.args[2]))

Base.show(io::IO, ĥ::AtomicOneBodyHamiltonian) =
    write(io, "ĥ($(ĥ.orbital))")

# * Diagonalization of one-body Hamiltonian

"""
    diagonalize_one_body(H, nev; method=:arnoldi_shift_invert, tol=1e-10, σ=-1)

Diagonalize the one-body Hamiltonian `H` and find the `nev` lowest
eigenpairs, using the specified diagonalization `method`; valid
choices are

- `:arnoldi` which performs the standard Krylov iteration looking for
  the eigenvalues with smallest real values,

- `:arnoldi_shift_invert` which performs the Krylov iteration but with
  the shifted and inverted matrix `(H - I*σ)⁻¹` looking for the
  eigenvalues with _largest_ real values,

- `:eigen` which uses Julia's built-in eigensolver.

`tol` sets the Krylov tolerance.
"""
function diagonalize_one_body(H::RadialOperator, nev::Int;
                              method::Symbol=:arnoldi_shift_invert, tol=1e-10, σ=-1,
                              verbosity=0, io=stdout)
    Hm = matrix(H)
    verbosity > 2 && println(io, "Diagonalizing via $(method)")
    if method == :arnoldi || method == :arnoldi_shift_invert
        A,target = method == :arnoldi ? (Hm,SR()) : (ShiftInvert(Hm, σ),LR())

        schur,history = partialschur(A, nev=nev, tol=tol, which=target)
        verbosity > 3 && println(io, history)

        θ = diag(schur.R)
        λ = method == :arnoldi ? θ : σ .+ inv.(θ)
        length(λ) < nev &&
            error("Could not converge the requested orbitals: $(history)")
        λ,schur.Q
    elseif method == :eigen
        ee = eigen(Hm)
        ee.values[1:nev],ee.vectors[:,1:nev]
    else
        throw(ArgumentError("Unknown diagonalization method $(method)"))
    end
end


text:
using GeoPhyInv
using Test
using Random

mgrid = [range(0.0, stop=10.,step=0.1), range(0.0, stop=30.,step=0.2)];
ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(10));

tgrid=range(0, stop=1.0, step=0.1);

data=Data(tgrid, ageom, [:P,:Vz]);

Random.randn!(data[3][:P]);

ext: jl


text:
using EIIndex
using Test
using LightGraphs

N = 1000; L = 10000; G = 3; S = 1000
g = erdos_renyi(N, L)
groups = rand(1:G, N)
cug = cugtest_edges(g, x -> ei_global(x, groups), S)
plot(cug, xlabel="EI-Index")

# @testset "EIIndex.jl" begin
#     # Write your tests here.
# end

max_stars_repo_name: jfb-h/EIIndex

max_issues_repo_name: jfb-h/EIIndex

max_forks_repo_name: jfb-h/EIIndex


text:
export softmax_crossentropy

doc"""
    softmax_crossentropy(p, x)

Cross-entropy function between p and ``softmax(x)``.
```math
f(x) = -\sum_{x} p(x) \log q(x)
```
where ``q = softmax(x)``

* p: Var of Vector{Int} or Matrix{Float}
* q: Var of Matrix{Float}

```julia
p = Var(rand(0:10,5))
x = Var(rand(Float32,10,5))
y = softmax_crossentropy(p, x)
```
"""
function softmax_crossentropy(p::Var, q::Var)
    logq = logsoftmax(q.data)
    ydata = softmax_crossentropy(p.data, logq)
    Var(ydata, ∇softmax_crossentropy!, (p,q,logq))
end
softmax_crossentropy(p::Node, q::Node) = Node(softmax_crossentropy, (p,q))

function softmax_crossentropy(p::Vector{Int}, logq::Matrix{T}) where T
    length(p) == size(logq,2) || throw("Length unmatch.")
    y = zeros(T, length(p))
    @inbounds for i = 1:length(y)
        if p[i] > 0
            y[i] = -logq[p[i],i]
        end
    end
    y
end

function softmax_crossentropy(p::Matrix{Int}, logq::Matrix{T}) where T
    size(p) == size(logq) || throw("Size mismatch.")
    y = Array{T}(size(p,2))
    @inbounds for j = 1:size(p,2)
        s = T(0)
        for i = 1:size(p,1)
            s += -p[i,j] * logq[i,j]
        end
        y[j] = s
    end
    y
end

function ∇softmax_crossentropy!(y::Var, p::Var, q::Var, logq)
    isnothing(q.grad) && return
    ∇softmax_crossentropy!(y.grad, p.data, q.grad, logq)
end

function ∇softmax_crossentropy!(gy::Vector{T}, p::Vector{Int}, gq::Matrix{T}, logq::Matrix{T}) where T
    @inbounds for j = 1:length(p)
        p[j] <= 0 && continue
        for i = 1:size(logq,1)
            delta = i == p[j] ? T(1) : T(0)
            gq[i,j] += gy[j] * (exp(logq[i,j]) - delta)
        end
    end
end

function ∇softmax_crossentropy!(gy::Vector{T}, p::Matrix{T}, gq::Matrix{T}, logq::Matrix{T}) where T
    @inbounds for j = 1:size(p,2)
        for i = 1:size(logq,1)
            gq[i,j] += gy[j] * (exp(logq[i,j]) - p[i,j])
        end
    end
end

@generated function softmax_crossentropy(p::CuVector{Cint}, logq::CuMatrix{T}) where T
    Ct = cstring(T)
    k = Kernel("""
    __global__ void softmax_crossentropy($Ct *y, int *p, $Ct *logq, int size1, int size2) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= size2) return;
        int qi = (p[idx]-1) + size1 * idx;
        y[idx] = p[idx] > 0 ? -logq[qi] : 0;
    }""")
    quote
        length(p) == size(logq,2) || throw("Length unmatch.")
        y = CuArray{T}(length(p))
        gdims, bdims = cudims(length(y))
        $k(gdims, bdims, pointer(y), pointer(p), pointer(logq), size(logq,1), size(logq,2))
        y
    end
end

@generated function softmax_crossentropy(p::CuMatrix{T}, logq::CuMatrix{T}) where T
    Ct = cstring(T)
    k = Kernel("""
    __global__ void softmax_crossentropy($Ct *y, $Ct *p, $Ct *logq, int length) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx < length) {
            y[idx] = -p[idx] * logq[idx];
        }
    }""")
    quote
        size(p) == size(logq) || throw("Length unmatch.")
        y = similar(p)
        gdims, bdims = cudims(length(y))
        $k(gdims, bdims, pointer(y), pointer(p), pointer(logq), length(y))
        vec(sum(y,1))
    end
end

@generated function ∇softmax_crossentropy!(gy::CuVector{T}, p::CuVector{Cint}, gq::CuMatrix{T}, logq::CuMatrix{T}) where T
    Ct = cstring(T)
    k = Kernel("""
    __global__ void softmax_crossentropy_grad($Ct *gy, int *p, Array<$Ct,2> gq, Array<$Ct,2> logq) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= logq.length()) return;

        int ndidxs[2];
        logq.ndindex(ndidxs, idx);
        int i = ndidxs[0];
        int j = ndidxs[1];
        if (p[j] > 0) {
            $Ct delta = (i == p[j]-1) ? 1 : 0;
            gq(i,j) += gy[j] * (exp(logq(i,j)) - delta);
        }
    }""")
    quote
        gdims, bdims = cudims(length(logq))
        $k(gdims, bdims, pointer(gy), pointer(p), gq, logq)
    end
end

@generated function ∇softmax_crossentropy!(gy::CuVector{T}, p::CuMatrix{T}, gq::CuMatrix{T}, logq::CuMatrix{T}) where T
    Ct = cstring(T)
    k = Kernel("""
    __global__ void softmax_crossentropy_grad($Ct *gy, $Ct *p, $Ct *gq, Array<$Ct,2> logq) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= logq.length()) return;

        int ndidxs[2];
        logq.ndindex(ndidxs, idx);
        int i = ndidxs[0];
        int j = ndidxs[1];
        gq[idx] += gy[j] * (exp(logq[idx]) - p[idx]);
    }""")
    quote
        gdims, bdims = cudims(length(logq))
        $k(gdims, bdims, rawpointer(gy), rawpointer(p), rawpointer(gq), logq)
    end
end


text:
const APL = AbstractPolynomialLike

function _expectation(μ::Measure{T, BT}, p::APL, f) where {T, BT}
    i = 1
    s = 0
    basis = MB.basis_covering_monomials(BT, monomials(p))
    for (c, m) in zip(coefficients(p, basis), basis)
        while i <= length(μ.basis) && m != μ.basis[i]
            i += 1
        end
        if i > length(μ.basis)
            error("The base function $m has a non-zero multiplier $c in $p,
                  but its expectation is not known in $μ")
        end
        s += f(μ.values[i], c)
        i += 1
    end
    s
end

"""
    MultivariateMoments.expectation(μ::AbstractMeasureLike, p::AbstractPolynomialLike)
    MultivariateMoments.expectation(p::AbstractPolynomialLike, μ::AbstractMeasureLike)

Computes the expectation ``\\mathbb{E}_{\\mu}[p]``.
"""
function expectation end

expectation(μ::Measure, p::APL) = _expectation(μ, p, (*))
expectation(p::APL, μ::Measure) = _expectation(μ, p, (a, b) -> b * a) # a and b may be noncommutative

"""
    dot(μ::AbstractMeasureLike, p::AbstractPolynomialLike)
    dot(p::AbstractPolynomialLike, μ::AbstractMeasureLike)

See [`expectation`](@ref)
"""
LinearAlgebra.dot(μ::AbstractMeasureLike, p::APL) = expectation(μ, p)
LinearAlgebra.dot(p::APL, μ::AbstractMeasureLike) = expectation(p, μ)


text:
#=
Definition of measures of cross-correlation and power-law correlation
=#
using MIDI, MusicManipulations
using Statistics, Random
using InformationMeasures
using DCCA

"""
    mutualinformation(x, y [, N]) → MI, MICI, center

Compute MI and it's confidence interval w.r.t the null hypothesis
(i.e no correlations), via bootstrap procedure (total of `N` times),
between `x` and `y`.

The confidence interval for the mutual information corresponding to
the null-hypothesis of no correlations. If the value of MI lies outside this
interval, then `x` and `y` are co-dependant.

Return MI, it's confidence interval and the center of the interval's distribution.
"""
function mutualinformation(MTD::Vector, vel::Vector, N = 10000)
    MTD = copy(MTD); vel = copy(vel)
    MI = round(get_mutual_information(MTD,vel), digits = 2)
    # direct estimation of error bars for MI is too complex.
    # test the null-hypothesis instead.
    # bootstraping procedure to get 95% confidence interval for MI under null-hypothesis.
    bootstrap = zeros(N)
    for i in 1:N
        bootstrap[i] = get_mutual_information(shuffle!(MTD), shuffle!(vel))
    end
    sort!(bootstrap)
    MiCI = (bootstrap[max(N÷40, 1)], bootstrap[N - N÷40])
    center = bootstrap[N÷2]
    return MI, MiCI, center
end

function normalized_mi(MTD::Vector, vel::Vector, N = 10000)
    mi, mici, cen = mutualinformation(MTD, vel, N)
    normalized_mi(mi, mici, cen)
end
function normalized_mi(mi::Real, mici::Tuple, cen::Real)
    mic = mi-cen
    if mic < 0
        return mic/(cen - mici[1])
    else
        return mic/(mici[2] - cen)
    end
end

"""
    pearson(x, y) → ρ, CI
Compute ρ, and it's 95% confidence interval between `x` and `y`.

Notice that `CI` is a 2-tuple (for low and upper bounds).
"""
function pearson(MTD::Vector,vel::Vector)
    r = cor(MTD,vel)
    #computing z transform and standard deviation for r'S error bars.
    #going to the z_transform space makes the distribution symmetric
    z_transform = 0.5*log((1+r)/(1-r))
    sigma = sqrt(1/(length(vel)-3))
    z_critical = 1.96*sigma
    #going back to the r space.
    rCI = (tanh(z_transform - z_critical), tanh(z_transform + z_critical))
    return r, rCI
end

"""
    bin_array(x, bin_size) → binned_array
Bin 'x' with a bin width of 'bin_size'.
"""
function bin_array(x::Array{Int,1}, bin_size::Int)
    new_bins = collect(unique(sort(x))[1]:bin_size:unique(sort(x))[end])
    binned_array = Int[]
    for value in x
        for index in 1:length(new_bins)-1
            if value >= new_bins[index] && value < new_bins[index+1]
                push!(binned_array,new_bins[index])
            end
            if value >= new_bins[end]
                push!(binned_array, new_bins[end])
                break
            end
        end
    end
    return binned_array
end

"""
    correlation_intervals(c)
Given container `c`, integrate with DrWatson's `produce_or_load` and
make a file that contains the Pearson and MI data.
"""
function correlation_intervals(c)
    mtds, vels = get_dataset(c)

    if get(c, :absolute, false)
        for i in 1:length(mtds)
            mtds[i] = abs.(mtds[i])
        end
    end

    # Pearson
    ρs = Float64[]
    ρus = Float64[]
    ρds = Float64[]
    # Mutual information stuff
    MIs = Float64[]
    centers = Float64[]
    MIus = Float64[]
    MIds = Float64[]

    for (m, v) in zip(mtds, vels)
        ρ, (d, u) = pearson(m, v)
        push!(ρs, ρ); push!(ρus, u - ρ); push!(ρds, ρ - d)

        MI, MICI, center = mutualinformation(m, v, 10000)
        push!(MIs, MI)
        push!(MIds, center - MICI[1])
        push!(MIus, MICI[2] - center)
        push!(centers, center)
    end
    return @dict( ρs, ρus, ρds, MIs, centers, MIus, MIds )
end

"""
    get_DCCA(x,y,start,stop,pts)
compute the rho_DCCA coefficients for the given window sizes, as well as a c
onfidence interval for the null hypothesis. Return two tuples :
- (points::Vector, values::Vector) which is the results of the DCCA analyses
  for the different time-scales contained in windows
- (l::Float64, u::Float64) the lower and upper bound for the 95% confidence
  interval
"""
function get_DCCA(x::Array{Float64,1},y::Array{Float64,1})
    windows, values = rhoDCCA(x,y)
    l, u = rhoDCCA_CI(x,y)
    return (windows, values), (l,u)
end

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
@enum(BRANCH_WEIGHT_STRATEGY, BW_INFEASIBLE, BW_INTERVAL_BRANCH,
                              BW_INTERVAL_BRANCH_REV, BW_INTERVAL_LP,
                              BW_INTERVAL_LP_REV)

Base.@kwdef mutable struct BranchOracle{T<:Real}
    strategy::BRANCH_WEIGHT_STRATEGY = BW_INTERVAL_LP
    𝛹m::Vector{T}     = T[]
    𝛹p::Vector{T}      = T[]
    δm::Vector{T}     = T[]
    δp::Vector{T}      = T[]
    ηm::Vector{T}     = T[]
    ηp::Vector{T}        = T[]
    score::Vector{T}         = T[]
    score_indx::Vector{Int}  = Int[]
    μ1::T                  = 0.1
    μ2::T                  = 1.3
    μ3::T                  = 0.8
    β::T                   = 0.05
    μ_score::T             = 0.15
    η_rel::T               = 3
end
function BranchOracle(T::Type{<:AbstractFloat}, n)
    BranchOracle{T}(𝛹m = ones(T,n),
                    𝛹p = ones(T,n),
                    δm = zeros(T,n),
                    δp = zeros(T,n),
                    ηm = zeros(T,n),
                    ηp = zeros(T,n)
                    )
end

function _variable_infeasibility(m, d::BranchOracle{T}, i) where T<:Real
    tsum = zero(T); tmin = Inf; tmax = -Inf
    for j in m._constr_dep[i]
        constraint_infeasiblity = m._constr_inf[j]
        tsum += constraint_infeasiblity
        if constraint_infeasiblity > tmax
            tmax = constraint_infeasiblity
        end
        if constraint_infeasiblity < tmin
            tmin = constraint_infeasiblity
        end
    end
    return d.μ1*tsum + d.μ2*tmin + d.μ3*tmax
end

function _set_δ!(::Val{BW_INFEASIBLE}, m::GlobalOptimizer, d::BranchOracle, i) where T<:Real
    _set_constraint_infeasibility!(m, i) # sets m._constr_inf
    v = _variable_infeasibility(m, d, i)
    m._variable_inf[i] = v
    d.δm[i] = v
    d.δp[i] = v
    return
end

function _set_δ!(::Val{BW_INTERVAL_BRANCH}, m::GlobalOptimizer, d::BranchOracle, i) where {T<:Real}
    l = _lower_bound(BranchVar, m, i)
    u = _upper_bound(BranchVar, m, i)
    d.δm[i] = isfinite(l) ? (xb - l) : _variable_infeasibility(m, d, i)
    d.δp[i]  = isfinite(u) ? (u - xb) : _variable_infeasibility(m, d, i)
    return
end

function _set_δ!(::Val{BW_INTERVAL_BRANCH_REV}, m::GlobalOptimizer, d::BranchOracle, i) where {T<:Real}
    l = _lower_bound(BranchVar, m, i)
    u = _upper_bound(BranchVar, m, i)
    d.δm[i] = isfinite(l) ? (u - xb) : _variable_infeasibility(m, d, i)
    d.δp[i] = isfinite(u) ? (xb - l) : _variable_infeasibility(m, d, i)
    return
end

function _set_δ!(::Val{BW_INTERVAL_LP}, m::GlobalOptimizer, d::BranchOracle, i) where {T<:Real}
    l = _lower_bound(BranchVar, m, i)
    u = _upper_bound(BranchVar, m, i)
    ρ = d.β*(u - l)
    xlp_adj = max(min(xlp, u - ρ), l + ρ)
    d.δm[i] = isfinite(l) ? (xlp_adj - l) : _variable_infeasibility(m, d, i)
    d.δp[i] = isfinite(u) ? (u - xlp_adj) : _variable_infeasibility(m, d, i)
    return
end

function _set_δ!(::Val{BW_INTERVAL_LP_REV}, m::GlobalOptimizer, d::BranchOracle, i) where {T<:Real}
    l = _lower_bound(BranchVar, m, i)
    u = _upper_bound(BranchVar, m, i)
    ρ = d.β*(u - l)
    xlp_adj = max(min(xlp, u - ρ), l + ρ)
    d.δm[i] = isfinite(l) ? (u - xlp_adj) : _variable_infeasibility(m, d, i)
    d.δp[i] = isfinite(u) ? (xlp_adj - l) : _variable_infeasibility(m, d, i)
    return
end

function _set_δ!(m::GlobalOptimizer, i)
    s = m._branch_oracle.strategy
    if s == BW_INTERVAL_LP
        _set_δ!(Val(BW_INTERVAL_LP), m, m._branch_oracle, i)
    elseif s == BW_INTERVAL_LP_REV
        _set_δ!(Val(BW_INTERVAL_LP_REV), m, m._branch_oracle, i)
    elseif s == BW_INTERVAL_BRANCH
        _set_δ!(Val(BW_INTERVAL_BRANCH), m, m._branch_oracle, i)
    elseif s == BW_INTERVAL_BRANCH_REV
        _set_δ!(Val(BW_INTERVAL_BRANCH_REV), m, m._branch_oracle, i)
    elseif s == BW_INFEASIBLE
        _set_δ!(Val(BW_INFEASIBLE), m, m._branch_oracle, i)
    end
    return nothing
end

function _pseudocost_branch!()
end

#=
Algorithm 3, Section 2.6 of Achterberg, Tobias, Thorsten Koch, and Alexander
Martin. "Branching rules revisited." Operations Research Letters 33.1 (2005):
42-54.
=#
function _reliability_branch!()
    MOI.set(m, MOI.InterationNumber(), k)
    _score!(m)
    sort_perm!(score_indx, score, rev = true)

    max_score = maximum(m._branch_score)
    no_improve_count = 0
    for i in unreliable
        if no_improve_count > no_improve_limit
            break
        end

        # solve subproblem lower
        obj_minus = MOI.get(m, MOI.ObjectiveValue()) - obj_val
        # solve subproblem upper
        obj_plus = MOI.get(m, MOI.ObjectiveValue()) - obj_val

        _set_pseudo_cost!(obj_plus, obj_minus, i)

        new_score = _score(obj_plus, obj_minus, d.μ_score)
        if new_score >= max_score
            no_improve_count += 1
        else
            m._branch_score[i] = new_score
            no_improve_count = 0
        end
    end
    MOI.set(m, MOI.InterationNumber(), k_orig)
end

_score(x::T, y, μ) where T<:Real = (one(T) - μ)*min(x, y) + max(x, y)
_score(d::BranchOracle{T}, i) where T<:Real = _score(𝛹m[i]*δm[i], 𝛹p[i]*δp[i], d.μ_score)
function _score!(m)
    map(i -> _score(m._branch_oracle, i), m._branch_score, m._branch_variable)
    return
end

function _select_branch_variable!(m)
    _score!(m)
    return argmax(m._branch_score)
end


text:
# module CheckArgument

"""
	参数检查 概率范围检查
	概率需要在0到1之间，否则抛出参数错误
"""
function check_probabilities(p::AbstractFloat)
    if p > 1 || p < 0
throw(ArgumentError("p should be between 0-1"))
end
end

"""
	参数检查 非负性检查
	如果检查项目不是非负数，抛出参数错误
"""
function check_is_nonnegative(nums::Real...)
    for n in nums
        if n < 0
            throw(ArgumentError("$n should be nonnegative"))
        end
    end
end

        """
    检测range是否为空
"""
function check_not_empty(ranges::AbstractRange...)
    for r in ranges
        if length(r) <= 0
            throw(ArgumentError("range should not be empty"))
        end
    end
end
# export check_is_nonnegative,check_probabilities
# end

max_stars_repo_name: mudssky/juliaPackages

max_issues_repo_name: mudssky/juliaPackages

max_forks_repo_name: mudssky/juliaPackages


text:
"""
    SamplingContext(rng, sampler, context)

Create a context that allows you to sample parameters with the `sampler` when running the model.
The `context` determines how the returned log density is computed when running the model.

See also: [`DefaultContext`](@ref), [`LikelihoodContext`](@ref), [`PriorContext`](@ref)
"""
struct SamplingContext{S<:AbstractSampler,C<:AbstractContext,R} <: AbstractContext
    rng::R
    sampler::S
    context::C
end

"""
    struct DefaultContext <: AbstractContext end

The `DefaultContext` is used by default to compute log the joint probability of the data
and parameters when running the model.
"""
struct DefaultContext <: AbstractContext end

"""
    struct PriorContext{Tvars} <: AbstractContext
        vars::Tvars
    end

The `PriorContext` enables the computation of the log prior of the parameters `vars` when
running the model.
"""
struct PriorContext{Tvars} <: AbstractContext
    vars::Tvars
end
PriorContext() = PriorContext(nothing)

"""
    struct LikelihoodContext{Tvars} <: AbstractContext
        vars::Tvars
    end

The `LikelihoodContext` enables the computation of the log likelihood of the parameters when
running the model. `vars` can be used to evaluate the log likelihood for specific values
of the model's parameters. If `vars` is `nothing`, the parameter values inside the `VarInfo` will be used by default.
"""
struct LikelihoodContext{Tvars} <: AbstractContext
    vars::Tvars
end
LikelihoodContext() = LikelihoodContext(nothing)

"""
    struct MiniBatchContext{Tctx, T} <: AbstractContext
        context::Tctx
        loglike_scalar::T
    end

The `MiniBatchContext` enables the computation of
`log(prior) + s * log(likelihood of a batch)` when running the model, where `s` is the
`loglike_scalar` field, typically equal to `the number of data points / batch size`.
This is useful in batch-based stochastic gradient descent algorithms to be optimizing
`log(prior) + log(likelihood of all the data points)` in the expectation.
"""
struct MiniBatchContext{Tctx,T} <: AbstractContext
    context::Tctx
    loglike_scalar::T
end
function MiniBatchContext(context=DefaultContext(); batch_size, npoints)
    return MiniBatchContext(context, npoints / batch_size)
end

"""
    PrefixContext{Prefix}(context)

Create a context that allows you to use the wrapped `context` when running the model and
adds the `Prefix` to all parameters.

This context is useful in nested models to ensure that the names of the parameters are
unique.

See also: [`@submodel`](@ref)
"""
struct PrefixContext{Prefix,C} <: AbstractContext
    context::C
end
function PrefixContext{Prefix}(context::AbstractContext) where {Prefix}
    return PrefixContext{Prefix,typeof(context)}(context)
end

const PREFIX_SEPARATOR = Symbol(".")

function PrefixContext{PrefixInner}(
    context::PrefixContext{PrefixOuter}
) where {PrefixInner,PrefixOuter}
    if @generated
        :(PrefixContext{$(QuoteNode(Symbol(PrefixOuter, PREFIX_SEPARATOR, PrefixInner)))}(
            context.context
        ))
    else
        PrefixContext{Symbol(PrefixOuter, PREFIX_SEPARATOR, PrefixInner)}(context.context)
    end
end

function prefix(::PrefixContext{Prefix}, vn::VarName{Sym}) where {Prefix,Sym}
    if @generated
        return :(VarName{$(QuoteNode(Symbol(Prefix, PREFIX_SEPARATOR, Sym)))}(vn.indexing))
    else
        VarName{Symbol(Prefix, PREFIX_SEPARATOR, Sym)}(vn.indexing)
    end
end


text:
#
# type EasyMessage
#   pts::Array{Float64,2}
#   bws::Array{Float64,1}
# end

mutable struct NBPMessage <: Singleton
  p::Dict{Int,EasyMessage}
end

mutable struct PotProd
    Xi::Int
    prev::Array{Float64,2}
    product::Array{Float64,2}
    potentials::Array{BallTreeDensity,1}
    potentialfac::Vector{AbstractString}
end
mutable struct CliqGibbsMC
    prods::Array{PotProd,1}
    lbls::Vector{Symbol}
    CliqGibbsMC() = new()
    CliqGibbsMC(a,b) = new(a,b)
end
mutable struct DebugCliqMCMC
    mcmc::Union{Void, Array{CliqGibbsMC,1}}
    outmsg::NBPMessage
    outmsglbls::Dict{Symbol, Int}
    priorprods::Vector{CliqGibbsMC} #Union{Void, Dict{Symbol, Vector{EasyMessage}}}
    DebugCliqMCMC() = new()
    DebugCliqMCMC(a,b,c,d) = new(a,b,c,d)
end

mutable struct UpReturnBPType
    upMsgs::NBPMessage
    dbgUp::DebugCliqMCMC
    IDvals::Dict{Int, EasyMessage} #Array{Float64,2}
end

mutable struct DownReturnBPType
    dwnMsg::NBPMessage
    dbgDwn::DebugCliqMCMC
    IDvals::Dict{Int,EasyMessage} #Array{Float64,2}
end

mutable struct ExploreTreeType
  fg::FactorGraph
  bt::BayesTree
  cliq::Graphs.ExVertex
  prnt::Any
  sendmsgs::Array{NBPMessage,1}
end

mutable struct MsgPassType
  fg::FactorGraph
  cliq::Graphs.ExVertex
  vid::Int
  msgs::Array{NBPMessage,1}
  N::Int
end

# type UpMsgPassType
#   fg::FactorGraph
#   cliq::Graphs.ExVertex
#   upMsgs::Array{NBPMessage,1}
#   N::Int
# end

#global pidx
pidx = 1
pidl = 1
pidA = 1
thxl = nprocs() > 4 ? floor(Int,nprocs()*0.333333333) : 1

# upploc to control processes done local to this machine and separated from other
# highly loaded processes. upploc() should be used for dispatching short and burst
# of high bottle neck computations. Use upp2() for general multiple dispatch.
function upploc()
    global pidl, thxl
    N = nprocs()
    pidl = (thxl > 1 ? ((pidl < thxl && pidl != 0 && thxl > 2) ? (pidl+1)%(thxl+1) : 2) : (N == 1 ? 1 : 2))
    return pidl
end

# upp2() may refer to processes on a different machine. Do not associate upploc()
# with processes on a separate computer -- this will be become more complicated when
# remote processes desire their own short fast 'local' burst computations. We
# don't want all the data traveling back and forth for shorter jobs
function upp2()
  global pidx, thxl
  N = nprocs()
  pidx = (N > 1 ? ((pidx < N && pidx != 0) ? (pidx+1)%(N+1) : thxl+1) : 1) #2 -- thxl+1
  return pidx
end

function uppA()
  global pidA
  N = nprocs()
  pidA = (N > 1 ? ((pidA < N && pidA != 0) ? (pidA+1)%(N+1) : 2) : 1) #2 -- thxl+1
  return pidA
end

function packFromIncomingDensities!(dens::Array{BallTreeDensity,1},
      wfac::Vector{AbstractString},
      vertid::Int,
      inmsgs::Array{NBPMessage,1} )
  #
  for m in inmsgs
    for idx in keys(m.p)
      if idx == vertid
        pdi = m.p[vertid]
        push!(dens, kde!(pdi) ) # kde!(pdi.pts, pdi.bws)
        push!(wfac, "msg")
      end
      # TODO -- we can inprove speed of search for inner loop
    end
  end
  nothing
end

# add all potentials associated with this clique and vertid to dens
function packFromLocalPotentials!(fgl::FactorGraph,
      dens::Vector{BallTreeDensity},
      wfac::Vector{AbstractString},
      cliq::Graphs.ExVertex,
      vertid::Int,
      N::Int,
      dbg::Bool=false )
  #
  for idfct in cliq.attributes["data"].potentials
    vert = getVert(fgl, idfct, api=localapi)
    data = getData(vert)
    # skip partials here, will be caught in packFromLocalPartials!
    if length( find(data.fncargvID .== vertid) ) >= 1 && !data.fnc.partial
      p = findRelatedFromPotential(fgl, vert, vertid, N, dbg )
      push!(dens, p)
      push!(wfac, vert.label)
    end
  end
  nothing
end

function packFromLocalPartials!(fgl::FactorGraph,
      partials::Dict{Int, Vector{BallTreeDensity}},
      cliq::Graphs.ExVertex,
      vertid::Int,
      N::Int,
      dbg::Bool=false)
  #

  for idfct in cliq.attributes["data"].potentials
    vert = getVert(fgl, idfct, api=localapi)
    data = getData(vert)
    if length( find(data.fncargvID .== vertid) ) >= 1 && data.fnc.partial
      p = findRelatedFromPotential(fgl, vert, vertid, N, dbg)
      pardims = data.fnc.usrfnc!.partial
      for dimnum in pardims
        if haskey(partials, dimnum)
          push!(partials[dimnum], marginal(p,[dimnum]))
        else
          partials[dimnum] = BallTreeDensity[marginal(p,[dimnum])]
        end
      end
    end
  end
  nothing
end

"""
    $(SIGNATURES)

Multiply different dimensions from partial constraints individually.
"""
function productpartials!(pGM::Array{Float64,2}, dummy::BallTreeDensity,
        partials::Dict{Int, Vector{BallTreeDensity}}  )

  #
  for (dimnum,pp) in partials
    dummy2 = marginal(dummy,[dimnum]) # kde!(rand(1,N),[1.0])
    if length(pp) > 1
      pGM[dimnum,:], = prodAppxMSGibbsS(dummy2, pp, Union{}, Union{}, 8)
    else
      pGM[dimnum,:] = getPoints(pp[1])
    end
  end
  nothing
end

"""
    $(SIGNATURES)

Multiply various full and partial dimension constraints.
"""
function prodmultiplefullpartials( dens::Vector{BallTreeDensity},
      partials::Dict{Int, Vector{BallTreeDensity}},
      Ndims::Int, N::Int  )
  #
  dummy = kde!(rand(Ndims,N),[1.0]); # TODO -- reuse memory rather than rand here
  pGM, = prodAppxMSGibbsS(dummy, dens, Union{}, Union{}, 8) #10
  for (dimnum,pp) in partials
    push!(pp, kde!(pGM[dimnum,:]))
  end
  productpartials!(pGM, dummy, partials)
  return pGM
end

"""
    $(SIGNATURES)

Multiply a single full and several partial dimension constraints.
"""
function prodmultipleonefullpartials( dens::Vector{BallTreeDensity},
      partials::Dict{Int, Vector{BallTreeDensity}},
      Ndims::Int, N::Int  )
  #
  dummy = kde!(rand(Ndims,N),[1.0]) # TODO -- reuse memory rather than rand here
  denspts = getPoints(dens[1])
  pGM = deepcopy(denspts)
  for (dimnum,pp) in partials
    push!(pp, kde!(pGM[dimnum,:]))
  end
  productpartials!(pGM, dummy, partials)
  return pGM
end

function productbelief(fg::FactorGraph,
      vertid::Int,
      dens::Vector{BallTreeDensity},
      partials::Dict{Int, Vector{BallTreeDensity}},
      N::Int;
      dbg::Bool=false )
  #

  pGM = Array{Float64}(0,0)
  lennonp, lenpart = length(dens), length(partials)
  if lennonp > 1
    Ndims = Ndim(dens[1])
    print("[$(lennonp)x$(lenpart)p,d$(Ndims),N$(N)],")
    pGM = prodmultiplefullpartials(dens, partials, Ndims, N)
  elseif lennonp == 1 && lenpart >= 1
    Ndims = Ndim(dens[1])
    print("[$(lennonp)x$(lenpart)p,d$(Ndims),N$(N)],")
    pGM = prodmultipleonefullpartials(dens, partials, Ndims, N)
  elseif lennonp == 0 && lenpart >= 1
    denspts = getVal(fg,vertid,api=localapi)
    Ndims = size(denspts,1)
    print("[$(lennonp)x$(lenpart)p,d$(Ndims),N$(N)],")
    dummy = kde!(rand(Ndims,N),[1.0]) # TODO -- reuse memory rather than rand here
    pGM = deepcopy(denspts)
    productpartials!(pGM, dummy, partials)
  # elseif lennonp == 0 && lenpart == 1
  #   print("[prtl]")
  #   pGM = deepcopy(getVal(fg,vertid,api=localapi) )
  #   for (dimnum,pp) in partials
  #     pGM[dimnum,:] = getPoints(pp)
  #   end
  elseif lennonp == 1 && lenpart == 0
    print("[drct]")
    pGM = getPoints(dens[1])
  else
    warn("Unknown density product on vertid=$(vertid), lennonp=$(lennonp), lenpart=$(lenpart)")
    pGM = Array{Float64,2}(0,1)
  end

  return pGM
end

function proposalbeliefs!(fgl::FactorGraph,
      destvertid::Int,
      factors::Vector{Graphs.ExVertex},
      dens::Vector{BallTreeDensity},
      partials::Dict{Int, Vector{BallTreeDensity}};
      N::Int=100,
      dbg::Bool=false )
  #
  for fct in factors
    data = getData(fct)
    p = findRelatedFromPotential(fgl, fct, destvertid, N, dbg)
    if data.fnc.partial   # partial density
      pardims = data.fnc.usrfnc!.partial
      for dimnum in pardims
        if haskey(partials, dimnum)
          push!(partials[dimnum], marginal(p,[dimnum]))
        else
          partials[dimnum] = BallTreeDensity[marginal(p,[dimnum])]
        end
      end
    else # full density
      push!(dens, p)
    end
  end
  nothing
end

function predictbelief(fgl::FactorGraph,
                       destvert::ExVertex,
                       factors::Vector{Graphs.ExVertex};
                       N::Int=0,
                       dbg::Bool=false )
  #
  destvertid = destvert.index
  dens = Array{BallTreeDensity,1}()
  partials = Dict{Int, Vector{BallTreeDensity}}()

  # determine number of particles to draw from the marginal
  nn = N != 0 ? N : size(getVal(destvert),2)

  # get proposal beliefs
  proposalbeliefs!(fgl, destvertid, factors, dens, partials, N=nn, dbg=dbg)

  # take the product
  pGM = productbelief(fgl, destvertid, dens, partials, nn, dbg=dbg )

  return pGM
end

function predictbelief(fgl::FactorGraph,
                       destvertsym::Symbol,
                       factorsyms::Vector{Symbol};
                       N::Int=0,
                       dbg::Bool=false,
                       api::DataLayerAPI=IncrementalInference.localapi  )
  #
  factors = Graphs.ExVertex[]
  for fsym in factorsyms
    push!(factors, getVert(fgl, fgl.fIDs[fsym], api=api))
  end
  vert = getVert(fgl, destvertsym, api=api)

  # determine the number of particles to draw from the marginal
  nn = N != 0 ? N : size(getVal(vert),2)

  # do the belief prediction
  predictbelief(fgl, vert, factors, N=nn, dbg=dbg)
end

function predictbelief(fgl::FactorGraph,
                       destvertsym::Symbol,
                       factorsyms::Colon;
                       N::Int=0,
                       dbg::Bool=false,
                       api::DataLayerAPI=IncrementalInference.localapi  )
  #
  predictbelief(fgl, destvertsym, ls(fgl, destvertsym, api=api), N=N, api=api, dbg=dbg )
end

function localProduct(fgl::FactorGraph,
                      sym::Symbol;
                      N::Int=100,
                      dbg::Bool=false,
                      api::DataLayerAPI=IncrementalInference.dlapi  )
  # TODO -- converge this function with predictbelief for this node
  # TODO -- update to use getVertId
  destvertid = fgl.IDs[sym] #destvert.index
  dens = Array{BallTreeDensity,1}()
  partials = Dict{Int, Vector{BallTreeDensity}}()
  lb = String[]
  fcts = Vector{Graphs.ExVertex}()
  cf = ls(fgl, sym, api=api)
  for f in cf
    vert = getVert(fgl,f, nt=:fnc, api=api)
    push!(fcts, vert)
    push!(lb, vert.label)
  end

  # get proposal beliefs
  proposalbeliefs!(fgl, destvertid, fcts, dens, partials, N=N, dbg=dbg)

  # take the product
  pGM = productbelief(fgl, destvertid, dens, partials, N, dbg=dbg )
  pp = kde!(pGM)

  return pp, dens, partials, lb
end
localProduct{T <: AbstractString}(fgl::FactorGraph, lbl::T; N::Int=100, dbg::Bool=false) = localProduct(fgl, Symbol(lbl), N=N, dbg=dbg)

"""
    initializeNode!(::FactorGraph, ::Symbol; N::Int=100, api::DataLayerAPI=dlapi)

Initialize the belief of a variable node in the factor graph struct.
"""
function initializeNode!(fgl::FactorGraph,
        sym::Symbol;
        N::Int=100,
        api::DataLayerAPI=IncrementalInference.dlapi )
  #

  vert = getVert(fgl, sym, api=api)
  # TODO -- this localapi is inconsistent, but get internal error due to problem with ls(fg, api=dlapi)
  belief,b,c,d  = localProduct(fgl, sym, api=localapi)
  pts = getPoints(belief)
  @show "initializing", sym, size(pts), Base.mean(pts,2), Base.std(pts,2)
  setVal!(vert, pts)
  api.updatevertex!(fgl, vert)

  nothing
end

function cliqGibbs(fg::FactorGraph,
      cliq::Graphs.ExVertex,
      vertid::Int,
      inmsgs::Array{NBPMessage,1},
      N::Int,
      dbg::Bool )
  #
  # several optimizations can be performed in this function TODO

  #consolidate NBPMessages and potentials
  dens = Array{BallTreeDensity,1}()
  partials = Dict{Int, Vector{BallTreeDensity}}()
  wfac = Vector{AbstractString}()
  packFromIncomingDensities!(dens, wfac, vertid, inmsgs)
  packFromLocalPotentials!(fg, dens, wfac, cliq, vertid, N)
  packFromLocalPartials!(fg, partials, cliq, vertid, N, dbg)

  potprod = !dbg ? nothing : PotProd(vertid, getVal(fg,vertid,api=localapi), Array{Float64}(0,0), dens, wfac)
  pGM = productbelief(fg, vertid, dens, partials, N, dbg=dbg )
  if dbg  potprod.product = pGM  end

  print(" ")
  return pGM, potprod
end

function fmcmc!(fgl::FactorGraph,
      cliq::Graphs.ExVertex,
      fmsgs::Vector{NBPMessage},
      IDs::Vector{Int},
      N::Int,
      MCMCIter::Int,
      dbg::Bool=false  )
  #
    println("---------- successive fnc approx ------------$(cliq.attributes["label"])")
    # repeat several iterations of functional Gibbs sampling for fixed point convergence
    if length(IDs) == 1
        MCMCIter=1
    end
    mcmcdbg = Array{CliqGibbsMC,1}()

    for iter in 1:MCMCIter
      # iterate through each of the variables, KL-divergence tolerence would be nice test here
      print("#$(iter)\t -- ")
      dbgvals = !dbg ? nothing : CliqGibbsMC([], Symbol[])
      for vertid in IDs
        # we'd like to do this more pre-emptive and then just execute -- just point and skip up only msgs
        densPts, potprod = cliqGibbs(fgl, cliq, vertid, fmsgs, N, dbg) #cliqGibbs(fg, cliq, vertid, fmsgs, N)
        if size(densPts,1)>0
          updvert = getVert(fgl, vertid, api=dlapi)
          setValKDE!(updvert, densPts)
          # Go update the datalayer TODO -- excessive for general case, could use local and update remote at end
          dlapi.updatevertex!(fgl, updvert)
          # fgl.v[vertid].attributes["val"] = densPts
          if dbg
            push!(dbgvals.prods, potprod)
            push!(dbgvals.lbls, Symbol(updvert.label))
          end
        end
      end
      !dbg ? nothing : push!(mcmcdbg, dbgvals)
      println("")
    end

    # populate dictionary for return NBPMessage in multiple dispatch
    # TODO -- change to EasyMessage dict
    d = Dict{Int,EasyMessage}() # Array{Float64,2}
    for vertid in IDs
      vert = dlapi.getvertex(fgl,vertid)
      pden = getKDE(vert)
      bws = vec(getBW(pden)[:,1])
      d[vertid] = EasyMessage(getVal(vert), bws)
      # d[vertid] = getVal(dlapi.getvertex(fgl,vertid)) # fgl.v[vertid]
    end
    println("fmcmc! -- finished on $(cliq.attributes["label"])")

    return mcmcdbg, d
end

function upPrepOutMsg!(d::Dict{Int,EasyMessage}, IDs::Array{Int,1}) #Array{Float64,2}
  print("Outgoing msg density on: ")
  len = length(IDs)
  m = NBPMessage(Dict{Int,EasyMessage}())
  for id in IDs
    m.p[id] = d[id]
  end
  return m
end

function upGibbsCliqueDensity(inp::ExploreTreeType, N::Int=200, dbg::Bool=false)
    print("up w $(length(inp.sendmsgs)) msgs")
    # Local mcmc over belief functions
    # this is so slow! TODO Can be ignored once we have partial working
    # loclfg = nprocs() < 2 ? deepcopy(inp.fg) : inp.fg

    # TODO -- some weirdness with: d,. = d = ., nothing
    mcmcdbg, d = Array{CliqGibbsMC,1}(), Dict{Int,EasyMessage}()

    priorprods = Vector{CliqGibbsMC}()

    if false
      IDS = [inp.cliq.attributes["data"].frontalIDs;inp.cliq.attributes["data"].conditIDs] #inp.cliq.attributes["frontalIDs"]
      mcmcdbg, d = fmcmc!(inp.fg, inp.cliq, inp.sendmsgs, IDS, N, 3, dbg)
    elseif true
      dummy, d = fmcmc!(inp.fg, inp.cliq, inp.sendmsgs, inp.cliq.attributes["data"].directFrtlMsgIDs, N, 1)
      if length(inp.cliq.attributes["data"].msgskipIDs) > 0
        dummy, dd = fmcmc!(inp.fg, inp.cliq, inp.sendmsgs, inp.cliq.attributes["data"].msgskipIDs, N, 1)
        for md in dd d[md[1]] = md[2]; end
      end
      # NOTE -- previous mistake, must iterate over directsvarIDs also
      if length(inp.cliq.attributes["data"].itervarIDs) > 0
        mcmcdbg, ddd = fmcmc!(inp.fg, inp.cliq, inp.sendmsgs, inp.cliq.attributes["data"].itervarIDs, N, 3, dbg)
        for md in ddd d[md[1]] = md[2]; end
      end
      if length(inp.cliq.attributes["data"].directPriorMsgIDs) > 0
        doids = setdiff(inp.cliq.attributes["data"].directPriorMsgIDs, inp.cliq.attributes["data"].msgskipIDs) # based on RYPKEMA2 example
        priorprods, dddd = fmcmc!(inp.fg, inp.cliq, inp.sendmsgs, doids, N, 1, dbg)
        for md in dddd d[md[1]] = md[2]; end
      end
    end

    #m = upPrepOutMsg!(inp.fg, inp.cliq, inp.sendmsgs, condids, N)
    m = upPrepOutMsg!(d, inp.cliq.attributes["data"].conditIDs)

    outmsglbl = Dict{Symbol, Int}()
    if dbg
      for (ke, va) in m.p
        outmsglbl[Symbol(inp.fg.g.vertices[ke].label)] = ke
      end
    end

    # Copy frontal variables back
    # for id in inp.cliq.attributes["frontalIDs"]
    #     inp.fg.v[id].attributes["val"] = loclfg.v[id].attributes["val"] # inp.
    # end
    # @show getVal(inp.fg.v[1])
    mdbg = !dbg ? DebugCliqMCMC() : DebugCliqMCMC(mcmcdbg, m, outmsglbl, priorprods)
    return UpReturnBPType(m, mdbg, d)
end
# else
#   rr = Array{Future,1}(4)
#   rr[1] = remotecall(upp2(), fmcmc!, inp.fg, inp.cliq, inp.sendmsgs, inp.cliq.attributes["directFrtlMsgIDs"], N, 1)
#   rr[2] = remotecall(upp2(), fmcmc!, inp.fg, inp.cliq, inp.sendmsgs, inp.cliq.attributes["msgskipIDs"], N, 1)
#   rr[3] = remotecall(upp2(), fmcmc!, inp.fg, inp.cliq, inp.sendmsgs, inp.cliq.attributes["itervarIDs"], N, 3)
#   rr[4] = remotecall(upp2(), fmcmc!, inp.fg, inp.cliq, inp.sendmsgs, inp.cliq.attributes["directvarIDs"], N, 1)
#   dummy, d = fetch(rr[1])
#   dummy, dd = fetch(rr[2])
#   mcmcdbg, ddd = fetch(rr[3])
#   dummy, dddd = fetch(rr[4])
#
#   for md in dd d[md[1]] = md[2]; end
#   for md in ddd d[md[1]] = md[2]; end
#   for md in dddd d[md[1]] = md[2]; end

function dwnPrepOutMsg(fg::FactorGraph, cliq::Graphs.ExVertex, dwnMsgs::Array{NBPMessage,1}, d::Dict{Int, EasyMessage}) #Array{Float64,2}
    # pack all downcoming conditionals in a dictionary too.
    if cliq.index != 1
      println("Dwn msg keys $(keys(dwnMsgs[1].p))")
    end # ignore root, now incoming dwn msg
    print("Outgoing msg density on: ")
    m = NBPMessage(Dict{Int,EasyMessage}())
    i = 0
    for vid in cliq.attributes["data"].frontalIDs
      m.p[vid] = deepcopy(d[vid]) # TODO -- not sure if deepcopy is required
      # outp = kde!(d[vid], "lcv") # need to find new down msg bandwidths
      # # i+=1
      # # outDens[i] = outp
      # bws = vec((getBW(outp))[:,1])
      # m.p[vid] = EasyMessage(  deepcopy(d[vid]) , bws  )
    end
    for cvid in cliq.attributes["data"].conditIDs
        i+=1
        # TODO -- convert to points only since kde replace by rkhs in future
        # outDens[i] = cdwndict[cvid]
        println("")
        # println("Looking for cvid=$(cvid)")
        m.p[cvid] = deepcopy(dwnMsgs[1].p[cvid]) # TODO -- maybe this can just be a union(,)
    end
    return m
end

function downGibbsCliqueDensity(fg::FactorGraph, cliq::Graphs.ExVertex, dwnMsgs::Array{NBPMessage,1}, N::Int=200, MCMCIter::Int=3, dbg::Bool=false)
    print("dwn")
    mcmcdbg, d = fmcmc!(fg, cliq, dwnMsgs, cliq.attributes["data"].frontalIDs, N, MCMCIter, dbg)
    m = dwnPrepOutMsg(fg, cliq, dwnMsgs, d)

    outmsglbl = Dict{Symbol, Int}()
    if dbg
      for (ke, va) in m.p
        outmsglbl[Symbol(fg.g.vertices[ke].label)] = ke
      end
    end

    mdbg = !dbg ? DebugCliqMCMC() : DebugCliqMCMC(mcmcdbg, m, outmsglbl, CliqGibbsMC[])
    return DownReturnBPType(m, mdbg, d)
end

function updateFGBT!(fg::FactorGraph, bt::BayesTree, cliqID::Int, ddt::DownReturnBPType; dbg::Bool=false)
    # if dlapi.cgEnabled
    #   return nothing
    # end
    cliq = bt.cliques[cliqID]
    if dbg
      cliq.attributes["debugDwn"] = deepcopy(ddt.dbgDwn)
    end
    for dat in ddt.IDvals
      #TODO -- should become an update call
        updvert = dlapi.getvertex(fg,dat[1])
        setValKDE!(updvert, deepcopy(dat[2])) # TODO -- not sure if deepcopy is required
        # updvert.attributes["latestEst"] = Base.mean(dat[2],2)
        dlapi.updatevertex!(fg, updvert, updateMAPest=true)
    end
    nothing
end

# TODO -- use Union{} for two types, rather than separate functions
function updateFGBT!(fg::FactorGraph, bt::BayesTree, cliqID::Int, urt::UpReturnBPType; dbg::Bool=false)
    # if dlapi.cgEnabled
    #   return nothing
    # end
    cliq = bt.cliques[cliqID]
    cliq = bt.cliques[cliqID]
    if dbg
      cliq.attributes["debug"] = deepcopy(urt.dbgUp)
    end
    for dat in urt.IDvals
      updvert = dlapi.getvertex(fg,dat[1])
      setValKDE!(updvert, deepcopy(dat[2])) # (fg.v[dat[1]], ## TODO -- not sure if deepcopy is required
      dlapi.updatevertex!(fg, updvert, updateMAPest=true)
    end
    println("updateFGBT! up -- finished updating $(cliq.attributes["label"])")
    nothing
end

# pass NBPMessages back down the tree -- pre order tree traversal
function downMsgPassingRecursive(inp::ExploreTreeType; N::Int=200, dbg::Bool=false)
    println("====================== Clique $(inp.cliq.attributes["label"]) =============================")

    mcmciter = inp.prnt != Union{} ? 3 : 0; # skip mcmc in root on dwn pass
    rDDT = downGibbsCliqueDensity(inp.fg, inp.cliq, inp.sendmsgs, N, mcmciter, dbg) #dwnMsg
    updateFGBT!(inp.fg, inp.bt, inp.cliq.index, rDDT, dbg=dbg)

    # rr = Array{Future,1}()
    pcs = procs()

    ddt=Union{}
    for child in out_neighbors(inp.cliq, inp.bt.bt)
        ett = ExploreTreeType(inp.fg, inp.bt, child, inp.cliq, [rDDT.dwnMsg])#inp.fg
        ddt = downMsgPassingRecursive( ett , N=N, dbg=dbg )
    end

    # return modifications to factorgraph to calling process
    return ddt
end

# post order tree traversal and build potential functions
function upMsgPassingRecursive(inp::ExploreTreeType; N::Int=200, dbg::Bool=false) #upmsgdict = Dict{Int, Array{Float64,2}}()
    println("Start Clique $(inp.cliq.attributes["label"]) =============================")
    childMsgs = Array{NBPMessage,1}()

    len = length(out_neighbors(inp.cliq, inp.bt.bt))
    for child in out_neighbors(inp.cliq, inp.bt.bt)
        ett = ExploreTreeType(inp.fg, inp.bt, child, inp.cliq, NBPMessage[]) # ,Union{})
        println("upMsgRec -- calling new recursive on $(ett.cliq.attributes["label"])")
        newmsgs = upMsgPassingRecursive(  ett, N=N, dbg=dbg ) # newmsgs
        println("upMsgRec -- finished with $(ett.cliq.attributes["label"]), w $(keys(newmsgs.p)))")
        push!(  childMsgs, newmsgs )
    end

    println("====================== Clique $(inp.cliq.attributes["label"]) =============================")
    ett = ExploreTreeType(inp.fg, inp.bt, inp.cliq, Union{}, childMsgs)

    urt = upGibbsCliqueDensity(ett, N, dbg) # upmsgdict
    updateFGBT!(inp.fg, inp.bt, inp.cliq.index, urt, dbg=dbg)
    println("End Clique $(inp.cliq.attributes["label"]) =============================")
    urt.upMsgs
end

function downGibbsCliqueDensity(inp::ExploreTreeType, N::Int=200, dbg::Bool=false)
  println("=================== Iter Clique $(inp.cliq.attributes["label"]) ===========================")
  mcmciter = inp.prnt != Union{} ? 3 : 0
  return downGibbsCliqueDensity(inp.fg, inp.cliq, inp.sendmsgs, N, mcmciter, dbg)
end

function prepDwnPreOrderStack!(bt::BayesTree, parentStack::Array{Graphs.ExVertex,1})
  # dwn message passing function
  nodedata = Union{}
  tempStack = Array{Graphs.ExVertex,1}()
  push!(tempStack, parentStack[1])

  while ( length(tempStack) != 0 || nodedata != Union{} )
    if nodedata != Union{}
      for child in out_neighbors(nodedata, bt.bt) #nodedata.cliq, nodedata.bt.bt
          ## TODO -- don't yet have the NBPMessage results to propagate belief
          ## should only construct ett during processing
          push!(parentStack, child) #ett
          push!(tempStack, child) #ett
      end
      nodedata = Union{} # its over but we still need to complete computation in each of the leaves of the tree
    else
      nodedata = tempStack[1]
      deleteat!(tempStack, 1)
    end
  end
  nothing
end

function findVertsAssocCliq(fgl::FactorGraph, cliq::Graphs.ExVertex)

  IDS = [cliq.attributes["data"].frontalIDs;cliq.attributes["data"].conditIDs] #inp.cliq.attributes["frontalIDs"]

  error("findVertsAssocCliq -- not completed yet")
  nothing
end

function partialExploreTreeType(pfg::FactorGraph, pbt::BayesTree, cliqCursor::Graphs.ExVertex, prnt, pmsgs::Array{NBPMessage,1})
    # println("starting pett")
    # TODO -- expand this to grab only partial subsection from the fg and bt data structures

    if length(pmsgs) < 1
      return ExploreTreeType(pfg, pbt, cliqCursor, prnt, NBPMessage[])
    else
      return ExploreTreeType(pfg, pbt, cliqCursor, prnt, pmsgs)
    end
    nothing
end

function dispatchNewDwnProc!(fg::FactorGraph,
      bt::BayesTree,
      parentStack::Array{Graphs.ExVertex,1},
      stkcnt::Int,
      refdict::Dict{Int,Future};
      N::Int=200,
      dbg::Bool=false  )
  #
  cliq = parentStack[stkcnt]
  while !haskey(refdict, cliq.index) # nodedata.cliq
    sleep(0.25)
  end

  rDDT = fetch(refdict[cliq.index]) #nodedata.cliq
  delete!(refdict,cliq.index) # nodedata

  if rDDT != Union{}
    updateFGBT!(fg, bt, cliq.index, rDDT, dbg=dbg)
  end

  emptr = BayesTree(Union{}, 0, Dict{Int,Graphs.ExVertex}(), Dict{String,Int}());

  for child in out_neighbors(cliq, bt.bt) # nodedata.cliq, nodedata.bt.bt
      haskey(refdict, child.index) ? error("dispatchNewDwnProc! -- why you already have dwnremoteref?") : nothing
      ett = partialExploreTreeType(fg, emptr, child, cliq, [rDDT.dwnMsg]) # bt
      refdict[child.index] = remotecall(downGibbsCliqueDensity, upp2() , ett, N) # Julia 0.5 swapped order
  end
  nothing
end

function processPreOrderStack!(fg::FactorGraph,
      bt::BayesTree,
      parentStack::Array{Graphs.ExVertex,1},
      refdict::Dict{Int,Future};
      N::Int=200,
      dbg::Bool=false  )
  #
    # dwn message passing function for iterative tree exploration
    stkcnt = 0

    @sync begin
      sendcnt = 1:length(parentStack) # separate memory for remote calls
      for i in 1:sendcnt[end]
          @async dispatchNewDwnProc!(fg, bt, parentStack, sendcnt[i], refdict, N=N, dbg=dbg) # stkcnt ##pidxI,nodedata
      end
    end
    nothing
end

function downMsgPassingIterative!(startett::ExploreTreeType; N::Int=200, dbg::Bool=false)
  # this is where we launch the downward iteration process from
  parentStack = Array{Graphs.ExVertex,1}()
  refdict = Dict{Int,Future}()

  # start at the given clique in the tree -- shouldn't have to be the root.
  pett = partialExploreTreeType(startett.fg, startett.bt, startett.cliq,
                                        startett.prnt, startett.sendmsgs)
  refdict[startett.cliq.index] = remotecall(downGibbsCliqueDensity, upp2(), pett, N)  # for Julia 0.5

  push!(parentStack, startett.cliq ) # r

  prepDwnPreOrderStack!(startett.bt, parentStack)
  processPreOrderStack!(startett.fg, startett.bt, parentStack, refdict, N=N, dbg=dbg)

  println("dwnward leftovers, $(keys(refdict))")
  nothing
end

function prepPostOrderUpPassStacks!(bt::BayesTree, parentStack::Array{Graphs.ExVertex,1}, childStack::Array{Graphs.ExVertex,1})
  # upward message passing preparation
  while ( length(parentStack) != 0 )
      #2.1 Pop a node from first stack and push it to second stack
      cliq = parentStack[end]
      deleteat!(parentStack, length(parentStack))
      push!(childStack, cliq )

      #2.2 Push left and right children of the popped node to first stack
      for child in out_neighbors(cliq, bt.bt) # nodedata.cliq, nodedata.bt.bt
          push!(parentStack, child )
      end
  end
  for child in childStack
    @show child.attributes["label"]
  end
  nothing
end

# for up message passing
function asyncProcessPostStacks!(fgl::FactorGraph,
      bt::BayesTree,
      chldstk::Vector{Graphs.ExVertex},
      stkcnt::Int,
      refdict::Dict{Int,Future};
      N::Int=200,
      dbg::Bool=false  )
  #
  if stkcnt == 0
    println("asyncProcessPostStacks! ERROR stkcnt=0")
    error("asyncProcessPostStacks! stkcnt=0")
  end
  cliq = chldstk[stkcnt]
  gomulti = true
  println("Start Clique $(cliq.attributes["label"]) =============================")
  childMsgs = Array{NBPMessage,1}()
  ur = nothing
  for child in out_neighbors(cliq, bt.bt)
      println("asyncProcessPostStacks -- $(stkcnt), cliq=$(cliq.attributes["label"]), start on child $(child.attributes["label"]) haskey=$(haskey(child.attributes, "remoteref"))")
        while !haskey(refdict, child.index)
          # println("Sleeping $(cliq.attributes["label"]) on lack of remoteref from $(child.attributes["label"])")
          # @show child.index, keys(refdict)
          sleep(0.25)
        end

      if gomulti
        ur = fetch(refdict[child.index])
      else
        ur = child.attributes["remoteref"]
      end
      updateFGBT!( fgl, bt, child.index, ur, dbg=dbg ) # deep copies happen in the update function
      #delete!(child.attributes, "remoteref")

      push!(childMsgs, ur.upMsgs)

  end
  println("====================== Clique $(cliq.attributes["label"]) =============================")
  emptr = BayesTree(Union{}, 0, Dict{Int,Graphs.ExVertex}(), Dict{String,Int}());
  pett = partialExploreTreeType(fgl, emptr, cliq, Union{}, childMsgs) # bt   # parent cliq pointer is not needed here, fix Graphs.jl first

  if haskey(cliq.attributes, "remoteref")
      println("asyncProcessPostStacks! -- WHY YOU ALREADY HAVE REMOTEREF?")
  end

  newprocid = upp2()
  if gomulti
    refdict[cliq.index] = remotecall(upGibbsCliqueDensity, newprocid, pett, N, dbg ) # swap order for Julia 0.5
  else
    # bad way to do non multi test
    cliq.attributes["remoteref"] = upGibbsCliqueDensity(pett, N, dbg)
  end

  # delete as late as possible, but could happen sooner
  for child in out_neighbors(cliq, bt.bt)
      # delete!(child.attributes, "remoteref")
      delete!(refdict, child.index)
  end

  println("End Clique $(cliq.attributes["label"]) =============================")
  nothing
end

# upward belief propagation message passing function
function processPostOrderStacks!(fg::FactorGraph,
      bt::BayesTree,
      childStack::Array{Graphs.ExVertex,1};
      N::Int=200,
      dbg::Bool=false  )
  #

  refdict = Dict{Int,Future}()

  stkcnt = length(childStack)
  @sync begin
    sendcnt = stkcnt:-1:1 # separate stable memory
    for i in 1:stkcnt
        @async asyncProcessPostStacks!(fg, bt, childStack, sendcnt[i], refdict, N=N, dbg=dbg) # deepcopy(stkcnt)
    end
  end
  println("processPostOrderStacks! -- THIS ONLY HAPPENS AFTER SYNC")
  # we still need to fetch the root node computational output
  if true
    ur = fetch(refdict[childStack[1].index])
  else
    ur = childStack[1].attributes["remoteref"]
  end
  # delete!(childStack[1].attributes, "remoteref") # childStack[1].cliq
  delete!(refdict, childStack[1].index)

  println("upward leftovers, $(keys(refdict))")

  updateFGBT!(fg, bt, childStack[1].index, ur, dbg=dbg ) # nodedata
  nothing
end

function upMsgPassingIterative!(startett::ExploreTreeType; N::Int=200, dbg::Bool=false)
  #http://www.geeksforgeeks.org/iterative-postorder-traversal/
  # this is where we launch the downward iteration process from
  parentStack = Array{Graphs.ExVertex,1}()
  childStack = Array{Graphs.ExVertex,1}()
  #Loop while first stack is not empty
  push!(parentStack, startett.cliq )
  # Starting at the root means we have a top down view of the tree
  prepPostOrderUpPassStacks!(startett.bt, parentStack, childStack)
  processPostOrderStacks!(startett.fg, startett.bt, childStack, N=N, dbg=dbg)
  nothing
end

function inferOverTree!(fgl::FactorGraph, bt::BayesTree; N::Int=200, dbg::Bool=false)
    println("Ensure all nodes are initialized")
    ensureAllInitialized!(fgl)
    println("Do multi-process inference over tree")
    cliq = bt.cliques[1]
    upMsgPassingIterative!(ExploreTreeType(fgl, bt, cliq, Union{}, NBPMessage[]),N=N, dbg=dbg);
    cliq = bt.cliques[1]
    downMsgPassingIterative!(ExploreTreeType(fgl, bt, cliq, Union{}, NBPMessage[]),N=N, dbg=dbg);
    nothing
end

function inferOverTreeR!(fgl::FactorGraph, bt::BayesTree; N::Int=200, dbg::Bool=false)
    println("Ensure all nodes are initialized")
    ensureAllInitialized!(fgl)
    println("Do recursive inference over tree")
    cliq = bt.cliques[1]
    upMsgPassingRecursive(ExploreTreeType(fgl, bt, cliq, Union{}, NBPMessage[]), N=N, dbg=dbg);
    cliq = bt.cliques[1]
    downMsgPassingRecursive(ExploreTreeType(fgl, bt, cliq, Union{}, NBPMessage[]), N=N, dbg=dbg);
    nothing
end


text:
# Copyright (c) 2020 California Institute of Technology (“Caltech”). U.S.
# Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   • Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#   • Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#   • Neither the name of Caltech nor its operating division, the Jet
#   Propulsion Laboratory, nor the names of its contributors may be used to
#   endorse or promote products derived from this software without specific
#   prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

module Cat2Anlm

export cat2amln

using SharedArrays
using Healpix
using ..HealpixHelpers
#using ..HealPy
using ..SeparableArrays
using ..Modes

function sortout(rθϕ, nside)
    # Sort by r. This helps reduce recomputations in the spline for gnl(r).
    p = sortperm(rθϕ[1,:])
    rθϕ = rθϕ[:,p]

    ## Sort by healpix pixel, r stays sorted within a pixel.
    ## Note: The idea here is that sorting by pixels avoids an awkward access
    ## pattern when creating the healpix map for a given ℓ. However, we don't do
    ## this, as it also thins out the galaxies and that leads to temporary
    ## variables being recomputed much more frequently in the spline for gnl(r).
    #θ = collect(rθϕ[2,:])
    #ϕ = collect(rθϕ[3,:])
    #pix = ang2pixRing(nside, θ, ϕ)
    #p = sortperm(pix)
    #rθϕ = rθϕ[:,p]

    r = collect(rθϕ[1,:])
    θ = collect(rθϕ[2,:])
    ϕ = collect(rθϕ[3,:])
    return r, θ, ϕ
end

########################### transform jl ##########################################

function transform_gnl(npix, pix, r, sphbesg_nl)
    Ngal = length(pix)
    map = fill(0.0, npix)
    for i=1:Ngal
        #@show pix[i],r[i]
        map[pix[i]] += sphbesg_nl(r[i])
    end
    return map
end

# pmu = pixel management unit
function make_pmu_pmupix(pix)
    pmu = unique(sort(pix))
    pmupix = Array{Int64,1}(indexin(pix, pmu))
    return pmu, pmupix
end

function transform_gnl_spmap(npix, pmu, pmupix, r, sphbesg_nl)
    Ngal = length(pmupix)
    spmap = fill(0.0, length(pmu))
    for i=1:Ngal
        spmap[pmupix[i]] += sphbesg_nl(r[i])
    end
    map = fill(0.0, npix)
    @. map[pmu] = spmap
    return map
end

function transform_jnl_binned_quadratic(npix, pix, k, l, r, rmax)
    # Note: r must be sorted
    Ngal = length(pix)
    map = fill(0.0, npix)
    # Note: want Δr to be much less than distance between zeros, or k*Δr ≪ π
    nr = 25 * ceil(Int, 2 * k * rmax)
    Δr = rmax / nr
    #@show Δr,1/k,nr
    idx = 1
    ksphbes_prev = sphericalbesselj(l, 0.0)
    for i=1:nr
        rmid = rmax * (i-0.5) / nr
        rhi = rmax * i / nr
        ksphbes_mid = sphericalbesselj(l, k * rmid)
        ksphbes_next = sphericalbesselj(l, k * rhi)
        a = ((ksphbes_prev - ksphbes_mid) + (ksphbes_next - ksphbes_mid)) / Δr^2
        b = (ksphbes_next - ksphbes_prev) / Δr
        c = ksphbes_mid
        while idx <= Ngal && r[idx] <= rhi
            x = r[idx] - rmid
            map[pix[idx]+1] += (a * x + b) * x + c
            idx += 1
        end
        ksphbes_prev = ksphbes_next
    end
    return map
end

function gen_jl_cache(lmax, qmax)
    # Note: We need Δq=kₘₐₓ⋅Δr to be much less than the distance between zeros,
    # or Δq ≪ π. We choose Δq = π / (2π * 4), or about 25 per zero. Then, the
    # number of bins we need is nᵣ = qmax / Δq + 1 = 8 * qmax + 1. Nah, let's
    # do much more until we get convergence for (n,ℓ)=(1,0) mode.
    nr = ceil(Int, 50 * qmax) + 1
    q = range(0.0, qmax, length=nr)
    jl = [sphericalbesselj.(l, q) for l=0:lmax]
    return q, jl
end

function transform_jnl_binned_quadratic_cached(npix, pix, k, l, r, rmax, jl_q, jl)
    # Note: r must be sorted
    Ngal = length(pix)
    map = fill(0.0, npix)
    mynr = findfirst(jl_q .>= k*rmax)
    idx = 1
    rhi = 0.0
    #test_kjl = Float64[]
    #@show nr,mynr,qmax,k,Δr
    for i=2:mynr-1
        rlo = jl_q[i-1] / k
        rmid = jl_q[i] / k
        rhi = jl_q[i+1] / k
        #Ngal_Δr = length(r[@. rlo < r <= rhi])
        #@show i,l,rlo,rmid,rhi,idx,Ngal_Δr
        ksphbes0 = jl[i-1]
        ksphbes1 = jl[i]
        ksphbes2 = jl[i+1]
        t0 = rlo - rmid
        t2 = rhi - rmid
        ksphbes_deriv0 = (ksphbes0 - ksphbes1) / t0
        ksphbes_deriv2 = (ksphbes2 - ksphbes1) / t2
        a = (ksphbes_deriv2 - ksphbes_deriv0) / (t2 - t0)
        b = (t2 * ksphbes_deriv0 - t0 * ksphbes_deriv2) / (t2 - t0)
        c = ksphbes1
        while idx <= Ngal && r[idx] <= rhi && r[idx] <= rmax
            t = r[idx] - rmid
            map[pix[idx]+1] += (a * t + b) * t + c
            #push!(test_kjl, (a * t + b) * t + c)
            idx += 1
        end
    end
    #@assert rhi >= rmax  # numerical round-off
    @assert rhi + eps(rhi) >= rmax
    #figure()
    #plot(r, k * sphericalbesselj.(l, k * r), label="Exact")
    #scatter(jl_q./k, k .* jl, label="Knots")
    #plot(r, test_kjl, label="Approx")
    #plot(r, test_kjl ./ (k * sphericalbesselj.(l, k*r)) .- 1, label="Error")
    #xlabel(L"r")
    #ylabel(L"k\,j_\ell(kr)")
    #legend()
    return map
end

####################### catalogue -> a_nlm #####################################
@doc raw"""
    cat2amln(rθϕ, amodes, nbar, win_rhat_ln)

Computes the spherical Fourier-Bessel decomposition coefficients from a
catalogue of sources. The number density is measured from the survey as $\bar n
= N_\mathrm{gals} / V_\mathrm{eff}$.

# Example
```julia-repl
julia> using SphericalFourierBesselDecompositions
julia> cat2amln(rθϕ, ...)
```
"""
function cat2amln(rθϕ, amodes, nbar, win_rhat_ln)
    r, θ, ϕ = sortout(rθϕ, amodes.nside)
    @show nbar length(r)
    sphbesg = amodes.basisfunctions
    knl = amodes.knl
    pix = ang2pixRing.(Ref(Resolution(amodes.nside)), θ, ϕ) # .+ 1  # python is 0-indexed
    pmu, pmupix = make_pmu_pmupix(pix)
    @show typeof(pmu) typeof(pmupix)
    npix = nside2npix(amodes.nside)
    ΔΩpix = 4π / npix
    anlm = fill(NaN+im*NaN, getnlmsize(amodes))
    #anlm = SharedArray{Complex{Float64}}(getnlmsize(amodes))
    #anlm .= NaN+NaN*im
    for n=1:amodes.nmax
        @show n,amodes.nmax,0:amodes.lmax_n[n]
        #@time @sync @distributed for l=0:amodes.lmax_n[n]
        @time for l=0:amodes.lmax_n[n]
            #@show n,l
            #map0 = transform_gnl(npix, pix, r, sphbesg.gnl[n,l+1])
            map1 = transform_gnl_spmap(npix, pmu, pmupix, r, sphbesg.gnl[n,l+1])
            #@time map2 = transform_jnl_binned_quadratic(npix, pix, knl[n,l+1], l, r, rmax)
            #map3 = transform_jnl_binned_quadratic_cached(npix, pix, knl[n,l+1], l, r, rmax, jl_q, jl[l+1])
            map = map1
            #@show map
            #@show extrema(map1 ./ map0 .- 1)
            #@show extrema(map2 ./ map0 .- 1)
            #@show extrema(map3 ./ map0 .- 1)
            #hp.mollview(map0, title="direct")
            #hp.mollview(map1, title="direct, splined")
            #hp.mollview(map2, title="quadratic")
            #hp.mollview(map3, title="quadratic cached")
            #hp.mollview(map0 .- map, title="map1 - map")
            #hp.mollview(map1 .- map, title="map1 - map")
            #hp.mollview(map2 .- map, title="map2 - map")
            #hp.mollview(map3 .- map, title="map3 - map")
            #readline()
            #close("all")
            #@show mean(map ./ map1),std(map ./ map1)
            #@show mean(map ./ map2),std(map ./ map2)
            #@show mean(map ./ map3),std(map ./ map3)

            map .*= 1 / (nbar * ΔΩpix)

            c = win_rhat_ln[:,l+1,n]
            #@show size(map) size(c)
            #@show map c mean(map) median(map)
            @. map = map - c
            #@show n,l,mean(map),median(map)
            #@show map

            # TODO: check manual implementation for a single ℓ, check libsharp,
            #       how does healpix do it?
            # TODO: start with small nside for small ℓ, dangerous for pixel window

            # HealPy:
            #alm = hp.map2alm(map, lmax=l, use_weights=true)
            #idx = hp.Alm.getidx.(l, l, 0:l) .+ 1  # python is 0-indexed

            # Healpix.jl:
            maphp = HealpixMap{Float64,Healpix.RingOrder}(map)
            alm = mymap2alm(maphp, lmax=l)
            idx = almIndex(alm, l, 0:l)

            baseidx = getidx(amodes, n, l, 0)
            @. anlm[baseidx:(baseidx+l)] = alm.alm[idx]
        end
    end
    @assert all(isfinite.(anlm))
    return anlm
end

end

# vim: set sw=4 et sts=4 :


text:
module Systems

export PHPSystem,Tube,Evaporator,Condenser,Liquid,Vapor,Wall,Mapping,PHPResult

# using ..Tools

"""
PHPSystem is a struct containing
    γ
    Hc
    He
    θc
    θe
    ω0
    ζ
    L dimensionless pipe total length
    Xc dimensionless condenser range
    Xe dimensionless evaporater range
"""

mutable struct Tube
    d::Float64
    L::Float64
    L2D::Float64
    angle::Float64
    closedornot::Bool
end

"""
PHPSystem is a struct containing
    γ
    Hc
    He
    θc
    θe
    ω0
    ζ
    L dimensionless pipe total length
    Xc dimensionless condenser range
    Xe dimensionless evaporater range
"""

mutable struct Evaporator
    He::Float64
    θe::Float64
    Xe::Array{Tuple{Float64,Float64},1}
    We::Array{Float64,1}
end

"""
PHPSystem is a struct containing
    γ
    Hc
    He
    θc
    θe
    ω0
    ζ
    L dimensionless pipe total length
    Xc dimensionless condenser range
    Xe dimensionless evaporater range
"""

mutable struct Condenser
    Hc::Float64
    θc::Float64
    Xc::Array{Tuple{Float64,Float64},1}
    Hwc::Float64
end

"""
PHPSystem is a struct containing
    γ
    Hc
    He
    θc
    θe
    ω0
    ζ
    L dimensionless pipe total length
    Xc dimensionless condenser range
    Xe dimensionless evaporater range
"""

mutable struct Liquid
    γ::Float64
    Hₗ::Float64
    ρ::Float64
    ω0::Array{Float64,1}
    ℘::Array{Float64,1}
    Xp::Array{Tuple{Float64,Float64},1}
    dXdt::Array{Tuple{Float64,Float64},1}
    Xarrays::Array{Array{Float64,1},1}
    θarrays::Array{Array{Float64,1},1}

end

"""
PHPSystem is a struct containing
    γ
    Hc
    He
    θc
    θe
    ω0
    ζ
    L dimensionless pipe total length
    Xc dimensionless condenser range
    Xe dimensionless evaporater range
"""

mutable struct Vapor
    γ::Float64
    P::Array{Float64,1}
    δ::Array{Float64,1}
end

"""
PHPSystem is a struct containing
    γ
    Hc
    He
    θc
    θe
    ω0
    ζ
    L dimensionless pipe total length
    Xc dimensionless condenser range
    Xe dimensionless evaporater range
"""

mutable struct Wall
    α::Float64
    Δt::Float64
    Xstations::Array{Float64,1}
    Xarray::Array{Float64,1}
    θarray::Array{Float64,1}
end

"""
PHPSystem is a struct containing
    γ
    Hc
    He
    θc
    θe
    ω0
    ζ
    L dimensionless pipe total length
    Xc dimensionless condenser range
    Xe dimensionless evaporater range
"""

mutable struct Mapping
    walltoliquid::Array{Tuple{Int64,Int64},1}
    liquidtowall::Array{Array{Int64,1},1}
end

"""
PHPSystem is a struct containing
    γ
    Hc
    He
    θc
    θe
    ω0
    ζ
    L dimensionless pipe total length
    Xc dimensionless condenser range
    Xe dimensionless evaporater range
"""
# add type names!
mutable struct PHPSystem
    tube::Tube
    evaporator::Evaporator
    condenser::Condenser
    liquid::Liquid
    vapor::Vapor
    wall::Wall
    mapping::Mapping
end

"""
PHPSystem is a struct containing
    γ
    Hc
    He
    θc
    θe
    ω0
    ζ
    L dimensionless pipe total length
    Xc dimensionless condenser range
    Xe dimensionless evaporater range
"""

struct PHPResult
    t
    Xp
    dXdt
    P
    θ
    M
end

end


text:
include("test_interface_contracts.jl")
include("test_skip_lists.jl")


text:
using LinearAlgebra
using Parameters, StatsBase

function cv_linear(X::AbstractArray,y::Vector,K::Int)
    n = length(y); m = round(Int,n/K)
    S = 0
    for j in 1:K
        test = j*m-m+1:j*m
        train = setdiff(1:n,test)
        β̂ = multiple_regression(X[train,:],y[train])
        ŷ = insert_ones(X[test,:])*β̂
        S += dot(y[test]-ŷ,y[test]-ŷ)
    end
    return S/n
end

function cv_fast(X::AbstractArray,y::Vector,K::Int)
    n = length(y); m = round(Int,n/K)
    X = insert_ones(X)
    H = X*((X'X)\X')
    e = (I(n) - H)*y
    S = 0
    @views for j in 1:K
        test = j*m-m+1:j*m;
        err = (I(m)-H[test,test]) \ e[test]
        S += dot(err,err)
    end
    return S/n
end

@with_kw mutable struct bt{T<:Number}
    original::T = 0.0
    bias::T = 0.0
    stderr::T = 0.0
end

function bootstrap(df::AbstractArray{T}, F::Function, r::Int) where T <:Number
    m,_ = size(df)
    original = F(df,1:m)
    u = Vector{T}(undef, r)
    for i in 1:r
        index = sample(1:m,m,replace=true)
        u[i] = F(df,index)
    end
    bias = mean(u) - original; stderr = std(u)
    result = bt()
    @pack! result = original, bias, stderr
    return result
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -11.08168115564476 <= q <= 12.629895565298366)

                     add_NL_constraint(m, :(gelu(-0.8013086439641905 + -0.577380120181807*gelu(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4])) + -0.9073468716494122*gelu(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4])) + -0.34890000179172853*gelu(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))) + gelu(-0.6400914085711236 + 0.5999574270201746*gelu(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4])) + -0.3302091884595031*gelu(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4])) + 0.9262651374928534*gelu(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))) + gelu(0.7171327203675486 + 0.9873737246227612*gelu(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4])) + -0.03055409136192777*gelu(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4])) + 0.20254351925716962*gelu(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using DifferentiableNAS
using Flux
using Flux: logitcrossentropy, onecold, onehotbatch, Optimiser
using Zygote: @nograd
using StatsBase: mean
using Parameters
using CUDA
using Distributions
using BSON
using Dates
using Plots
include("../CIFAR10.jl")
include("../training_utils.jl")

argparams = trial_params()

m = DARTSModel(track_acts = true, num_cells = argparams.num_cells, channels = argparams.channels) |> gpu

optimiser_α = Optimiser(WeightDecay(1e-3),ADAM(3e-4,(0.5,0.999)))
optimiser_w = Optimiser(WeightDecay(3e-4),Momentum(0.025, 0.9))

train, val = get_processed_data(argparams.val_split, argparams.batchsize, argparams.trainval_fraction)
test = get_test_data(argparams.test_fraction)

function (hist::historiessml)()
    push!(hist.normal_αs_sm, softmax.(copy(m.normal_αs)) |> cpu)
    push!(hist.reduce_αs_sm, softmax.(copy(m.reduce_αs)) |> cpu)
    push!(hist.activations, copy(m.activations.currentacts) |> cpu)
    CUDA.reclaim()
    GC.gc()
end
histepoch = historiessml()
histbatch = historiessml()
losses = [0.0, 0.0]

base_folder = prepare_folder("acts")

cbepoch = CbAll(CUDA.reclaim, histepoch, save_progress, CUDA.reclaim)
cbbatch = CbAll(CUDA.reclaim, histbatch, CUDA.reclaim)

acts = activationpre(loss, m, val)
for epoch in 1:argparams.epochs
    @show epoch
    Activationtrain1st!(loss, m, train, val, optimiser_α, optimiser_w, acts, losses)
end

max_stars_repo_name: huangyxi/DifferentiableNAS

max_issues_repo_name: huangyxi/DifferentiableNAS

max_forks_repo_name: huangyxi/DifferentiableNAS


text:
using SatelliteSpots
using Documenter

DocMeta.setdocmeta!(SatelliteSpots, :DocTestSetup, :(using SatelliteSpots); recursive=true)

makedocs(;
    modules=[SatelliteSpots],
    authors="mileslucas <mdlucas@hawaii.edu> and contributors",
    repo="https://github.com/JuliaHCI/SatelliteSpots.jl/blob/{commit}{path}#{line}",
    sitename="SatelliteSpots.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaHCI.github.io/SatelliteSpots.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaHCI/SatelliteSpots.jl",
    devbranch="main",
)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c7de68ac-906c-4959-82ac-5fc8fe4870a5
using PlotlyJS, WebIO, Luxor, PlutoUI

# ╔═╡ 40415ac0-db9f-11eb-3ca4-1d76bfbbc79a
md"""
# Learn basic Pluto stuff
"""

# ╔═╡ 2e20d719-0f6e-46d3-89be-c73a56a6e44e
md"""
Plotly scatter plot
"""

# ╔═╡ 8e614340-d98d-4946-8da9-657594b5009d
@bind T3 Slider(6:16, show_value=true)

# ╔═╡ 2ab468f0-adce-4d2b-9f6f-bffe253f000b
begin
	l = Layout(
        title="Penguins",
        xaxis=attr(title="fish", showgrid=true),
        yaxis=attr(title="Weight", showgrid=true),
        legend=attr(x=0.1, y=0.1),
		width = 10,
		height = 10,
    )
    trace1 = scatter(;x=1:4, y=[10, 15, 13, 17], mode="markers")
    trace2 = scatter(;x=2:5, y=[16, 5, 11, 9], mode="lines")
    trace3 = scatter(;x=1:4, y=[12, 9, T3, 12], mode="lines+markers")
	# uncomment to display in Electron
	#plt = plot([trace1, trace2, trace3], l)
	#display(plt)

	# Display in-place
    p = Plot([trace1, trace2, trace3], l)   # for display in-place
end

# ╔═╡ 15c7a365-41cb-4ef8-b0d4-1eeeff3e15c7
@bind power_level html"<input type='range'>"

# ╔═╡ c4ce098b-dd0c-43a6-9380-d12da68e8d65
println("Power Level: $power_level")

# ╔═╡ 46085fcf-be98-44f9-872e-823e0513ff40
@bind S Slider(10:100, show_value=true)

# ╔═╡ 49448abd-bd53-4a60-9f7d-bf625dcbc177
@drawsvg begin
	background("slateblue4")
	for (k,i) in pairs(between.(O + (0,S), O - (0,S), range(0,1,length=3)))
		sethue([Luxor.julia_red, Luxor.julia_purple, Luxor.julia_green][k])
		Luxor.circle(i, S, :fill)
	end
end 100 100

# ╔═╡ Cell order:
# ╠═c7de68ac-906c-4959-82ac-5fc8fe4870a5
# ╠═40415ac0-db9f-11eb-3ca4-1d76bfbbc79a
# ╠═2e20d719-0f6e-46d3-89be-c73a56a6e44e
# ╠═2ab468f0-adce-4d2b-9f6f-bffe253f000b
# ╠═8e614340-d98d-4946-8da9-657594b5009d
# ╠═15c7a365-41cb-4ef8-b0d4-1eeeff3e15c7
# ╠═c4ce098b-dd0c-43a6-9380-d12da68e8d65
# ╠═46085fcf-be98-44f9-872e-823e0513ff40
# ╠═49448abd-bd53-4a60-9f7d-bf625dcbc177


text:
module SpringCollab2020StrawberryCollectionField

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/StrawberryCollectionField" StrawberryCollectionField(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    delayBetweenBerries::Bool=true, includeGoldens::Bool=false)

const placements = Ahorn.PlacementDict(
    "Strawberry Collection Field (Spring Collab 2020)" => Ahorn.EntityPlacement(
        StrawberryCollectionField,
        "rectangle"
    )
)

end


text:
mutable struct ThreadsCache{A, B, C}
    comm::A
    world_sz::B
    rank::C
end

function ThreadsCache(comm::TPI.Comm=TPI.Comm())
    return ThreadsCache(comm,
    TPI.Comm_size(comm),
    TPI.Comm_rank(comm))
end

function ThreadsCaches(n::Int)
    cms = TPI.Comms(n)
    return [ThreadsCache(c) for c=cms]
end

parallel_execution_cache(::ParallelThreaded) = ThreadsCache()
parallel_execution_cache(tc::ThreadsCache) = ThreadsCache(tc.comm)

num_workers(tc::ThreadsCache) = tc.world_sz

function workers_mean!(target, data, par_type::ThreadsCache)
    mean!(target, data)
    TPI.Allreduce!(target, target, +, par_type.comm)
    target ./= par_type.world_sz

    return target
end

function workers_mean!(target, par_type::ThreadsCache)
    TPI.Allreduce!(target, +, par_type.comm)
    target ./= par_type.world_sz

    return target
end

function workers_mean(data, par_type::ThreadsCache)
    data_m = mean(data)
    target = TPI.Allreduce(data_m, +, par_type.comm)

    return target/par_type.world_sz
end

function workers_sum!(target, par_type::ThreadsCache)
    TPI.Allreduce!(target, +, par_type.comm)
    return target
end

function all_bcast_from_root!(data, par_type::ThreadsCache)
    TPI.Bcast!(data, 1, par_type.comm)
    return data
end

function worker_local_seed(seed, par_cache::ThreadsCache)
    rng = Random.MersenneTwister(seed)
    local_seeds = rand(rng, eltype(seed), par_cache.world_sz)
    TPI.Bcast!(local_seeds, 1, par_cache.comm)
    return local_seeds[Threads.threadid()]
end

function worker_block(iter, par_cache::ThreadsCache)
    return _iterator_blocks(iter, Threads.threadid(), par_cache.world_sz)
end

my_block(stuff, par_cache::ThreadsCache) = stuff[Threads.threadid()]

function worker_allgatherv!(data, counts, par_cache::ThreadsCache)
    TPI.Allgatherv!(data, counts, par_cache.comm)
    return data
end


text:
guo_instances = DataDep("GCSPET_Guo",
            "This package is tyring to automatically install the original instances of the GCSPET. The files will be installed to the folder: ",
            "https://www.researchgate.net/profile/Peng_Guo9/publication/322538542_2017GCSPET_Instances/data/5a5ef5b0aca272d4a3e02dd4/2017GCSPET-Instances.rar")

instancedir(key) = @datadep_str key

function instancelist(key)
    dir = instancedir(key)
    if key == "GCSPET_Guo"
            small = readdir(joinpath(dir, "SGCSPET_Instances"))
            large = readdir(joinpath(dir, "LGCSPET_Instances"))
            return vcat(small, large)
    else
            return readdir(dir)
    end
end

function instancepath(key, name)
    dir = instancedir(key)
    if key == "GCSPET_Guo"
        if parse(Int, first(split(name, "_"))) < 30
            ifile = joinpath(dir, "SGCSPET_Instances", name)
        else
            ifile = joinpath(dir, "LGCSPET_Instances", name)
        end
    else
        ifile = joinpath(dir, name)
    end

    return ifile
end

const small_guo = "SGCSPET"
const large_guo = "LGCSPET"


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

struct TagsApi <: SwaggerApi
    client::Swagger.Client
end

"""
Creates a predefined tag name.
This operation allows adding a name to the list of predefined tag names for the given subscription. A tag name can have a maximum of 512 characters and is case-insensitive. Tag names cannot have the following prefixes which are reserved for Azure use: 'microsoft', 'azure', 'windows'.
Param: tagName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: TagDetails
"""
function _swaggerinternal_tagsCreateOrUpdate(_api::TagsApi, tagName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", TagDetails, "/subscriptions/{subscriptionId}/tagNames/{tagName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "tagName", tagName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function tagsCreateOrUpdate(_api::TagsApi, tagName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsCreateOrUpdate(_api, tagName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function tagsCreateOrUpdate(_api::TagsApi, response_stream::Channel, tagName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsCreateOrUpdate(_api, tagName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Creates or updates the entire set of tags on a resource or subscription.
This operation allows adding or replacing the entire set of tags on the specified resource or subscription. The specified entity can have a maximum of 50 tags.
Param: scope::String (required)
Param: api_version::String (required)
Param: parameters::TagsResource (required)
Return: TagsResource
"""
function _swaggerinternal_tagsCreateOrUpdateAtScope(_api::TagsApi, scope::String, api_version::String, parameters; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", TagsResource, "/{scope}/providers/Microsoft.Resources/tags/default", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "scope", scope)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function tagsCreateOrUpdateAtScope(_api::TagsApi, scope::String, api_version::String, parameters; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsCreateOrUpdateAtScope(_api, scope, api_version, parameters; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function tagsCreateOrUpdateAtScope(_api::TagsApi, response_stream::Channel, scope::String, api_version::String, parameters; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsCreateOrUpdateAtScope(_api, scope, api_version, parameters; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Creates a predefined value for a predefined tag name.
This operation allows adding a value to the list of predefined values for an existing predefined tag name. A tag value can have a maximum of 256 characters.
Param: tagName::String (required)
Param: tagValue::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: TagValue
"""
function _swaggerinternal_tagsCreateOrUpdateValue(_api::TagsApi, tagName::String, tagValue::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", TagValue, "/subscriptions/{subscriptionId}/tagNames/{tagName}/tagValues/{tagValue}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "tagName", tagName)  # type String
    Swagger.set_param(_ctx.path, "tagValue", tagValue)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function tagsCreateOrUpdateValue(_api::TagsApi, tagName::String, tagValue::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsCreateOrUpdateValue(_api, tagName, tagValue, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function tagsCreateOrUpdateValue(_api::TagsApi, response_stream::Channel, tagName::String, tagValue::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsCreateOrUpdateValue(_api, tagName, tagValue, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Deletes a predefined tag name.
This operation allows deleting a name from the list of predefined tag names for the given subscription. The name being deleted must not be in use as a tag name for any resource. All predefined values for the given name must have already been deleted.
Param: tagName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_tagsDelete(_api::TagsApi, tagName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/tagNames/{tagName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "tagName", tagName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function tagsDelete(_api::TagsApi, tagName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsDelete(_api, tagName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function tagsDelete(_api::TagsApi, response_stream::Channel, tagName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsDelete(_api, tagName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Deletes the entire set of tags on a resource or subscription.

Param: scope::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_tagsDeleteAtScope(_api::TagsApi, scope::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/{scope}/providers/Microsoft.Resources/tags/default", ["azure_auth"])
    Swagger.set_param(_ctx.path, "scope", scope)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function tagsDeleteAtScope(_api::TagsApi, scope::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsDeleteAtScope(_api, scope, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function tagsDeleteAtScope(_api::TagsApi, response_stream::Channel, scope::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsDeleteAtScope(_api, scope, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Deletes a predefined tag value for a predefined tag name.
This operation allows deleting a value from the list of predefined values for an existing predefined tag name. The value being deleted must not be in use as a tag value for the given tag name for any resource.
Param: tagName::String (required)
Param: tagValue::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_tagsDeleteValue(_api::TagsApi, tagName::String, tagValue::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/tagNames/{tagName}/tagValues/{tagValue}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "tagName", tagName)  # type String
    Swagger.set_param(_ctx.path, "tagValue", tagValue)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function tagsDeleteValue(_api::TagsApi, tagName::String, tagValue::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsDeleteValue(_api, tagName, tagValue, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function tagsDeleteValue(_api::TagsApi, response_stream::Channel, tagName::String, tagValue::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsDeleteValue(_api, tagName, tagValue, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Gets the entire set of tags on a resource or subscription.

Param: scope::String (required)
Param: api_version::String (required)
Return: TagsResource
"""
function _swaggerinternal_tagsGetAtScope(_api::TagsApi, scope::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", TagsResource, "/{scope}/providers/Microsoft.Resources/tags/default", ["azure_auth"])
    Swagger.set_param(_ctx.path, "scope", scope)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function tagsGetAtScope(_api::TagsApi, scope::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsGetAtScope(_api, scope, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function tagsGetAtScope(_api::TagsApi, response_stream::Channel, scope::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsGetAtScope(_api, scope, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Gets a summary of tag usage under the subscription.
This operation performs a union of predefined tags, resource tags, resource group tags and subscription tags, and returns a summary of usage for each tag name and value under the given subscription. In case of a large number of tags, this operation may return a previously cached result.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: TagsListResult
"""
function _swaggerinternal_tagsList(_api::TagsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", TagsListResult, "/subscriptions/{subscriptionId}/tagNames", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function tagsList(_api::TagsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsList(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function tagsList(_api::TagsApi, response_stream::Channel, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsList(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Selectively updates the set of tags on a resource or subscription.
This operation allows replacing, merging or selectively deleting tags on the specified resource or subscription. The specified entity can have a maximum of 50 tags at the end of the operation. The 'replace' option replaces the entire set of existing tags with a new set. The 'merge' option allows adding tags with new names and updating the values of tags with existing names. The 'delete' option allows selectively deleting tags based on given names or name/value pairs.
Param: scope::String (required)
Param: api_version::String (required)
Param: parameters::TagsPatchResource (required)
Return: TagsResource
"""
function _swaggerinternal_tagsUpdateAtScope(_api::TagsApi, scope::String, api_version::String, parameters; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", TagsResource, "/{scope}/providers/Microsoft.Resources/tags/default", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "scope", scope)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function tagsUpdateAtScope(_api::TagsApi, scope::String, api_version::String, parameters; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsUpdateAtScope(_api, scope, api_version, parameters; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function tagsUpdateAtScope(_api::TagsApi, response_stream::Channel, scope::String, api_version::String, parameters; _mediaType=nothing)
    _ctx = _swaggerinternal_tagsUpdateAtScope(_api, scope, api_version, parameters; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export tagsCreateOrUpdate, tagsCreateOrUpdateAtScope, tagsCreateOrUpdateValue, tagsDelete, tagsDeleteAtScope, tagsDeleteValue, tagsGetAtScope, tagsList, tagsUpdateAtScope


text:
module MPSGE

import JuMP, Complementarity, MacroTools
import JuMP: value, set_value, Containers
import JuMP.Containers: DenseAxisArray

export add!, Model, solve!, algebraic_version
export Sector, Commodity, Consumer, Production, DemandFunction, Endowment, Input, Output, Parameter, Demand
export value, set_value, set_fixed!
export @parameter, @sector, @commodity, @consumer, @production, @demand

include("model.jl")
include("macros.jl")
include("build/build_helpers.jl")
include("build/build_variables.jl")
include("build/build_implicitconstraints.jl")
include("build/build_zeroprofit.jl")
include("build/build_marketclearance.jl")
include("build/build_incomebalance.jl")
include("build/build_startvalues_bounds.jl")
include("build/build.jl")
include("algebraic_wrapper.jl")

end