text:
module NuclearRepulsion

using TypesParticles

export calcnuclrepul

"""
    Use: calcnuclrepul(system:SystemOfAtoms)

Input:
Results:

Description: returns the total energy from the nuclear-nuclear repulsion. This is a direct
sum calculation, as opposed to Ewald summation or particle mesh solvers.

TODO: return energy in units of eV or Ha

""" function calcnuclrepul(system::SystemOfAtoms)
    energy = 0.00e0;
    if system.natoms == 1
        return energy
    else
        for i=1:system.natoms
            atomi = system.atoms[i]
            for j=i+1:system.natoms;
                atomj = system.atoms[j];
                mag_rij = sqrt((atomj.x-atomi.x)^2 +
                              (atomj.y-atomi.y)^2 +
                              (atomj.z-atomi.z)^2);
                energy = energy + atomi.Z*atomj.Z/mag_rij;
            end
        end #system.natoms
        return energy
    end
end #get_nucl_nucl

end #NuclearRepulsion


text:
@testset "Parallel.Random Vertex Cover" begin

    g3 = StarGraph(5)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g3)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            @test (length(c)== 2 && (c[1] == 1 || c[2] == 1))
        end
    end

    g4 = CompleteGraph(5)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g4)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            @test length(c)== 4 #All except one vertex
        end
    end

    g5 = PathGraph(5)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g5)
            c = @inferred(Parallel.vertex_cover(g, 4, RandomVertexCover(); parallel=parallel))
            sort!(c)
            @test (c == [1, 2, 3, 4] || c == [1, 2, 4, 5] || c == [2, 3, 4, 5])
        end
    end
end

ext: jl


text:
"""
    estimating_function_template(nobs::Function,
                                 ef_contribution::Function)

[Composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) for defining an `estimating_function_template`.

Arguments
===
+ `nobs`: a function of `data` that computes the number of observations of the particular data type,
+ `observation`: a function of `data` and an indexing value `i` which extracts the `i`th observation from the data
+ `ef_contribution`: a function of the parameters `theta`, the `data` and the observation index `i` that returns a vector of length `length(theta)`.

Result
===
An `estimating_function_template` object with fields `nobs` and `obj_contributions`.
"""
struct estimating_function_template
    nobs::Function
    observation::Function
    ef_contribution::Function
end

"""
    estimating_function(theta::Vector,
                        data::Any,
                        template::estimating_function_template,
                        br::Bool = false,
                        concentrate::Vector{Int64} = Vector{Int64}())

Evaluate a vector of estimating functions at `theta` by adding up all contributions in `data`, according to an [`estimating_function_template`](@ref).

Arguments
===
+ `theta`: a `Vector` of parameter values at which to evaluate the estimating functions
+ `data`: typically an object of [composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) with all the data required to compute the `estimating_function`.
+ `template`: an [`estimating_function_template`](@ref) object.
+ `br`: a `Bool`. If `false` (default), the estimating functions is constructed by adding up all contributions in
`data`, according to [`estimating_function_template`](@ref), before it is evaluated at `theta`. If `true` then the empirical bias-reducing adjustments in [Kosmidis & Lunardon, 2020](http://arxiv.org/abs/2001.03786) are computed and added to the estimating functions.
+ `concentrate`: a `Vector{Int64}`; if specified, empirical bias-reducing adjustments are added only to the subset of estimating functions indexed by `concentrate`. The default is to add empirical bias-reducing adjustments to all estimating functions.

Result
===
A `Vector`.

Details
===
`data` can be used to pass additional constants other than the actual data to the objective.
"""
function estimating_function(theta::Vector,
                             data::Any,
                             template::estimating_function_template,
                             br::Bool=false,
                             concentrate::Vector{Int64}=Vector{Int64}())
    p = length(theta)
    n_obs = template.nobs(data)
    ef = zeros(p)
    for i in 1:n_obs
        obs_i = template.observation(data, i)
        ef += template.ef_contribution(theta, obs_i)
    end
    if (br)
        quants = ef_quantities(theta, data, template, br, concentrate)
        ef + quants[1]
    else
        ef
    end
end

"""
    get_estimating_function(data::Any,
                            template::estimating_function_template,
                            br::Bool = false,
                            concentrate::Vector{Int64} = Vector{Int64}(),
                            regularizer::Any = Vector{Int64}())

Construct the estimating functions by adding up all contributions in the `data` according to [`estimating_function_template`](@ref).

Arguments
===
+ `data`: typically an object of [composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) with all the data required to compute the `estimating_function`.
+ `template`: an [`estimating_function_template`](@ref) object.
+ `br`: a `Bool`. If `false` (default), the estimating functions is constructed by adding up all contributions in
`data`, according to [`estimating_function_template`](@ref), before it is evaluated at `theta`. If `true` then the empirical bias-reducing adjustments in [Kosmidis & Lunardon, 2020](http://arxiv.org/abs/2001.03786) are computed and added to the estimating functions.
+ `concentrate`: a `Vector{Int64}`; if specified, empirical bias-reducing adjustments are added only to the subset of estimating functions indexed by `concentrate`. The default is to add empirical bias-reducing adjustments to all estimating functions.
+ `regularizer`: a function of `theta` and `data` returning a `Vector` of dimension equal to the number of the estimating functions, which is added to the (bias-reducing) estimating function; the default value will result in no regularization.

Result
===
An in-place function that stores the value of the estimating functions inferred from `template`, in a preallocated vector passed as its first argument, ready to be used withing `NLsolve.nlsolve`. This is the in-place version of [`estimating_function`](@ref) with the extra `regularizer` argument.
"""
function get_estimating_function(data::Any,
                                 template::estimating_function_template,
                                 br::Bool=false,
                                 concentrate::Vector{Int64}=Vector{Int64}(),
                                 regularizer::Any=Vector{Int64}())
    ## regulizer here has different type and default than fit, because
    ## the dimension of theta cannot be inferred
    has_regularizer = typeof(regularizer) <: Function
    function (F, theta::Vector)
        if has_regularizer
            out = estimating_function(theta, data, template, br, concentrate) + regularizer(theta, data)
        else
            out = estimating_function(theta, data, template, br, concentrate)
        end
        for i in 1:length(out)
            F[i] = out[i]
        end
    end
end

function ef_quantities(theta::Vector,
                       data::Any,
                       template::estimating_function_template,
                       adjustment::Bool=false,
                       concentrate::Vector{Int64}=Vector{Int64}())
    estfun = (pars, obs) -> template.ef_contribution(pars, obs)

    function ja_obs(eta::Vector{T}, obs) where {T <: Any}
        x = (eta, obs)
        f = (params, obs_i) -> estfun(params, obs_i)
        if !(haskey(JCACHE, T))
            jtape = ReverseDiff.compile(ReverseDiff.JacobianTape(f, x))
            JCACHE[T] = (
                jtape,
                (
                    zeros(T, (length(eta), length(eta))),
                    zeros(T, (length(obs), length(eta)))
                )
            )
        end
        jtape, y = JCACHE[T]
        return ReverseDiff.jacobian!(y, jtape, x)[1]
    end

    p = length(theta)
    n_obs = template.nobs(data)
    psi = zeros(p)
    emat = zeros(p, p)
    jmat = zeros(p, p)

    if adjustment
        u = (theta, obs) -> ForwardDiff.jacobian(eta -> ja_obs(eta, obs), theta)
        psi2 = Vector(undef, p)
        for j in 1:p
            psi2[j] = zeros(p, p)
        end
        umat = zeros(p * p, p)
        for i in 1:n_obs
            obs = template.observation(data, i)
            cpsi = estfun(theta, obs)
            psi += cpsi
            emat += cpsi * cpsi'
            jaco = ja_obs(theta, obs)
            jmat += -jaco
            umat += u(theta, obs)
            for j in 1:p
                psi2[j] += jaco[j, :] * cpsi'
            end
        end
    else
        for i in 1:n_obs
            obs = template.observation(data, i)
            cpsi = estfun(theta, obs)
            psi += cpsi
            emat += cpsi * cpsi'
            jmat += - ja_obs(theta, obs)
        end
    end

    jmat_inv = try
        inv(jmat)
    catch
        fill(NaN, p, p)
    end
    vcov = jmat_inv * (emat * jmat_inv')

    if adjustment
        Afun(j::Int64) = -tr(jmat_inv * psi2[j] + vcov * umat[j:p:(p * p - p + j), :] / 2)
        A = map(Afun, 1:p)
        ## if concentrate then redefine A
        if length(concentrate) > 0
            if any((concentrate .> p) .| (concentrate .< 1))
                error(concentrate, " should be a vector of integers between ", 1, " and ", p)
            else
                ist = concentrate
                nce = deleteat!(collect(1:p), concentrate)
                A_ist = A[ist]
                A_nce = A[nce]
                A = vcat(A_ist + inv(jmat_inv[ist, ist]) * jmat_inv[ist, nce] * A_nce,
                         zeros(length(nce)))
            end
        end
        [A, jmat_inv, emat, psi]
    else
        [vcov, jmat_inv, emat, psi]
    end
end


text:
@reexport module _Algorithm

export @forcartesian

# https://github.com/timholy/Cartesian.jl/blob/master/src/compat.jl
macro forcartesian(sym, sz, ex)
    idim = gensym()
    N = gensym()
    sz1 = gensym()
    isdone = gensym()
    quote
        if !(isempty($(esc(sz))) || prod($(esc(sz))) == 0)
            $N = length($(esc(sz)))
            $sz1 = $(esc(sz))[1]
            $isdone = false
            $(esc(sym)) = ones(Int, $N)
            while !$isdone
                $(esc(ex))
                if ($(esc(sym))[1] += 1) > $sz1
                    $idim = 1
                    while $(esc(sym))[$idim] > $(esc(sz))[$idim] && $idim < $N
                        $(esc(sym))[$idim] = 1
                        $idim += 1
                        $(esc(sym))[$idim] += 1
                    end
                    $isdone = $(esc(sym))[$N] > $(esc(sz))[$N]
                end
            end
        end
    end
end

include("binarysearch.jl")
include("cycledetection.jl")

end


text:
@doc raw"""
    flowpp(p0, ps, k, m, T,δ=0.02)

compute $-\frac{1}{\pi}\tilde{\partial_k}\Im I_{1, k}(p)$

To be noticed that, `flowpp` doesn't contains $\tilde{\partial_k}\mathcal{F}_4$,
we will consider it separately.

At $p_0=2E_{\pi,k}$, `flowpp` has a $\delta$ function contribution, we use a rectangle
function with width $\delta$ to approximate.
"""
function flowpp(p0, ps, k, m, T,δ=0.02)
    if k > ps / 2
        if p0 >= Epi(k + ps, m) + Epi(k, m)
            return 0.0
        elseif 2 * Epi(k, m) < p0 < Epi(k + ps, m) + Epi(k, m)
            return ppfun(p0, ps, k, Epi(k, m), T)
        elseif 2 * Epi(k - δ, m) <= p0 <= 2 * Epi(k, m)
            return -peak(p0, ps, m, T,δ)
        elseif p0 < 2 * Epi(k - δ, m)
            return 0.0
        end
    elseif k <= ps / 2
        if p0 >= Epi(k + ps, m) + Epi(k, m)
            return 0.0
        elseif Epi(ps, m) + Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
            return ppfun(p0, ps, k, Epi(k, m), T)
        elseif Epi(k, m) + Epi(ps - k, m) <= p0 < Epi(ps, m) + Epi(k, m)
            return ppfun(p0, ps, k, Epi(k, m), T)
        elseif p0 < Epi(k, m) + Epi(ps - k, m)
            return 0.0
        end
    end
end

@doc raw"""
    flowpm(p0, ps, k, m, T)

compute $-\frac{1}{\pi}\tilde{\partial_k}\Im I_{2, k}(p)$
"""
function flowpm(p0, ps, k, m, T)
    if k > ps / 2
        if p0 > Epi(k + ps, m) - Epi(k, m)
            return 0.0
        elseif p0 <= Epi(k + ps, m) - Epi(k, m)
            return pmfun(p0, ps, k, Epi(k, m), T)
        end
    elseif k <= ps / 2
        if p0 > Epi(k + ps, m) - Epi(k, m)
            return 0.0
        elseif Epi(k - ps, m) - Epi(k, m) < p0 <= Epi(k + ps, m) - Epi(k, m)
            return pmfun(p0, ps, k, Epi(k, m), T)
        elseif p0 <= Epi(k - ps, m) - Epi(k, m)
            return 0.0
        end
    end
end

@doc raw"""
    flowpp_intcostheqs(p0, ps, qsmax, k, m, T)

compute
```math
\begin{aligned}
&\int_0^{qsmax}\!\!dq_s q_s^2\int_{-1}^{1}\!\!d\cos\theta \tilde{\partial_k}F_1\left(\sqrt{p_s^2+q_s^2+2*p_s*q_s\cos\theta}\right)\\
&=\int_0^{qsmax}\!\!dq_s q_s^2\int_{-1}^{1}\!\!d\cos\theta \tilde{\partial_k}F_1\left(\sqrt{p_s^2+q_s^2-2*p_s*q_s\cos\theta}\right)
\end{aligned}
```

"""
function flowpp_intcostheqs(p0, ps, qsmax, k, m, T)
    if p0 <= 2 * Epi(k, m)
        return 0.0
    elseif p0 > 2 * Epi(k, m)
        Ek2 = sqrt(k^2 + p0 * (-2 * sqrt(k^2 + m) + p0))
        Ek = sqrt(k^2 + m)
        if ps >= k + Ek2 + qsmax
            return 0.0
        elseif k + Ek2 <= ps < k + Ek2 + qsmax
            return ppfuncostheps1(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 + k - qsmax <= ps < k + Ek2
            return ppfuncostheps2(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 <= ps < Ek2 + k - qsmax
            return ppfuncostheps3(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 - k + qsmax <= ps < Ek2
            return ppfuncostheps4(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 - k <= ps < Ek2 - k + qsmax && qsmax <= ps
            return ppfuncostheps5(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 - k <= ps < qsmax && ps >= -Ek2 + k + qsmax
            return ppfuncostheps6(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 - k <= ps < qsmax && ps < -Ek2 + k + qsmax
            return ppfuncostheps7(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 - k - qsmax < ps < Ek2 - k &&
               ps <= (Ek2 - k) / 2 &&
               Ek2 - k + ps <= qsmax
            return ppfuncostheps8(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 - k - qsmax < ps < Ek2 - k &&
               ps <= (Ek2 - k) / 2 &&
               Ek2 - k + ps > qsmax
            return ppfuncostheps9(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 - k - qsmax < ps < Ek2 - k &&
               ps > (Ek2 - k) / 2 &&
               Ek2 - k + ps <= qsmax
            return ppfuncostheps10(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif Ek2 - k - qsmax < ps < Ek2 - k &&
               ps > (Ek2 - k) / 2 &&
               Ek2 - k + ps > qsmax
            return ppfuncostheps11(p0, ps, qsmax, k, Ek, Ek2, T)
        elseif ps <= Ek2 - k - qsmax
            return 0.0
        end
    end
end

@doc raw"""
    flowpm_intcostheqs(p0, ps, qsmax, k, m, T)

compute
```math
\begin{aligned}
&\int_0^{qsmax}\!\!dq_s q_s^2\int_{-1}^{1}\!\!d\cos\theta \tilde{\partial_k}F_2\left(\sqrt{p_s^2+q_s^2+2p_sq_s\cos\theta}\right)\\
&=\int_0^{qsmax}\!\!dq_s q_s^2\int_{-1}^{1}\!\!d\cos\theta \tilde{\partial_k}F_2\left(\sqrt{p_s^2+q_s^2-2*p_s*q_s\cos\theta}\right)
\end{aligned}
```
"""
function flowpm_intcostheqs(p0, ps, qsmax, k, m, T)
    Ek1 = sqrt(k^2 + p0 * (2 * sqrt(k^2 + m) + p0))
    Ek = sqrt(k^2 + m)
    # println(ps-qsmax,",",Ek1 - k)
    if ps > Ek1 + k + qsmax
        return 0.0
    elseif Ek1 + k < ps < Ek1 + k + qsmax
        #check
        return pmfuncostheps6(p0, ps, qsmax, k, Ek, Ek1, T)
    elseif Ek1 + k - qsmax < ps < Ek1 + k
        return pmfuncostheps8(p0, ps, qsmax, k, Ek, Ek1, T)
    elseif Ek1 < ps < Ek1 + k - qsmax
        return pmfuncostheps7(p0, ps, qsmax, k, Ek, Ek1, T)
    elseif Ek1 - k + qsmax <= ps < Ek1
        return pmfuncostheps9(p0, ps, qsmax, k, Ek, Ek1, T)
    elseif Ek1 - k <= ps < Ek1 - k + qsmax
        if p0 >= -sqrt(k^2 + m) + sqrt(m + (k + qsmax)^2)
            return pmfuncostheps10(p0, ps, qsmax, k, Ek, Ek1, T)
        elseif p0 < -sqrt(k^2 + m) + sqrt(m + (k + qsmax)^2)
            if ps <= qsmax / 2
                return pmfuncostheps11(p0, ps, qsmax, k, Ek, Ek1, T)
            elseif qsmax >= ps > qsmax / 2 && qsmax - ps >= Ek1 - k
                return pmfuncostheps11(p0, ps, qsmax, k, Ek, Ek1, T)
            elseif qsmax >= ps > qsmax / 2 && qsmax - ps < Ek1 - k
                return pmfuncostheps12(p0, ps, qsmax, k, Ek, Ek1, T)
            elseif ps > qsmax
                return pmfuncostheps12(p0, ps, qsmax, k, Ek, Ek1, T)
                # elseif ps > qsmax &&ps-qsmax > Ek1 - k
                #     println("su2")
                #     return pmfuncostheps13(p0, ps, qsmax, k, Ek, Ek1, T)
            end
        end
    elseif ps < Ek1 - k
        if ps <= -k + Ek1 - qsmax
            return 0.0
        elseif ps > -k + Ek1 - qsmax && ps >= qsmax
            return pmfuncostheps14(p0, ps, qsmax, k, Ek, Ek1, T)
        elseif ps > -k + Ek1 - qsmax &&
               ps < qsmax &&
               ps <= (Ek1 - k) / 2 &&
               Ek1 - k + ps <= qsmax
            return pmfuncostheps15(p0, ps, qsmax, k, Ek, Ek1, T)
        elseif ps > -k + Ek1 - qsmax &&
               ps < qsmax &&
               ps <= (Ek1 - k) / 2 &&
               Ek1 - k + ps > qsmax
            return pmfuncostheps16(p0, ps, qsmax, k, Ek, Ek1, T)
        elseif ps > -k + Ek1 - qsmax &&
               ps < qsmax &&
               ps > (Ek1 - k) / 2 &&
               Ek1 - k + ps <= qsmax
            return pmfuncostheps17(p0, ps, qsmax, k, Ek, Ek1, T)
        elseif ps > -k + Ek1 - qsmax &&
               ps < qsmax &&
               ps > (Ek1 - k) / 2 &&
               Ek1 - k + ps > qsmax
            return pmfuncostheps18(p0, ps, qsmax, k, Ek, Ek1, T)
        end
    end
end

#integrate delta function in F1 , we integrate out qs, cos(θ) and k' this will be used in Im part calculation
function deltasumkfix(p0, ps, k, T, Npi, IRScale, UVScale, mfun, lamfun)
    #find the location of k0 where p0==2Epi(k,m)
    deltaf(x) = 2 * Epi(x, mfun(x)) - p0
    if deltaf(IRScale) * deltaf(UVScale) >= 0
        return 0.0
    else
        k0 = find_zero(deltaf, (IRScale, UVScale))
        # println(" k0=", k0)
        # k0 should lies between k~Λ, when k>k0 you will get 0
        if k > k0
            return 0.0
            #δ function only appears in p<2k so we have the following division
        elseif k <= k0
            if k <= 2 * k0 - ps
                if ps > k
                    return (
                        k^3 *
                        k0 *
                        (k^2 + 5 * ps * (-2 * k0 + ps)) *
                        coth(p0 / (4 * T)) *
                        lamfun(k0)^2
                    ) / (
                        30 *
                        p0 *
                        pi^2 *
                        ps *
                        abs(2 * k0 + derivative(mfun, k0))
                    )
                elseif ps <= k
                    return -1 / 120 *(
                        k0 *
                        (5 * k^3 * (-3 * k + 8 * k0) - 10 * k^2 * ps^2 + ps^4) *
                        coth(p0 / (4 * T)) *
                        lamfun(k0)^2
                    ) / (p0 * pi^2 * abs(2 * k0 + derivative(mfun, k0)))
                end
            elseif max(ps - k, 0.0) > 2 * k0
                return 0.0
            elseif 2*k < ps < 2 * k0 + k
                # println("locate end")
                return (k0*(k+2*k0-ps)^3*(6*k0^2+(k-ps)*(4*k+ps)-k0*(9*k+ps))*coth(p0/(4*T))*lamfun(k0)^2)/(240*p0*pi^2*ps*abs(2*k0+derivative(mfun,k0)))
            elseif 2 * k0 - k<ps < 2 * k0
                return (k0*(k+2*k0-ps)^3*(6*k0^2+(k-ps)*(4*k+ps)-k0*(9*k+ps))*coth(p0/(4*T))*lamfun(k0)^2)/(240*p0*pi^2*ps*abs(2*k0+derivative(mfun,k0)))
            end
        end
    end
end

function deltasumkfixSmooth(p0, ps, k, T, Npi, IRScale, UVScale, mfun, lamfun,slop=10.0)
    #find the location of k0 where p0==2Epi(k,m)
    deltaf(x) = 2 * Epi(x, mfun(x)) - p0
    if deltaf(IRScale) * deltaf(UVScale) >= 0
        return 0.0
    else
        k0 = find_zero(deltaf, (IRScale, UVScale))
        # println(" k0=", k0)
        # k0 should lies between k~Λ, when k>k0 you will get 0
        if k <= 2 * k0 - ps || k>k0
            if ps > k
                return 0.5*(tanh(slop*(k0-k))+1)*(
                    k^3 *
                    k0 *
                    (k^2 + 5 * ps * (-2 * k0 + ps)) *
                    coth(p0 / (4 * T)) *
                    lamfun(k0)^2
                ) / (
                    30 *
                    p0 *
                    pi^2 *
                    ps *
                    abs(2 * k0 + derivative(mfun, k0))
                )
            elseif ps <= k
                return -1 / 120 *0.5*(tanh(slop*(k0-k))+1)* (
                    k0 *
                    (5 * k^3 * (-3 * k + 8 * k0) - 10 * k^2 * ps^2 + ps^4) *
                    coth(p0 / (4 * T)) *
                    lamfun(k0)^2
                ) / (p0 * pi^2 * abs(2 * k0 + derivative(mfun, k0)))
            end
        elseif max(ps - k, 0.0) > 2 * k0
            return 0.0
        elseif 2 * k0 - k < ps < 2 * k0 + k
            # println("locate end")
            return 0.5*(tanh(slop*(k0-k))+1)*(k0*(k+2*k0-ps)^3*(6*k0^2+(k-ps)*(4*k+ps)-k0*(9*k+ps))*coth(p0/(4*T))*lamfun(k0)^2)/(240*p0*pi^2*ps*abs(2*k0+derivative(mfun,k0)))
        end
    end
end

function deltasumkfix_dp0(p0, ps, k, T, Npi, IRScale, UVScale, mfun, lamfun)
    #find the location of k0 where p0==2Epi(k,m)
    deltaf(x) = 2 * Epi(x, mfun(x)) - p0
    if deltaf(IRScale) * deltaf(UVScale) >= 0
        return 0.0
    else
        k0 = find_zero(deltaf, (IRScale, UVScale))
        # println(" k0=", k0)
        # k0 should lies between k~Λ, when k>k0 you will get 0
        if k > k0
            return 0.0
            #δ function only appears in p<2k so we have the following division
        elseif k <= k0
            if ps <= 2 * k0 - k
                if ps > k
                    return -1 / 120 * (
                        k^3 *
                        lamfun(k0) *
                        (
                            4 *
                            k0 *
                            (k^2 + 5 * ps * (-2 * k0 + ps)) *
                            T *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            coth(p0 / (4 * T)) *
                            (2 * k0 + derivative(mfun, k0)) *
                            lamfun(k0) +
                            k0 *
                            p0 *
                            (k^2 + 5 * ps * (-2 * k0 + ps)) *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            csch(p0 / (4 * T))^2 *
                            (2 * k0 + derivative(mfun, k0)) *
                            lamfun(k0) -
                            8 *
                            k0 *
                            p0 *
                            (k^2 + 5 * ps * (-2 * k0 + ps)) *
                            T *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            coth(p0 / (4 * T)) *
                            derivative(lamfun, k0) *
                            sqrt(k0^2 + mfun(k0)) +
                            40 *
                            k0 *
                            p0 *
                            ps *
                            T *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            coth(p0 / (4 * T)) *
                            lamfun(k0) *
                            sqrt(k0^2 + mfun(k0)) -
                            4 *
                            p0 *
                            (k^2 + 5 * ps * (-2 * k0 + ps)) *
                            T *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            coth(p0 / (4 * T)) *
                            lamfun(k0) *
                            sqrt(k0^2 + mfun(k0)) +
                            4 *
                            k0 *
                            p0 *
                            (k^2 + 5 * ps * (-2 * k0 + ps)) *
                            T *
                            coth(p0 / (4 * T)) *
                            dabs(2 * k0 + derivative(mfun, k0)) *
                            (2 + derivative(mfun, k0, nu = 2)) *
                            lamfun(k0) *
                            sqrt(k0^2 + mfun(k0))
                        )
                    ) / (
                        p0^2 *
                        pi^2 *
                        ps *
                        T *
                        abs(2 * k0 + derivative(mfun, k0))^2 *
                        (2 * k0 + derivative(mfun, k0))
                    )
                elseif ps <= k
                    return (
                        lamfun(k0) * (
                            4 *
                            k0 *
                            (
                                -15 * k^4 + 40 * k^3 * k0 - 10 * k^2 * ps^2 +
                                ps^4
                            ) *
                            T *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            coth(p0 / (4 * T)) *
                            (2 * k0 + derivative(mfun, k0)) *
                            lamfun(k0) +
                            k0 *
                            p0 *
                            (
                                -15 * k^4 + 40 * k^3 * k0 - 10 * k^2 * ps^2 +
                                ps^4
                            ) *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            csch(p0 / (4 * T))^2 *
                            (2 * k0 + derivative(mfun, k0)) *
                            lamfun(k0) -
                            8 *
                            k0 *
                            p0 *
                            (
                                -15 * k^4 + 40 * k^3 * k0 - 10 * k^2 * ps^2 +
                                ps^4
                            ) *
                            T *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            coth(p0 / (4 * T)) *
                            derivative(lamfun, k0) *
                            sqrt(k0^2 + mfun(k0)) -
                            160 *
                            k^3 *
                            k0 *
                            p0 *
                            T *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            coth(p0 / (4 * T)) *
                            lamfun(k0) *
                            sqrt(k0^2 + mfun(k0)) -
                            4 *
                            p0 *
                            (
                                -15 * k^4 + 40 * k^3 * k0 - 10 * k^2 * ps^2 +
                                ps^4
                            ) *
                            T *
                            abs(2 * k0 + derivative(mfun, k0)) *
                            coth(p0 / (4 * T)) *
                            lamfun(k0) *
                            sqrt(k0^2 + mfun(k0)) +
                            4 *
                            k0 *
                            p0 *
                            (
                                -15 * k^4 + 40 * k^3 * k0 - 10 * k^2 * ps^2 +
                                ps^4
                            ) *
                            T *
                            coth(p0 / (4 * T)) *
                            dabs(2 * k0 + derivative(mfun, k0)) *
                            (2 + derivative(mfun, k0, nu = 2)) *
                            lamfun(k0) *
                            sqrt(k0^2 + mfun(k0))
                        )
                    ) / (
                        480 *
                        p0^2 *
                        pi^2 *
                        T *
                        abs(2 * k0 + derivative(mfun, k0))^2 *
                        (2 * k0 + derivative(mfun, k0))
                    )
                end
            elseif max(ps - k, 0.0) > 2 * k0
                return 0.0
            elseif 2 * k0 - k < ps < 2 * k0 + k
                # println("locate end")
                return (
                    (k + 2 * k0 - ps)^2 *
                    lamfun(k0) *
                    (
                        -4 *
                        k0 *
                        (k + 2 * k0 - ps) *
                        (
                            6 * k0^2 + (k - ps) * (4 * k + ps) -
                            k0 * (9 * k + ps)
                        ) *
                        T *
                        abs(2 * k0 + derivative(mfun, k0)) *
                        coth(p0 / (4 * T)) *
                        (2 * k0 + derivative(mfun, k0)) *
                        lamfun(k0) -
                        k0 *
                        p0 *
                        (k + 2 * k0 - ps) *
                        (
                            6 * k0^2 + (k - ps) * (4 * k + ps) -
                            k0 * (9 * k + ps)
                        ) *
                        abs(2 * k0 + derivative(mfun, k0)) *
                        csch(p0 / (4 * T))^2 *
                        (2 * k0 + derivative(mfun, k0)) *
                        lamfun(k0) +
                        8 *
                        k0 *
                        p0 *
                        (k + 2 * k0 - ps) *
                        (
                            6 * k0^2 + (k - ps) * (4 * k + ps) -
                            k0 * (9 * k + ps)
                        ) *
                        T *
                        abs(2 * k0 + derivative(mfun, k0)) *
                        coth(p0 / (4 * T)) *
                        derivative(lamfun, k0) *
                        sqrt(k0^2 + mfun(k0)) +
                        4 *
                        k0 *
                        p0 *
                        (k + 2 * k0 - ps) *
                        (-9 * k + 12 * k0 - ps) *
                        T *
                        abs(2 * k0 + derivative(mfun, k0)) *
                        coth(p0 / (4 * T)) *
                        lamfun(k0) *
                        sqrt(k0^2 + mfun(k0)) +
                        24 *
                        k0 *
                        p0 *
                        (
                            6 * k0^2 + (k - ps) * (4 * k + ps) -
                            k0 * (9 * k + ps)
                        ) *
                        T *
                        abs(2 * k0 + derivative(mfun, k0)) *
                        coth(p0 / (4 * T)) *
                        lamfun(k0) *
                        sqrt(k0^2 + mfun(k0)) +
                        4 *
                        p0 *
                        (k + 2 * k0 - ps) *
                        (
                            6 * k0^2 + (k - ps) * (4 * k + ps) -
                            k0 * (9 * k + ps)
                        ) *
                        T *
                        abs(2 * k0 + derivative(mfun, k0)) *
                        coth(p0 / (4 * T)) *
                        lamfun(k0) *
                        sqrt(k0^2 + mfun(k0)) -
                        4 *
                        k0 *
                        p0 *
                        (k + 2 * k0 - ps) *
                        (
                            4 * k^2 + 6 * k0^2 - k0 * ps - ps^2 -
                            3 * k * (3 * k0 + ps)
                        ) *
                        T *
                        coth(p0 / (4 * T)) *
                        dabs(2 * k0 + derivative(mfun, k0)) *
                        (2 + derivative(mfun, k0, nu = 2)) *
                        lamfun(k0) *
                        sqrt(k0^2 + mfun(k0))
                    )
                ) / (
                    960 *
                    p0^2 *
                    pi^2 *
                    ps *
                    T *
                    abs(2 * k0 + derivative(mfun, k0))^2 *
                    (2 * k0 + derivative(mfun, k0))
                )
            end
        end
    end
end

function deltasumkfix_dp0_deltafun(
    p0,
    ps,
    k0,
    k,
    T,
    Npi,
    IRScale,
    UVScale,
    mfun,
    lamfun,
)
    if k <= 2 * k0 - ps
        if ps > k
            return (
                k^3 *
                k0 *
                (k^2 + 5 * ps * (-2 * k0 + ps)) *
                coth(p0 / (4 * T)) *
                lamfun(k0)^2
            ) / (30 * p0 * pi^2 * ps * abs(2 * k0 + derivative(mfun, k0)))
        elseif ps <= k
            return -1 / 120 * (
                k0 *
                (5 * k^3 * (-3 * k + 8 * k0) - 10 * k^2 * ps^2 + ps^4) *
                coth(p0 / (4 * T)) *
                lamfun(k0)^2
            ) / (p0 * pi^2 * abs(2 * k0 + derivative(mfun, k0)))
        end
    elseif max(ps - k, 0.0) > 2 * k0
        return 0.0
    elseif 2 * k0 - k < ps < 2 * k0 + k
        # println("locate end")
        return (
            k0 *
            (k + 2 * k0 - ps)^3 *
            (6 * k0^2 + (k - ps) * (4 * k + ps) - k0 * (9 * k + ps)) *
            coth(p0 / (4 * T)) *
            lamfun(k0)^2
        ) / (240 * p0 * pi^2 * ps * abs(2 * k0 + derivative(mfun, k0)))
    end
end

#integrate type1 delta function in F1 , we integrate out qs, cos(θ) and q0' this will be used in Re part calculation
function delta1_intcosthqs(p0, ps, qsmax, k, m, T,δk=0.02)
    #δ function only appears in p<2k so we have the following division
    if ps >= 2 * k + qsmax
        return 0.0
        #for Principal Value integration we ommit k0-δ ~ k0+δ
    elseif 2 * k + qsmax > ps >= 2 * k
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return -1 / 120 * (
                k *
                (2 * k - ps + qsmax)^3 *
                (
                    6 * k^2 - (ps - qsmax) * (ps + 4 * qsmax) -
                    k * (ps + 9 * qsmax)
                ) *
                coth(Epi(k, m) / (2 * T)) *
                (p0 + Epi(k, m))
            ) / (p0^2 * pi^2 * ps * Epi(k, m) * (p0 + 2 * Epi(k, m)))
        else
            return (
                k *
                (2 * k - ps + qsmax)^3 *
                (
                    6 * k^2 - (ps - qsmax) * (ps + 4 * qsmax) -
                    k * (ps + 9 * qsmax)
                ) *
                coth(Epi(k, m) / (2 * T))
            ) / (120 * pi^2 * ps * Epi(k, m) * (-p0^2 + 4 * Epi(k, m)^2))
        end
    elseif 2 * k > ps >= 2 * k - qsmax
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return (
                k *
                qsmax^2 *
                (
                    10 * (-2 * k + ps)^2 * (k + ps) +
                    20 * (2 * k - ps) * ps * qsmax +
                    15 * (-k + ps) * qsmax^2 - 4 * qsmax^3
                ) *
                coth(Epi(k, m) / (2 * T)) *
                (p0 + Epi(k, m))
            ) / (120 * p0^2 * pi^2 * ps * Epi(k, m) * (p0 + 2 * Epi(k, m)))
        else
            return -1 / 120 * (
                k *
                qsmax^2 *
                (
                    -10 * (-2 * k + ps)^2 * (k + ps) +
                    20 * ps * (-2 * k + ps) * qsmax +
                    15 * (k - ps) * qsmax^2 +
                    4 * qsmax^3
                ) *
                coth(Epi(k, m) / (2 * T))
            ) / (pi^2 * ps * Epi(k, m) * (p0^2 - 4 * Epi(k, m)^2))
        end
    elseif 2 * k - qsmax > ps >= qsmax
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return (
                k *
                qsmax^3 *
                (10 * k * ps - 5 * ps^2 - qsmax^2) *
                coth(Epi(k, m) / (2 * T)) *
                (p0 + Epi(k, m))
            ) / (15 * p0^2 * pi^2 * ps * Epi(k, m) * (p0 + 2 * Epi(k, m)))
        else
            return (
                k *
                qsmax^3 *
                (5 * ps * (-2 * k + ps) + qsmax^2) *
                coth(Epi(k, m) / (2 * T))
            ) / (15 * pi^2 * ps * Epi(k, m) * (-p0^2 + 4 * Epi(k, m)^2))
        end
    elseif ps < qsmax
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return (
                k *
                (
                    ps^4 - 10 * ps^2 * qsmax^2 +
                    5 * (8 * k - 3 * qsmax) * qsmax^3
                ) *
                coth(Epi(k, m) / (2 * T)) *
                (p0 + Epi(k, m))
            ) / (60 * p0^2 * pi^2 * Epi(k, m) * (p0 + 2 * Epi(k, m)))
        else
            return (
                k *
                (
                    ps^4 - 10 * ps^2 * qsmax^2 +
                    5 * (8 * k - 3 * qsmax) * qsmax^3
                ) *
                coth(Epi(k, m) / (2 * T))
            ) / (60 * pi^2 * Epi(k, m) * (p0^2 - 4 * Epi(k, m)^2))
        end
    end
end

function delta1_intcosthqs_Pv(p0, ps, qsmax, k, m, T)
    #δ function only appears in p<2k so we have the following division
    if ps >= 2 * k + qsmax
        return 0.0
        #for Principal Value integration we ommit k0-δ ~ k0+δ
    elseif 2 * k + qsmax > ps >= 2 * k
        return (
            k *
            (2 * k - ps + qsmax)^3 *
            (
                6 * k^2 - (ps - qsmax) * (ps + 4 * qsmax) -
                k * (ps + 9 * qsmax)
            ) *
            coth(Epi(k, m) / (2 * T))
        ) / (120 * pi^2 * ps * Epi(k, m))
    elseif 2 * k > ps >= 2 * k - qsmax
        return 1 / 120 * (
            k *
            qsmax^2 *
            (
                -10 * (-2 * k + ps)^2 * (k + ps) +
                20 * ps * (-2 * k + ps) * qsmax +
                15 * (k - ps) * qsmax^2 +
                4 * qsmax^3
            ) *
            coth(Epi(k, m) / (2 * T))
        ) / (pi^2 * ps * Epi(k, m))
    elseif 2 * k - qsmax > ps >= qsmax
        return (
            k *
            qsmax^3 *
            (5 * ps * (-2 * k + ps) + qsmax^2) *
            coth(Epi(k, m) / (2 * T))
        ) / (15 * pi^2 * ps * Epi(k, m))
    elseif ps < qsmax
        return -(
            k *
            (
                ps^4 - 10 * ps^2 * qsmax^2 +
                5 * (8 * k - 3 * qsmax) * qsmax^3
            ) *
            coth(Epi(k, m) / (2 * T))
        ) / (60 * pi^2 * Epi(k, m))
    end
end

#integrate type2 delta function in F1 , we integrate out qs, cos(θ) and q0' this will be used in Re part calculation
function delta2_intcosthqs(p0, ps, qsmax, k, m, T,δk=0.02)
    #δ function only appears in p<2k so we have the following division
    Ek=Epi(k,m)
    if ps >= 2 * k + qsmax
        return 0.0
        #for Principal Value integration we ommit k0-δ ~ k0+δ
    elseif 2 * k + qsmax > ps >= 2 * k
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return delta2funcosthqs1(p0, ps, qsmax, k, Ek, T)
        end
    elseif 2 * k > ps >= 2 * k - qsmax
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return delta2funcosthqs2(p0, ps, qsmax, k, Ek, T)
        end
    elseif 2 * k - qsmax > ps >= qsmax
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return delta2funcosthqs3(p0, ps, qsmax, k, Ek, T)
        end
    elseif ps < qsmax
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return delta2funcosthqs4(p0, ps, qsmax, k, Ek, T)
        end
    end
end

function delta2Smooth_intcosthqs(p0, ps, qsmax, k, m, T,δk=0.0,δ=0.01)
    #δ function only appears in p<2k so we have the following division
    Ek=Epi(k,m)
    if ps >= 2 * k + qsmax
        return 0.0
        #for Principal Value integration we ommit k0-δ ~ k0+δ
    elseif 2 * k + qsmax > ps >= 2 * k
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return delta2Smoothfuncosthqs1(p0, ps, qsmax, k, Ek, T,δ)
        end
    elseif 2 * k > ps >= 2 * k - qsmax
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return delta2Smoothfuncosthqs2(p0, ps, qsmax, k, Ek, T,δ)
        end
    elseif 2 * k - qsmax > ps >= qsmax
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return delta2Smoothfuncosthqs3(p0, ps, qsmax, k, Ek, T,δ)
        end
    elseif ps < qsmax
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return delta2Smoothfuncosthqs4(p0, ps, qsmax, k, Ek, T,δ)
        end
    end
end

function F4_intcosthqs(p0, ps, qsmax, k, m, T)
    #δ function only appears in p<2k so we have the following division
    Ek=Epi(k,m)
    if ps >= 2 * k + qsmax
        return 0.0
    elseif 2 * k + qsmax > ps >= 2 * k
        return F4funcosthqs1(p0, ps, qsmax, k, Ek, T)
    elseif 2 * k > ps >= 2 * k - qsmax
        return F4funcosthqs2(p0, ps, qsmax, k, Ek, T)
    elseif 2 * k - qsmax > ps >= qsmax
        return F4funcosthqs3(p0, ps, qsmax, k, Ek, T)
    elseif ps < qsmax
        return F4funcosthqs4(p0, ps, qsmax, k, Ek, T)
    end
end

function F4tilde_intcosthqs(p0, ps, k, m, T,δk=0.02)
    Ek=Epi(k,m)
    if ps >= 2 * k + k
        return 0.0
        #for Principal Value integration we ommit k0-δ ~ k0+δ
    elseif 2 * k + k > ps >= 2 * k
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return F4tildefuncosthqs1(p0, ps, k, Ek, T)
        end
    elseif 2 * k > ps >= 2 * k - k
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return F4tildefuncosthqs2(p0, ps, k, Ek, T)
        end
    elseif ps < k
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return F4tildefuncosthqs4(p0, ps, k, Ek, T)
        end
    end
end

function F4tildeSmooth_intcosthqs(p0, ps,qsmax, k, m, T,δk=0.1,δ=0.05)
    Ek=Epi(k,m)
    if ps >= 2 * k + k
        return 0.0
        #for Principal Value integration we ommit k0-δ ~ k0+δ
    elseif 2 * k + k > ps >= 2 * k
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return F4tildeSmoothfuncosthqs1(p0, ps,qsmax, k, Ek, T,δ)
        end
    elseif 2 * k > ps >= 2 * k - k
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return F4tildeSmoothfuncosthqs2(p0, ps, qsmax,k, Ek, T,δ)
        end
    elseif ps < k
        if p0 + δk > 2 * Epi(k, m) > p0 - δk
            return 0.0
        else
            return F4tildeSmoothfuncosthqs4(p0, ps, qsmax,k, Ek, T,δ)
        end
    end
end

function F4tildek_intcosthqs(p0, ps, k, m, T, UVScale, δk = 0.02)
    F4tilde_intcosthqs(p0, ps, k, m, T, δk) -
    F4tilde_intcosthqs(p0, ps, UVScale, m, T, δk)
end

function F4tildekSmooth_intcosthqs(p0, ps, k, m, T, UVScale, δk=0.1,δ=0.05)
    F4tildeSmooth_intcosthqs(p0, ps,k, k, m, T,δk,δ) -
    F4tildeSmooth_intcosthqs(p0, ps,k, UVScale, m, T,δk,δ)
end


text:
using RouletteWheels
using FactCheck
using Distributions
using Compat

const algos = [LinearWalk, BisectingSearch, StochasticAcceptance]

function in_ci(wheel, n)
    ps = collect(wheel)
    ps /= sum(ps)
    tallies = rand_tally(wheel, n)

    for (p, tally) in zip(ps, tallies)
        (l_b, u_b) = quantile(Binomial(n, p), [0.00001, 0.99999])
        if u_b < tally < l_b
            return false
        end
    end

    return true
end

facts("A RouletteWheel") do
    context("iterates over the underlying frequency") do
        for algo in algos
            @fact collect(algo([1,9,9,4])) --> [1, 9, 9, 4]
        end
    end

    context("samples in proportion to the given frequency") do
        freqs = [5,4,9,2,4]
        for algo in algos
            @fact in_ci(algo(freqs), 10000) --> true
        end
    end

    context("samples in proportion to the given proportion") do
        ps = [5,4,9,2,4]
        ps /= sum(ps)
        for algo in algos
            @fact in_ci(algo(ps), 10000) --> true
        end
    end

    context("allows for push! given a normalize!") do
        freqs = [5,4,9,2,4]
        for algo in algos
            wheel = algo(freqs)

            push!(wheel, 1)
            @fact length(wheel) --> 6
            normalize!(wheel)
            @fact wheel[6] --> 1

            @fact in_ci(wheel, 10000) --> true
        end
    end
end

facts("The WheelFromDict wrapper") do
    context("samples underlying values instead of an index") do
        d = @compat Dict{Symbol, Int}(:red => 1, :green => 5, :blue => 20)
        wheel = WheelFromDict(d)
        @fact length(wheel) --> 3

        for (k, v) in wheel
            @fact k ∈ keys(d) --> true
            @fact v --> d[k]
        end

        sampled_d = rand_dict(wheel, 10000)
        @fact sampled_d[:red] < sampled_d[:green] < sampled_d[:blue] --> true

        wheel[:gold] = 90
        normalize!(wheel)
        sampled_d = rand_dict(wheel, 10000)
        @fact sampled_d[:red] < sampled_d[:green] < sampled_d[:blue] --> true
        @fact sampled_d[:blue] < sampled_d[:gold] --> true
    end
end

facts("The select_fastest function") do
    context("selects the fastest wheel over some proportions") do
        fastest, all_timings = select_fastest(1:3) do wheel
            for _ in 1:1000
                rand(wheel)
            end
        end
        @fact length(all_timings) --> 3
        @fact all_timings[fastest] --> minimum(values(all_timings))
    end
end


text:
mutable struct Callbacks
    cbs::Vector
    runner::CallbackRunner
    graph::SimpleDiGraph
    initialized::Bool
end
Callbacks(cbs, runner = LinearRunner()) = Callbacks(cbs, runner, callbackgraph(cbs), false)

"""
    mutable struct BatchState

Stores data of the last processed batch.

# Fields

- `xs`
- `ys`
- `ŷs`
- `loss`
- `grads`

(!) If used in callbacks, some fields may be `nothing` as
they are reset after every step.
"""
@with_kw mutable struct BatchState
    xs = nothing
    ys = nothing
    ŷs = nothing
    loss = nothing
    grads = nothing
end

mutable struct Learner
    model
    data
    optimizer
    lossfn
    params
    batch::BatchState
    callbacks::Callbacks
    cbstate::PropDict
end

"""
    Learner(model, data, optimizer, lossfn, [callbacks...; kwargs...])

Holds and coordinates all state of the training. `model` is trained by
optimizing `lossfn` with `optimizer` on `data`.

## Arguments

- `model`
- `data`: Tuple of data iterators in the order `(traindata, valdata, [testdata])`.
    Must be iterable and return tuples of `(xs, ys)`
- `lossfn`: Function with signature `lossfn(model(x), y) -> Number`
- `optimizer`
- `callbacks...`: Any other unnamed arguments are callbacks

## Keyword arguments

- `usedefaultcallbacks = true`: Whether to add some basic callbacks. Included
    are [`Metrics`](#), [`Recorder`](#), [`ProgressPrinter`](#),
    [`StopOnNaNLoss`](#), and [`MetricsPrinter`](#).
- `cbrunner = LinearRunner()`: Callback runner to use.

## Fields

*(Use this as a reference when implementing callbacks)*

- `model`, `optimizer`, and `lossfn` are stored as passed in
- `data` is a `NamedTuple` of `(training = ..., validation = ..., test = ...)`.
    Some values might be `nothing` if you didn't pass in multiple data iterators.
- `params`: an instance of `model`'s parameters of type `Flux.Params`
- `batch`: State of the current batch, including:

    - `batch.xs`: model inputs
    - `batch.ys`: target outputs
    - `batch.ŷs`: model outputs, i.e. `model(xs)`
    - `batch.loss`: batch loss, i.e. `lossfn(ŷs, ys)`
    - `batch.gs`: batch gradients, instance of `Zygote.Grads`

    (!) Note: Depending on the progress of the step, some fields may be `nothing`,
    e.g. the `gs` before the backward pass.
- `cbstate::`[`PropDict`](#): Special state container that callbacks can
    save state to for other callbacks. Its keys depend on what callbacks
    are being used. See the [custom callbacks guide](../docs/callbacks/custom.md)
    for more info.
"""
function Learner(
        model, data, optimizer, lossfn, callbacks::Vararg{<:Callback};
        usedefaultcallbacks = true, cbrunner = LinearRunner()
    )
    callbacks = collect(Callback, callbacks)

    if usedefaultcallbacks
        for cb in defaultcallbacks()
            if !any(typeof(cb) .== typeof.(callbacks))
                push!(callbacks, cb)
            end
        end
    end

    return Learner(
        model, dataiters(data), optimizer, lossfn,
        Flux.params(model),
        BatchState(),
        Callbacks(callbacks, cbrunner),
        PropDict())
end

Base.show(io::IO, learner::Learner) = print(io, "Learner()")

defaultcallbacks()::Vector{AbstractCallback} = [
    ProgressPrinter(),
    MetricsPrinter(),
    StopOnNaNLoss(),
    Recorder(),
    Metrics(),
]

#  Callback handling
handle(event, learner, phase) = handle(learner.callbacks.runner, event, phase, learner)

# Other

getdataiter(phase::AbstractTrainingPhase, learner) = learner.data.training
getdataiter(phase::ValidationPhase, learner) = learner.data.validation

function model!(learner, model)
    learner.model = model
    learner.params = Flux.params(model)
end

numsteps(learner::Protected, phase) = numsteps(getfield(learner, :data), phase)
numsteps(learner, phase) = length(getdataiter(phase, learner))

dataiters(train, val = nothing, test = nothing) = (training = train, validation = val, test = test)
dataiters(t::Tuple) = dataiters(t...)
dataiters(t::NamedTuple) = keys(t) == (:training, :validation, :test) ? t : error("Wrong keys.")

function setcallbacks!(learner, callbacks)
    learner.callbacks = Callbacks(callbacks)
end

function addcallback!(learner, callback)
    learner.callbacks = Callbacks(vcat(learner.callbacks.cbs, callback))
end


text:
"""
    anymap(f, xs)

Like `map`, but force to create `Vector{Any}`.
"""
anymap(f, xs) = Any[f(x) for x in xs]

"""
    default_julia_pass(ir::IRCode, sv::OptimizationState)

The default julia optimization pass.
"""
function default_julia_pass(ir::IRCode, sv::OptimizationState)
    ir = compact!(ir)
    ir = ssa_inlining_pass!(ir, ir.linetable, sv.inlining, sv.src.propagate_inbounds)
    ir = compact!(ir)

    @static if VERSION < v"1.8-"
        ir = Core.Compiler.getfield_elim_pass!(ir)
    else
        ir = Core.Compiler.sroa_pass!(ir)
    end

    ir = adce_pass!(ir)
    ir = type_lift_pass!(ir)
    ir = compact!(ir)
    if JLOptions().debug_level == 2
        verify_ir(ir)
        verify_linetable(ir.linetable)
    end
    return ir
end

"""
    no_pass(ir::IRCode, ::OptimizationState)

No pass.
"""
no_pass(ir::IRCode, ::OptimizationState) = ir

function replace_from_perm(stmt, perm)
    @match stmt begin
        SSAValue(id) => SSAValue(findfirst(isequal(id), perm))
        Expr(head, args...) => Expr(head, map(x->replace_from_perm(x, perm), args)...)
        _ => stmt
    end
end

function permute_stmt(e, perm::Vector{Int})
    @switch e begin
        @case SSAValue(id)
            return SSAValue(findfirst(isequal(id), perm))
        @case ::Expr
            return replace_from_perm(e, perm)
        @case GotoIfNot(cond, dest)
            cond = permute_stmt(cond, perm)
            # NOTE:
            # dest in IRCode refers to #BB
            # don't update it
            # dest = findfirst(isequal(dest), perm)
            return GotoIfNot(cond, dest)
        @case GotoNode(label)
            return GotoNode(findfirst(isequal(label), perm))
        @case ReturnNode(val)
            return ReturnNode(permute_stmt(val, perm))
        @case _
            return e
    end
end

"""
    permute_stmts!(ir::IRCode, perm::Vector{Int})

Permute statements according to `perm`.
"""
function permute_stmts!(ir::IRCode, perm::Vector{Int})
    inst = Any[permute_stmt(ir.stmts.inst[v], perm) for v in perm]
    copyto!(ir.stmts.inst, inst)
    permute!(ir.stmts.flag, perm)
    permute!(ir.stmts.line, perm)
    permute!(ir.stmts.type, perm)
    permute!(ir.stmts.flag, perm)
    return ir
end

@static if VERSION < v"1.8-"
    argtypes(sig::Signature) = sig.atypes
else
    argtypes(sig::Signature) = sig.argtypes
end

function is_allconst(sig::Signature)
    allconst = true
    for atype in argtypes(sig)
        if !isa(atype, Const)
            allconst = false
            break
        end
    end
    return allconst
end

function is_arg_allconst(ir, arg)
    if arg isa Argument
        return false
    elseif arg isa SSAValue
        return is_arg_allconst(ir, ir.stmts[arg.id][:inst])
    elseif !is_inlineable_constant(arg) && !isa(arg, QuoteNode)
        return false
    end
    return true
end

function is_const_call_inlineable(sig::Signature)
    is_allconst(sig) || return false
    f, ft, atypes = sig.f, sig.ft, argtypes(sig)

    if isa(f, IntrinsicFunction) && is_pure_intrinsic_infer(f) && intrinsic_nothrow(f, atypes[2:end])
        return true
    end

    if isa(f, Builtin) && (f === Core.tuple || f === Core.getfield)
        return true
    end
    return false
end

function unwrap_arg(ir, arg)
    if arg isa QuoteNode
        return arg.value
    elseif arg isa SSAValue
        return unwrap_arg(ir, ir.stmts[arg.id][:inst])
    else
        return arg
    end
end

"""
    inline_const!(ir::IRCode)

This performs constant propagation on `IRCode` so after the constant propagation
during abstract interpretation, we can force inline constant values in `IRCode`.
"""
function inline_const!(ir::IRCode)
    for i in 1:length(ir.stmts)
        stmt = ir.stmts[i][:inst]
        @switch stmt begin
            @case GlobalRef(mod, name)
                t = Core.Compiler.abstract_eval_global(mod, name)
                t isa Const || continue
                ir.stmts[i][:inst] = quoted(t.val)
                ir.stmts[i][:type] = t
            @case Expr(:call, f, args...)
                new_stmt = Expr(:call, f, map(eval_global, args)...)
                sig = Core.Compiler.call_sig(ir, new_stmt)
                sig === nothing && continue
                if is_const_call_inlineable(sig)
                    fargs = anymap(x::Const -> x.val, argtypes(sig)[2:end])
                    val = sig.f(fargs...)
                    ir.stmts[i][:inst] = quoted(val)
                    ir.stmts[i][:type] = Const(val)
                else
                    ir.stmts[i][:inst] = new_stmt
                end
            @case Expr(:invoke, mi, f, args...)
                ir.stmts[i][:inst] = Expr(:invoke, mi, f, map(eval_global, args)...)
            @case Expr(:new, t, args...)
                allconst = all(x->is_arg_allconst(ir, x), args)
                allconst && isconcretetype(t) && !ismutabletype(t) || continue
                args = anymap(arg->unwrap_arg(ir, arg), args)
                val = ccall(:jl_new_structv, Any, (Any, Ptr{Cvoid}, UInt32), t, args, length(args))
                ir.stmts[i][:inst] = quoted(val)
                ir.stmts[i][:type] = Const(val)
            @case _
                nothing
        end
    end
    return ir
end

function eval_global(x)
    @switch x begin
        @case GlobalRef(mod, name)
            t = Core.Compiler.abstract_eval_global(mod, name)
            t isa Const || return x
            return quoted(t.val)
        @case _
            return x
    end
end

"""
    const_invoke!(f, ir::IRCode, ref::GlobalRef)

Replace the function invoke `Expr(:invoke, _, ref, args...)` with
`f(args...)` if its arguments `args` are all constant.
"""
function const_invoke!(f, ir::IRCode, ref::GlobalRef)
    for i in 1:length(ir.stmts)
        stmt = ir.stmts[i][:inst]

        @switch stmt begin
            @case Expr(:invoke, _, &ref, args...)
                if all(x->is_arg_allconst(ir, x), args)
                    args = anymap(x->unwrap_arg(ir, x), args)
                    val = f(args...)

                    ir.stmts[i][:inst] = quoted(val)
                    ir.stmts[i][:type] = Const(val)
                end
            @case _
                nothing
        end
    end
    return ir
end

function contains_const_invoke(ir::IRCode, ref::GlobalRef)
    for i in 1:length(ir.stmts)
        stmt = ir.stmts[i][:inst]

        @switch stmt begin
            @case Expr(:invoke, _, &ref, args...)
                all(x->is_arg_allconst(ir, x), args) && return true
            @case _
                nothing
        end
    end
    return false
end


text:
#=

  Mr Greenguest puzzle (fancy dress) in ConstraintSolver.jl

  Problem (and LPL) code in

  http://diuflx71.unifr.ch/lpl/GetModel?name=/demo/demo2

  """
  (** Mr. Greenfan wants to give a dress party where the male guests
   * must wear green dresses. The following rules are given:
   * 1 If someone wears a green tie he has to wear a green shirt.
   * 2 A guest may only wear green socks and a green shirt
   *   if he wears a green tie or a green hat.
   * 3 A guest wearing a green shirt or a green hat or who does
   *   not wear green socks must wear a green tie.
   * 4 A guest who is not dressed according to rules 1-3 must
   *   pay a $11 entrance fee.
   * Mr Greenguest wants to participate but owns only a green shirt
   * (otherwise he would have to pay one for $9). He could buy
   * a green tie for $10, a green hat (used) for $2 and green socks
   * for $12.
   * What is the cheapest solution for Mr Greenguest to participate?
   *)
  """

  The solution:
    t = 1
    h = 0
    r = 1
    s = 0
    n = 0
    cost = 10

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function fancy(print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    @variable(model, x[1:5], Bin)
    T,H,R,S,N = x
    @variable(model, 0 <= cost <= 10000, Int)

    #  This is a straight translation from the LPL
    #
    # ( ((T ==1) => (R ==1)) #\/ N ==1)  % Picat constraint
    @constraint(model, T := { R ==1 || N ==1 } )

    # ( ((S ==1 \/ R==1) => (T==1 \/ H==1)) \/ N==1)
    b1 = @variable(model, binary=true)
    @constraint(model, b1 := {S ==1 || R==1})
    b2 = @variable(model, binary=true)
    @constraint(model, b2 := {N==1})
    @constraint(model, b1 => {T ==1 || H ==1 } )
    @constraint(model, b1 + b2 >= 1)

    # ( ((R ==1 \/ H==1 \/ (S!=1)) => T==1) \/ N==1)
    b3 = @variable(model, binary=true)
    @constraint(model, b3 := {R ==1 || H==1 || S!=1})
    @constraint(model, b1 + b3 >= 1)

    @constraint(model, cost == 10*T + 2*H + 12*S + 11*N)

    @objective(model,Min,cost)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                cost_val = convert.(Integer,JuMP.value.(cost; result=sol))
                println("x:$x_val  cost:$cost_val")
                println("t: $(x_val[1])")
                println("h: $(x_val[2])")
                println("r: $(x_val[3])")
                println("s: $(x_val[4])")
                println("n: $(x_val[5])")
                println("cost: $cost_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

@time fancy()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text: println(x = 42)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Unzip: Unzip
using Documenter: doctest

doctest(Unzip)


text:
# Quadrotor in Maze
T = Float64

# model
model = Dynamics.quadrotor
# model = Dynamics.quadrotor_euler
model_d = rk3(model)
n = model.n; m = model.m
q0 = [1.;0.;0.;0.] # unit quaternion

x0 = zeros(T,n)
x0[1:3] = [0.; 0.; 10.]
x0[4:7] = q0

xf = zero(x0)
xf[1:3] = [0.;60.; 10.]
xf[4:7] = q0;

# cost
Q = (1.0e-3)*Diagonal(I,n)
Q[4:7,4:7] = (1.0e-2)*Diagonal(I,4)
R = (1.0e-4)*Diagonal(I,m)
Qf = 1000.0*Diagonal(I,n)

u_min = 0.
u_max = 50.
x_max = Inf*ones(model.n)
x_min = -Inf*ones(model.n)

x_max[1:3] = [25.0; Inf; 20]
x_min[1:3] = [-25.0; -Inf; 0.]
bnd3 = BoundConstraint(n,m,u_min=u_min)

xf_no_quat_U = copy(xf)
xf_no_quat_L = copy(xf)
xf_no_quat_U[4:7] .= Inf
xf_no_quat_L[4:7] .= -Inf
xf_no_quat_U[8:10] .= 0.
xf_no_quat_L[8:10] .= 0.
bnd_xf = BoundConstraint(n,m,x_min=xf_no_quat_L,x_max=xf_no_quat_U)

N = 101 # number of knot points
tf = 5.0
dt = tf/(N-1) # total time

U_hover = [0.5*9.81/4.0*ones(m) for k = 1:N-1] # initial hovering control trajectory
obj = LQRObjective(Q, R, Qf, xf, N) # objective with same stagewise costs

quadrotor = Problem(model_d, obj, x0=x0, xf=xf, N=N, dt=dt)
initial_controls!(quadrotor,U_hover); # initialize problem with controls

quadrotor.constraints[1] += bnd3
for k = 2:N-1
    quadrotor.constraints[k] += bnd3
end
quadrotor.constraints[N] += bnd_xf

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to Date and Time conversion.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] http://aa.usno.navy.mil/data/docs/JulianDate.php
#   [2] https://www.ietf.org/timezones/data/leap-seconds.list
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/time/time.jl
# ========================

# Function get_Δat
# ----------------

@testset "Function get_Δat"  begin
    # Leap seconds values obtained from [2].
    leap_secs = [2272060800	10	1 "Jan" 1972;
                 2303683200	12	1 "Jan" 1973;
                 2287785600	11	1 "Jul" 1972;
                 2335219200	13	1 "Jan" 1974;
                 2366755200	14	1 "Jan" 1975;
                 2398291200	15	1 "Jan" 1976;
                 2429913600	16	1 "Jan" 1977;
                 2461449600	17	1 "Jan" 1978;
                 2492985600	18	1 "Jan" 1979;
                 2524521600	19	1 "Jan" 1980;
                 2571782400	20	1 "Jul" 1981;
                 2603318400	21	1 "Jul" 1982;
                 2634854400	22	1 "Jul" 1983;
                 2698012800	23	1 "Jul" 1985;
                 2776982400	24	1 "Jan" 1988;
                 2840140800	25	1 "Jan" 1990;
                 2871676800	26	1 "Jan" 1991;
                 2918937600	27	1 "Jul" 1992;
                 2950473600	28	1 "Jul" 1993;
                 2982009600	29	1 "Jul" 1994;
                 3029443200	30	1 "Jan" 1996;
                 3076704000	31	1 "Jul" 1997;
                 3124137600	32	1 "Jan" 1999;
                 3345062400	33	1 "Jan" 2006;
                 3439756800	34	1 "Jan" 2009;
                 3550089600	35	1 "Jul" 2012;
                 3644697600	36	1 "Jul" 2015;
                 3692217600	37	1 "Jan" 2017;
                ]

    for i = 1:size(leap_secs,1)
        ΔAT   = leap_secs[i,2]
        day   = leap_secs[i,3]
        month = begin
            if leap_secs[i,4] == "Jan"
                1
            elseif leap_secs[i,4] == "Jul"
                7
            else
                error("Invalid array `leap_secs`.")
            end
        end
        year  = leap_secs[i,5]

        @test get_Δat(date_to_jd(year,month,day,0,0,0))         == ΔAT

        # One second before the entry in `leap_secs`, we must have `ΔAT-1` leap
        # seconds. However, this is not true for the very first line.
        (i == 1) && continue
        @test get_Δat(date_to_jd(year,month,day,0,0,0)-1/86400) == ΔAT-1
    end
end

# Functions jd_ut1_to_utc and jd_utc_to_ut1
# -----------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-7: Calculating Dynamical Time [1].
#
# The following date:
#
#   Mountain Standard Time Zone: May 14, 2004, 10:43
#
# is converted into the following one using ΔUT1 = -0.463326 s:
#
#   UT1: May 14, 2004, 16:42:59.5367
#
# Scenario 02
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction [1].
#
# The following UTC date:
#
#   UTC: April 6, 2004, 07:51:28.386009
#
# is converted into the following one using ΔUT1 = -0.4399619 s:
#
#   UT1: April 6, 2004, 07:51:27.946047
#
################################################################################

@testset "Functions jd_ut1_to_utc and jd_utc_to_ut1" begin

    ## Scenario 01
    ## ===========

    ΔUT1 = -0.463326

    ## jd_utc_to_ut1
    ## -----------

    # At the mentioned date, Mountain Standard Time is 6h behind UTC.
    JD_UTC = date_to_jd(2004, 5, 14, 10+6, 43, 0)
    JD_UT1 = jd_utc_to_ut1(JD_UTC, ΔUT1)

    (year, month, day, hour, minute, second) = jd_to_date(JD_UT1)

    @test year   == 2004
    @test month  == 5
    @test day    == 14
    @test hour   == 16
    @test minute == 42
    @test second  ≈ 59.5367 atol=1e-4

    ## jd_ut1_to_utc
    ## -------------

    JD_UT1 = date_to_jd(2004, 5, 14, 16, 42, 59.5367)
    JD_UTC = jd_ut1_to_utc(JD_UT1, ΔUT1)

    (year, month, day, hour, minute, second) = jd_to_date(JD_UTC)

    @test year   == 2004
    @test month  == 5
    @test day    == 14
    @test hour   == 10+6
    @test minute == 43
    @test second  ≈ 0.0000 atol=1e-4

    ## Scenario 02
    ## ===========

    ΔUT1 = -0.4399619

    ## jd_utc_to_ut1
    ## -------------

    JD_UTC = date_to_jd(2004, 4, 6, 07, 51, 28.386009)
    JD_UT1 = jd_utc_to_ut1(JD_UTC, ΔUT1)

    (year, month, day, hour, minute, second) = jd_to_date(JD_UT1)

    @test year   == 2004
    @test month  == 4
    @test day    == 6
    @test hour   == 7
    @test minute == 51
    @test second  ≈ 27.946047 atol=1e-4

    ## jd_ut1_to_utc
    ## -------------

    JD_UT1 = date_to_jd(2004, 4, 6, 07, 51, 27.946047)
    JD_UTC = jd_ut1_to_utc(JD_UT1, ΔUT1)

    (year, month, day, hour, minute, second) = jd_to_date(JD_UTC)

    @test year   == 2004
    @test month  == 4
    @test day    == 6
    @test hour   == 7
    @test minute == 51
    @test second  ≈ 28.386009 atol=1e-4
end

# Functions jd_tt_to_utc and jd_utc_to_tt
# ---------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-7: Calculating Dynamical Time [1].
#
# The following date:
#
#   Mountain Standard Time Zone: May 14, 2004, 10:43
#
# is converted into:
#
#   TT: May 14, 2004, 16:44:04.1840
#
# Scenario 02
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction [1].
#
# The following UTC date:
#
#   UTC: April 6, 2004, 07:51:28.386009
#
# is converted into:
#
#   JD [TT]: 2453101.828154745
#
################################################################################

@testset "Functions jd_tt_to_utc and jd_utc_to_tt" begin

    ## Scenario 01
    ## ===========

    ## jd_utc_to_tt
    ## ------------

    # At the mentioned date, Mountain Standard Time is 6h behind UTC.
    JD_UTC = date_to_jd(2004, 5, 14, 10+6, 43, 0)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    (year, month, day, hour, minute, second) = jd_to_date(JD_TT)

    @test year   == 2004
    @test month  == 5
    @test day    == 14
    @test hour   == 16
    @test minute == 44
    @test second  ≈ 04.1840 atol=1e-4

    ## jd_tt_to_utc
    ## ------------

    JD_TT  = date_to_jd(2004, 5, 14, 16, 44, 04.1840)
    JD_UTC = jd_tt_to_utc(JD_TT)

    (year, month, day, hour, minute, second) = jd_to_date(JD_UTC)

    @test year   == 2004
    @test month  == 5
    @test day    == 14
    @test hour   == 10+6

    # TODO: Fix this rounding problem.
    if minute == 42
        @test minute == 42
        @test second  ≈ 60.0000 atol=1e-4
    else
        @test minute == 43
        @test second  ≈ 0.0000 atol=1e-4
    end

    ## Scenario 02
    ## ===========

    ## jd_utc_to_tt
    ## ------------

    # At the mentioned date, Mountain Standard Time is 6h behind UTC.
    JD_UTC = date_to_jd(2004, 4, 6, 07, 51, 28.386009)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    @test JD_TT ≈ 2453101.828154745 atol=1e-9

    ## jd_tt_to_utc
    ## ------------

    JD_TT  = 2453101.828154745
    JD_UTC = jd_tt_to_utc(JD_TT)

    (year, month, day, hour, minute, second) = jd_to_date(JD_UTC)

    @test year   == 2004
    @test month  == 4
    @test day    == 6
    @test hour   == 07
    @test minute == 51
    @test second  ≈ 28.386009 atol=1e-4

end

# File: ./src/time/julian_day.jl
# ==============================

# Function jd_to_gmst
# -------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Data obtained from [1]:
#
#   ╔═════════════════════╦════════════════╗
#   ║    Gregorian Day    ║   Julian Day   ║
#   ╠═════════════════════╬════════════════╣
#   ║ 1986-06-19 21:35:22 ║ 2446601.399560 ║
#   ║ 1987-05-19 04:00:00 ║ 2446934.666667 ║
#   ║ 2018-04-16 20:19:37 ║ 2458225.346956 ║
#   ║ 1822-09-07 12:00:00 ║ 2386781.000000 ║
#   ║ 1900-01-01 00:00:00 ║ 2415020.500000 ║
#   ║ 2000-01-01 00:00:00 ║ 2451544.500000 ║
#   ║ 2013-10-19 22:00:00 ║ 2456585.416667 ║
#   ╚═════════════════════╩════════════════╝
#
# Data obtained from SatelliteToolbox:
#
#   ╔═════════════════════╦════════════════╗
#   ║    Gregorian Day    ║   Julian Day   ║
#   ╠═════════════════════╬════════════════╣
#   ║ 1986-06-19 21:35:22 ║ 2446601.399560 ║
#   ║ 1987-05-19 04:00:00 ║ 2446934.666667 ║
#   ║ 2018-04-16 20:19:37 ║ 2458225.346956 ║
#   ║ 1822-09-07 12:00:00 ║ 2386781.000000 ║
#   ║ 1900-01-01 00:00:00 ║ 2415020.500000 ║
#   ║ 2000-01-01 00:00:00 ║ 2451544.500000 ║
#   ║ 2013-10-19 22:00:00 ║ 2456585.416667 ║
#   ╚═════════════════════╩════════════════╝
#
################################################################################

@testset "Function jd_to_gmst" begin
    @test date_to_jd(1986, 06, 19, 21, 35, 22) ≈ 2446601.399560 atol=1e-6
    @test date_to_jd(1987, 05, 19, 04, 00, 00) ≈ 2446934.666667 atol=1e-6
    @test date_to_jd(2018, 04, 16, 20, 19, 37) ≈ 2458225.346956 atol=1e-6
    @test date_to_jd(1822, 09, 07, 12, 00, 00) ≈ 2386781.000000 atol=1e-6
    @test date_to_jd(1900, 01, 01, 00, 00, 00) ≈ 2415020.500000 atol=1e-6
    @test date_to_jd(2000, 01, 01, 00, 00, 00) ≈ 2451544.500000 atol=1e-6
    @test date_to_jd(2013, 10, 19, 22, 00, 00) ≈ 2456585.416667 atol=1e-6
end

# Function jd_to_date
# -------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Data obtained from [1]:
#
#   ╔═════════════════════╦════════════════╗
#   ║    Gregorian Day    ║   Julian Day   ║
#   ╠═════════════════════╬════════════════╣
#   ║ 1986-06-19 21:35:22 ║ 2446601.399560 ║
#   ║ 1987-05-19 04:00:00 ║ 2446934.666667 ║
#   ║ 2018-04-16 20:19:37 ║ 2458225.346956 ║
#   ║ 1822-09-07 12:00:00 ║ 2386781.000000 ║
#   ║ 1900-01-01 00:00:00 ║ 2415020.500000 ║
#   ║ 2000-01-01 00:00:00 ║ 2451544.500000 ║
#   ║ 2013-10-19 22:00:00 ║ 2456585.416667 ║
#   ╚═════════════════════╩════════════════╝
#
# Data obtained from SatelliteToolbox:
#
#   ╔═════════════════════╦════════════════╗
#   ║    Gregorian Day    ║   Julian Day   ║
#   ╠═════════════════════╬════════════════╣
#   ║ 1986-06-19 21:35:22 ║ 2446601.399560 ║
#   ║ 1987-05-19 04:00:00 ║ 2446934.666667 ║
#   ║ 2018-04-16 20:19:37 ║ 2458225.346956 ║
#   ║ 1822-09-07 12:00:00 ║ 2386781.000000 ║
#   ║ 1900-01-01 00:00:00 ║ 2415020.500000 ║
#   ║ 2000-01-01 00:00:00 ║ 2451544.500000 ║
#   ║ 2013-10-19 22:00:00 ║ 2456585.416667 ║
#   ╚═════════════════════╩════════════════╝
#
################################################################################

@testset "Function jd_to_date" begin
    @test jd_to_date(Int,2446601.399560) == (1986,06,19,21,35,22)
    @test jd_to_date(Int,2446934.666667) == (1987,05,19,04,00,00)
    @test jd_to_date(Int,2458225.346956) == (2018,04,16,20,19,37)
    @test jd_to_date(Int,2386781.000000) == (1822,09,07,12,00,00)
    @test jd_to_date(Int,2415020.500000) == (1900,01,01,00,00,00)
    @test jd_to_date(Int,2451544.500000) == (2000,01,01,00,00,00)
    @test jd_to_date(Int,2456585.416667) == (2013,10,19,22,00,00)
end

# Function date_to_jd
# -------------------

# Issue #19
# ==============================================================================

@testset "Issue #19" begin
    year_i = rand(1900:2100)

    year, month, day, hour, min, sec = jd_to_date(date_to_jd(year_i,12, 31, 23, 59, 59))

    @test year  == year_i
    @test month == 12
    @test day   == 31
    @test hour  == 23
    @test min   == 59
    @test sec    ≈ 59 rtol = 1e-5

    year, month, day, hour, min, sec = jd_to_date(date_to_jd(year_i, 1, 1, 0, 0, 0))

    @test year  == year_i
    @test month == 1
    @test day   == 1
    @test hour  == 0
    @test min   == 0
    @test sec    ≈ 0 rtol = 1e-5

    date_time_1 = jd_to_date(DateTime, date_to_jd(year_i, 12, 31, 23, 59, 59))
    date_time_2 = DateTime(year_i, 12, 31, 23, 59, 59)

    @test date_time_1 === date_time_2
end

# Issue #30
# ==============================================================================

@testset "Issue #30" begin
    year, month, day, h, m, s = jd_to_date(Int, 2.45754095833333e6)

    @test year  == 2016
    @test month == 6
    @test day   == 1
    @test h     == 11
    @test m     == 0
    @test s     == 0
end

# Function date_to_jd(::DateTime)
# ==============================================================================

@testset "date_to_jd with DateTime" begin
    JDdatetime = date_to_jd(DateTime(2020, 8, 14, 12, 4, 1))
    JDnums     = date_to_jd(2020, 8, 14, 12, 4, 1)
    @test JDdatetime == JDnums

    JDdatetime_ms = date_to_jd(DateTime(2020, 8, 14, 12, 4, 1, 500))
    JDnums_ms     = date_to_jd(2020, 8, 14, 12, 4, 1.5)
    JDdatetime_back_ms = jd_to_date(DateTime,JDnums_ms)
    @test JDdatetime_ms ≈ JDnums_ms atol=6e-9
    @test JDdatetime_back_ms == DateTime(2020, 8, 14, 12, 4, 1, 500)
end


text:
export run_sim, run_sim_total_recombine

function selection_flip(par::PopRates; var_sites = [], prob = 0.1)
    # Generates a new parameter struct with variable sites flipped at random with probability = prob
    n = length(var_sites) # length of vectors with selection
    k = rand(Binomial(n,prob))
    inds = var_sites[StatsBase.sample(1:length(var_sites),k, replace=false)]
    β1 = β1_vector(par.βf)
    β1[inds] .*= -1
    return PopRates(par, β1 = β1)
end

function run_sim(pop::PopState, par::PopRates, timepoints;
	var_sites = [], # vector of indices that are flipping
	flip_prob = 0.0) # probability that they flip at every time point.
    df = initialize_stats() # empty data frame for holding the stats
    for t in timepoints
        run_until!(pop,par,t) # evolve the poulation state until pop.time > t.
        record_stats!(df, pop, par)  # record the statistics in the data_frame
        par = selection_flip(par; var_sites = var_sites, prob = flip_prob) # possibly stochastically change the environment.
    end
    return df # return the df, now heavy with juicy statistics for you to plot or analyze
end

# default definition when pop is provided
function run_sim(par::PopRates, timepoints; kwdargs...) # population loci set to wright equilibrium,
	pop = initialize_pop(par)
	run_sim(pop, par, timepoints; kwdargs...)
end

"""
bottleneck removes all individuals for which survives = false
"""
function bottleneck!(pop::PopState, survives::Function)
    pop.individuals = pop.individuals[survives.(pop.individuals)]
end

function shuffle_col!(data)
    vdata = Align(data,1)
    for i = 1:size(vdata, 1)
        shuffle!(@view vdata[i,:])
    end
    return
end

function total_recombine(pop::PopState)
    shuffle_col!(pop.individuals)
end

function run_sim_total_recombine(par::PopRates, timepoints; kwdargs...) # population loci set to wright equilibrium,
	pop = initialize_pop(par)
	run_sim_total_recombine(pop, par, timepoints; kwdargs...)
end

function run_sim_total_recombine(pop::PopState, par::PopRates, timepoints;
	var_sites = [], # vector of indices that are flipping
	flip_prob = 0.0) # probability that they flip at every time point.
    df = initialize_stats() # empty data frame for holding the stats
    for t in timepoints
        run_until!(pop,par,t) # evolve the poulation state until pop.time > t.
        record_stats!(df, pop, par)  # record the statistics in the data_frame
        total_recombine(pop) # get near limit of infinite recombination.
        par = selection_flip(par; var_sites = var_sites, prob = flip_prob) # possibly stochastically change the environment.
    end
    return df # return the df, now heavy with juicy statistics for you to plot or analyze
end

max_stars_repo_name: StatPhysBio/Tomoko

max_issues_repo_name: StatPhysBio/Tomoko

max_forks_repo_name: StatPhysBio/Tomoko


text:
using Test

using SumOfSquares
const CEG = SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph

function fill_in_3_nodes_test(G, x, y, z)
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 0
    @test CEG.fill_in(G, z) == 0
    CEG.add_edge!(G, x, y)
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 0
    @test CEG.fill_in(G, z) == 0
    CEG.add_edge!(G, y, z)
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 1
    @test CEG.fill_in(G, z) == 0
    @testset "Disable y" begin
        H = copy(G)
        CEG.disable_node!(H, y)
        @test CEG.fill_in(H, x) == 0
        @test CEG.fill_in(H, z) == 0
        CEG.add_edge!(H, z, x)
        @test CEG.fill_in(H, x) == 0
        @test CEG.fill_in(H, z) == 0
    end
    CEG.add_edge!(G, z, x)
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 0
    @test CEG.fill_in(G, z) == 0
end

function fill_in_4_nodes_test(G, w, x, y, z)
    @test CEG.fill_in(G, w) == 0
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 0
    @test CEG.fill_in(G, z) == 0
    CEG.add_edge!(G, w, x)
    @test CEG.fill_in(G, w) == 0
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 0
    @test CEG.fill_in(G, z) == 0
    CEG.add_edge!(G, y, z)
    @test CEG.fill_in(G, w) == 0
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 0
    @test CEG.fill_in(G, z) == 0
    CEG.add_edge!(G, w, y)
    @test CEG.fill_in(G, w) == 1
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 1
    @test CEG.fill_in(G, z) == 0
    @testset "wz first" begin
        H = copy(G)
        CEG.add_edge!(H, w, z)
        @test CEG.fill_in(H, w) == 2
        @test CEG.fill_in(H, x) == 0
        @test CEG.fill_in(H, y) == 0
        @test CEG.fill_in(H, z) == 0
        @testset "Disable w" begin
            I = copy(H)
            CEG.add_edge!(I, x, z)
            @test CEG.fill_in(I, x) == 0
            @test CEG.fill_in(I, y) == 0
            @test CEG.fill_in(I, z) == 1
            CEG.add_edge!(I, x, y)
            @test CEG.fill_in(I, x) == 0
            @test CEG.fill_in(I, y) == 0
            @test CEG.fill_in(I, z) == 0
        end
        CEG.add_edge!(H, x, z)
        @test CEG.fill_in(H, w) == 1
        @test CEG.fill_in(H, x) == 0
        @test CEG.fill_in(H, y) == 0
        @test CEG.fill_in(H, z) == 1
        CEG.add_edge!(H, x, y)
        @test CEG.fill_in(H, w) == 0
        @test CEG.fill_in(H, x) == 0
        @test CEG.fill_in(H, y) == 0
        @test CEG.fill_in(H, z) == 0
    end
    CEG.add_edge!(G, x, z)
    @test CEG.fill_in(G, w) == 1
    @test CEG.fill_in(G, x) == 1
    @test CEG.fill_in(G, y) == 1
    @test CEG.fill_in(G, z) == 1
    @testset "xy then wz" begin
        H = copy(G)
        CEG.add_edge!(H, x, y)
        @test CEG.fill_in(H, w) == 0
        @test CEG.fill_in(H, x) == 1
        @test CEG.fill_in(H, y) == 1
        @test CEG.fill_in(H, z) == 0
        CEG.add_edge!(H, w, z)
        @test CEG.fill_in(H, w) == 0
        @test CEG.fill_in(H, x) == 0
        @test CEG.fill_in(H, y) == 0
        @test CEG.fill_in(H, z) == 0
    end
    CEG.add_edge!(G, w, z)
    @test CEG.fill_in(G, w) == 1
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 0
    @test CEG.fill_in(G, z) == 1
    @testset "Disable w" begin
        H = copy(G)
        @test CEG.fill_in(H, x) == 0
        @test CEG.fill_in(H, y) == 0
        @test CEG.fill_in(H, z) == 1
        CEG.add_edge!(H, x, y)
        @test CEG.fill_in(H, x) == 0
        @test CEG.fill_in(H, y) == 0
        @test CEG.fill_in(H, z) == 0
    end
    CEG.add_edge!(G, x, y)
    @test CEG.fill_in(G, w) == 0
    @test CEG.fill_in(G, x) == 0
    @test CEG.fill_in(G, y) == 0
    @test CEG.fill_in(G, z) == 0
end

@testset "Chordal Extensions" begin
    @testset "CEG.LabelledGraph" begin
        @test CEG.LabelledGraph{Int}() isa CEG.LabelledGraph{Int}
        @test_throws ErrorException CEG.LabelledGraph()
        G = CEG.LabelledGraph{Int}()
        @test CEG.add_node!(G, 1) == 1
        @test CEG.add_edge!(G, (1, 2)) == (1, 2)
        @test CEG.add_edge!.(G, [(1, 2), (2, 3)]) isa Vector{Tuple{Int, Int}}
        @test CEG.add_clique!(G, [1, 2, 3]) == nothing
    end

    @testset "show" begin
        G = CEG.LabelledGraph{Symbol}()
        @test sprint(show, G) == "Symbol-LabelledGraph with\nNodes:\nEdges:\n"
        CEG.add_edge!(G, :x, :y)
        @test sprint(show, G) isa String
    end

    @testset "sub_sets" begin
        G = CEG.LabelledGraph{Symbol}()
        CEG.add_clique!(G, [:x, :y, :z])
        @test CEG.neighbors(G.graph, 1) == Set([2, 3])
    end

    @testset "fill_in with 3 nodes" begin
        G = CEG.Graph()
        x = CEG.add_node!(G)
        y = CEG.add_node!(G)
        z = CEG.add_node!(G)
        fill_in_3_nodes_test(copy(G), x, y, z)
        fill_in_3_nodes_test(CEG.FillInCache(copy(G)), x, y, z)
    end

    @testset "fill_in with 4 nodes" begin
        G = CEG.Graph()
        w = CEG.add_node!(G)
        x = CEG.add_node!(G)
        y = CEG.add_node!(G)
        z = CEG.add_node!(G)
        fill_in_4_nodes_test(copy(G), w, x, y, z)
        fill_in_4_nodes_test(CEG.FillInCache(copy(G)), w, x, y, z)
    end

    @testset "cliques" begin
        G = CEG.LabelledGraph{Symbol}()
        CEG.add_edge!(G, :x, :y)
        CEG.add_edge!(G, :y, :z)
        @test CEG.fill_in.(G.graph, [1, 2, 3]) == [0, 1, 0]
        @test !CEG.is_clique(G.graph, [1, 2, 3])
        @test CEG.is_clique(G.graph, [1, 2])
    end

    @testset "chordal" begin
        G = CEG.LabelledGraph{Symbol}()
        CEG.add_edge!.(G, [(:x, :y), (:y, :z)])
        H, cliques = CEG.chordal_extension(G, CEG.GreedyFillIn())
        @test H.int2n == G.int2n
        @test H.graph.neighbors == G.graph.neighbors
        @test cliques == [[:y, :z], [:x, :y]]

        G = CEG.LabelledGraph{Int}()
        CEG.add_edge!.(G, [(1, 2), (1, 3), (3, 4), (2, 4)])
        H, cliques = CEG.chordal_extension(G, CEG.GreedyFillIn())
        @test length.(cliques) == [3, 3]
    end

    @testset "Add edges between two `chordal_extension` calls" begin
        G = CEG.LabelledGraph{Int}()
        CEG.add_node!.(G, [1, 2, 3, 4, 5])
        CEG.add_edge!(G, 1, 2)
        CEG.add_edge!(G, 2, 3)
        CEG.add_edge!(G, 3, 4)
        CEG.add_edge!(G, 1, 4)

        H, cliques = CEG.chordal_extension(G, CEG.GreedyFillIn())
        @test cliques == [[2, 3, 4], [1, 2, 4], [5]]
        @test isempty(G.graph.disabled)
        @test isempty(H.graph.disabled)

        CEG.add_edge!(H, 1, 5)
        I, cliques = CEG.chordal_extension(H, CEG.GreedyFillIn())
        @test cliques == [[1, 2, 4], [2, 3, 4], [1, 5]]
        @test isempty(G.graph.disabled)
        @test isempty(H.graph.disabled)
        @test isempty(I.graph.disabled)
    end

end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
# It can be seen that the number, 125874, and its double, 251748, contain
# exactly the same digits, but in a different order.
#
# Find the smallest positive integer, x, such that 2x, 3x, 4x, 5x, and 6x,
# contain the same digits.

using ProjectEulerSolutions

# Check sets of digits to see if they are equivalent, but the trick is to cut
# down the search range from 100 to 166, 1000 to 1666, etc.  This becomes much
# simpler with the use of sort and the built-in digits function.
function p052solution_builtins(n::Integer=3)::Integer

    for factor_power in 2:n
        # Number must start with a 1
        x_start = 10^factor_power
        # Upper bound must be the same number of digits
        x_end = div(10^(factor_power + 1), 6)
        for x in x_start:x_end
            digit_set = sort(digits(x * 2))
            if digit_set == sort(digits(x * 3)) &&
               digit_set == sort(digits(x * 4)) &&
               digit_set == sort(digits(x * 5)) &&
               digit_set == sort(digits(x * 6))
                return x
            end
        end
    end
    return -1
end

# Explicitly handles the digit checking rather than using built-ins.  Faster,
# but more complicated.
function p052solution_explicit(n::Integer=3)::Integer
    for n_digits in 3:n
        # Number must start with a 1
        x_start = 10^(n_digits - 1)
        # Upper bound must be the same number of digits
        x_end = div(10^n_digits, 6)
        digs = zeros(Int8, 10)
        for x in x_start:x_end
            fill!(digs, 0)
            broke_early = false
            for m in 1:5
                num = x * (m + 1)
                for j in 1:n_digits
                    num2 = (num % 10) + 1
                    digs[num2] += 1
                    # Break early if one of the multipliers doesn't match
                    if digs[num2] != m
                        broke_early = true
                        break
                    end
                    num = div(num, 10)
                end
                if broke_early  # Inner loop broke early
                    break
                end
            end
            if ! broke_early
                return x
            end
        end
    end
    return -1
end

p052 = Problems.Problem(Dict("Built-ins" => p052solution_builtins,
                             "Explicit" => p052solution_explicit))

Problems.benchmark(p052, 7)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


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

using Markdown
using InteractiveUtils

# ╔═╡ 8d2a2140-4d6b-11eb-1ed0-5dfeba2603af
begin
	using Pkg
    Pkg.activate(pwd())
end

# ╔═╡ b9452d10-4d6b-11eb-1074-83294cc14eec
begin
	using CSV
	using GraphIO
	using PlutoUI
	using GraphPlot
	using DataFrames
	using LightGraphs
	using SimpleWeightedGraphs
end

# ╔═╡ a4778db2-4d6b-11eb-0e97-c368743a47ea
md"""
# Data Mining
## Semester Project - Clustering by Density
### Aadam - CS1945
"""

# ╔═╡ b046df70-4e5a-11eb-350f-d1abb711d409
md"""
## Utility Functions
"""

# ╔═╡ 822e6af0-4d73-11eb-32e2-216f0814a3c0
# Convert string keys into integer mappings
function str_int(array)
	d = Dict{String, Int}()
	for (i, s) in enumerate(array)
		push!(d, s => i)
	end
	d
end

# ╔═╡ f0bd31ee-4d77-11eb-2834-85e4981d8680
# Return key of dictionary containing value 'i'
function get_key(i, d)
	for (k,v) in d
		if v==i
			return k
		end
	end
end

# ╔═╡ f60d9610-4d75-11eb-002e-6f23a07f66ad
# Return a cleaned DataFrame from a CSV file
function get_df(path)
	df = CSV.read(path, DataFrame, datarow=1)
	d = str_int(union(df.Column1, df.Column2))
	df.sources = [d[i] for i in df.Column1]
	df.destinations = [d[i] for i in df.Column2]
	rename!(df, :Column3 => :weights)
	return df, d
end

# ╔═╡ 097a9b52-4e37-11eb-2afc-0102d9e2c69d
# Neighbors of a cluster are the nodes connected to any node of cluster but
# not part of the cluster.
function get_neighbors(g, c, excluded)
	neighbors_set = Set()

	for n in c
		union!(neighbors_set, neighbors(g, n))
	end

	return collect(setdiff(setdiff(neighbors_set, c), excluded))
end

# ╔═╡ 30c41780-4e5b-11eb-0508-3b2e824a0607
md"""
The priority is determined based on two measures
1. the sum of the weights of the edges between a neighbor and the cluster,
2. the number of edges between a neighbor and the cluster.
"""

# ╔═╡ d4c26d30-4e3f-11eb-07d2-2d722ea4c878
function count_edges(g, cluster, n)
	count = 0
	for c in cluster
		count += has_edge(g, c, n)
	end
	return count
end

# ╔═╡ 95d86970-4e45-11eb-0929-692e7fab4639
# Remove nodes from the graph that are already in the cluster
function rem_nodes(g, cluster)
	for c in cluster
		rem_vertex!(g, c)
	end
	g
end

# ╔═╡ 7b666270-4e5b-11eb-1f65-b5eafbc71680
md"""
**Density** is given as:

$d = \frac{2 \times \left|E\right|}{\left|N\right| \times (\left|N\right| - 1)}$
"""

# ╔═╡ 3df86fd0-4e3f-11eb-3a9c-1b73eabf5568
function density(N, E)
	return (2 * E)/(N * (N - 1))
end

# ╔═╡ c063ab80-4e5b-11eb-283f-6749471fae2d
md"""
To determine whether node is part of periphery we use **cluster property** `cp` of node. It is given as:

$cp = \frac{\left|E_c\right|}{d \times \left|N_c\right|}$
"""

# ╔═╡ ab04f4e0-4e49-11eb-19be-49390f4285e3
function periphery(N, E, d)
	return E/(d * N)
end

# ╔═╡ eba1e8b0-4e57-11eb-1cee-252cd18c350c
md"""
## Algorithm
- Start with a single node as a cluster
- Grow cluster by adding nodes from neighbors one by one
- Continue expanding cluster as long as following 2 conditions are satisfied
  - Density of cluster > = d’
  - Node is in periphery of cluster
- Remove Cluster from graph
- Apply same procedure to rest of nodes to find other clusters in graph.

---

We will follow the algorithm as shown below:
$(LocalResource("algo.png"))
"""

# ╔═╡ ba80451e-4d6b-11eb-0949-6de777cd3780
filename = joinpath("data", "PPI-I.txt")

# ╔═╡ 5655df00-4d76-11eb-01d6-9d21ce3aa823
df, d = get_df(filename)

# ╔═╡ baf68690-4d6b-11eb-3151-b5e91d478309
g = SimpleWeightedGraph(
	Vector(df.sources),
	Vector(df.destinations),
	Vector(df.weights)
)

# ╔═╡ ebca8580-4e3a-11eb-2a8b-711f8bc91692
function get_priorities(graph, cluster, neighbors, edge_weights)
	m1 = Vector()
	m2 = Vector()

	for n in neighbors
		weight_sum = 0
		edge_count = 0

		for c in cluster
			weight_sum += edge_weights[n, c]
			edge_count += has_edge(g, n, c)
		end
		push!(m1, weight_sum)
		push!(m2, edge_count)
	end
	return m1, m2
end

# ╔═╡ 843b4660-4d7d-11eb-387c-6ddb73c11e73
function density_clustering(G, dʼ, cpʼ)
	g = deepcopy(G)
	clusters = Vector()

	# Δ(g) returns highest node degree (LightGraphs)
	while (Δ(g) != 0)

		edge_weights = adjacency_matrix(g) * adjacency_matrix(g)
		for i in 1:size(edge_weights)[1]
			edge_weights[i, i] = 0
		end
		node_weights = degree_matrix(g)

		highest_nw = findall(node_weights .== maximum(node_weights))[1]
		starting_node = highest_nw != 0 ? highest_nw[1] : Δ(g)

		cluster = Vector()
		push!(cluster, starting_node)
		cluster_edges = 0
		excluded = Set()

		while true
			cluster_neighbors = get_neighbors(g, cluster, excluded)
			m1, m2 = get_priorities(g, cluster, cluster_neighbors, edge_weights)
			sorted_priorities = sortperm(m1)

			if (isempty(sorted_priorities))
				break
			end

			@label addNode
			highest_priority = pop!(sorted_priorities)
			highest_node = cluster_neighbors[highest_priority]

			push!(cluster, highest_node)
			cluster_edges += count_edges(g, cluster, highest_node)

			d = density(length(cluster), cluster_edges)
			cp = periphery(length(cluster), count_edges(g, cluster, highest_node), d)

			println("Density: $d and CP: $cp")

			if (d < dʼ && cp < cpʼ)
				push!(excluded, pop!(cluster))
				if (!isempty(sorted_priorities))
					println("Remaining Neighbors: $(length(sorted_priorities))")
					@goto addNode
				else
					println("No neighbors left")
					break
				end
			end
		end

		push!(clusters, cluster)
		println("Added $(length(clusters)) clusters.")
		g = rem_nodes(g, cluster)
	end
	[c for c in clusters if length(c) > 2]
end

# ╔═╡ 845006e0-4d7d-11eb-1fe7-1d59f31b1ca6
clusters = density_clustering(g, 0.6, 0.6)

# ╔═╡ d9f94720-4e43-11eb-063f-49c9038aa982
# gplot(g)

# ╔═╡ 822fb7b0-4e69-11eb-3cd4-07f9188ca131
md"""
### Clusters
"""

# ╔═╡ 2710db30-4e6d-11eb-2441-b7d19fe74697
with_terminal() do
    for (i, c) in enumerate(clusters)
		println("Cluster $i: $(convert.(Int, c))")
		println()
	end
end

# ╔═╡ e1bb91a0-4e6d-11eb-1214-09f5dcf2e4a2
md"""
#### Save to File
"""

# ╔═╡ 7f986070-4e68-11eb-16a3-93a29ab63135
output_file = joinpath("output", "clusters.txt")

# ╔═╡ da256032-4e43-11eb-3690-c186bee53f59
open(output_file, "w") do io
	for c in clusters
		println(io, convert.(Int, c))
		println(io)
	end
end

# ╔═╡ 409381b0-4e6e-11eb-1654-21a675b7eb7e
md"""
Save clusters with their original node labels.
"""

# ╔═╡ f1324f20-4e6d-11eb-21fc-e3461700c517
output_file_names = joinpath("output", "named_clusters.txt")

# ╔═╡ 06fcec20-4e6e-11eb-3a44-f754f81f3290
open(output_file_names, "w") do io
	for c in clusters
		println(io, [get_key(i, d) for i in c])
		println(io)
	end
end

# ╔═╡ Cell order:
# ╟─8d2a2140-4d6b-11eb-1ed0-5dfeba2603af
# ╟─a4778db2-4d6b-11eb-0e97-c368743a47ea
# ╠═b9452d10-4d6b-11eb-1074-83294cc14eec
# ╟─b046df70-4e5a-11eb-350f-d1abb711d409
# ╠═822e6af0-4d73-11eb-32e2-216f0814a3c0
# ╠═f0bd31ee-4d77-11eb-2834-85e4981d8680
# ╠═f60d9610-4d75-11eb-002e-6f23a07f66ad
# ╠═097a9b52-4e37-11eb-2afc-0102d9e2c69d
# ╟─30c41780-4e5b-11eb-0508-3b2e824a0607
# ╠═ebca8580-4e3a-11eb-2a8b-711f8bc91692
# ╠═d4c26d30-4e3f-11eb-07d2-2d722ea4c878
# ╠═95d86970-4e45-11eb-0929-692e7fab4639
# ╟─7b666270-4e5b-11eb-1f65-b5eafbc71680
# ╠═3df86fd0-4e3f-11eb-3a9c-1b73eabf5568
# ╟─c063ab80-4e5b-11eb-283f-6749471fae2d
# ╠═ab04f4e0-4e49-11eb-19be-49390f4285e3
# ╟─eba1e8b0-4e57-11eb-1cee-252cd18c350c
# ╠═843b4660-4d7d-11eb-387c-6ddb73c11e73
# ╠═ba80451e-4d6b-11eb-0949-6de777cd3780
# ╠═5655df00-4d76-11eb-01d6-9d21ce3aa823
# ╠═baf68690-4d6b-11eb-3151-b5e91d478309
# ╠═845006e0-4d7d-11eb-1fe7-1d59f31b1ca6
# ╠═d9f94720-4e43-11eb-063f-49c9038aa982
# ╟─822fb7b0-4e69-11eb-3cd4-07f9188ca131
# ╟─2710db30-4e6d-11eb-2441-b7d19fe74697
# ╟─e1bb91a0-4e6d-11eb-1214-09f5dcf2e4a2
# ╠═7f986070-4e68-11eb-16a3-93a29ab63135
# ╠═da256032-4e43-11eb-3690-c186bee53f59
# ╟─409381b0-4e6e-11eb-1654-21a675b7eb7e
# ╠═f1324f20-4e6d-11eb-21fc-e3461700c517
# ╠═06fcec20-4e6e-11eb-3a44-f754f81f3290


text:
# workspace()
using JuLearn,Plots,JLD
plotlyjs()
@Eval lr=1e-2 batchsize=100

include("predict.jl")
include("dataProvider.jl")

inputBase,outputBase=dataProvider()
T=10
net=rnn([:d,:h,:h,:l],[1 1;1 10;10 1;1 1],T,:regression)

# backprop
err=backprop(net,inputBase,outputBase,nIter=30000)
p=plot(err);display(p.o)

# prediction
input=inputBase[:,100+(0:T÷2-1),1]
target=outputBase[:,100+(0:T-1),1]'
prediction=predict(net,input)
prediction=map(i->prediction[i][1],1:T)
p=plot([target prediction],label=["target" "prediction"])
display(p.o)


text:
@show VERSION
VERSION < v"0.4" && exit(1)

if isdefined(Base, :bloop) && !isempty(methods(abs))
    @show abs(bloop)
end

a, b, c = 1, 2, 3
vars = filter(x -> eval(x) isa Integer, names(Main))
println("Integer variables: ", join(vars, ", "), ".")
println("Sum of integers in the global scope: ", sum(eval.(vars)), ".")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode